因子分析（Factor analysis, FA） - ricket-sjtu/bi028 GitHub Wiki

1. 什么是因子分析（Factor analysis）？

因子分析（FA）就是用少量的不可观测的潜在“因子”，对观察变量及其协方差结构进行建模的一种方法。

1.1 FA与PCA

乍看之下，FA与PCA颇为相似。然而，实际上两者是相反的：

  • PCA的PC是观察变量的线性函数；
  • FA的观察变量是Factor的线性组合

1.2 问题

  1. 什么是因子载荷（factor loading）？特异性方差（specific variances）？共同度（communalities）？
  2. 如何应用主成份和最大似然估计方法估计因子模型的参数？
  3. 因子旋转（factor rotation）？旋转后因子载荷（rotated factor loading）？

2. 模型

2.1 回归模型（一般要求$m << p$）：

$$ \begin{array}{llc} X_1 & = & \mu_1 + l_{11}f_1 + l_{12}f_2 + \dots + l_{1m}f_m + \epsilon_1\ X_2 & = & \mu_2 + l_{21}f_1 + l_{22}f_2 + \dots + l_{2m}f_m + \epsilon_2 \ & & \vdots \ X_p & = & \mu_p + l_{p1}f_1 + l_{p2}f_2 + \dots + l_{pm}f_m + \epsilon_p \end{array} $$

2.2 因子载荷：

$$ \mathbf{L} = \left(\begin{array}{cccc} l_{11}& l_{12}& \dots & l_{1m}\ l_{21} & l_{22} & \dots & l_{2m}\ \vdots & \vdots & & \vdots \ l_{p1} & l_{p2} & \dots & l_{pm} \end{array}\right) = \text{matrix of factor loadings} $$

2.3 特异性因子：

$$ \mathbf{\epsilon} = \left(\begin{array}{c} \epsilon_1\ \epsilon_2\ \vdots\ \epsilon_p \end{array}\right) = \text{vector of specific factors} $$

2.4 表达式的矩阵形式：

$$ \textbf{X} = \mu + \textbf{Lf}+\mathbf{\epsilon} $$

3. 模型有哪些假设？

3.1 针对均值的假设

  1. 特异性因子的均值为0：$E(\epsilon_i) = 0$
  2. 共同因子的均值也为0：$E(f_i) = 0$

3.2 针对方差的假设

  1. 共同因子的方差为1：$\text{var}(f_i) = 1$
  2. 特异性因子的方差：$\text{var}(\epsilon_i) = \psi_i$，其中$\psi_i$称为特异性方差（specific variance）

3.3 针对相关性的假设

  1. 共同因子相互不相关：$\text{cov}(f_i, f_j) = 0$
  2. 不同特异性因子互不相关：$\text{cov}(\epsilon_i, \epsilon_j) = 0$
  3. 特异性因子和共同因子互不相关：$\text{cov}(f_i, \epsilon_j) = 0$

由此，我们可以得到

  • 第$i$个特征（观察变量）的方差为：$\sigma^2_i = \text{var}(X_i) = \sum_{j=1}^{m}l^2_{ij}+\psi_i$；
  • 其中$\sum_{j=1}^{m}l^2_{ij}$称为第$i$个特征的共性系数（communality）；
  • 第$i$个和第$j$个特征的协方差为：$\sigma_{ij}= \text{cov}(X_i, X_j) = \sum_{k=1}^{m}l_{ik}l_{jk}$；
  • 特征$i$和因子$j$的协方差为：$\text{cov}(X_i, f_j) = l_{ij}$；
  • 模型的方差-协方差矩阵为：$\Sigma = \mathbf{LL'} + \Psi$，其中 $$ \mathbf{\Psi} = \left(\begin{array}{cccc} \psi_1 & 0 & \dots & 0 \ 0 & \psi_2 & \dots & 0 \ \vdots & \vdots & \ddots & \vdots\ 0 & 0 & \dots & \psi_p \end{array}\right) $$

3.4 因子分析前的假设检验

在对数据进行因子分析前，必须进行KMO检验和Bartlett检验：

  • KMO检验统计值大于0.5；
  • Bartlett检验的P值小于0.05 才能进行因子分析。

KMO检验

KMO（Kaiser-Meyer-Olkin)检验统计量是用于比较变量间简单相关系数（simple correlation）和偏相关系数（partial correlation）的指标:

  • 当所有变量间的简单相关系数平方和远远大于偏相关系数平方和时，KMO值接近1，意味着变量间的相关性越强，原有变量越适合进行因子分析；
  • 当所有变量间的简单相关系数平方和接近0时，KMO值接近0，意味着变量间的相关性越弱，原有变量越不适合作因子分析。

