因素分析

本文的閱讀等級:高級

因素分析 (factor analysis) 是統計學中一種多變量分析法。因素分析與主成分分析具有一些相同的概念與技巧,但兩者的建模推理方向相反。假設可量測的隨機向量 \mathbf{x}=(x_1,\ldots,x_p)^T 服從一個未知的機率分布 p(\mathbf{x}),期望值為 \hbox{E}[\mathbf{x}]=\boldsymbol{\mu}=(\mu_1,\ldots,\mu_p)^T,共變異數矩陣為 \hbox{cov}[\mathbf{x}]=\Sigma=[\sigma_{ij}]1\le i,j\le p。主成分分析的主要功用是降維 (dimension reduction),我們從原始的變數 x_1,\ldots,x_p 構築一組新變數 z_1,\ldots,z_k1\le k<p。具體地說,低維隨機向量 \mathbf{z}=(z_1,\ldots,z_k)^T 由離差 (deviation) \mathbf{x}-\boldsymbol{\mu} 的線性映射產生:

\displaystyle  \mathbf{z}=W^T(\mathbf{x}-\boldsymbol{\mu})

其中 W 是一個 p\times k 階矩陣滿足 W^TW=I_k (見“主成分分析”)。在因素分析,我們設想隨機向量 \mathbf{x} 的資料生成模型 (generative model) 如下:

\displaystyle   \mathbf{x}=\boldsymbol{\mu}+F\mathbf{z}+\boldsymbol{\epsilon}

其中 \mathbf{z}=(z_1,\ldots,z_k)^T 是一組無法量測的隱藏變數,稱為隱藏因素 (hidden factor)、共同因素 (common factor) 或簡稱因素,F 是一個 p\times k 階變換矩陣[1]\boldsymbol{\epsilon}=(\epsilon_1,\ldots,\epsilon_p)^T 是代表雜音的隨機向量。本文討論的問題包括:

  • 因素分析如何描述多隨機變數的產生?
  • 如何估計因素分析的模型參數?
  • 因素分析如何解釋隱藏因素的涵義?
  • 因素分析如何應用於降維?
  • 因素分析與主成分分析有哪些相同與相異的性質?

模型

給定可量測的隨機變數 x_1,\ldots,x_p,因素分析假設每個隨機變數 x_i 的離差 x_i-\mu_i 為因素 z_1,\ldots,z_k1\le k<p,的線性組合加上殘差項 \epsilon_i,如下:

\displaystyle \begin{aligned}  x_1-\mu_1&=f_{11}z_1+\cdots+f_{1k}z_k+\epsilon_1,\\  &\vdots\\  x_p-\mu_p&=f_{p1}z_1+\cdots+f_{pk}z_k+\epsilon_p,  \end{aligned}

或以向量─矩陣形式表示為

\mathbf{x}-\boldsymbol{\mu}=F\mathbf{z}+\boldsymbol{\epsilon}

其中 F=[f_{ij}] 是一個 p\times k 階常數矩陣,稱為因素負荷 (factor loading),描述因素 \mathbf{z} 至多變數 \mathbf{x} 的變換。你不妨想像 z_1,\ldots,z_k 代表未知的信號源,將 k 個信號同時送入 p 個通道,每個通道各自有增強或衰減信號的方式 (線性組合),通道末端再引入雜音 (殘差) 就是我們接收到的變數 x_1,\ldots,x_p

為便利分析,假設因素 \mathbf{z}=(z_1,\ldots,z_k)^T 與殘差 \boldsymbol{\epsilon}=(\epsilon_1,\ldots,\epsilon_p)^T 具有底下的機率性質:

  1. \boldsymbol{\mu}=\mathbf{0},若不成立,直接從 \mathbf{x} 減去 \hbox{E}[\mathbf{x}]
  2. \hbox{E}[\mathbf{z}]=\mathbf{0}\hbox{cov}[\mathbf{z}]=I,即每個信號的強度相同 (\hbox{var}[z_j]=1) 且信號彼此獨立 (\hbox{cov}[z_i,z_j]=0i\neq j)。
  3. \hbox{E}[\boldsymbol{\epsilon}]=\mathbf{0}\hbox{cov}[\boldsymbol{\epsilon}]=\Psi=\hbox{diag}(\psi_1,\ldots,\psi_p),即雜音的強度為 \hbox{var}[\epsilon_i]=\psi_i 且雜音彼此獨立 (\hbox{cov}[\epsilon_i,\epsilon_j]=0i\neq j)。
  4. \hbox{cov}[\epsilon_i,z_j]=0i=1,\ldots,pj=1,\ldots,k,即雜音與信號完全無關。

因素分析的模型參數包括 p\times k 階因素負荷 F 與殘差的變異數 \Psi=\hbox{diag}(\psi_1,\ldots,\psi_p),合計 p(k+1) 個參數值。根據上述假設,你可以計算隨機向量 \mathbf{x} 的期望值與共變異數矩陣。期望值 \hbox{E}[\cdot] 是線性算子,因此

