代码收藏家技术教程 2022-09-21

数学建模常用算法—因子分析

今天跟大家聊一下，在面对评价类题型时候，比较常用的算法—因子分析。

模型的含义

因子分析的基本目的就是用少数几个因子去描述许多指标或因素之间的联系，即将相关比较密切的几个变量归在同一类中，每一类变量就成为一个因子，以较少的几个因子反映原资料的大部分信息。运用这种研究技术，我们可以方便地找出影响消费者购买、消费以及满意度的主要因素是哪些，以及它们的影响力。运用这种研究技术，我们还可以为市场细分做前期分析。

模型的建立与求解

1.相关性检验

因子分析法的使用前提是原始变量存在相关性。这很好理解，我要把一些原始表量归纳为一类，就需要它们代表同一种属性。万一它们互不相关、自成一派，那根本没办法归纳，无法降维，也就无法进行因子分析。

在stata中进行一般使用巴特利特球形检验和KMO检验。

巴特利特球形检验用于检验相关阵是否是单位阵，即各变量是否独立。它是以变量的相关系数矩阵为出发点，零假设：相关系数矩阵是一个单位阵。如果巴特利球形检验的统计计量数值较大，且对应的相伴概率值小于用户给定的显著性水平，则应该拒绝零假设；反之，则不能拒绝零假设，认为相关系数矩阵可能是一个单位阵，不适合做因子分析。若假设不能被否定，则说明这些变量间可能各自独立提供一些信息，缺少公因子。一般实证中，p值小于0.05就可以进行下去。

KMO统计值是通过比较各变量间简单相关系数和偏相关系数的大小判断变量间的相关性。相关性强时，偏相关系数远小于简单相关系数，KMO值接近1。一般来说，KMO>0.9非常适合因子分析；0.8＜KMO＜0.9适合；0.7以上尚可，0.6时效果很差，0.5以下不适宜作因子分析。不过实证中能达到0.6真的已经非常非常可了，所以KMO大于0.5就可以试着做下去。

在stata14.0中，输入命令factortest x1 x2 x2….xn，就可以得到两个检验值。如下图，方框内为命令语句，圆圈内分别问巴特利特球形检验p值和KMO检验值。

（合格的KMO值）

在做适用性检验时，我出现过KMO统计值小于0.5的情况，如下图。可能有2个原因：一是原始数据的异常值对结果有影响；二是原始变量不适合做因子分析。（当然原因二是我极度十分非常不愿意接受的，因为这意味着我需要另外再去找变量。。。）

（不合格的KMO值）

1.针对异常值，只需要对原始数据进行缩尾处理即可。

（1）安装winsor2命令：ssc install winsor2 replace

（2）执行winsor2命令：winsor2 x1 x2…xn, replace cuts(1 99)

这个命令的意思是，对变量进行上下1%的替换缩尾处理，也就是说，如果一个样本变量值大于变量的99分位数，则该样本的值被强制指定为99分位数的值；如果一个样本变量值小于变量的1分位数，则该样本的值被强制指定为1分位数的值。

2.针对原始变量不适合做因子分析，需要检验相各变量间的关性，然后将相关性很低的变量剔除。命令为pwcorr x1 x…xn。像这个里面的g31 g41和g51，跟其他变量之间的相关性都非常低，所以应该剔除，再去寻找其他变量。

（相关系数矩阵）

2.构造因子变量

我们应该选择哪些因子来解释原始变量？解释力度如何？这种力度可以接受吗？这些问题都是在构造因子变量时需要解决的。不过归根结底都是关于累计贡献率——怎么看累计贡献率 & 累计贡献率够不够。

在看累计贡献率时，stata是默认取特征根大于1的因子作为主成分，这表示主成分的解释力度大于直接使用1个原始变量。别忘了，因子分析的本质就是降维，如果提取的主成分特征根小于1，就不能达到化繁为简的目的，那还不如直接用原始变量。

一般认为，累计贡献率达到80%是比较好的，说明主成分包含了原始变量80%以上的信息。不过在实证中，很多人似乎都在65%左右徘徊，如果可以自圆其说的话，其实也可以接受。

输入命令factor x1 x2…xn pcf，就可以得到主成分分析的结果，如下图。4表示提取了4个主成分，因为只有前4个因子的特征根大于1，最后一列表示累计方差贡献率。这表示，提取的4个主成分可以解释原始变量中61.29%的信息。

（提取主成分）

在做因子分析时，我遇到了2个问题：一是提取出的主成分比我预想的多一个，二是累计贡献率不够。

1.主成分多一个。如下图，其他年份都只有4个主成分，而这一年的却有5个主成分。在去掉最后一个因子，累计贡献率不会太低的情况下，可以指定只提取4个主成分。命令为factor x1 x2…xn， factors(4) pcf，在之后的分析中就可以只用4个因子。

（5个主成分）

（4个主成分）

2.累计贡献率不够。我之前做过一版原始数据的因子分析，提取出主成分以后累计贡献率不足55%。这时候就需要看一下到底是哪些原始变量在捣乱，可以通过因子载荷矩阵来看，如下图。一般来说，反常因子的值要小于0.6，像下图中的g6m1达到了0.87，可以将它剔除，这样整体的累计贡献率会提高。

（剔除反常因子）

3.因子旋转

本来做完主成分分析，我们就可以知道哪些因子解释了哪些原始变量，并且知道解释力度是多少。但还是不够清晰，我们需要让因子走向极端。比如，f1对g1的解释力度是0.6，对g2的解释力度是0.5，我们怎么知道f1到底更侧重哪个属性呢？因子旋转要做的事情就是明确这种态度，如果让F1对应0.9的g1、对应0.01的g2，这样就会让结论更清晰。输入命令rotate，各因子对应的的原始变量就会非常清晰了。

如下图，第一个是旋转前的因子载荷矩阵，f1似乎对g30，g40，g50的解释力度都还不错，但是它具体代表什么属性恨不清晰；第二张图是旋转以后的因子载荷矩阵，我们可以很清晰地看出来f1在解释g30和g40。

（旋转前因子载荷矩阵）

（旋转后因子载荷矩阵）

4.计算因子得分

在进行因子旋转以后，输入命令predict f1 f2 f3，就可以得到各因子的得分。

（因子得分预测值）

这时候要构建综合得分函数，还需要知道各因子的权重。非常简单，以各因子的方差贡献率分别除以累计贡献率，得出各因子的权重。下图是stata旋转结果和用Excel制作的权数。输入命令

gen Fm1=0.3610*f1m1+0.3535*f2m1+0.2856*f3m1即可计算每个样本的综合得分。

（各因子累计贡献率）

（各因子对应权数）

5.实证检验

因为我这个是配对样本，看了下很多相关论文用的是均值t检验，所以我很开心地输入了ttest F0 ==F1，检验结果如下图，一把辛酸泪。。。我在论坛上看到一个帖子也遇到了同我一样的问题：配对样本均值t检验的结果p值全部为1。

（配对样本t检验）

后来我查了很多资料发现，配对样本均值t检验的使用前提很严格，需要样本呈正态分布。所以我对自己的样本数据做了下正态性检验，输入命令sktest Fm1 F0 F1 F2 F3，结果如下图：p值全部小于0.05，不是正态分布。后来我对样本做了一些处理：倒数、平方、立方等等，颠七倒八的整了一堆，但都没办法实现正态分布。这种情况下不能直接用参数检验，所以得用非参数检验。