1. 前言
主成分分析是一种数据降维的技巧，能将大量相关变量转化为一组很少的不想管变量，这些不相关变量称之为主成分，因子分析是一系列用来发现一组变量的潜在结构的方法，通过寻找一组更小的，潜在的结构解释观测到的，显式的变量间的关系.
主成分分析是观测变量的线性组合，形成线性组合的全中都是通过最大化个主成分所解释的方差来获取，同时要保证各主成分各不相关.
因子分析是被当做是观测变量的结构基础或是原因，而不是线性组合，代表观测变量方差的误差无法用因子解释，相关的因子和误差无法直接观测，但可通过变量间的相关关系退到得到.
相关经验法则：因子分许需要5—10倍与变量数的样本数才可用.
2. R中的主成分和因子分析
R中的基础函数princomp()函数和factanal()函数，也可采用psych包中的函数
psych包中的因子分析函数
|
函数 |
描述 |
|
principal() |
含多种可选的方差旋转阀的主成分分析 |
|
fa() |
可用主轴，最小残差，加权最小平方或最大似然法估计的因子分析 |
|
fa.parallel() |
含平行分析的碎石图 |
|
factor.plot() |
绘制因子分析或主成分分析的结果 |
|
fa.diagram() |
绘制因子分析或主成分的载荷矩阵 |
|
scree() |
因子分析和主成分分析的碎石图 |
采用这两种方法需要一些步骤才可获得最终结果，常见的步骤如下：
（1）：数据预处理：主成分分析和因子分析都根据观测变量间的相关性来退到结果，可将原始数据矩阵或者相关系数矩阵输入到principao()和fa()函数中，若输入的是初始矩阵，则相关系数矩阵会自动计算，但是要确保数据中没有存在缺失值.
（2）：选择因子模型，判断是因子分析还是主成分分析，若是因子分析，则要选择一种估计因子模型的方法，例如有最大似然估计法
（3）：判断要选择的主成分/因子数目
（4）：选择主成分/因子
（5）：旋转主成分/因子
（6）:解释结果
（7）：计算主成分或因子得分
3.主成分分析
主成分是讲一组较少的不想管变量替代大量相关变量，同时尽可能保留初始变量的信息.采用数据集USJudgeRatings中美国高等法官法官的评分进行主成分分析.用较少的变量来替代总结这11个变量，首先需要进行的是判断主成分的个数，准则如下：
(1):根据先验经验和理论知识判断主成分输
(2):跟要要解释的变量方差的累计值的阈值来判断所需主成分数
(3):通过检查变量间K*K的相关系数矩阵来判断保留的主成分数
较常见是基于特征值的方法，每个主成分都与相关系数矩阵的特征值相关联，第一主成分与最大的特征值相关联，第二主成分与第二大的特征值相关联,KH准则建议保留特征值大于1的主成分，碎石图是绘制了特征值与主成分数的图形，在图形变化最大处之上的主成分皆可保留，最后，仍需进行模拟，依据初始矩阵大小的随机数据矩阵来判断要提取的特征值.若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值，该主成分可以保留，该方法称为平行分析.
利用fa.parallel()函数，同时对三种特征值判别准则进行评价.
fa.parallel(USJudgeRatings[, -1], fa ="pc", n.iter = 100
+ show.legend = FALSE, main = "Scree plot with parallelanalysis")
由图展示了基于观测特征值的碎石检验（即含有X的线段），根据100个随机数据矩阵退到出来的特征值均值（虚线）及大于1的特征值准则.三种准则表明选择一个胡成分可保留数据集的大部分信息，下步为使用principal()函数提取主函数.其格式为：
principai(r,nfactors=,rotate=,scores=)
其中r是相关系数矩阵或者原始矩阵，nffactors是主成分数（默认值为1），rotate是指旋转方法（默认方法是最大方差旋转（varimax）），scores设定是否需要计算主成分的得分，默认值是否.n.iter是指随机矩阵的个数
library(psych)
pc<-principal(USJudgeRatings[,-1],nfactors=1,rotate="varimax",scores=T)
pc
Principal Components Analysis
Call: principal(r = USJudgeRatings[, -1],nfactors = 1, rotate = "varimax",
scores = T)
Standardized loadings (pattern matrix) basedupon correlation matrix
PC1 h2 u2 com
INTG 0.92 0.84 0.1565 1
DMNR 0.91 0.83 0.1663 1
DILG 0.97 0.94 0.0613 1
CFMG 0.96 0.93 0.0720 1
DECI 0.96 0.92 0.0763 1
PREP 0.98 0.97 0.0299 1
FAMI 0.98 0.95 0.0469 1
ORAL 1.00 0.99 0.0091 1
WRIT 0.99 0.98 0.0196 1
PHYS 0.89 0.80 0.2013 1
RTEN 0.99 0.97 0.0275 1
PC1
SS loadings 10.13
Proportion Var 0.92
Mean item complexity = 1
Test of the hypothesis that 1 component issufficient.
The root mean square of the residuals (RMSR)is 0.04
withthe empirical chi square 6.21 with prob < 1
Fit based upon off diagonal values = 1
输入的变量没有最后一个变量的原始数据，并制定获取一个未旋转的主成分，PC栏即是包含了成分载荷，即是指观测变量与主成分的相关系数，若提取的主成分不知一个，那么会有pc2,pc3等栏.
成分载荷即表示主成分的与变量的相关关系，由上述结果可得，第一主成分与每个变量高度相关.
h2指公因子方差，即每个主成分对每个变量的方差解释程度.
u2值得是成分唯一性，即方差无法被主成分解释的比例.
SSloading行包含了与主成分关联的特征值，指特定主成分关联后标准化后的方差值.Prosprtion Var行表示每个主成分对整个数据的解释程度.本例中，0.92即表示第一主成分解释了整个数据集的0.92.
第二个例子是Harman23.cor看女生身体的8个身体指标.
fa.parallel(Harman23.cor$cov,n.obs=302,fa="pc",n.iter=100,show.legend= F,main="scree plot with parallel analysis")
n.obs指代的是样本量的大小，采用的数据是Harman23.cor的协方差矩阵.
由碎石图可看，本例取得了两个主成分.
Principal Components Analysis
Call: principal(r = Harman23.cor$cov,nfactors = 2, rotate = "none")
Standardized loadings (pattern matrix) basedupon correlation matrix
PC1 PC2 h2 u2 com
height 0.86 -0.37 0.88 0.123 1.4
arm.span 0.84 -0.44 0.90 0.097 1.5
forearm 0.81 -0.46 0.87 0.128 1.6
lower.leg 0.84 -0.40 0.86 0.139 1.4
weight 0.76 0.52 0.85 0.150 1.8
bitro.diameter 0.67 0.53 0.74 0.261 1.9
chest.girth 0.62 0.58 0.72 0.283 2.0
chest.width 0.67 0.42 0.62 0.375 1.7
PC1 PC2
SS loadings 4.67 1.77
Proportion Var 0.58 0.22
Cumulative Var 0.58 0.81
Proportion Explained 0.73 0.27
Cumulative Proportion 0.73 1.00
Mean item complexity = 1.7
Test of the hypothesis that 2 components aresufficient.
The root mean square of the residuals (RMSR)is 0.05
Fit based upon off diagonal values = 0.99
由PC1，PC2解释了0.58，第二主成分解释了0.22，二者共解释了88%的方差.
2.3主成分旋转
旋转的方法分为正交旋转（使成分间间不相关）和斜交旋转（是成分间变得相关）
rc<-principal(Harman23.cor$cov,nfactors=2,rotate="varimax")
Principal Components Analysis
Call: principal(r = Harman23.cor$cov,nfactors = 2, rotate = "varimax")
RC1 RC2 h2 u2 com
height 0.90 0.25 0.88 0.123 1.2
arm.span 0.93 0.19 0.90 0.097 1.1
forearm 0.92 0.16 0.87 0.128 1.1
lower.leg 0.90 0.22 0.86 0.139 1.1
weight 0.26 0.88 0.85 0.150 1.2
bitro.diameter 0.19 0.84 0.74 0.261 1.1
chest.girth 0.11 0.84 0.72 0.283 1.0
chest.width 0.26 0.75 0.62 0.375 1.2
RC1 RC2
SS loadings 3.52 2.92
Proportion Var 0.44 0.37
Cumulative Var 0.44 0.81
Proportion Explained 0.55 0.45
Cumulative Proportion 0.55 1.00
Mean item complexity = 1.1
Test of the hypothesis that 2 components aresufficient.
The root mean square of the residuals (RMSR)is 0.05
Fit based upon off diagonal values = 0.99
此处采用的是方差极大旋转法，可见第一主成分分别解释前四个变量，第二个主成分解释了后四个变量，总共解释方差解释了81%.
2.4获取主成分
利用Principal，获得每个对象在该主成分的得分
pc<-principal(USJudgeRatings[,-1],nfactoes=1,rotate="none")
head(pc$scores)
PC1
AARONSON,L.H. -0.1857981
ALEXANDER,J.M. 0.7469865
ARMENTANO,A.J. 0.0704772
BERDON,R.I. 1.1358765
BRACKEN,J.J. -2.1586211
BURNS,E.B. 0.7669406
此段代码可得不同的对象在某个主成分上的得分，也可看到不同的变量与其评分之间的相关系数.
cor(USJudgeRatings$CONT,pc$scores)
PC1
[1,] -0.008815895
可得，与该变量与该得分的关系较低.
对于不同变量与成分，可用如下得到不同成分的得分系数.根据下式的矩阵，可得到不同变量在不同成分上的得分，且都基于山体指标已经标准化的前提.(均值为0，方差为1)，但是倘若是要探索可解释的因子，采用的是因子分析.
round(unclass(rc$weights),2)
RC1 RC2
height 0.28 -0.05
arm.span 0.30 -0.08
forearm 0.30 -0.09
lower.leg 0.28 -0.06
weight -0.06 0.33
bitro.diameter -0.08 0.32
chest.girth -0.10 0.34
chest.width -0.04 0.27
3. 探索性因子分析
因子分析是发现数据中较少，无法观测变量，来解释可观测变量的相关性，这些无法观测到的变量称为因子.
再次采用的是能力数据的因子分析.
cov2cor()函数将协方差矩阵给转化为相关系数矩阵，不过要求数据集不能有缺失值.
covariance<-ability.cov$cov
correlations<-cov2cor(covariance)
correlations
general picture blocks maze reading vocab
general 1.00 0.47 0.55 0.34 0.58 0.51
picture 0.47 1.00 0.57 0.19 0.26 0.24
blocks 0.55 0.57 1.00 0.45 0.35 0.36
maze 0.34 0.19 0.45 1.00 0.18 0.22
reading 0.58 0.26 0.35 0.18 1.00 0.79
vocab 0.51 0.24 0.36 0.22 0.79 1.00
计算完相关系数矩阵后，下步工作为判断需要提取多少个因子.
3.1判断提取的共民族数
用fa.parallel()函数判断提取的公因子数
3.2提取公因子
决定提取两个因子时，采用fa()函数获得相应的结果，格式如下：
fa(r,nfactors=,n.obs=,rotate=,scores=,fm=)
其中，r是指相关系数矩阵或者原始矩阵，nfactor设定身体去的因子数目，且默认值为1，n.obs为观测数目，当输入的矩阵为相关系数矩阵时需要填写该项，rotate值得是旋转的方法，在因子分析中默认的是互变异数最小法，scores为是否计算因子得分，fm为设定因子话方法，且系统默认值为极小残差法.
与主成分分析法不同，提取公因子的方法较多，包括最大似然法（ml），主轴迭代法（pa），加权最小二乘法（wls），广义加权最小二乘法（gls）,和最小残差法，用较多为极大似然法，但有时不会收敛.
fa<-fa(correlations,nfactors=2,rotate="none",fm="pa")
fa
Factor Analysis using method = pa
Call: fa(r = correlations, nfactors = 2,rotate = "none", fm = "pa")
Standardized loadings (pattern matrix) basedupon correlation matrix
PA1 PA2 h2 u2 com
general 0.75 0.07 0.57 0.432 1.0
picture 0.52 0.32 0.38 0.623 1.7
blocks 0.75 0.52 0.83 0.166 1.8
maze 0.39 0.22 0.20 0.798 1.6
reading 0.81 -0.51 0.91 0.089 1.7
vocab 0.73 -0.39 0.69 0.313 1.5
PA1 PA2
SS loadings 2.75 0.83
Proportion Var 0.46 0.14
Cumulative Var 0.46 0.60
Proportion Explained 0.77 0.23
Cumulative Proportion 0.77 1.00
Mean item complexity = 1.5
Test of the hypothesis that 2 factors aresufficient.
The degrees of freedom for the null modelare 15 and the objective function was 2.5
The degrees of freedom for the model are4 and the objective function was 0.07
The root mean square of the residuals (RMSR)is 0.03
The df corrected root mean square of theresiduals is 0.06
Fit based upon off diagonal values = 0.99
Measures of factor score adequacy
PA1 PA2
Correlation of scores with factors 0.96 0.92
Multiple R square of scores with factors 0.93 0.84
Minimum correlation of possible factorscores 0.86 0.68
可得，两个因子解释了6个心理学测验的60%的方差.
3.3因子旋转
尝试不同的旋转方法来探索因子分析.
fa.varimax<-fa(correlations,nfactors=2,rotate="varimax",fm="pa")
fa.varimax
Factor Analysis using method = pa
Call: fa(r = correlations, nfactors = 2,rotate = "varimax", fm = "pa")
Standardized loadings (pattern matrix) basedupon correlation matrix
PA1 PA2 h2 u2 com
general 0.49 0.57 0.57 0.432 2.0
picture 0.16 0.59 0.38 0.623 1.1
blocks 0.18 0.89 0.83 0.166 1.1
maze 0.13 0.43 0.20 0.798 1.2
reading 0.93 0.20 0.91 0.089 1.1
vocab 0.80 0.23 0.69 0.313 1.2
PA1 PA2
SS loadings 1.83 1.75
Proportion Var 0.30 0.29
Cumulative Var 0.30 0.60
Proportion Explained 0.51 0.49
Cumulative Proportion 0.51 1.00
Mean item complexity = 1.3
可见不同的变量在不同的因素有着差别较大的载荷，正交旋转阀是强制让两个因子不想管，而允许两个变量相关则采用斜交旋转法，例如promax.
fa.promax <-fa(correlations, nfactors = 2, rotate = "promax", fm ="pa")
fa.promax
Factor Analysis using method = pa
Call: fa(r = correlations, nfactors = 2,rotate = "promax", fm = "pa")
Warning: A Heywood case was detected.
Standardized loadings (pattern matrix) basedupon correlation matrix
PA1 PA2 h2 u2 com
general 0.37 0.48 0.57 0.432 1.9
picture -0.03 0.63 0.38 0.623 1.0
blocks -0.10 0.97 0.83 0.166 1.0
maze 0.00 0.45 0.20 0.798 1.0
reading 1.00 -0.09 0.91 0.089 1.0
vocab 0.84 -0.01 0.69 0.313 1.0
PA1 PA2
SS loadings 1.83 1.75
Proportion Var 0.30 0.29
Cumulative Var 0.30 0.60
Proportion Explained 0.51 0.49
Cumulative Proportion 0.51 1.00
Withfactor correlations of
PA1 PA2
PA1 1.00 0.55
PA2 0.55 1.00
探讨正交旋转和斜交旋转的不同：正交旋转侧重于因子结构矩阵（变量与因子的相关系数），斜交矩阵会考虑，因子结构矩阵，因子模式矩阵（标准化回归系数矩阵，其列出了因子预测变量的权重），因子关联矩阵.（因子相关系数矩阵），在斜交旋转中，采用的因子模式矩阵，是标准化的回归系数，而并非相关系数，而矩阵的列用来对因子进行命名.而因子关联矩阵中因子相关数达到了0.55，若相关数较低，就采取正交旋转.采用如下函数可用来计算因子载荷矩阵，其中F=p*phi，F即是因子载荷矩阵，P为因子模式矩阵，Phi是因子关联矩阵.
fsm<- function(oblique) {
if (class(oblique)[2]=="fa" & is.null(oblique$Phi)) {
warning("Object doesn't look like oblique EFA")
} else {
P <- unclass(oblique$loading)#显示变量内部数据
F <- P %*% oblique$Phi
colnames(F) <- c("PA1", "PA2")
return(F)
}
}
fsm(fa.promax)
PA1 PA2
general 0.64 0.69
picture 0.32 0.61
blocks 0.43 0.91
maze 0.25 0.45
reading 0.95 0.46
vocab 0.83 0.45
此即为因子载荷矩阵，与正交旋转的因子载荷矩阵相比，斜交的载荷矩阵列的噪声较大，这是因为允许潜在因子相关，随斜交方法更复杂，但是模型将更符合真实数据.
使用factor.plot()或者fa.diagram()函数，可以绘制正交或者斜交结果的图形.
factor.plot(fa.promax,labels=rownames(fa.promax$loadings))
fa.diagram(fa.promax)
3.4因子得分
在fa()函数中添加score=T,即可计算因子得分了.