《R语言实战》自学笔记71-主成分和因子分析

第14章主成分和因子分析

主成分分析
主成分分析（(Principal Component Analysis，PCA）是一种数据降维技巧，它能将大量相关变量转化为一组很少的不相关变量，这些无关变量称为主成分（原来变量的线性组合）。整体思想就是化繁为简，抓住问题关键，也就是降维思想。
主成分分析法是通过恰当的数学变换，使新变量——主成分成为原变量的线性组合，并选取少数几个在变差总信息量中比例较大的主成分来分析事物的一种方法。主成分在变差信息量中的比例越大，它在综合评价中的作用就越大。

因子分析
探索性因子分析法（Exploratory Factor Analysis，EFA）是一系列用来发现一组变量的潜在结构的方法。它通过寻找一组更小的、潜在的或隐藏的结构来解释已观测到的、显式的变量间的关系。

PCA与EFA模型间的区别
参见图14-1。主成分（PC1和PC2）是观测变量（X1到X5）的线性组合。形成线性组合的权重都是通过最大化各主成分所解释的方差来获得，同时还要保证个主成分间不相关。相反，因子（F1和F2）被当做是观测变量的结构基础或“原因”，而不是它们的线性组合。

image.png

14.1 R中的主成分和因子分析

R的基础安装包提供了PCA和EFA的函数，分别为princomp()和factanal()。
最常见的分析步骤
(1)数据预处理。PCA和EFA都根据观测变量间的相关性来推导结果。用户可以输入原始数据矩阵或者相关系数矩阵到principal()和fa()函数中。若输入初始数据，相关系数矩阵将会被自动计算，在计算前请确保数据中没有缺失值。
(2)选择因子模型。判断是PCA（数据降维）还是EFA（发现潜在结构）更符合你的研究目标。如果选择EFA方法，你还需要选择一种估计因子模型的方法（如最大似然估计）。
(3)判断要选择的主成分/因子数目。
(4)选择主成分/因子。
(5)旋转主成分/因子。
(6)解释结果。
(7)计算主成分或因子得分。

image.png

14.2 主成分分析

PCA的目标是用一组较少的不相关变量代替大量相关变量，同时尽可能保留初始变量的信息，这些推导所得的变量称为主成分，它们是观测变量的线性组合。如第一主成分为：
$PC_1=a_1X_1=a_2X_2+.....+a_kX_k$
它是k个观测变量的加权组合，对初始变量集的方差解释性最大。第二主成分也是初始变量的线性组合，对方差的解释性排第二，同时与第一主成分正交（不相关）。后面每一个主成分都最大化它对方差的解释程度，同时与之前所有的主成分都正交。理论上来说，你可以选取与变量数相同的主成分，但从实用的角度来看，我们都希望能用较少的主成分来近似全变量集。

主成分与原始变量之间的关系
（1）主成分保留了原始变量绝大多数信息。
（2）主成分的个数大大少于原始变量的数目。
（3）各个主成分之间互不相关。
（4）每个主成分都是原始变量的线性组合。

数据集USJudgeRatings包含了律师对美国高等法院法官的评分。数据框包含43个观测，12个变量。

USJudgeRatings # 查看数据集。
##                 CONT INTG DMNR DILG CFMG DECI PREP FAMI ORAL WRIT PHYS RTEN
## AARONSON,L.H.    5.7  7.9  7.7  7.3  7.1  7.4  7.1  7.1  7.1  7.0  8.3  7.8
## ALEXANDER,J.M.   6.8  8.9  8.8  8.5  7.8  8.1  8.0  8.0  7.8  7.9  8.5  8.7
## ARMENTANO,A.J.   7.2  8.1  7.8  7.8  7.5  7.6  7.5  7.5  7.3  7.4  7.9  7.8
## BERDON,R.I.      6.8  8.8  8.5  8.8  8.3  8.5  8.7  8.7  8.4  8.5  8.8  8.7
## BRACKEN,J.J.     7.3  6.4  4.3  6.5  6.0  6.2  5.7  5.7  5.1  5.3  5.5  4.8
## BURNS,E.B.       6.2  8.8  8.7  8.5  7.9  8.0  8.1  8.0  8.0  8.0  8.6  8.6
## CALLAHAN,R.J.   10.6  9.0  8.9  8.7  8.5  8.5  8.5  8.5  8.6  8.4  9.1  9.0
## COHEN,S.S.       7.0  5.9  4.9  5.1  5.4  5.9  4.8  5.1  4.7  4.9  6.8  5.0
## DALY,J.J.        7.3  8.9  8.9  8.7  8.6  8.5  8.4  8.4  8.4  8.5  8.8  8.8
## DANNEHY,J.F.     8.2  7.9  6.7  8.1  7.9  8.0  7.9  8.1  7.7  7.8  8.5  7.9
## DEAN,H.H.        7.0  8.0  7.6  7.4  7.3  7.5  7.1  7.2  7.1  7.2  8.4  7.7
## DEVITA,H.J.      6.5  8.0  7.6  7.2  7.0  7.1  6.9  7.0  7.0  7.1  6.9  7.2
## DRISCOLL,P.J.    6.7  8.6  8.2  6.8  6.9  6.6  7.1  7.3  7.2  7.2  8.1  7.7
## GRILLO,A.E.      7.0  7.5  6.4  6.8  6.5  7.0  6.6  6.8  6.3  6.6  6.2  6.5
## HADDEN,W.L.JR.   6.5  8.1  8.0  8.0  7.9  8.0  7.9  7.8  7.8  7.8  8.4  8.0
## HAMILL,E.C.      7.3  8.0  7.4  7.7  7.3  7.3  7.3  7.2  7.1  7.2  8.0  7.6
## HEALEY.A.H.      8.0  7.6  6.6  7.2  6.5  6.5  6.8  6.7  6.4  6.5  6.9  6.7
## HULL,T.C.        7.7  7.7  6.7  7.5  7.4  7.5  7.1  7.3  7.1  7.3  8.1  7.4
## LEVINE,I.        8.3  8.2  7.4  7.8  7.7  7.7  7.7  7.8  7.5  7.6  8.0  8.0
## LEVISTER,R.L.    9.6  6.9  5.7  6.6  6.9  6.6  6.2  6.0  5.8  5.8  7.2  6.0
## MARTIN,L.F.      7.1  8.2  7.7  7.1  6.6  6.6  6.7  6.7  6.8  6.8  7.5  7.3
## MCGRATH,J.F.     7.6  7.3  6.9  6.8  6.7  6.8  6.4  6.3  6.3  6.3  7.4  6.6
## MIGNONE,A.F.     6.6  7.4  6.2  6.2  5.4  5.7  5.8  5.9  5.2  5.8  4.7  5.2
## MISSAL,H.M.      6.2  8.3  8.1  7.7  7.4  7.3  7.3  7.3  7.2  7.3  7.8  7.6
## MULVEY,H.M.      7.5  8.7  8.5  8.6  8.5  8.4  8.5  8.5  8.4  8.4  8.7  8.7
## NARUK,H.J.       7.8  8.9  8.7  8.9  8.7  8.8  8.9  9.0  8.8  8.9  9.0  9.0
## O'BRIEN,F.J.     7.1  8.5  8.3  8.0  7.9  7.9  7.8  7.8  7.8  7.7  8.3  8.2
## O'SULLIVAN,T.J.  7.5  9.0  8.9  8.7  8.4  8.5  8.4  8.3  8.3  8.3  8.8  8.7
## PASKEY,L.        7.5  8.1  7.7  8.2  8.0  8.1  8.2  8.4  8.0  8.1  8.4  8.1
## RUBINOW,J.E.     7.1  9.2  9.0  9.0  8.4  8.6  9.1  9.1  8.9  9.0  8.9  9.2
## SADEN.G.A.       6.6  7.4  6.9  8.4  8.0  7.9  8.2  8.4  7.7  7.9  8.4  7.5
## SATANIELLO,A.G.  8.4  8.0  7.9  7.9  7.8  7.8  7.6  7.4  7.4  7.4  8.1  7.9
## SHEA,D.M.        6.9  8.5  7.8  8.5  8.1  8.2  8.4  8.5  8.1  8.3  8.7  8.3
## SHEA,J.F.JR.     7.3  8.9  8.8  8.7  8.4  8.5  8.5  8.5  8.4  8.4  8.8  8.8
## SIDOR,W.J.       7.7  6.2  5.1  5.6  5.6  5.9  5.6  5.6  5.3  5.5  6.3  5.3
## SPEZIALE,J.A.    8.5  8.3  8.1  8.3  8.4  8.2  8.2  8.1  7.9  8.0  8.0  8.2
## SPONZO,M.J.      6.9  8.3  8.0  8.1  7.9  7.9  7.9  7.7  7.6  7.7  8.1  8.0
## STAPLETON,J.F.   6.5  8.2  7.7  7.8  7.6  7.7  7.7  7.7  7.5  7.6  8.5  7.7
## TESTO,R.J.       8.3  7.3  7.0  6.8  7.0  7.1  6.7  6.7  6.7  6.7  8.0  7.0
## TIERNEY,W.L.JR.  8.3  8.2  7.8  8.3  8.4  8.3  7.7  7.6  7.5  7.7  8.1  7.9
## WALL,R.A.        9.0  7.0  5.9  7.0  7.0  7.2  6.9  6.9  6.5  6.6  7.6  6.6
## WRIGHT,D.B.      7.1  8.4  8.4  7.7  7.5  7.7  7.8  8.2  8.0  8.1  8.3  8.1
## ZARRILLI,K.J.    8.6  7.4  7.0  7.5  7.5  7.7  7.4  7.2  6.9  7.0  7.8  7.1
str(USJudgeRatings) # 查看数据集结构。
## 'data.frame':    43 obs. of  12 variables:
##  $ CONT: num  5.7 6.8 7.2 6.8 7.3 6.2 10.6 7 7.3 8.2 ...
##  $ INTG: num  7.9 8.9 8.1 8.8 6.4 8.8 9 5.9 8.9 7.9 ...
##  $ DMNR: num  7.7 8.8 7.8 8.5 4.3 8.7 8.9 4.9 8.9 6.7 ...
##  $ DILG: num  7.3 8.5 7.8 8.8 6.5 8.5 8.7 5.1 8.7 8.1 ...
##  $ CFMG: num  7.1 7.8 7.5 8.3 6 7.9 8.5 5.4 8.6 7.9 ...
##  $ DECI: num  7.4 8.1 7.6 8.5 6.2 8 8.5 5.9 8.5 8 ...
##  $ PREP: num  7.1 8 7.5 8.7 5.7 8.1 8.5 4.8 8.4 7.9 ...
##  $ FAMI: num  7.1 8 7.5 8.7 5.7 8 8.5 5.1 8.4 8.1 ...
##  $ ORAL: num  7.1 7.8 7.3 8.4 5.1 8 8.6 4.7 8.4 7.7 ...
##  $ WRIT: num  7 7.9 7.4 8.5 5.3 8 8.4 4.9 8.5 7.8 ...
##  $ PHYS: num  8.3 8.5 7.9 8.8 5.5 8.6 9.1 6.8 8.8 8.5 ...
##  $ RTEN: num  7.8 8.7 7.8 8.7 4.8 8.6 9 5 8.8 7.9 ...

image.png

14.2.1 判断主成分的个数

用来判断PCA中需要多少个主成分的准则：
根据先验经验和理论知识判断主成分数；
根据要解释变量方差的积累值的阈值来判断需要的主成分数；
通过检查变量间k × k的相关系数矩阵来判断保留的主成分数。
最常见的是基于特征值的方法。每个主成分都与相关系数矩阵的特征值相关联，第一主成分与最大的特征值相关联，第二主成分与第二大的特征值相关联，依此类推。
Kaiser-Harris准则建议保留特征值大于1的主成分，特征值小于1的成分所解释的方差比包含在单个变量中的方差更少。Cattell碎石检验则绘制了特征值与主成分数的图形。这类图形可以清晰地展示图形弯曲状况，在图形变化最大处之上的主成分都可保留。最后，你还可以进行模拟，依据与初始矩阵相同大小的随机数据矩阵来判断要提取的特征值。若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值，那么该主成分可以保留。该方法称作平行分析。

library(psych) # 调用psych包。
fa.parallel(USJudgeRatings[,-1], fa = "pc", n.iter = 100, show.legend = FALSE, main = "Scree plot with parallel analysis") # 碎石图判断主成分个数。
abline(h=1,lwd=1,col="green") # 添加特征值准则线。

image.png

图形解读：线段和x符号组成的图（蓝色线）：特征值曲线；
红色虚线：根据100个随机数据矩阵推导出来的平均特征值曲线；
绿色实线：特征值准则线（即：y=1的水平线）
判别标准：特征值大于平均特征值，且大于y=1的特征值准则线，被认为是可保留的主成分。根据判别标准，保留1个主成分即可。

fa.parallel函数学习
fa.parallel(data,n.obs=,fa=”pc”/”both”,n.iter=100,show.legend=T/F)
data：原始数据数据框；
n.obs：当data是相关系数矩阵时，给出原始数据（非原始变量）个数，data是原始数据矩阵时忽略此参数；
fa：“pc”为仅计算主成分，“fa”为因子分析，“both”为计算主成分及因子；
n.iter：模拟平行分析次数；
show.legend：显示图例。

14.2.2 提取主成分

principal(r, nfactors = , rotate = , scores = )

r：相关系数矩阵或原始数据矩阵；
nfactors：设定主成分数（默认为1）；
rotate：指定旋转的方法，默认最大方差旋转（varimax）。
scores：设定是否需要计算主成分得分（默认不需要）。

library(psych) # 调用psych包。
pc <- principal(USJudgeRatings[,-1], nfactors = 1) # 提取1个主成分。
pc # 返回结果。
## Principal Components Analysis
## Call: principal(r = USJudgeRatings[, -1], nfactors = 1)
## Standardized loadings (pattern matrix) based upon 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 is sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.04 
##  with the empirical chi square  6.21  with prob <  1 
## 
## Fit based upon off diagonal values = 1

PC1栏包含了成分载荷，指观测变量与主成分的相关系数。如果提取不止一个主成分，那么还将会有PC2、PC3等栏。成分载荷（component loadings）可用来解释主成分的含义，解释主成分与各变量的相关程度。
h2栏为成分公因子方差，即主成分对每个变量的方差解释度。
u2栏为成分唯一性，即方差无法被主成分解释的部分（1-h2）。
SS loadings包含了与主成分相关联的特征值，其含义是与特定主成分相关联的标准化后的方差值，即可以通过它来看90%的方差可以被多少个成分解释，从而选出主成分（即可使用nfactors=原始变量个数来把所有特征值查出，当然也可以直接通过eigen函数对它的相关矩阵进行查特征值）。
Proportion Var表示每个主成分对整个数据集的解释程度。
Cumulative Var表示各主成分解释程度之和。
Proportion Explained及Cumulative Proportion分别为按现有总解释方差百分比划分主成分及其累积百分比。

结果解读：第一主成分（PC1）与每个变量都高度相关，也就是说，它是一个可用来进行一般性评价的维度。ORAL变量99.1%的方差都可以被PC1来解释，仅仅有0.91%的方差不能被PC1解释。第一主成分解释了11个变量92%的方差。

head(Harman23.cor) # 查看数据集Harman23.cor 
## $cov
##                height arm.span forearm lower.leg weight bitro.diameter
## height          1.000    0.846   0.805     0.859  0.473          0.398
## arm.span        0.846    1.000   0.881     0.826  0.376          0.326
## forearm         0.805    0.881   1.000     0.801  0.380          0.319
## lower.leg       0.859    0.826   0.801     1.000  0.436          0.329
## weight          0.473    0.376   0.380     0.436  1.000          0.762
## bitro.diameter  0.398    0.326   0.319     0.329  0.762          1.000
## chest.girth     0.301    0.277   0.237     0.327  0.730          0.583
## chest.width     0.382    0.415   0.345     0.365  0.629          0.577
##                chest.girth chest.width
## height               0.301       0.382
## arm.span             0.277       0.415
## forearm              0.237       0.345
## lower.leg            0.327       0.365
## weight               0.730       0.629
## bitro.diameter       0.583       0.577
## chest.girth          1.000       0.539
## chest.width          0.539       1.000
## 
## $center
## [1] 0 0 0 0 0 0 0 0
## 
## $n.obs
## [1] 305
library(psych) # 调用psych包。
fa.parallel(Harman23.cor$cov, n.obs=302, fa="pc", n.iter=100, show.legend=FALSE, main="Scree plot with parallel analysis") # 判定主成分数量。

image.png

结果解读：通过碎石图可以判定选择的主成分个数为2个。

library(psych) # 调用psych包。
pc1 <- principal(Harman23.cor$cov, nfactors=2, rotate="none") # 提取2个主成分。
pc1 # 返回提取主成分的结果。
## Principal Components Analysis
## Call: principal(r = Harman23.cor$cov, nfactors = 2, rotate = "none")
## Standardized loadings (pattern matrix) based upon 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 are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.05 
## 
## Fit based upon off diagonal values = 0.99

结果解读：从结果Proportion Var： 0.58和0.22可以判定，第一主成分解释了身体测量指标58%的方差，而第二主成分解释了22%，两者总共解释了81%的方差。对于高度变量，两者则共解释了其88%的方差。

14.2.3 主成分旋转

旋转是一系列将成分载荷阵变得更容易解释的数学方法，它们尽可能地对成分去噪。旋转方法有两种：使选择的成分保持不相关（正交旋转），和让它们变得相关（斜交旋转）。旋转方法也会依据去噪定义的不同而不同。最流行的正交旋转是方差极大旋转，它试图对载荷阵的列进行去噪，使得每个成分只是由一组有限的变量来解释（即载荷阵每列只有少数几个很大的载荷，其他都是很小的载荷）。结果列表中列的名字都从PC变成了RC，以表示成分被旋转。

library(psych) # 调用psych包。
rc <- principal(Harman23.cor$cov, nfactors=2, rotate="varimax") # 主成分数判定，采用旋转。
rc # 返回结果。
## Principal Components Analysis
## Call: principal(r = Harman23.cor$cov, nfactors = 2, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                 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 are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.05 
## 
## Fit based upon off diagonal values = 0.99

14.2.4 获取主成分得分

当scores = TRUE时，主成分得分存储在principal()函数返回对象的scores元素中。

library(psych) # 调用psych包。
pc2 <-principal(USJudgeRatings[,-1], nfactors=1, score=TRUE) # 获取成分得分。
head(pc2$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

rc <- principal(Harman23.cor$cov, nfactors = 2, rotate = "varimax") # 获取主成分得分的系数。
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

14.3 探索性因子分析

如果你的目标是寻求可解释观测变量的潜在隐含变量，可使用因子分析。
EFA的目标是通过发掘隐藏在数据下的一组较少的、更为基本的无法观测的变量，来解释一
组可观测变量的相关性。这些虚拟的、无法观测的变量称作因子。（每个因子被认为可解释多个
观测变量间共有的方差，因此准确来说，它们应该称作公共因子。）
$X_i=a_1F_1+a_2F_2+a_3F_3....+a_pF_p+U_i$

其中 $X_i$ 是第i个可观测变量（i = 1…k）， $F_j$ 是公共因子（j = 1…p），并且p<k。 $U_i$ 是 $X_i$ 变量独有的部分（无法被公共因子解释）。 $a_i$ 可认为是每个因子对复合而成的可观测变量的贡献值。

ability.cov # 查看数据集。
## $cov
##         general picture  blocks   maze reading   vocab
## general  24.641   5.991  33.520  6.023  20.755  29.701
## picture   5.991   6.700  18.137  1.782   4.936   7.204
## blocks   33.520  18.137 149.831 19.424  31.430  50.753
## maze      6.023   1.782  19.424 12.711   4.757   9.075
## reading  20.755   4.936  31.430  4.757  52.604  66.762
## vocab    29.701   7.204  50.753  9.075  66.762 135.292
## 
## $center
## [1] 0 0 0 0 0 0
## 
## $n.obs
## [1] 112
options(digits = 2) # 设置数值显示的小数点位数。
covariances <- ability.cov$cov # 提取协方差矩阵的cov。
correlations <- cov2cor(covariances) # 将协方差矩阵转化为相关系数矩阵。
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

14.3.1 判断需提取的公共因子数

library(psych) # 调用psych包。
covariances <- ability.cov$cov # 提取协方差矩阵的cov。
correlations <- cov2cor(covariances) # 将协方差矩阵转化为相关系数矩阵。
fa.parallel(correlations,n.obs = 112,fa = "both",n.iter = 100, main = "Scree plots with paralled analysis") # 判断要提取的因子数。

image.png

碎石检验的前两个特征值（三角形）都在拐角处之上，并且大于基于100次模拟数据矩阵的特征值均值。对于EFA，Kaiser-Harris准则的特征值数大于0，而不是1。
结果解读：PCA结果建议提取一个或者两个成分，EFA建议提取两个因子。

14.3.2 提取公共因子

fa(r, nfactors=, n.obs=, rotate=, scores=, fm=)
 r是相关系数矩阵或者原始数据矩阵；
 nfactors设定提取的因子数（默认为1）；
 n.obs是观测数（输入相关系数矩阵时需要填写）；
 rotate设定旋转的方法（默认互变异数最小法）；
 scores设定是否计算因子得分（默认不计算）；
 fm设定因子化方法（默认极小残差法）。
与PCA不同，提取公共因子的方法很多，包括最大似然法（ml）、主轴迭代法（pa）、加权最小二乘法（wls）、广义加权最小二乘法（gls）和最小残差法（minres）。统计学家青睐使用最大似然法，因为它有良好的统计性质。

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) based upon 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 are sufficient.
## 
## The degrees of freedom for the null model are  15  and the objective function was  2.5
## The degrees of freedom for the model are 4  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 the residuals is  0.06 
## 
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy             
##                                                    PA1  PA2
## Correlation of (regression) scores with factors   0.96 0.92
## Multiple R square of scores with factors          0.93 0.84
## Minimum correlation of possible factor scores     0.86 0.68

结果解读：两个因子的Proportion Var分别为0.46和0.14，两个因子解释了六个心理学测试60%的方差。

14.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) based upon 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
## Test of the hypothesis that 2 factors are sufficient.
## 
## The degrees of freedom for the null model are  15  and the objective function was  2.5
## The degrees of freedom for the model are 4  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 the residuals is  0.06 
## 
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy             
##                                                    PA1  PA2
## Correlation of (regression) scores with factors   0.96 0.92
## Multiple R square of scores with factors          0.91 0.85
## Minimum correlation of possible factor scores     0.82 0.71

