`vignettes/mediation.Rmd`

`mediation.Rmd`

This vignettes demonstrates the `mediation()`

-function. Before we start, we fit some models, including a mediation-object from the *mediation*-package, which we use for comparison with *brms* and *rstanarm*.

library(bayestestR) library(mediation) library(brms) library(rstanarm) # load sample data data(jobs) set.seed(123) # linear models, for mediation analysis b1 <- lm(job_seek ~ treat + econ_hard + sex + age, data = jobs) b2 <- lm(depress2 ~ treat + job_seek + econ_hard + sex + age, data = jobs) # mediation analysis, for comparison with brms m1 <- mediate(b1, b2, sims = 1000, treat = "treat", mediator = "job_seek")

# Fit Bayesian mediation model in brms f1 <- bf(job_seek ~ treat + econ_hard + sex + age) f2 <- bf(depress2 ~ treat + job_seek + econ_hard + sex + age) m2 <- brm(f1 + f2 + set_rescor(FALSE), data = jobs, cores = 4)

# Fit Bayesian mediation model in rstanarm m3 <- stan_mvmer( list(job_seek ~ treat + econ_hard + sex + age + (1 | occp), depress2 ~ treat + job_seek + econ_hard + sex + age + (1 | occp)), data = jobs, cores = 4, refresh = 0 )

`mediation()`

is a summary function, especially for mediation analysis, i.e. for multivariate response models with casual mediation effects.

In the models `m2`

and `m3`

, `treat`

is the treatment effect and `job_seek`

is the mediator effect. For the *brms* model (`m2`

), `f1`

describes the mediator model and `f2`

describes the outcome model. This is similar for the *rstanarm* model.

`mediation()`

returns a data frame with information on the *direct effect* (median value of posterior samples from treatment of the outcome model), *mediator effect* (median value of posterior samples from mediator of the outcome model), *indirect effect* (median value of the multiplication of the posterior samples from mediator of the outcome model and the posterior samples from treatment of the mediation model) and the *total effect* (median value of sums of posterior samples used for the direct and indirect effect). The *proportion mediated* is the indirect effect divided by the total effect.

The simplest call just needs the model-object.

# for brms mediation(m2) #> # Causal Mediation Analysis for Stan Model #> #> Treatment: treat #> Mediator : job_seek #> Response : depress2 #> #> Effect | Estimate | 89% ETI #> ---------------------------------------------------- #> Direct Effect (ADE) | -0.040 | [-0.110, 0.031] #> Indirect Effect (ACME) | -0.015 | [-0.036, 0.004] #> Mediator Effect | -0.240 | [-0.285, -0.195] #> Total Effect | -0.055 | [-0.129, 0.018] #> #> Proportion mediated: 28.14% [-71.11%, 127.40%] # for rstanarm mediation(m3) #> # Causal Mediation Analysis for Stan Model #> #> Treatment: treat #> Mediator : job_seek #> Response : depress2 #> #> Effect | Estimate | 89% ETI #> ---------------------------------------------------- #> Direct Effect (ADE) | -0.040 | [-0.111, 0.031] #> Indirect Effect (ACME) | -0.018 | [-0.037, 0.002] #> Mediator Effect | -0.241 | [-0.286, -0.197] #> Total Effect | -0.057 | [-0.130, 0.017] #> #> Proportion mediated: 30.59% [-75.65%, 136.82%]

Typically, `mediation()`

finds the treatment and mediator variables automatically. If this does not work, use the `treatment`

and `mediator`

arguments to specify the related variable names. For all values, the 89% credible intervals are calculated by default. Use `ci`

to calculate a different interval.

Here is a comparison with the *mediation* package. Note that the `summary()`

-output of the *mediation* package shows the indirect effect first, followed by the direct effect.

summary(m1) #> #> Causal Mediation Analysis #> #> Quasi-Bayesian Confidence Intervals #> #> Estimate 95% CI Lower 95% CI Upper p-value #> ACME -0.0157 -0.0387 0.01 0.19 #> ADE -0.0438 -0.1315 0.04 0.35 #> Total Effect -0.0595 -0.1530 0.02 0.21 #> Prop. Mediated 0.2137 -2.0277 2.70 0.32 #> #> Sample Size Used: 899 #> #> #> Simulations: 1000 mediation(m2, ci = .95) #> # Causal Mediation Analysis for Stan Model #> #> Treatment: treat #> Mediator : job_seek #> Response : depress2 #> #> Effect | Estimate | 95% ETI #> ---------------------------------------------------- #> Direct Effect (ADE) | -0.040 | [-0.124, 0.046] #> Indirect Effect (ACME) | -0.015 | [-0.041, 0.008] #> Mediator Effect | -0.240 | [-0.294, -0.185] #> Total Effect | -0.055 | [-0.145, 0.034] #> #> Proportion mediated: 28.14% [-181.46%, 237.75%] mediation(m3, ci = .95) #> # Causal Mediation Analysis for Stan Model #> #> Treatment: treat #> Mediator : job_seek #> Response : depress2 #> #> Effect | Estimate | 95% ETI #> ---------------------------------------------------- #> Direct Effect (ADE) | -0.040 | [-0.129, 0.048] #> Indirect Effect (ACME) | -0.018 | [-0.042, 0.006] #> Mediator Effect | -0.241 | [-0.296, -0.187] #> Total Effect | -0.057 | [-0.151, 0.033] #> #> Proportion mediated: 30.59% [-221.09%, 282.26%]

If you want to calculate mean instead of median values from the posterior samples, use the `centrality`

-argument. Furthermore, there is a `print()`

-method, which allows to print more digits.

m <- mediation(m2, centrality = "mean", ci = .95) print(m, digits = 4) #> # Causal Mediation Analysis for Stan Model #> #> Treatment: treat #> Mediator : job_seek #> Response : depress2 #> #> Effect | Estimate | 95% ETI #> ------------------------------------------------------ #> Direct Effect (ADE) | -0.0395 | [-0.1237, 0.0456] #> Indirect Effect (ACME) | -0.0158 | [-0.0405, 0.0083] #> Mediator Effect | -0.2401 | [-0.2944, -0.1846] #> Total Effect | -0.0553 | [-0.1454, 0.0341] #> #> Proportion mediated: 28.60% [-181.01%, 238.20%]

As you can see, the results are similar to what the *mediation* package produces for non-Bayesian models.