R語(yǔ)言因子分析

因子模型： X=μ + A*F* + ε
其中F=[(f1,f2,…,fm)]^T為公共因子向量，[ε=(ε1,ε2,…,εp)]^T為特殊因子向量，A=[(aij)]^(p×m)為因子載荷矩陣。

I.參數(shù)估計(jì)

為了建立因子模型，需要要得到因子載荷矩陣A=[(aij)]^(p×m)和特殊方差矩陣D=diag(σ1^2,σ2^2,…,σp^2)這兩個(gè)參數(shù)的估計(jì)。
常用的參數(shù)估計(jì)方法有如下三種：主成分法、主因子法和極大似然法。

接下來(lái)會(huì)分別介紹以上三種方法具體方法，和綜合三種方法的一個(gè)簡(jiǎn)便寫法。

例. 12項(xiàng)智力指標(biāo)的因子分析

研究者收集了40名學(xué)生的12項(xiàng)智力指標(biāo)，分別為常識(shí)（x1）、類同（x2）、計(jì)算（x3）、詞匯（x4）、理解（x5）、數(shù)字廣度（x6）、常填圖（x7）、圖片排列（x8）、積木（x9）、拼圖（x10）、譯碼（x11）和迷津（x12）。將原始數(shù)據(jù)經(jīng)過(guò)標(biāo)準(zhǔn)化處理后，計(jì)算其相關(guān)系數(shù)矩陣，結(jié)果列在下表中。取m=2，試進(jìn)行因子分析

#輸入相關(guān)矩陣的數(shù)值
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))
#生成相關(guān)系數(shù)矩陣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.主成分法
（需要設(shè)置的參數(shù)是R，因子個(gè)數(shù)m，后面會(huì)講到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)
}

函數(shù)輸入值S是樣本方差陣或相關(guān)矩陣，m是主因子的個(gè)數(shù)，函數(shù)的輸出值是列表形式，其內(nèi)容有估計(jì)參數(shù)的辦法（主成分法），因子載荷（loadings），共性方差和特殊方差，以及因子F對(duì)變量X的貢獻(xiàn)、貢獻(xiàn)率和累積貢獻(xiàn)率。

#調(diào)用因子分析主成分法的函數(shù)
source("factor.analy1.R")
#顯示結(jié)果.估計(jì)參數(shù)的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對(duì)變量X的貢獻(xiàn)、貢獻(xiàn)率和累積貢獻(xiàn)率
fa1<-factor.analy1(R, m=2); fa1
#協(xié)方差陣S的近似公式，誤差平方和Q(m)     (近似公式為E=S-A*A^T-D)
E1 <- R-fa1$loadings %*% t(fa1$loadings)-diag(fa1$var[,2])
sum(E1^2)
因子個(gè)數(shù)m的選取

#求特征值，對(duì)其求和
eigen(cor(R))
sum(eigen(cor(R))$values)
#選取滿足 m個(gè)λ累加/所有λ累加 >= P0 的最小m，P0一般取[0.7,1)
(5.561644e+00 + 1.676901e+00 + 1.434965e+00) / sum(eigen(cor(R))$values)
#可以取m=3
#下面檢驗(yàn)是否此時(shí)Q(m)最小
fa11 <- factor.analy1(R, m=3); fa11
#協(xié)方差陣S的近似公式，誤差平方和Q(m)     (近似公式為E=S-A*A^T-D)
E11 <- R-fa11$loadings %*% t(fa11$loadings)-diag(fa11$var[,2])
sum(E11^2)

結(jié)果看到，sum(E1^2)=1.060286 > sum(E11^2)=0.9550174。說(shuō)明公因子個(gè)數(shù)m選擇適當(dāng)時(shí)，近似公式S的誤差平方和Q(m)更優(yōu)

2.主因子法

（需要設(shè)置的參數(shù)是R，因子個(gè)數(shù)m，特殊方差的估計(jì)值d，m值選取參考主成分法，d選取方法后面會(huì)講到）