结果解读：阅读和词汇在第一因子上载荷较大，画图、积木图案和迷宫在第二因子上载荷较大，非语言的普通智力测量在两个因子上载荷较为平均，这表明存在一个语言智力因子和一个非语言智力因子。

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")
## Standardized loadings (pattern matrix) based upon 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
## 
##  With factor correlations of 
##      PA1  PA2
## PA1 1.00 0.55
## PA2 0.55 1.00
## 
## Mean item complexity =  1.2
## Test of the hypothesis that 2 factors are sufficient.
## 
## The degrees of freedom for the null model are  15  and the objective function was  2.5
## The degrees of freedom for the model are 4  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 the residuals is  0.06 
## 
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy             
##                                                    PA1  PA2
## Correlation of (regression) scores with factors   0.97 0.94
## Multiple R square of scores with factors          0.93 0.88
## Minimum correlation of possible factor scores     0.86 0.77

正交旋转和斜交旋转的不同之处。
对于正交旋转，因子分析的重点在于因子结构矩阵（变量与因子的相关系数），而对于斜交旋转，因子分析会考虑三个矩阵：因子结构矩阵、因子模式矩阵和因子关联矩阵。
因子模式矩阵即标准化的回归系数矩阵。它列出了因子预测变量的权重。因子关联矩阵即因子相关系数矩阵。

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。

