R語言因子分析
自己整理編寫的R語言常用資料分析模型的模板，原檔案為Rmd格式，直接複製貼上過來，作為個人學習筆記儲存和分享。部分參考薛毅的《統計建模與R軟體》和《R語言實戰》
因子模型： X=μ + A*F* + ε
其中F=[(f1,f2,…,fm)]^T為公共因子向量，[ε=(ε1,ε2,…,εp)]^T為特殊因子向量，A=[(aij)]^(p×m)為因子載荷矩陣。
I.引數估計
為了建立因子模型，需要要得到因子載荷矩陣A=[(aij)]^(p×m)和特殊方差矩陣D=diag(σ1^2,σ2^2,…,σp^2)這兩個引數的估計。
常用的引數估計方法有如下三種：主成分法、主因子法和極大似然法。
接下來會分別介紹以上三種方法具體方法，和綜合三種方法的一個簡便寫法。
例. 12項智力指標的因子分析
研究者收集了40名學生的12項智力指標，分別為常識（x1）、類同（x2）、計算（x3）、詞彙（x4）、理解（x5）、數字廣度（x6）、常填圖（x7）、圖片排列（x8）、積木（x9）、拼圖（x10）、譯碼（x11）和迷津（x12）。將原始資料經過標準化處理後，計算其相關係數矩陣，結果列在下表中。取m=2，試進行因子分析
#輸入相關矩陣的數值
x <- c(
1.000,
0.6904 ,1.000,
0.4115 ,0.4511, 1.000,
0.4580, 0.7068, 0.4018 1.主成分法
（需要設定的引數是R，因子個數m，後面會講到m如何選取）
下面給出主成分法的R程式（factor.analy1.R）
factor.analy1<-function(S, m){
p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
eig<-eigen(S)
for (i in 1:m)
A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
h<-diag(A%*%t(A))
rowname<-c("SS loadings", "Proportion Var", "Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Principal Component Method")
list(method=method, loadings=A,
var=cbind(common=h, spcific=diag_S-h), B=B)
}函式輸入值S是樣本方差陣或相關矩陣，m是主因子的個數，函式的輸出值是列表形式，其內容有估計引數的辦法（主成分法），因子載荷（loadings），共性方差和特殊方差，以及因子F對變數X的貢獻、貢獻率和累積貢獻率。
#呼叫因子分析主成分法的函式
source("factor.analy1.R")
#顯示結果.估計引數的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對變數X的貢獻、貢獻率和累積貢獻率
fa1<-factor.analy1(R, m=2); fa1
#協方差陣S的近似公式，誤差平方和Q(m) (近似公式為E=S-A*A^T-D)
E1 <- R-fa1$loadings %*% t(fa1$loadings)-diag(fa1$var[,2])
sum(E1^2)
因子個數m的選取
#求特徵值，對其求和
eigen(cor(R))
sum(eigen(cor(R))$values)
#選取滿足 m個λ累加/所有λ累加 >= P0 的最小m，P0一般取[0.7,1)
(5.561644e+00 + 1.676901e+00 + 1.434965e+00) / sum(eigen(cor(R))$values)
#可以取m=3
#下面檢驗是否此時Q(m)最小
fa11 <- factor.analy1(R, m=3); fa11
#協方差陣S的近似公式，誤差平方和Q(m) (近似公式為E=S-A*A^T-D)
E11 <- R-fa11$loadings %*% t(fa11$loadings)-diag(fa11$var[,2])
sum(E11^2)
結果看到，sum(E1^2)=1.060286 > sum(E11^2)=0.9550174。說明公因子個數m選擇適當時，近似公式S的誤差平方和Q(m)更優
2.主因子法
（需要設定的引數是R，因子個數m，特殊方差的估計值d，m值選取參考主成分法，d選取方法後面會講到）
按照主因子法的思想編寫相應的R程式：（factor.analy2.R）
factor.analy2<-function(R, m, d){
p<-nrow(R); diag_R<-diag(R); sum_rank<-sum(diag_R)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
kmax=20; k<-1; h <- diag_R-d
repeat{
diag(R)<- h; h1<-h; eig<-eigen(R)
for (i in 1:m)
A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
h<-diag(A %*% t(A))
if ((sqrt(sum((h-h1)^2))<1e-4)|k==kmax) break
k<-k+1
}
rowname<-c("SS loadings", "Proportion Var", "Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Principal Factor Method")
list(method=method, loadings=A,
var=cbind(common=h, spcific=diag_R-h), B=B, iterative=k)
} 函式輸入值R是樣本方差陣或相關矩陣，m是主因子的個數，d是特殊方差的估計值，函式的輸出值是列表形式，其內容有估計引數的辦法（主因子法），因子載荷（loadings），共性方差和特殊方差，以及因子F對變數X的貢獻、貢獻率和累積貢獻率，以及求解的迭代次數。
相同資料，相關係數矩陣R，取公因子個m=2，特殊方差的估計值為0
#輸入特殊方差var$spcific估計值，可以全部取0，下面會介紹怎麼取合適的特殊方差估計值
d<-c(0,0,0,0,0,0,0,0,0,0,0,0)
#呼叫呼叫因子分析主因子法的函式
source("factor.analy2.R")
#顯示結果.估計引數的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對變數X的貢獻、貢獻率和累積貢獻率，iterative-迭代次數
fa2<-factor.analy2(R, m=3, d); fa2
#近似公式S的誤差平方和Q(m)
E2<- R-fa2$loadings %*% t(fa2$loadings)-diag(fa2$var[,2])
sum(E2^2)
用了13次迭代得到穩定解，再計算Q(m)
sum(E2^2)=0.3141111，優於主成分法
特殊方差估計值σi^2的常用選取方法
## 1.σi^2 = 1/rii,其中rii為R的逆矩陣的第i個對角線元素，此時Q(m)=sum(E21^2)
#R的逆矩陣R^-1
solve(R)
#取其對角線值，再求倒數
1 / diag(solve(R))
#將剛才的結果作為特殊方差估計值，我們來驗證是否Q(m)會更優
d1 <- c(0.4113202,0.2159605,0.5974511,0.3610979,0.3659987,0.4522035,0.4673815,0.7639169,0.4743578,0.6385381,0.6627739,0.5743706)
#呼叫呼叫因子分析主因子法的函式
source("factor.analy2.R")
#顯示結果.估計引數的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對變數X的貢獻、貢獻率和累積貢獻率，iterative-迭代次數
fa21 <- factor.analy2(R, m=3, d1); fa21
#近似公式S的誤差平方和Q(m)
E21 <- R-fa21$loadings %*% t(fa21$loadings)-diag(fa21$var[,2])
sum(E21^2)
## 2.σi^2 = 1-hi^2,其中hi^2=max(j/i) |rij|
## 3.σi^2 = 1-hi^2,其中hi^2=1,此時σi^2全取0，此時Q(m)=sum(E2^2)
#### 這裡R為相關矩陣，對角線元素全為1，其餘元素都為0-1間的小數，所以方法2.和3.在這裡是一樣的sum(E21^2) = 0.3106186 < sum(E2^2)=0.3141111，證明特殊方差估計值的選取方法，1.要優於2.、3.
3.極大似然法
（需要設定的引數是R，因子個數m，特殊方差的估計值d，m值選取參考主成分法，d值選取參考主因子法）
按照極大似然法的思想編寫相應的R程式：（factor.analy3.R）
factor.analy3<-function(S, m, d){
p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
kmax=20; k<-1
repeat{
d1<-d; d2<-1/sqrt(d); eig<-eigen(S * (d2 %o% d2))
for (i in 1:m)
A[,i]<-sqrt(eig$values[i]-1)*eig$vectors[,i]
A<-diag(sqrt(d)) %*% A
d<-diag(S-A%*%t(A))
if ((sqrt(sum((d-d1)^2))<1e-4)|k==kmax) break
k<-k+1
}
rowname<-c("SS loadings","Proportion Var","Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Maximum Likelihood Method")
list(method=method, loadings=A,
var=cbind(common=diag_S-d, spcific=d),B=B,iterative=k)
} 函式輸入值R是樣本方差陣或相關矩陣，m是主因子的個數，d是特殊方差的估計值，函式的輸出值是列表形式，其內容有估計引數的辦法（主因子法），因子載荷（loadings），共性方差和特殊方差，以及因子F對變數X的貢獻、貢獻率和累積貢獻率，以及求解的迭代次數。
相同資料，相關係數矩陣R，取公因子個m=2，特殊方差的估計值為:
#輸入特殊方差var$spcific估計值（用上例中方法1.的結果d1）
d1 <- c(0.4113202,0.2159605,0.5974511,0.3610979,0.3659987,0.4522035,0.4673815,0.7639169,0.4743578,0.6385381,0.6627739,0.5743706)
#呼叫呼叫因子分析極大似然法的函式
source("factor.analy3.R")
#顯示結果.估計引數的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對變數X的貢獻、貢獻率和累積貢獻率，iterative-迭代次數
fa3 <- factor.analy3(R, m=3, d1); fa3
#近似公式S的誤差平方和Q(m)
E3 <- R-fa3$loadings %*% t(fa3$loadings)-diag(fa3$var[,2])
sum(E3^2)
sum(E3^2) = 0.3412492
4.綜合以上三種方法
（method=“xxx”）
將上述3種方法結合在一起，並考慮主成分估計中介紹的因子個數m的選取方法，和在主因子法中介紹的特殊方差初始估計方法，編寫相應的R程式
factor.analy.R
用一條函式，通過改變引數method=“xxx” ， 可以更方便對比三種方法的結果
factor.analy<-function(S, m=0,
d=1/diag(solve(S)), method="likelihood"){
if (m==0){
p<-nrow(S); eig<-eigen(S)
sum_eig<-sum(diag(S))
for (i in 1:p){
if (sum(eig$values[1:i])/sum_eig>0.70){
m<-i; break
}
}
}
source("factor.analy1.R")
source("factor.analy2.R")
source("factor.analy3.R")
switch(method,
princomp=factor.analy1(S, m), #method=“princomp”時輸入S,m=i兩個引數
factor=factor.analy2(S, m, d), #method=“factor”時輸入S,m=i,d=c(x,..,x)三個引數
likelihood=factor.analy3(S, m, d) #method=“likehood”時輸入S,m=i,d=c(x,..,x)三個引數
)
}函式輸入樣本方差矩陣S或樣本相關矩陣R。因子個數m（預設值由貢獻率計算出m值）。特殊方差的初始估計d（預設值為^σi方 = 1/rii）
計算因子載荷的方法，method=princomp採用主成分法，method=factor採用主因子法，method=likelihood（預設值）採用極大似然法
函式輸出就是採用前面介紹的三種方法的輸出格式。
#使用factor.analy.R的例項：
source("factor.analy.R")
fa4 <- factor.analy(S=R,m=3,method = "princomp") ; fa4
#近似公式S的誤差平方和Q(m)
E4 <- R-fa4$loadings %*% t(fa4$loadings)-diag(fa4$var[,2])
sum(E4^2) #可以看到，這裡E4算出的Q(m)與E11算出的Q(m)是相同的
II.方差最大的正交旋轉
某醫院為了合理評價該院各月的醫療工作質量，蒐集了3年有關門診人次、出院人數、病床利用率、病床週轉次數、平均住院天數、治癒好轉率、病死率、診斷符合率、搶救成功率等9個指標資料，試採用因子分析法，探討其綜合評價指標體系。
使方差最大的因子載荷矩陣
先用三種方法之一計算的因子載荷估計矩陣，再用varimax()函式得到方差最大的因子載荷矩陣
#匯入原始資料
hospital <- read.csv("hospital.csv",header=T)
#生成hospital表格的相關係數矩陣R
R <- cor(hospital)
for (i in 1:9){
for (j in 1:i){
R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]
}
}
#呼叫因子分析的特殊方差初始估計方法
source("factor.analy.R")
#以princomp方法為例
fa<-factor.analy(R, m=2, method="princomp")
vm1<-varimax(fa$loadings, normalize = F); vm1因子分析的計算函式
事實上，在R軟體中，提供了作因子分析計算的函式–factanal()函式，它可以從樣本資料、樣本的方差矩陣和相關矩陣出發對資料作因子分析，並可直接給出方差最大的載荷因子矩陣。
#顯示factanal()函式的幫助頁面，引數設定問題
?factanal()
#取公因子個數m=2，選用II中例子裡的相關係數矩陣R，利用factanal函式得到fa結果
fa <- factanal(factors = 4,covmat = R)
#或者不用相關係數矩陣R，直接用csv格式檔案：fa <- factanal(X=~.,factors=2,data=hospital)
#顯示結果
fa在上述資訊中，call表示呼叫函式的方法，uniquenesses是特殊方差，loadings是因子載荷矩陣，其中Factor1，Factor2是因子，X1，X2，…,X9是對應的變數，SS loadings是公共因子對變數X的總方差貢獻，Proportion Var是方差貢獻率，Cumulative Var是累積方差貢獻率。
IV.因子得分
迴歸法和加權最小二乘法
##匯入原始資料
hospital <- read.csv("hospital.csv",header=T)
#相關矩陣特徵值
eigen(cor(R))$values
sum(eigen(cor(R))$values[1:3])/sum(eigen(cor(R))$values)
#前3個因子的累積貢獻率達到0.8134434,接下來選取因子個數為3
#不同方法計算因子得分
fa_1<-factanal(~., factors=3, data=hospital, scores="Bartlett") #加權最小二乘法
fa_2<-factanal(~., factors=3, data=hospital, scores="regression") #迴歸法
fa_1;fa_2
#畫出各組在第a、第b公共因子下的散點圖
plot(fa$scores[, 1:2], type="n"); text(fa$scores[,1], fa$scores[,2]) #第一、第二公共因子下的散點圖
plot(fa$scores[, c(1,3)], type="n"); text(fa$scores[,1], fa$scores[,3]) #第一、第三公共因子下的散點圖上面是採用迴歸法，也可以使用加權最小二乘法來畫圖。由散點圖，可以直觀選出偏向哪個公共因子的組別。
根據選項factors=4的設定，3個潛在因子被保留，前3個因子的累積貢獻率達到81.3%，上式為全部變數在3個潛在因子F1-F3上的因子載荷矩陣
例如：x1由4個因子表達的式子為：
x1=0.447*F1 + 0.519*F2 + -0.101*F4
從矩陣上看，因子1在多數原始指標上均有較大的載荷，因子2在x1（門診人次）、x2（出院人數）、x3（病床利用率）、x4（病床週轉次數）上有較大的載荷，因子3在x2（出院人數）、x5（平均住院天數）、x6（治癒好轉率）上有較大的載荷。除因子1可以認定為綜合因子外，其他3個因子意義不明顯。
因子旋轉
fa_1<-factanal(~., factors=3, data=hospital, scores="Bartlett") #加權最小二乘法
vm <- varimax(fa_1$loadings,normalize = F) ; vm經過因子旋轉處理，3個潛在因子在9個原始指標上的因子載荷矩陣如上表所示。
對該因子載荷進行分析，可看出：因子F1在x1（門診人次）、x2（出院人數）、x5（平均住院天數）、x8（診斷符合率）、x9（搶救成功率）上因子載荷較大；F2在x3（病床利用率）、x4（病床週轉次數）上的因子載荷較大；F3在x6（治癒好轉率）、x7（病死率）上的因子載荷較大
我們可以推出：因子F1反映了該醫院醫療工作質量各方面的情況，為綜合因子；F2反映了病床利用情況；F3反映了醫療水平的高低
將旋轉後的因子載荷與主成分分析的因子載荷矩陣比較可知：因子旋轉後，除F1的因子載荷仍分佈多數指標上外，其他2個因子的載荷明顯地集中到少數指標上，說明旋轉對因子載荷起到明顯的分離作用，使得各因子解釋的變數更加清晰。