基于R语言的因子分析（附代码）

一、因子分析数学推导

        用少数几个不可观测的随机变量（因子）去描述许多随机变量之间的协方差关系。主义因子分析可以视作主成分分析的一种推广，但因子是不可观察且不必是相互正交的变量。其基本思想是：根据相关性大小将变量分组，每组变量代表一个基本结构，反映观测到的相关性。R型因子分析可以用来研究变量之间的相关关系，Q型因子分析研究样品之间的相关关系，二者知识形式上的不同数学处理上是一样的，本文以R型因子为例展开说明。

1.1正交因子模型

        假定P维期望为,协方差矩阵为的随机变量X线性地依赖于少数几个不可观测的随机变量(公因子)和附加的方差源(特殊因子，误差)，那么因子模型可表作

  

其中为数学期望，为第i个变量在第j个因子上的载荷。

当以下条件成立时，我们称该模型为正交因子模型（F为公因子构成矩阵）：且

可推导知，

也就是Xi有m个公因子和一个特殊因子构成。公因子的贡献总和称为共同度，特殊因子的方差贡献为特殊度。注意到，在非退化的情况下模型是不唯一的（缺少约束时模型是不可识别的）。

1.2 因子模型的估计

        建立因子模型首先要估计载荷矩阵和特殊方差，常用的估计方法有主成分发、主因子法和最大似然法。

1.2.1 主成分法

        设为的特征值，为对应的标准正交化特征向量。这里可以对随机变量的协方差矩阵进行分解，并且假定特殊因子是可以在分解过程中忽略的。在实际中，常用样本的协方差矩阵来代替，考虑到变量的量纲差别，往往需要将数据标准化，这样求得的样本方差-协方差矩阵就是原来数据的相关系数矩阵R，并对R进行以上操作。而后根据各特征值的占比和阈值进行比较取定公因子。

1.2.2 主因子法

        作R的约相关矩阵并对此进行估计，具体过程不在文本赘述，其和主成分法的不同是：特殊方差汗的估计在不同步骤进行。如果我们希望求得拟合程度更高的节，可以采用迭代的方法，即将上一步中的特殊方差作为下一步的初始估计。

1.2.3极大似然法

   在模型假定下计算似然函数，记其极大似然估计为

注意对极大似然解，当因子数增加时，愿挨银子的估计载荷以及对x的贡献度发生改变，区别于前两种方法。

1.3因子正交旋转

        满足方差结构的因子模型并不唯一，模型的公因子和载荷矩阵不唯一，并且公因子和因子载荷矩阵作正交变换后并不改变公共度，乘因子载荷的正交变换和伴随的因子正交变换为因子正交旋转，从而得到更简单的结构。

1.4因子得分

        因子得分即对公因子的估计

1.4.1加权最小二乘法

        先将因子模型改写为后经过变化并由最小二乘法可得的估计。是集中通常用它们的某种估计来代替。

1.4.2回归法

        在正交因子模型中假设服从元正态分布，用回归预测的方法估计。

1.4.3综合因子得分

        以各因子的方差贡献率为权重，由各因子的下行组合得到综合评价指标函数：

二、案例分析

####导入数据 安装相应库###############
setwd('C:/Users/myq20170530/data')
data<-read.csv('eg6.1.csv',header = 1)
data<-scale(data)##数据标准化处理
cor(data)  ##观察各个变量之间的相关程度

install.packages('psych')##笔者mvstats库无法下载故而用psych代替
library('psych')
install.packages('GPArotation')
library('GPArotation')