fsm(fa.promax) # 通过fsm函数获取变量和因子间的相关系数。
##          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.promax, labels=rownames(fa.promax$loadings)) # 绘制正交或者斜交结果的图形。

image.png

图形解读：词汇和阅读在第一个因子（PA1）上载荷较大，而积木图案、画图和迷宫在第二个因子（PA2）上载荷较大。普通智力测验在两个因子上较为平均。

fa.diagram(fa.promax, simple=FALSE) # 两因子斜交旋转结果图。

image.png

14.3.4 因子得分

fa.promax$weights # 通过二因子斜交旋转法获得用来计算因子得分的权重。
##           PA1   PA2
## general 0.078 0.211
## picture 0.020 0.090
## blocks  0.037 0.702
## maze    0.027 0.035
## reading 0.743 0.030
## vocab   0.177 0.036

与可精确计算的主成分得分不同，因子得分只是估计得到的。它的估计方法有多种，fa()函数使用的是回归方法。

14.3.5 其他与 EFA 相关的包

R包含了其他许多对因子分析非常有用的软件包。FactoMineR包不仅提供了PCA和EFA方法，还包含潜变量模型。它有许多此处我们并没考虑的参数选项，比如数值型变量和类别型变量的使用方法。FAiR包使用遗传算法来估计因子分析模型，它增强了模型参数估计能力，能够处理不等式的约束条件，GPArotation包则提供了许多因子旋转方法。最后，还有nFactors包，它提供了用来判断因子数目的许多复杂方法。

