用R做验证性因素分析（CFA）

用R做因素分析是因为我那个Mplus出问题了。

因素分析是潜变量分析的一种，目的是用潜变量（latent variables）来解释外显变量（manifest variables）的共同变异。验证性因素分析（Confirmatory Factor Analysis, CFA）用于检验研究者对于外显变量关系之间的理论模型，在编制测验的研究中很常用。潜变量就是不能直接测量的变量，或者叫因素、因子（factor）。外显变量是能直接获得测量值的变量，最常见的外显变量就是测验得分。如果我们认为某几个外显变量的测量值能够反映出一个潜变量的水平，就把这几个外显变量称为这个潜变量的指标（indicator），比如智力测验的分数就可以作为智力的指标。还有一对概念叫做内生变量（endogenous variables）和外生变量（exogenous variable）。内生变量的变异来源在模型内（模型中的其他变量），外生变量的变异来源在模型外。在因素分析中一般将外显变量作为内生变量（变异来源于潜变量），将潜变量作为外生变量。在其他潜变量建模中可能会出现外生外显变量、内生外显变量、外生潜变量和内生潜变量。

在R中进行验证性因素分析主要用lavaan这个包，但是lavaan中没有绘制路径图的函数，所以还需要semPlot包。接下来的示例将用八项SPSS焦虑问卷（SPSS Anxiety Questionnaire 8, SAQ-8）的数据来进行演示。这个问卷主要测量学生对于统计学的焦虑，共有8个项目，计分方式是5点评分。

I dream that Pearson is attacking me with correlation coefficients.

I don’t understand statistics.

I have little experience with computers.

Computers are useful only for playing games.

I have never been good at mathematics.

My friends are better at statistics than me.

All computers hate me.

I did badly at mathematics at school.

这是数据集的下载地址。数据的格式是sav。

装载包、导入数据

R的默认功能中没有读取sav文件的函数，所以我们需要rio包中的import()。这个函数基本上能够将所有格式的数据文件读取为数据框，而且自动生成列名，非常方便。dplyr是一个用来进行数据初步处理的包，其中定义了管道运算符%>%，作用是将左边的对象传递给右边的函数作为参数。管道运算符能提高代码的可读性，直接从左向右看就可以，不用一层一层拆套娃。比如以下两行代码的作用完全相同，都是先计算数据框saq的相关矩阵再将所有的相关系数保留两位有效数字。

saq %>% cor() %>% round(2)
round(cor(saq),2)

这个例子只用到了两个函数所以效果不够明显，然而在进行数据整理的时候可能需要连续进行处理缺失值、筛选变量、创建新变量、修改变量名、给个案分组、生成分组摘要等一系列操作，使用管道运算就很有必要了。

装载包

library(lavaan)
library(semPlot)
library(rio)
library(dplyr)

导入数据

saq = import("SAQ8.sav") %>% na.omit() #删掉有缺失值的行

检验数据是否适合因素分析

事实上大多数研究都不做这一步。

简单直观的方法就是看相关矩阵。

> (cor = saq %>% cor() %>% round(2)) #给赋值语句加上括号就会自动输出赋值对象
      q01   q02   q03   q04   q05   q06   q07   q08
q01  1.00 -0.10 -0.34  0.44  0.40  0.22  0.31  0.33
q02 -0.10  1.00  0.32 -0.11 -0.12 -0.07 -0.16 -0.05
q03 -0.34  0.32  1.00 -0.38 -0.31 -0.23 -0.38 -0.26
q04  0.44 -0.11 -0.38  1.00  0.40  0.28  0.41  0.35
q05  0.40 -0.12 -0.31  0.40  1.00  0.26  0.34  0.27
q06  0.22 -0.07 -0.23  0.28  0.26  1.00  0.51  0.22
q07  0.31 -0.16 -0.38  0.41  0.34  0.51  1.00  0.30
q08  0.33 -0.05 -0.26  0.35  0.27  0.22  0.30  1.00
> cor[cor!=1] %>% abs() %>% range() #|r|的最大值和最小值，先去掉了主对角线上的1
[1] 0.05 0.51

一般来说相关系数的绝对值达到0.3就认为是中等相关。粗略地看一下可以认为项目分数之间存在相关。

