多元统计分析——数据降维——因子分析（FA）

xia ge tou lia

已于 2023-06-29 22:52:17 修改

阅读量2.5w

点赞数 34

分类专栏：机器学习 sklearn fa-kit 文章标签： python 机器学习人工智能

于 2020-06-22 22:53:09 首次发布

本文链接：https://blog.csdn.net/huangguohui_123/article/details/105869253

版权

机器学习同时被 3 个专栏收录

17 篇文章 39 订阅

订阅专栏

sklearn

8 篇文章 1 订阅

订阅专栏

fa-kit

1 篇文章 1 订阅

订阅专栏

一、因子分析简介

1、定义

1904年，英国心理学家CharlesSpearman研究了33名学生在古典语、法语和英语三门成绩，三门成绩的相关性系数如下：

三门成绩的高度相关会不会是由于它们三个成绩的背后有一个共同的因素，来决定这三门成绩的？比如语言能力？

对于个原始变量 $Y_{1},...,Y_{p}$ 来说，那些高度相关的变量很可能会遵循一个共同的潜在结构——或可称之为公共因子 (Common factor） 。简单的说就是：公共因子是用一个共同的因素来刻画几个高度相关的变量。

然而，这些“公共因子”通常是无法观测的，故称为潜变量 (latentvariables)。这在心理学、社会学及行为科学等学科中非常常见，比如
“智力”和“社会阶层”。

因子分析（Factor analysis）旨在提出因子模型（Factor model）来研究如何用几个公共因子，记作 $F_{1},...,F_{m}$ ，通常 m<p ，来刻画原始变量之间的相关性。

2、正交因子模型

Charles Spearman基于学生3门语言成绩的数据提出了单因子模型（Single factor model）：

$Y_{1}=l_{1}F+\varepsilon _{1}$

$Y_{2}=l_{2}F+\varepsilon _{2}$

$Y_{3}=l_{3}F+\varepsilon _{3}$

其中代表公共因子（Common factor）， $\varepsilon _{j}$ 代表特殊因子（Specific factor），即代表 $Y_{j}$ 的特殊部分； $l_{j}$ 代表系数/载荷（Loading），即来说明公共因子对 $Y_{j}$ 的解释力。

当然，大多数时候一个公共因子是不够的，错综复杂的变量可能需要多个公共因子刻画，这就是我们将要学习的正交因子模型（Orthogonal factor model）

假设可观测随机向量 $y=(Y_1,...,Y_{p})'$ 的均值为 $\mu$ ，协方差矩阵为 $\Sigma$ 。正交因子模型假定线性依赖于个不可观测公共因子 $f=(F_{1},...,F_{m})'$ 和个不可观测的特殊因子 (specific factors) $\varepsilon =(\varepsilon _{1},...,\varepsilon _{p})'$ ，通常 m<p ：

$Y_{1}-\mu _{1}=l_{11}F_{1}+l_{12}F_{2}+..+l_{1m}F_{m}+\varepsilon _{1}$

$Y_{2}-\mu _{2}=l_{21}F_{1}+l_{22}F_{2}+..+l_{2m}F_{m}+\varepsilon _{2}$

$Y_{p}-\mu _{p}=l_{p1}F_{1}+l_{p2}F_{2}+..+l_{pm}F_{m}+\varepsilon _{p}$

系数 $l_{jk}$ 称为第个变量在第个因子上的载荷(loading)，体现了该公共因子对此变量的解释力。

使用矩阵记号，上述模型可写为

$y_{p\times 1}-\mu _{p\times 1}=L_{p\times m}f_{m\times 1}+\varepsilon _{p\times 1}$ ，其中 $L_{p\times m}$ 即称作载荷矩阵(loading matrix)。

我们知道公共因子是虚拟的，不可观测的，所以理论上，我们可以给其任意的假设，为了分析方便，我们给它一种最简单的假设，其中和 $\varepsilon$ 假设满足：

① $E(f)=0_{m\times 1},COV(f)=E(ff')=I_{m\times m}$ （单位阵，即假设是正交的， $F_{j}$ 的方差假设为1， $F_{j}$ 和 $F_{k}(j\neq k)$ 之间假设是不相关的)

② $E(\varepsilon )=0_{p\times 1},COV(\varepsilon )=E(\varepsilon \varepsilon ')=\Psi _{p\times p}=\begin{bmatrix} \Psi _{1} & 0 & ... & 0\\ 0 & \Psi _{1} & ... & 0\\ . & . & & .\\ . &. & &. \\ 0& 0 & 0 & \Psi _{p} \end{bmatrix}$ （对角阵）

③而且和 $\varepsilon$ 是无关的： $COV(\varepsilon ,f)=0_{p\times m}$

有了这些假设，我们就可以回头来看，载荷矩阵到底是如何解释原始变量和公共因子的关系的？

已知模型如下：

矩阵形式模型： $y-\mu=Lf+\varepsilon$

元素形式模型： $Y_{j}-\mu _{j}=l_{j1}F_{1}+l_{j2}F_{2}+..+l_{jm}F_{m}+\varepsilon _{j},j=1,...,p$

我们可以计算原始变量和公共因子的协方差矩阵：

$\begin{align}COV(y,f) &=COV(Lf+\varepsilon ,f) \\& =COV(Lf,f)+COV(\varepsilon,f) \\&=E(Lff')+0_{p\times m} \\&=L \end{align}$

综上，我们得到结论： COV(y,f)=L ，即刻画的是和之间的协方差，因此载荷 $l_{jk}$ 测度了第个变量与第个公共因子之间的关联 $cov(Y_{j},F_{k})=l_{jk}$ 。

原始变量的协方差矩阵

$\begin{align}\Sigma =COV(y)&=COV(Lf+\varepsilon)\\&=COV(Lf)+COV(\varepsilon)+COV(Lf,\varepsilon)\\&=LL'+\Psi \end{align}$

其中 $\Sigma$ 中每一个对角线元素

$\begin{align}\sigma _{jj}=var(Y_{j}) &=var(l_{j1}F_{1}+l_{j2}F_{2}+..+l_{jm}F_{m}+\varepsilon _{j})\\&=var(l_{j1}F_{1})+var(l_{j2}F_{2})+..+var(l_{jm}F_{m})+var(\varepsilon _{j})\\&=l_{j1}^{2}+...+l_{jm}^{2}+\Psi _{j}\\&=h_{j}^{2}+\Psi _{j} \end{align}$

其中 $h_{j}^{2}=l_{j1}^{2}+...+l_{jm}^{2}$

$h_{j}^{2}$ 体现了共同性(communality)，指由个公共因子贡献的方差。

$\Psi _{j}$ 称为唯一性(uniqueness)或个体方差(specific variance)，指无法由公共因子贡献的方差部分。

我们也可以看 $Y_{j}$ 和 $Y_{k}$ 之间的关系： $\sigma _{jk}=cov(Y_{j},Y_{k})=l_{j1}l_{k1}+..+l_{jm}l_{km},j\neq k$ 。

3、因子分析的计算过程

估计因子载荷矩阵
进行因子旋转
估计公共因子（因子得分）

二、载荷矩阵的估计方法

1、主成分法

1.1、基本方法

已知矩阵形式模型： $y-\mu=Lf+\varepsilon$

如何估计载荷矩阵？

上面得到的因子模型的协方差矩阵分解： $\Sigma =L_{p\times m}L_{p\times m}'+\Psi$

基于 $\Sigma$ 的谱分解： $\Sigma =\sum ^{p}_{j=1}\lambda _{j}e_{j}e'_{j}=\Lambda_{p\times p} \Lambda_{p\times p} '$ 其中 $\Lambda =(\sqrt{\lambda _{1}}e_{1},...,\sqrt{\lambda _{p}}e_{p})$ ，而 $\sqrt{\lambda _{j}}e_{j}$ 是 $\Sigma$ 的特征值-特征向量。

对比以下两个式子：

$\Sigma =L_{p\times m}L_{p\times m}'+\Psi$

$\Sigma =\sum ^{p}_{j=1}\lambda _{j}e_{j}e'_{j}=\Lambda_{p\times p} \Lambda_{p\times p} '$

当最后 p-m 个特征值很小的时候，我们可以忽略掉 $\Psi$ ，此时： $\Sigma\approx L_{p\times m}L_{p\times m}'$ ，此时我们定义： $L_{p\times m} =(\sqrt{\lambda _{1}}e_{1},...,\sqrt{\lambda _{m}}e_{m})$ ，即为 $\Lambda_{p\times p}$ 的前列；得到了 $L_{p\times m}$ ，我们就可以算出 $\Psi$ ： $\Psi=diag(\Psi _{1},...,\Psi _{p}),\Psi _{j}=\sigma _{jj}-\sum _{k=1}^{m}l^{2}_{jk},j=1,...,p$ ，即理解为： $\Psi$ 对角阵的每一个元素就是的方差没有被 $L_{p\times m}$ 刻画的部分。

在样本层面， $\Sigma$ 可用样本协方差矩阵代替。这种估计方法称为主成分法 (Principal component method)，或主成分解（Principal component solution）

定理（主成分解）

$(\widehat{\lambda} _{j},\widehat{e}_{j}),j=1,...,p$ 为样本协方差矩阵的特征值-标准化特征向量，其中 $\widehat{\lambda} _{1}\geq ...\geq \widehat{\lambda} _{p}$ 。记为公共因子的数目。那么因子载荷矩阵估计的主成分解为 $\widehat{L} =(\sqrt{\widehat{\lambda} _{1}}\widehat{e}_{1},...,\sqrt{\widehat{\lambda} _{m}}\widehat{e}_{m})$ 。估计的个体方差是 $S-\widehat{L}\widehat{L}'$ 的对角元素，即 $\widehat{\Psi}=diag(\widehat{\Psi} _{1},...,\widehat{\Psi} _{p})$ ，其中 $\widehat{\Psi} _{j}=s _{jj}-\sum _{k=1}^{m}\widehat{l}^{2}_{jk},j=1,...,p$ （ $\widehat{\Psi} _{j}$ 代表个体方差， $s _{jj}$ 代表总方差， $\sum _{k=1}^{m}\widehat{l}^{2}_{jk}$ 代表共同度 $\widehat{h}_{j}^{2}=\widehat{l}_{j1}^{2}+...+\widehat{l}_{jm}^{2}$ ）。

【小贴士】

的对角元素等于 $\widehat{L}\widehat{L}'+\widehat{\Psi }$ 的对角元素。然而，的非对角元通常不能由 $\widehat{L}\widehat{L}'+\widehat{\Psi }$ 复原。

第个的因子的载荷估计 $\sqrt{\widehat{\lambda} _{k}}\widehat{e}_{k}$ 与第个主成分的系数 $\widehat{e}_{k}$ 成比例（因子载荷估计其实就是加权的主成分）。

1.2、因子的方差（信息量）贡献

$Y_{1}-\mu _{1}=l_{11}F_{1}+l_{12}F_{2}+..+l_{1m}F_{m}+\varepsilon _{1}$

$Y_{2}-\mu _{2}=l_{21}F_{1}+l_{22}F_{2}+..+l_{2m}F_{m}+\varepsilon _{2}$

$Y_{p}-\mu _{p}=l_{p1}F_{1}+l_{p2}F_{2}+..+l_{pm}F_{m}+\varepsilon _{p}$

第个因子对总方差 $s_{11}+s_{22}+...+s_{pp}=tr(S)$ 的贡献为 $\sum ^{p}_{j=1}\widehat{l}^{2}_{jk}=\widehat{\lambda }_{k}$ （主成分中，第个主成分的方差贡献也为 $\widehat{\lambda }_{k}$ ，要注意的是主成分当中， $\widehat{\lambda }_{k}$ 代表的是第个主成分的方差；但是在因子分析当中，各个公共因子的方差我们假设为1， $\widehat{\lambda }_{k}$ 代表的是同一个公共因子（例如 $F_{2}$ ）对 $Y_{j},j=1,...,p$ 的方差贡献之和 $\sum ^{p}_{j=1}\widehat{l}^{2}_{j2}=\widehat{l}_{12}^{2}+\widehat{l}_{22}^{2}+...+\widehat{l}_{p2}^{2}$ ），贡献比例为 $\frac{\widehat{\lambda }_{k}}{tr(S)},k=1,...,m$ 。

以第2个因子 $F_{2}$ 为例，我们知道 $Y_{j},j=1,...,p$ 由对应的公共因子 $F_{2}$ 贡献的方差是 $l_{j2}^{2},j=1,...,p$ ，则总方差 $s_{11}+s_{22}+...+s_{pp}=tr(S)$ 也可以看做是 $Y_{j},j=1,...,p$ 的方差之和，则第2个因子对总方差的贡献为 $\sum ^{p}_{j=1}\widehat{l}^{2}_{j2}=\widehat{l}_{12}^{2}+\widehat{l}_{22}^{2}+...+\widehat{l}_{p2}^{2}$ ，即为载荷矩阵的第2列的平方和，载荷矩阵的第2列的估计为 $\sqrt{\lambda _{2}}e_{2}$ ， $(\sqrt{\lambda _{2}}e_{2})^{2}=(\sqrt{\lambda _{2}})^{2}e_{2}e'_{2}=\lambda _{2}$ 。

1.3、基于标准化变量的主成分法

正如主成分分析中提到的，当变量的度量单位/量纲相差很多时，建议使用标准化后的变量（相关系数矩阵）进行因子分析，这等价于基于样本相关系数矩阵使用主成分法，即将换成。

比更常用，是大多数软件包的默认设置。

基于标准化变量的主成分分析详见《机器学习——数据降维——主成分分析（PCA）和奇异值分解（SVD）》第一章第4节。

1.4、因子的个数的选择

心理学，社会科学或者其他行为科学的研究者们可能基于他们的具体领域或经验来指定因子的个数。

如果没有先验信息（例如前面对于各科成绩，我们事先确定一个公共因子——语言能力），如何确定？

$S \approx LL'+\Psi$

$S -(LL'+\Psi)$ 所有元素的平方和 $\leq \sum ^{p}_{k=m+1}\widehat{\lambda }^{2}_{k}$

因此，可以参考主成分分析，评估被忽略的特征值的贡献，即可以使用百分比截点、平均截点和碎石图。主成分个数的选择详见《机器学习——数据降维——主成分分析（PCA）和奇异值分解（SVD）》第一章第5节。

1.5、案例

一个12岁的女孩对她身边的7个人进行9分制评分。评分基于五个维度进行的，分别是“友好(kind)”、“聪明(intelligent)”、“快乐(happy)”、“受人喜爱(likeable)”和“公正( just)”：

1.5.1、计算相关系数矩阵

corr=data.corr()
corr

输出：

根据相关系数的大小，变量可能分为两组：{0, 2, 3} 和 {1, 4}，因此我们期望通过2个因子解释变量间的相关性。

1.5.2、求载荷矩阵

eigVals,eigVects=np.linalg.eig(np.mat(corr)) #协方差矩阵的特征值和特征向量
Loading=np.sqrt(eigVals)*np.array(eigVects) #主成分乘上权重，λ的开平方，注意：此时eigVects中一列是代表着一个特征向量
Loading

输出：

array([[ 0.9694553 ,  0.23114797, -0.07845059,         nan, -0.02392768],
       [ 0.5194021 , -0.80694535,  0.28034884,         nan, -0.02156826],
       [ 0.78451739,  0.58724124,  0.16754548,         nan,  0.10774369],
       [ 0.97087045,  0.20994906, -0.03935088,         nan, -0.10855173],
       [ 0.70396444, -0.66692688, -0.23125766,         nan,  0.07850152]])

我们使用 m=2 个因子，基于相关系数矩阵的主成分法得到的载荷矩阵估计（ m=2 ）、共同度及个体方差如下：

L_=pd.DataFrame(Loading[:,:2],columns=['lj1','lj2'],index=['Kind','Intelligent','Happy','Likeabel','Just'])
L_['共同度']=(Loading[:,:2]*Loading[:,:2]).sum(axis=1)
L_['个体方差']=1-L_['共同度']
L_

输出：

由 $l_{j1}$ 和 $l_{j2}$ 可知这两个公共因子是如何决定原始变量的，如 $Intelligent=0.519l_{j1}-0.807l_{j2}$ ， $Happy=0.785l_{j1}+0.587l_{j2}$ ， $Just=0.704l_{j1}-0.667l_{j2}$ ，三个原始变量在两个公共因子的表现（权重）基本都不低，这样就导致一个问题：公共因子的可解释性很差。在后面的章节我们会介绍一种方法——因子旋转，来解决此类问题。

1.5.3、方差贡献和方差贡献率

eigVals[:2].sum()/eigVals.sum()

输出：

0.9613695373197961

通过两个公共因子能刻画96%的方差，通过以上碎石图能够更加直观的看出：

1.5.4、验算：根据 $\Sigma =L_{p\times m}L_{p\times m}'+\Psi$ ，用两个公共因子能否复原原来的协方差矩阵 $\Sigma$ （或者相关系数矩阵）？

pd.DataFrame(np.dot(Loading[:,:2],(Loading[:,:2].T))+np.diag(1-(Loading[:,:2]*Loading[:,:2]).sum(axis=1)))

输出：

计算过程： $L_{p\times m}L_{p\times m}'+\Psi=\begin{pmatrix} 0.969 & 0.231 \\ 0.519&-0.807 \\ 0.785 & 0.587\\ 0.971&0.210 \\ 0.704 &-0.667 \end{pmatrix}\begin{pmatrix} 0.969 & 0.519 &0.785 &0.971&0.704\\0.231&-0.807 & 0.587 &0.210 &-0.667 \end{pmatrix}+\begin{pmatrix} 0.007& 0& 0& 0& 0 \\ 0& 0.79& 0& 0& 0 \\ 0& 0& 0.40& 0& 0\\0& 0& 0& 0.13& 0 \\ 0& 0& 0& 0& 0.60 \end{pmatrix}= \begin{pmatrix} 1.000& 0.321& 0.897& 0.992& 0.525 \\ 0.321& 1.000& -0.062& 0.338& 0.906 \\ 0.897& -0.062& 1.000& 0.887& 0.158\\0.992& 0.338& 0.887& 1.000& 0.540 \\ 0.525& 0.906& 0.158& 0.540& 1.000 \end{pmatrix}$ 对比原来的相关系数矩阵：

可以看出，很好的复原了原来的相关性系数矩阵，由两个公共因子刻画原始变量是足够的。

2、主因子法

$S \approx LL'+\Psi$

在主成分方法中，我们忽略 $\Psi$ ，对或进行谱分解。

而主因子法/主轴法(principal factor/axis method)，则使用一个初始的估计值 $\widehat{\Psi}^{(0)}$ ，从而对 $S-\widehat{\Psi}^{(0)}$ 或 $R-\widehat{\Psi}^{(0)}$ 使用与主成分方法相同
的操作进行因子分析。

我们可以迭代地更新估计值 $\widehat{\Psi}^{(0)}$ 以及 $S-\widehat{\Psi}^{(0)}$ 或 $R-\widehat{\Psi}^{(0)}$ 的分解，直到收敛。

大致步骤如下：

给定一个初始的估计值 $\widehat{\Psi}^{(0)}$ ，常用的 $\widehat{\Psi}^{(0)}$ 初值：使用时，初值为 $1/diag(diag(S^{-1}))$ ；使用时，初值为 $1/diag(diag(R^{-1}))$ 。
对 $S-\widehat{\Psi}^{(0)}$ 或 $R-\widehat{\Psi}^{(0)}$ 使用与主成分方法相同的操作（谱分解）进行因子分析，得到载荷矩阵 $\widehat{L}$ 的估计。
求 $S-\widehat{L}\widehat{L}'$ 的对角元素，即 $\widehat{\Psi}=diag(\widehat{\Psi} _{1},...,\widehat{\Psi} _{p})$ ，其中 $\widehat{\Psi} _{j}=s _{jj}-\sum _{k=1}^{m}\widehat{l}^{2}_{jk},j=1,...,p$ ，得到新的个体方差估计值 $\widehat{\Psi}^{(1)}$ 。
重复2和3步骤，使得 $\widehat{\Psi}$ 越来越接近真值（收敛）。

缺陷：

因为 $S-\widehat{\Psi}^{(0)}$ 或 $R-\widehat{\Psi}^{(0)}$ 有可能不是正定的，这个方法有可能得到负的特征值，使得结果较难解释（负的方差），如下图。

下表对由主成分和主因子法得到的载荷估计进行比较：

【小贴士】

性质：正定矩阵的特征值都是正的。

3、极大似然法

当样本 $f_{1},...,f_{n}$ 和 $\varepsilon _{1},...,\varepsilon _{n}$ 服从联合正态分布， $y _{1},...,y _{n}$ 也因此独立同分布地服从 $N_{p}(\mu ,\Sigma )$ 分布时（其中 $\Sigma =LL'+\Psi$ ）， $\Sigma$ 的似然函数可以如下表示：

$L(\Sigma )\propto |\Sigma |^{-n/2}exp\left [ -\frac{1}{2} tr\left \{ \Sigma ^{-1} \left ( \sum ^{n}_{i=1}(y_{i}-\overline{y})(y_{i}-\overline{y})' \right )\right \} \right ]$

其中 $\widehat{\mu }_{MLE}=\overline{y}$ （ $\widehat{\mu }_{MLE}$ 表示 $\widehat{\mu }$ 的极大似然估计）已经代入上式，将 $\Sigma =LL'+\Psi$ 代入上式，然后极大化这一个方程，和 $\Psi$ 的极大似然估计就可以通过迭代计算得到。

三、因子旋转

1、因子及载荷的不唯一性

$y-\mu=Lf+\varepsilon$

因子是不可观测的，我们可以将其看作虚拟的，所以其实不唯一的。故以上因子模型等价于 $y-\mu=L^{*}f^{*}+\varepsilon=(LT)(T'f)=L(TT')f$ ，其中 $L^{*}=LT$ ， $f^{*}=T'f$ ，是满足 TT'=T'T=I 的任意正交矩阵。

【小贴士】

任意正交阵乘以其自身转置的结果，为单位阵。

新的载荷矩阵 $L^{*}$ 仍然满足 $\Sigma =LL'+\Psi=L^{*}L^{*}'+\Psi$ 。

新的 $f^{*}$ 仍然满足 $E(f^{*})=0_{m\times 1},COV(f^{*})=I_{m\times m},COV(f^{*},\varepsilon )=0_{p\times m}$ 。

所以因子及其载荷矩阵并不唯一，可以按照任意的正交矩阵提供的方向旋转。这种不唯一性为“因子旋转”提供了理论基础。

我们的最终目标是寻找使因子及载荷结构更简单、解释更清晰的旋转方向。

2、直觉理解

$Y_{1}-\mu _{1}=l_{11}F_{1}+l_{12}F_{2}+..+l_{1m}F_{m}+\varepsilon _{1}$

$Y_{2}-\mu _{2}=l_{21}F_{1}+l_{22}F_{2}+..+l_{2m}F_{m}+\varepsilon _{2}$

$Y_{p}-\mu _{p}=l_{p1}F_{1}+l_{p2}F_{2}+..+l_{pm}F_{m}+\varepsilon _{p}$

我们旋转因子 $f^{*}=T'f$ 就相当于旋转载荷矩阵 $L^{*}=LT$ ，其本质上是在做同一件事情。我们对处理载荷矩阵更加熟悉，载荷矩阵代表原始变量和公共因子之间的关系（协方差/相关系数，取决于使用协方差矩阵还是用相关系数矩阵进行谱分解）。

我们希望得到的是，每个原始变量都由某个因子主要决定（对应载荷数值很大），而与其他因子关系不大（对应载荷绝对值接近0，则此公共因子对原始变量没有作用），即我们想要得到一个结果尽量简单（稀疏）的载荷矩阵，便于我们解释。

从几何的角度， $\widehat{L}$ 的第行的载荷构成了原始变量 $Y_{j}$ 在因子/载荷空间的坐标。所以因子旋转的目标，是让坐标轴靠近尽可能多的点。

回到12岁女孩调查的案例：

上图中的载荷矩阵我们可以投射到如下的载荷（因子）空间当中的5个点。载荷矩阵的每一行可以表示为一个点的坐标。需要注意的是，图中的直角三角形的，根据勾股定理的 $\widehat{h}_{j}^{2}=\widehat{l}^{2}_{j1}+\widehat{l}^{2}_{j2}$ ，即对应的载荷平方和，等于共同度，共同度，即点到原点的距离平方。

总结：因子旋转的目的是得到更加好解释、定义和理解的因子方向（因为因子是虚拟的）。即通过“旋转坐标轴”的方式，尽可能使新的坐标轴靠近更多的点，这边“旋转坐标轴”有两种方式：正交旋转和斜交旋转。

3、正交旋转

正交旋转(Orthogonal rotation)：

原来垂直的坐标轴经过旋转后仍保持垂直。
角度和距离都保持不变。
共同度也不变。
点的相对位置也维持原状。
只有参考系改变了。

3.1、图像法

如果 m=2 ，我们可以通过观察因子载荷坐标系决定如何旋转坐标轴。

我们选择一个角度 $\phi$ 来让坐标轴尽可能靠近更多的点。新的旋转载荷 $(\widehat{l}^{*}_{j1},\widehat{l}^{*}_{j2})$ 可以通过 $\widehat{L}^{*}=\widehat{L}T$ （旋转坐标轴）得到，其中 $T=\begin{pmatrix} cos\phi & -sin\phi\\ sin\phi& cos\phi \end{pmatrix}$ 。

矩阵旋转的角度：顺时针旋转则为负值，逆时针旋转则为正值，更多关于旋转矩阵的内容，详见《线性代数——线性变换——旋转矩阵（泰勒公式、虚数、欧拉公式）》。

回到12岁女孩的例子：取 $\phi=-55^{\circ}$ ， $T=\begin{pmatrix} 0.574 &0.819 \\ -0.819&0.574 \end{pmatrix}$ 。

通过 $\phi=-55^{\circ}$ 的旋转，我们得到了下面的载荷矩阵：

$\widehat{L}^{*}=\widehat{L}T=\begin{pmatrix} 0.969 & 0.231 \\ 0.519&-0.807 \\ 0.785 & 0.587\\ 0.971&0.210 \\ 0.704 &-0.667 \end{pmatrix}\begin{pmatrix} 0.574 &0.819 \\ -0.819&0.574 \end{pmatrix}=\begin{pmatrix}0.367 & 0.927 \\ 0.959& -0.037 \\ -0.031 &0.980\\ 0.385&0.916 \\ 0.950 & 0.194\end{pmatrix}$

旋转后的载荷与原载荷对比：旋转之后的载荷，有一些很大（下图红色框），有一些接近于0（下图蓝色框），无限接近于0的我们可以忽略，这样就会使结果更加好解释、简洁。

旋转后的载荷很容易解释：

第一个因子主要由intelligence和just构成，可解释为人的逻辑思维能力。
第二个因子与变量kind, happy, likeable高度相关，可描述为亲和力。

3.2、最大方差法——最常用的正交因子旋转

$Y_{1}-\mu _{1}=l_{11}F_{1}+l_{12}F_{2}+..+l_{1m}F_{m}+\varepsilon _{1}$

$Y_{2}-\mu _{2}=l_{21}F_{1}+l_{22}F_{2}+..+l_{2m}F_{m}+\varepsilon _{2}$

$Y_{p}-\mu _{p}=l_{p1}F_{1}+l_{p2}F_{2}+..+l_{pm}F_{m}+\varepsilon _{p}$

使用图像法来正交旋转局限于 m=2 ，因为只有二维的空间，才能在图像上呈现出来。

对 m>2 ，有更多解析的方法已经被提出。最常用的为最大方差法（Varimax method），也就是寻找能够最大化载荷矩阵中每一列载荷平方的方差的旋转载荷。

最大方差法（Varimax method）原理：

对，我们也可以用最大方差法，设因子载荷矩阵 $L=\begin{bmatrix} l_{11}& l_{12}\\ l_{21}& l_{22}\\ . &. \\ . &. \\ l_{p1}& l_{p2} \end{bmatrix}$ ，正交阵 $T=\begin{pmatrix} cos\phi & -sin\phi\\ sin\phi& cos\phi \end{pmatrix}$ ，旋转过后的载荷矩阵设为 $L^{*}=\begin{bmatrix} l^{*}_{11}& l^{*}_{12}\\ l^{*}_{21}& l^{*}_{22}\\ . &. \\ . &. \\ l^{*}_{p1}& l^{*}_{p2} \end{bmatrix}$ ，我们要做的是，选定合适的 $\phi$ 可以令旋转过后的载荷矩阵 $L^{*}$ 的列向量的方差之和（ $S_{1}+...+S_{m}=\sum ^{m}_{i=1}\frac{(\sum ^{p}_{j=1}(l^{*}_{ij}-\overline{l^{*}_{i}})^{2})}{p}=\frac{\sum ^{m}_{i=1}\sum ^{p}_{j=1}(l^{*}_{ij}-\overline{l^{*}_{i}})^{2}}{p}$ ，）达到最大。
当载荷矩阵的列时，我们可以逐次对每两个因子进行上述旋转，直至个因子全部配对旋转，即需要旋转 $C_{m}^{2}$ 次，全部旋转完毕算作一轮循环。一般来说，一轮循环是不够的，需要经过多轮旋转，直至方差无法再增大为止。

回到12岁女孩的数据，最大方差法得到的旋转因子载荷与图像法得到的旋转十分接近：

4、斜交旋转

斜交旋转(Oblique rotation)：

不要求轴保持垂直。
因此旋转更加自由。
也更容易让轴靠近更多的点。

回想12岁女孩的例子，如果旋转后的坐标轴允许不再垂直（即斜交旋转），代表 $F_{2}^{*}$ 的轴可以更加靠近第1个和第4个变量对应的点。

正交旋转的好处是：角度和距离都保持不变。共同度也不变。

不像正交旋转中使用正交矩阵，斜交旋转使用一个更一般的非奇异变换矩阵（可逆即可）来得到 $f^{*}=Q'f$ ，那么 $COV(f^{*})=Q'IQ=Q'Q\neq I$ ，因此新的因子之间是相关的，不是正交的。

由于距离和角度不再保持不变， $f^{*}$ 的共同度与的也不同。

斜交旋转的好处是，当不要求坐标轴相互垂直时，旋转后的坐标轴更容易“穿过”多数坐标点。

例：我们收集了25个家庭的第一胎和第二胎婴儿头围测量数据（头长和头宽），如下图，右图为相关系数矩阵：

直观的选择2个公共因子，分别做正交旋转（下图左）和斜交旋转（下图右）：

与最大方差旋转载荷（正交旋转）相比，斜交旋转载荷的结构更加简洁（更加0-1两极化），但解释起来完全一样，即 $F_{1}$ 主要解释第二胎的头长和头宽， $F_{2}$ 主要解释第一胎的头长和头宽。

斜交旋转得到的两个坐标轴的夹角为38°，两个旋转后的因子之间的相关性则为0.79（可以由 Q'Q 的非对角元得到，或由 $cos38^{\circ}=0.79$ 得到）。故我们总结一个结论：两坐标轴间的夹角小于45°，可只选择一个因子。

四、估计因子得分

1、因子得分

有时研究者则希望得到因子得分（Factor score）， $\widehat{f}_{i}=(\widehat{F}_{i1},...,\widehat{F}_{im}),i=1,...,n$ ，即每个因子在不同个体上的取值，从而我们可以

查看不同个体的因子表现情况（比如智力问题）
将因子得分作为“观测值”后续进行其他分析，如分类（降维之后再分类）

2、因子得分的估计——最小二乘法

我们先回顾回归模型： $y=X\omega +\varepsilon$

在标准线性回归中我们需要找到是误差最小的, 即预测的值与真实的值之间的差值。

f(w)=(y-Xw )'(y-Xw )

使用矩阵表示将会是的求解和程序更为简单:

f(w)=(y-Xw)'(y-Xw)=y'y-y'Xw-w'X'y+w'X'Xw

将 f(w) 对求导可得:

$\frac{\partial f(w)}{\partial w}=-X'y-X'y+2X'Xw$

【注意】

矩阵求导性质：

1）、 $\frac{\partial Ax}{\partial x}=A^{T}$

2）、 $\frac{\partial x^{T}A}{\partial x}=A$

使其等于0，便可得到：

0=-X'y-X'y+2X'Xw

求得： $\widehat{w}=(X'X)^{-1}X'y$

类比回归模型，我们在以下的因子模型当中

因子模型： $y_{i}-\mu=Lf_{i}+\varepsilon _{i}$

前面我们已经将载荷矩阵估计出来了，现在目标是估计： $(L'L)^{-1}L'(y-\mu )$ 。

由于 $\varepsilon_{i}$ 的方差不相同，我们可以通过加权最小二乘法来估计： $(L'\Psi ^{-1}L)^{-1}L'\Psi ^{-1}(y-\mu )$ 。

【小贴士】

加权最小二乘法的残差平方和为 $f(f)=(y-\mu -Lf )'\Psi ^{-1}(y-\mu -Lf )$ ，其中 $\Psi= \begin{bmatrix} \Psi _{1} & 0 & ... & 0\\ 0 & \Psi _{1} & ... & 0\\ . & . & & .\\ . &. & &. \\ 0& 0 & 0 & \Psi _{p} \end{bmatrix}$ ，和标准最小二乘法的差别在于： $\Psi^{-1}$ 使得具有较小的方差的样本具有较大权重，经过这般调整，使其方差相同。最后求得 $f=(L'\Psi ^{-1}L)^{-1}L'\Psi ^{-1}(y-\mu )$ ，详细加权最小二乘估计可见《加权最小二乘法与局部加权线性回归》

将已估计的参数 $\widehat{L}$ ， $\widehat{\Psi }$ 和 $\widehat{\mu }=\overline{y}$ （样本均值）代入上式： $\widehat{f}_{i}=(\widehat{L}'\widehat{\Psi} ^{-1}\widehat{L})^{-1}\widehat{L}'\widehat{\Psi} ^{-1}(y_{i}-\overline{y} ),i=1,...,n$ ，最后得到的是一个 $n\times m$ 维的矩阵。

3、因子得分的估计——回归法

基于正态假设：假设和服从

$\begin{pmatrix} f \\ y-\mu \end{pmatrix}\sim N_{m+p}(0,\begin{bmatrix} I&L\\ L'& \Sigma =LL'+\Psi \end{bmatrix})$ ，即服从联合正态分布。

则在给定时的条件分布为正态分布，均值为 $E(f|y)=L'\Sigma^{-1}(y-\mu )$ ，于是我们就可以利用给定的来估计。

我们就可以得到第个因子得分的估计： $\widehat{f}_{i}=\widehat{L}'(\widehat{L}\widehat{L}'+\widehat{\Psi})^{-1}(y_{i}-\overline{y} ),i=1,...,n$ ，或者 $\widehat{f}_{i}=\widehat{L}'S^{-1}(y_{i}-\overline{y} ),i=1,...,n$

【小贴士】子向量的条件正态性：

假设 $y\sim N_{p}(\mu ,\Sigma )$ ，现在对、 $\mu$ 、 $\Sigma$ 以第个元素为界进行分割如下： $y=\begin{pmatrix} y_{1}\\ y_{2} \end{pmatrix}$ ， $\mu=\begin{pmatrix} \mu _{1}\\ \mu _{2} \end{pmatrix}$ ， $\Sigma =\begin{pmatrix} \Sigma _{11}&\Sigma _{12}\\ \Sigma _{21}&\Sigma _{22}\end{pmatrix}$ ，这里 $\Sigma _{11}$ 是 $r\times r$ 的，并且 $\left | \Sigma _{22} \right |>0$ 。

于是在给定 $y_{2}$ 时， $y_{1}$ 的条件分布仍然是多元正态： $y_{1}|y_{2}\sim N_{r}(\mu _{1}+\Sigma _{12}\Sigma^{-1} _{22}(y_{2}-\mu _{2}),\Sigma _{11}-\Sigma _{12}\Sigma ^{-1}_{22}\Sigma _{21})$ 。

注意：

$E(y_{1}|y_{2})$ 是关于 $y_{2}$ 的线性方程，同时 $COV(y_{1}|y_{2})$ 不依赖 $y_{2}$ 。

五、spss示例

例：下表记录了52位学生6门功课的考试分数的部分数据，案例数据集下载：案例数据集《多元统计分析-数据降维-因子分析（FA）》， $Y_{1}$ 到 $Y_{6}$ 依次表示数学、物理、化学、语文、历史、英语成绩，将数据导入spss软件当中，部分数据如下：

1、适用性检验

因子分析的前提是原始维度具有相关性。因此，适用性检验就是检验原始维度间是否具有相关性，如果没有相关性，则不适合做因子分析。具体操作步骤如下

选择菜单：分析——>降维——>因子分析
设置变量：在弹出的“因子分析”对话框中，将 $Y_{1}$ 到 $Y_{6}$ 变量选入“变量"对话框中。
设置选项：单击”描述“按钮，选择”原始分析结果“和”KMO和Bartlett的球形度检验（K)"。

于是，输出KMO和Bartlett的球形度检验结果，见下表。KMO=0.804＞0.7，Bartlett 的球形度检验Sig.=0.000＜0.05，通过因子分析适用性检验，说明本案例适合做因子分析。

球形检验原理：球形检验主要是用于检验数据的分布，以及各个变量间的独立情况。按照理想情况，如果我们有一个变量，那么所有的数据都在一条线上。如果有两个完全独立的变量，则所有的数据在两条垂直的线上。如果有三条完全独立的变量，则所有的数据在三条相互垂直的线上。如果有个变量，那所有的数据就会在条相互垂直的线上，在每个变量取值范围大致相等的情况下（常见于各种调查问卷的题目），所有的数据分布就像在一个球形体里面。想象一下万剑穿心的情形，大抵就是那个样子。如果不对数据分布进行球形检验，在做因素分析的时候就会违背因素分析的假设——各个变量在一定程度上相互独立。在spss中,KMO检验用于检查变量间的相关性和偏相关性，取值在0~1之前。KMO统计量越接近于1，变量间的相关性越强，偏相关性越弱，因子分析的效果越好。实际分析中，KMO统计量在0.7以上时效果比较好；当KMO统计量在0.5以下，此时不适合应用因子分析法，应考虑重新设计变量结构或者采用其他统计分析方法。

2、因子提取（主成分法）

按上面的适用性检验操作，还会输出公因子方差和解释的总方差，用于因子提取。首先，从公因子方差（下表）可知，提取的因子对这6个变量的解释度均超过70%，部分超过80%，说明所提取的因子对原始维度具有一定的解释力。

其次，从解释的总方差（见下表）可知，按特征值＞1的标准（下表中“元件”即为“因子”，初始特征值中的“总计”即为各因子的特征值），应提取前2个因子，从下表可以看出前2个因子的累积方差贡献率为82.873%，表明这些因子能够解释总体信息量的82.873%。

3、因子旋转

按上方操作，还会输出成分矩阵（将下表），用于判断是否需要因子旋转。

上表中的数值称为因子载荷，表示因子对6个维度信息的解释程度。

从上表可知，因子1解释了“Y5”这个维度81.0%的信息，而因子2解释了该维度46.6%的信息，81.0%与46.6%数值相近，表明因子1和因子2都具有该维度的特征，具有相关性。同理，该表还显示出因子1和因子2都具有“Y6”维度的特征。如前所述，因子分析的目的就是剔除相关性，使各个因子具有差异化的特征，而目前的成分矩阵没有达到既定效果，因此需要进行因子旋转。

为什么因子旋转能使各因子具有差异化的特征呢？以因子1和因子2为例，这两个因子可以被分别看成是一个平面直角坐标系的横纵坐标轴，而“Y5”维度就是该坐标系中的一个点，该点在横纵坐标轴上的投影（即该点的横纵坐标值）就是该维度在因子1和因子2上的因子载荷。目前该点的横纵坐标值相近（81.0%和46.6%），但如果调整坐标轴的夹角或该坐标系与水平面的夹角（即因子旋转），该点的横纵坐标值必然跟着改变，当调整到该点的横纵坐标值相差悬殊，例如横坐标值很大、纵坐标值很小时，则因子1具有该维度的特征，因子2不具有（几乎）该维度的特征，从而实现了两个因子特征的差异化。

如何进行因子旋转呢？

在上面操作的基础上，按下图所示分别设置 “ 旋转”和“选项”选项。其中“大方差法”是常用的“旋转”方法，若效果不好，则可尝试其他方法；在“选项”中可以设置不显示绝对值过小的因子载荷（本例设置为0.4），这样更容易观察因子特征。

按上图所示进行设置后，得到旋转成分矩阵（见下表），该矩阵中每个维度上仅有1个因子具有较大（＞40%）的因子载荷，即2个因子在6个维度上的特征被明显区隔开，2个因子具有差异化的特征。

4、因子命名

旋转后的成分矩阵有效区分了2个因子的维度特征，为直观呈现各因子的维度特征，方便对因子命名，这里根据旋转成分矩阵（见上图）整理出下表。

	因子1	因子2
维度特征（因子载荷值）	语文（0.878）	数学（0.838）
	历史（0.918）	物理（0.783）
	英语（0.926）	化学（0.896）

旋转过后的公共因子，第一个公共因子主要决定后面三个变量（ $Y_{4}$ ， $Y_{5}$ ， $Y_{6}$ ，即语文、英语、历史），对应的可以概括为“文科因子” ，第二个公共因子主要决定前面三个变量（ $Y_{1}$ ， $Y_{2}$ ， $Y_{3}$ ，即数学、物理、化学），对应的可以概括为“理科因子” 。即我们通过因子分析，可以将我们的变量分成两类，具体分类如下：

5、计算因子得分

因子得分指每条记录（记每一个学生）在所提取因子上的得分，得分越高，表明该条记录（即该学生）越具有该因子的特征。因子得分为后续设置因子变量奠定了基础。计算因子得分的具体操作如下：

在上面操作的基础上，设置“得分”选项，在“因子分析：因子得分”对话框中选择“保存为变量”、“回归”和“显示因子得分系数矩阵”，单击“继续”按钮，返回到“因子分析”对话框中，单击“确定”按钮完成设置（见下图）。

由于选择了“显示因子得分矩阵系数”选项，因此在 SPSS的输出窗口中输出了成分得分系数矩阵。

由于选择了“保存为变量”选项，因此在SPSS的数据视图中生成了因子得分，作为新变量保存在数据文件中（见下图）。这意味着每一个学生都有了2个因子的得分，这为下一步确定学生的因子类型奠定了基础。

6、设置因子变量

从上图可知，每个学生（即每条记录）都有2个因子的得分，根据得分大小，可以判断出这些学生的因子类型。例如，第1个学生（即第1条记录）的因子1至因子2的得分分别为0.66294和 -0.68470，显然因子1得分高，因此，第1个客户的因子变量值为1，表示该客户属于因子1的类型，即为“偏文科”学生。

那么如何设置因子变量呢？具体操作步骤如下：

选择菜单：转换——>计算变量
设置选项：“目标变量”选项，输入“因子类别”。变量值“1”的设置方法（“如果”选项：选择“如果个案满足条件则包括（F)"；在被激活的条件框中设置公式：FAC1_1=MAX(FAC1_1,FAC2_1)。（MAX是最大值函数。该公式表示，如果某条记录满足因子1得分即FAC1_1最大的条件。公式中的FAC1_1,FAC2_1从变量框选入，MAX从函数组中的”统计量“选入。）单击”继续“按钮，回到”计算变量“对话框。”数字表达式“选项：输入”1“（1表示因子1））；变量值“2”的设置方法（“如果”选项：选择“如果个案满足条件则包括（F)"；在被激活的条件框中设置公式：FAC2_1=MAX(FAC1_1,FAC2_1)。单击”继续“按钮，回到”计算变量“对话框。”数字表达式“选项：输入”2“（2表示因子2））。单击”确定“按钮，完成设置。

于是得到输出结果：变量名为“因子类别”，变量值为自然数1～2，其中1表示“文科因子”，2表示“理科因子”。如下：

在变量视图中设置因子变量的值标签，这样就完成了“因子类别”变量的设置与编码。

接下来就可以将使用“因子类别”变量替代反映原始数据的4个维度作为例如聚类分析的细分维度之一了。

六、python示例

1、python适用于因子分析的包

sklearn当中也有因子分析的库sklearn.decomposition.FactorAnalysis，但是是采用极大似然估计的方式计算结果，不能旋转。不能旋转的因子分析对原始维度缺少一定的解释力，并且因子间可能存在一定的相关性，达不到因子分析的既定效果。

故我们可以用factor_analyzer.FactorAnalyzer:既可做因子分析也能做因子的旋转，格式如下：FactorAnalyzer(rotation=None, n_factors=n, method='principal')

2、案例解析——因子分析

案例：回到第5节中考试成绩的例子。52位学生6门功课的考试分数： $Y_{1}$ 到 $Y_{6}$ 依次表示数学、物理、化学、语文、历史、英语成绩。

同案例的主成分分析（PCA）过程可见《机器学习——数据降维——主成分分析（PCA）和奇异值分解（SVD）》

目的：对比于主成分分析（PCA），因子分析可以更好的解释这两科的差异性，可以由“文科”和“理科”来刻画，或者是否有其他更合理解释角度？

2.1、计算相关系数矩阵

corr=data.corr()
corr

输出：

从相关性系数矩阵上看，前三个变量（ $Y_{1}$ ， $Y_{2}$ ， $Y_{3}$ ）更加相关，后三个变量（ $Y_{4}$ ， $Y_{5}$ ， $Y_{6}$ ）更加相关，我们暂定用 m=2 个公共因子来刻画，待估计的因子模型初步定为： $Y_{j}-\mu _{j}=l_{j1}F_{1}+l_{j2}F_{2}+\varepsilon _{j},j=1,...,6$ 。

2.2、估计载荷矩阵

2.2.1、极大似然法

from factor_analyzer import FactorAnalyzer
fa = FactorAnalyzer(rotation=None,  #是否旋转因子
                    n_factors=2, #取两个因子
                    method='ml'  #极大似然估计
                   )
data_=fa.fit_transform(data)
fa.loadings_  #载荷矩阵

输出：

array([[-0.67553954,  0.5618195 ],
       [-0.59919972,  0.42715296],
       [-0.48652622,  0.65618063],
       [ 0.91689514,  0.1039815 ],
       [ 0.85579329,  0.23917605],
       [ 0.88275542,  0.26616829]])

取2个公共因子的载荷矩阵（极大似然估计）如上。

fa.get_factor_variance()

输出：

(array([3.40443772, 1.06753672]),
 array([0.56740629, 0.17792279]),
 array([0.56740629, 0.74532907]))

第一行表示前两个因子对总体方差的贡献，第二行表示方差贡献的占比，第三行表示方差贡献的累积占比，我们可以看出，前两个因子的方差贡献率在74.53%。

2.2.2、主成分法

from factor_analyzer import FactorAnalyzer
fa = FactorAnalyzer(rotation=None,  #是否旋转因子
                    n_factors=2, #取两个因子
                    method='principal'  #主成分法
                   )
data_=fa.fit_transform(data)
fa.loadings_  #载荷矩阵

输出：

array([[-0.79365794,  0.42244881],
       [-0.73419111,  0.40079553],
       [-0.63972828,  0.63215697],
       [ 0.88829277,  0.31294912],
       [ 0.81009848,  0.46605162],
       [ 0.82849201,  0.45674661]])

取2个公共因子的载荷矩阵（主成分法）如上，注意：以上是未经过旋转的载荷矩阵。

fa.get_factor_variance()

输出：

(array([3.70990438, 1.26248122]),
 array([0.6183174 , 0.21041354]),
 array([0.6183174 , 0.82873093]))

第一行表示前两个因子对总体方差的贡献，第二行表示方差贡献的占比，第三行表示方差贡献的累积占比，我们可以看出，前两个因子的方差贡献率在82.87%。比上面的极大释然估计的方差贡献率（74.53%）要高的，极大释然估计是基于正态分布的，也许原始样本分布并不是很接近“正态分布”。所以我们后续的结果倾向于继续使用主成分法来讨论。

2.3、计算旋转因子载荷——基于主成分法

from factor_analyzer import FactorAnalyzer
fa = FactorAnalyzer(rotation='varimax',  #最大方差法因子旋转
                    n_factors=2, #取两个因子
                    method='principal'  #主成分法
                   )
data_=fa.fit_transform(data)
fa.loadings_  #载荷矩阵

输出：

取2个公共因子的载荷矩阵（主成分法）如上。

fa.get_factor_variance()

输出：

(array([2.66060044, 2.31178516]),
 array([0.44343341, 0.38529753]),
 array([0.44343341, 0.82873093]))

第一行表示前两个因子对总体方差的贡献，第二行表示方差贡献的占比，第三行表示方差贡献的累积占比，我们可以看出，前两个因子的方差贡献率在82.87%。与未旋转的方差贡献率相比没有变化。

但是两个因子的载荷与未旋转之前相比差了很多：

两个公共因子，我们可以画图展现出来更加直观：

import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

fig=plt.figure(figsize=(10,4)) #表示绘制图形的画板尺寸为6*4.5；
ax=fig.add_subplot(1,2,1)

loading=fa_1.loadings_ #未旋转的因子载荷

plt.scatter(loading[:,[0]],loading[:,[1]])
plt.xlabel('第一个公共因子')
plt.ylabel('第二个公共因子')
for i in range(len(loading)):
    plt.text(loading[i,0],loading[i,1],'Y%s'%(i+1),ha='left',va='top',fontsize=12,rotation=0,alpha=50) #columns.index(i)返回下标,#ha=‘right'表示点在注释右边，va='bottom'表示点在注释底部，alpha表示透明程度
plt.title('未旋转的因子载荷')
    
ax=fig.add_subplot(1,2,2)
loading=fa_2.loadings_ #旋转的因子载荷

plt.scatter(loading[:,[0]],loading[:,[1]])
plt.xlabel('第一个公共因子')
plt.ylabel('第二个公共因子')
for i in range(len(loading)):
    plt.text(loading[i,0],loading[i,1],'Y%s'%(i+1),ha='left',va='top',fontsize=12,rotation=0,alpha=50) #columns.index(i)返回下标,#ha=‘right'表示点在注释右边，va='bottom'表示点在注释底部，alpha表示透明程度
plt.title('旋转过后的因子载荷')

输出：

从图中我们可以看到，旋转过后的因子载荷矩阵，更加两级分化（更靠近轴和轴），解释起来将更加清晰。

2.4、计算样本因子得分——基于主成分法

即计算所有的学生在这两个因子上面的表现情况：

from factor_analyzer import FactorAnalyzer
fa_2 = FactorAnalyzer(rotation='varimax',  #最大方差法因子旋转
                    n_factors=2, #取两个因子
                    method='principal'  #主成分法
                   )
data_=fa_2.fit_transform(data)
data_

输出：

array([[ 0.66680653, -0.69388904],
       [-1.08617442, -0.15723761],
       [-1.61684799, -1.90160396],
       [-0.72920389,  0.15382713],
       [-1.76907501, -1.13891434],
       [-1.85801275,  0.64436162],
       [-1.31915844,  1.26530921],
       [-1.74117705, -1.23097618],
       [ 0.95725666, -0.6764035 ],
       [-0.84648864,  0.35102727],
       [-0.90511521, -0.1115654 ],
       [ 0.57900914, -0.24869373],
       [-0.30651906,  1.18778613],
       [ 0.58064712,  0.56696106],
       [-0.77862214, -0.8278795 ],
       [-0.60582854,  0.58630962],
       [ 0.06228757,  0.58790371],
       [ 1.26041827,  0.09990526],
       [ 1.6165137 , -0.36917334],
       [ 0.89554779,  0.99821778],
       [ 1.37610116,  0.20272853],
       [-0.05046723,  1.42552427],
       [ 0.12385714,  1.64671775],
       [-0.10702467,  0.24068394],
       [-0.73672797,  1.32010354],
       [ 0.90379062,  1.71032496],
       [-0.11252153, -0.69081664],
       [-0.1599815 , -1.738235  ],
       [ 0.36325364, -0.35529494],
       [ 2.16545793, -1.06019214],
       [ 0.5120452 , -1.75445471],
       [ 0.24420505, -1.47332172],
       [ 1.41023729,  1.35395876],
       [-0.79261961,  0.83837722],
       [-0.03566869,  0.40380043],
       [-0.79881006, -0.63984564],
       [-0.13342206, -1.64587618],
       [-1.21747732,  0.92526294],
       [ 0.35162225, -0.20707309],
       [ 0.02153877,  1.42305047],
       [ 0.78520054, -0.03781809],
       [-0.34213278, -0.17680972],
       [ 0.05145438,  1.1367464 ],
       [-1.08960807, -0.70504028],
       [-0.9456918 ,  2.06015672],
       [-0.56815626, -1.00460889],
       [ 0.58396812,  0.36319028],
       [-0.1713918 , -1.01788887],
       [ 2.3945996 , -0.90107041],
       [ 0.09178151, -0.84650009],
       [ 1.61876748,  0.2735249 ],
       [ 1.20755702, -0.1545769 ]])

我们已知：第一个公共因子可以概括为“文科因子” ，第二个公共因子可以概括为“理科因子”。所以，例如第一个学生：则表现为文科好，理科差的情况，第二个学生则是文理都相对较差。注意：因子都是正向去决定变量的，正值则代表在因子的表现上较好，负值则代表在因子的表现上较差。

我们可以将所有学生在因子上的表现画在图上：

fig=plt.figure(figsize=(10,10)) #表示绘制图形的画板尺寸为6*4.5；
plt.scatter(data_[:,[0]],data_[:,[1]])
plt.xlabel('文科公共因子')
plt.ylabel('理科公共因子')
#画横纵坐标参考线
plt.hlines(y=0,xmin=-2,xmax=2.5,linestyles='dashdot')
plt.vlines(x=0,ymin=-2,ymax=2.5,linestyles='dashdot')

#画样本点标签
for i in range(len(data_)):
    plt.text(data_[i,0],data_[i,1],'Y%s'%(i+1),ha='left',va='top',fontsize=10,rotation=0,alpha=50) #columns.index(i)返回下标,#ha=‘right'表示点在注释右边，va='bottom'表示点在注释底部，alpha表示透明程度

#画原始参考坐标系
for i in range(len(loading)):
    plt.annotate('Y%s'%(i+1)   #箭头的文字
                 ,xy=(0,0)  #箭头的终点
                 ,xytext=(loading[i,0]*2,loading[i,1]*2)  #起点
                 ,color='r'  #文字的类型
                 ,arrowprops=dict(arrowstyle="<-",color='red',connectionstyle="arc3")  #箭头的样式
                ) 

#打标题
plt.title('因子得分')

输出：

我们可以将各个同学在因子空间上的得分分成四个象限：第一象限即可以表示为“学霸”类型，第二和第四象限可以表示为“偏科”类型，第三象限即为要“好好努力”类型。

同样我们可以将原始变量空间表现在我们的因子载荷空间上（图中的红色坐标系）。我们发现： $Y_{1}$ ， $Y_{2}$ ， $Y_{3}$ 主要由第二个因子决定， $Y_{4}$ ， $Y_{5}$ ， $Y_{6}$ 主要由第一个因子决定。

七、因子分析与主成分分析的比较

区别：

因子分析通常指是一种模型；而主成分分析不涉及模型，是寻找综合指标刻画数据差异性的降维方法。
因子分析旨在通过公共因子解释原始变量间的相关性；而主成分分析旨在通过综合指标解释个体之间的差异性。
在因子分析中，原变量表示为因子的线性组合；而在主成分分析中，主成分则是原变量的线性组合。
因子分析的估计不唯一，可以通过因子旋转寻找最容易解释的公共因子；而主成分的构造是唯一的。
因子分析需要构造因子模型，着重要求新变量具有实际的意义，能解释原始变量间的内在结构。主成分分析仅仅是变量变换，强调新变量贡献了多大比例的方差，不关心新变量是否有明确的实际意义。

联系：

两者都是降维和信息浓缩的方法。
生成的新变量均代表了原始变量的大部分信息且互相独立，都可以用于后续的回归分析、判别分析、聚类分析等等。

八、因子分析总结

因子分析可以很好地满足维度分析的需求；
对于没有业务经验的数据分析人员来讲，是通过观察每个原始变量在因子上的权重绝对值来给因子取名成的。而对于业务知识丰富的数据分析人员，已经对变量的分类有一个预判，并通过进行不同的变量转换（标准化）方式和旋转方式使得预判别为同一组的原始变量在共同的因子上权重绝对值最大化。所以因子分析的要点在于选择变量转换方式。
因子分析作为维度分析的手段，是构造合理的聚类模型的有效步骤。