应用多元统计分析-因子分析
因子分析是主成分分析的推广和拓展，它是多元统计分析中降维的一种方法，是一种用来分析隐藏在表面现象背后的因子作用的一类统计模型。因子分析研究相关阵或协差阵的内部依赖关系。研究的是是否存在更少的因子，可以解释原始变量与因子之间的相关关系。
1 因子模型
首先，我们来看一个简单而直观的例子，来了解因子分析到底是怎样的一种方法。
为了解学生的学习能力，观测了n个学生的p个科目的成绩，用\(X_1,…X_p\)表示p个科目，\(X_i=(x_{i1},…x_{ip})\)表示第i个学生的p个科目的成绩，现在我们想知道的是，哪些因素对学生的学习能力影响是比较大的。
如果我们假设学生的学习能力仅由智能高低和其他影响因素决定，那么就是说，各个科目可以看作是由两部分组成：\[X_i=a_if+\epsilon_i,i=1,2…p\] 其中\(f\)是对所有的\(X_i\)都起作用的公共因子，它是表示智能高低的因子；系数\(a_i\)称为因子载荷；\(\epsilon_i\)是科目\(X_i\)特有的特殊因子，这就是一个最简单的因子模型。
接着呢，我们可以把简单因子模型推广到多个因子的情况。就是说，学生的学习能力由很多因素决定，比如说记忆因子、计算因子、灵活因子等。分别计这些因子为\(f_1,…f_m\),则可以得到\[X_i=a_{i1}f_1+a_{i2}f_2+…+a_{im}f_m+\epsilon_i,i=1,2,…,p\] 用这m个不可观测且互不相关的公共因子\(f_1,…f_m\)和一个特殊因子\(\epsilon\)来描述原始可测的相关变量解释分析学生的学习能力。它们的系数\(a_{i1},…1_{ip}\)称为因子载荷，表示第i个科目在m个方面的表现，这就是一个因子模型。
严格的数学表述为：设\(X=(X_1,…X_P)^T\)是可观测的随机变量，且\[E(X)=\mu=(\mu_1,…\mu_p)^T,var(X)=\Sigma=(\sigma_{ij})_{p\times p}\] 因子分析的一般模型为\[\begin{cases}X_1-\mu_1=a_{11}f_1+a_{12}f_2+…+a_{1m}f_m+\epsilon_1\\X_2-\mu_2=a_{21}f_1+a_{22}f_2+…+a_{2m}f_m+\epsilon_2\\…\\X_p-\mu_p=a_{p1}f_1+a_{p2}f_2+…+a_{pm}f_m+\epsilon_p\end{cases}\]
其中，\(f_1,f_2,…,f_m\)为 公共因子，\(\epsilon_1,\epsilon_2,…\epsilon_p\)为特殊因子，它们都是不可观测的随机变量。公共因子出现在每一个原始变量的表达式中，可以理解为原始变量共同具有的公共因素，每个公共因子一般至少对两个原始变量有作用，否则它将归入特殊因子。上述模型的矩阵形式为\[X=\mu+AF+\epsilon\]
因子分析的主要应用有两个方面，一是寻求基本结构，简化观测系统，将具有错综复杂关系的对象（变量或样本）综合为少数几个因子，以再现因子与原始变量之间的内在联系；二是用于分类，对于p个变量或n个样本进行分类。
在拿到一个数据集后，很多人都很疑惑，不知道因子分析是对特征变量分析还是对样本分析。其实呢，因子分析既可以对特征变量分析（R型因子分析），也可以对样本分析（Q型因子分析）。R型因子分析研究变量之间的相关关系，通过对变量的相关阵或协方差内部结构的研究，找出控制所有变量的几个公共因子。Q型因子分析研究样本之间的相关关系，通过对样本的相似矩阵内部结构的研究找出控制样本的几个因素。
值得注意得到是，公共因子的个数一般有两种确定方法，一是根据实际问题的意义或专业理论知识来确定；二是用确定主成分个数的原则来确定。
2 参数估计
这里的参数估计主要指的是公共因子载荷矩阵和特殊因子的方差矩阵的估计。估计方法一般有三种，分别是主成分法和主因子法。
2.1 主成分法
设样本相关阵S的特征值为\(\lambda_1\geq\lambda_2,…\geq\lambda_p\)，相应的单位正交特征向量为\(l_1,l_2,…l_p\)，当最后p-m个特征值较小时，则样本相关阵可以近似分解为\[S\approx\lambda_1l_1l_1’+…+\lambda_ml_ml_m’=AA’+D=(\sqrt{\lambda_1}l_1,…,\sqrt{\lambda_m}l_m)(\sqrt{\lambda_1}l_1,…,\sqrt{\lambda_m}l_m)’\\+diag(\sigma_1^2,\sigma_2^2,…\sigma_p^2)\] 式中给出的\(\hat{A}\)和\(\hat{D}\)即为因子模型的一个解。
2.2 主因子法
从相关阵出发，设\(R=AA’+D\)，则\(R-D=AA’=R^*\)称为简约相关阵。计算简约相关阵的特征值和单位正交特征向量，并且选取前m个特征值\(\lambda_1^*\geq\lambda_2^*,…\geq\lambda_p^*\) ，相对应的单位正交特征向量为\(l_1^*,l_2^*,…l_p^*\)，则有近似分解式为:\[R\approx\hat{A}\hat{A}’\]其中，\[\hat{A}=(\sqrt{\hat{\lambda_1^*}}\hat{l_1^*},\sqrt{\hat{\lambda_2^*}}\hat{l_2^*},…\sqrt{\hat{\lambda_m^*}}\hat{l_m^*})\]\[\hat{\sigma_i}=1-(\hat{h_i})^2\]\[\hat{D}=diag(\hat{\sigma_1},\hat{\sigma_2},…\hat{\sigma_p})\] 则\(\hat{A}\)和\(\hat{D}\)为因子的模型的一个解。如果希望求得的解更加精确，则可以利用迭代的方法，直到解稳定为止。
2.3 极大似然法
该方法具体思想主要是运用了求参数极大似然估计的思想，这种思想我们在参数估计章节中已具体讲到。以下是极大似然法的步骤：
首先，给定初始矩阵\(D_0\)，计算\(D_0^{1/2}SD_0^{1/2}\)的特征值和单位正交特征向量,\(\Theta_0\)为前m个特征值组成的对角阵，\(L_0\)为相应的前m个单位正交向量为列向量的矩阵。得到\(A_1=D_0L_0(\Theta_0-I_m)^{1/2}\),\(D_1=diag(S-A_1A_1′)\);然后再按上述方法得到\(A_2\)和\(D_2\)，直到得到的解稳定为止。
3 因子旋转
因子分析的目的不仅是求出公共因子，更主要的是应该知道每个公共因子实际意义。什么叫做实际意义？我们知道，因子分析模型中的公共因子的系数\(a_{ij}\)表示的是第i个变量与第j个公共因子的相关关系，其绝对值越大表示相关的密切程度越高。
那么现在我们想知道的是，对每一个因子来说，受其影响最大的变量有哪些？受其影响最小的变量有哪些，如果知道了这些，我们就可以很轻松地解释每一个因子的价值体现在哪里。那么，我们就想通过因子载荷矩阵的变换使得各个公共因子的典型代表变量突出一些，也就是使放大载荷大的值，缩小载荷小的值，但同时我们要保证这些载荷在同一度量下。采取的方法就是对载荷矩阵进行正交变换，使得因子载荷的每一列元素的平方向0或1两级转化，通常采用正交旋转得到方差最大的载荷矩阵。
数学表述为：设因子模型是\(X=AF+\epsilon\)，其中\(F\)为公因子向量，对\(F\)施行正交变换，令\(Z=\Gamma^TF\)(\(\Gamma\)为任一m阶正交矩阵)，则\[X=A\Gamma Z+\epsilon\]
4 因子得分
当我们得到公共因子和因子载荷后，就应当反过来考察每一个样本，比如我们最开始给出的例子，在得到一个智能因子后，应该考察所有同学在这个记忆因子、计算因子、灵活因子这3个因子上的得分情况，以便于我们挑选出优秀的同学和需要帮助的同学。估计因子得分常用的两种方法为加权最小二乘法和回归方法。
4.1 加权最小二乘法
根据因子模型的假设条件，我们考虑加权最小二乘函数\[g(F)=(x-AF)’D^{-1}(x-AF)\] 求\(F\)的估计值\(\hat{F}\),使得\[g(\hat{F})=min(g(F))\]则\[\hat{F}=(A’D^{-1}A)^{-1}A’D^{-1}x\]这就是因子得分的加权最小二乘估计。
4.2 回归法
\[x=AF+\varepsilon\]\[F\sim N(0,I_m),\varepsilon\sim N_p(0,D)\]\[F和\varepsilon相互独立\]则通过统计知识可以得到回归法计算出来的因子得分为\(\hat{F}=A’\Sigma^{-1}x\)。
5 R语言实现
在R语言中，用factanal()函数作因子分析，其中是运用极大似然法进行参数估计，函数的基本调用格式为
factanal(x，factors,data,...score=c("none", "regression", "Bartlett"))例如，有48名应聘者，公司为这48名应聘者的15项指标标分，用因子分析方法对这15项指标做因子分析，在因子分析中选取5个因子。
factor<-read.csv("E:/BJintern/factor.csv") # 读取数据
head(factor)#显示数据集前6行信息## FL APP AA LA SC LC HON SMS EXP DRV AMB GSP POT KJ SUTT
## 1 6 7 2 5 8 7 8 8 3 8 9 7 5 7 10
## 2 9 10 5 8 10 9 9 10 5 9 9 8 8 8 10
## 3 7 8 3 6 9 8 9 7 4 9 9 8 6 8 10
## 4 5 6 8 5 6 5 9 2 8 4 5 8 7 6 5
## 5 6 8 8 8 4 4 9 5 8 5 5 8 8 7 7
## 6 7 7 7 6 8 7 10 5 9 6 5 8 6 6 6fa=factanal(~.,factor=5,data=factor,score="regression") # 对所有变量进行因子分析
fa#显示模型各项信息##
## Call:
## factanal(x = ~., factors = 5, data = factor, scores = "regression")
##
## Uniquenesses:
## FL APP AA LA SC LC HON SMS EXP DRV AMB GSP
## 0.439 0.597 0.509 0.197 0.118 0.005 0.292 0.140 0.365 0.223 0.098 0.119
## POT KJ SUTT
## 0.084 0.005 0.267
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4 Factor5
## FL 0.127 0.722 0.102 -0.117
## APP 0.451 0.134 0.270 0.206 0.258
## AA 0.129 0.686
## LA 0.222 0.246 0.827
## SC 0.917 0.167
## LC 0.851 0.125 0.279 -0.420
## HON 0.228 -0.220 0.777
## SMS 0.880 0.266 0.111
## EXP 0.773 0.171
## DRV 0.754 0.393 0.199 0.114
## AMB 0.909 0.187 0.112 0.165
## GSP 0.783 0.295 0.354 0.148 -0.181
## POT 0.717 0.362 0.446 0.267
## KJ 0.418 0.399 0.563 -0.585
## SUTT 0.351 0.764 0.148
##
## Factor1 Factor2 Factor3 Factor4 Factor5
## SS loadings 5.490 2.507 2.188 1.028 0.331
## Proportion Var 0.366 0.167 0.146 0.069 0.022
## Cumulative Var 0.366 0.533 0.679 0.748 0.770
##
## Test of the hypothesis that 5 factors are sufficient.
## The chi square statistic is 60.97 on 40 degrees of freedom.
## The p-value is 0.0179从结果可以看出，第一公因子主要表现应聘者的外露能力，第二公因子主要表现应聘者的经验，第三公因子主要表现求职者是否讨人喜欢，第四公因子主要表现为求职者的专业能力，第五公因子主要表现为求职者的外貌。
为直观起见，画出这48名应聘者在第1、第2公共因子下的散点图
plot(fa$scores[,1:2],type="n")
text(fa$scores[,1],fa$scores[,2])
我们已经知道第1公共因子主要表现求职者外露能力，第2公共因子主要表现求职者的经验。公司可以选择两者得分都比较高的应聘者。如39，40，7，8等应聘者。如偏重外露能力，则选取第1公共因子得分较大的应聘者。如偏重经验，则可以考虑第2公共因子得分较大的应聘者。公司也可以根据情况，依据此法分析其他公共因子。
6 参考文献
[1]薛毅. 统计建模与R软件[M]. 清华大学出版社, 2007.