如果用统计方法来检验的话就是进行Bartlett球形度检验（Bartlett's test of sphericity）和KMO检验（ Kaiser-Meyer-Olkin Measure of Sampling Adequacy）。Bartlett检验在方差分析中被用来检验方差齐性，其实际作用是比较一个相关矩阵和单位矩阵（主对角线元素是1，其余元素是0）是否有显著差异，有显著差异就说明变量之间显著相关。KMO检验则是比较变量之间的偏相关系数和相关系数的差异，KMO统计量能够反映可能由共同因子造成的变异的比例，越接近1越好，一般大于0.8就非常足够了。

The Kaiser-Meyer-Olkin Measure of Sampling Adequacy is a statistic that indicates the proportion of variance in your variables that might be caused by underlying factors. High values (close to 1.0) generally indicate that a factor analysis may be useful with your data. If the value is less than 0.50, the results of the factor analysis probably won't be very useful.

Bartlett's test of sphericity tests the hypothesis that your correlation matrix is an identity matrix, which would indicate that your variables are unrelated and therefore unsuitable for structure detection. Small values (less than 0.05) of the significance level indicate that a factor analysis may be useful with your data.

来源：IBN Knowledge Center

以上两个检验可以用psych包中的corr.test.bartlett()和KMO()这两个函数来完成。这两个函数的对象都是相关矩阵。

> cortest.bartlett(cor,n = nrow(saq)) #n是样本容量
$chisq
[1] 4159.883

$p.value
[1] 0

$df
[1] 28
> KMO(cor)
Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor)
Overall MSA =  0.82 
MSA for each item = 
 q01  q02  q03  q04  q05  q06  q07  q08 
0.84 0.67 0.82 0.85 0.87 0.76 0.78 0.88

Bartlett球形度检验结果显著，KMO统计量Overall MSA = 0.82，说明数据适合因素分析。

自由度和模型识别

自由度决定了模型能否识别，也就是模型中的参数能不能计算估计值。因子分析模型的自由度等于数据提供的信息数减去模型中需要估计的参数数量。数据中的信息数量和模型要估计的参数数量的关系就像是一次方程的数量和未知数的数量的关系，要求解n个未知数需要至少n个不同的方程，要估计n个参数也需要至少n个信息。模型识别情况有以下三种：

自由度小于0，模型识别不足（under-identified），参数无法估计
自由度等于0，模型完全识别（just-identified）,参数可以估计
自由度大于0，模型过度识别（over-identified），参数可以估计

接下来说明这两个数量是怎么计算的。

因子分析的起点是协方差矩阵，或者相关矩阵（相关系数是标准化的协方差）。

      q01   q02   q03   q04   q05   q06   q07   q08
q01  1.00 -0.10 -0.34  0.44  0.40  0.22  0.31  0.33
q02 -0.10  1.00  0.32 -0.11 -0.12 -0.07 -0.16 -0.05
q03 -0.34  0.32  1.00 -0.38 -0.31 -0.23 -0.38 -0.26
q04  0.44 -0.11 -0.38  1.00  0.40  0.28  0.41  0.35
q05  0.40 -0.12 -0.31  0.40  1.00  0.26  0.34  0.27
q06  0.22 -0.07 -0.23  0.28  0.26  1.00  0.51  0.22
q07  0.31 -0.16 -0.38  0.41  0.34  0.51  1.00  0.30
q08  0.33 -0.05 -0.26  0.35  0.27  0.22  0.30  1.00

因为主对角线两侧的值是对称的，所以K个变量的相关矩阵中共有 $K(K+1)$ 个不同的值，这就是数据能够提供的用来建模的信息数量。

然后我们来看模型中有多少个需要估计的参数值。以下模型中f1是因子，q1、q2、q3分别是f1的三个指标，是残差（不能被模型解释的部分，相当于回归分析中的）。路径图中一般用矩形表示外显变量，用圆形表示潜变量。

CFA-单因素3指标路径图

模型中的参数可以分为以下几类

因素负荷（factor loading）：因子和指标之间的红色箭头，用 $\lambda$ 表示
残差和指标的相关：残差和指标之间的红色箭头，固定值为1
因子之间的协方差：因子间的蓝色双向箭头，用 $\psi$ 表示，因子与自身的协方差就是因子的方差
残差方差：残差的蓝色双向箭头，用 $\theta$ 表示
残差之间的协方差：因为假定残差之间无相关，所以都为0，在图中没有标示

以上参数中需要估计的有因素负荷因子方差和残差方差，总共7个。那么问题来了，我们只有 $3×4/2=6$ 个信息，所以这个模型识别不足。解决方法有两种：

marker method 将第一个因素负荷固定为1，这也是lavaan的默认方法
variance standardization method 将因子的方差固定为1，可以估计所有的因素负荷

