因子分析——數學原理及R語言代碼

正交因子分析

  • 目的
  • 數學原理
    • 參數估計方法
      • 主成分法
      • 主因子法
      • 極大似然法
    • 因子旋轉
    • 模型檢驗
    • 因子得分
      • 加權最小二乘法
      • 回歸法
  • 代碼實現
    • 注意事項
    • 例子
  • 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?的表達式中。


上述因子分析模型具有如下性質:

  1. ‘ Σ = A A T + D ‘ `\Sigma=AA^T+D` ‘Σ=AAT+D

  2. ‘ 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+=μ?+A?F+ε?

  3. 因子載荷不唯一;

  4. ‘ 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?

  5. ‘ h i 2 = ∑ j = 1 n a i j 2 ‘ `h_i^2=\sum\limits_{j=1}^{n}a_{ij}^2` hi2?=j=1n?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=1n?aij2?+σi2?=hi2?+σi2?,i=1,2,,m

  6. ‘ g j 2 = ∑ i = 1 m a i j 2 ‘ `g_j^2=\sum\limits_{i=1}^{m}a_{ij}^2` gj2?=i=1m?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=1m?Var(Xi?)=j=1n?gj2?+i=1n?σ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=1m?Var(Xi?)?=tr[Cov(X)]=tr(AAT+D)=i=1m?j=1n?aij2?+i=1n?σi2?=j=1n?i=1m?aij2?+i=1n?σi2?=j=1n?gj2?+i=1n?σ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=1m?λ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=1m?λi?li?liT?i=1n?λ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=1n?λi?)/(i=1m?λ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=1n?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=1m?(dij2??dˉj?)2=m1?i=1m?(dij2??p1?i=1m?dij2?)=m1? ?i=1m?dij4??mm21?(i=1m?dij2?)2 ?=m21? ?mi=1m?dij4??m1?(i=1m?dij2?)2 ?=m21? ?mi=1m?hi4?aij4???m1?(i=1m?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=1n?Vj?=m21?? ? ??j=1n? ?mi=1m?hi4?aij4???m1?(i=1m?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=1p?ln{(2π)2m?(detΣ)21?1?e?21?tr[(Xi??μ)(Xi??μ)TΣ?1]}=i=1p?{?2m?ln(2π)?21?ln(detΣ)?21?tr[(Xi??μ)(Xi??μ)TΣ?1]}=?2p?{mln(2π)+ln(detΣ)+p1?i=1p?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) XNm?(μ+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=1m?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=1m?bjk?Xk?+εj?)=k=1m?σ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, ...)
  1. x:formula或矩陣或數據框。formula時要用data參數來指定數據來源,因為因子分析是對變量的考察,所以formula應形如~X_1+X_2+X_3~左側為空,右端指定觀測變量。若為矩陣或數據框,其行數應為樣本數,列數應為觀測變量數。若covmat被指定,則不應指定x
  2. factors:指定潛在因子個數,應為正整數。
  3. datax為formula時指定數據來源。
  4. covmat:指定觀測變量協方差矩陣,被指定時不需要傳遞x, data
  5. n.obs:樣本數,指定covmat時需要明確樣本數。
  6. subset, na.action:選擇樣本子集及缺失值處理
  7. start:特殊方差迭代時的初始值,無需深究,除非你懂了極大似然估計法估計參數。
  8. scores:指定因子得分計算方式,"none"表示不計算因子得分。
  9. rotation:因子旋轉,值只能為"varimax"none,前者表示使用varimax函數來進行正交旋轉,后者表示不進行因子旋轉。

注意事項

  1. factanal()函數只提供正交因子分析,斜交因子分析請另外找包。
  2. 若只提供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 
  1. Uniquenesses即為特殊方差的值
  2. loadings表示因子載荷,其中較小的因子載荷會顯示空值
  3. 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=1m?Var(Xi?),Cumulative Var為累計方差貢獻率。
  4. 最底下是似然比檢驗的結果,若p值小于0.05,則應增大潛在因子數。

Reference

1.薛毅,統計建模與R軟件

本文來自互聯網用戶投稿,該文觀點僅代表作者本人,不代表本站立場。本站僅提供信息存儲空間服務,不擁有所有權,不承擔相關法律責任。
如若轉載,請注明出處:http://www.pswp.cn/bicheng/80733.shtml
繁體地址,請注明出處:http://hk.pswp.cn/bicheng/80733.shtml
英文地址,請注明出處:http://en.pswp.cn/bicheng/80733.shtml

如若內容造成侵權/違法違規/事實不符,請聯系多彩編程網進行投訴反饋email:809451989@qq.com,一經查實,立即刪除!

相關文章

ACL訪問控制列表:access-list 10 permit 192.168.10.1

ACL訪問控制列表 標準ACL語法 1. 創建ACL access-list <編號> <動作> <源IP> <通配符掩碼> // 編號范圍 1-99 // 動作&#xff1a;permit 允許 、 deny 拒絕2. 示例 //允許192.168.1.0/24g整個網絡,0.0.0.255 反掩碼 access-list 10 permit 192.1…

解決社區錄音應用橫屏狀態下,錄音后無法播放的bug

最近看到社區有小伙伴反映&#xff0c;社區錄音應用橫屏時&#xff0c;錄音后無法播放的問題。現分享解決辦法。 社區錄音應用的來源&#xff1a;https://gitee.com/openharmony/applications_app_samples/tree/OpenHarmony-5.0.2-Release/code/SystemFeature/Media/Recorder …

每周靶點分享:Angptl3、IgE、ADAM9及文獻分享:抗體的多樣性和特異性以及結構的新見解

本期精選了《脂質代謝的關鍵調控者Angptl3》《T細胞活化抑制因子VISTA靶點》《文獻分享&#xff1a;雙特異性抗體重輕鏈配對設計》三篇文章。以下為各研究內容的概述&#xff1a; 1. 脂質代謝的關鍵調控者Angptl3 血管生成素相關蛋白3&#xff08;Angptl3&#xff09;是血管生…

保持Word中插入圖片的清晰度

大家有沒有遇到這個問題&#xff0c;原本繪制的高清晰度圖片&#xff0c;插入word后就變模糊了。先說原因&#xff0c;word默認啟動了自動壓縮圖片功能&#xff0c;分享一下如何關閉這項功能&#xff0c;保持Word中插入圖片的清晰度。 ①在Word文檔中&#xff0c;點擊左上角的…

Datawhale AI春訓營 day

待補充 2025星火杯應用賽入門應用 創空間 魔搭社區 {"default": {"system": "你是星火大模型&#xff0c;一個由科大訊飛研發的人工智能助手。請用簡潔、專業、友好的方式回答問題。","description": "默認系統提示詞"}…

項目全棧實戰-基于智能體、工作流、API模塊化Docker集成的創業分析平臺

目錄 思維導圖 前置知識 Docker是什么&#xff1f; Docker的核心概念&#xff1a; Docker在本項目中的作用 1. 環境隔離與一致性 2. 簡化部署流程 3. 資源管理與擴展性 4. 服務整合與通信 5. 版本控制和回滾 6. 開發與生產環境一致性 總結 前端 1.小程序 2.web …

正則表達式實用指南:原理、場景、優化與引擎對比

正則表達式實用指南&#xff1a;原理、場景、優化與引擎對比 正則表達式&#xff08;Regular Expression&#xff0c;簡稱 regex 或 regexp&#xff09;是程序員處理文本數據時不可或缺的“瑞士軍刀”。無論是表單校驗、日志分析、數據清洗&#xff0c;還是敏感信息脫敏&#…

OSCP - Hack The Box - Sau

主要知識點 CVE-2023-27163漏洞利用systemd提權 具體步驟 執行nmap掃描&#xff0c;可以先看一下55555端口 Nmap scan report for 10.10.11.224 Host is up (0.58s latency). Not shown: 65531 closed tcp ports (reset) PORT STATE SERVICE VERSION 22/tcp o…

5.1.1 WPF中Command使用介紹

WPF 的命令系統是一種強大的輸入處理機制,它比傳統的事件處理更加靈活和可重用,特別適合 MVVM (Model, View, ViewModel)模式開發。 一、命令系統核心概念 1.命令系統基本元素: 命令(Command): 即ICommand類,使用最多的是RoutedCommand,也可以自己繼承ICommand使用自定…

Dagster Pipes系列-2:增強外部腳本與Dagster的交互能力

在現代數據工程中&#xff0c;自動化和監控是確保數據管道高效運行的關鍵因素。Dagster作為一款強大的數據編排工具&#xff0c;提供了多種方式來實現這些目標。本文將深入探討如何使用Dagster Pipes修改外部代碼&#xff0c;以實現日志記錄、結構化元數據報告以及資產檢查等功…

C++類和對象進階 —— 與數據結構的結合

&#x1f381;個人主頁&#xff1a;工藤新一 &#x1f50d;系列專欄&#xff1a;C面向對象&#xff08;類和對象篇&#xff09; &#x1f31f;心中的天空之城&#xff0c;終會照亮我前方的路 &#x1f389;歡迎大家點贊&#x1f44d;評論&#x1f4dd;收藏?文章 文章目錄 […

Java中進階并發編程

第一章、并發編程的挑戰 并發和并行&#xff1a;指多線程或多進程 線程的本質&#xff1a;操作系統能夠進行運算調度的最小單位&#xff0c;是進程&#xff08;Process&#xff09;中的實際工作單元 進程的本質&#xff1a;操作系統進行資源分配和調度的基本單位&#xff0c…

《 指針變量類型與內存訪問:揭秘背后的奧秘》

&#x1f680;個人主頁&#xff1a;BabyZZの秘密日記 &#x1f4d6;收入專欄&#xff1a;C語言 &#x1f30d;文章目入 一、指針變量類型的基本概念二、指針類型與內存訪問字節數的關系&#xff08;一&#xff09;整型指針&#xff08;二&#xff09;字符型指針&#xff08;三&…

mapbox進階,使用mapbox-plugins插件加載餅狀圖

????? 主頁: gis分享者 ????? 感謝各位大佬 點贊?? 收藏? 留言?? 加關注?! ????? 收錄于專欄:mapbox 從入門到精通 文章目錄 一、??前言1.1 ??mapboxgl.Map 地圖對象1.1 ??mapboxgl.Map style屬性二、??使用mapbox-plugins插件加載餅狀圖1. ?…

GraphicLayer與BusineDataLayer層級控制

補充說明&#xff1a; 當參與層級控制的元素是點型元素時&#xff0c;是無法參與ZIndex層級控制的&#xff0c;此時可以換個解決方案 1.給不同的高度值實現&#xff0c;元素間的層級控制覆蓋 import * as mars3d from "mars3d"export let map // mars3d.Map三維地…

uniapp 百家云直播插件打包失敗

打包錯誤日志 Android自有證書 打包失敗 錯誤日志: https://app.liuyingyong.cn/build/errorLog/cf41a610-effe-11ef-88db-05262d4c3e5d原因&#xff1a;需要導入插件依賴 依賴地址&#xff1a;https://ext.dcloud.net.cn/plugin?id16289 百家云直播插件地址 直播插…

【C++】”如虎添翼“:模板初階

泛型編程&#xff1a; C中一種使用模板來實現代碼重用和類型安全的編程范式。它允許程序員編寫與數據類型無關的代碼&#xff0c;從而可以用相同的代碼邏輯處理不同的數據類型。模板是泛型編程的基礎 模板分為兩類&#xff1a; 函數模板&#xff1a;代表了一個函數家族&#x…

十五、多態與虛函數

十五、多態與虛函數 15.1 引言 面向對象編程的基本特征&#xff1a;數據抽象&#xff08;封裝&#xff09;、繼承、多態基于對象&#xff1a;我們創建類和對象&#xff0c;并向這些對象發送消息多態&#xff08;Polymorphism&#xff09;&#xff1a;指的是相同的接口、不同的…

點云特征提取的兩大經典范式:Voxel-based 與 Pillar-based

點云特征提取的兩大經典范式&#xff1a;Voxel-based 與 Pillar-based 在點云處理領域&#xff0c;尤其是針對 3D 目標檢測任務&#xff0c;特征提取是核心環節之一。目前&#xff0c;Voxel-based&#xff08;體素化&#xff09;和 Pillar-based&#xff08;柱狀化&#xff09…

前蘋果首席設計官回顧了其在蘋果的設計生涯、公司文化、標志性產品的背后故事

每周跟蹤AI熱點新聞動向和震撼發展 想要探索生成式人工智能的前沿進展嗎&#xff1f;訂閱我們的簡報&#xff0c;深入解析最新的技術突破、實際應用案例和未來的趨勢。與全球數同行一同&#xff0c;從行業內部的深度分析和實用指南中受益。不要錯過這個機會&#xff0c;成為AI領…