In this example, the Monte Carlo method is used to generate confidence intervals for the indirect effect in a simple mediation model where variable X has an effect on variable Y, through a mediating variable M. This mediating or indirect effect is a product of path coefficients from the fitted model. The Simple Mediation Model

library(semmcci)
library(lavaan)

## Data

n <- 1000
a <- 0.50
b <- 0.50
cp <- 0.25
s2_em <- 1 - a^2
s2_ey <- 1 - cp^2 - a^2 * b^2 - b^2 * s2_em - 2 * cp * a * b
em <- rnorm(n = n, mean = 0, sd = sqrt(s2_em))
ey <- rnorm(n = n, mean = 0, sd = sqrt(s2_ey))
X <- rnorm(n = n)
M <- a * X + em
Y <- cp * X + b * M + ey
df <- data.frame(X, M, Y)

## Model Specification

The indirect effect is defined by the product of the slopes of paths X to M labeled as a and M to Y labeled as b. In this example, we are interested in the confidence intervals of indirect defined as the product of a and b using the := operator in the lavaan model syntax.

model <- "
Y ~ cp * X + b * M
M ~ a * X
indirect := a * b
direct := cp
total := cp + (a * b)
"

## Model Fitting

We can now fit the model using the sem() function from lavaan.

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

## Monte Carlo Confidence Intervals

The fit lavaan object can then be passed to the MC() function 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.2333 0.0264 20000 0.1465 0.1651 0.1817 0.2846 0.3012 0.3165
#> b        0.5082 0.0272 20000 0.4161 0.4385 0.4551 0.5611 0.5786 0.5988
#> a        0.4820 0.0264 20000 0.3959 0.4140 0.4299 0.5337 0.5495 0.5642
#> Y~~Y     0.5462 0.0244 20000 0.4660 0.4832 0.4979 0.5944 0.6095 0.6263
#> M~~M     0.7527 0.0339 20000 0.6441 0.6655 0.6858 0.8187 0.8387 0.8620
#> indirect 0.2449 0.0188 20000 0.1885 0.1983 0.2094 0.2831 0.2964 0.3105
#> direct   0.2333 0.0264 20000 0.1465 0.1651 0.1817 0.2846 0.3012 0.3165
#> total    0.4782 0.0264 20000 0.3911 0.4102 0.4262 0.5295 0.5457 0.5629

## 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.2422 0.0268 20000 0.1535 0.1722 0.1893 0.2945 0.3101 0.3294
#> b        0.5123 0.0246 20000 0.4289 0.4494 0.4639 0.5607 0.5764 0.5937
#> a        0.4963 0.0239 20000 0.4157 0.4338 0.4488 0.5427 0.5575 0.5754
#> Y~~Y     0.5558 0.0235 20000 0.4805 0.4951 0.5095 0.6018 0.6164 0.6334
#> M~~M     0.7537 0.0237 20000 0.6690 0.6892 0.7055 0.7986 0.8118 0.8272
#> X~~X     1.0000 0.0000 20000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
#> indirect 0.2542 0.0176 20000 0.1974 0.2101 0.2206 0.2892 0.3005 0.3154
#> direct   0.2422 0.0268 20000 0.1535 0.1722 0.1893 0.2945 0.3101 0.3294
#> total    0.4964 0.0239 20000 0.4159 0.4342 0.4485 0.5418 0.5564 0.5712

## Standardized Monte Carlo Confidence Intervals - An Alternative Approach

In this example, confidence intervals for the standardized indirect effect are generated by specifying the standardized indirect effect as a derived parameter using the := operator. The standardized indirect effect in a simple mediation model involves paths $$a$$ and $$b$$, and the standard deviations of $$X$$ and $$Y$$. It is given by

$\begin{equation} a b \frac{s_X}{s_Y} \end{equation}$

where

$\begin{equation} s_X = \sqrt{ s_{X}^{2} } \end{equation}$

and

$\begin{equation} s_Y = \sqrt{ c^{\prime 2} s_{X}^{2} + a^2 b^2 s_{X}^{2} + b^2 s_{e_{M}}^{2} + 2 c^{\prime} b a s_{X}^{2} + s_{e_{Y}}^{2} } \end{equation}$

where

$$s_{e_{Y}}^{2}$$ and $$s_{e_{M}}^{2}$$ are the residual variances in the regression equations.

The standardized indirect effect can be defined using the := operator and the named parameters in the model.

model <- "
Y ~ cp * X + b * M
M ~ a * X
X ~~ s2_X * X
M ~~ s2_em * M
Y ~~ s2_ey * Y
indirect_std := a * b * (
sqrt(s2_X) / sqrt(
(
cp^2 * s2_X + a^2 * b^2 * s2_X
) + (
b^2 * s2_em
) + (
2 * cp * b * a * s2_X
) + s2_ey
)
)
"
fit <- sem(data = df, model = model, fixed.x = FALSE)

The row indirect_std corresponds to the confidence intervals for the standardized indirect effect.

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.2333 0.0263 20000 0.1480 0.1648 0.1810 0.2838 0.3016 0.3204
#> b            0.5082 0.0270 20000 0.4154 0.4387 0.4553 0.5613 0.5769 0.5920
#> a            0.4820 0.0264 20000 0.3974 0.4155 0.4309 0.5338 0.5499 0.5733
#> s2_X         1.0590 0.0472 20000 0.9108 0.9390 0.9675 1.1531 1.1797 1.2133
#> s2_em        0.7527 0.0339 20000 0.6404 0.6655 0.6865 0.8192 0.8401 0.8665
#> s2_ey        0.5462 0.0243 20000 0.4661 0.4840 0.4987 0.5940 0.6086 0.6258
#> indirect_std 0.2542 0.0174 20000 0.2011 0.2104 0.2204 0.2884 0.3000 0.3128