1.数学原理
1.1数学模型
1.2正交因子模型假设
注意:下面的推导都是基于这一假设。因此,这里的模型都是属于正交因子模型。
1.3正交因子模型的协方差结构
1.4各类方差贡献的介绍
在1.3正交因子模型的协方差结构中,我们介绍了“方差贡献”,下面给出关于“方差贡献”更详细的介绍。
1.5因子表示具有不唯一性
利用此性质,我们可以选择不同的正交矩阵,做以下操作:
从而获得容易解释的载荷矩阵和因子。
这一操作被称为“因子旋转”。
因子旋转的类型
因子载荷的正交变换和伴随的因子正交变换为因子正交旋转。
进一步地, 可以修改正交因子模型, 允许共同因子之间相关, 称这样的变换为因子斜交旋转(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)
- 注意:此法得到的特征值没有按从大到小的顺序排列,需要自行排序
- 排序的目的:后面需要利用特征值来求累积贡献率,然后利用累积贡献率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)
参考资料
9 因子分析 | 多元统计分析讲义 (pku.edu.cn)
Python数学建模算法与应用 (司守奎,孙玺菁主编)
标签:sort,Python,步骤,np,因子,因子分析,factor,print,rotation From: https://blog.csdn.net/Yhw20040823/article/details/141341218