【R语言】主成分分析与因子分析

一、主成分分析

主成分分析（Principal Component Analysis, PCA）是一种常用的无监督数据降维技术，广泛应用于统计学、数据科学和机器学习等领域。它通过正交化线性变换将（高维）原始数据投影到一个新的坐标系（低维空间），同时保留数据的主要变异信息，使得新坐标系中的第一个坐标轴（主成分1）上的方差最大，第二个坐标轴（主成分2）上的方差次大，并且各个坐标轴之间相互正交（即不相关）。PCA 的目标是通过保留数据的主要变异方向来减少数据的维度，同时尽可能少地丢失信息。

psych扩展包能提供很丰富和有用的函数，它的输出结果也更接近商业统计分析软件，如SAS和SPSS等，所使用的函数为principal()函数。

principal()函数

# 准备数据（以 iris 数据集为例）
head(iris)
# 去除标签列，只保留数值列
iris_data <- iris[, 1:4]
head(iris_data)

 principal()函数的nfactors用来指定各种主成分，iris_data里面有4种变量，所以指定为4。

library(psych)
pcal_iris <- principal(iris_data, nfactors=4)
pcal_iris

从上面的结果可知，前三个成分（RC1+RC3+RC2）的累计解释比例（Cumulative Proportion）为99%，前两个成分（RC1+RC3）的累计解释比例（Cumulative Proportion）为72%，所以选择这两个或三个主要成分是合适的。 

pcal_iris1 <- principal(iris_data, nfactors=3)
pcal_iris1

从上面结果可知，当nfactors为3时，累计比例能近似达到100%，说明用三个主成分代替原有变量的信息是足够的。 

二、因子分析

因子分析（Factor Analysis）是一种降维技术，它用于研究变量之间的内在关系，试图通过少数几个“潜在变量”（或称“因子”）来解释多个观测变量之间的相关性，其核心思想是降维和结构探测。这些潜在变量是不可直接观测的，但可以通过它们对观测变量的影响来推断。因子分析在心理学、社会学、市场研究、生物学等多个领域都有广泛应用。

根据是否已知潜在结构，可以将因子分析分为探索性因子分析（Exploratory Factor Analysis，EFA）和验证性因子分析（Confirmatory Factor Analysis，CFA）。

1、注意事项

在进行因子分析，需要注意以下两个方面。

保证有足够的样本量。一般认为样本量小于50时，不适合做因子分析；样本量至少在100以上；样本量达到1000时，效果会比较好。而且，样本量的选择还受原始变量数量的影响，一般样本量至少是原始变量数量的5倍以上，10倍更好。

原始变量之间应该要有足够的相关性。如果所有或者大部分原始变量是相互独立或者相关系数都小于0.3，则不能从中提取公共因子，即数据不适合进行因子分析。

原始变量之间的相关性可以使用KMO（Kaiser-Meyer-Oklin）检验或Bartlett球形检验。

因子分析可以使用psych扩展包中的fa()函数。

fa(r, nfactors, n.obs, rotate, scores, fm)

1. r：相关系数矩阵或者原始数据矩阵；
2. nfactor：因子数，默认为1；
3. n.obs：观测数，当r为相关系数矩阵时需手动输入；
4. rotate：设定因子旋转的方法，默认为promax（斜交旋转）；还有 varimax（正交旋转）旋转有助于更好地解释因子结构。
5. scores：是否计算因子得分，默认为FALSE，且要求r为原始数据矩阵；
6. fm：因子提取方法，默认为minres，一般选择ml（最大似然法）。

2、分析步骤

下面使用R语言中的内置数据集mtcars做演示：

2.1 准备数据

library(psych)
library(GPArotation) # 支持因子旋转
head(mtcars)
data <- mtcars[,c("mpg", "disp", "hp", "drat", "wt", "qsec")]
head(data)

2.2 检查数据的适用性

检验相关性：如果变量间相关性较低（绝对值<0.3），可能不适合因子分析

cor_data <- cor(data)
cor_data

从下面结果可知，绝大多数的绝对值都大于0.3。 

 KMO检验：KMO值>0.6表示数据适合做因子分析。

KMO(data)

从下面结果可知，Overall MSA=0.76>0.6。 

 Bartlett球形检验： p<0.05时适合做因子分析。

cortest.bartlett(cor_data, n = nrow(data)) 

从下面结果可知，p=1.332068e-30 < 0.05。 

 2.3 确定因子数量

特征值（Eigenvalues）：碎石图，选择特征值>1的因子

eigen_values <- eigen(cor_data)$values
plot(eigen_values, type = "b", main = "Scree Plot")

从以下结果可知，只能选择2个因子。 

 平行分析（Parallel Analysis）：

fa.parallel(cor_data, n.obs=nrow(data), fm="ml", fa = "fa", n.iter=100)

从以下碎石图结果可知，建议取2个因子（虚线上面的小三角形个数） 

2.4 进行因子分析

使用正交旋转。

#   数据为相关矩阵 
#   nfactors: 因子数
#   最大迭代次数为100次
#   rotate: 旋转方法（"varimax"正交旋转，"oblimin"斜交旋转）
#   fm: 因子提取方法（"pa"主成分，"ml"极大似然）
result <- fa(r = cor_data, nfactors = 2, n.obs = nrow(data), n.iter=100, rotate = "varimax", fm = "ml")
result

从以下结果可知，2个因子一共解释了100%的变异。 

使用斜交旋转。

result <- fa(r = cor_data, nfactors = 2, n.obs = nrow(data), n.iter=100, rotate = "promax", fm = "ml")
result

 2.5 可视化结果

# 绘制因子载荷图
fa.diagram(result)