\displaystyle   \hbox{E}[\mathbf{x}]=\hbox{E}[F\mathbf{z}+\boldsymbol{\epsilon}]=F\hbox{E}[\mathbf{z}]+\hbox{E}[\boldsymbol{\epsilon}]=\mathbf{0}

使用此結果與前述假設,\mathbf{x} 的共變異數矩陣 \Sigma=[\sigma_{ij}] 計算如下 (見“共變異數矩陣的性質”):

\displaystyle  \begin{aligned}  \Sigma&=\hbox{cov}[\mathbf{x}]=\hbox{E}[\mathbf{x}\mathbf{x}^T]\\  &=\hbox{E}\left[(F\mathbf{z}+\boldsymbol{\epsilon})(F\mathbf{z}+\boldsymbol{\epsilon})^T\right]\\  &=\hbox{E}\left[(F\mathbf{z}\mathbf{z}^TF^T+\boldsymbol{\epsilon}\mathbf{z}^TF^T+F\mathbf{z}\boldsymbol{\epsilon}^T+\boldsymbol{\epsilon}\boldsymbol{\epsilon}^T\right]\\  &=F\hbox{E}[\mathbf{z}\mathbf{z}^T]F^T+\hbox{E}[\boldsymbol{\epsilon}\mathbf{z}^T]F^T+F\hbox{E}[\mathbf{z}\boldsymbol{\epsilon}^T]+\hbox{E}[\boldsymbol{\epsilon}\boldsymbol{\epsilon}^T]\\  &=F\hbox{cov}[\mathbf{z}]F^T+\hbox{cov}[\boldsymbol{\epsilon}]\\  &=FF^T+\Psi.  \end{aligned}

乘開上式可得 x_i 的變異數

\displaystyle  \hbox{var}[x_i]=\sigma_{ii}=\sum_{j=1}^kf_{ij}^2+\psi_i

其中包含兩個部分:(FF^T)_{ii}=\sum_{j=1}^kf_{ij}^2 稱為 x_i 的共同性 (communality),意指透過共同的因素,x_i 與其他變數共享的變異成分;\psi_i 稱為 x_i 的特殊變異數 (specific variance),意指 x_i 自己獨有不與其他變數共享的變異成分。因為 \Psi 是對角矩陣,x_ix_j 的共變異數為 (FF^T)_{ij},即

\displaystyle  \hbox{cov}[x_i,x_j]=\sigma_{ij}=\sum_{l=1}^kf_{il}f_{jl}

因素負荷 F 的元 f_{ij} 是變數 x_i 與因素 z_j 的共變異數,如下:

\displaystyle  \begin{aligned}  \hbox{cov}[x_i,z_j]&=\hbox{cov}\left[\sum_{l=1}^kf_{il}z_l+\epsilon_i,z_j\right]=\hbox{E}\left[\left(\sum_{l=1}^kf_{il}z_l+\epsilon_i\right)z_j\right]\\  &=\sum_{l=1}^kf_{il}\hbox{E}[z_lz_j]=\sum_{l=1}^kf_{il}\,\hbox{cov}[z_l,z_j]=f_{ij}.  \end{aligned}

因素分析並不具有唯一的因素負荷 F。令 Q 是一個 k\times k 階正交矩陣 (orthogonal matrix),Q^TQ=QQ^T=I。將上式嵌入因素分析模型,

\mathbf{x}=FQQ^T\mathbf{z}+\boldsymbol{\epsilon}=F'\mathbf{z}'+\boldsymbol{\epsilon}

其中 F'=FQ 是對應新因素 \mathbf{z}'=Q^T\mathbf{z} 的因素負荷。不難驗證 \mathbf{z}' 滿足前述假設條件:\hbox{E}[\mathbf{z}']=\mathbf{0}\hbox{cov}[\mathbf{z}']=I\hbox{cov}[\epsilon_i,z'_j]=0i=1,\ldots,pj=1,\ldots,k。共變異數矩陣不會因選擇的因素負荷而改變,

\Sigma=FF^T+\Psi=FQQ^TF^T+\Psi=F'(F')^T+\Psi

就解釋隨機向量 \mathbf{x} 的共變異數矩陣而言,因素 \mathbf{z} 配合負載 F 等價於因素 \mathbf{z}' 配合負載 F'。要排除因素分析的非唯一性,你可以加入一些額外條件,譬如,條件1:F^TF 是對角矩陣,或條件2:F^T\Psi^{-1}F 是對角矩陣。

參數估計

假設你有一筆樣本大小為 n 的數據 \{\mathbf{x}_1,\ldots,\mathbf{x}_n\},每個數據點視為一個向量 \mathbf{x}_l=(x_{l1},\ldots,x_{lp})^T\in\mathbb{R}^p。定義 n\times p 階數據矩陣

\displaystyle   X=\begin{bmatrix}  \mathbf{x}_1^T\\  \vdots\\  \mathbf{x}_n^T  \end{bmatrix}=\begin{bmatrix}  x_{11}&\cdots&x_{1p}\\  \vdots&\ddots&\vdots\\  x_{n1}&\cdots&x_{np}  \end{bmatrix}

