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)