基于R lavaan 进行SEM分析总结

前期对SEM模型的构建以及问卷题项设计工作已经完成，本部分从问卷数据收集以及处理完成后开始分析。

一理论部分

1.1验证性因子分析

验证性因子分析CFA的主要目的在于进行效度验证，它与探索性因子分析有所不同。
探索性因子分析主要是为了找出影响观测变量的因子个数，以及各个因子和各个观测变量之间的相关程度，推断模型的内部结果。而验证性因子分析的主要是决定事前定义因子的模型拟合实际数据的能力，以试图检验观测变量的因子个数和因子载荷是否与基于预先建立的理论的预期一致。具体区别可以点击探索性因子分析与验证性因子分析的区别参考

验证性因子分析主要工作有：信度与效度的分析

1.信度分析
信度，即可靠性，是指问卷的设计的可信程度。通常使用的衡量指标是Alpha系数，一般Alpha系数在0.9以上，可靠性比较好，在0.8-0.9是可接受的；在0.7-0.8则需要进行调整题项
2.效度分析：
一般常常使用的效度是聚合效度和区分效度
（1）聚合效度/收敛效度： 一般强调属于同一因子（指标）下的测量项落在同一因子下的有效程度。常常使用AVE和CR这两个指标进行分析
AVE，即平均提炼方差，反映了潜变量或同一因子下的各观测变量对该潜变量的平均差异解释力，AVE达到0.5以上为可以接受的水平
CR，即潜变量的组成信度，用来分析潜变量的各观测变量间的一致性，CR达到0.7以上表明组成信度较好；
有时还会要求每个潜变量的测量题项因子载荷系数（factor loading）值大于0.7
（2）区分效度/区别效度 ：强调不同因子下的题项确实不在同一因子下面的可靠程度。区分效度分析一般使用AVE根号值和相关分析结果进行对比，如果每个因子的AVE根号值均大于“该因子与其它因子的相关系数最大值”，则说明具有良好的区别效度。如下图来源于网络的一张图片，对角线洋红色的数据为因子的AVE根号值，其余数据为各个因子相关系数，因子1的AVE根号值是大于它与因子2和因子3，因子4的相关系数的，其余也是，可以看出具有良好的区别效度。

1.2 模型拟合分析

模型拟合分析常常使用共同方法偏差(CMV)的方式来进行测量。共同方法偏差，是指由于测量外部的某些因素导致数据出现集中的偏差。测量偏差可能是由于测量方法本身，测量模型本身，或者环境造成的差异。
对于模型拟合的指标有很多，包括如下：
1.χ2卡方拟合指数 ，原假设是模型协方差阵等于样本协方差阵。如果模型拟合的好，卡方值应该不显著。在这种情况下，数据拟合不好的模型被拒绝。
2.RMR 残差均方根 ，RMR 是样本方差和协方差减去对应估计的方差和协方差的平方和，再取平均值的平方根。一般RMR越小，拟合越好，理论上要小于0.08。
3.RMSEA 近似误差均方根，该值也是越小越好，RMSE理论上小于0.06，也可以放宽到0.08。
4.CFI 比较拟合指数，值在0-1之间。CFI 值越大越好，一般在0.9以上模型可接受
5.GFI 拟合优度指数，值在0-1之间，一般GFI 应该等于或大于0.90。
6.PGFI 简效拟合优度指数。它是简效比率(PRATIO，独立模式的自由度与内定模式的自由度的比率)乘以GFI。 PGFI 应该等于或大于0.90，越接近1越好。
7.PNFI 简效拟合优度指数，等于PRATIO乘以 NFI。 PNFI应该等于或大于0.90，越接近1越好。
8.NFI 规范拟合指数，变化范围在0和1间，当为1的时候标识完全拟合。一般NFI 小于0.90 表示需要重新设置模型。越接近1越好。
9.TLI Tucker-Lewis 系数， 也叫做Bentler-Bonett 非规范拟合指数 (NNFI)。TLI接近1表示拟合良好。
10.Chisqare/df卡方值与自由度的比值，该值越小越好，一般要小于2，放宽到3也是可以接受的。

由于指标众多，也有很多取舍，但是常常使用的重要参考指标为：Chisqare/df，RMSEA ，CFI

二使用R进行SEM模型简单分析

2.1需要的包：

#进行lavaan分析的包
library(lavaan)
#进行信效度辅助分析的包
library(semTools)
#进行sem模型画图的包
library(semPlot)
#进行导入sav文件的包
library(Hmisc)

2.2 录入数据与模型

#录入问卷数据
basic_data=spss.get("data.sav")
#录入构建的模型
basic_model<- readLines("modle.lav")

录入数据：也可以不读取sav格式数据，可以使用R自带的方式读取csv数据，格式如下：

basic_data<-read.csv("data.csv")

模型设计： 模型的简单语法为：

=~：用来定义潜变量
~ ：用来定义回归关系方程
~~：用来定义协方差

例如建立如下图所示的模型：

代码如下