其中 x_{li} 表示第 i 個隨機變數 x_i 的第 l 個量測值,i=1,\ldots,pl=1,\ldots,n。我們假設樣本已經去除平均數,\sum_{l=1}^n\mathbf{x}_l=\mathbf{0}。共變異數矩陣 \Sigma 通常以樣本共變異數矩陣估計:

\displaystyle  S=\frac{1}{n-1}\sum_{i=1}^n\mathbf{x}_i\mathbf{x}_i^T=\frac{1}{n-1}X^TX

接下來的問題是從給定的樣本共變異數矩陣 S=[s_{ij}] 求出因素分析的參數估計 \hat{F}\hat{\Psi} 使得 S\approx \hat{F}\hat{F}^T+\hat{\Psi}。主因素分析 (principal factor analysis) 與最大似然估計 (maximum likelihood estimation) 是兩種常見的估計法。

主因素分析

主因素分析運用主成分分析方法計算。我們預先估計 \hat{\Psi}=\hbox{diag}(\hat{\psi}_1,\ldots,\hat{\psi}_p)。估計特殊變異數 \hat{\psi}_i 的一個合理常用方式是對變數 x_i 建立線性回歸,以其餘的 p-1 個變數當作預測變數,如下:

\displaystyle  x_{li}=b_1x_{l1}+\cdots+b_{i-1}x_{l,i-1}+b_{i+1}x_{l,i+1}+\cdots+b_px_{lp}+e_{li},~~l=1,\ldots,n,

其中 e_{li} 為殘差。利用最小平方法求出參數 \{b_1,\ldots,b_p\}\setminus \{b_i\},可得最小的殘差平方和 \sum_{l=1}^n\hat{e}_{li}^2,再令 \hat{\psi}_i=\frac{1}{n-1}\sum_{l=1}^n\hat{e}_{li}^2

定義簡約 (reduced) 共變異數矩陣 P=S-\hat{\Psi},接下來要找出 \hat{F} 使得 P\approx \hat{F}\hat{F}^T。因為 P 是一個實對稱矩陣,寫出正交對角化表達式 (見“可對角化矩陣的譜分解”)

\displaystyle  P=UDU^T=\begin{bmatrix}  \mathbf{u}_1&\cdots&\mathbf{u}_p  \end{bmatrix}\begin{bmatrix}  d_1&&\\  &\ddots&\\  &&d_p  \end{bmatrix}\begin{bmatrix}  \mathbf{u}_1^T\\  \vdots\\  \mathbf{u}_p^T  \end{bmatrix}=\sum_{i=1}^pd_i\mathbf{u}_i\mathbf{u}_i^T

其中 U=\begin{bmatrix}  \mathbf{u}_1&\cdots&\mathbf{u}_p  \end{bmatrix}p\times p 階正交矩陣,U^T=U^{-1},且 D=\hbox{diag}(d_1,\ldots,d_p)d_1\ge \cdots\ge d_p。如果 d_1,\ldots,d_k>0 且其餘的 d_i 接近零,則

\displaystyle  P\approx\sum_{i=1}^kd_i\mathbf{u}_i\mathbf{u}_i^T=U_kD_kU_k^T=(U_kD_k^{1/2})(U_kD_k^{1/2})^T

上式定義了 U_k=\begin{bmatrix}  \mathbf{u}_1&\cdots&\mathbf{u}_k  \end{bmatrix}D_k=\hbox{diag}(d_1,\ldots,d_k)。因此,因素負荷的估計值為

\hat{F}=[\hat{f}_{ij}]=U_kD_k^{1/2}

使用這個結果,我們可以重新估計特殊變異數:

\displaystyle  \hat{\psi}_i=s_{ii}-\sum_{j=1}^k\hat{f}_{ij}^2,~~i=1,\ldots,p

重複上面兩個步驟交替計算 \hat{F}\hat{\Psi},直到滿足收斂條件才停止。因為 U_k^TU_k=I_k\hat{F} 滿足前述條件1,

\displaystyle  \hat{F}^T\hat{F}=(U_kD_k^{1/2})^T(U_kD_k^{1/2})=D_k^{1/2}U_k^TU_kD_k^{1/2}=D_k

如果特殊變異矩陣不存在,\Psi=0,因素分析模型 \mathbf{x}=F\mathbf{z} 退化為主成分分析 \mathbf{z}'=W^T\mathbf{x}。請注意,因素分析的因素 \mathbf{z} 其實與主成分 \mathbf{z}' 有不同的數值範圍,因為 W^TW=I_k,意思是 W 有單範正交 (orthonormal) 的行向量 (column vector)。從樣本共變異數矩陣的正交對角化表達式 S=V\Lambda V^T,其中 V^T=V^{-1}\Lambda=\hbox{diag}(\lambda_1,\ldots,\lambda_p)\lambda_1\ge\cdots\ge \lambda_p\ge 0,推得因素負荷估計 \hat{F}=V_k\Lambda^{1/2}_k,這裡 V_k\Lambda_k^{1/2} 的定義方式同 U_kD_k^{1/2}。因此,\hat{F}^T\hat{F}=\Lambda_k。若 \Lambda_k 是可逆的,則

