Python之因子分析详细步骤

1.数学原理

1.1数学模型

1.2正交因子模型假设

注意：下面的推导都是基于这一假设。因此，这里的模型都是属于正交因子模型。

1.3正交因子模型的协方差结构

1.4各类方差贡献的介绍

在1.3正交因子模型的协方差结构中，我们介绍了“方差贡献”，下面给出关于“方差贡献”更详细的介绍。

1.5因子表示具有不唯一性

利用此性质，我们可以选择不同的正交矩阵 $T_{m\times m}$ ，做以下操作：

从而获得容易解释的载荷矩阵和因子。

这一操作被称为“因子旋转”。

因子旋转的类型

因子载荷的正交变换和伴随的因子正交变换为因子正交旋转。

进一步地，可以修改正交因子模型，允许共同因子之间相关，称这样的变换为因子斜交旋转(bolique rotation)。如果中的矩阵不是正交阵，则中的各个公共因子可能相关，但公共因子的方差仍要求等于1。

因子分析中常用的正交旋转和斜交旋转方法：

varimax旋转：正交旋转方法，要求每个因子仅有少数绝对值较大的载荷，多数载荷接近0。通过迭代最小化载荷系数的一个二次函数实现。结果使得每个因子仅与少数原始变量有强相关，而与其余原始变量近似不相关。
quartimax旋转：正交旋转方法，要求每个原始变量仅与一个公共因子强相关，而与其余公共因子近似不相关。不如varimax接受程度高。
oblimin旋转：斜交旋转方法，追求载荷阵中0尽可能多，设置公共因子之间的相关系数在较小值。
promax旋转：斜交旋转方法，追求载荷阵中0尽可能多，方法是取正交旋转得到的载荷阵元素的幂次。要求幂次尽可能低，公共因子间的相关性尽可能低。

1.6计算因子得分

在1.1数学模型中，我们建立了如下模型：

因子分析中，首先要得到因子载荷阵，对因子进行解释。

其次，可以从数据中估计公共因子的取值（注意公共因子在模型中是不可观测的，或者说是未知的）。公共因子的取值的估计称为因子得分。

注意：在模型中，m<p，因此无法得到公共因子的精确值（即：因子得分），只能估计。

常使用加权最小二乘法估计因子得分。

加权最小二乘法估计因子得分

得到因子得分后，就可以利用因子得分进行其他分析。

2.Python代码实现

关键的库：factor_analyzer、scipy

下面用这个例子来说明Python中实现因子分析的步骤

2.0导入库和数据

import numpy as np
from numpy.linalg import inv
import pandas as pd
from scipy.stats import zscore
from  factor_analyzer import FactorAnalyzer as FA
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity  # 用于Bartlett's球状检验
from factor_analyzer.factor_analyzer import calculate_kmo  # 用于KMO检验# 导入数据
data=np.loadtxt("data12_5_1.txt")
print(data)

2.1数据标准化

zscore()

# 数据标准化
data_norm=zscore(data,ddof=1)

对应的数学公式：

2.2Bartlett's球状检验和KMO检验

Bartlett's球状检验:计算巴特利特P值
若P值<0.05,则说明变量间有相关性,则可因子分析
KMO检验:计算KMO值
KMO值介于0~1之间，越接近1，说明变量间相关性越强，偏相关性越弱，因子分析的效果越好。
若KMO值>0.6,说明变量间有相关性,则可因子分析。

# Bartlett's球状检验:计算巴特利特P值
chi_square_value,p_value=calculate_bartlett_sphericity(data_norm)
print("Bartlett's球状检验:\n","卡方值:",chi_square_value,"\nP值(若<0.05,则说明变量间有相关性,则可因子分析):",p_value)# KMO检验:计算KMO值
kmo_all,kmo_model=calculate_kmo(data_norm)
print("\nKMO(KMO值介于0~1之间，越接近1，说明变量间相关性越强，偏相关性越弱，因子分析的效果越好。\n""若KMO值>0.6,说明变量间有相关性,则可因子分析):",kmo_model)

2.3求相关系数矩阵

np.corrcoef()

# 求相关系数矩阵
r=np.corrcoef(data_norm.T)

对应的数学公式及分析：

2.4求相关系数矩阵的特征值、特征向量，并排序

法一：np.linalg.eig(r)

