因子分析——数学原理及R语言代码

因子分析

FactorAnalysis的目的是从多个高度相关的观测变量中提取出少数几个LatentFactor，这些因子代表了变量背后的共通结构，从而实现降维并提升可解释性。

假设对一组学生进行了以下六门课程的测试：语文、英语、数学、物理、化学、生物，发现语文和英语成绩之间高度相关，数学、物理、化学、生物也彼此高度相关。此时可以猜测：这些成绩可能是由两个更基本的“能力”决定的，比如语言能力和理科能力。通过因子分析就可以提取出这两个潜在因子，并发现语文和英语主要由“语言能力”因子决定，理科四门主要由“理科能力”因子解释。这样就可以用两个因子有效地概括了六个变量的结构，同时让模型更易解释、更简洁。

设$‘\mathbf{X}‘$是一个可观测的$‘m‘$维随机向量，$‘E(\mathbf{X})=\mathbf{\mu },\mathrm{Cov}(\mathbf{X})=\Sigma =({\sigma }_{ij})‘$。因子分析的数学模型为：
X = μ + A F + ε { E ⁡ ( F ) = 0 ,    Cov ⁡ ( F ) = I n E ⁡ ( ε ) = 0 ,    Cov ⁡ ( ε ) = D = diag ⁡ { σ 1 2 , … , σ m 2 } Cov ⁡ ( F , ε ) = 0 \begin{gathered} \mathbf{X}=\boldsymbol{\mu}+AF+\varepsilon \\ \begin{cases} \operatorname{E}(F)=\mathbf{0},\;\operatorname{Cov}(F)=I_n \\ \operatorname{E}(\varepsilon)=\mathbf{0},\;\operatorname{Cov}(\varepsilon)=D=\operatorname{diag}\{\sigma_1^2,\dots,\sigma_m^2\} \\ \operatorname{Cov}(F,\varepsilon)=\mathbf{0} \end{cases} \end{gathered} X=μ+AF+ε⎩ ⎨ ⎧​E(F)=0,Cov(F)=In​E(ε)=0,Cov(ε)=D=diag{σ12​,…,σm2​}Cov(F,ε)=0​​
其中$‘F=(f_1,\dots ,f_n{)}^T‘$是不可观测的$‘n‘$维随机变量，$‘\epsilon ‘$是不可观测的$‘m‘$维随机变量，分别称$‘F‘$和$‘\epsilon ‘$为CommonFactor和SpecificFactor。$‘A=(a_{ij})‘$是一个非随机矩阵，$‘a_{ij}‘$表示公共因子$‘f_j‘$、随机变量$‘{\mathbf{X}}_i‘$的因子载荷。$‘a_{1j},a_{2j},\dots ,a_{ij}‘$中至少有两个不为$‘0‘$，否则可将$‘f_i‘$并入到$‘{\epsilon }_i‘$中去；$‘{\epsilon }_i‘$也仅出现在$‘{\mathbf{X}}_i‘$的表达式中。

上述因子分析模型具有如下性质：

1. $‘\Sigma =A{A}^T+D‘$；

2. 若$‘{\mathbf{X}}^{\ast }=C\mathbf{X}‘$，则有：
   $$\mathbf{Y}=C\mathbf{\mu }+CAF+C\epsilon ={\mathbf{\mu }}^{\ast }+A^{\ast }F+{\epsilon }^{\ast }$$

3. 因子载荷不唯一；

4. $‘\mathrm{Cov}(\mathbf{X},F)=A‘$，即$‘\mathrm{Cov}({\mathbf{X}}_i,F_j)=a_{ij}‘$

5. 令$‘h_i^2=\sum _{j=1}^n{a}_{ij}^2‘$，则有：
   $$\mathrm{Var}({\mathbf{X}}_i)={\sigma }_{ii}=\sum _{j=1}^n{a}_{ij}^2+{\sigma }_i^2=h_i^2+{\sigma }_i^2,i=1,2,\dots ,m$$

6. 令$‘g_j^2=\sum _{i=1}^m{a}_{ij}^2‘$，则有：
   $$\sum _{i=1}^m\mathrm{Var}({\mathbf{X}}_i)=\sum _{j=1}^n{g}_j^2+\sum _{i=1}^n{\sigma }_i^2$$

Proof. (1)由[prop:CovMat](3)(4)(5)可得：
$$\begin{array}{rl}\Sigma & =\mathrm{Cov}(\mathbf{X})=\mathrm{Cov}(\mathbf{\mu }+AF+\epsilon ,\mathbf{\mu }+AF+\epsilon )\\ & =\mathrm{Cov}(\mathbf{\mu },\mathbf{\mu }+AF+\epsilon )+\mathrm{Cov}(AF,\mathbf{\mu }+AF+\epsilon )+\mathrm{Cov}(\epsilon ,\mathbf{\mu }+AF+\epsilon )\\ & =\mathrm{Cov}(AF,\mathbf{\mu })+\mathrm{Cov}(AF)+\mathrm{Cov}(AF,\epsilon )+\mathrm{Cov}(\epsilon ,\mathbf{\mu })+\mathrm{Cov}(\epsilon ,AF)+\mathrm{Cov}(\epsilon )\\ & =A\mathrm{Cov}(F)A^T+A\mathrm{Cov}(F,\epsilon )+\mathrm{Cov}(\epsilon ,F)A^T+D\\ & =A{A}^T+D\end{array}$$

(2)显然。

(3)取正交矩阵$‘Q‘$，令$‘A^{\ast }=AQ‘$，$‘F^{\ast }=Q^{T}F‘$，则依然有：
$$E(F^{\ast })=Q^{T}E(F)=0,\mathrm{Cov}(F^{\ast })=Q^T\mathrm{Cov}(F)Q=I_n,\mathbf{X}=\mathbf{\mu }+A^{\ast }{F}^{\ast }+\epsilon$$

(4)由[prop:CovMat](3)(4)(5)可得：
$$\mathrm{Cov}(\mathbf{X},F)=\mathrm{Cov}(\mathbf{\mu }+AF+\epsilon ,F)=\mathrm{Cov}(\mathbf{\mu },F)+\mathrm{Cov}(AF,F)+\mathrm{Cov}(\epsilon ,F)=A$$

(5)由(1)即可得到结论。

(6)由(1)可得：
$$\begin{array}{rl}\sum _{i=1}^m\mathrm{Var}({\mathbf{X}}_i)& =\mathrm{tr}[\mathrm{Cov}(\mathbf{X})]=\mathrm{tr}(A{A}^T+D)=\sum _{i=1}^m\sum _{j=1}^n{a}_{ij}^2+\sum _{i=1}^n{\sigma }_i^2\\ & =\sum _{j=1}^n\sum _{i=1}^m{a}_{ij}^2+\sum _{i=1}^n{\sigma }_i^2=\sum _{j=1}^n{g}_j^2+\sum _{i=1}^n{\sigma }_i^2\end{array}$$
◻

称$‘h_i^2‘$为变量$‘{\mathbf{X}}_i‘$的CommonVariance，它反映了公共因子对$‘{\mathbf{X}}_i‘$的方差贡献度。称$‘{\sigma }_i^2‘$为$‘{\mathbf{X}}_i‘$的SpecificVariance，它反映了特殊因子$‘{\epsilon }_i‘$对$‘{\mathbf{X}}_i‘$的方差贡献度。$‘g_j^2‘$可视为公共因子$‘f_j‘$对$‘{\mathbf{X}}_1,\dots ,{\mathbf{X}}_m‘$的总方差贡献度。

参数估计方法

主成分法

设观测变量$‘\mathbf{X}‘$的协方差矩阵$‘\Sigma ‘$，它的特征值从大到小依次为$‘{\lambda }_1,\dots ,{\lambda }_m‘$，对应的单位正交特征向量分别为$‘l_1,\dots ,l_m‘$。于是$‘\Sigma ‘$有分解式：
$$\Sigma =\left(\begin{array}{cccc}{l}_1& l_2& \cdots & l_m\end{array}\right)\left(\begin{array}{cccc}{\lambda }_1& 0& \cdots & 0\\ 0& {\lambda }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\lambda }_m\end{array}\right)\left(\begin{array}{c}{l}_1^T\\ l_2^T\\ ⋮\\ l_m^T\end{array}\right)=\sum _{i=1}^m{\lambda }_i{l}_i{l}_i^T$$
由[prop:CovMat](2)和[theo:PositiveSemidefinite](3)的第五条可知$‘{\lambda }_m⩾0‘$。当最后$‘m-n‘$个特征值较小时，$‘\Sigma ‘$有如下近似：
$$\Sigma =\sum _{i=1}^m{\lambda }_i{l}_i{l}_i^T\approx \sum _{i=1}^n{\lambda }_i{l}_i{l}_i^T+D=A{A}^T+D$$
其中：
A = ( λ 1 l 1 ⋯ λ n l n ) ,    D = diag ⁡ { σ 1 2 , … , σ m 2 } ,    σ i 2 = σ i i − h i 2 A= \begin{pmatrix} \sqrt{\lambda_1}l_1 & \cdots & \sqrt{\lambda_n}l_n \end{pmatrix},\; D=\operatorname{diag}\{\sigma_1^2,\dots,\sigma_m^2\},\; \sigma_i^2=\sigma_{ii}-h_i^2 A=(λ1​ ​l1​​⋯​λn​ ​ln​​),D=diag{σ12​,…,σm2​},σi2​=σii​−hi2​
与PCA一样，一般通过使$‘\left(\sum _{i=1}^n{\lambda }_i\right)/\left(\sum _{i=1}^m{\lambda }_i\right)‘$大于一定比例来选择$‘n‘$的具体值。

