---
title: "High Dimensional Bayesian Mediation Analysis in R"
author: "Alexander Rix"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{High Dimensional Bayesian Mediation Analysis in R}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```
# bama

`bama` is a Bayesian inference method that uses continuous shrinkage priors for
high-dimensional mediation analysis, developed by Song et al (2018).
`bama` provides estimates for the regression coefficients as well as
the posterior inclusion probability for ranking mediators.

## Installation
You can install `bama` from CRAN
```{r, eval = FALSE}
install.packages("bama")

```

or from github via `devtools`
```{r, eval = FALSE}
# install.packages(devtools)
devtools::install_github("umich-cphds/bama", built_opts = c())
```
`bama` requires the R packages `Rcpp` and `RcppArmadillo`, so you may want to
install / update them before downloading. If you decide to install `bama` from
source (eg github), you will need a C++ compiler that supports C++11. On Windows
this can accomplished by installing
[Rtools](https://cran.r-project.org/bin/windows/Rtools/), and
[Xcode](https://apps.apple.com/us/app/xcode/id497799835?mt=12) on MacOS.

## Example problem
`bama` contains a semi-synthetic example data set, `bama.data` that is used in
this example. `bama.data` contains a continuous response `y` and a continuous
exposure `a` that is mediated by 100 mediators, `m[1:100]`.

```{r}
library(bama)
# print just the first 10 columns
head(bama.data[,1:10])
```

The mediators have an internal correlation structure that is based off the
covariance matrix from the Multi-Ethnic Study of Atherosclerosis (MESA) data.
However, `bama` does not model internal correlation between mediators.
Instead, `bama` employs continuous Bayesian shrinkage priors to select mediators
and assumes that all the potential mediators contribute small effects
in mediating the exposure-outcome relationship, but only a small proportion of
mediators exhibit large effects.

We use no adjustment covariates in this example, so we just include the
intercept. Also, in a real world situation, it may be beneficial to normalize
the input data.

```{r}

Y <- bama.data$y
A <- bama.data$a

# grab the mediators from the example data.frame
M <- as.matrix(bama.data[, paste0("m", 1:100)], nrow(bama.data))

# We just include the intercept term in this example.
C <- matrix(1, nrow(M), 1)

# Initial guesses for coefficients
beta.m  <- rep(0, ncol(M))
alpha.a <- rep(0, ncol(M))

set.seed(12345)
# It is recommended to pick a larger number for burnin.
bama.out <- bama(Y, A, M, C, beta.m, alpha.a, burnin = 3000, ndraws = 100)

# rank the mediators by PIP and return the highest 10
summary(bama.out, rank = T)[1:10,]
```

Here, the summary function calculates the posterior inclusion probability (PIP)
`r1 = r3 = 1 | Data`, and rank each mediator by its PIP. Each mediator's
estimate and 95\% credible interval is also calculated.

## Reference
Yanyi Song, Xiang Zhou et al. Bayesian Shrinkage Estimation of
High Dimensional Causal Mediation Effects in Omics Studies.
bioRxiv [10.1101/467399](https://doi.org/10.1101/467399)

## Contact
If you would like to report a bug, ask questions, or suggest something, please
e-mail Alexander Rix at `alexrix@umich.edu`.
