R语言做结构方程模型入门

这是用R做结构方程的案例, 最下方有视频教程.

目的

这篇教程目的是介绍如何使用R语言的lavaan包来做结构方程模型，
如果你不清楚什么是结构方程（SEM）， 建议你先看这个介绍视频。
当然， 我们会在教程的前面简要介绍结构方程的概念， 以便于我们教程的完整性。

案例介绍

这是一个研究学生的学业成绩影响因素的研究，
目前数据已经采集技术， 你可以在这里下载这个数据,
这个数据有9个观测指标（9个变量）：Motivation（冬季）, Harmony, Stability, Negative Parental Psychology, SES, Verbal IQ, Reading, Arithmetic and Spelling（动机、和谐、稳定、消极的父母心理、SES、语言智力、阅读、算术和拼写）。研究员假设三个潜变量：Adjustment, Risk, Achievement， 他们的测量指标如下， 其中包含了变量名及其解释：

Adjustment
- motiv Motivation
- harm Harmony
- stabi Stability

Risk
- ppsych (Negative) Parental Psychology
- ses SES
- verbal Verbal IQ

Achievement
- read Reading
- arith Arithmetic
- spell Spelling

模型表示

对于任何结构方程研究， 你最好事先绘制理论模型图， 下面就是我们的假定模型图：

根据理论模型， 如果用数学符号表示， 我们可以绘制具有数学意义的模型图：

我们先来解释一下图中的符号的意义：

• observed variable: 显变量或者叫观测指标， 就是在你的数据中有一列数据可以代表该变量， 图中用矩形表示
• latent variable: 潜变量， 数据中不存在， 但是可以在模型中构建出来，图中用椭圆表示
• exogenous variable: 外生变量， 可以认为是自变量， 就是在模型中没有变量预测它， 它用于预测其他变量， 它既可以是显变量$x$也可以是潜变量$\xi$
• endogenous variable: 内生变量， 可以认为是因变量， 只要有箭头指向它， 它就可以被成为是内生变量， 即可以是显变量$y$也可以是潜变量$\eta$
• measurement model: 测量模型， 多个显变量和一个潜变量构成的结构， 比如图中$\xi_1$与$x_1 x_2 x_3$构成的结构
• factor: 因子， 一个潜变量， 由多个显变量所定义， 例如图中的$\eta_1 \eta_2 \xi_1$
• loading: 因子载荷， 因子和观测指标之间的连线， 可以用字母$\lambda$表示
• structural model: 由自变量和因变量构成的结构， 比如$\eta_1 \eta_2 \xi_1$及其连线构成的结构
• regression path: 回归路径， 自变量和因变量之间的连线， 用字母$\gamma\$

R代码

有了数据和模型， 下一步就是如何用R来表示模型和估计模型。

模型的表示

R中最常用的SEM包就是lavaan , 它定义了一些模型表示的方法, 列在这里:

• ~ 波浪号, 表示观测变量的回归关系, (e.g., y ~ x)
• =~ 等号波浪号, 用于表示测量模型, (e.g., f =~ q + r + s), f表示潜变量, q r s 表示测量指标
• ~~ 协方差 (e.g., x ~~ x)
• ~1 截距或者均值, 比如x ~ 1估计了变量x的均值
• * 固定参数, 比如f =~ 1*q意思是f对q的因子载荷固定为1
• NA* 自由参数或者载荷 (e.g., f =~ NA*q)
• a* 参数标签, 比如f =~ a*q, 意思是f对q的因子载荷是a, 在设定模型constrants时非常有用

参考这个模型图, 我们有三个潜变量, 每个潜变量有三个测量指标, 测量模型的表示为:

1
2
3

adjust =~ motiv + harm + stabi
risk =~ verbal + ses + ppsych
achieve =~ read + arith + spell

结构方程需要你清楚模型中有多少个因变量, 因为有多少因变量, 就可以写几行回归方程, 识别因变量最好的方法就是看有几个变量被箭头指向,
根据模型图, 我们得出两个回归方程, 结构模型的表示:

1
2

adjust ~ risk 
achieve ~ adjust + risk

R 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# m6b就是模型代码, 它在R中数据类型是字符串
m6b <- '
# measurement model
adjust =~ motiv + harm + stabi
risk =~ verbal + ses + ppsych
achieve =~ read + arith + spell
# regressions
adjust ~ risk 
achieve ~ adjust + risk
'

# 使用SEM函数, 第一个参数是模型代码, 第二个参数是数据(DataFrame)
fit6b <- sem(m6b, data=dat)

# 使用summary输出结果
summary(fit6b, standardized=TRUE, fit.measures=TRUE)

实战

安装和加载库

1
2

install.packages("lavaan", dependencies=TRUE)
library(lavaan)

下载数据

这是一个网络数据, 你可以下载到本地, 点这里,
你也可以直接加载到R中.

1

dat <- read.csv("https://stats.idre.ucla.edu/wp-content/uploads/2021/02/worland5.csv")

1

head(dat) # 查看数据前几行

根据前面测量模型的代码, 我们知道, 各个变量的意义,

• adjust的测量指标是: motiv + harm + stabi
• risk 的测量指标是: verbal + ses + ppsych
• achieve 的测量指标是: read + arith + spell

估计模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# m6b就是模型代码, 它在R中数据类型是字符串
m6b <- '
# measurement model
adjust =~ motiv + harm + stabi
risk =~ verbal + ses + ppsych
achieve =~ read + arith + spell
# regressions
adjust ~ risk 
achieve ~ adjust + risk
'

# 使用SEM函数, 第一个参数是模型代码, 第二个参数是数据(DataFrame)
fit6b <- sem(m6b, data=dat)

# 使用summary输出结果
summary(fit6b, standardized=TRUE, fit.measures=TRUE)

展示模型的概览:

1

show(fit6b)

• 获取模型拟合参数

1
2
3

fits = fitMeasures(fit6b)

fits['rmsea']

1

fits[c('chisq', 'df',  'pvalue', 'cfi' ,'rmsea', 'srmr', 'aic')]

计算中介效应和总效应

我们注意到, 这是一个中介模型, 我们想要检验中介效应以及总效应是否显著, 可以使用下面的代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# m6b就是模型代码, 它在R中数据类型是字符串
m6b <- '
# measurement model
adjust =~ motiv + harm + stabi
risk =~ verbal + ses + ppsych
achieve =~ read + arith + spell
# regressions
adjust ~ a*risk # 设置参数a
achieve ~ b*adjust + c*risk # 设置参数b和c

# 中介效应
indirect := a * b
total := a*b + c
'

# 使用SEM函数, 第一个参数是模型代码, 第二个参数是数据(DataFrame)
fit6b <- sem(m6b, data=dat)

# 使用summary输出结果
summary(fit6b, standardized=TRUE, fit.measures=TRUE)

视频教程

注意
统计咨询请加QQ 2726725926, 微信 shujufenxidaizuo, SPSS统计咨询是收费的, 不论什么模型都可以, 只限制于1个研究内.
跟我学统计可以代做分析, 每单几百元不等.
本文由jupyter notebook转换而来, 您可以在这里下载notebook
可以在微博上@mlln-cn向我免费题问
请记住我的网址: mlln.cn 或者 jupyter.cn