主因子法

令$‘A{A}^T=\Sigma -D‘$。取$‘{\hat{\sigma }}_1^2,\dots ,{\hat{\sigma }}_m^2‘$为特殊方差的合理初始估计（(1)全零，(2)取$‘\underset{j\ne i}{\max }{\sigma }_{ij}‘$），则有：
$$\widehat{A{A}^T}=\left(\begin{array}{cccc}{\sigma }_{11}-{\hat{\sigma }}_1^2& {\sigma }_{12}& \cdots & {\sigma }_{1m}\\ {\sigma }_{21}& {\sigma }_{22}-{\hat{\sigma }}_2^2& \cdots & {\sigma }_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ {\sigma }_{m1}& {\sigma }_{m2}& \cdots & {\sigma }_{mm}-{\hat{\sigma }}_m^2\end{array}\right)$$
取$‘\widehat{A{A}^T}‘$前$‘n‘$个大于$‘0‘$的特征值，从大到小依次为$‘{\hat{\lambda }}_1,\dots ,{\hat{\lambda }}_n‘$，对应的单位正交特征向量为$‘{\hat{l}}_1,\dots ,{\hat{l}}_n‘$，则有近似的：
$$\hat{A}=\left(\begin{array}{ccc}\sqrt{{\hat{\lambda }}_1}{\hat{l}}_1& \cdots & \sqrt{{\hat{\lambda }}_n}{\hat{l}}_n\end{array}\right)$$
令$‘{\hat{\sigma }}_i^2={\sigma }_{ii}-{\hat{h}}_i^2‘$，继续上面的迭代过程以得到稳定的近似解。

