#《多元统计分析及R语言建模》#
#Multivariate Statistical Analysis and Modeling for R Language#
#王斌会 著#
#暨南大学出版社#
第一章 多元统计分析概述
多元统计分析的内容：
图示法，线性相关分析，线性回归分析，广义与一般线性模型，判别分析，聚类分析，主成分分析，因子分析，对应分析，典型相关分析，多维标度法，综合评价法
第二章 多元数据的数学表达及R应用
x1 <- c(171,175,159,155,152,158,154,164,168,166,159,164)
x2 <- c(57,64,41,38,35,44,41,51,57,49,47,46)
2.1矩阵
rbind(x1,x2) #按行合并形成矩阵
cbind(x1,x2) #按列合并形成矩阵
matrix(data = x1,nrow=3,ncol=4,byrow=F,dimnames=NULL) #，neow为行数，ncol为列数，byrow为是否按行进行
a <- matrix(1:12,nrow = 3,ncol = 4)
b <- matrix(1:12,nrow = 4,ncol = 3)
a%*%b #求矩阵积
diag(a) #取对角线元素
diag(diag(a)) #用对角线元素创建对角矩阵
diag(3) #创建3阶单位阵
solve(A,b) #求Ax=b,b缺省为单位阵，所以不写b的话,可以求A的逆矩阵
eigen(A)#求特征值
chol(A) #Choleskey分解，即A=p’p
svd(A) #奇异值分解，即A=UDV’ ，U'U=V'V=I
qr(A) #QR分解
kronecker(A,B) #kronecker积
dim(A) #矩阵维数
nrow(A) #行数
ncol(A) #列数
rowSums(A) #行列求和，平均
rowMeans(A)
colSums(A)
colMeans(A)
apply(A,1,sum) #A为矩阵，1为行2为列，sum里填想好要的函数操作
2.2数据框
数据框是矩阵的推广
data.frame('身高'=heights,'体重'=weights) 建立数据框同时命名
2.3读取数据
mydata <- read.table("clipboard",header = T)
mydata <- read.table("text.txt",header = T)
mydata <- read.csv("text.csv",header = T)
数据很大，最好不要读入，而使用链接
library(RODBC)
Rcode <- odbcConnectExcel("Rcode.xls") #odbcConnectExcel只认32位系统,且不认xlsx，xlsx用odbcConnectExcel2007
sqlFetch(Rcode,"data") #显示data表格的数据
close(Rcode) #关闭连接
2.4因素分析
data <- read.table("clipboard",header = T)
attach(data) #将命令绑定在data上，可以不用$符号而直接使用列名。 detach()解除绑定
2.4.1 单因素分析
table(年龄)
barplot(table(年龄),col = 1:7)
pie(table(结果))
2.4.2 双因素分析
table(年龄,性别) #二维列联表
barplot(table(年龄,性别),beside= T,col = 1:7)
barplot(table(性别,年龄),beside= T,col = 1:7)
2.4.3 三因素分析
table(年龄,性别,结果)
ftable(年龄,性别,结果)
table(性别,年龄,结果)
ftable(性别,年龄,结果)
ft <- ftable(性别,结果,年龄)
rowSums(ft)
colSums(ft)
sum(ft)
2.5案例分析
Case1 <- read.table("clipboard",header = T)
summary(Case1)
T1 <- table(地区)
T1
barplot(T1)
f <- hist(月收入)
boxplot(月收入~性别)
t.test(月收入~性别)
T2 <- table(性别,观点)
T2
barplot(T2,beside=T)
detach(Case1)
第三章 多元数据的直观表示及R使用
数据：RstatM.xls,d3.1
X <- read.table("clipboard",header = T)
barplot(apply(X,1,mean),las = 3) #按行做均值条形图。las表示坐标是竖着
barplot(apply(X,2,mean)) #按列做均值条形图
barplot(apply(X[2:8],2,mean)) #去掉第一列，按列做均值条形图
barplot(apply(X,2,median)) #按列做中位数条形图
boxplot(X) #箱线图
boxplot(X,horizontal = T) #箱线图,横着
stars(X,full = T,key.loc = c(13,1.5),draw.segments = F) #星象图,full为是否半圆，key.loc 为图例，draw.segments为全圆半圆
调和线图：
library(mvstats)
多元数据实战：
Case2 <- read.table("clipboard",header = T)
RM <- rowMeans(Case2)
barplot(RM,las=3) #行均值
stars(Case2,draw.segments = T,key.loc = c(8,1.6))
第四章 多元相关与回归分析及R使用
书P62
4.1.两变量线性相关系数
cor(x1,x2,method = 'pearson') #两变量间的相关系数，method有‘pearson’‘kendall’‘spearman’，结果大于零，说明成正的相关关系
相关系数的假设检验：
cor.test(x1,x2,alternative='two.sided') #alternative有'two.sided（双侧）','less（左侧）''greater（右侧）'
#Pearson's product-moment correlation
#data: Case2[, 3] and Case2[, 6]
#t = 3.4792, df = 19, p-value = 0.002511
#alternative hypothesis: true correlation is not equal to 0
#95 percent confidence interval:
# 0.262950 0.831572
#sample estimates:
# cor
#0.6238251
# p-value = 0.002511 <0.05 ，证明线性关系。
4.2 一元线性回归分析
yx <- read.table("clipboard",header = T)
fm <- lm(y~x,data = yx)
plot(y~x,data = yx)
abline(fm)
anova(fm) #模型方差分析,看P值
summary(fm) #看P值T值R方
4.3 多元线性回归分析
y=b0+b1*x1+b2*x2+......
每个b：偏回归系数
yX <- read.table("clipboard",header = T)
(fm <- lm(y~x1+x2+x3+x4,data=yX))
anova(fm) #模型方差分析,看P值,发现P<0.0001，模型有意义
summary(fm) #看P值T值R方，x1和x3P值大于0.5，说明x1和x3对y值没有显著地影响
4.4 多元线性相关分析
4.4.1矩阵相关分析
cor(yX)
pairs(yX) #多元数据散点图
4.4.2复相关分析
即R方
4.5 一元非线性回归分析
4.5.1曲线回归
4.5.2可直线化的曲线类型
4.5.1.1多项式曲线
4.5.1.2对数函数
一元曲线回归具体应用
x <- c(1.5,2.8,4.5,7.5,10.5,13.5,15.1,16.5,19.5,22.5,24.5,26.5)
y <- c(7.0,5.5,4.6,3.6,2.9,2.7,2.5,2.4,2.2,2.1,1.9,1.8)
plot(x,y)
1.直线回归
lm.1 <- lm(y~x)
summary(lm.1)$coef #模型系数
summary(lm.1)$r.sq #模型系数(r方),0.79效果不佳
plot(x,y)
abline(lm.1)
2.多项式回归，y=a+bx+cx^2
x1=x
x2=x^2
lm.2 <- lm(y~x1+x2)
summary(lm.2)$coef #模型系数
summary(lm.2)$r.sq #模型系数(r方),0.95还不错
plot(x,y)
lines(x,fitted(lm.2))
3.对数拟合 y=a+bLOGx
lm.log=lm(y~log(x))
summary(lm.log)$coef #模型系数
summary(lm.log)$r.sq #模型系数(r方),0.98
plot(x,y)
lines(x,fitted(lm.log))
4.指数拟合 y=ae^bx
lm.exp=lm(log(y)~x)
summary(lm.exp)$coef #模型系数
summary(lm.exp)$r.sq #模型系数(r方),0.91
plot(x,y)
lines(x,exp(fitted(lm.exp)))
5.幂函数 y=ax^b
总结
lm.1 <- lm(y~x)
x1=x
x2=x^2
lm.2 <- lm(y~x1+x2)
lm.log <- lm(y~log(x))
lm.exp <- lm(log(y)~x)
lm.pow <- lm(log(y)~log(x))
t <- c(summary(lm.1)$r.sq,summary(lm.2)$r.sq,summary(lm.log)$r.sq,summary(lm.exp)$r.sq,summary(lm.pow)$r.sq)
a <- which(t==max(t))
if(a==1){print('线性回归，r.sq=');print(summary(lm.1)$r.sq);plot(x,y);abline(lm.1)}
if(a==2){print('多项式回归，r.sq=');summary(lm.2)$r.sq;plot(x,y);lines(x,fitted(lm.2)) }
if(a==3){print('对数拟合，r.sq=');summary(lm.log)$r.sq;plot(x,y);lines(x,fitted(lm.log))}
if(a==4){print('指数拟合，r.sq=');summary(lm.exp)$r.sq;plot(x,y);lines(x,exp(fitted(lm.exp)))}
if(a==5){print('幂函数拟合，r.sq=');summary(lm.pow)$r.sq;plot(x,y);lines(x,exp(fitted(lm.pow)))}
4.5 多元非线性回归分析
d <- read.table("clipboard",header = T)
model <- nls(Y~A0*(exp(m*t))*(L^a)*(K^b),data=d,start=list(A0=0.45,m=0,a=0.5,b=0.5))
model
多元回归分析例题
case3 <- read.table("clipboard",header = T)
cor(case3)
plot(case3) 相当于pairs(case3)
第五章 广义线性模型与一般线性模型(模型的选用看Ppt的xy的形式)
5.1广义线性模型
5.1.1logistic模型
y变量为0-1时，判断哪些变量会影响y
d <- read.table("clipboard",header = T)
logit.glm <- glm(y~x1+x2+x3,famili=binomial,data=d) #二项分布
summary(logit.glm) #看结果，发现存在未通过的变量，所以做筛选
logit.step <- step(logit.glm,direction="both") #逐步回归，删除无关变量
summary(logit.step) #看结果
pre1 <- predict(logit.step,data.frame(x1=1)) #预测x1为1时的y结果
p1 <- exp(pre1)/(1+exp(pre1)) #计算y概率
pre2 <- predict(logit.step,data.frame(x1=0)) #预测x1为1时的y结果
p2 <- exp(pre2)/(1+exp(pre2)) #计算y概率
c(p1,p2) #0.32 0.65 x1=0时y=1的概率是x1=1时y=1的概率的两倍
5.1.2对数线性模型
d <- read.table("clipboard",header = T)
log.glm <- glm(y~x1+x2,family = poisson(link=log),data=d)
summary(log.glm) #查看p1,p2,小于0.01说明有影响，越小说明对y影响越大
5.2 一般线性模型
可以建立模型，考虑交叉效应
第六章 判别分析
确定型判别 Fisher判别
线性判别分析例：
install.packages("MASS")
library(MASS)
d <- read.table("clipboard",header = T)
plot(d$x1,d$x2)
text(d$x1,d$x2,d$G,adj=-0.5) #x坐标，y坐标，内容，位置调整
ld <- lda(d$G~d$x1+d$x2) #线性判别模型
z <- predict(ld) #预测
newG <- z$class
cbind(d$G,z$x,newG) #结果
tab <- table(d$G,newG) #混淆矩阵
predict(ld,data.frame(x1=8.1,x2=2.0),data=d6.1)
距离判别法例(2类)：
d <- read.table("clipboard",header = T)
attach(d)
plot(Q,C)
text(Q,C,G,adj=-0.5)
plot(Q,P)
text(Q,P,G,adj=-0.5)
library(mvstats)
discrim.dist(cbind(Q,C,P),as.factor(G)) #马氏距离判别模型
discrim.dist(cbind(Q,C,P),as.factor(G),c(8.0,7.5,65)) #新样本判断
detach(d)
多总体距离判别(多类):
d <- read.table("clipboard",header = T)
attach(d)
plot(Q,C);text(Q,C,G,adj=-0.5,cex=0.75)
plot(Q,P);text(Q,P,G,adj=-0.5,cex=0.75)
plot(C,P);text(C,P,G,adj=-0.5,cex=0.75)
异方差距离判别：
D <- discrim.dist(cbind(Q,C,P),as.factor(G))
cbind(G,D=t(D[[1]]),t(D[[2]])) #结果
table(D$G,D$newG) #1判错
sum(diag(prop.table(table(G,newG)))) #正确率
同方差距离判别：
D <- discrim.dist(cbind(Q,C,P),as.factor(G),var.equal=T)
cbind(G,D=t(D[[1]]),t(D[[2]])) #结果
table(D$G,D$newG) #2判错
sum(diag(prop.table(table(G,newG)))) #正确率0.9
线性判别(等方差)：
library(MASS)
ld <- lda(G~Q+C+P)
z <- predict(ld) #预测
newG <- z$class
cbind(G,z$x,newG) #结果
table(G,newG) #混淆矩阵
sum(diag(prop.table(table(G,newG)))) #正确率0.9
detach(d)
概率型判别 Bayes判别
先验概率相等，等价于Fisher线性判别分析
ld1 <- lda(G~Q+C+P,piror=c((1,1,1)/3))
z1 <- predict(ld1)
cbind(G,z1$x,z1$class)
table(G,z1$class)
z1$post #给出分类的概率，这与Fisher判别不同
先验概率不相等
ld1 <- lda(G~Q+C+P,piror=c((5,8,7)/20))
z2 <- predict(ld2)
cbind(G,z2$x,z2$class)
table(G,z2$class)
第七章 聚类分析
Q型聚类：纵向样本数分类
样本的统计量叫 距离
R型聚类：横向决策变量分类，降维
变量的统计量叫 相似系数
系统聚类
例子
library(mvstats)
d <- read.table("clipboard",header = T) #7.2
plot(d)
H.clust(d,"euclidean","single",plot = T) #最短距离法
H.clust(d,"euclidean","complete",plot = T) #最长距离法
H.clust(d,"euclidean","median",plot = T) #中间距离法
H.clust(d,"euclidean","average",plot = T) #类平均法
H.clust(d,"euclidean","centroid",plot = T) #重心法
H.clust(d,"euclidean","ward",plot = T) #ward法
kmeans聚类(快速聚类)(需要给定类数)
1.做数据
x1 <- matrix(rnorm(1000,mean=0,sd=0.3),ncol=10) #中心在0
x2 <- matrix(rnorm(1000,mean=1,sd=0.3),ncol=10) #中心在1
x <- rbind(x1,x2)
H.clust(x,"euclidean","complete") #系统聚类
cl <- kmeans(x,2) #kmeans聚类聚类
cl #看结果
pch1 <-rep("1",100)
pch2 <-rep("2",100)
plot(x,col=cl$cluster,pch=c(pch1,pch2),cex=0.7)
points(cl$centers,col=3,pch="*",cex=3)
第八章 主成分分析
变量降维方法
例
X <- read.table("clipboard",header = T)
cor(X)
PCA <- princomp(X,cor = T) #主成分分析
PCA
summary(PCA)
PCA$loadings #主成分载荷,将累计载荷大于0.8的，选定为主成分
screeplot(PCA,type = "lines")
PCA$scores #主成分得分
通过看主成分载荷矩阵的绝对值大小，可以了解到选定的主成分，代表了原变量的意义
library(mvstats)
princomp.rank(PCA,m=2,plot=T) #主成分排名
第九章 因子分析
原变量之间的联系与区别
例
X <- read.table("clipboard",header = T)
cor(X)
(FA0 <- factanal(X,3,rot = "none")) #极大似然法因子分析
library(mvstats)
(Fac <- factpc(X,3)) #主成分法因子分析
(Fa1 <- factanal(X,3,rot = "varimax")) #varimax法旋转因子分析
#Call:
#factanal(x = X, factors = 3, rotation = "varimax")
#Uniquenesses:
# x1 x2 x3 x4 x5 x6
#0.005 0.005 0.005 0.271 0.005 0.548
#Loadings:
# Factor1 Factor2 Factor3
#x1 0.983 0.155 Factor1代表盈利能力
#x2 0.985 0.142 Factor2代表偿债能力
#x3 -0.990 -0.124 Factor3代表企业发展
#x4 0.127 0.844
#x5 0.293 0.953
#x6 0.210 0.631
#
# Factor1 Factor2 Factor3
#SS loadings 1.998 1.800 1.367
#Proportion Var 0.333 0.300 0.228
#Cumulative Var 0.333 0.633 0.861
#
#The degrees of freedom for the model is 0 and the fit was 1.1422
因子得分：
Fa1 <- factanal(X,3,scores = "regression" ) #回归估计发的极大似然法因子分析,如果数据为多元正态分布时，用这个。否则用下面的主成分法进行计算。
Fa1$scores
Fac1 <- factpc(X,3,scores = "regression" ) #回归估计发的主成分法因子分析
Fac1$scores
因子得分排名：
factanal.rank(Fa1,plot = T)
第十章 对应分析
样本，变量之间的关系
x <- read.table("clipboard",header = T)
library(MASS)
ca2 <- corresp(x,nf=2) #对应分析
ca2
biplot(ca2)
abline(v=0,h=0,lty=3)
第十一章 典型相关关系
两组变量之间的相互关系
例1
X <- read.table("clipboard",header = T) #d11.1
(R <- cor(X))
R11 <- R[1:3,1:3]
R12 <- R[1:3,4:6]
R21 <- R[4:6,1:3]
R22 <- R[4:6,4:6]
A <- solve(R11)%*%R12%*%solve(R22)%*%R21
ev <- eigen(A)$values
ev #典型相关系数
xy <- scale(X)
ca <- cancor(xy[ ,1:3],xy[ ,4:6]) #典型相关分析
ca$cor #典型相关系数
ca$xcoef #第一组变量典型载荷
ca$ycoef #第二组变量典型载荷
library(mvstats)
cancor.test(xy[ ,1:3],xy[ ,4:6],plot = T) 典型相关分析及作图
例2
d <- read.table("clipboard",header = T) #d11.2
library(mvstats)
cancor.test(d[ ,1:4],d[ ,5:10],plot = T)
第十二章 多维标度法
例
D <- read.table("clipboard",header = T) #12.1
library(MASS)
D <- as.matrix(D)
fit <- isoMDS(D,k=2)
fit
x <- fit$points[ ,1]
y <- fit$points[ ,2]
plot(x,y,type="n")
text(x,y,labels=row.names(D))
第十三章 综合评价方法
用一个指标说明被评价对象的总体水平
13.1综合评分法
加权平均分
B1 <- read.table("clipboard",header = T) #13.1
B1_z <- z_data(B1)
B1_z
Si <- apply(B1_z,mean)
cbind(B1_z,Si)
cbind(Si=Si,ri,rank(-Si))
13.2层次分析法
B1 <- c(1,4,5,3,6,7,1/4,1,2,1/2,3,4,1/5,1/2,1,1/3,2,3,1/3,2,3,1,4,5,1/6,1/3,1/2,1/4,1,2,1/7,1/4,1/3,1/5,1/2,1)
B1_W <- weight(B1) #B1权重
B1_W
CI_CR(B1) #一致性检验
S_rank(B1_Z，B1——W) #按B1得到的综合得分及排名
打赏
微信扫一扫，打赏作者吧～