\begin{aligned}  \mathbf{x}=\hat{F}\mathbf{z}&\Rightarrow \hat{F}^T\mathbf{x}=\hat{F}^T\hat{F}\mathbf{z}\\  &\Rightarrow \Lambda_k^{1/2}V_k^T\mathbf{x}=\Lambda_k\mathbf{z}\\  &\Rightarrow V_k^T\mathbf{x}=\Lambda_k^{1/2}\mathbf{z}.  \end{aligned}

最後一個式子表明主成分分析的變換矩陣為 W=V_k,滿足 W^TW=V_k^TV_k=I_k,主成分與因素的關係為 \mathbf{z}'=\Lambda_k^{1/2}\mathbf{z},即有 (見“共變異數矩陣的性質”)

\hbox{cov}[\mathbf{z'}]=\hbox{cov}[\Lambda_k^{1/2}\mathbf{z}]  =\Lambda_k^{1/2}\hbox{cov}[\mathbf{z}]\Lambda_k^{1/2}=\Lambda_k^{1/2}I_k\Lambda_k^{1/2}=\Lambda_k

主成分保留了最大的變異總量 \hbox{trace}(\hbox{cov}[\mathbf{z'}])=\hbox{trace}\Lambda_k=\lambda_1+\cdots+\lambda_k。因此,共同因素 \mathbf{z} 即為標準化的主成分 \mathbf{z}'z_j=z'_j/\sqrt{\lambda_j}。附帶一提,利用數據矩陣 X 的奇異值分解可求得 V\Lambda (見“主成分分析與低秩矩陣近似”)。

如果特殊變異矩陣是一個純量矩陣,\Psi=\psi I\psi>0,我們稱此模型為機率主成分分析 (probabilistic PCA),則

\displaystyle  \begin{aligned}  S&=V\Lambda V^T\\  &=V\begin{bmatrix}  \lambda_1-\psi&&\\  &\ddots&\\  &&\lambda_p-{\psi}  \end{bmatrix}V^T+V\begin{bmatrix}  {\psi}&&\\  &\ddots&\\  &&\psi\end{bmatrix}V^T\\  &\approx V\begin{bmatrix}  \Lambda_k-{\psi}I_k&0\\  0&0  \end{bmatrix}V^T+{\psi}I_p\\  &=V_k(\Lambda_k-{\psi}I_k)V_k^T+{\psi}I_p.  \end{aligned}

我們選擇估計值 \hat{\psi} 使得主對角元的誤差平方和 \sum_{i=k+1}^p(\lambda_i-{\psi})^2 為最小,即

\displaystyle  \hat{\psi}=\frac{1}{p-k}\sum_{i=k+1}^p\lambda_i

因素負荷則為

\hat{F}=V_k(\Lambda_k-\hat{\psi}I_k)^{1/2}=V_k\,\hbox{diag}\left(\sqrt{\lambda_1-\hat{\psi}},\ldots,\sqrt{\lambda_k-\hat{\psi}}\right)

最大似然估計

最大似然估計必須知道隨機向量 \mathbf{x} 的機率分布 p(\mathbf{x})。我們假設每一個數據點 \mathbf{x}_i 服從常態分布,\mathbf{x}_i\sim\mathcal{N}(\mathbf{0},\Sigma),其中 \Sigma=FF^T+\Psi,似然函數 \mathcal{L} 取對數 (見“多變量常態分布的最大似然估計”) 並充分利用跡數性質 (見“跡數的性質與應用”),可得

\displaystyle   \begin{aligned}  \log\mathcal{L}&=\log\prod_{i=1}^n\mathcal{N}(\mathbf{x}_i|\mathbf{0},\Sigma)\\  &=\log\prod_{i=1}^n(2\pi)^{-p/2}|\Sigma|^{-1/2}\exp\left\{-\frac{1}{2}\mathbf{x}_i^T\Sigma^{-1}\mathbf{x}_i\right\}\\  &=\sum_{i=1}^n\left(-\frac{p}{2}\log(2\pi)-\frac{1}{2}\log\vert \Sigma \vert-\frac{1}{2}\mathbf{x}_i^T\Sigma^{-1}\mathbf{x}_i\right)\\  &=-\frac{np}{2}\log(2\pi)-\frac{n}{2} \log\vert \Sigma\vert-\frac{1}{2}\sum_{i=1}^n\hbox{trace}\left(\mathbf{x}_i^T\Sigma^{-1}\mathbf{x}_i\right)\\  &=-\frac{np}{2}\log(2\pi)-\frac{n}{2} \log\vert \Sigma\vert-\frac{1}{2}\hbox{trace}\left(\sum_{i=1}^n\mathbf{x}_i\mathbf{x}_i^T\Sigma^{-1}\right)\\  &=-\frac{np}{2}\log(2\pi)-\frac{n}{2} \log\vert \Sigma\vert-\frac{n-1}{2}\hbox{trace}\left(S\Sigma^{-1}\right),  \end{aligned}

其中 |\Sigma|=\det \Sigma。接著求出 \hat{\Psi}\hat{F} 使 \log\mathcal{L} 有最大值。為簡化記號,令 \Sigma^{-1}=[\rho_{ij}]。計算 \log\mathcal{L}\psi_i1\le i\le p,的偏導數 (推導見[2]),

\displaystyle  \begin{aligned}  \frac{\partial \log\mathcal{L}}{\partial \psi_{i}}&=-\frac{n}{2}\frac{\partial \log|\Sigma|}{\partial\psi_i}-\frac{n-1}{2}\frac{\partial \,\hbox{trace}(S\Sigma^{-1})}{\partial \psi_i}\\  &=-\frac{n}{2}\rho_{ii}+\frac{n-1}{2}\sum_{l=1}^p\sum_{j=1}^p\rho_{il}s_{lj}\rho_{ji}.  \end{aligned}

設上式為零,以矩陣形式表示為

\displaystyle  \hbox{diag}(\Sigma^{-1})=\frac{n-1}{n}\hbox{diag}(\Sigma^{-1}S\Sigma^{-1})=\hbox{diag}(\Sigma^{-1}\hat{S}\Sigma^{-1})

或寫成

\hbox{diag}(\Sigma^{-1}(\Sigma-\hat{S})\Sigma^{-1})=0

上面定義了有偏樣本共變異數矩陣 \hat{S}=[\hat{s}_{ij}]=\frac{n-1}{n}S\hbox{diag}(Y) 表示 Y 的主對角元構成的對角矩陣。計算 \log\mathcal{L}f_{ij} 的偏導數,1\le i,j\le p,如下 (推導見[3]):

\displaystyle  \begin{aligned}  \frac{\partial \log\mathcal{L}}{\partial f_{ij}}&=-\frac{n}{2}\frac{\partial \log|\Sigma|}{\partial f_{ij}}-\frac{n-1}{2}\frac{\partial \,\hbox{trace}(S\Sigma^{-1})}{\partial f_{ij}}\\  &=-n\sum_{l=1}^p\rho_{il}f_{lj}+n\sum_{l=1}^p\sum_{m=1}^p\sum_{r=1}^p\rho_{il}\hat{s}_{lm}\rho_{mr}f_{rj}.  \end{aligned}

設上式等於零,用矩陣形式表示為

\displaystyle  \Sigma^{-1}F=\Sigma^{-1}\hat{S}\Sigma^{-1}F

左乘 \Sigma,可得

\displaystyle  F=\hat{S}\Sigma^{-1}F

接著我們將估計值必須滿足的兩個必要條件 \hbox{diag}(\Sigma^{-1}(\Sigma-\hat{S})\Sigma^{-1})=0F=\hat{S}\Sigma^{-1}F 中的 \Sigma^{-1} 消去。令 \Gamma=F^T\Psi^{-1}F。寫出

\Sigma \Psi^{-1}F=(FF^T+\Psi)\Psi^{-1}F=F\Gamma+F=F(\Gamma+I)

上式左乘 \Sigma^{-1},可得 \Psi^{-1}F=\Sigma^{-1}F(\Gamma+I)。合併以上結果,

F(\Gamma+I)=\hat{S}\Sigma^{-1}F(\Gamma+I)=\hat{S}\Psi^{-1}F

或寫成

(\hat{S}-\Psi)\Psi^{-1}F=F\Gamma

此即 \frac{\partial L}{\partial f_{ij}}=0 給出的估計值所要滿足的必要條件。當上式成立時,我們可以證明[4]

\displaystyle  \Sigma^{-1}(\Sigma-\hat{S})\Sigma^{-1}=\Psi^{-1}(\Sigma-\hat{S})\Psi^{-1}

因此,\hbox{diag}(\Psi^{-1}(\Sigma-\hat{S})\Psi^{-1})=\hbox{diag}(\Psi^{-1}(FF^T+\Psi-\hat{S})\Psi^{-1})=0。但 \Psi^{-1} 是對角矩陣,此式等價於

\displaystyle  \hbox{diag}(FF^T+\Psi)=\hbox{diag}(\hat{S})

此即 \frac{\partial L}{\partial \psi_i}=0 給出的估計值要滿足的必要條件。另外,我們還加入條件2:\Gamma=F^T\Psi^{-1}F 是對角矩陣。由於不存在代數解,你可以採用迭代法,如梯度下降法或牛頓法找出最大似然估計 \hat{F}\hat{\Psi}

旋轉

因素分析有一個特別的應用:解釋隱藏因素的意義,並將可量測變數分群 (clustering)。為便於度量,假設每個隨機變數 x_i 已被標準化,\hbox{E}[x_i]=0\hbox{var}[x_i]=1i=1,\ldots,p。換句話說,共變異數矩陣 \Sigma 被相關係數矩陣 R 取代。前面提到因素分析沒有唯一的因素負荷。令 A=FQ,其中 F 是從數據樣本獲得的因素負荷估計,Q 是一個旋轉矩陣,即 Q 是一個正交矩陣且 \det Q=1。這時候,新因素變成 \mathbf{z}'=\mathbf{z}Q^T,我們仍然保留因素之間的無關性[5]。請注意,旋轉並不會改變因素分析解的結構,它改變的是解的描述方式。理論上,你可以自由地選擇旋轉矩陣 Q 直到對因素獲得滿意的解釋為止。在一般的應用場合,我們希望找到一個旋轉矩陣 Q 使得因素負荷 A 最大化某個目標函數致使變數間存在簡單的結構。量化簡單變數結構的方式很多,下面介紹兩種方法:varimax 旋轉與 quartimax 旋轉。

varimax 旋轉

因為 R=AA^T+\Psi,因素負荷 A(i,j)a_{ij}1\le i\le p1\le j\le k,為變數 x_i 與因素 z_j 的相關係數,滿足 -1\le a_{ij}\le 1 (見“相關係數”)。所謂變數間有簡單的結構是指:第一,每個變數至多僅與一個因素有高相關性;第二,所有的 a_{ij}^2 接近 10。為了實現這個想法,varimax 聚焦於 A 的行 (column)。定義 a^2_{ij} 的第 j 個行變異數為

\displaystyle  v_j=\frac{1}{p}\sum_{i=1}^p(a_{ij}^2)^2-\left(\frac{1}{p}\sum_{i=1}^pa_{ij}^2\right)^2

varimax 程序要找到一個旋轉矩陣 Q 使得所有的行變異數總和有最大值:

\displaystyle  v=\sum_{j=1}^kv_j=\frac{1}{p}\sum_{j=1}^k\left(\sum_{i=1}^pa_{ij}^4-\frac{1}{p}\left(\sum_{i=1}^pa_{ij}^2\right)^2\right)

對於某個 z_j,當 a_{ij}^2 趨於 1,最大化 v 將迫使其餘的 a_{il}^2 趨於零,l\neq i,換句話說,變數 x_i 至多僅與一個因素 z_j 有相關性。據此,與同一因素 z_j 有高相關性的所有變數 x_l 可以視為一組,這樣就達到變數分群的目的。

quartimax 旋轉

不同於 varimax 專注於因素負荷 A 的行,quartimax 旋轉考慮 A 的列 (row)。因為 AA^T=FF^T,變數 i 的共同性 \sum_{j=1}^ka_{ij}^2=(AA^T)_{ii}──透過共同的因素與其他變數共享的變異成分──不會因旋轉 Q 而改變,也就是說 \sum_{j=1}^ka_{ij}^2 是一個常數。再者,所有變數的共同性平方和 \sum_{i=1}^p(\sum_{j=1}^ka_{ij}^2)^2 也是一個常數,展開如下:

\displaystyle  \sum_{i=1}^p\left(\sum_{j=1}^ka_{ij}^2\right)^2=\sum_{i=1}^p\sum_{j=1}^ka_{ij}^4+\sum_{i=1}^p\left(  \sum_{j=1}^k\sum_{l\neq j}a_{ij}^2a_{il}^2\right)

如果變數有簡單的結構,a^2_{ij}a^2_{il}j\neq l,不可同時趨於 1,故上式第二項應該越小越好,換言之,第一項要越大越好。quartimax 程序即在找到一個旋轉矩陣 Q 使下式有最大值:

\displaystyle  q=\sum_{i=1}^p\sum_{j=1}^ka_{ij}^4

因素分數

如同主成分分析,因素分析也可用於數據資料的降維。假設你已經得知因素負荷 F 與特殊共變異數矩陣 \Psi。如何求得對應數據點 \mathbf{x}_i\in\mathbb{R}^p 的因素分數 (factor score) \mathbf{z}_i\in\mathbb{R}^ki=1,\ldots,n?底下介紹兩個方法:確定法與機率法。

確定法

第一個方法把因素分數 \mathbf{z}_i 看成待決定的模型參數,根據 \mathbf{x}_i=F\mathbf{z}_i+\boldsymbol{\epsilon},假設 \mathbf{x}_i\sim\mathcal{N}(F\mathbf{z}_i,\Psi)。似然函數的對數為

\displaystyle  \log\mathcal{L}=-\frac{p}{2}\log(2\pi)-\frac{1}{2}\log\vert \Psi \vert-\frac{1}{2}(\mathbf{x}_i-F\mathbf{z}_i)^T\Psi^{-1}(\mathbf{x}_i-F\mathbf{z}_i)

求導可得

\displaystyle  \frac{\partial \log\mathcal{L}}{\partial \mathbf{z}_i}=-F^T\Psi^{-1}F\mathbf{z}_i+F^T\Psi^{-1}\mathbf{x}_i

設上式等於零,解出最大似然估計

\displaystyle  \hat{\mathbf{z}}_i=\left(F^T\Psi^{-1}F\right)^{-1}F^T\Psi^{-1}\mathbf{x}_i=\Gamma^{-1}F^T\Psi^{-1}\mathbf{x}_i,~~~i=1,\ldots,n

或以矩陣表示為

\hat{Z}=\begin{bmatrix}  \hat{\mathbf{z}}_1^T\\  \vdots\\  \hat{\mathbf{z}}_n^T  \end{bmatrix}=\begin{bmatrix}  \hat{\mathbf{x}}_1^T\\  \vdots\\  \hat{\mathbf{x}}_n^T  \end{bmatrix}\Psi^{-1}F\left(F^T\Psi^{-1}F\right)^{-1}=X\Psi^{-1}F\Gamma^{-1}

機率法

另一個方法將因素分數 \mathbf{z}_i 當作隨機向量,\mathbf{z}_i\sim\mathcal{N}(\mathbf{0},I_k)。因此,聯合機率分布為

\begin{bmatrix}  \mathbf{z}_i\\  \mathbf{x}_i  \end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}  \mathbf{0}\\  \mathbf{0}  \end{bmatrix},\begin{bmatrix}  I_k&F^T\\  F&\Sigma  \end{bmatrix}\right)

由此知 (見“多變量常態分布”,(8))

p(\mathbf{z}_i|\mathbf{x}_i)\sim\mathcal{N}\left(F^T\Sigma^{-1}\mathbf{x}_i,I_k-F^T\Sigma^{-1}F\right)

條件密度函數的期望值是一個自然的估計,使用 \Psi^{-1}F=\Sigma^{-1}F(\Gamma+I)

\hat{\mathbf{z}}_i=F^T\Sigma^{-1}\mathbf{x}_i=(\Gamma+I)^{-1}F^T\Psi^{-1}\mathbf{x}_i,~~~i=1,\ldots,n

或以矩陣表示為

\hat{Z}=X\Sigma^{-1}F=X\Psi^{-1}F(\Gamma+I)^{-1}

因素分析與主成分分析比較

因素分析與主成分分析的相同性質:

  1. 兩者都被歸於探索資料分析分法 (exploratory data analysis)。
  2. 兩者皆可應用於降維。如果特殊變異數 \psi_i 很小,兩者給出相近的因素分數 (標準化後);如果特殊變異數為零,因素分析與主成分分析有相同的結果。

因素分析與主成分分析的相異性質:

  1. 因素分析假設一個資料產生模型,主成分分析沒有任何假設。
  2. 因素分析強調因素至可量測變數的變換,主成分分析強調可量測變數至主成分的變換。
  3. 因素分析的因素分數計算複雜,主成分分析的主成分分數計算相對簡單。
  4. 因素分析 (配合適當的旋轉) 的負載矩陣可用於解釋變數間的關係,主成分分析則不具備此功能。

註解
[1] 統計學的慣例是以希臘字母表示機率模型的參數,許多文獻以 \Lambda 表示因素負荷。考量線性代數多以 \Lambda 表示特徵值矩陣,故本文改用英文字母 F
[2] 令 \Sigma^{-1}=[\rho_{ij}]。使用鏈式法則,Jacobi 公式 d|A|=\hbox{trace}((\hbox{adj}A)(dA)) (見“跡數的性質與應用”),以及逆矩陣公式 A^{-1}=\frac{1}{|A|}\hbox{adj}A (見“伴隨矩陣”),可得

\displaystyle  \begin{aligned}  \frac{\partial \log|\Sigma|}{\partial \psi_i}&=\frac{1}{|\Sigma|}\frac{\partial |\Sigma|}{\partial\psi_i}=\frac{1}{|\Sigma|}\hbox{trace}\left((\hbox{adj}\Sigma)\frac{\partial\Sigma}{\partial \psi_i}\right)\\  &=\hbox{trace}\left(\frac{1}{|\Sigma|}(\hbox{adj}\Sigma)\frac{\partial \Psi}{\partial \psi_i}\right)=\hbox{trace}\left(\Sigma^{-1}E_{ii}\right)\\  &=\rho_{ii},  \end{aligned}

其中 \frac{\partial\Psi}{\partial\psi_i}=E_{ii}(i,i) 元等於 1,其餘元皆為零。因為 \Sigma \Sigma^{-1}=I,對於任意參數 \theta,下列導數成立 (見“矩陣導數”,MS-6):

\displaystyle  \frac{\partial \Sigma^{-1}}{\partial \theta}=-\Sigma^{-1}\frac{\partial \Sigma}{\partial \theta}\Sigma^{-1}

使用跡數循環不變性,

\displaystyle  \begin{aligned}  \frac{\partial\,\hbox{trace}(S\Sigma^{-1})}{\partial \psi_i}&=\hbox{trace}\left(S\frac{\partial\Sigma^{-1}}{\partial \psi_i}\right)=-\hbox{trace}\left(S\Sigma^{-1}\frac{\partial \Sigma}{\partial \psi_i}\Sigma^{-1}\right)\\  &=-\hbox{trace}\left(S\Sigma^{-1}E_{ii}\Sigma^{-1}\right)=-\hbox{trace}\left(\Sigma^{-1}S\Sigma^{-1}E_{ii}\right)\\  &=-(\Sigma^{-1}S\Sigma^{-1})_{ii}=-\sum_{l=1}^p\sum_{j=1}^p\rho_{il}s_{lj}\rho_{ji}.  \end{aligned}

[3] 類似 [2] 的計算方式,

\displaystyle  \begin{aligned}  \frac{\partial \log|\Sigma|}{\partial f_{ij}}&=\frac{1}{|\Sigma|}\frac{\partial |\Sigma|}{\partial f_{ij}}=\frac{1}{|\Sigma|}\hbox{trace}\left((\hbox{adj}\Sigma)\frac{\partial \Sigma}{\partial f_{ij}}\right)\\  &=\hbox{trace}\left(\Sigma^{-1}\frac{\partial (FF^T)}{\partial f_{ij}}\right)=\hbox{trace}\left(\Sigma^{-1}\left(\frac{\partial F}{\partial f_{ij}}F^T+F\frac{\partial F^T}{\partial f_{ij}}\right)\right)\\  &=\hbox{trace}\left(\Sigma^{-1}E_{ij}F^T\right)+\hbox{trace}\left(\Sigma^{-1}FE_{ij}^T\right)\\  &=\hbox{trace}\left(FE^T_{ij}\Sigma^{-1}\right)+\hbox{trace}\left(\Sigma^{-1}FE_{ij}^T\right)\\  &=2\,\hbox{trace}\left(\Sigma^{-1}FE_{ij}^T\right)=2\left(\Sigma^{-1}F\right)_{ij}=2\sum_{l=1}^p\rho_{il}f_{lj},  \end{aligned}

其中 \frac{\partial F}{\partial f_{ij}}=E_{ij}(i,j) 元等於 1,其餘元為零。使用上式結果,

\displaystyle  \begin{aligned}  \frac{\partial\,\hbox{trace}(S\Sigma^{-1})}{\partial f_{ij}}&=\hbox{trace}\left(S\frac{\partial\Sigma^{-1}}{\partial f_{ij}}\right)=-\hbox{trace}\left(S\Sigma^{-1}\frac{\partial \Sigma}{\partial f_{ij}}\Sigma^{-1}\right)\\  &=-2\,\hbox{trace}\left(S\Sigma^{-1}FE_{ij}^T\Sigma^{-1}\right)=-2\,\hbox{trace}\left(\Sigma^{-1}S\Sigma^{-1}FE_{ij}^T\right)\\  &=-2(\Sigma^{-1}S\Sigma^{-1}F)_{ij}=-2\sum_{l=1}^p\sum_{m=1}^p\sum_{r=1}^p\rho_{il}s_{lm}\rho_{mr}f_{rj}.  \end{aligned}

[4] 式 \Psi^{-1}(\Sigma-\hat{S})\Psi^{-1} 左右同時乘 \Sigma=FF^T+\Psi

\begin{aligned}  \Sigma\Psi^{-1}(\Sigma-\hat{S})\Psi^{-1}\Sigma&=(FF^T+\Psi)\Psi^{-1}(FF^T+\Psi-\hat{S})\Psi^{-1}(FF^T+\Psi)\\  &=FF^T\Psi^{-1}(FF^T+\Psi-\hat{S})\Psi^{-1}(FF^T+\Psi)+(FF^T+\Psi-\hat{S})\Psi^{-1}(FF^T+\Psi)\\  &=(FF^T+\Psi-\hat{S})\Psi^{-1}(FF^T+\Psi)\\  &=(FF^T+\Psi-\hat{S})\Psi^{-1}FF^T+FF^T+\Psi-\hat{S}\\  &=FF^T+\Psi-\hat{S}\\  &=\Sigma-\hat{S},  \end{aligned}

因為

\displaystyle  \begin{aligned}  FF^T\Psi^{-1}(FF^T+\Psi-\hat{S})&=F\Gamma F^T+FF^T-FF^T\Psi^{-1}\hat{S}\\  &=F\left((\Gamma+I)F^T-F^T\Psi^{-1}\hat{S}\right)\\  &=F\left(F(\Gamma+I)-\hat{S}\Psi^{-1}F\right)^T\\  &=0,  \end{aligned}

上式取轉置,(FF^T+\Psi-\hat{S})\Psi^{-1}FF^T=0

[5] 除了直交轉軸 (orthogonal rotation),另有一種方式稱為斜交轉軸 (oblique rotation),但我們將失去因素之間的無關性。

廣告
本篇發表於 機器學習 並標籤為 , , , 。將永久鏈結加入書籤。

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s