Input: 协方差矩阵 $‘\Sigma ‘$，初始特殊方差估计
$‘{\hat{\sigma }}_1^2,\dots ,{\hat{\sigma }}_m^2‘$，目标因子数 $‘n‘$
Output: 因子载荷矩阵估计 $‘\hat{A}‘$，特殊方差估计
$‘{\hat{\sigma }}_i^2‘$

初始化 $‘{\hat{\sigma }}_i^2‘$ 为合理值 构造矩阵
$‘\widehat{A{A}^T}=\Sigma -\mathrm{diag}({\hat{\sigma }}_1^2,\dots ,{\hat{\sigma }}_m^2)‘$
对 $‘\widehat{A{A}^T}‘$ 做特征值分解，得到部分特征值
$‘{\hat{\lambda }}_1⩾\cdots ⩾{\hat{\lambda }}_n‘$，及对应单位正交特征向量
$‘{\hat{l}}_1,\dots ,{\hat{l}}_n‘$ 构造因子载荷矩阵估计：
$‘\hat{A}=({\hat{a}}_{ij})=\left(\begin{array}{ccc}\sqrt{{\hat{\lambda }}_1}{\hat{l}}_1& \cdots & \sqrt{{\hat{\lambda }}_n}{\hat{l}}_n\end{array}\right)‘$ 令
$‘{\hat{h}}_i^2=\sum _{j=1}^n{\hat{a}}_{ij}^2‘$，更新
$‘{\hat{\sigma }}_i^2={\sigma }_{ii}-{\hat{h}}_i^2,i=1,2,\dots ,m‘$

