正交因子分析
- 目的
- 數學原理
- 參數估計方法
- 主成分法
- 主因子法
- 極大似然法
- 因子旋轉
- 模型檢驗
- 因子得分
- 加權最小二乘法
- 回歸法
- 代碼實現
- 注意事項
- 例子
- Reference
目的
FactorAnalysis的目的是從多個高度相關的觀測變量中提取出少數幾個LatentFactor,這些因子代表了變量背后的共通結構,從而實現降維并提升可解釋性。
假設對一組學生進行了以下六門課程的測試:語文、英語、數學、物理、化學、生物,發現語文和英語成績之間高度相關,數學、物理、化學、生物也彼此高度相關。此時可以猜測:這些成績可能是由兩個更基本的“能力”決定的,比如語言能力和理科能力。通過因子分析就可以提取出這兩個潛在因子,并發現語文和英語主要由“語言能力”因子決定,理科四門主要由“理科能力”因子解釋。這樣就可以用兩個因子有效地概括了六個變量的結構,同時讓模型更易解釋、更簡潔。
正交因子分析認為因子間是正交的,即不相關的,斜交因子分析中因子間可以是相關的,本文只探討正交因子分析。
數學原理
設 ‘ X ‘ `\mathbf{X}` ‘X‘是一個可觀測的 ‘ m ‘ `m` ‘m‘維隨機向量, ‘ E ? ( X ) = μ , Cov ? ( X ) = Σ = ( σ i j ) ‘ `\operatorname{E}(\mathbf{X})=\boldsymbol{\mu},\;\operatorname{Cov}(\mathbf{X})=\Sigma=(\sigma_{ij})` ‘E(X)=μ,Cov(X)=Σ=(σij?)‘。正交因子分析的數學模型為:
X = μ + A F + ε { E ? ( F ) = 0 , Cov ? ( F ) = I n E ? ( ε ) = 0 , Cov ? ( ε ) = D = diag ? { σ 1 2 , … , σ m 2 } Cov ? ( F , ε ) = 0 \begin{gathered} \mathbf{X}=\boldsymbol{\mu}+AF+\varepsilon \\ \begin{cases} \operatorname{E}(F)=\mathbf{0},\;\operatorname{Cov}(F)=I_n \\ \operatorname{E}(\varepsilon)=\mathbf{0},\;\operatorname{Cov}(\varepsilon)=D=\operatorname{diag}\{\sigma_1^2,\dots,\sigma_m^2\} \\ \operatorname{Cov}(F,\varepsilon)=\mathbf{0} \end{cases} \end{gathered} X=μ+AF+ε? ? ??E(F)=0,Cov(F)=In?E(ε)=0,Cov(ε)=D=diag{σ12?,…,σm2?}Cov(F,ε)=0??
其中 ‘ F = ( f 1 , … , f n ) T ‘ `F=(f_1,\dots,f_n)^T` ‘F=(f1?,…,fn?)T‘是不可觀測的 ‘ n ‘ `n` ‘n‘維隨機變量, ‘ ε ‘ `\varepsilon` ‘ε‘是不可觀測的 ‘ m ‘ `m` ‘m‘維隨機變量,分別稱 ‘ F ‘ `F` ‘F‘和 ‘ ε ‘ `\varepsilon` ‘ε‘為CommonFactor和SpecificFactor。 ‘ A = ( a i j ) ‘ `A=(a_{ij})` ‘A=(aij?)‘是一個非隨機矩陣, ‘ a i j ‘ `a_{ij}` ‘aij?‘表示公共因子 ‘ f j ‘ `f_j` ‘fj?‘、隨機變量 ‘ X i ‘ `\mathbf{X}_i` ‘Xi?‘的因子載荷。 ‘ a 1 j , a 2 j , … , a i j ‘ `a_{1j},a_{2j},\dots,a_{ij}` ‘a1j?,a2j?,…,aij?‘中至少有兩個不為 ‘ 0 ‘ `0` ‘0‘,否則可將 ‘ f i ‘ `f_i` ‘fi?‘并入到 ‘ ε i ‘ `\varepsilon_i` ‘εi?‘中去; ‘ ε i ‘ `\varepsilon_i` ‘εi?‘也僅出現在 ‘ X i ‘ `\mathbf{X}_i` ‘Xi?‘的表達式中。
上述因子分析模型具有如下性質:
-
‘ Σ = A A T + D ‘ `\Sigma=AA^T+D` ‘Σ=AAT+D‘;
-
若 ‘ X ? = C X ‘ `\mathbf{X}^*=C\mathbf{X}` ‘X?=CX‘,則有:
Y = C μ + C A F + C ε = μ ? + A ? F + ε ? \mathbf{Y}=C\boldsymbol{\mu}+CAF+C\varepsilon=\boldsymbol{\mu}^*+A^*F+\varepsilon^* Y=Cμ+CAF+Cε=μ?+A?F+ε? -
因子載荷不唯一;
-
‘ Cov ? ( X , F ) = A ‘ `\operatorname{Cov}(\mathbf{X},F)=A` ‘Cov(X,F)=A‘,即 ‘ Cov ? ( X i , F j ) = a i j ‘ `\operatorname{Cov}(\mathbf{X}_i,F_j)=a_{ij}` ‘Cov(Xi?,Fj?)=aij?‘
-
令 ‘ h i 2 = ∑ j = 1 n a i j 2 ‘ `h_i^2=\sum\limits_{j=1}^{n}a_{ij}^2` ‘hi2?=j=1∑n?aij2?‘,則有:
Var ? ( X i ) = σ i i = ∑ j = 1 n a i j 2 + σ i 2 = h i 2 + σ i 2 , i = 1 , 2 , … , m \operatorname{Var}(\mathbf{X}_i)=\sigma_{ii}=\sum_{j=1}^{n}a_{ij}^2+\sigma_i^2=h_i^2+\sigma_i^2,\;i=1,2,\dots,m Var(Xi?)=σii?=j=1∑n?aij2?+σi2?=hi2?+σi2?,i=1,2,…,m -
令 ‘ g j 2 = ∑ i = 1 m a i j 2 ‘ `g_j^2=\sum\limits_{i=1}^{m}a_{ij}^2` ‘gj2?=i=1∑m?aij2?‘,則有:
∑ i = 1 m Var ? ( X i ) = ∑ j = 1 n g j 2 + ∑ i = 1 n σ i 2 \sum_{i=1}^{m}\operatorname{Var}(\mathbf{X}_i)=\sum_{j=1}^{n}g_j^2+\sum_{i=1}^{n}\sigma_i^2 i=1∑m?Var(Xi?)=j=1∑n?gj2?+i=1∑n?σi2?
Proof. (1)由[prop:CovMat](3)(4)(5)可得:
Σ = Cov ? ( X ) = Cov ? ( μ + A F + ε , μ + A F + ε ) = Cov ? ( μ , μ + A F + ε ) + Cov ? ( A F , μ + A F + ε ) + Cov ? ( ε , μ + A F + ε ) = Cov ? ( A F , μ ) + Cov ? ( A F ) + Cov ? ( A F , ε ) + Cov ? ( ε , μ ) + Cov ? ( ε , A F ) + Cov ? ( ε ) = A Cov ? ( F ) A T + A Cov ? ( F , ε ) + Cov ? ( ε , F ) A T + D = A A T + D \begin{aligned} \Sigma&=\operatorname{Cov}(\mathbf{X})=\operatorname{Cov}(\boldsymbol{\mu}+AF+\varepsilon,\boldsymbol{\mu}+AF+\varepsilon) \\ &=\operatorname{Cov}(\boldsymbol{\mu},\boldsymbol{\mu}+AF+\varepsilon)+\operatorname{Cov}(AF,\boldsymbol{\mu}+AF+\varepsilon)+\operatorname{Cov}(\varepsilon,\boldsymbol{\mu}+AF+\varepsilon) \\ &=\operatorname{Cov}(AF,\boldsymbol{\mu})+\operatorname{Cov}(AF)+\operatorname{Cov}(AF,\varepsilon)+\operatorname{Cov}(\mathbf{\varepsilon},\boldsymbol{\mu})+\operatorname{Cov}(\varepsilon,AF)+\operatorname{Cov}(\varepsilon) \\ &=A\operatorname{Cov}(F)A^T+A\operatorname{Cov}(F,\varepsilon)+\operatorname{Cov}(\varepsilon,F)A^T+D \\ &=AA^T+D \end{aligned} Σ?=Cov(X)=Cov(μ+AF+ε,μ+AF+ε)=Cov(μ,μ+AF+ε)+Cov(AF,μ+AF+ε)+Cov(ε,μ+AF+ε)=Cov(AF,μ)+Cov(AF)+Cov(AF,ε)+Cov(ε,μ)+Cov(ε,AF)+Cov(ε)=ACov(F)AT+ACov(F,ε)+Cov(ε,F)AT+D=AAT+D?
(2)顯然。
(3)取正交矩陣 ‘ Q ‘ `Q` ‘Q‘,令 ‘ A ? = A Q ‘ `A^*=AQ` ‘A?=AQ‘, ‘ F ? = Q T F ‘ `F^*=Q^TF` ‘F?=QTF‘,則依然有:
E ? ( F ? ) = Q T E ? ( F ) = 0 , Cov ? ( F ? ) = Q T Cov ? ( F ) Q = I n , X = μ + A ? F ? + ε \operatorname{E}(F^*)=Q^T\operatorname{E}(F)=\mathbf{0},\;\operatorname{Cov}(F^*)=Q^T\operatorname{Cov}(F)Q=I_n,\;\mathbf{X}=\boldsymbol{\mu}+A^*F^*+\varepsilon E(F?)=QTE(F)=0,Cov(F?)=QTCov(F)Q=In?,X=μ+A?F?+ε
(4)由[prop:CovMat](3)(4)(5)可得:
Cov ? ( X , F ) = Cov ? ( μ + A F + ε , F ) = Cov ? ( μ , F ) + Cov ? ( A F , F ) + Cov ? ( ε , F ) = A \operatorname{Cov}(\mathbf{X},F)=\operatorname{Cov}(\boldsymbol{\mu}+AF+\varepsilon,F)=\operatorname{Cov}(\boldsymbol{\mu},F)+\operatorname{Cov}(AF,F)+\operatorname{Cov}(\varepsilon,F)=A Cov(X,F)=Cov(μ+AF+ε,F)=Cov(μ,F)+Cov(AF,F)+Cov(ε,F)=A
(5)由(1)即可得到結論。
(6)由(1)可得:
∑ i = 1 m Var ? ( X i ) = tr ? [ Cov ? ( X ) ] = tr ? ( A A T + D ) = ∑ i = 1 m ∑ j = 1 n a i j 2 + ∑ i = 1 n σ i 2 = ∑ j = 1 n ∑ i = 1 m a i j 2 + ∑ i = 1 n σ i 2 = ∑ j = 1 n g j 2 + ∑ i = 1 n σ i 2 \begin{aligned} \sum_{i=1}^{m}\operatorname{Var}(\mathbf{X}_i)&=\operatorname{tr}[\operatorname{Cov}(\mathbf{X})]=\operatorname{tr}(AA^T+D)=\sum_{i=1}^{m}\sum_{j=1}^{n}a_{ij}^2+\sum_{i=1}^{n}\sigma_i^2 \\ &=\sum_{j=1}^{n}\sum_{i=1}^{m}a_{ij}^2+\sum_{i=1}^{n}\sigma_i^2=\sum_{j=1}^{n}g_j^2+\sum_{i=1}^{n}\sigma_i^2 \end{aligned} i=1∑m?Var(Xi?)?=tr[Cov(X)]=tr(AAT+D)=i=1∑m?j=1∑n?aij2?+i=1∑n?σi2?=j=1∑n?i=1∑m?aij2?+i=1∑n?σi2?=j=1∑n?gj2?+i=1∑n?σi2??
?
稱 ‘ h i 2 ‘ `h_i^2` ‘hi2?‘為變量 ‘ X i ‘ `\mathbf{X}_i` ‘Xi?‘的CommonVariance,它反映了公共因子對 ‘ X i ‘ `\mathbf{X}_i` ‘Xi?‘的方差貢獻度。稱 ‘ σ i 2 ‘ `\sigma_i^2` ‘σi2?‘為 ‘ X i ‘ `\mathbf{X}_i` ‘Xi?‘的SpecificVariance,它反映了特殊因子 ‘ ε i ‘ `\varepsilon_i` ‘εi?‘對 ‘ X i ‘ `\mathbf{X}_i` ‘Xi?‘的方差貢獻度。 ‘ g j 2 ‘ `g_j^2` ‘gj2?‘可視為公共因子 ‘ f j ‘ `f_j` ‘fj?‘對 ‘ X 1 , … , X m ‘ `\mathbf{X}_1,\dots,\mathbf{X}_m` ‘X1?,…,Xm?‘的總方差貢獻度。
參數估計方法
主成分法
設觀測變量 ‘ X ‘ `\mathbf{X}` ‘X‘的協方差矩陣 ‘ Σ ‘ `\Sigma` ‘Σ‘,它的特征值從大到小依次為 ‘ λ 1 , … , λ m ‘ `\lambda_1,\dots,\lambda_m` ‘λ1?,…,λm?‘,對應的單位正交特征向量分別為 ‘ l 1 , … , l m ‘ `l_1,\dots,l_m` ‘l1?,…,lm?‘。于是 ‘ Σ ‘ `\Sigma` ‘Σ‘有分解式:
Σ = ( l 1 l 2 ? l m ) ( λ 1 0 ? 0 0 λ 2 ? 0 ? ? ? ? 0 0 ? λ m ) ( l 1 T l 2 T ? l m T ) = ∑ i = 1 m λ i l i l i T \Sigma= \begin{pmatrix} l_1 & l_2 & \cdots &l_m \end{pmatrix} \begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_m \end{pmatrix} \begin{pmatrix} l_1^T \\ l_2^T \\ \vdots \\ l_m^T \end{pmatrix} =\sum_{i=1}^{m}\lambda_il_il_i^T Σ=(l1??l2????lm??) ?λ1?0?0?0λ2??0??????00?λm?? ? ?l1T?l2T??lmT?? ?=i=1∑m?λi?li?liT?
由[prop:CovMat](2)和[theo:PositiveSemidefinite](3)的第五條可知 ‘ λ m ? 0 ‘ `\lambda_m\geqslant0` ‘λm??0‘。當最后 ‘ m ? n ‘ `m-n` ‘m?n‘個特征值較小時, ‘ Σ ‘ `\Sigma` ‘Σ‘有如下近似:
Σ = ∑ i = 1 m λ i l i l i T ≈ ∑ i = 1 n λ i l i l i T + D ^ = A ^ A ^ T + D ^ \Sigma=\sum_{i=1}^{m}\lambda_il_il_i^T\approx\sum_{i=1}^{n}\lambda_il_il_i^T+\hat{D}=\hat{A}\hat{A}^T+\hat{D} Σ=i=1∑m?λi?li?liT?≈i=1∑n?λi?li?liT?+D^=A^A^T+D^
其中:
A ^ = ( λ 1 l 1 ? λ n l n ) , D ^ = diag ? ( Σ ? A ^ A ^ T ) \hat{A}= \begin{pmatrix} \sqrt{\lambda_1}l_1 & \cdots & \sqrt{\lambda_n}l_n \end{pmatrix},\; \hat{D}=\operatorname{diag}(\Sigma-\hat{A}\hat{A}^T) A^=(λ1??l1????λn??ln??),D^=diag(Σ?A^A^T)
與PCA一樣,一般通過使 ‘ ( ∑ i = 1 n λ i ) / ( ∑ i = 1 m λ i ) ‘ `\left(\sum\limits_{i=1}^{n}\lambda_i\right)/\left(\sum\limits_{i=1}^{m}\lambda_i\right)` ‘(i=1∑n?λi?)/(i=1∑m?λi?)‘大于一定比例來選擇 ‘ n ‘ `n` ‘n‘的具體值。
主因子法
令 ‘ A A T = Σ ? D ‘ `AA^T=\Sigma-D` ‘AAT=Σ?D‘。取 ‘ σ ^ 1 2 , … , σ ^ m 2 ‘ `\hat{\sigma}_1^2,\dots,\hat{\sigma}_m^2` ‘σ^12?,…,σ^m2?‘為特殊方差的合理初始估計((1)全零,(2)取 ‘ max ? j ≠ i σ i j ‘ `\max\limits_{j\ne i}\sigma_{ij}` ‘j=imax?σij?‘),則有:
A A T ^ = ( σ 11 ? σ ^ 1 2 σ 12 ? σ 1 m σ 21 σ 22 ? σ ^ 2 2 ? σ 2 m ? ? ? ? σ m 1 σ m 2 ? σ m m ? σ ^ m 2 ) \widehat{AA^T}= \begin{pmatrix} \sigma_{11}-\hat{\sigma}_1^2 & \sigma_{12} & \cdots & \sigma_{1m} \\ \sigma_{21} & \sigma_{22}-\hat{\sigma}_2^2 & \cdots & \sigma_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{m1} & \sigma_{m2} & \cdots & \sigma_{mm}-\hat{\sigma}_m^2 \end{pmatrix} AAT = ?σ11??σ^12?σ21??σm1??σ12?σ22??σ^22??σm2???????σ1m?σ2m??σmm??σ^m2?? ?
取 ‘ A A T ^ ‘ `\widehat{AA^T}` ‘AAT ‘前 ‘ n ‘ `n` ‘n‘個大于 ‘ 0 ‘ `0` ‘0‘的特征值,從大到小依次為 ‘ λ ^ 1 , … , λ ^ n ‘ `\hat{\lambda}_1,\dots,\hat{\lambda}_n` ‘λ^1?,…,λ^n?‘,對應的單位正交特征向量為 ‘ l ^ 1 , … , l ^ n ‘ `\hat{l}_1,\dots,\hat{l}_n` ‘l^1?,…,l^n?‘,則有近似的:
A ^ = ( λ ^ 1 l ^ 1 ? λ ^ n l ^ n ) \hat{A}= \begin{pmatrix} \sqrt{\hat{\lambda}_1}\hat{l}_1 & \cdots & \sqrt{\hat{\lambda}_n}\hat{l}_n \end{pmatrix} A^=(λ^1??l^1????λ^n??l^n??)
令 ‘ σ ^ i 2 = σ i i ? h ^ i 2 ‘ `\hat{\sigma}_i^2=\sigma_{ii}-\hat{h}_i^2` ‘σ^i2?=σii??h^i2?‘,繼續上面的迭代過程以得到穩定的近似解。
Input: 協方差矩陣 ‘ Σ ‘ `\Sigma` ‘Σ‘,初始特殊方差估計
‘ σ ^ 1 2 , … , σ ^ m 2 ‘ `\hat{\sigma}^2_1, \ldots, \hat{\sigma}^2_m` ‘σ^12?,…,σ^m2?‘,目標因子數 ‘ n ‘ `n` ‘n‘
Output: 因子載荷矩陣估計 ‘ A ^ ‘ `\hat{A}` ‘A^‘,特殊方差估計
‘ σ ^ i 2 ‘ `\hat{\sigma}_i^2` ‘σ^i2?‘
初始化 ‘ σ ^ i 2 ‘ `\hat{\sigma}_i^2` ‘σ^i2?‘ 為合理值 構造矩陣
‘ A A T ^ = Σ ? diag ? ( σ ^ 1 2 , … , σ ^ m 2 ) ‘ `\widehat{AA^T} = \Sigma - \operatorname{diag}(\hat{\sigma}_1^2, \ldots, \hat{\sigma}_m^2)` ‘AAT =Σ?diag(σ^12?,…,σ^m2?)‘
對 ‘ A A T ^ ‘ `\widehat{AA^T}` ‘AAT ‘ 做特征值分解,得到部分特征值
‘ λ ^ 1 ? ? ? λ ^ n ‘ `\hat{\lambda}_1 \geqslant \cdots \geqslant \hat{\lambda}_n` ‘λ^1????λ^n?‘,及對應單位正交特征向量
‘ l ^ 1 , … , l ^ n ‘ `\hat{l}_1, \ldots, \hat{l}_n` ‘l^1?,…,l^n?‘ 構造因子載荷矩陣估計:
‘ A ^ = ( a ^ i j ) = ( λ ^ 1 l ^ 1 ? λ ^ n l ^ n ) ‘ `\hat{A}=(\hat{a}_{ij}) = \begin{pmatrix} \sqrt{\hat{\lambda}_1} \hat{l}_1 & \cdots & \sqrt{\hat{\lambda}_n} \hat{l}_n \end{pmatrix}` ‘A^=(a^ij?)=(λ^1??l^1????λ^n??l^n??)‘ 令
‘ h ^ i 2 = ∑ j = 1 n a ^ i j 2 ‘ `\hat{h}_i^2 = \sum\limits_{j=1}^n \hat{a}_{ij}^2` ‘h^i2?=j=1∑n?a^ij2?‘,更新
‘ σ ^ i 2 = σ i i ? h ^ i 2 , i = 1 , 2 , … , m ‘ `\hat{\sigma}_i^2 = \sigma_{ii} - \hat{h}_i^2,\;i=1,2,\dots,m` ‘σ^i2?=σii??h^i2?,i=1,2,…,m‘
極大似然法
推導太過困難,可以參考
J?reskog, K. G. (1963). Statistical Estimation in Factor Analysis. Almqvist and Wicksell.
Lawley, D. N. and Maxwell, A. E. (1971). Factor Analysis as a Statistical Method. Second edition. Butterworths.
因子旋轉
為了提高因子的可解釋性,我們希望每個因子對觀測變量的影響是集中且明顯的,即一個因子主要對少數幾個變量有顯著影響,對其余變量幾乎沒有作用。這種結構反映在因子載荷矩陣 ‘ A ‘ `A` ‘A‘上即為 ‘ A ‘ `A` ‘A‘每一列的元素 ‘ a i j , i = 1 , 2 , … , m ‘ `a_{ij},\;i=1,2,\dots,m` ‘aij?,i=1,2,…,m‘不是均勻地分布在中間水平,而是趨于兩極分化:其絕對值要么接近于 ‘ 0 ‘ `0` ‘0‘,要么較大。這樣可以使得每個因子更容易被識別和解釋——因為它只與一小組變量高度相關。這種結構等價于希望載荷矩陣 ‘ A ‘ `A` ‘A‘的每一列具有稀疏性,從而便于賦予因子明確的語義標簽。
由[prop:FactorAnalysis](3)可知在初步求得因子載荷矩陣 ‘ A ‘ `A` ‘A‘后,可以使用一個正交矩陣右乘 ‘ A ‘ `A` ‘A‘,此時仍能得到一個因子模型。使用正交矩陣來右乘 ‘ A ‘ `A` ‘A‘相當于是對因子 ‘ F ‘ `F` ‘F‘進行旋轉變換,我們可以通過不斷旋轉 ‘ F ‘ `F` ‘F‘來得到更加稀疏的因子載荷矩陣,從而提高因子的可解釋性。
如何旋轉?怎么衡量旋轉后因子載荷矩陣的優良性?
令:
d i j 2 = a i j 2 h i 2 , i = 1 , 2 , … , m , j = 1 , 2 , … , n d_{ij}^2=\frac{a_{ij}^2}{h_i^2},\quad i=1,2,\dots,m,\;j=1,2,\dots,n dij2?=hi2?aij2??,i=1,2,…,m,j=1,2,…,n
‘ d i j 2 ‘ `d_{ij}^2` ‘dij2?‘衡量了因子 ‘ j ‘ `j` ‘j‘對觀測變量 ‘ X i ‘ `\mathbf{X}_i` ‘Xi?‘的影響,且消除了 ‘ a i j ‘ `a_{ij}` ‘aij?‘的正負號帶來的差異和各觀測變量在因子載荷大小上的不同帶來的差異。定義第 ‘ j ‘ `j` ‘j‘列 ‘ p ‘ `p` ‘p‘個數據 ‘ d i j 2 , i = 1 , 2 , … , m ‘ `d_{ij}^2,\;i=1,2,\dots,m` ‘dij2?,i=1,2,…,m‘的方差為:
V j = 1 m ∑ i = 1 m ( d i j 2 ? d ˉ j ) 2 = 1 m ∑ i = 1 m ( d i j 2 ? 1 p ∑ i = 1 m d i j 2 ) = 1 m [ ∑ i = 1 m d i j 4 ? m 1 m 2 ( ∑ i = 1 m d i j 2 ) 2 ] = 1 m 2 [ m ∑ i = 1 m d i j 4 ? 1 m ( ∑ i = 1 m d i j 2 ) 2 ] = 1 m 2 [ m ∑ i = 1 m a i j 4 h i 4 ? 1 m ( ∑ i = 1 m a i j 2 h i 2 ) 2 ] \begin{aligned} V_j&=\frac{1}{m}\sum_{i=1}^{m}(d_{ij}^2-\bar{d}_j)^2=\frac{1}{m}\sum_{i=1}^{m}\left(d_{ij}^2-\frac{1}{p}\sum_{i=1}^{m}d_{ij}^2\right) \\ &=\frac{1}{m}\left[\sum_{i=1}^{m}d_{ij}^4-m\frac{1}{m^2}\left(\sum_{i=1}^{m}d_{ij}^2\right)^2\right] \\ &=\frac{1}{m^2}\left[m\sum_{i=1}^{m}d_{ij}^4-\frac{1}{m}\left(\sum_{i=1}^{m}d_{ij}^2\right)^2\right] \\ &=\frac{1}{m^2}\left[m\sum_{i=1}^{m}\frac{a_{ij}^4}{h_i^4}-\frac{1}{m}\left(\sum_{i=1}^{m}\frac{a_{ij}^2}{h_i^2}\right)^2\right] \end{aligned} Vj??=m1?i=1∑m?(dij2??dˉj?)2=m1?i=1∑m?(dij2??p1?i=1∑m?dij2?)=m1? ?i=1∑m?dij4??mm21?(i=1∑m?dij2?)2 ?=m21? ?mi=1∑m?dij4??m1?(i=1∑m?dij2?)2 ?=m21? ?mi=1∑m?hi4?aij4???m1?(i=1∑m?hi2?aij2??)2 ??
若 ‘ V j ‘ `V_j` ‘Vj?‘越大,則第 ‘ j ‘ `j` ‘j‘個因子對觀測變量的影響越集中。定義因子載荷矩陣 ‘ A ‘ `A` ‘A‘的方差為:
V = ∑ j = 1 n V j = 1 m 2 { ∑ j = 1 n [ m ∑ i = 1 m a i j 4 h i 4 ? 1 m ( ∑ i = 1 m a i j 2 h i 2 ) 2 ] } V=\sum_{j=1}^{n}V_j=\frac{1}{m^2}\left\{\sum_{j=1}^{n}\left[m\sum_{i=1}^{m}\frac{a_{ij}^4}{h_i^4}-\frac{1}{m}\left(\sum_{i=1}^{m}\frac{a_{ij}^2}{h_i^2}\right)^2\right]\right\} V=j=1∑n?Vj?=m21?? ? ??j=1∑n? ?mi=1∑m?hi4?aij4???m1?(i=1∑m?hi2?aij2??)2 ?? ? ??
若 ‘ V ‘ `V` ‘V‘越大,則表明因子對觀測變量的影響越集中。
綜上,我們只需使得旋轉后得到的因子載荷矩陣 ‘ A ‘ `A` ‘A‘的方差 ‘ V ‘ `V` ‘V‘達到最大即可。
模型檢驗
由上面的討論可以看出,潛在因子的數目是一個超參數,也是一個非常重要的參數,我們該如何選擇呢?有沒有什么辦法能夠確定這一超參數的值?
在正態性假設下(仍需假設 X \mathbf{X} X是 n n n維正態隨機向量),我們可以對求解后的因子分析模型進行似然比檢驗。
設樣本數為 p p p,分別為 X 1 , X 2 , … , X p \mathbf{X_1},\mathbf{X_2},\dots,\mathbf{X_p} X1?,X2?,…,Xp?,都獨立同分布于 N ? m ( μ , Σ ) \operatorname{N}_m(\boldsymbol{\mu},\Sigma) Nm?(μ,Σ)。構建似然比檢驗假設:
H 0 : Σ = A A T + D , H 1 : Σ 為其它任一正定矩陣 \begin{equation*} H_0:\Sigma=AA^T+D,\quad H_1:\Sigma\text{為其它任一正定矩陣} \end{equation*} H0?:Σ=AAT+D,H1?:Σ為其它任一正定矩陣?
由\cref{def:MultiNormalPDF2}可得此時備擇假設下的對數似然函數為:
L 1 = ∑ i = 1 p ln ? { 1 ( 2 π ) m 2 ( det ? Σ ) 1 2 e ? 1 2 tr ? [ ( X i ? μ ) ( X i ? μ ) T Σ ? 1 ] } = ∑ i = 1 p { ? m 2 ln ? ( 2 π ) ? 1 2 ln ? ( det ? Σ ) ? 1 2 tr ? [ ( X i ? μ ) ( X i ? μ ) T Σ ? 1 ] } = ? p 2 { m ln ? ( 2 π ) + ln ? ( det ? Σ ) + 1 p ∑ i = 1 p tr ? [ ( X i ? μ ) ( X i ? μ ) T Σ ? 1 ] } \begin{align*} L_1&=\sum_{i=1}^{p}\ln\left\{\frac{1}{(2\pi)^{\frac{m}{2}}(\det\Sigma)^{\frac{1}{2}}}e^{-\frac{1}{2}\operatorname{tr}[(\mathbf{X_i}-\boldsymbol{\mu})(\mathbf{X_i}-\boldsymbol{\mu})^T\Sigma^{-1}]}\right\} \\ &=\sum_{i=1}^{p}\left\{-\frac{m}{2}\ln(2\pi)-\frac{1}{2}\ln(\det\Sigma)-\frac{1}{2}\operatorname{tr}[(\mathbf{X_i}-\boldsymbol{\mu})(\mathbf{X_i}-\boldsymbol{\mu})^T\Sigma^{-1}]\right\} \\ &=-\frac{p}{2}\left\{m\ln(2\pi)+\ln(\det\Sigma)+\frac{1}{p}\sum_{i=1}^{p}\operatorname{tr}[(\mathbf{X_i}-\boldsymbol{\mu})(\mathbf{X_i}-\boldsymbol{\mu})^T\Sigma^{-1}]\right\} \end{align*} L1??=i=1∑p?ln{(2π)2m?(detΣ)21?1?e?21?tr[(Xi??μ)(Xi??μ)TΣ?1]}=i=1∑p?{?2m?ln(2π)?21?ln(detΣ)?21?tr[(Xi??μ)(Xi??μ)TΣ?1]}=?2p?{mln(2π)+ln(detΣ)+p1?i=1∑p?tr[(Xi??μ)(Xi??μ)TΣ?1]}?
上式可以化為(樣本因子分析時需要注意使用協方差矩陣的無偏估計,系數應取 p ? 1 p-1 p?1):
L 1 ( Σ ) = ? p 2 [ m ln ? ( 2 π ) + ln ? ( det ? Σ ) + tr ? ( S Σ ? 1 ) ] \begin{equation*} L_1(\Sigma)=-\frac{p}{2}\left[m\ln(2\pi)+\ln(\det\Sigma)+\operatorname{tr}(S\Sigma^{-1})\right] \end{equation*} L1?(Σ)=?2p?[mln(2π)+ln(detΣ)+tr(SΣ?1)]?
其中 S S S為樣本協方差陣。這個似然函數在 Σ = S \Sigma=S Σ=S時取最大值,于是:
L 1 = ? p 2 [ m ln ? ( 2 π ) + ln ? ( det ? S ) + p ] \begin{equation*} L_1=-\frac{p}{2}[m\ln(2\pi)+\ln(\det S)+p] \end{equation*} L1?=?2p?[mln(2π)+ln(detS)+p]?
同理,此時原假設下的似然函數值為:
L 2 = ? p 2 [ m ln ? ( 2 π ) + ln ? ( det ? Σ ^ ) + tr ? ( S Σ ^ ) ] \begin{equation*} L_2=-\frac{p}{2}\left[m\ln(2\pi)+\ln(\det\hat{\Sigma})+\operatorname{tr}(S\hat{\Sigma})\right] \end{equation*} L2?=?2p?[mln(2π)+ln(detΣ^)+tr(SΣ^)]?
其中 Σ ^ = A ^ A ^ T ? D ^ \hat{\Sigma}=\hat{A}\hat{A}^T-\hat{D} Σ^=A^A^T?D^。
由似然比檢驗原理可知:
? 2 [ L 2 ( Σ ) ? L 1 ( Σ ) ] = ~ χ d f 2 \begin{equation*} -2[L_2(\Sigma)-L_1(\Sigma)]=\sim\chi^2_{df} \end{equation*} ?2[L2?(Σ)?L1?(Σ)]=~χdf2??
當
p [ ln ? ( det ? Σ ^ ) + tr ? ( S Σ ^ ) ? ln ? ( det ? S ) ? p ] > χ 0.95 2 ( d f ) \begin{equation*} p[\ln(\det\hat{\Sigma})+\operatorname{tr}(S\hat{\Sigma})-\ln(\det S)-p]>\chi^2_{0.95}(df) \end{equation*} p[ln(detΣ^)+tr(SΣ^)?ln(detS)?p]>χ0.952?(df)?
時應拒絕原假設,即 n n n個因子不足以解釋數據,應增大因子個數。其中 χ 0.95 2 ( d f ) \chi^2_{0.95}(df) χ0.952?(df)為分布的 0.95 0.95 0.95分位數。
不介紹自由度的計算,有點復雜,R語言factanal函數中使用的自由度計算公式為dof <- 0.5 * ((p - factors)^2 - p - factors)
,其中p
為變量數,factors
為潛在因子數。
因子得分
在擬合得到因子載荷矩陣后,我們可以反過來求解各樣本因子的取值,這樣一來就可以根據因子值去進行進一步的分析。例如在開頭的例子里,我們可以得到每個學生語言能力與理科能力的值,進而可以進行分類或選擇。
因子得分主要有兩種計算方式。
加權最小二乘法
考慮加權最小二乘函數:
φ ( F ) = ( X ? μ ? A F ) T D ? 1 ( X ? μ ? A F ) \begin{equation*} \varphi(F)=(\mathbf{X}-\boldsymbol{\mu}-AF)^TD^{-1}(\mathbf{X}-\boldsymbol{\mu}-AF) \end{equation*} φ(F)=(X?μ?AF)TD?1(X?μ?AF)?
求:
F ^ = arg ? min ? φ ( F ) \begin{equation*} \hat{F}=\arg\min\varphi(F) \end{equation*} F^=argminφ(F)?
由極值的必要條件得到:
? φ ( F ) ? F = ? 2 A T D ? 1 ( X ? μ ? A F ) = 0 A T D ? 1 ( X ? μ ) = A T D ? 1 A F F = ( A T D ? 1 A ) ? 1 A T D ? 1 ( X ? μ ) \begin{gather*} \frac{\partial\varphi(F)}{\partial F}=-2A^TD^{-1}(\mathbf{X}-\boldsymbol{\mu}-AF)=0 \\ A^TD^{-1}(\mathbf{X}-\boldsymbol{\mu})=A^TD^{-1}AF \\ F=(A^TD^{-1}A)^{-1}A^TD^{-1}(\mathbf{X}-\boldsymbol{\mu}) \end{gather*} ?F?φ(F)?=?2ATD?1(X?μ?AF)=0ATD?1(X?μ)=ATD?1AFF=(ATD?1A)?1ATD?1(X?μ)?
需要注意 A T D ? 1 A A^TD^{-1}A ATD?1A的可逆性。\par
若認為 X ~ N ? m ( μ + A F , D ) \mathbf{X}\sim\operatorname{N}_m(\boldsymbol{\mu}+AF,D) X~Nm?(μ+AF,D),則上述解得的 F F F也是極大似然估計的結果。
稱加權最小二乘法得到的因子得分為Bartlett因子得分。
從求解過程可以看出,該方法實際上是對特殊方差更大的變量施以更寬容的殘差值。
回歸法
設:
f j = ∑ i = 1 m b j i X i + ε j , Cov ? ( X i , ε j ) = 0 , j = 1 , 2 , … , n \begin{equation*} f_j=\sum_{i=1}^{m}b_{ji}\mathbf{X}_i+\varepsilon_j,\;\operatorname{Cov}(\mathbf{X}_i,\varepsilon_j)=0,\quad j=1,2,\dots,n \end{equation*} fj?=i=1∑m?bji?Xi?+εj?,Cov(Xi?,εj?)=0,j=1,2,…,n?
由\cref{prop:FactorAnalysis}(4)和\cref{prop:CovMat}(3)(5)可知:
a i j = Cov ? ( X i , f j ) = Cov ? ( X i , ∑ k = 1 m b j k X k + ε j ) = ∑ k = 1 m σ i k b j k \begin{equation*} a_{ij}=\operatorname{Cov}(\mathbf{X}_i,f_j)=\operatorname{Cov}\left(\mathbf{X}_i,\sum_{k=1}^{m}b_{jk}\mathbf{X}_k+\varepsilon_j\right)=\sum_{k=1}^{m}\sigma_{ik}b_{jk} \end{equation*} aij?=Cov(Xi?,fj?)=Cov(Xi?,k=1∑m?bjk?Xk?+εj?)=k=1∑m?σik?bjk??
令 B = ( b i j ) B=(b_{ij}) B=(bij?),則有:
A = Σ B T \begin{equation*} A=\Sigma B^T \end{equation*} A=ΣBT?
于是 B = A T Σ ? 1 B=A^T\Sigma^{-1} B=ATΣ?1,需要注意 Σ \Sigma Σ的可逆性。回歸法的因子得分即為:
F = A T Σ ? 1 X \begin{equation*} F=A^T\Sigma^{-1}\mathbf{X} \end{equation*} F=ATΣ?1X?
代碼實現
R語言中使用Factanal函數進行因子分析,注意它使用極大似然估計法進行參數估計。
factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,subset, na.action, start = NULL,scores = c("none", "regression", "Bartlett"),rotation = "varimax", control = NULL, ...)
- x:formula或矩陣或數據框。formula時要用data參數來指定數據來源,因為因子分析是對變量的考察,所以formula應形如
~X_1+X_2+X_3
,~
左側為空,右端指定觀測變量。若為矩陣或數據框,其行數應為樣本數,列數應為觀測變量數。若covmat
被指定,則不應指定x
。 - factors:指定潛在因子個數,應為正整數。
- data:
x
為formula時指定數據來源。 - covmat:指定觀測變量協方差矩陣,被指定時不需要傳遞
x, data
。 - n.obs:樣本數,指定
covmat
時需要明確樣本數。 - subset, na.action:選擇樣本子集及缺失值處理
- start:特殊方差迭代時的初始值,無需深究,除非你懂了極大似然估計法估計參數。
- scores:指定因子得分計算方式,
"none"
表示不計算因子得分。 - rotation:因子旋轉,值只能為
"varimax"
或none
,前者表示使用varimax
函數來進行正交旋轉,后者表示不進行因子旋轉。
注意事項
factanal()
函數只提供正交因子分析,斜交因子分析請另外找包。- 若只提供
covmat
,無法進行似然比檢驗,結果將會輸出擬合殘差平方和與自由度,但如果此時提供n.obs
,則會輸出似然比檢驗的結果。
例子
x<-c(1.000,0.923, 1.000,0.841, 0.851, 1.000,0.756, 0.807, 0.870, 1.000,0.700, 0.775, 0.835, 0.918, 1.000,0.619, 0.695, 0.779, 0.864, 0.928, 1.000,0.633, 0.697, 0.787, 0.869, 0.935, 0.975, 1.000,0.520, 0.596, 0.705, 0.806, 0.866, 0.932, 0.943, 1.000)
names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")
R<-matrix(0, nrow=8, ncol=8, dimnames=list(names, names))
for (i in 1:8){for (j in 1:i){R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]}
}
factanal(factors = 3, covmat = R, n.obs = 20)Call:
factanal(factors = 3, covmat = R, n.obs = 20)Uniquenesses:X1 X2 X3 X4 X5 X6 X7 X8
0.005 0.106 0.144 0.080 0.057 0.038 0.010 0.090 Loadings:Factor1 Factor2 Factor3
X1 0.284 0.956
X2 0.395 0.851 0.119
X3 0.543 0.723 0.196
X4 0.682 0.595 0.318
X5 0.797 0.501 0.240
X6 0.900 0.382
X7 0.914 0.391
X8 0.913 0.273 Factor1 Factor2 Factor3
SS loadings 4.110 3.139 0.221
Proportion Var 0.514 0.392 0.028
Cumulative Var 0.514 0.906 0.934Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 0.79 on 7 degrees of freedom.
The p-value is 0.998
- Uniquenesses即為特殊方差的值
- loadings表示因子載荷,其中較小的因子載荷會顯示空值
- SS loadings為 g j 2 g_j^2 gj2?,Proportion Var為因子對所有觀測變量總方差的貢獻率,即 g j 2 / ∑ i = 1 m Var ? ( X i ) g_j^2/\sum\limits_{i=1}^m\operatorname{Var}(\mathbf{X}_i) gj2?/i=1∑m?Var(Xi?),Cumulative Var為累計方差貢獻率。
- 最底下是似然比檢驗的結果,若p值小于0.05,則應增大潛在因子數。
Reference
1.薛毅,統計建模與R軟件