Introduction to BayesGmed
In BayesGmed: Bayesian Causal Mediation Analysis using 'Stan'

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Motivation

This vignette introduces the R package BayesGmed, a package designed for a Bayesian causal mediation analysis in Stan, a probabilistic programming language for Bayesian inference. BayesGmed uses a parametric approach for the specification of the outcome and mediator model, and uses the g-formula approach for estimation. In addition to the estimation of causal mediation effects, the package also allows researchers to conduct sensitivity analysis.

Getting started

You can install from Github via devtools

devtools::install_gitlab(belayb/BayesGmed)

Example

BayesGmed comes with an example data, example_data, which contains a binary outcome $Y$, a single binary mediator $M$, a binary exposure $A$, and two numeric confounding variables $Z = (Z_1, Z_2)$. The first 6 rows of the data is visualized below.

library(BayesGmed)

head(example_data)

The aim in this example data is to estimate the direct and indirect effect of $A$ on $Y$ adjusting for $Z$. To do this, we may proced as follow.

[logit( P(Y_{i} = 1|A_{i},M_{i},\mathbf{Z}{i}) ) = \alpha{0} + \mathbf{\alpha}{Z}^{'}\mathbf{Z}{i} + \alpha_{A}A_{i} + \alpha_{M}M_{i},]

[logit(P(M_{i} = 1| A_{i},\mathbf{Z}{i} ) ) = \beta{0} + \mathbf{\beta}{Z}^{'}\mathbf{Z}{i} + \beta_{A}A_{i}.]

To complete the model specification, we assume the following priors for the model parameters.

$$ \begin{aligned} \beta &\sim MVN(location_m, scale_m)\ \alpha &\sim MVN(location_y, scale_y) \end{aligned} $$

We then need to specify the parameters of the prior distrbutions or assume a hyper-prior. BayesGmed currently only allow fixed values and these values can be passed as part of the model call and if not the following defult values will be used

P <- 3 # number of covariates plus the intercept term 
priors <- list(scale_m = 2.5*diag(P+1), 
               scale_y = 2.5*diag(P+2),
               location_m = rep(0, P+1), 
               location_y = rep(0, P+2))

The model can then be fitted as follow

fit1 <- bayesgmed(outcome = "Y", mediator =  "M", treat = "A", covariates = c("Z1", "Z2"), dist.y = "binary",
dist.m = "binary", link.y = "logit", link.m = "logit", data = example_data)

bayesgmed_summary(fit1)