knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
hdbm is a Bayesian inference method that uses continuous shrinkage priors for
high-dimensional mediation analysis, developed by Song et al (2018).
hdbm provides estimates for the regression coefficients as well as
the posterior inclusion probability for ranking mediators.
You can install
hdbm from CRAN
or from github via
# install.packages(devtools) devtools::install_github("umich-cphds/hdbm", built_opts = c())
hdbm requires the R packages
RcppArmadillo, so you may want to
install / update them before downloading. If you decide to install
source (eg github), you will need a C++ compiler that supports C++11. On Windows
this can accomplished by installing
Xcode on MacOS.
hdbm contains a semi-synthetic example data set,
hdbm.data that is used in
hdbm.data contains a continuous response
y and a continuous
a that is mediated by 100 mediators,
library(hdbm) # print just the first 10 columns head(hdbm.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.
hdbm does not model internal correlation between mediators.
hdbm 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.
Y <- hdbm.data$y A <- hdbm.data$a # grab the mediators from the example data.frame M <- as.matrix(hdbm.data[, paste0("m", 1:100)], nrow(hdbm.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. hdbm.out <- hdbm(Y, A, M, C, C, beta.m, alpha.a, burnin = 1000, ndraws = 100) # Which mediators are active? active <- which(colSums(hdbm.out$r1 * hdbm.out$r3) > 100 / 2) colnames(M)[active]
Here, we calculate the posterior inclusion probability
r1 = r3 = 1 | Data,
and classify a mediator as active if its posterior probability is greater than
Yanyi Song, Xiang Zhou et al. Bayesian Shrinkage Estimation of High Dimensional Causal Mediation Effects in Omics Studies. bioRxiv 10.1101/467399
If you would like to report a bug, ask questions, or suggest something, please
e-mail Alexander Rix at
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.