14.4 其他潜变量模型

14.5 小结

image.png

实战练习

主成分分析

1.数据导入
数据结构：对10株玉米进行了生物学性状考察，考察指标有株高，穗位，茎粗，穗长，秃顶，穗粗，穗行数，行粒数。

df20 <- read.table(file = "D:/Documents/R wd/df20.csv", header = T, sep = ",") # 数据导入。
df20 # 查看数据。
##    株号 株高.cm. 穗位.cm. 茎粗.cm. 穗长.cm. 秃顶.cm. 穗粗.cm. 穗行数.行.
## 1     1      237       90       14       21      2.5       52         18
## 2     2      233       88       16       19      3.0       44         18
## 3     3      229       84       15       20      2.0       50         12
## 4     4      245       80       17       20      0.0       47         16
## 5     5      230       70       14       15      4.0       46         16
## 6     6      215       70       13       18      4.0       46         16
## 7     7      210       84       11       18      6.0       48         14
## 8     8      208       85       12       18      3.0       47         14
## 9     9      229       90       15       18      2.5       46         14
## 10   10      232       87       17       18      2.0       47         16
##    行粒数.粒.
## 1          34
## 2          35
## 3          37
## 4          36
## 5          35
## 6          23
## 7          23
## 8          24
## 9          32
## 10         31

