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如何选取）
下面给出主成分法的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个因子的载荷明显地集中到少数指标上，说明旋转对因子载荷起到明显的分离作用，使得各因子解释的变量更加清晰。