---
title: "Summary of Mediation Analysis using Bayesian Regression Models"
author: "Daniel Lüdecke"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
---

<!--
vignette: >
  %\VignetteIndexEntry{Summary of Mediation Analysis using Bayesian Regression Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
-->

```{r, SETTINGS-knitr, echo = FALSE, warning = FALSE, message = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, 
  comment = "#>", 
  dev = "png", 
  fig.width = 7, 
  fig.height = 5, 
  message = FALSE, 
  warning = FALSE
)

options(width = 800)

if (!requireNamespace("mediation", quietly = TRUE) ||
    !requireNamespace("httr", quietly = TRUE) ||
    !requireNamespace("brms", quietly = TRUE) ||
    !requireNamespace("insight", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
}
```

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_.

```{r}
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")
```
```{r eval=FALSE}
# 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)
```
```{r echo=FALSE}
m2 <- insight::download_model("brms_mv_6")
```

`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.

```{r, message=TRUE}
mediation(m2)
```

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.

```{r}
summary(m1)

mediation(m2, prob = .95)
```

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.

```{r, message=TRUE}
mediation(m2, typical = "mean", prob = .95) %>% print(digits = 4)
```

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

# References

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