判断主成分数量

library(psych) # 调用psych包。
df20.cor <- cor(df20[,-1]) # 计算相关矩阵。
fa.parallel(df20[,-1], fa = "pc", n.iter = 100, show.legend = FALSE, main = "Scree plot with parallel analysis") # 碎石检验判断主成分个数。
abline(h=1,lty=1,lwd=2,col="green") # 添加特征值准则线。

image.png

结果解读：选择2个主成分即可保留样本大量信息。

3.提取主成分

df20.df <- principal(df20[,-1], nfactors = 2, score = T, rotate = "varimax") # 提取2个主成分。
df20.df # 返回结果。
## Principal Components Analysis
## Call: principal(r = df20[, -1], nfactors = 2, rotate = "varimax", scores = T)
## Standardized loadings (pattern matrix) based upon correlation matrix
##              RC1   RC2   h2    u2 com
## 株高.cm.    0.95  0.17 0.93 0.066 1.1
## 穗位.cm.    0.12  0.72 0.53 0.471 1.1
## 茎粗.cm.    0.96  0.01 0.93 0.073 1.0
## 穗长.cm.    0.27  0.84 0.79 0.213 1.2
## 秃顶.cm.   -0.81 -0.32 0.76 0.242 1.3
## 穗粗.cm.   -0.14  0.80 0.65 0.349 1.1
## 穗行数.行.  0.46 -0.18 0.24 0.758 1.3
## 行粒数.粒.  0.85  0.19 0.75 0.248 1.1
## 
##                        RC1  RC2
## SS loadings           3.51 2.07
## Proportion Var        0.44 0.26
## Cumulative Var        0.44 0.70
## Proportion Explained  0.63 0.37
## Cumulative Proportion 0.63 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.1 
##  with the empirical chi square  5.9  with prob <  0.95 
## 
## Fit based upon off diagonal values = 0.95