KMO > 0.9：非常适合；KMO>0.8：适合；KMO>0.7：一般；KMO>0.6：不太适合；KMO<0.5：极不适合。

library(MASS)
kmo <- function( data ){
  X <- cor(as.matrix(data))
  iX <- ginv(X)
  S2 <- diag(diag((iX^-1)))
  AIS <- S2%*%iX%*%S2                      # anti-image covariance matrix
  IS <- X+AIS-2*S2                         # image covariance matrix
  Dai <- sqrt(diag(diag(AIS)))
  IR <- ginv(Dai)%*%IS%*%ginv(Dai)         # image correlation matrix
  AIR <- ginv(Dai)%*%AIS%*%ginv(Dai)       # anti-image correlation matrix
  a <- apply((AIR - diag(diag(AIR)))^2, 2, sum)
  AA <- sum(a)
  b <- apply((X - diag(nrow(X)))^2, 2, sum)
  BB <- sum(b)
  MSA <- b/(b+a)                        # indiv. measures of sampling adequacy

  AIR <- AIR-diag(nrow(AIR))+diag(MSA)  # Examine the anti-image of the
                                        # correlation matrix. That is the
                                        # negative of the partial correlations,
                                        # partialling out all other variables.

  kmo <- BB/(AA+BB)                     # overall KMO statistic

  # Reporting the conclusion
    if (kmo >= 0.00 && kmo < 0.50){
      test <- 'The KMO test yields a degree of common variance
unacceptable for FA.'
    } else if (kmo >= 0.50 && kmo < 0.60){
      test <- 'The KMO test yields a degree of common variance miserable.'
    } else if (kmo >= 0.60 && kmo < 0.70){
      test <- 'The KMO test yields a degree of common variance mediocre.'
    } else if (kmo >= 0.70 && kmo < 0.80){
      test <- 'The KMO test yields a degree of common variance middling.'
    } else if (kmo >= 0.80 && kmo < 0.90){
      test <- 'The KMO test yields a degree of common variance meritorious.'
    } else {
      test <- 'The KMO test yields a degree of common variance marvelous.'
    }

    ans <- list(  overall = kmo,
                  report = test,
                  individual = MSA,
                  AIS = AIS,
                  AIR = AIR )
    return(ans)

}    # end of kmo()

Bartlett检验

Bartlett检验的统计值是通过计算相关系数矩阵的行列式（determinant）得到的。如果统计值很大，则对应$p-val$小于0.05，则认为该矩阵不可能是单位矩阵，也就是说变量之间存在相关性，适合进行因子分析；反之，则表明该矩阵接近单位矩阵，变量之间不存在相关性，不适合进行因子分析。

psych::cortest.bartlett(data)

4. 因子分析方法

如何确定因子数目

确定因子数目必须遵循下列几个原则：

  1. 累积比例（累积解释的方差）至少为80%；
  2. 最小的特征值不小于1；
  3. 碎石图（scree plot）中的拐点