#定义四个潜变量
f1=~x11+x12+x13
f2=~x21+x22+x23
f3=~x31+x32+x33
f4=~x41+x42+x43
#定义f3和f4之间的关系
f3~f1+f2
f4~f1+f2+f3
#f3，f2分别与f1的协方差，f4与f2的协方差
f3~~f1
f2~~f1
f4~~f2

模型可以写在代码中：一般格式为：

basic_model<-'
f1=~x11+x12+x13
f2=~x21+x22+x23
'

如果模型比较大，可以写到记事本中，保存为lav格式的文件读取也行，读取方法如本部分前面的代码
有关lavaan模型建模的详情可参考lavaan教程或者翻译版中文lavaan简明教程

2.3 使用lavaan进行分析

lavaan是SEM模型分析的核心包，在上述模型构建好之后，主要进行CFA分析，对信效度、模型拟合指数进行相关分析，可以使用cfa函数，也可以使用sem函数，分析结果是类似的。
1.分析代码

#进行验证性因子分析
cfa_fit<-cfa(basic_model,basic_data)
#计算AVE,CR等值（信效度分析）
reliability(cfa_fit)
# 模型拟合结果，fit.measures可显示模型对应指标，这里选择了“all"就会把所有的指标参数都会罗列出来
fitmeasures(cfa_fit,fit.measures="all",baseline.modle=null)
#显示总体概要结果，standardized=TRUE参数显示标准化结果
summary(cfa_fit,standardized=TRUE)
#如果想显示所有标准化结果，可使用如下语句
standardizedSolution(cfa_fit)
#进行路径分析画图
semPaths(cfa_fit,whatLables="std", style="lisrel",nCharNodes=0,nCharEdges=0,title=TRUE, layout="spring",edg.lable.cex=1)

在分析过程中并没有找到区别效度分析的函数，可能是个人对lavaan学习不太深入。这时可以使用两种方法做区别效度：
方法1： 使用smartPLS模型分析数据，在分析结果中，discriminant validity就是区别效度的结果，smartPLS的操作可参考smart PLS入门教程
方法2： 对lavaan的分析结果进行后期excel或者R分析的加工，利用得出的协方差和方差可以得出两两因子之间的相关系数，并自行编辑到结果表格中，相关系数的计算公式为：

其中Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差，本部分不再进行赘述

2.4 分析结果

1.信效度分析结果
代码如下：

#计算AVE,CR等值（信效度分析）
reliability(cfa_fit)

运行结果如下：

          f1        f2        f3        f4     total
alpha  0.8593373 0.8738389 0.8830737 0.9025800 0.9503724
omega  0.8611799 0.8735066 0.8865331 0.9021715 0.9632151
omega2 0.8611799 0.8735066 0.8865331 0.9021715 0.9632151
omega3 0.8626579 0.8721693 0.8903044 0.8984987 0.9621632
avevar 0.6744191 0.6974425 0.7229420 0.7553726 0.7147706

在上述结果中omega值为CR，其中alpha值以及AVE值都已列举出来，可以看出AVE基本都大于0.5，alpha值在0.8以上,CR值也在0.8以上，说明模型信效度良好
2.模型拟合结果分析
代码运行结果为：

  npar                fmin               chisq                  df              pvalue 161.000               4.573            3009.020            1270.000               0.000 baseline.chisq         baseline.df     baseline.pvalue                 cfi                 tli 20840.944            1378.000               0.000               0.911               0.903 nnfi                 rfi                 nfi                pnfi                 ifi 0.903               0.843               0.856               0.789               0.911 rni                logl   unrestricted.logl                 aic                 bic 0.911          -22405.964          -20901.454           45133.928           45745.094 ntotal                bic2               rmsea      rmsea.ci.lower      rmsea.ci.upper 329.000           45234.403               0.065               0.062               0.067 rmsea.pvalue                 rmr          rmr_nomean                srmr        srmr_bentler 0.000               0.114               0.114               0.050               0.050
srmr_bentler_nomean                crmr         crmr_nomean          srmr_mplus   srmr_mplus_nomean 0.050               0.051               0.051               0.050               0.050 cn_05               cn_01                 gfi                agfi                pgfi 149.046             152.999               0.731               0.696               0.648 mfi                ecvi 0.071              10.125

在上述指标中，CFI已经达到0.911，符合模型的预期，RMSEA的值为0.065，低于0.08，模型的拟合效果还是很理想的，chisqurar/df=2009.02/1270<3,拟合结果还是比较好的
3 模型概述
模型代码如下：

#显示总体概要结果
summary(cfa_fit,standardized=TRUE)