按照主因子法的思想編寫相應(yīng)的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)
}  

函數(shù)輸入值R是樣本方差陣或相關(guān)矩陣，m是主因子的個(gè)數(shù)，d是特殊方差的估計(jì)值，函數(shù)的輸出值是列表形式，其內(nèi)容有估計(jì)參數(shù)的辦法（主因子法），因子載荷（loadings），共性方差和特殊方差，以及因子F對(duì)變量X的貢獻(xiàn)、貢獻(xiàn)率和累積貢獻(xiàn)率，以及求解的迭代次數(shù)。

相同數(shù)據(jù)，相關(guān)系數(shù)矩陣R，取公因子個(gè)m=2，特殊方差的估計(jì)值為0
#輸入特殊方差var$spcific估計(jì)值，可以全部取0，下面會(huì)介紹怎么取合適的特殊方差估計(jì)值
d<-c(0,0,0,0,0,0,0,0,0,0,0,0)
#調(diào)用調(diào)用因子分析主因子法的函數(shù)
source("factor.analy2.R")
#顯示結(jié)果.估計(jì)參數(shù)的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對(duì)變量X的貢獻(xiàn)、貢獻(xiàn)率和累積貢獻(xiàn)率，iterative-迭代次數(shù)
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次迭代得到穩(wěn)定解，再計(jì)算Q(m)
sum(E2^2)=0.3141111，優(yōu)于主成分法

特殊方差估計(jì)值σi^2的常用選取方法

##  1.σi^2 = 1/rii,其中rii為R的逆矩陣的第i個(gè)對(duì)角線元素，此時(shí)Q(m)=sum(E21^2)
#R的逆矩陣R^-1
solve(R)
#取其對(duì)角線值，再求倒數(shù)
1 / diag(solve(R))  
#將剛才的結(jié)果作為特殊方差估計(jì)值，我們來(lái)驗(yàn)證是否Q(m)會(huì)更優(yōu)
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)
#調(diào)用調(diào)用因子分析主因子法的函數(shù)
source("factor.analy2.R")
#顯示結(jié)果.估計(jì)參數(shù)的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對(duì)變量X的貢獻(xiàn)、貢獻(xiàn)率和累積貢獻(xiàn)率，iterative-迭代次數(shù)
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,此時(shí)σi^2全取0，此時(shí)Q(m)=sum(E2^2)
#### 這里R為相關(guān)矩陣，對(duì)角線元素全為1，其余元素都為0-1間的小數(shù)，所以方法2.和3.在這里是一樣的
sum(E21^2) = 0.3106186 < sum(E2^2)=0.3141111，證明特殊方差估計(jì)值的選取方法，1.要優(yōu)于2.、3.

3.極大似然法

（需要設(shè)置的參數(shù)是R，因子個(gè)數(shù)m，特殊方差的估計(jì)值d，m值選取參考主成分法，d值選取參考主因子法）

按照極大似然法的思想編寫相應(yīng)的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)
}  

函數(shù)輸入值R是樣本方差陣或相關(guān)矩陣，m是主因子的個(gè)數(shù)，d是特殊方差的估計(jì)值，函數(shù)的輸出值是列表形式，其內(nèi)容有估計(jì)參數(shù)的辦法（主因子法），因子載荷（loadings），共性方差和特殊方差，以及因子F對(duì)變量X的貢獻(xiàn)、貢獻(xiàn)率和累積貢獻(xiàn)率，以及求解的迭代次數(shù)。

相同數(shù)據(jù)，相關(guān)系數(shù)矩陣R，取公因子個(gè)m=2，特殊方差的估計(jì)值為:

