R 语言 主成分分析

1.R中的主成分和因子分析

R的基础安装包中提供了PCA和EFA的函数，分别为princomp （）和factanal（）

psych包中有用的因子分析函数

函数	描述
principal（）	含多种可选的方差放置方法的主成分分析
fa（）	可用主轴、最小残差、加权最小平方或最大似然法估计的因子分析
fa.parallel（）	含平等分析的碎石图
factor.plot（）	绘制因子分析或主成分分析的结果
fa.diagram（）	绘制因子分析或主成分分析的载荷矩阵
scree（）	因子分析和主成分分析的碎石图

PCA/EFA 分析流程：

（1）数据预处理；PCA和EFA都是根据观测变量间的相关性来推导结果。用户可以输入原始数据矩阵或相关系数矩阵列到principal()和fa（）函数中，若输出初始结果，相关系数矩阵将会被自动计算，在计算前请确保数据中没有缺失值；

（2）选择因子分析模型。判断是PCA（数据降维）还是EFA（发现潜在结构）更符合你的分析目标。若选择EFA方法时，还需要选择一种估计因子模型的方法（如最大似然估计）。

（3）判断要选择的主成分/因子数目；

（4）选择主成分/因子；

（5）旋转主成分/因子；

（6）解释结果；

（7）计算主成分或因子得分。

#主成分分析：psych包的：princomp()#因子分析：factanal()
library(psych)
#例1（根据原始变量做主成分分析）
View(USJudgeRatings)
#主成分个数选择（图形方法选择主成分个数）
fa.parallel(USJudgeRatings[,-1], fa="pc", n.iter=100,
            show.legend=FALSE, main="Scree plot with parallel analysis")
#主成分结果
pc <- principal(USJudgeRatings[,-1], nfactors=1,scores=T)
#每个调查对象在该主成分上的得分
head(pc$scores)
#律师与法官的接触频数与法官评分间的相关系数
cor(USJudgeRatings$CONT, pc$score)
cor(USJudgeRatings)


#例2（根据相关关系图做主成分分析）
View(Harman23.cor)
#主成分个数选择
fa.parallel(Harman23.cor$cov, n.obs=302, fa="pc", n.iter=100,
            show.legend=FALSE, main="Scree plot with parallel analysis")
#主成分结果
principal(Harman23.cor$cov, nfactors=2, rotate="none")
#主成分旋转(方差极大旋转)
principal(Harman23.cor$cov, nfactors=2, rotate="varimax")
#获取主成分得分系数
rc <- principal(Harman23.cor$cov, nfactors=2, rotate="varimax")
round(unclass(rc$weights), 2)
PC1 = 0.28*height + 0.30*arm.span + 0.30*forearm + 0.29*lower.leg -
      0.06*weight - 0.08*bitro.diameter - 0.10*chest.girth -
      0.04*chest.width
PC2 = -0.05*height - 0.08*arm.span - 0.09*forearm - 0.06*lower.leg +
       0.33*weight + 0.32*bitro.diameter + 0.34*chest.girth +
       0.27*chest.width

#例3
#默认生成所有主成分查看主成分的贡献率，然后选择
pr <- princomp(USJudgeRatings[,-1],cor = T)
options(digits = 4)
summary(pr,loadings = T)#列出主成分对应原始变量的系数
#简单画碎石图的函数
screeplot(pr,type = "line",main = "碎石图")
#用biplot绘制主成分的散点图
pr <- princomp(USJudgeRatings[,-1],cor = T)
biplot(pr,choices = 1:2)

#R 中labdsv包
library(labdsv)
#
pac <- pca(USJudgeRatings[,-1],dim = 4,cor = T)
summary(pac)
#查看载荷系数
loadings.pca(pac)
#绘制碎石图以及累计方差图
#分割图形区域
op <- par(mfrow= c(1,2))
varplot.pca(pac)
Hit Return to Continue