结果解读：主成分1可解释44%的方差，主成分2解释了26%的方差，合计解释了70%的方差。

4.获取主成分得分

dft <- round(unclass(df20.df$weights),2) # 获取主成分得分。
dft # 返回结果。
##              RC1   RC2
## 株高.cm.    0.27 -0.01
## 穗位.cm.   -0.04  0.36
## 茎粗.cm.    0.29 -0.10
## 穗长.cm.   -0.01  0.41
## 秃顶.cm.   -0.21 -0.08
## 穗粗.cm.   -0.13  0.43
## 穗行数.行.  0.16 -0.15
## 行粒数.粒.  0.24  0.01

5.主成分方程

PC1 = 0.27 $\times$ 株高 - 0.04 $\times$ 穗位 + 0.29 $\times$ 茎粗 - 0.01 $\times$ 穗长 - 0.21 $\times$ 秃顶 - 0.13 $\times$ 穗粗 + 0.16 $\times$ 穗行数 + 0.24 $\times$ 行粒数

PC2 = -0.01 $\times$ 株高 + 0.36 $\times$ 穗位 - 0.10 $\times$ 茎粗 + 0.41 $\times$ 穗长 - 0.08 $\times$ 秃顶 + 0.43 $\times$ 穗粗 - 0.15 $\times$ 穗行数 + 0.01 $\times$ 行粒数

