当前位置：物联沃-IOTWORD物联网 > 技术教程 > Python之因子分析详细步骤
代码收藏家代码收藏家 技术教程

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检验

    1. Bartlett's球状检验:计算巴特利特P值
             若P值<0.05,则说明变量间有相关性,则可因子分析
    2. 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)

  • 注意：此法得到的特征值没有按从大到小的顺序排列，需要自行排序
  • 排序的目的：后面需要利用特征值来求累积贡献率，然后利用累积贡献率\geqslant85%来挑选前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数学建模算法与应用 (司守奎，孙玺菁主编)

    作者：whY的笔记

    物联沃分享整理
    物联沃-IOTWORD物联网 » Python之因子分析详细步骤

    代码收藏家 普通

    分享到：

    相关推荐

    发表回复