因子旋转

为了提高因子的可解释性，我们希望每个因子对观测变量的影响是集中且明显的，即一个因子主要对少数几个变量有显著影响，对其余变量几乎没有作用。这种结构反映在因子载荷矩阵$‘A‘$上即为$‘A‘$每一列的元素$‘a_{ij},i=1,2,\dots ,m‘$不是均匀地分布在中间水平，而是趋于两极分化：其绝对值要么接近于$‘0‘$，要么较大。这样可以使得每个因子更容易被识别和解释——因为它只与一小组变量高度相关。这种结构等价于希望载荷矩阵$‘A‘$的每一列具有稀疏性，从而便于赋予因子明确的语义标签。

由[prop:FactorAnalysis](3)可知在初步求得因子载荷矩阵$‘A‘$后，可以使用一个正交矩阵右乘$‘A‘$，此时仍能得到一个因子模型。使用正交矩阵来右乘$‘A‘$相当于是对因子$‘F‘$进行旋转变换，我们可以通过不断旋转$‘F‘$来得到更加稀疏的因子载荷矩阵，从而提高因子的可解释性。

如何旋转？怎么衡量旋转后因子载荷矩阵的优良性？

令：
$$d_{ij}^2=\frac{a_{ij}^2}{h_i^2},i=1,2,\dots ,m,j=1,2,\dots ,n$$
$‘d_{ij}^2‘$衡量了因子$‘j‘$对观测变量$‘{\mathbf{X}}_i‘$的影响，且消除了$‘a_{ij}‘$的正负号带来的差异和各观测变量在因子载荷大小上的不同带来的差异。定义第$‘j‘$列$‘p‘$个数据$‘d_{ij}^2,i=1,2,\dots ,m‘$的方差为：
$$\begin{array}{rl}{V}_j& =\frac{1}{m}\sum _{i=1}^m(d_{ij}^2-{\bar{d}}_j{)}^2=\frac{1}{m}\sum _{i=1}^m\left(d_{ij}^2-\frac{1}{p}\sum _{i=1}^m{d}_{ij}^2\right)\\ & =\frac{1}{m}\left[\sum _{i=1}^m{d}_{ij}^4-m\frac{1}{m^2}{\left(\sum _{i=1}^m{d}_{ij}^2\right)}^2\right]\\ & =\frac{1}{m^2}\left[m\sum _{i=1}^m{d}_{ij}^4-\frac{1}{m}{\left(\sum _{i=1}^m{d}_{ij}^2\right)}^2\right]\\ & =\frac{1}{m^2}\left[m\sum _{i=1}^m\frac{a_{ij}^4}{h_i^4}-\frac{1}{m}{\left(\sum _{i=1}^m\frac{a_{ij}^2}{h_i^2}\right)}^2\right]\end{array}$$
若$‘V_j‘$越大，则第$‘j‘$个因子对观测变量的影响越集中。定义因子载荷矩阵$‘A‘$的方差为：
$$V=\sum _{j=1}^n{V}_j=\frac{1}{m^2}\left\{\sum _{j=1}^n\left[m\sum _{i=1}^m\frac{a_{ij}^4}{h_i^4}-\frac{1}{m}{\left(\sum _{i=1}^m\frac{a_{ij}^2}{h_i^2}\right)}^2\right]\right\}$$
若$‘V‘$越大，则表明因子对观测变量的影响越集中。

综上，我们只需使得旋转后得到的因子载荷矩阵$‘A‘$的方差$‘V‘$达到最大即可。

代码实现

R语言中使用Factanal函数进行因子分析，注意它使用极大似然估计法进行参数估计。

Reference

1.薛毅，统计建模与R软件