4.1 主成份法（principal component method）

  1. 输入矩阵$\mathbf{X} \in \mathbb{R}^{n \times p}$，其中每个分量$X_i = \left(\begin{array}{c}X_{i1}\ X_{i2}\ \vdots \ X_{ip}\end{array}\right), i=1,2,\ldots, n$
  2. 定义方差-协方差矩阵：$\textbf{S} = \frac{1}{n-1}\sum_{i=1}^{n}\mathbf{(X_i - \bar{x})(X_i - \bar{x})^T}$
  3. 计算方差-协方差矩阵的特征值（$\hat{\lambda}_1, \hat{\lambda}_2, \dots, \hat{\lambda}_p$）和特征向量（$\hat{\mathbf{e}}_1, \hat{\mathbf{e}}_2, \dots, \hat{\mathbf{e}}_p$）
  4. 方差-协方差矩阵的谱分解（spectral decomposition） $$ \Sigma = \sum_{i=1}^{p}\lambda_i \mathbf{e_ie_i^T} \cong \sum_{i=1}^{m}\lambda_i \mathbf{e_ie_i^T} = \left(\begin{array}{cccc}\sqrt{\lambda_1}\mathbf{e_1} & \sqrt{\lambda_2}\mathbf{e_2} & \dots & \sqrt{\lambda_m}\mathbf{e_m}\end{array}\right) \left(\begin{array}{c}\sqrt{\lambda_1}\mathbf{e_1^T}\ \sqrt{\lambda_2}\mathbf{e_2^T}\ \vdots\ \sqrt{\lambda_m}\mathbf{e_m^T}\end{array}\right) = \mathbf{LL^T} $$ 也就是将特征分解从$p$维压缩到$m$维，其余$p-m$维被忽略，这里$m$为因子个数
  5. 因子载荷（factor loading）就可以这样估计：$\hat{l}{ij} = \hat{e}{ji}\sqrt{\hat{\lambda}_j}$
  6. 相应的，特异性方差这样计算：$\mathbf{\Psi} = \Sigma - \mathbf{LL^T}$，其中$\mathbf{\Psi}$的对角线元素就代表了特异性方差： $$ \hat{\Psi}i = s_i^2 - \sum{j=1}^m \hat{l}{ij}^2 = s_i^2 - \sum{j=1}^m \lambda_i \hat{e}_{ji}^2 $$
  7. 其中$\hat{h}i^2 = \sum{j=1}^m \hat{l}_{ij}^2$定义为第$i$个变量的共同度（communality，数值上相当于多重线性回归的$R^2$，也就是变量的方差有多大比例能被几个因子所解释）

总结

  1. 因子载荷矩阵（factor loading matrix）的每一行的平方和是对应变量的共同度（communality），反映的是该变量对公共因子（common factor）的依赖程度，数值越大，表明变量对因子的依赖程度越高；相应的，特异性因子（特异性方差）也就越小
  2. 因子载荷矩阵（factor loading matrix）的每一列的平方和反映的是对应因子对变量的贡献率，值越大则贡献率越高

4.2 最大似然估计法（maximum likelihood estimation method）

4.2.1 MLE的假设

最大似然估计方法假设数据来自独立采样的样本，符合均值为$\mu$且方差-协方差矩阵符合下面形式的： $$ \Sigma = \mathbf{LL^T +\Psi} $$ 的多元正态分布（multivariate normal distribution），其中

  • $\mathbf{L}$是因子载荷矩阵（factor loading matrix）；
  • $\Psi$是特异性方差对角矩阵（diagonal matrix of specific variances）

4.2.2 MLE的目标

对于具有$n$个输入观察值$\mathbf{X_1},\mathbf{X_2}, \dots, \mathbf{X_n}$的输入矩阵来说，最大似然估计需要做的就是

  • 估计均值$\hat{\mu}$；
  • 估计因子载荷$\hat{\mathbf{L}}$；
  • 估计特异性方差矩阵$\hat{\Psi}$

估计是通过最大化对数似然函数（log-likelihood）实现： $$ \ell(\mathbf{\mu, L, \Psi}) = - \frac{np}{2}\log{2\pi}- \frac{n}{2}\log{|\mathbf{LL^T + \Psi}|} - \frac{1}{2}\sum_{i=1}^{n}\mathbf{(X_i-\mu)^T(LL^T+\Psi)(X_i-\mu)} $$

