SEM in R by lavaan
吴欣益
2019-10-27
library(tidyverse)
library(lavaan)
library(semPlot)ways
- 先与一般线性回归做比较看看结果
- 完整的做了一次SEM，并输出结果
- 另外一些要求下的输出值
首先验证了只有两个变量做回归时与线性回归结果同
ols <- read_csv("./ols.csv")
ols %>% count(wa4)ols1 <- ols %>%
mutate(
wa4_1 = case_when(wa4 == 1 ~ 3,
wa4 == 0 ~ 1)
)
lm(wd2_years ~ wa4_1,data = ols1)
#>
#> Call:
#> lm(formula = wd2_years ~ wa4_1, data = ols1)
#>
#> Coefficients:
#> (Intercept) wa4_1
#> 15.1364 0.4615
model <- '
# latent variables
# F1 =~ wz301+ care_ave + wk6
# max_edu =~ wf701 + wh9 +F1
# finc_per =~ wf701 + F1 + wh9
# wf701 =~ F1 + wh9
# regressions
wd2_years ~ wa4_1
# residual covariances
# max_edu ~~ finc_per
'
fit <- sem(model,data=ols1)
summary(fit)
#> lavaan 0.6-5 ended normally after 13 iterations
#>
#> Estimator ML
#> Optimization method NLMINB
#> Number of free parameters 2
#>
#> Number of observations 3768
#>
#> Model Test User Model:
#>
#> Test statistic 0.000
#> Degrees of freedom 0
#>
#> Parameter Estimates:
#>
#> Information Expected
#> Information saturated (h1) model Structured
#> Standard errors Standard
#>
#> Regressions:
#> Estimate Std.Err z-value P(>|z|)
#> wd2_years ~
#> wa4_1 0.462 0.056 8.175 0.000
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> .wd2_years 8.093 0.186 43.405 0.000正式做SEM
接下来导入要做模型的数据
city_data <- read_csv("./city_data.csv")可以做CFA，（潜变量之间没回归），也可以做SEM，
#先将准备做的模模型关系写出来
#latent variablen那一块是潜变量和观测变量之家的关系用 =~，定义潜变量
#regression那一块是潜变量之间的关系用 ~，回归方程组，可以有潜变量
#residual covariances那是观测变量之间拉相关 用 ~~ ，方差和协方差
#只有截距项的回归方程用~1, 比如y1 ~ 1
model <- '
# latent variables
f1 =~ wz301+ care_ave + wk6
# regressions
wh9 ~ f1 + max_edu + wf701 +finc_per
f1 ~ max_edu + wf701 + finc_per
wf701 ~ max_edu + finc_per
# residual covariances
#截距项 ##########################lavaan包默认cfa()和sem()函数将潜变量截距固定为0，否则模型不能被估计
'
fit <- sem(model,data=city_data)#sem()和cfa()这连个参数结果相同，以后可能不同一般常看的结果
- summary()：结果概览
- parameterEstimates(fit)：参数估计
- standardizedSolution()：标准化参数估计
- fitted()：拟合协方差矩阵和均值向量
- resid()：拟合模型残差，[观测协方差矩阵/均值向量]和[隐含协方差矩阵/均值向量]的差
- fitMeasures()：返回所有拟合度量指标
- fitMeasures(fit, “cfi”) # 提取某个具体指标，例如cfi
- inspect()：查看拟合的lavaan对象具体信息
summary(fit)#前6行为头部，包含版本号，收敛情况，迭代次数，观测数，用来计算参数的估计量,模型检验统计量，自由度和相关的p值
#> lavaan 0.6-5 ended normally after 67 iterations
#>
#> Estimator ML
#> Optimization method NLMINB
#> Number of free parameters 17
#>
#> Number of observations 854
#>
#> Model Test User Model:
#>
#> Test statistic 31.643
#> Degrees of freedom 8
#> P-value (Chi-square) 0.000
#>
#> Parameter Estimates:
#>
#> Information Expected
#> Information saturated (h1) model Structured
#> Standard errors Standard
#>
#> Latent Variables:
#> Estimate Std.Err z-value P(>|z|)
#> f1 =~
#> wz301 1.000
#> care_ave 1.159 0.186 6.243 0.000
#> wk6 2.030 0.328 6.195 0.000
#>
#> Regressions:
#> Estimate Std.Err z-value P(>|z|)
#> wh9 ~
#> f1 7.953 2.072 3.838 0.000
#> max_edu -0.049 0.081 -0.608 0.543
#> wf701 0.038 0.022 1.718 0.086
#> finc_per -0.035 0.089 -0.390 0.697
#> f1 ~
#> max_edu 0.037 0.005 6.983 0.000
#> wf701 0.007 0.002 4.055 0.000
#> finc_per 0.001 0.009 0.150 0.881
#> wf701 ~
#> max_edu 0.187 0.078 2.396 0.017
#> finc_per -0.296 0.181 -1.636 0.102
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> .wz301 0.398 0.022 18.144 0.000
#> .care_ave 0.458 0.026 17.608 0.000
#> .wk6 1.460 0.082 17.749 0.000
#> .wh9 8.156 0.809 10.084 0.000
#> .wf701 65.890 3.189 20.664 0.000
#> .f1 0.037 0.011 3.288 0.001
Est <- parameterEstimates(fit, ci = FALSE, standardized = TRUE)#est为未标准化的，std.all为标准化后的
subset(Est, lhs == "f1") #anova(fit,othermodel) 这里如果有其他模型，可以写入相比较
fitMeasures(fit, c("cfi", "rmsea"))#返回检验模型的那些值
#> cfi rmsea
#> 0.931 0.059
#inspect(fit)
#fitted(fit)
# coef(fit)把结构图画出来
semPaths(fit,whatLabels="std.all",style="lisrel",
residuals=TRUE)
如果模型很长，可以将模型(单引号之间的内容)储存入myModel.lav的txt文档中，读取
#myModel <- readLines("/mydirectory/myModel.lav")其他，写在后面
#对于一个对应4个指标的潜变量，lavaan默认将第一个指标的因子载荷固定为1，其他指标为自由
#但如果有一个很好的理由来让所有因子载荷都固定为1，可以照如下做法：
#f =~ y1 + 1*y2 + 1*y3 + 1*y4
#如果想将某对变量的相关系数设为0，需要为其添加一个协方差公式并将参数设为0.在下面的语法中，wh9，max_edu，和care_ave 设为0。
#若希望固定f1的因子为一个单位，所以不需要设定第一个指标wz301为1，因此我们给wz301乘以NA使wz301的因子载荷自由且未知，整个模型如下：
model <- '
# latent variables
f1 =~ NA*wz301+ care_ave + wk6
# regressions
wh9 ~ f1 + max_edu + wf701 +finc_per
f1 ~ max_edu + wf701 + finc_per
wf701 ~ max_edu + finc_per
# residual covariances
#wh9 ~~ 0*care_ave
#max_edu ~~ 0*care_ave
#截距项
'
fit <- sem(model,data=city_data)
#这时出来的结果图就不会有虚线，
semPaths(fit,whatLabels="std.all",style="lisrel",
residuals=TRUE)
有些情况我们会预设变量参数相等，可以用相同的label达到效果，也可以用equal()函数
model <- '
# latent variables
f1 =~ NA*wz301+ care_ave + equal("visual=~care_ave")*wk6
# regressions
wh9 ~ f1 + max_edu + wf701 +finc_per
f1 ~ max_edu + wf701 + finc_per
wf701 ~ max_edu + finc_per
# residual covariances
#wh9 ~~ 0*care_ave
#max_edu ~~ 0*care_ave
#截距项
'
fit <- sem(model,data=city_data)
semPaths(fit,whatLabels="std.all",style="lisrel",
residuals=TRUE) ## 若对回归系数有所约束
#本来结果为
set.seed(1234)
Data <- data.frame(y = rnorm(100),
x1 = rnorm(100),
x2 = rnorm(100),
x3 = rnorm(100))
model <- '
y ~ b1*x1 + b2*x2 + b3*x3'
fit <- sem(model, data = Data)
coef(fit)
#> b1 b2 b3 y~~y
#> -0.052 0.084 0.139 0.970
#想要为b1增加两个约束，
model.constr <- '# model with labeled parameters
y ~ b1*x1 + b2*x2 + b3*x3
# constraints
b1 == (b2 + b3)^2
b1 > exp(b2 + b3)'
fit <- sem(model.constr, data = Data)
coef(fit)
#> b1 b2 b3 y~~y
#> 0.495 -0.405 -0.299 1.610- 对截距有所约束
three-factor model
y1 =~ x1 + x2 + x3 y2 =~ x4 + x5 + x6 y3 =~ x7 + x8 + x9 ### intercepts with fixed values x1 + x2 + x3 + x4 ~ 0.5*1
另外还有分组，进行多组分析，默认情况下，所有组会拟合相同的模型
潜变量模型的另一个重要种类是潜增长曲线模型。增长模型常用于分析纵向和发展数据
默认为列表状态删除，在进行统计量的计算时，把含有缺失值的记录删除
可以通过在summary()函数中增加参数modindices = TRUE，来得出修正指标