注意：此法得到的特征值没有按从大到小的顺序排列，需要自行排序
排序的目的：后面需要利用特征值来求累积贡献率，然后利用累积贡献率 $\geqslant$ 85%来挑选前n个公共因子。只有这样，才能保证挑选出来公共因子，既能代表较多的信息，又数量较少。

# 求相关系数矩阵的特征值、特征向量
eigval,eigvec=np.linalg.eig(r)
print("特征值:\n",eigval)
print("特征向量:\n",eigvec)# 对特征值、特征向量按从大到小排序
address_sort=np.argsort(-eigval)
result_sort=np.zeros(len(eigval))
result_sort[address_sort]=np.arange(0,len(eigval))
result_sort=[int(i) for i in result_sort]
print("特征值从大到小的顺序",result_sort)eigval=eigval[result_sort]
eigvec=eigvec[:,result_sort]
print("排序后的特征值:\n",eigval)
print("排序后的特征向量:\n",eigvec)

法二：建立一个不旋转的因子分析模型，再利用FA.get_eigenvalues()

注意：此法会得到2组特征值，但只有第一组是需要的，本人不太清楚是为什么

优点：得到的特征值已经按从大到小的顺序排列

fa0=FA(n_factors=len(r),rotation=None)
fa0.fit(data_norm)
val1,val2=fa0.get_eigenvalues()  # 注意：这里会得到2组特征值，但只有第一组是需要的
print("fa0的特征值:\n",val1)

2.5求累积贡献率

contibute_ratio=eigval/sum(eigval)
print("贡献率:\n",contibute_ratio)
print("累积贡献率:\n",np.cumsum(contibute_ratio))

对应的数学公式及分析：

2.6选公共因子数量

选取累积贡献率达到85%以上的前n个公共因子

'''由于前3个公共因子的累计贡献率已经超过0.85，
因此认为用3个公共因子 就能较好地进行评价'''
n=3

2.7构建并训练模型

FA(n_factors,rotation)

# 构建模型
fa=FA(n_factors=n,rotation="varimax")
'''FA
参数解读：
n_factors:公共因子数
rotation:因子旋转方式
有以下的选择：
varimax (orthogonal rotation):正交旋转方法，方差最大化
promax (oblique rotation):斜交旋转方法，追求载荷阵中0尽可能多，方法是取正交旋转得到的载荷阵元素的幂次。要求幂次尽可能低，公共因子间的相关性尽可能低。
oblimin (oblique rotation)斜交旋转方法，追求载荷阵中0尽可能多，设置公共因子之间的相关系数在较小值。
oblimax (orthogonal rotation)
quartimin (oblique rotation)
quartimax (orthogonal rotation)正交旋转方法，要求每个原始变量仅与一个公共因子强相关，而与其余公共因子近似不相关。 
equamax (orthogonal rotation)
'''
fa.fit(data_norm)  # 训练模型

2.8因子载荷矩阵

FA.loadings_

factor_load=fa.loadings_
print("因子载荷矩阵:\n",factor_load)

对应的数学公式及分析：

2.9求各类方差贡献

# 各因子对总方差贡献
factor_contribute=np.sum(factor_load**2,axis=0)
print("各因子对总方差贡献:\n",factor_contribute)# 共性方差
same_variance=np.sum(factor_load**2,axis=1)# 特殊方差
specific_variance=1-same_variance  # 如果数据已经标准化，则 共性方差+特殊方差=1
print("特殊方差:\n",specific_variance)

2.10求因子得分

法一：自己用矩阵运算

# 求因子得分
ss=inv(np.diag(specific_variance))  # 因子得分函数的一部分
factor_score=data_norm@ss@factor_load@inv(factor_load.T@ss@factor_load)
print("因子得分:\n",factor_score)

法二：FA.transform()

# 计算因子得分
factor_score=fa.transform(data_norm)
print("因子得分:\n",factor_score)

注意：两种方法求出的因子得分有一些差别，具体原因不明。

对应的数学公式及分析：

至此，因子分析的步骤已经完成！

下面用因子得分进行评价

2.11计算评价得分

# 计算评价得分
evaluate=factor_score@factor_contribute/sum(factor_contribute)  # 以 各因子对总方差贡献 为权重，求因子得分的加权平均，即可得到评价得分
print("评价得分:\n",evaluate)

对应的数学公式及分析：

2.12按评价得分排序

address_sort=np.argsort(-evaluate)
result_sort=np.zeros(17)
result_sort[address_sort]=np.arange(1,18)
print("排名次序:\n",result_sort)