4.2.3 MLE的问题与解决方案

  • 上面的对数似然函数的优化目标函数没有闭合解，必须通过数值迭代方法进行求解；
  • 数值方法的解非唯一，例如可以通过旋转得到不同的解；
  • 为了保证解的唯一性，可以加入约束：$\mathbf{L'\Psi^{-1}L}$为对角矩阵；
  • 数值解不稳定

4.2.4 拟合优度检验（goodness of fit）

由于以上方法中，模型假设 $$ \mathbf{\Sigma = LL^T + \Psi} $$ 因此，拟合优度检验引入Bartlett校正的似然比检验统计值（LRT）： $$ X^2 = \left(n-1-\frac{2p+4m-5}{6}\right)\log \frac{|\mathbf{\hat{L}\hat{L}^T}+\mathbf{\hat{\Psi}}|}{|\hat{\mathbf{\Sigma}}|} $$

  • $n-1-\frac{2p+4m-5}{6}$称为Bartlett校正常数（Bartlett-corrected constant）；
  • 分子$|\mathbf{\hat{L}\hat{L}^T}+\mathbf{\hat{\Psi}}|$对应的是上面模型的似然值，这里$|\cdot|$是行列式；
  • 分母$|\hat{\mathbf{\Sigma}}|$对应的无模型的似然值$\hat{\mathbf{\Sigma}} = \frac{n-1}{n}\mathbf{S}$，$S$是样本方差-协方差矩阵；
  • 当零假设成立时，$\mathbf{X}^2 \sim \chi^2_{\frac{(p-m)^2-p-m}{2}}$
  • 如果模型拟合得很好，则$\mathbf{X}^2$值会很小；反之，$\mathbf{X}^2$值会很大

If $p<0.05$, we conclude that the relationships among the variables is not adequately described by the factor model. This suggests that we do not have the correct model.

The only real remedy that we can apply is to increase the number $m$, number of factors, until an adequate fit is achieved.

However, $m$ must satisfy $p(m+1) < p(p+1)/2$.

5. 因子旋转（factor rotation）

由于因子载荷矩阵$\mathbf{L}$反映的是变量与公共因子（common factors）之间的关系，因此出于模型理解（interpretation）的需要，我们希望在$\mathbf{L}$的每一行上都是稀疏的，但是情况往往并非如此。

由于因子模型$\mathbf{X = \mu + LF + \epsilon}$中，包含均值（mean vector）$\mu$，因子载荷矩阵（factor loading matrix）$\mathbf{L}$，公共因子（common factors）$\mathbf{F}$，以及特异性因子（specific factors）$\epsilon$。

对于这个模型，必然存在另一个旋转模型$\mathbf{X = \mu + LF + \epsilon}$，其中$\mathbf{L* = LT}, \mathbf{F* = T^T F}$，且$\mathbf{T}$为正交矩阵满足$\mathbf{T^T T = T T^T = I}$。因此，我们的目标就是确定一个正交矩阵$\mathbf{T}$，使产生的模型容易被理解。

5.1 Varimax因子旋转

Varimax旋转是最常用的因子旋转方法之一：

  1. 计算标准化载荷：$\tilde{l}^_{ij}= \hat{l}^_{ij}/\hat{h}_i$
  2. 选择旋转矩阵$\mathbf{T}$，使得下面的目标函数最大化： $$ V = \frac{1}{p}\sum_{j=1}^{m}\left{\sum_{i=1}^{p}(\tilde{l}^{ij})^4 - \frac{1}{p}\left(\sum{i=1}^{p}(\tilde{l}^_{ij})^2 \right)^2 \right} $$ 这其实就是所有因子的标准化载荷（standardized loading）的样本方差的期望值。

6. 因子分值（factor scores）的计算

因子分值的计算与PCA中主成份分值的计算类似。对于因子模型： $$ \mathbf{Y_i = \mu + Lf_i + \epsilon_i}; i = 1,2,\dots, n, $$ 计算因子分值，也就是为每个观察样本估计 $$ \mathbf{f_1, f_2, \dots, f_n} $$

6.1 计算方法

常用的计算方法主要有三种：

  • 最小二乘法（ordinary least squares）；
  • 加权最小二乘法（weighted least squares）；
  • 回归法（regression）

6.2 最小二乘法（OLS）

  • 最小化残差平方和： $$ \sum_{j=1}^{p}\epsilon^2_{ij} = \sum_{j=1}^{p}(y_{ij}-\mu_j-l_{j1}f_1 - l_{j2}f_2 - \dots - l_{jm}f_m)^2 = (\mathbf{Y_i - \mu - Lf_i})^T (\mathbf{Y_i - \mu - Lf_i}) $$
  • 最小二乘法的理论解为： $$ \mathbf{\hat{f}_i = (L^T L)^{-1}L^T(Y_i-\mu)} $$
  • 实际解为： $$ \mathbf{\hat{f}_i = \left(\hat{L}^T \hat{L}\right)^{-1}\hat{L}^T (Y_i-\bar{y})} = \left(\begin{array}{c} \frac{1}{\sqrt{\hat{\lambda}_1}}\mathbf{\hat{e}_1^T(Y_i-\bar{y})}\ \frac{1}{\sqrt{\hat{\lambda}_2}}\mathbf{\hat{e}_2^T(Y_i-\bar{y})}\ \vdots \ \frac{1}{\sqrt{\hat{\lambda}_m}}\mathbf{\hat{e}_m^T(Y_i-\bar{y})}\end{array}\right) $$ 这里$\hat{\lambda}_i}, \hat{e}_i, i=1,\dots,m$为前$m$个特征值和特征向量。