#輸入特殊方差var$spcific估計(jì)值（用上例中方法1.的結(jié)果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)
#調(diào)用調(diào)用因子分析極大似然法的函數(shù)
source("factor.analy3.R")
#顯示結(jié)果.估計(jì)參數(shù)的方法為主成分法，loadings-因子載荷，var-共性方差和特殊方差，以及B-因子F對(duì)變量X的貢獻(xiàn)、貢獻(xiàn)率和累積貢獻(xiàn)率，iterative-迭代次數(shù)
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種方法結(jié)合在一起，并考慮主成分估計(jì)中介紹的因子個(gè)數(shù)m的選取方法，和在主因子法中介紹的特殊方差初始估計(jì)方法，編寫相應(yīng)的R程序
factor.analy.R

用一條函數(shù)，通過(guò)改變參數(shù)method=“xxx” ， 可以更方便對(duì)比三種方法的結(jié)果

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”時(shí)輸入S,m=i兩個(gè)參數(shù)
             factor=factor.analy2(S, m, d),    #method=“factor”時(shí)輸入S,m=i,d=c(x,..,x)三個(gè)參數(shù)
             likelihood=factor.analy3(S, m, d) #method=“l(fā)ikehood”時(shí)輸入S,m=i,d=c(x,..,x)三個(gè)參數(shù)
          )
}

函數(shù)輸入樣本方差矩陣S或樣本相關(guān)矩陣R。因子個(gè)數(shù)m（缺省值由貢獻(xiàn)率計(jì)算出m值）。特殊方差的初始估計(jì)d（缺省值為^σi方 = 1/rii）
計(jì)算因子載荷的方法，method=princomp采用主成分法，method=factor采用主因子法，method=likelihood（缺省值）采用極大似然法
函數(shù)輸出就是采用前面介紹的三種方法的輸出格式。

#使用factor.analy.R的實(shí)例：
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.方差最大的正交旋轉(zhuǎn)

某醫(yī)院為了合理評(píng)價(jià)該院各月的醫(yī)療工作質(zhì)量，搜集了3年有關(guān)門診人次、出院人數(shù)、病床利用率、病床周轉(zhuǎn)次數(shù)、平均住院天數(shù)、治愈好轉(zhuǎn)率、病死率、診斷符合率、搶救成功率等9個(gè)指標(biāo)數(shù)據(jù)，試采用因子分析法，探討其綜合評(píng)價(jià)指標(biāo)體系。

使方差最大的因子載荷矩陣

先用三種方法之一計(jì)算的因子載荷估計(jì)矩陣，再用varimax()函數(shù)得到方差最大的因子載荷矩陣

#導(dǎo)入原始數(shù)據(jù)
hospital <- read.csv("hospital.csv",header=T)
#生成hospital表格的相關(guān)系數(shù)矩陣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]
   }
}
#調(diào)用因子分析的特殊方差初始估計(jì)方法
source("factor.analy.R")
#以princomp方法為例
fa<-factor.analy(R, m=2, method="princomp")
vm1<-varimax(fa$loadings, normalize = F); vm1

因子分析的計(jì)算函數(shù)

事實(shí)上，在R軟件中，提供了作因子分析計(jì)算的函數(shù)–factanal()函數(shù)，它可以從樣本數(shù)據(jù)、樣本的方差矩陣和相關(guān)矩陣出發(fā)對(duì)數(shù)據(jù)作因子分析，并可直接給出方差最大的載荷因子矩陣。

#顯示factanal()函數(shù)的幫助頁(yè)面，參數(shù)設(shè)置問(wèn)題
?factanal()
#取公因子個(gè)數(shù)m=2，選用II中例子里的相關(guān)系數(shù)矩陣R，利用factanal函數(shù)得到fa結(jié)果
fa <- factanal(factors = 4,covmat = R)
#或者不用相關(guān)系數(shù)矩陣R，直接用csv格式文件：fa <- factanal(X=~.,factors=2,data=hospital)
#顯示結(jié)果
fa

在上述信息中，call表示調(diào)用函數(shù)的方法，uniquenesses是特殊方差，loadings是因子載荷矩陣，其中Factor1，F(xiàn)actor2是因子，X1，X2，…,X9是對(duì)應(yīng)的變量，SS loadings是公共因子對(duì)變量X的總方差貢獻(xiàn)，Proportion Var是方差貢獻(xiàn)率，Cumulative Var是累積方差貢獻(xiàn)率。
IV.因子得分

