PCA和EFA
PCA中，透過對多元資料進行降維，將大量相關變數轉化為一組很少的不相關變數，這些無關變數稱為主成分（Principal Component，PC），是所有觀測變數的線性組合。形成線性組合的權重都是透過最大化各主成分所解釋的方差來獲得，同時還要保證各主成分間不相關（彼此正交）。例如，使用PCA將30種環境變數（之間可能存在大量相關或冗餘）轉化為幾個無關的主成分，並儘可能保留原始資料集的資訊。
與PCA相比，EFA是一系列用來發現一組變數的潛在結構的方法，它透過尋找一組更小的、潛在的或隱藏的結構來解釋已觀測到的、顯式的變數間的關係。EFA的目標是透過發掘隱藏在資料下的一組較少的、更為基本的無法觀測的變數，來解釋一組可觀測變數的相關性，這些虛擬的、無法觀測的變數稱為因子（Factor，F），因子被當作觀測變數的結構基礎或原因，而不是代表它們的線性組合。例如，使用EFA探索30種環境變數的相互關係，可用三個潛在因子（地理位置、氣候環境、生物因素）來表示。
每個因子被認為可解釋多個觀測變數間共有的方差，因此它們常被稱為公共因子。代表觀測變數方差的誤差（error，e）無法用因子來解釋。因子和誤差均無法直接測量，透過變數間的相互關係推導得到。
圖注：主成分分析和因子分析模式圖，展示了可觀測變數（X1到X5）、主成分（PC1、PC2）、因子（F1、F2）和誤差（e1到e5）。
R包psych的EFA
接下來展示R語言執行EFA的方法，以psych包的方法為例。
資料集
模擬生成了包含200個物件（行）6個變數（列）的資料集（X），接下來使用EFA探索變數間的關係。
注：真實的用於因子分析的資料中，不可能只對這麼少的變數進行因子分析，此處模擬資料集僅為操作方便。
#模擬資料
library(mvtnorm)
N <- 200
P <- 6
Q <- 2
Lambda <- matrix(c(0.7, -0.4, 0.8, 0, -0.2, 0.9, -0.3, 0.4, 0.3, 0.7, -0.8, 0.1), nrow = P, ncol = Q, byrow = TRUE)
set.seed(123)
Kf <- diag(Q)
mu <- c(5, 15)
FF <- rmvnorm(N, mean = mu, sigma = Kf)
E <- rmvnorm(N, mean = rep(0, P), sigma = diag(P))
X <- FF %*% t(Lambda) + E
也可首先計算6個變數間的相關性。
#計算變數間的相關矩陣，以 pearson 相關係數為例
corMat <- cor(X, method = 'pearson')
head(corMat)
#相關圖
library(corrplot)
corrplot(corMat, method = 'number', number.cex = 0.8, diag = FALSE, tl.cex = 0.8)
corrplot(corMat, add = TRUE, type = 'upper', method = 'pie', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
判斷需提取的公共因子數量
EFA要求指定提取的公共因子數量，可在執行EFA前作個評估。
library(psych)
#確定最佳因子數量，詳情 ?fa.parallel
#輸入變數間的相關矩陣，並手動輸入原始變數集的物件數量（n.obs）
fa.parallel(corMat, n.obs = N, fa = 'both', n.iter = 100)
#或者直接使用原始變數集
fa.parallel(X, fa = 'both', n.iter = 100)
顯示提取兩個因子為宜。
碎石檢驗的前兩個特徵值（藍線三角形）都在拐角處之上，並且大於基於100次隨機模擬資料矩陣的特徵值均值（紅色虛線）。
並且，根據Kaiser-Harris準則的特徵值數大於0，也推薦選擇兩個因子。（注意不是根據Kaiser-Harris>1的標準，Kaiser-Harris>1是對PCA而言的，如果期望在PCA中評估選擇多少個主成分合適，則可使用該標準）
提取公共因子及因子旋轉
執行EFA。上述評估提示選擇兩個因子為宜，在該步將它作為引數輸入。
#EFA 分析，詳情 ?fa
#輸入變數間的相關矩陣，nfactors 指定提取的因子數量，並手動輸入原始變數集的物件數量（n.obs）
#rotate 設定旋轉的方法，fm 設定因子化方法
fa_varimax <- fa(r = corMat, nfactors = 2, n.obs = N, rotate = 'varimax', fm = 'pa')
#或者也可直接使用原始變數集
fa_varimax <- fa(r = X, nfactors = 2, rotate = 'varimax', fm = 'pa')
fa_varimax
EFA中提取因子的方法很多，如最大似然法（ml）、主軸迭代法（pa）、加權最小二乘法（wls）、廣義加權最小二乘法（gls）、最小殘差法（minres）等。這裡選擇了主軸迭代法。
因子旋轉是一系列將因子載荷陣變得更容易解釋的數學方法，它們儘可能地對因子去噪。旋轉方法有兩種：使選擇的因子保持不相關（正交旋轉），和讓它們變得相關（斜交旋轉）。最流行的正交旋轉是方差極大旋轉（式中引數rotate = 'varimax'），它試圖對載荷陣的列進行去噪，使每個因子只由一組有限的變數來解釋（即載荷陣每列只有少數幾個很大的載荷，其它都是很小的載荷）。
因子模式矩陣（標準化的迴歸係數矩陣，列出了因子預測變數的權重）顯示，前兩個因子解釋了資料集的39%的方差。
對於該模擬的資料集，變數1、2在第一因子上載荷較大，變數3在第二因子上載荷較大。
使用正交旋轉將人為地強制兩個因子不相關。如果允許因子間存在相關，則可使用引數rotate = 'promax'。
#上述展示了正交旋轉，以下是斜交旋轉
fa_promax <- fa(r = corMat, nfactors = 2, n.obs = N, rotate = 'promax', fm = 'pa')
fa_promax
因子模式矩陣顯示，前兩個因子解釋了資料集的39%的方差。
因子關聯矩陣顯示兩個因子之間存在負相關關係。在因子之間存在較大相關關係時，因子載荷陣會出現較明顯的噪聲（相較於正交旋轉而言，因為允許潛在因子相關），雖然斜交方法更為複雜，但更符合真實資料。如果因子間的關聯很低，推薦更換為正交旋轉方法簡化問題。
#斜交旋轉中，因子結構矩陣（載荷陣）使用 F=P*Phi 獲得
#F 為因子載荷陣，P 為因子模式矩陣，Phi 為因子關聯矩陣
fsm <- function(oblique) {
P <- unclass(oblique$loading)
F <- P %*% oblique$Phi
colnames(F) <- c('PA1', 'PA2')
F
}
fsm(fa_promax)
根據因子結構矩陣，變數1、2在第一因子上載荷較大，變數3在第二因子上載荷較大。因為這個模擬資料的兩個因子之間相關程度很小，所以斜交旋轉的載荷陣和正交旋轉相似。
繪製正交或斜交結果的圖形。
#視覺化
par(mfrow = c(2, 2))
factor.plot(fa_varimax, title = '正交旋轉')
fa.diagram(fa_varimax, simple = TRUE, main = '正交旋轉')
factor.plot(fa_promax, title = '斜交旋轉')
fa.diagram(fa_promax, simple = FALSE, main = '斜交旋轉')
可以據此評估哪些變數在哪些因子上的載荷更大，以確定哪些因子主要代表了哪些變數間的相互關係。
如果出現了未能被已知因子所代表的變數（例如，本模擬資料集的變數4），則表明還存在潛在的因子。如果有必要，則可以在計算過程中增添提取因子的數量。
因子得分
儘管EFA中不怎麼關注因子得分，如果期望獲得，需使用原始變數矩陣（而非變數間的相關矩陣）計算。
與精確計算的主成分得分不同，因子得分只是估計得到的。它的估計方法有很多種，fa()預設使用迴歸方法計算（計算時新增引數score='regression'）。
#因子得分，以斜交旋轉結果為例
fa_promax <- fa(r = X, nfactors = 2, rotate = 'promax', fm = 'pa', score = 'regression')
head(fa_promax$scores)
#得分系數（標準化的迴歸權重）
head(fa_promax$weight)
其它常用於EFA的R包
FactoMineR包不僅提供了PCA和EFA方法，還包含潛變數模型，並提供了分別針對數值型變數和分類變數的計算方法。FAiR包使用遺傳演算法估計因子分析模型。它增強了模型引數估計能力，能夠處理不等式的約束條件。GPArotatio包則提供了許多因子旋轉方法。nFactors包提供了用來判斷因子數目的許多複雜方法。
參考資料