6.3 加权最小二乘法（WLS）

  • WLS与OLS的区别是在目标函数中加入了特异性方差的倒数作为权重： $$ \sum_{j=1}^{p}\frac{\epsilon^2_{ij}}{\Psi_j} = \sum_{j=1}^{p}\frac{(y_{ij}-\mu_j - l_{j1}f_1 - l_{j2}f_2 -\dots - l_{jm}f_m)^2}{\Psi} = \mathbf{(Y_i-\mu-Lf_i)^T \Psi^{-1}(Y_i-\mu-Lf_i)} $$
  • 理论解： $$ \mathbf{\hat{f}_i = (L^T \Psi^{-1}L)^{-1}L^T \Psi^{-1}(Y_i-\mu)} $$
  • 实际解： $$ \mathbf{\hat{f}_i = (\hat{L}^T \hat{\Psi}^{-1}\hat{L})^{-1}\hat{L}^T \hat{\Psi}^{-1}(Y_i-\bar{y})} $$

6.4 回归法（Regression method）

  • 如果因子分析的因子载荷矩阵是通过最大似然估计得到的，那么应该采用回归法计算因子得分；
  • 联合概率分布： $$ \left(\begin{array}{c} \mathbf{Y_i} \ \mathbf{f_i} \end{array}\right) \sim N \left[\left( \begin{array}{c} \mathbf{\mu} \ 0 \end{array}\right), \left( \begin{array}{cc} \mathbf{LL^T+\Psi} & \mathbf{L} \ \mathbf{L^T} & \mathbf{I} \end{array}\right)\right] $$
  • 给定数据$\mathbf{Y_i}$的公共因子得分$\mathbf{f}_i$的条件期望为： $$ E(\mathbf{f_i | Y_i}) = \mathbf{L^T(LL^T+\Psi)^{-1}(Y_i-\mu)} $$
  • 其估计值$\mathbf{\hat{f}_i$可写成： $$ \mathbf{\hat{f}_i = \hat{L}^T \left(\hat{L}\hat{L}^T + \hat{\Psi}\right)^{-1}(Y_i-\bar{y})} $$

7. 总结

  1. 因子载荷矩阵的意义是什么？
  2. 如何用主成份法或最大似然法估计因子载荷矩阵和特异性方差？
  3. 共同度（communality）是如何反映因子模型的充分性的？
  4. 如何用似然比检验（LRT）评价因子模型的拟合优度？
  5. 因子旋转能够增加因子模型的实际意义，常用的因子旋转方法为Varimax。
  6. 因子得分的计算方法有哪些？分别适用于哪些特殊场景？