plot(df20.df) # 主成分分析可视化

image.png

图形解读：此图反映了变量与主成分的关系，三个蓝点对应的RC2值较高，点上的标号2，4，6对应变量名穗位，穗长，穗粗，说明第2主成分主要解释了这些变量，与这些变量相关性强；黑点分别对应株高，茎粗，穗行数，行粒数，说明第一主成分与这些变量相关性强，第一主成分主要解释的也是这些变量，而5号点秃顶对于两个主成分均没有显示好的相关性。

因子分析

判断因子数量

library(psych) # 调用psych包。
df20.cor <- cor(df20[,-1]) # 计算相关矩阵。
fa.parallel(df20[,-1], fa = "fa", n.iter = 100, show.legend = FALSE, main = "Scree plot with parallel analysis") # 碎石检验判断主成分个数。
abline(h=0,lty=1,lwd=2,col="green") # 添加特征值准则线。

image.png

图解：可以看到需要提取4个因子。

2.提取因子

df20.df1 <- fa(df20.cor, nfactors = 4, rotate = "promax",fm="ml") # 最大似然法提取4个因子。
df20.df1 # 返回结果。
## Factor Analysis using method =  ml
## Call: fa(r = df20.cor, nfactors = 4, rotate = "promax", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##              ML1   ML2   ML3   ML4   h2     u2 com
## 株高.cm.    0.97 -0.06  0.19  0.14 1.00 0.0050 1.1
## 穗位.cm.   -0.04  0.60 -0.02  0.02 0.35 0.6467 1.0
## 茎粗.cm.    0.87  0.09  0.03 -0.32 1.00 0.0049 1.3
## 穗长.cm.   -0.01  0.88  0.12  0.31 0.98 0.0169 1.3
## 秃顶.cm.   -0.68 -0.39  0.14  0.19 0.85 0.1493 1.9
## 穗粗.cm.    0.04  0.18 -0.06  0.73 0.63 0.3709 1.1
## 穗行数.行.  0.02  0.05  0.97 -0.09 0.97 0.0327 1.0
## 行粒数.粒.  1.03 -0.23 -0.07  0.26 0.87 0.1321 1.2
## 
##                        ML1  ML2  ML3  ML4
## SS loadings           3.24 1.45 1.04 0.92
## Proportion Var        0.41 0.18 0.13 0.11
## Cumulative Var        0.41 0.59 0.72 0.83
## Proportion Explained  0.49 0.22 0.16 0.14
## Cumulative Proportion 0.49 0.71 0.86 1.00
## 
##  With factor correlations of 
##       ML1   ML2   ML3   ML4
## ML1  1.00  0.43  0.24 -0.09
## ML2  0.43  1.00 -0.01  0.21
## ML3  0.24 -0.01  1.00 -0.04
## ML4 -0.09  0.21 -0.04  1.00
## 
## Mean item complexity =  1.3
## Test of the hypothesis that 4 factors are sufficient.
## 
## The degrees of freedom for the null model are  28  and the objective function was  7.5
## The degrees of freedom for the model are 2  and the objective function was  0.34 
## 
## The root mean square of the residuals (RMSR) is  0.03 
## The df corrected root mean square of the residuals is  0.11 
## 
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1  ML2  ML3  ML4
## Correlation of (regression) scores with factors   1.00 0.99 0.98 0.98
## Multiple R square of scores with factors          1.00 0.98 0.97 0.97
## Minimum correlation of possible factor scores     0.99 0.96 0.94 0.93

