1、起源
因子分析最早由英國心理學家C.Spearman發表了第一篇有關因子分析的文章《對智力測驗得分進行統計分析》，從中提出的：他發現學生的英語、法語和古典語成績非常有相關性，他認為這三門課程背后有一個共同的因素驅動，最后將這個因素定義為“語言能力”。由此解開了因子分析的序幕。
2、基本思想
因子分析通過研究眾多變量間的內部依賴關系，探求觀測數據中的基本結構，並用少數幾個假想變量（因子）來表示基本的數據結構。原始變量是可觀測的顯在變量，而假想變量是不可觀測的潛在變量，稱為因子。
如在企業品牌形象研究中，消費者可通過24個指標構成的評價體系，對商場的方方面面進行優劣評價，但是消費者會主要關心3個放米娜，即商店的環境、商店的服務和商品的價格。因子分析可以通過24個變量，找出反應商店環境、商店服務水平和商品價格的3個潛在因子，對商店進行綜合評價。
$ X{i} = u{_i} + a{_i}{_1}F{_1} + a{_i}{_2}F{_2} + a{_i}{_3}F{_3} + e{_i} $
$ F{_1} 、F{_2}、F{_3} $就是因子，不被包含的部分 \(e{_i}\) 叫特殊因子。
3、因子分析特點
-
因子變量的數量 << 原有的指標變量數量。
-
因子變量不是對原有變量的取舍，而是根據原有變量的信息進行重新組合，能夠反映原有變量大部分的信息。
-
因子變量不存在線性相關關系。
-
因子變量具有命名解釋性，即該變量是對某些原始變量信息的綜合與反應。
4、算法用途
-
降維，減少分析變量的數量
-
分類，將內部具有關聯變量/樣本划分一類
5、分析步驟
-
a. 選擇分析變量，並進行標准化處理
-
b. 計算變量間的相關系數矩陣、相關系數矩陣的特征值、特征向量
-
c. 充分性檢驗：KMO和巴特萊特球度檢驗，驗證變量是否適合做因子分析。
KMO(Kaiser-Meyer-Olkin)檢驗
KMO檢驗統計量是用於比較變量間簡單相關系數矩陣和偏相關系數的指標。數學定義為：
是否適合因子分析：KMO在0.9以上非常適合；0.8表示適合，0.7表示一般，0.6表示不太適合；0.5以下表示極不適合。
巴特萊特球度檢驗
該檢驗以原有變量的相關系數矩陣為出發點，其零假設\(H_0\)是：相關系數矩陣為單位矩陣，即相關系數矩陣主對角元素均為1，非主對角元素均為0（即原始變量之間無相關關系）。
依據檢驗統計量服從卡方分布，如果統計量卡方值較大且對應的sig值小於給定的顯著性水平\(\alpha\)時，零假設不成立。即說明相關系數矩陣不太可能是單位矩陣，變量之間存在相關關系，適合做因子分析。
-
d. 提取公因子：只取方差>1或特征值>1（方差小於1的因子其貢獻可能很小。）；因子的累積方差貢獻率達到80%。
-
e. 因子旋轉：提取因子的實際意義更容易解釋
正交旋轉和斜交旋轉，其中主要使用正交旋轉的方差最大旋轉法。
- f. 計算因子得分
6、應用實例
因子分析用到的庫 factor_analyzer
pip install factor_analyzer
6.1 數據處理
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
# 導入數據
df = pd.read_csv(r'C:\Users\Desktop\bfi.csv')
df.head()
# 刪除無關列
df.drop(["gender", "education", "age", "Unnamed: 0"], axis=1, inplace=True)
# 查看是否有缺失值
df.isnull().sum()
可以看到數據存在一些缺失值，需要對缺失值做刪除處理
# 刪除缺失值
df.dropna(inplace=True)
df.shape
處理完缺失值，樣本數據大小為2436個樣本 * 25個變量。
6.2 充分性檢驗
# 導入所需庫
from factor_analyzer import FactorAnalyzer
from factor_analyzer.factor_analyzer import calculate_kmo
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
# 因子分析可靠性檢驗
kmo_all, kmo_model = calculate_kmo(df)
chi_square_value, p_value = calculate_bartlett_sphericity(df)
print("kmo_all:", kmo_all, end="\n\n")
print("kmo_model:", kmo_model, end="\n\n")
print("chi_square_value:", chi_square_value, end="\n\n")
print("p_value:", p_value)
KMO檢驗表明KMO統計量值為0.85，表明數據適合進行因子分析；同時巴特萊特球形檢驗p值為0，拒絕原假設（相關系數矩陣為單位矩陣），變量之間存在相關關系，適合做因子分析。
6.3 提取公因子
進行探索性因子分析，確定提取公因子個數。計算相關系數矩陣特征根和特征向量
# 探索性因子分析
fa = FactorAnalyzer(25, rotation=None)
fa.fit(df)
# 相關系數矩陣的特征根和特征向量
ev, v = fa.get_eigenvalues()
ev, v
# 根據特征根>1發現，可提取6個公因子
繪制碎石圖進一步確定因子個數
# 繪制碎石圖，選擇因子數
plt.scatter(range(1,df.shape[1]+1),ev)
plt.plot(range(1,df.shape[1]+1),ev)
plt.title('Scree Plot')
plt.xlabel('Factors')
plt.ylabel('Eigenvalue')
plt.grid()
plt.show()
6.4 因子旋轉
建立明確的因子分析模型，並對因子進行方差最大化的正交旋轉。
# 建立因子分析模型
fa_six = FactorAnalyzer(6, rotation="varimax")
fa_six.fit(df)
# 輸出因子的載荷
fa_six.loadings_
# pd.DataFrame(fa_six.loadings_, index=df.columns)
結果看起來並不直觀，無法看出哪些因子對變量的解釋程度較高，對數據進行可視化展示結果。
import seaborn as sns
df_cm = pd.DataFrame(np.abs(fa_six.loadings_), index=df.columns)
plt.figure(figsize = (14,14))
ax = sns.heatmap(df_cm, annot=True, cmap="BuPu")
# 設置y軸的字體的大小
ax.yaxis.set_tick_params(labelsize=15)
plt.title('Factor Analysis', fontsize='xx-large')
# 設置y軸標簽
plt.ylabel('Personality items', fontsize='xx-large')
# 保存圖片
# plt.savefig(r'C:\Users\Desktop\factorAnalysis.png', dpi=500)
由上圖，因子6對所有變量都沒有高負荷，也不容易解釋，需要調整因子個數，選擇5個公因子。重復上面的步驟：
# 建立因子分析模型，設置公因子個數為5
fa_five = FactorAnalyzer(5, rotation="varimax")
fa_five.fit(df)
import seaborn as sns
df_cm = pd.DataFrame(np.abs(fa_five.loadings_), index=df.columns)
plt.figure(figsize = (14,14))
ax = sns.heatmap(df_cm, annot=True, cmap="BuPu")
# 設置y軸的字體的大小
ax.yaxis.set_tick_params(labelsize=15)
plt.title('Factor Analysis', fontsize='xx-large')
# 設置y軸標簽
plt.ylabel('Personality items', fontsize='xx-large')
根據上圖結果：
因子1在變量（N1, N2, N3, N4, N5）上具有較高載荷，可定義因子1為神經質因子。
因子2在變量（E1, E2, E3, E4, E5）上具有較高載荷，可定義因子2為外向型因子。
因子3在變量（C1, C2, C3, C4, C5）上具有較高載荷，可定義因子3為盡責性因子。
因子4在變量（A1, A2, A3, A4, A5）上具有較高載荷，可定義因子4為認同性因子。
因子5在變量（O1, O2, O3, O4, O5）上具有較高載荷，可定義因子5為開放性因子。
# 方差累計貢獻
fa_v = fa_five.get_factor_variance()
fa_dt = pd.DataFrame({
"特征根": fa_v[0],
"方差貢獻率": fa_v[1],
"方差累計貢獻率": fa_v[2]
})
fa_dt
5個因子的特征值之和占特征值總和的42.36%，也可以說5個因子解釋了全部變量的42.36%的信息。
6.5 計算因子得分
# 計算因子得分
score = fa_five.transform(df)
score
# 計算綜合得分
x = score @ (fa_v[1])
result = pd.DataFrame(x, columns=["綜合得分"], index=df.index)
result.sort_values(by="綜合得分", ascending=False, inplace=True)
result
完整代碼：
import numpy as np
import pandas as pd
from factor_analyzer import FactorAnalyzer
from factor_analyzer.factor_analyzer import calculate_kmo
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
import matplotlib.pyplot as plt
%matplotlib inline
# 導入數據
df = pd.read_csv(r'C:\Users\wy\Desktop\bfi.csv')
df.head()
# 刪除無關列
df.drop(["gender", "education", "age", "Unnamed: 0"], axis=1, inplace=True)
# 查看是否有缺失值
df.isnull().sum()
# 刪除缺失值
df.dropna(inplace=True)
# 因子分析可靠性檢驗
kmo_all, kmo_model = calculate_kmo(df)
chi_square_value, p_value = calculate_bartlett_sphericity(df)
print("kmo_all:", kmo_all, end="\n\n")
print("kmo_model:", kmo_model, end="\n\n")
print("chi_square_value:", chi_square_value, end="\n\n")
print("p_value:", p_value)
# 探索性因子分析
fa = FactorAnalyzer(25, rotation=None)
fa.fit(df)
# 相關系數矩陣的特征根和特征向量
ev, v = fa.get_eigenvalues()
ev, v
# 根據特征根>1發現，可提取6個公因子
# 繪制碎石圖，選擇因子數
plt.scatter(range(1,df.shape[1]+1),ev)
plt.plot(range(1,df.shape[1]+1),ev)
plt.title('Scree Plot')
plt.xlabel('Factors')
plt.ylabel('Eigenvalue')
plt.grid()
plt.show()
# 建立因子分析模型
fa_six = FactorAnalyzer(6, rotation="varimax")
fa_six.fit(df)
# 輸出因子的載荷
fa_six.loadings_
# 建立因子分析模型
fa_five = FactorAnalyzer(5, rotation="varimax") # 根據6個公因子的熱力圖發現Factor6在每個變量上均無載荷。故調整為5個公因子。
fa_five.fit(df)
import seaborn as sns
df_cm = pd.DataFrame(np.abs(fa_five.loadings_), index=df.columns)
plt.figure(figsize = (14,14))
ax = sns.heatmap(df_cm, annot=True, cmap="BuPu")
# 設置y軸的字體的大小
ax.yaxis.set_tick_params(labelsize=15)
plt.title('Factor Analysis', fontsize='xx-large')
# 設置y軸標簽
plt.ylabel('Personality items', fontsize='xx-large')
# 方差累計貢獻
fa_v = fa_five.get_factor_variance()
fa_dt = pd.DataFrame({
"特征根": fa_v[0],
"方差貢獻率": fa_v[1],
"方差累計貢獻率": fa_v[2]
})
fa_dt
# 計算因子得分
score = fa_five.transform(df)
score
# 計算綜合得分
x = score @ (fa_v[1])
result = pd.DataFrame(x, columns=["綜合得分"], index=df.index)
result.sort_values(by="綜合得分", ascending=False, inplace=True)
result
參考：