执行结果如下：

 Estimator                                         ML#使用最大似然估计法Optimization method                           NLMINBNumber of free parameters                         30Number of observations                           329Model Test User Model:Test statistic                               184.131Degrees of freedom                                48P-value (Chi-square)                           0.000Parameter Estimates:Information                                 ExpectedInformation saturated (h1) model          StructuredStandard errors                             StandardLatent Variables:#显示载荷因子情况，可以看出基本都大于0.7Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.allf1 =~                                                                 x11               1.000                               1.091    0.814x12               1.004    0.065   15.528    0.000    1.095    0.776x13               1.093    0.060   18.066    0.000    1.193    0.874f2 =~                                                                 x21               1.000                               1.239    0.880x22               0.958    0.055   17.414    0.000    1.187    0.778x23               0.911    0.044   20.536    0.000    1.129    0.857f3 =~                                                                 x31               1.000                               1.195    0.794x32               1.006    0.057   17.560    0.000    1.203    0.846x33               1.099    0.056   19.462    0.000    1.313    0.911f4 =~                                                                 x41               1.000                               1.400    0.856x42               0.969    0.046   21.081    0.000    1.356    0.879x43               0.853    0.041   20.936    0.000    1.194    0.876Regressions:#模型路径系数，可以看出拟合的双尾检验有些不太显著Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.allf3 ~                                                                  f1                0.247    0.135    1.824    0.068    0.225    0.225f2                0.614    0.123    5.008    0.000    0.636    0.636f4 ~                                                                  f1               -0.117    0.126   -0.933    0.351   -0.092   -0.092f2                0.075    0.128    0.589    0.556    0.067    0.067f3                1.136    0.106   10.709    0.000    0.970    0.970Covariances:#协方差结果Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.allf1 ~~                                                                 f2                1.201    0.120   10.026    0.000    0.888    0.888Variances:#方差结果Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all.x11               0.607    0.060   10.036    0.000    0.607    0.338.x12               0.792    0.074   10.697    0.000    0.792    0.398.x13               0.441    0.054    8.184    0.000    0.441    0.237.x21               0.448    0.051    8.710    0.000    0.448    0.226.x22               0.920    0.083   11.046    0.000    0.920    0.395.x23               0.462    0.048    9.537    0.000    0.462    0.266.x31               0.836    0.073   11.394    0.000    0.836    0.369.x32               0.575    0.054   10.677    0.000    0.575    0.284.x33               0.351    0.042    8.347    0.000    0.351    0.169.x41               0.718    0.070   10.279    0.000    0.718    0.268.x42               0.541    0.056    9.619    0.000    0.541    0.227.x43               0.434    0.045    9.731    0.000    0.434    0.233f1                1.191    0.137    8.678    0.000    1.000    1.000f2                1.534    0.155    9.885    0.000    1.000    1.000.f3                0.414    0.061    6.732    0.000    0.290    0.290.f4                0.176    0.050    3.485    0.000    0.090    0.090

4.进行画图
进行画图代码如下：

semPaths(cfa_fit,what="std", style="lisrel",nCharNodes=0,nCharEdges=0,title=TRUE, layout="spring",edg.lable.cex=1,fade=FALSE)

代码运行的结果为：

从得出的图形可以看出，使用semPath函数得出的模型并不像理想的那样美观，如果潜变量特别多的话，图形会变得非常糟糕，这时建议还是使用专业软件进行制图，最简单的，也可以借用ppt手画吧。具体的使用教程可以参考semPlot教程

三代码分析过程出现的问题

在使用semPlot包出现了两个问题：
安装完成后不能使用：检查结果是出现底层的支持包未安装，如果出现这种情况，可以根据提示把支持的包先安装好，然后再重新安装semPlot包
安装加载完成后不能调用函数：提示不能识别支持文件还有路径，具体错误很长，类似如下：

Error in gzfile(file, "wb") : 无法打开链结
此外: Warning messages:
1: In qgraph::qgraph(Edgelist, labels = nLab, bidirectional = Bidir,  :The following arguments are not documented and likely not arguments of qgraph and thus ignored: edg.lable.cex

此时最可能的原因是你在安装包的存储路径中含有中文路径，导致R不能识别。由于我的包安装在用户目录下不能更改路径名称，所以采用了强大的网友的解决帮助方法，在使用该语句前，添加一句话，dev.new()，这里新建一个画布，不用把生成的图像存储到目标路径中，所以就不用读取目录路径，详情可以参考不能画图的解决方法

dev.new()
semPaths(cfa_fit,what="std", style="lisrel",nCharNodes=0,nCharEdges=0,title=TRUE, layout="spring",edg.lable.cex=1,fade=FALSE)

但是最后semPlot包使用的频率也不太高
附录完整代码：

#引入包
library(lavaan)
library(semTools)
library(semPlot)
library(Hmisc)
library(semPlot)#录入数据
basic_data=spss.get("data.sav")
#录入模型
basic_model<- readLines("modle.lav")#进行验证性因子分析
cfa_fit<-cfa(basic_model,basic_data)
#计算AVE,CR等值（信效度分析）
reliability(cfa_fit)
# 模型拟合指数
fitmeasures(cfa_fit,fit.measures="all",baseline.modle=null)
#显示总体概要结果
summary(cfa_fit,standardized=TRUE)
#进行路径分析画图
dev.new()
semPaths(cfa_fit,whatLables="std", style="lisrel",nCharNodes=0,nCharEdges=0,title=TRUE, layout="spring",edg.lable.cex=1)
dev.off()