2.13代码汇总

import numpy as np
from numpy.linalg import inv
import pandas as pd
from scipy.stats import zscore
from  factor_analyzer import FactorAnalyzer as FA
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity  # 用于Bartlett's球状检验
from factor_analyzer.factor_analyzer import calculate_kmo  # 用于KMO检验# 导入数据
data=np.loadtxt("data12_5_1.txt")
print(data)# 数据标准化
data_norm=zscore(data,ddof=1)# Bartlett's球状检验:计算巴特利特P值
chi_square_value,p_value=calculate_bartlett_sphericity(data_norm)
print("Bartlett's球状检验:\n","卡方值:",chi_square_value,"\nP值(若<0.05,则说明变量间有相关性,则可因子分析):",p_value)# KMO检验:计算KMO值
kmo_all,kmo_model=calculate_kmo(data_norm)
print("\nKMO(KMO值介于0~1之间，越接近1，说明变量间相关性越强，偏相关性越弱，因子分析的效果越好。\n""若KMO值>0.6,说明变量间有相关性,则可因子分析):",kmo_model)# 求相关系数矩阵
r=np.corrcoef(data_norm.T)# 求相关系数矩阵的特征值、特征向量
eigval,eigvec=np.linalg.eig(r)
print("特征值:\n",eigval)
print("特征向量:\n",eigvec)# 对特征值、特征向量按从大到小排序
address_sort=np.argsort(-eigval)
result_sort=np.zeros(len(eigval))
result_sort[address_sort]=np.arange(0,len(eigval))
result_sort=[int(i) for i in result_sort]
print("特征值从大到小的顺序",result_sort)eigval=eigval[result_sort]
eigvec=eigvec[:,result_sort]
print("排序后的特征值:\n",eigval)
print("排序后的特征向量:\n",eigvec)contibute_ratio=eigval/sum(eigval)
print("贡献率:\n",contibute_ratio)
print("累积贡献率:\n",np.cumsum(contibute_ratio))# 用此法也能获得相关系数矩阵的特征值
fa0=FA(n_factors=len(r),rotation=None)
fa0.fit(data_norm)
val1,val2=fa0.get_eigenvalues()  # 注意：这里会得到2组特征值，但只有第一组是需要的?
print("fa0的特征值:\n",val1)'''由于前3个公共因子的累计贡献率已经超过0.85，因此认为用3个公共因子 就能较好地进行评价'''
n=3
# 构建模型
fa=FA(n_factors=n,rotation="varimax")
'''FA
参数解读：
n_factors:公共因子数
rotation:因子旋转方式
有以下的选择：
varimax (orthogonal rotation):正交旋转方法，方差最大化
promax (oblique rotation):斜交旋转方法，追求载荷阵中0尽可能多，方法是取正交旋转得到的载荷阵元素的幂次。要求幂次尽可能低，公共因子间的相关性尽可能低。
oblimin (oblique rotation)斜交旋转方法，追求载荷阵中0尽可能多，设置公共因子之间的相关系数在较小值。
oblimax (orthogonal rotation)
quartimin (oblique rotation)
quartimax (orthogonal rotation)正交旋转方法，要求每个原始变量仅与一个公共因子强相关，而与其余公共因子近似不相关。 
equamax (orthogonal rotation)
'''
fa.fit(data_norm)  # 训练模型factor_load=fa.loadings_
print("因子载荷矩阵:\n",factor_load)# 各因子对总方差贡献
factor_contribute=np.sum(factor_load**2,axis=0)
print("各因子对总方差贡献:\n",factor_contribute)# 共性方差
same_variance=np.sum(factor_load**2,axis=1)
# 特殊方差
specific_variance=1-same_variance  # 如果数据已经标准化，则 共性方差+特殊方差=1
print("特殊方差:\n",specific_variance)# 求因子得分
ss=inv(np.diag(specific_variance))  # 因子得分函数的一部分
factor_score=data_norm@ss@factor_load@inv(factor_load.T@ss@factor_load)
print("因子得分:\n",factor_score)
print("因子得分2:\n",fa.transform(data_norm))# 计算评价得分
evaluate=factor_score@factor_contribute
print("评价得分:\n",evaluate)address_sort=np.argsort(-evaluate)
result_sort=np.zeros(17)
result_sort[address_sort]=np.arange(1,18)
print("排名次序:\n",result_sort)