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.000,
   0.5535, 0.6620, 0.4122, 0.7119, 1.000,
   0.3923, 0.6317, 0.4520, 0.4583, 0.5299, 1.000,
   0.1415, 0.3009, 0.2025, 0.2665, 0.2480, 0.1590, 1.000,
   0.0077, 0.0344, 0.1855, 0.1065, 0.0003, 0.1100, 0.3595, 1.000,
   0.2385, 0.3523, 0.3646, 0.3644, 0.3388, 0.3982, 0.5004, 0.3314, 1.000,
   0.0333, 0.1726, 0.1311, 0.1757, 0.1998, 0.0342, 0.5758, 0.1420, 0.2808, 1.000,
   0.0898, 0.3878, 0.2041, 0.3191, 0.3186, 0.2914, 0.2537, 0.2025, 0.3971, 0.1468, 1.000,
   0.2215, 0.2427, 0.4124, 0.2169, 0.1459, 0.0985, 0.4222, 0.2156, 0.5016, 0.2286, 0.0776, 1.000)
names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12")
R<-matrix(0, nrow=12, ncol=12, dimnames=list(names, names))
#生成相關系數矩陣R
for (i in 1:12){
   for (j in 1:i){
      R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]
   }
}

1.主成分法
（需要設置的參數是R，因子個數m，后面會講到m如何選?。?br /> 下面給出主成分法的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個因子的載荷明顯地集中到少數指標上，說明旋轉對因子載荷起到明顯的分離作用，使得各因子解釋的變量更加清晰。