6.4 因子分析
因子分析是一种降维的方法，通过探寻变量的潜在结构，将其归纳为较少的因子，从而实现降维的目的。
6.4.1 正交因子模型
设X为可观测的p维随机向量，其均值为\(\mu\)，协方差矩阵为\(\Sigma\)，则正交因子模型可被表示为
\[ \begin{array}{c} X=AF+\mu+\varepsilon \\ E(F)=0, \; Cov(F)=I_m \\ E(\varepsilon)=0, \; Cov(\varepsilon)=\Psi=diag(\psi_1,...,\psi_2) \\ Cov(F,\varepsilon)=0 \end{array} \tag{6.48} \]
可通过中心化处理而消去\(\mu\)
其中\(A\)为\(p \times m\)的因子载荷矩阵，\(F=(F_1, ...,F_m)'\)为公共因子向量，\(\varepsilon=(\varepsilon_1,...,\varepsilon_p)'\)为特殊因子向量。
则X的协差阵可表示为
\[ \Sigma = Cov(X)=Cov(AF+\mu+\varepsilon)=ACov(F)A'+Cov(\varepsilon)=AA'+\Psi \tag{6.49} \]
若X已进行标准化处理，则
\[ R=AA'+\Psi \tag{6.50} \]
由于\(\Psi\)是对角阵，则\(\Sigma\)和\(R\)的非对角线元素全部由因子载荷矩阵A决定
- 因子载荷矩阵的统计意义
\[ \begin{aligned} Cov(X,F)&=E[X-E(X)][F-E(F)]' \\ &=E[(X-\mu)F'] \\ &= E[(AF+\varepsilon)F'] \\ &= AE(FF')+E(\varepsilon F') \\ &= A \end{aligned} \tag{6.51} \]
由此可知，对于原始数据而言，因子载荷矩阵A记录了X与F之间的协方差，即\(a_{ij}=Cov(X_i,F_j)\)。
若对原始数据X进行标准化处理得到\(X^*\)，则
\[ \rho_{ij}=\frac{Cov(X_i,F_j)}{\sqrt{Var(X_i)}\sqrt{Var(F_j)}}=Cov(X_i, F_j)=a_{ij} \tag{6.52} \]
- 变量共同度的统计意义
由式(6.48)可得
\[ Var(X_i)=\sigma_{ii}=\sum_{j=1}^m a_{ij}^2+\psi_i:=h_i+\psi_i \tag{6.53} \]
其中\(h_i\)表现为全部公共因子对\(X_i\)的变异程度的解释，反映了\(X_i\)对公共因子的依赖程度，称为变量共同度，即因子载荷矩阵A的行元素平方和。
- 公共因子方差贡献的统计意义
\[ \sum_{i=1}^p Var(X_i) = \sum_{i=1}^p a_{i1}^2Var(F_1)+...+\sum_{i=1}^pa_{im}^2Var(F_m)+\sum_{i=1}^pVar(\varepsilon):=\sum_{j=1}^m q_j+\sum_{i=1}^p \psi_i \tag{6.54} \]
其中\(q_j=\sum_{i=1}^p a_{ij}^2\)，是因子载荷矩阵的列元素之和，反映了第j个公共因子对变量X的贡献及其相对重要性。进一步地，称\(\frac{q_j}{\sum_{i=1}^pVar(X_i)}\)为公共因子\(F_j\)所解释的总方差的比例。
性质：
因子分析具有尺度不变性
即对随机向量X的各个分量进行尺度放缩\(X^*=CX, \; C=diag(c_1,...,c_p)\)后的随机向量\(X^*\)仍满足正交因子模型的条件，且公共因子\(F\)不变。
\(CX=C\mu+CAF+C\varepsilon\)，可将\(CA\)视作\(A^*\)因此公共因子\(F\)不变，变的是载荷矩阵
因子载荷矩阵不唯一
对于\(AF\)，可以有任一正交矩阵\(\Gamma\)，使得\(AF=(A\Gamma)(\Gamma'F)\)，其中可将\(A\Gamma\)视作\(A^*\)，将\(\Gamma'F\)视作\(F^*\)，进行因子旋转后仍保持正交因子模型的条件。
正交矩阵的作用相当于进行旋转
6.4.2 参数估计
在参数估计前，得先判断变量之间是否存在足够强的相关关系来进行因子分析。
KMO检验
KMO检验用于检查变量间的相关性和偏相关性。KMO统计量的取值越接近1，表明变量间的相关性越强，偏相关性越弱，因子分析的效果越好。KMO统计量的值一般在0.7以上比较适合做因子分析。
R语言：
KMO()
Bartlett球形检验
Bartlett球形检验的原假设是变量的相关阵是单位阵，当拒绝原假设时，则说明变量间具有相关性，因子分析有效。
R语言：
psych::cortest.bartlett()
6.4.2.1 主成分法
令\(\Sigma\)的特征值为\(\lambda_1 \geq \lambda_2 \geq ... \geq \lambda_p \geq 0\)，对\(\Sigma\)进行谱分解
\[ \begin{aligned} \Sigma &= P\Lambda P' \\ &= \sum_{i=1}^p \lambda_ip_ip_i' \\ &= \sum_{i=1}^m \lambda_ip_ip_i'+\sum_{i=m+1}^p \lambda_ip_ip_i' \\ &\approx \begin{pmatrix} \sqrt{\lambda_1}p_1 & ... & \sqrt{\lambda_m}p_m \end{pmatrix} \begin{pmatrix} \sqrt{\lambda_1}p_1' \\ \vdots \\ \sqrt{\lambda_m}p_m' \end{pmatrix} + \begin{pmatrix} \psi_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \psi_p \end{pmatrix} \\ &= \hat A \hat A'+\hat \Psi \end{aligned} \tag{6.55} \]
对\(\Sigma\)进行谱分解，取特征值较大的前m个特征值及对应的特征向量构成因子载荷矩阵A，忽视掉其余较小的特征值及对应的特征向量。其中的\(\psi_i\)则通过\(\hat \psi_i = \sigma_{ii}-\sum_{j=1}^m a_{ij}^2\)计算得到。
此时因子载荷矩阵中的每一列都是\(\sqrt{\lambda_i}p_i\)。特别的，当\(\lambda_{m+1}=...=\lambda_p=0\)时，此时有\(\Sigma=\hat A \hat A'\)，表明\(\hat A=P\Lambda^{\frac{1}{2}}\)就是真实的载荷矩阵
当然也可以对变量进行标准化后利用主成分法进行求解。此时公共因子\(F_j\)的贡献为\(q_j=\sum_{i=1}^p \hat a_{ij}^2=||\sqrt{\lambda_j}p_j||^2 = \lambda_j\)
m的选择可参考主成分分析
6.4.2.2 主因子法
由于因子具有尺度不变性，对X进行标准化后得到
\[ R=AA'+\Psi \tag{6.56} \]
若\(R\)和\(\Psi\)都已知，则称\(R^*=R-\Psi=AA'\)为“约相关矩阵”。
其中\(R\)可表示为
\[ \begin{cases} \rho_{ij}=a_{i1}a_{j1}+a_{i2}a_{j2}+...+a_{im}a_{jm} \quad &1 \leq i \neq j \leq q \\ \rho_{ii} = a_{i1}^2 + a_{i2}^2 + ...+a_{im}^2+\psi_i = h_i+\psi_i=1 \quad &i=1,2,...,p \end{cases} \tag{6.57} \]
\(R^*\)可表示为
\[ \begin{cases} \rho_{ij}=a_{i1}a_{j1}+a_{i2}a_{j2}+...+a_{im}a_{jm} \quad &1 \leq i \neq j \leq q \\ \rho_{ii} = a_{i1}^2 + a_{i2}^2 + ...+a_{im}^2= h_i \quad &i=1,2,...,p \end{cases} \tag{6.58} \]
关键点为\(h_i=1-\psi_i\)
设\(R^*\)的特征值为\(\lambda_1^* \geq \lambda_2^* \geq ... \geq \lambda_p^* \geq 0\)，对应的特征向量为\(p_i^*\)，此时可对\(R^*\)进行谱分解，得到
\[ R^*=R-\Psi=\hat A \hat A' \approx \begin{pmatrix} \sqrt{\lambda_1^*}p_1^* & ... & \sqrt{\lambda_m^*}p_m^* \end{pmatrix} \begin{pmatrix} \sqrt{\lambda_1^*}p_1^{*'} \\ \vdots \\ \sqrt{\lambda_m^*}p_m^{*'} \end{pmatrix} \tag{6.59} \]
此时\(\hat A = \begin{pmatrix} \sqrt{\lambda_1^*}p_1^* & ... & \sqrt{\lambda_m^*}p_m^* \end{pmatrix}\)，特殊因子方差的估计为\(\hat \psi_i = 1-\sum_{j=1}^m \hat a_{ij}^2\)
主因子法的前提是得知道\(\Psi\)，但在实际中\(\Psi\)一般未知，因此可根据样本进行估计。也可根据迭代算法，根据初始值\(\Psi^{(0)}\)开始迭代，得到\(\hat A^{(1)}\)，再根据\(\hat A^{(1)}\)得到\(\Psi^{(1)}\)，不断重复，直至收敛
6.4.2.3 极大似然法
假设公共因子\(F \sim N_m(0,I)\)，特殊因子\(\varepsilon\sim N_p(0,\Psi)\)，且相互独立，则\(X \sim N_p(\mu, \Sigma)\)。构造似然函数如下所示
\[ L(\mu, \Sigma)=(2\pi)^{-\frac{np}{2}}|\Sigma|^{-\frac{n}{2}}\exp[-\frac{1}{2}\sum_{i=1}^n(X_{(i)}-\mu)'\Sigma^{-1}(X_{(i)}-\mu)] \tag{6.60} \]
其中\(\mu\)的极大似然估计为\(\hat \mu = \bar X\)，回代可得
\[ \begin{aligned} L(\Sigma)&=(2\pi)^{-\frac{np}{2}}|\Sigma|^{-\frac{n}{2}}\exp[-\frac{1}{2}\sum_{i=1}^n(X_{(i)}-\bar X)'\Sigma^{-1}(X_{(i)}-\bar X)] \\ &= (2\pi)^{-\frac{np}{2}}|\Sigma|^{-\frac{n}{2}}\exp[-\frac{1}{2}tr(\sum_{i=1}^n(X_{(i)}-\bar X)'\Sigma^{-1}(X_{(i)}-\bar X))] \\ &= (2\pi)^{-\frac{np}{2}}|\Sigma|^{-\frac{n}{2}}\exp[-\frac{1}{2}tr(\sum_{i=1}^n(X_{(i)}-\bar X)(X_{(i)}-\bar X)'\Sigma^{-1})] \\ &= (2\pi)^{-\frac{np}{2}}|\Sigma|^{-\frac{n}{2}}\exp[-\frac{1}{2}tr(D\Sigma^{-1})] \end{aligned} \tag{6.61} \]
其中\(D=\sum_{i=1}^n(X_{(i)}-\bar X)(X_{(i)}-\bar X)'\)为样本离差阵。
根据\(\Sigma = AA'+\Psi\)，则
\[ L(A,\Psi)=(2\pi)^{-\frac{np}{2}}|AA'+\Psi|^{-\frac{n}{2}}\exp[-\frac{1}{2}tr(D(AA'+\Psi)^{-1})] \tag{6.62} \]
求解\(\hat A\)和\(\hat \Psi\)使得对数似然函数\(\ln L(A,\Psi)\)达到最大即可。