总之就是通过固定一个参数值来减少一个需要估计的参数，这样需要估计的参数数量变成6个，模型变成了完全识别。这两种方法的另一个作用是为因素负荷指定单位，也是进行CFA的必要条件。要求单因子模型中每个因素的指标不应少于3个，多于1个因子时每个因子至少有2个指标。

因为标准的CFA模型假定每个指标只在一个因子上有负荷，一般都是过度识别模型。而饱和模型假设所有变量都存在相关（所有指标在所有因子上都有负荷），所以对数据的拟合是完全的，但是对于了解变量之间的关系没有帮助。

建模和拟合度指标

首先我们看一下探索性因素分析（exploratory factor analysis, EFA）的结果来初步确定因子结构。EFA也是个大坑，就不讲了。

SAQ8 EFA

从这个路径图我们可以看出：

有3个因子，因子1对应项目1、4、5、8，因子2对应项目6、7，因子3对应项目2、3
因子之间存在相关

根据这两点我们可以开始定义CFA模型了。

在lavaan中定义模型主要用到以下几个符号：

=~用来定义潜变量，比如f1 =~ y1 + y2 +y3代表因子f1的指标是y1、y2和y3
~ 回归运算符，相当于回归方程中的等号，左侧是因变量，右侧是自变量，比如f1~f2+f3
~~用来表示方差和协方差，比如f1~~f2表示f1和f2的协方差，f1~~f1表示f1因子的方差
~1为变量指定截距

在lavaan中要将模型描述保存为一个字符串变量（用一对单引号括起来），再将这个变量和数据传递给cfa()函数就可以了。

模型描述如下

#因子之间的相关可以省略，因为cfa()默认使用marker method，将每个因子的第一个负荷固定为1而且因子之间相关
#如果理论模型中因子相互独立则需要使用variance standardization method
CfaModel = '
f1 =~ q01 + q04 + q05 + q08
f2 =~ q06 + q07
f3 =~ q02 + q03
'
cfafit = cfa(model = CfaModel,data = saq)
summary(cfafit,
        fit.measures = T, #显示模型拟合指标
        standardized = T) #显示标准化的估计值
dev.new() #新建一个图形窗口
#绘制路径图
semPaths(cfafit,what = "std", #显示标准化的估计值，如果要显示原始估计值则 what = "par"
         rotation = 2, #将潜变量置于左侧，显变量置于右侧
         edge.color = "black", #箭头和线条颜色设为黑色
         esize = 0.5, #箭头的粗细
         edge.label.cex = 1, #所有值的字号
         exoVar = F #不显示外生变量的方差)

以下是结果

> summary(cfafit, fit.measures = T, standardized = T)
lavaan 0.6-5 ended normally after 49 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of free parameters                         19
                                                      
  Number of observations                          2571
                                                      
Model Test User Model:
                                                      
  Test statistic                                58.271
  Degrees of freedom                                17
  P-value (Chi-square)                           0.000

Model Test Baseline Model:

  Test statistic                              4164.572
  Degrees of freedom                                28
  P-value                                        0.000

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.990
  Tucker-Lewis Index (TLI)                       0.984

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)             -26381.599
  Loglikelihood unrestricted model (H1)     -26352.464
                                                      
  Akaike (AIC)                               52801.198
  Bayesian (BIC)                             52912.387
  Sample-size adjusted Bayesian (BIC)        52852.018

Root Mean Square Error of Approximation:

  RMSEA                                          0.031
  90 Percent confidence interval - lower         0.022
  90 Percent confidence interval - upper         0.040
  P-value RMSEA <= 0.05                          1.000

Standardized Root Mean Square Residual:

  SRMR                                           0.018

