1. 簡介
- 因子分析是一種研究觀測變數變動的共同原因和特殊原因, 從而達到簡化變數結構目的的多元統計方法.
- 因子分析模型是主成分分析的推廣, 也是利用降維的思想, 將複雜的原始變數歸結為少數幾個綜合因子的一種多變數統計分析方法.
1.1 應用
- 尋求變數的基本結構, 簡化變數系統.
- 用於分類, 根據因子得分值, 在因子軸所構成的空間中將變數或者樣本進行分類 (能夠分析樣品間差異的原因).
1.2 型別
- R型因子分析: 研究變數之間的相關關係.
- Q型因子分析: 研究樣本之間的相關關係.
2. 因子分析模型
- 因子分析模型主要涉及矩陣的相關運算.
- 在日常分析中使用最多的就是 R 型因子分析, 下面也將主要介紹 R 型因子的分析模型, 可以對比 Q 型因子分析模型加強對模型的理解.
2.1 因子分析的數學模型
2.1.1 R 型因子分析模型
- 概述:
- R型因子分析是將每一個變數都表示成公共因子的線性函式與特殊因子之和, 即
\[X_i = a_{i1}F_1 + a_{i2}F_2 + \cdots + a_{im}F_m + \epsilon_i, \quad (i=1,2,\cdots, p) \tag{2.1.1-1} \]上式中的\(F_1,F_2,\cdots,F_m\)稱為公共因子, \(\epsilon_i\)稱為\(X_i\)的特殊因子。該模型可用矩陣表示為：\[X = AF + \epsilon \]\[\begin{bmatrix} X_1 \\ X_2 \\ \vdots \\ X_p \end{bmatrix} = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_p \end{bmatrix} + A_{p \times m}\begin{bmatrix} F_1 \\ F_2 \\ \vdots \\ F_m \end{bmatrix} \tag{2.1.1-2} \]- 構造模型滿足:
- \(m \le p\)
- \(Cov(F, \epsilon) = 0\), 即公共因子與特殊因子是不相關的
- \(D_F = D(F) = I_m\), 即各個公共因子不相關且方差為1
- 各個特殊因子不相關, 方差不要求相等
\[D_{\epsilon} = D(\epsilon) = \begin{bmatrix} \sigma_1^2 & & & 0 \\ & \sigma_2^2 & & \\ & & \ddots & \\ 0 & & & \sigma_p^2 \end{bmatrix}\]\[D_{\epsilon} = D(\epsilon) = \begin{bmatrix} \sigma_1^2 & & & 0 \\ & \sigma_2^2 & & \\ & & \ddots & \\ 0 & & & \sigma_p^2 \end{bmatrix}\] - 公共因子(潛在因子)是不可觀測變數且只存在於某種理論意義之中, 可理解為在高維空間中的互相垂直的m個座標軸。雖然潛在變數不能直接測得, 但它一定與某些可測變數有著某種程度的關聯。
2.1.2 Q 型因子分析 (因子得分)
- 概述:
- 類似地, Q 型因子分析數學模型可表示為:
\[X_i = a_{i1}F_1 + a_{i2}F_2 + \cdots + a_{im}F_m + \epsilon_i, \quad (i=1,2,\cdots, n) \tag{2.1.1-3} \]- Q型因子分析與R型因子分析模型的差體現在\(X_1, X_2, \cdots, X_n\)表示的是n個樣品。
2.2 主成分分析與因子分析的異同
- 相同點:
- R 型或 Q 型因子分析都用公因子 F 代替 X, 一般要求\(m<p, m<n\), 因此因子分析與主成分分析一樣, 也是一種降低變數維度數的方法.
- 因子分析求解過程同主成分分析類似, 也是從一個協方差陣出發的.
- 區別:
- 主成分分析的數學模型本質上是一種線性變化, 將原始座標變換到變異程度大的方向, 突出資料變異的方向, 歸納重要資訊。
- 因子分析從本質上看是從顯在變數去 "提煉" 潛在因子的過程。並且因子的形式也不是唯一確定的。一般來說, 作為 "自變數" 的因子\(F_1,F_2,\cdots,F_m\)是不可直接觀測的.
2.3 因子載荷陣
2.3.1 因子載荷陣不唯一的原因
- 變數X的協差陣\(\Sigma\)的分解式為
\[\begin{aligned}
D(X) & = D(AF + \epsilon) = E[(AF + \epsilon)(AF + \epsilon)'] \\
& = AE(FF')A' + AE(F\epsilon') + E(\epsilon F')A' + E(\epsilon\epsilon') \\
& = AD(F)A' + D(\epsilon) \\
& = AA' + D(\epsilon)
\end{aligned}
\]
如果X為標準化的隨機向量, 則\(\Sigma\)就是相關矩陣\(R = (\rho_{ij}), 即R = AA' + D_{\epsilon}\)
- 對於\(m \times m\)的正交矩陣T, 令\(A^* = AT, F^* = T'F\), 模型可表示為: \(X = A^*F^* + \epsilon\)
由於\[Cov(F^*, \epsilon) = E(F^*\epsilon') = T'E(F\epsilon') = 0 \\ D(F^*) = T'D(F)T = T'T = I_{m \times m}\]所以仍滿足模型的條件
同樣\(\Sigma\)也可以分解為\(\Sigma = A^*A^{*'} + D_{\epsilon}\)
2.3.2 因子載荷陣的統計意義
- 因子載荷\(a_{ij}\)的統計意義
- 對於因子模型
\[X_i = a_{i1}F_1 + a_{i2}F_2 + \cdots + a_{ij}F_j + \cdots + a_{im}F_m + \epsilon_i \quad (i=1,2,\cdots,p) \]- 可得\(X_i與F_j\)的協方差為:
\[\begin{aligned} Cov(X_i,F_j) & = Cov(\sum_{k=1}^ma_{ik}F_k + \epsilon, F_j) \\ & = Cov(\sum_{k=1}^ma_{ik}F_k, F_j) + Cov(\epsilon_i, F_j) \\ & = a_{ij} \end{aligned}\]- 若對\(X_i\)做了標準化處理, \(X_i\)的標準差為1, 且\(F_j\)的標準差為1, 有
\[r_{X_i,F_j} = \frac{Cov(X_i,F_j)}{\sqrt{D(X_i)}\sqrt{D(F_j)}} = Cov(X_i,F_j) = a_{ij} \]- 那麼, 對於標準化後的\(X_i\), \(a_{ij}\)是\(X_i\)與\(F_j\)的相關係數, 表示\(X_i\)依賴\(F_j\)的分量(比重)。因此統計學上也稱其為權。
- 歷史原因, 心理學家將它叫做載荷, 表示第i個變數在第j個公共因子上的載荷, 反應了第i個變數在第j個公共因子上的相對重要性.
- 變數共同度的統計意義
- 由因子模型, 知
\[\begin{aligned} D(X_i) & = a_{i1}^2D(F_1) + a_{i2}^2D(F_2) + \cdots + a_{im}^2D(F_m) + D(\epsilon_i) \\ & = a_{i1}^2 + a_{i2}^2 + \cdots + a_{im}^2 + D(\epsilon_i) \\ & = h_i^2 + \sigma_i^2 \end{aligned}\]- 變數\(X_i\)的共同度為: \(h_i^2 = \sum_{j=1}^ma_{ij}^2 \quad i=1,2,\cdots,p\)
- 如果對\(X_i\)作了標準化處理, 有
\[1 = h_i^2 + \sigma_i^2
\]
- 公因子\(F_i\)的方差貢獻率\(g_j^2\)的統計意義
- 設因子載荷矩陣為\(A\), 稱第j列元素的平方和, 即:
\[g_j^2 = \sum_{i=1}^pa_{ij}^2 \quad j=1,2,\cdots,m \]為公共因子\(F_j\)對X的貢獻, 即\(g_j^2\)表示同一個公共因子\(F_j\)對各變數所提供的方差貢獻之總和, 反映該因子對全部變數的總影響- 它是衡量每一個公共因子相對重要性的一個尺度。公共因子\(F_j\)的相對重要性可用\(g_j^2/h\)來衡量, 它稱為公共因子\(F_j\)對X的貢獻率.
- 因子分析的目的就是要由原始隨機向量的協差陣\(\Sigma\)或相關陣\(R\)求出因子分析模型的解, 即求出載荷陣\(A\)和特性方差陣\(D_\epsilon\), 並給公因子賦予有實際背景的解釋.
3. 因子載荷矩陣求解
3.1 主成分法
- 概述:
- 在進行因子分析之前先對資料進行一次主成分分析, 然後把前幾個主成分作為未旋轉的公共因子。
- 該方法所得的特殊因子\(\epsilon_1,\epsilon_2,\cdots,\epsilon_p\)之間並不相互獨立, 不完全符合因子模型的假設前提。
- 當共同度較大時, 特殊因子所起的作用較小, 特殊因子之間的相關性所帶來的影響幾乎可以忽略。
3.1.1 求解方法
- 假定從相關陣出發求解主成分, 設有p個變數, 則可以找出p個主成分\(Y_1,Y_2,\cdots,Y_p\)(按由大到小順序排列), 則主成分與原始變數之間存在如下關係式:\[\begin{cases} Y_1 = \gamma_{11}X_1 + \gamma_{12}X_2 + \cdots + \gamma_{1p}X_p \\ Y_2 = \gamma_{21}X_1 + \gamma_{22}X_2 + \cdots + \gamma_{2p}X_p \\ \cdots\cdots \\ Y_p = \gamma_{p1}X_1 + \gamma_{p2}X_2 + \cdots + \gamma_{pp}X_p \\ \end{cases} \tag{3.1.1-1} \]式中, \(\gamma_ij\)為隨機變數X的相關矩陣的特徵根所對應的特徵向量的分量, 因為特徵向量之間彼此正交, 從X到Y的轉換關係式可逆的, 則可得出由Y到X的轉換關係為:\[ \begin{cases} X_1 = \gamma_{11}Y_1 + \gamma_{21}Y_2 + \cdots + \gamma_{p1}Y_p \\ X_2 = \gamma_{12}Y_1 + \gamma_{22}Y_2 + \cdots + \gamma_{p2}Y_p \\ \cdots\cdots \\ X_p = \gamma_{1p}Y_1 + \gamma_{2p}Y_2 + \cdots + \gamma_{pp}Y_p \\ \end{cases} \tag{3.1.1-2} \]對上面每一等式只保留前m個主成分而把後面的部分用\(\epsilon_i\)代替, 則式(3.1.1-2)轉化為:\[\begin{cases} X_1 = \gamma_{11}Y_1 + \gamma_{21}Y_2 + \cdots + \gamma_{m1}Y_m + \epsilon_1 \\ X_2 = \gamma_{12}Y_1 + \gamma_{22}Y_2 + \cdots + \gamma_{m2}Y_m + \epsilon_2 \\ \cdots\cdots \\ X_p = \gamma_{1p}Y_1 + \gamma_{2p}Y_2 + \cdots + \gamma_{mp}Y_m + \epsilon_p \\ \end{cases} \tag{3.1.1-3} \]該式與因子模型(2.1.1-2)相一致, 並且\(Y_i=(i=1,2,\cdots,m)\)之間相互獨立。將\(Y_i\)轉化成合適的公共因子, 只需將主成分\(Y_i\)變成方差為1的變數, 即除以其標準差(特徵根的平方根), 則可令\(F_i = Y_i/\sqrt{\lambda_i}, \quad a_{ij} = \gamma_{ji}\), 式(3.1.1-3)變為:\[\begin{cases} X_1 = a_{11}F_1 + a_{12}F_2 + \cdots + a_{1m}F_m + \epsilon_1 \\ X_2 = a_{21}F_1 + a_{22}F_2 + \cdots + a_{2m}F_m + \epsilon_2 \\ \cdots\cdots \\ X_p = a_{p1}F_1 + a_{p2}F_2 + \cdots + a_{pm}F_m + \epsilon_p \end{cases} \tag{3.1.1-4} \]這與因子模型完全一致, 這樣就得到了載荷矩陣\(A\)和一組初始公共因子(未旋轉)。
- 一般設\(\lambda_1,\lambda_2,\cdots,\lambda_p \quad (\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_p)\)為樣本相關矩陣R的特徵根, \(\gamma_1,\gamma_2,\cdots,\gamma_p\)為對應的標準化正交化特徵向量。設\(m < p\), 則因子載荷矩陣A的一個解為:\[\widehat{A} = (\sqrt{\lambda_1}\gamma_1, \sqrt{\lambda_2}\gamma_2,\cdots,\sqrt{\lambda_m}\gamma_m) \tag{3.1.1-5} \]共同度的估計為:\[\widehat{h}_i^2 = \widehat{a}_{i1}^2 + \widehat{a}_{i2}^2 + \cdots + \widehat{a}_{im}^2 \tag{3.1.1-6} \]
- 公因子的數目m, 當用主成分進行因子分析時, 可以借鑑確定主成分個數的準則。一般具體問題具體分析, 總之要使所選取的公共因子能夠合理地描述原始變數相關陣的結構, 同時要有利於因子模型的解釋。
3.2 主軸因子法
- 概述:
- 求解思路與主成分發類似, 均從分析矩陣的結構入手, 不同的地方在於, 主成分發是在所有的p個主成分能解釋標準化原始變數所有方差的基礎上進行分析的, 而主軸因子法中, 假定m個公共因子只能解釋原始變數的部分方差, 利用公共因子方差(或共同度)來代替相關矩陣主對角線上的元素1, 並以新得到的這個矩陣(稱為調整相關矩陣)為出發點, 對其分別求解特徵根與特徵向量並得到因子解。
3.2.1 求解方法
- 在因子模型(2.1.1-1)中, 不難得到如下關於X的相關矩陣R的關係式:\[R = AA' + \Sigma_\epsilon \]式中A為因子載荷矩陣, \(\Sigma_\epsilon\)為一對角陣, 其對角元素為相應特殊因子的方差。則稱\(R^* = R -\Sigma_\epsilon = AA'\)為調整相關矩陣, 顯然\(R^*\)的主對角元素不再是1, 而是共同度\(h_i^2\)。分別求解\(R^*\)的特徵值與標準正交特徵向量, 進而求出因子載荷矩陣A。此時, \(R^*\)有m個正的特徵值。設\(\lambda_1^* \ge \lambda_2^* \ge \cdots \ge \lambda_m^*\)為\(R^*\)的特徵根, \(\gamma_1^*, \gamma_2^*, \cdots, \gamma_m^*\)為對應的標準正交話特徵向量。m < p, 則因子載荷矩陣A的一個主軸因子解為:\[\widehat{A} = (\sqrt{\lambda_1^*}\gamma_1^*, \sqrt{\lambda_2^*}\gamma_2^*,\cdots,\sqrt{\lambda_m^*}\gamma_m^*) \tag{3.1.1.7} \]
- 注意到, 上面的分析是以首先得到調整相關矩陣\(R^*\)為基礎的, 而實際上, \(R^*\)與共同度(或相對的剩餘方差)都是未知的, 需要先進行估計。一般先給出一個初始估計, 然後估計出載荷矩陣A後再給出較好的共同度或剩餘方差的估計。初始估計的方法有很多, 可嘗試對原始變數先進行一次主成分分析, 給出初始估計值。
3.3 極大似然法
- 如果假定公共因子F和特殊因子\(\epsilon\)服從正態分佈, 則能夠得到因子載荷和特殊因子方差的極大似然估計。設\(X_1,X_2,\cdots,X_p\)為來自正態總體\(N(\mu,\Sigma)\)的隨機樣本, 其中\(\Sigma = AA' + \Sigma_\epsilon\)。從似然函式的理論知\[L(\mu, \Sigma) = \frac{1}{(2\pi)^{np/2}|\Sigma|^{n/2}}e^{-1/2tr\{\Sigma^{-1}[\sum_{j=1}^n(X_j - \bar{X})(X_j - \bar{X})' + n(\bar{X} - \mu)(\bar{X} - \mu)']\}} \tag{3.3.1-1} \]它通過\(\Sigma\)依賴於A和\(\Sigma_\epsilon\), 但式(3.3.1)不能唯一確定A, 為此, 新增如下條件:\[A'\Sigma_\epsilon^{-1}A = \Lambda \tag{3.3.1-2} \]其中\(\Lambda\)是一個對角陣, 用數值極大化的方法可以得到極大似然估計\(\widehat{A}\)和\(\widehat{\Sigma_\epsilon}\)。極大似然估計\(\widehat{A}, \widehat{\Sigma_\epsilon}和\widehat{\mu} = \bar{X}\), 將使\(\widehat{A}'\widehat{\Sigma}_\epsilon^{-1}\widehat{A}\)為對角陣, 使公式(3.3.1-2)達到最大。
4. 公因子重要性的分析
4.1 因子旋轉
- 概述
- 不管用何種方法確定初始因子載荷矩陣A, 它們都不是唯一的。
- 因所得的初始因子解，各主因子的典型代表變數不是很突出，容易使因子的意義含糊不清，不便於對實際問題進行分析，出於這種考慮，可以對初始公共因進行線性組合，即進行因子旋轉。
- 分類:
- 正交旋轉: 由初始載荷矩陣A右乘一正交陣得到, 經過正交旋轉得到的新公共因子仍然保持彼此獨立的性質
- 斜交旋轉: 放棄了因子之間彼此獨立的限制, 形式更簡潔, 其實際意義更容易理解。
- 說明：
- 對於一個具體問題做因子旋轉，有時需要多次才能得到滿意效果。每一次旋轉後，矩陣各列平方的相對方差之和總會比上一次有所增加。如此繼續下去，當總方差的改變不大時，就可以停止旋轉，這樣就得到新的一組公共因子及相應的因子載荷矩陣，使得其各列元素的相對方差之和最大。
4.2 因子得分概述
- 因子得分就是公共因子\(F_1,F_2,\cdots,F_m\)在每一個樣品點上的得分。
4.2.1 計算公式
- 因子模型中，公共因子的個數少於原始變數的個數，且是不可觀測的隱變數，載荷矩陣A不可逆, 因而不能直接求得公共因子用原始變數表示的精確線性組合。可採用迴歸的思想求出線性組合係數的估計值，即建立如下以公共因子為因變數，原始變數為自變數的迴歸方程：\[F_j = \beta_{j1}X_1 + \beta_{j2}X_2 + \cdots + \beta_{jp}X_p, \quad j=1,2,\cdots,m \tag{4.2.1-1} \]此處因為原始變數與公共因子變數均為標準化變數, 所以迴歸模型中不存在常數項。在最小二乘意義下, 可以得到F的估計值:\[\widehat{F} = A'R^{-1}X \tag{4.2.1-2} \]式中, A為因子載荷矩陣, R為原始變數的相關陣, X為原始變數向量。
5. 因子分析步驟與邏輯框圖
- 瞭解完因子分析的數學模型及因子載荷矩陣的求解方法後, 下面歸納一下因子分析的一般步驟.
5.1 步驟
5.1.1 選取原始變數進行標準化
- 根據研究的問題選取合適的原始變數, 並對變數做無量綱化, 通常有區間縮放, 標準化和歸一化等方法.
5.1.2 適用性驗證
- 並非所有的資料都是適合因子分析的, 因此在因子分析之前需要做適用性檢驗, 即分析變數之間的相關性.
- 共有三種檢驗方法:
- 計算相關係數矩陣:
- 滿足大多數簡單相關係數大於 0.3 即適合做因子分析.
- Bartlett 球度檢驗
- 用於檢驗相關陣中各變數間的相關性，是否為單位陣，即檢驗各個變數是否各自獨立.
- KMO Kaiser-Meyer-Olkin) 檢驗統計量: 用於比較變數間簡單相關係數和偏相關係數的指標
- 取值在0和1之間, 0.9以上表示非常適合; 0.8表示適合; 0.7表示一般; 0.6 表示不太適合; 0.5 以下表示極不適合.
- 計算相關係數矩陣:
5.1.3 求解初始公共因子及因子載荷矩陣
- 計算特徵值和特徵向量
- 計算累計貢獻率大於 80% 的公共因子個數
5.1.4 因子旋轉
- 包括正交旋轉和斜交旋轉, 一般使用較多的是正交旋轉方法以獲得相互獨立的因子.
5.1.5 計算因子得分
5.1.6 根據因子得分值進行進一步分析
5.2 邏輯框圖
6. python示例
6.1 factor_analyzer 包
- API:
- FactorAnalyzer(rotation=None, n_factors=n, method='principal')
- 引數:
rotation: 旋轉的方式 (包括 None: 不旋轉, 'varimax': 最大方差法,'promax': 最優斜交旋轉 )n_factors: 公因子的數量;method: 因子分析的方法 ( 包括 'minres' : 最小殘差因子法,'principal' : 主成分分析法 )
6.2 示例
- 導包
import pandas as pd
import numpy as np
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity, calculate_kmo
from factor_analyzer import FactorAnalyzer
import seaborn as sns
import matplotlib.pyplot as plt
- 讀取樣本資料
plt.style.use({'figure.figsize':(20, 8)})
score_df = pd.read_csv('src_data/factor/test02.csv')
- 對資料進行標準化, 並進行適用性驗證
std_score_df = (score_df - score_df.mean()) / score_df.std()
# 計算皮爾森相關係數
corr_score_df = std_score_df.corr()
- 繪製熱圖, 檢視變數之間的相關性
# 從熱圖明細看出變數間具有相關性
sns.heatmap(corr_score_df, cmap='Blues', annot=True)
- KMO Kaiser-Meyer-Olkin) 檢驗統計量
# KMO Kaiser-Meyer-Olkin) 檢驗統計量
kmo_per_variable, kmo_total = calculate_kmo(std_score_df)
# kmo_total = 0.8148048571126616 適合因子分析
- 計算特徵值和特徵向量
eig_value, eig_vector = nlg.eig(corr_score_df)
#利用變數名和特徵值建立一個 DataFrame
eig_df = pd.DataFrame()
eig_df['names'] = std_score_df.columns
eig_df['eig_value'] = eig_value
eig_df.sort_values(by='eig_value', ascending=False, inplace=True)
eig_df.set_index('names').plot(grid=True)
- 繪圖確定初始公共因子個數
cum_eig_df = pd.DataFrame()
cum_eig_df['factors'] = pd.Series(range(1, 11))
cum_eig_df['cum_eig_value_ratio'] = np.cumsum(eig_df['eig_value'].values) / np.sum(eig_df['eig_value'].values)
# 從累計折線圖中可看出當因子個數 5 時, 累計貢獻率大於 0.8
cum_eig_df.set_index('factors').plot(grid=True)
- 因子分析
# 因子旋轉使用正交旋轉最大方差法, 因子載荷矩陣使用主成分分析法
rotation_fa = FactorAnalyzer(5, method='principal', rotation='varimax')
# 求解載荷矩陣
rota_factor_df = pd.DataFrame(np.abs(rotation_fa.loadings_), index=std_score_df.columns)
# 繪製載荷矩陣熱圖
sns.heatmap(rota_factor_df, annot=True, cmap='Blues')
# 至此該樣本中提取出 5 個綜合因子
7. 參考
- 多元統計分析-何曉群編著.中國人民大學出版社
- factor_analyzer 包