回歸法和加權(quán)最小二乘法

##導(dǎo)入原始數(shù)據(jù)
hospital <- read.csv("hospital.csv",header=T)
#相關(guān)矩陣特征值
eigen(cor(R))$values
sum(eigen(cor(R))$values[1:3])/sum(eigen(cor(R))$values)
#前3個(gè)因子的累積貢獻(xiàn)率達(dá)到0.8134434,接下來(lái)選取因子個(gè)數(shù)為3
#不同方法計(jì)算因子得分
fa_1<-factanal(~., factors=3, data=hospital, scores="Bartlett")   #加權(quán)最小二乘法
fa_2<-factanal(~., factors=3, data=hospital, scores="regression") #回歸法
fa_1;fa_2
#畫出各組在第a、第b公共因子下的散點(diǎn)圖
plot(fa$scores[, 1:2], type="n"); text(fa$scores[,1], fa$scores[,2])    #第一、第二公共因子下的散點(diǎn)圖
plot(fa$scores[, c(1,3)], type="n"); text(fa$scores[,1], fa$scores[,3]) #第一、第三公共因子下的散點(diǎn)圖

上面是采用回歸法，也可以使用加權(quán)最小二乘法來(lái)畫圖。由散點(diǎn)圖，可以直觀選出偏向哪個(gè)公共因子的組別。

根據(jù)選項(xiàng)factors=4的設(shè)定，3個(gè)潛在因子被保留，前3個(gè)因子的累積貢獻(xiàn)率達(dá)到81.3%，上式為全部變量在3個(gè)潛在因子F1-F3上的因子載荷矩陣

例如：x1由4個(gè)因子表達(dá)的式子為：

x1=0.447*F1 + 0.519*F2 + -0.101*F4

從矩陣上看，因子1在多數(shù)原始指標(biāo)上均有較大的載荷，因子2在x1（門診人次）、x2（出院人數(shù)）、x3（病床利用率）、x4（病床周轉(zhuǎn)次數(shù)）上有較大的載荷，因子3在x2（出院人數(shù)）、x5（平均住院天數(shù)）、x6（治愈好轉(zhuǎn)率）上有較大的載荷。除因子1可以認(rèn)定為綜合因子外，其他3個(gè)因子意義不明顯。

因子旋轉(zhuǎn)
fa_1<-factanal(~., factors=3, data=hospital, scores="Bartlett")   #加權(quán)最小二乘法
vm <- varimax(fa_1$loadings,normalize = F) ; vm

經(jīng)過(guò)因子旋轉(zhuǎn)處理，3個(gè)潛在因子在9個(gè)原始指標(biāo)上的因子載荷矩陣如上表所示。

對(duì)該因子載荷進(jìn)行分析，可看出：因子F1在x1（門診人次）、x2（出院人數(shù)）、x5（平均住院天數(shù)）、x8（診斷符合率）、x9（搶救成功率）上因子載荷較大；F2在x3（病床利用率）、x4（病床周轉(zhuǎn)次數(shù)）上的因子載荷較大；F3在x6（治愈好轉(zhuǎn)率）、x7（病死率）上的因子載荷較大

我們可以推出：因子F1反映了該醫(yī)院醫(yī)療工作質(zhì)量各方面的情況，為綜合因子；F2反映了病床利用情況；F3反映了醫(yī)療水平的高低

將旋轉(zhuǎn)后的因子載荷與主成分分析的因子載荷矩陣比較可知：因子旋轉(zhuǎn)后，除F1的因子載荷仍分布多數(shù)指標(biāo)上外，其他2個(gè)因子的載荷明顯地集中到少數(shù)指標(biāo)上，說(shuō)明旋轉(zhuǎn)對(duì)因子載荷起到明顯的分離作用，使得各因子解釋的變量更加清晰。