结果解读：因子1到4解释了80%的方差。

3.获取因子得分

df20.df1$weights # 返回得分。
##                 ML1     ML2     ML3     ML4
## 株高.cm.    0.63292 -0.5131  0.2709  1.4293
## 穗位.cm.   -0.00017  0.0188 -0.0030 -0.0029
## 茎粗.cm.    0.38640  0.6016 -0.2073 -1.6872
## 穗长.cm.    0.01902  0.9569 -0.0025  0.3491
## 秃顶.cm.   -0.01549 -0.0600  0.0467  0.0328
## 穗粗.cm.    0.00723  0.0018 -0.0174  0.0642
## 穗行数.行. -0.12547 -0.0308  0.9377 -0.2525
## 行粒数.粒.  0.03649 -0.0440 -0.0546  0.0999

可视化

factor.plot(df20.df1,labels=rownames(df20.df1$loadings)) # 可视化。
fa.diagram(df20.df1,simple = TRUE) # 因子分析可视化

image.png

图解：可以看出，因子1和因子2的相关系数为0.4，行粒数，株高，茎粗，秃顶在因子1的载荷较大，穗长，穗位在因子2上的载荷较大；因子3只有穗行数相关，因子4只有穗粗相关。

参考资料：

《R语言实战》（中文版），人民邮电出版社，2013.
如何理解主成分分析法（PCA），https://zhuanlan.zhihu.com/p/170398464
主成分分析法，https://blog.csdn.net/weixin_43914260/article/details/99585202

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 214,951评论 6赞 497
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 91,606评论 3赞 389
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 160,601评论 0赞 350
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 57,478评论 1赞 288
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 66,565评论 6赞 386
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 50,587评论 1赞 293
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,590评论 3赞 414
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,337评论 0赞 270
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,785评论 1赞 307
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,096评论 2赞 330
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,273评论 1赞 344
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,935评论 5赞 339
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,578评论 3赞 322
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,199评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,440评论 1赞 268
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,163评论 2赞 366
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,133评论 2赞 352

《R语言实战》自学笔记71-主成分和因子分析

第14章 主成分和因子分析

14.1 R中的主成分和因子分析

14.2 主成分分析

14.2.1 判断主成分的个数

14.2.2 提取主成分

14.2.3 主成分旋转

14.2.4 获取主成分得分

14.3 探索性因子分析

14.3.1 判断需提取的公共因子数

14.3.2 提取公共因子

14.3.3 因子旋转

14.3.4 因子得分

14.3.5 其他与 EFA 相关的包

14.4 其他潜变量模型

14.5 小结

实战练习

推荐阅读更多精彩内容

第14章主成分和因子分析