# Summary of Mediation Analysis using Bayesian Regression Models

#### Daniel Lüdecke

#### 2020-02-06

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.

# References

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