Parameter Estimates:

  Information                                 Expected
  Information saturated (h1) model          Structured
  Standard errors                             Standard

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  f1 =~                                                                 
    q01               1.000                               0.519    0.628
    q04               1.278    0.051   25.027    0.000    0.664    0.700
    q05               1.100    0.049   22.665    0.000    0.571    0.592
    q08               0.841    0.042   19.979    0.000    0.437    0.501
  f2 =~                                                                 
    q06               1.000                               0.662    0.590
    q07               1.449    0.073   19.982    0.000    0.959    0.870
  f3 =~                                                                 
    q02               1.000                               0.278    0.327
    q03               3.758    0.428    8.774    0.000    1.046    0.973

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  f1 ~~                                                                 
    f2                0.222    0.015   14.813    0.000    0.645    0.645
    f3               -0.079    0.010   -7.956    0.000   -0.548   -0.548
  f2 ~~                                                                 
    f3               -0.082    0.011   -7.518    0.000   -0.443   -0.443

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .q01               0.415    0.015   28.517    0.000    0.415    0.606
   .q04               0.459    0.018   24.952    0.000    0.459    0.510
   .q05               0.604    0.020   29.777    0.000    0.604    0.649
   .q08               0.570    0.018   32.133    0.000    0.570    0.749
   .q06               0.820    0.030   27.571    0.000    0.820    0.652
   .q07               0.294    0.041    7.217    0.000    0.294    0.242
   .q02               0.647    0.020   32.967    0.000    0.647    0.893
   .q03               0.061    0.109    0.559    0.576    0.061    0.053
    f1                0.270    0.018   15.350    0.000    1.000    1.000
    f2                0.438    0.033   13.392    0.000    1.000    1.000
    f3                0.077    0.012    6.501    0.000    1.000    1.000

结果中要注意的有以下几个内容

自由参数数量和自由度

自由参数数量（Number of free parameters）也就是模型中要估计的参数数量，这个模型中共有8个因素负荷、8个残差方差、3个因子间的协方差和3个因子方差，其中有3个因素负荷（每个因子的第一个）被固定为1，所以自由参数数量是 $8+8+3-3=19$ 。

数据提供的信息总量是 $8×9/2=36$ ，因此模型自由度 $df=36-19=17$ 。

卡方检验结果

在Model Test User Model中，由于卡方检验对于样本容量很敏感，因此样本容量大于400（在社会心理学研究中很常见的情况）时一般都会得到显著的结果。卡方检验结果对于模型评价基本没有价值，但是论文中一般都会报告。

相对拟合指标

评价模型拟合度的相对指标主要是CFI（Comparative Fit Index）和TLI（Tucker Lewis Index），要解释这两个指数的意义我们要先知道什么是基线模型（Baseline Model）。基线模型中没有任何的协方差，也就是对与变量之间的关系没有假设。也就是这样：

SAQ8 baseline model

因为模型中没有假定指标的关系，所以路径图中没有潜变量。可以认为基线模型是拟合最差的模型，而相对拟合指标就是的模型相比基线模型的改进程度。饱和模型的相对拟合指标都是1，所以CFI和TLI值越接近1说明拟合度越好，一般大于0.9就认为模型拟合良好。UCLA的CFA教程认为，因为CFI和TLI高度相关，在报告结果时应该只报告其中一个。

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.990
  Tucker-Lewis Index (TLI)                       0.984

绝对拟合指标

RMSEA（Root Mean Square Error of Approximation）反映了模型误差的大小（一种非常不准确但是直观的解释）。一般认为说明拟合良好，说明拟合很差（不可接受）。

Root Mean Square Error of Approximation:

  RMSEA                                          0.031
  90 Percent confidence interval - lower         0.022
  90 Percent confidence interval - upper         0.040
  P-value RMSEA <= 0.05                          1.000

因素负荷

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  f1 =~                                                                 
    q01               1.000                               0.519    0.628
    q04               1.278    0.051   25.027    0.000    0.664    0.700
    q05               1.100    0.049   22.665    0.000    0.571    0.592
    q08               0.841    0.042   19.979    0.000    0.437    0.501
  f2 =~                                                                 
    q06               1.000                               0.662    0.590
    q07               1.449    0.073   19.982    0.000    0.959    0.870
  f3 =~                                                                 
    q02               1.000                               0.278    0.327
    q03               3.758    0.428    8.774    0.000    1.046    0.973

Estimate一列是用原始数据计算出的因素负荷，可以看到每个因子的第一个负荷被固定为1。Std.lv一列是将回归方程的自变量（X）标准化之后得到的回归系数，Std.all则是将回归方程中的自变量和因变量都标准化之后得到的回归系数。如何根据因素负荷来删除或保留项目没有统一的标准，一般要结合理论和数据来决定。

路径图

SAQ8 CFA

图中每个箭头上的数字是标准化的系数，系数越大则线条和字的颜色越深（可以看到q03的残差方差接近透明）。虚线箭头代表marker method中这个因素负荷被固定为1。

实际上报告CFA结果并没有统一的标准，一般都会报告卡方检验结果、CFI、TLI和RMSEA等拟合指标、因素负荷以及路径图。

用R做验证性数据分析