####用主成分分析和极大似然法进行因子分析##########
fa(data,nfactors=2,fm='ml',rotate = 'none')
>Factor Analysis using method =  ml
Call: fa(r = data, nfactors = 2, rotate = "none", fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
     ML1  ML2   h2   u2 com
x1 -0.68 0.56 0.77 0.23 1.9
x2 -0.60 0.43 0.54 0.46 1.8
x3 -0.49 0.66 0.67 0.33 1.8
x4  0.92 0.10 0.85 0.15 1.0
x5  0.86 0.24 0.79 0.21 1.2
x6  0.88 0.27 0.85 0.15 1.2

                       ML1  ML2
SS loadings           3.40 1.07
Proportion Var        0.57 0.18
Cumulative Var        0.57 0.75
Proportion Explained  0.76 0.24
Cumulative Proportion 0.76 1.00

Mean item complexity =  1.5
Test of the hypothesis that 2 factors are sufficient.

df null model =  15  with the objective function =  4.28 with Chi Square =  205.97
df of  the model are 4  and the objective function was  0.08 

The root mean square of the residuals (RMSR) is  0.02 
The df corrected root mean square of the residuals is  0.04 

The harmonic n.obs is  52 with the empirical chi square  0.62  with prob <  0.96 
The total n.obs was  52  with Likelihood Chi Square =  3.64  with prob <  0.46 

Tucker Lewis Index of factoring reliability =  1.007
RMSEA index =  0  and the 90 % confidence intervals are  0 0.203
BIC =  -12.17
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                                   ML1  ML2
Correlation of (regression) scores with factors   0.97 0.89
Multiple R square of scores with factors          0.95 0.80
Minimum correlation of possible factor scores     0.89 0.59


> fac1=principal(data,2,rotate='none')
> fac1
Principal Components Analysis
Call: principal(r = data, nfactors = 2, rotate = "none")
Standardized loadings (pattern matrix) based upon correlation matrix
     PC1  PC2   h2   u2 com
x1 -0.79 0.42 0.81 0.19 1.5
x2 -0.73 0.40 0.70 0.30 1.5
x3 -0.64 0.63 0.81 0.19 2.0
x4  0.89 0.31 0.89 0.11 1.2
x5  0.81 0.47 0.87 0.13 1.6
x6  0.83 0.46 0.90 0.10 1.6

                       PC1  PC2
SS loadings           3.71 1.26
Proportion Var        0.62 0.21
Cumulative Var        0.62 0.83
Proportion Explained  0.75 0.25
Cumulative Proportion 0.75 1.00

Mean item complexity =  1.6
Test of the hypothesis that 2 components are sufficient.

The root mean square of the residuals (RMSR) is  0.06 
 with the empirical chi square  5.96  with prob <  0.2 

Fit based upon off diagonal values = 0.99


fac2=principal(data,2,rotate = 'varimax')
> fac2
Principal Components Analysis
Call: principal(r = data, nfactors = 2, rotate = "varimax")
Standardized loadings (pattern matrix) based upon correlation matrix
     RC1   RC2   h2   u2 com
x1 -0.32  0.84 0.81 0.19 1.3
x2 -0.29  0.78 0.70 0.30 1.3
x3 -0.07  0.90 0.81 0.19 1.0
x4  0.88 -0.35 0.89 0.11 1.3
x5  0.92 -0.18 0.87 0.13 1.1
x6  0.93 -0.20 0.90 0.10 1.1

                       RC1  RC2
SS loadings           2.66 2.31
Proportion Var        0.44 0.39
Cumulative Var        0.44 0.83
Proportion Explained  0.54 0.46
Cumulative Proportion 0.54 1.00

Mean item complexity =  1.2
Test of the hypothesis that 2 components are sufficient.

The root mean square of the residuals (RMSR) is  0.06 
 with the empirical chi square  5.96  with prob <  0.2 

Fit based upon off diagonal values = 0.99

###对结果进行可视化处理######################
plot(fac2$loadings,xlabel='Factor1',ylabel='Factor2')#输出因子载荷图

text(fac2$loadings,paste=("x",1:12,sep=""),cex=1.5)
fac2.plotdata<-fac2$scores
rownames(fac2.plotdata)<-unlist(name)
plot.text(fac2.plotdata)#因子得分图


plot(fac2$loadings,xlab="Factor1",ylab="Factor2")#各因子得分图和原坐标在因子上的方向图

 

 三、结果分析

        由因子载荷图可以看出后三个随机变量离第一个因子较近，前三者离第二个因子较近。

        从上述极大似然法和主成分法得到的因子分析结果可以看出：当m取2时极大似然法累计贡献率为74%而主成分法的累计贡献率为87%，说明主成分法效果比极大似然法好。

      从因子正交旋转的结果可以看出累计贡献率为82.87%。可见作正交旋转后因子的含义更清晰。

      本文内容主要来源于费宇先生所著《多元统计分析——基于R》一文，笔者才疏学浅，若有不妥谬误之处，还望温和批评指正。