This vignettes demontrates the `mediation()`

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

```
library(sjstats)
library(mediation)
library(brms)
# 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
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)
```

`mediation()`

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

In the model *m2*, *treat* is the treatment effect, *job_seek* is the mediator effect, *f1* describes the mediator model and *f2* describes the outcome 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.

```
mediation(m2)
#>
#> # Causal Mediation Analysis for Stan Model
#>
#> Treatment: treat
#> Mediator: job_seek
#> Response: depress2
#>
#> Estimate HDI (90%)
#> Direct effect: -0.04 [-0.11 0.03]
#> Indirect effect: -0.02 [-0.04 0.00]
#> Total effect: -0.05 [-0.13 0.02]
#>
#> Proportion mediated: 28.14% [-79.57% 135.86%]
```

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 90% HDIs are calculated by default. Use `prob`

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, prob = .95)
#>
#> # Causal Mediation Analysis for Stan Model
#>
#> Treatment: treat
#> Mediator: job_seek
#> Response: depress2
#>
#> Estimate HDI (95%)
#> Direct effect: -0.04 [-0.12 0.04]
#> Indirect effect: -0.02 [-0.04 0.01]
#> Total effect: -0.05 [-0.15 0.03]
#>
#> Proportion mediated: 28.14% [-178.65% 234.94%]
```

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

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

-method, which allows to print more digits.

```
mediation(m2, typical = "mean", prob = .95) %>% print(digits = 4)
#>
#> # Causal Mediation Analysis for Stan Model
#>
#> Treatment: treat
#> Mediator: job_seek
#> Response: depress2
#>
#> Estimate HDI (95%)
#> Direct effect: -0.0395 [-0.1244 0.0450]
#> Indirect effect: -0.0158 [-0.0400 0.0086]
#> Total effect: -0.0553 [-0.1482 0.0302]
#>
#> Proportion mediated: 28.5975% [-178.1953% 235.3902%]
```

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

Bürkner, P. C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1-28