Skip to contents

In this example, the Monte Carlo method is used to generate confidence intervals for the indirect effects in a serial mediation model with two mediators where X is the predictor, M1 is the first mediator, M2 is the second mediator, and Y is the dependent variable.

The Serial Mediation Model

Data

n <- 1000
X <- rnorm(n = n)
M1 <- 0.50 * X + rnorm(n = n)
M2 <- 0.15 * X + 0.50 * M1 + rnorm(n = n)
Y <- 0.10 * X + 0.05 * M1 + 0.50 * M2 + rnorm(n = n)
df <- data.frame(X, M1, M2, Y)

Model Specification

We can define several indirect effects in this example:

  • \(X \rightarrow M1 \rightarrow M2 \rightarrow Y\)
  • \(X \rightarrow M1 \rightarrow M2\)
  • \(X \rightarrow M1 \rightarrow Y\)
  • \(M1 \rightarrow M2 \rightarrow Y\)

These indirect effects are defined using the := operator in the lavaan model syntax.

model <- "
  Y ~ cp * X + b1 * M1 + b2 * M2
  M2 ~ a2 * X + k * M1
  M1 ~ a1 * X
  # X -> M1 -> M2 -> Y
  a1kb2 := a1 * k * b2
  # X -> M1 -> M2
  a1k := a1 * k
  # X -> M1 -> Y
  a2b2 := a2 * b2
  # M1 -> M2 -> Y
  kb2 := k * b2
"

Model Fitting

fit <- sem(data = df, model = model)

Monte Carlo Confidence Intervals

The fit lavaan object can then be passed to the MC() function from semmcci to generate Monte Carlo confidence intervals.

MC(fit, R = 20000L, alpha = c(0.001, 0.01, 0.05))
#> Monte Carlo Confidence Intervals
#>           est     se     R   0.05%    0.5%   2.5%  97.5%  99.5% 99.95%
#> cp     0.0868 0.0355 20000 -0.0298 -0.0053 0.0169 0.1554 0.1778 0.2017
#> b1     0.1190 0.0350 20000  0.0087  0.0290 0.0513 0.1880 0.2095 0.2301
#> b2     0.4809 0.0305 20000  0.3807  0.4022 0.4215 0.5404 0.5587 0.5832
#> a2     0.1373 0.0364 20000  0.0152  0.0447 0.0660 0.2092 0.2325 0.2580
#> k      0.4817 0.0331 20000  0.3739  0.3974 0.4166 0.5469 0.5672 0.5883
#> a1     0.5098 0.0310 20000  0.4082  0.4295 0.4483 0.5700 0.5895 0.6108
#> Y~~Y   0.9744 0.0439 20000  0.8304  0.8602 0.8884 1.0602 1.0882 1.1224
#> M2~~M2 1.0581 0.0477 20000  0.9039  0.9366 0.9653 1.1522 1.1817 1.2121
#> M1~~M1 0.9712 0.0431 20000  0.8311  0.8600 0.8862 1.0553 1.0815 1.1052
#> a1kb2  0.1181 0.0132 20000  0.0795  0.0860 0.0931 0.1452 0.1543 0.1631
#> a1k    0.2456 0.0226 20000  0.1773  0.1894 0.2022 0.2908 0.3057 0.3202
#> a2b2   0.0660 0.0180 20000  0.0068  0.0210 0.0315 0.1022 0.1137 0.1269
#> kb2    0.2317 0.0217 20000  0.1654  0.1789 0.1902 0.2751 0.2897 0.3110

Standardized Monte Carlo Confidence Intervals

Standardized Monte Carlo Confidence intervals can be generated by passing the result of the MC() function to the MCStd() function.

Note: We recommend setting fixed.x = FALSE when generating standardized estimates and confidence intervals to model the variances and covariances of the exogenous observed variables if they are assumed to be random. If fixed.x = TRUE, which is the default setting in lavaan, MC() will fix the variances and the covariances of the exogenous observed variables to the sample values.

fit <- sem(data = df, model = model, fixed.x = FALSE)
unstd <- MC(fit, R = 20000L, alpha = c(0.001, 0.01, 0.05))
MCStd(unstd)
#> Standardized Monte Carlo Confidence Intervals
#>           est     se     R   0.05%    0.5%   2.5%  97.5%  99.5% 99.95%
#> cp     0.0723 0.0291 20000 -0.0229 -0.0036 0.0146 0.1291 0.1456 0.1640
#> b1     0.1098 0.0318 20000  0.0075  0.0299 0.0485 0.1721 0.1915 0.2203
#> b2     0.4779 0.0273 20000  0.3866  0.4051 0.4232 0.5302 0.5475 0.5661
#> a2     0.1151 0.0307 20000  0.0108  0.0336 0.0556 0.1752 0.1935 0.2162
#> k      0.4470 0.0282 20000  0.3510  0.3711 0.3911 0.5014 0.5179 0.5364
#> a1     0.4602 0.0249 20000  0.3747  0.3930 0.4097 0.5081 0.5220 0.5407
#> Y~~Y   0.6725 0.0242 20000  0.5922  0.6084 0.6230 0.7181 0.7313 0.7433
#> M2~~M2 0.7396 0.0238 20000  0.6565  0.6759 0.6911 0.7846 0.7978 0.8134
#> M1~~M1 0.7882 0.0229 20000  0.7076  0.7275 0.7419 0.8322 0.8456 0.8596
#> X~~X   1.0000 0.0000 20000  1.0000  1.0000 1.0000 1.0000 1.0000 1.0000
#> a1kb2  0.0983 0.0103 20000  0.0677  0.0739 0.0787 0.1191 0.1262 0.1344
#> a1k    0.2057 0.0175 20000  0.1527  0.1621 0.1717 0.2409 0.2515 0.2649
#> a2b2   0.0550 0.0150 20000  0.0048  0.0160 0.0262 0.0848 0.0950 0.1036
#> kb2    0.2136 0.0186 20000  0.1564  0.1666 0.1778 0.2503 0.2635 0.2755