Chong Xing, Center for Research Methods and Data Analysis, University of Kansas <cxing@ku.edu>
Please visit http://crmda.ku.edu/guides
Keywords: Linear Regression, Path Analysis, Mediation, R, lavaan
2019 January 24
Abstract
This guide outlines how fit two path models in R using the lavaan package. The first example shows how to specify and estimate an indirect effect (or mediation) model using lavaan with bootstrapped standard errors. The second example demonstrates how to fit the same model for two groups.Load the lavaan package at the beginning of the session
library(lavaan)
This is lavaan 0.6-3
lavaan is BETA software! Please report any bugs.
The data file is read in, columns are named, and missing values are specified.
dat <- read.csv("../../data/job_placement.csv", header = FALSE)
colnames(dat) <- c("id", "wjcalc", "wjspl", "wratspl", "wratcalc",
"waiscalc", "waisspl", "edlevel", "newschl",
"suspend", "expelled", "haveld", "gender", "age")
dat[dat == 99999] <- NA
dat$gender <- as.factor(dat$gender)
levels(dat$gender) <- c("Male", "Female")
This model contains one outcome variable “wjspl”, one predictor variable “age”, and one mediator variable “edlevel”. The hypothesis being tested is that the influence of age on wjspl can be explained (or partially explained) by edlevel.
model.path <-
' wjspl ~ c*age ## the c path - the direct path from IV to DV
wjspl ~ b*edlevel ## the b path - the direct path from mediator to DV
edlevel ~ a*age ## the a path - the direct path from IV to mediator
ab := a*b ## using := operator to create ab (the indirect effect)
total := c + a*b ## using := operator to create total effect
'
This fits the model using the sem function and returns an output summary. Bootstrapped standard errors (1000 bootstrapped draws) are requested.
output.path <- sem(model = model.path, data = dat, se = "bootstrap",
bootstrap = 1000)
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in lavaan(slotOptions = lavoptions, slotParTable = lavpartable, :
lavaan WARNING: the optimizer warns that a solution has NOT been found!
Warning in bootstrap.internal(object = NULL, lavmodel. = lavmodel,
lavsamplestats. = lavsamplestats, : lavaan WARNING: only 990 bootstrap
draws were successful
summary(output.path)
lavaan 0.6-3 ended normally after 24 iterations
Optimization method NLMINB
Number of free parameters 5
Used Total
Number of observations 320 322
Estimator ML
Model Fit Test Statistic 0.000
Degrees of freedom 0
Minimum Function Value 0.0000000000000
Parameter Estimates:
Standard Errors Bootstrap
Number of requested bootstrap draws 1000
Number of successful bootstrap draws 990
Regressions:
Estimate Std.Err z-value P(>|z|)
wjspl ~
age (c) 0.440 0.243 1.811 0.070
edlevel (b) 1.226 0.364 3.367 0.001
edlevel ~
age (a) 0.269 0.037 7.366 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.wjspl 46.890 3.502 13.389 0.000
.edlevel 1.360 0.104 13.130 0.000
Defined Parameters:
Estimate Std.Err z-value P(>|z|)
ab 0.330 0.106 3.113 0.002
total 0.770 0.218 3.522 0.000
The same path model can be estimated for multiple groups (e.g., gender).
model.path.2 <-
' wjspl ~ c(b1, b2)*edlevel + c(c1, c2)*age ## separate b and c paths for two groups
edlevel ~ c(a1, a2)*age ## separate a paths for two groups
a1b1 := a1*b1 ## the indirect effect for group 1
a2b2 := a2*b2 ## the indirect effect for group 2
total1 := c1 + a1*b1 ## the total effect for group 1
total2 := c2 + a2*b2 ## the total effect for group 2
'
This fits the model using the sem function and returns an output summary. Bootstrapped standard errors (1000 bootstrapped draws) are requested. The indirect and the total effects are estimated simultaneously for two groups.
output.path.2 <- sem(model = model.path.2, data = dat, se = "bootstrap",
bootstrap = 1000, group = "gender")
summary(output.path.2)
lavaan 0.6-3 ended normally after 60 iterations
Optimization method NLMINB
Number of free parameters 14
Number of observations per group
Female 101
Male 219
Estimator ML
Model Fit Test Statistic 0.000
Degrees of freedom 0
Minimum Function Value 0.0000000000000
Chi-square for each group:
Female 0.000
Male 0.000
Parameter Estimates:
Standard Errors Bootstrap
Number of requested bootstrap draws 1000
Number of successful bootstrap draws 1000
Group 1 [Female]:
Regressions:
Estimate Std.Err z-value P(>|z|)
wjspl ~
edlevel (b1) 1.810 0.569 3.184 0.001
age (c1) 0.471 0.402 1.172 0.241
edlevel ~
age (a1) 0.265 0.057 4.644 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
.wjspl 13.339 6.757 1.974 0.048
.edlevel 5.678 1.135 5.001 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.wjspl 41.686 5.149 8.096 0.000
.edlevel 1.426 0.171 8.355 0.000
Group 2 [Male]:
Regressions:
Estimate Std.Err z-value P(>|z|)
wjspl ~
edlevel (b2) 1.031 0.427 2.417 0.016
age (c2) 0.385 0.317 1.214 0.225
edlevel ~
age (a2) 0.272 0.045 6.087 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
.wjspl 22.293 5.386 4.139 0.000
.edlevel 5.852 0.893 6.551 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.wjspl 48.374 4.388 11.023 0.000
.edlevel 1.299 0.124 10.474 0.000
Defined Parameters:
Estimate Std.Err z-value P(>|z|)
a1b1 0.479 0.191 2.502 0.012
a2b2 0.280 0.120 2.336 0.020
total1 0.951 0.358 2.657 0.008
total2 0.665 0.283 2.350 0.019
R version 3.5.1 (2018-07-02)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lavaan_0.6-3 stationery_0.98.5.7
loaded via a namespace (and not attached):
[1] Rcpp_1.0.0 digest_0.6.18 plyr_1.8.4 xtable_1.8-3
[5] magrittr_1.5 stats4_3.5.1 evaluate_0.12 zip_1.0.0
[9] stringi_1.2.4 pbivnorm_0.6.0 openxlsx_4.1.0 rmarkdown_1.11
[13] tools_3.5.1 stringr_1.3.1 foreign_0.8-71 kutils_1.59
[17] yaml_2.2.0 xfun_0.4 compiler_3.5.1 mnormt_1.5-5
[21] htmltools_0.3.6 knitr_1.21
Available under Created Commons license 3.0