## ---- 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("lavaan", quietly = TRUE) ||
    !requireNamespace("brms", quietly = TRUE) ||
    !requireNamespace("rstanarm", quietly = TRUE) ||
    !requireNamespace("insight", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
}

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
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")

## ----eval=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # 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)
#  

## ----echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
m2 <- insight::download_model("brms_mv_6")

## ----eval=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  # 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
#  )

## ----echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
m3 <- insight::download_model("stanmvreg_2")

## ---- message=TRUE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# for brms
mediation(m2)

# for rstanarm
mediation(m3)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
summary(m1)

mediation(m2, ci = .95)

mediation(m3, ci = .95)

## ---- message=TRUE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
m <- mediation(m2, centrality = "mean", ci = .95)
print(m, digits = 4)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(lavaan)
data(jobs)
set.seed(1234)

model <- ' # direct effects
             depress2 ~ c1*treat + c2*econ_hard + c3*sex + c4*age + b*job_seek
             
           # mediation
             job_seek ~ a1*treat + a2*econ_hard + a3*sex + a4*age
             
           # indirect effects (a*b)
             indirect_treat := a1*b
             indirect_econ_hard := a2*b
             indirect_sex := a3*b
             indirect_age := a4*b
             
           # total effects
             total_treat := c1 + (a1*b)
             total_econ_hard := c2 + (a2*b)
             total_sex := c3 + (a3*b)
             total_age := c4 + (a4*b)
         '
m4 <- sem(model, data = jobs)
summary(m4)

# just to have the numbers right at hand and you don't need to scroll up
mediation(m2, ci = .95)