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 devtools
hdbm
requires the R packages Rcpp
and
RcppArmadillo
, so you may want to install / update them
before downloading. If you decide to install hdbm
from
source (eg github), you will need a C++ compiler that supports C++11. On
Windows this can accomplished by installing Rtools, and Xcode
on MacOS.
hdbm
contains a semi-synthetic example data set,
hdbm.data
that is used in this example.
hdbm.data
contains a continuous response y
and
a continuous exposure a
that is mediated by 100 mediators,
m[1:100]
.
library(hdbm)
# print just the first 10 columns
head(hdbm.data[,1:10])
#> y a m1 m2 m3 m4
#> 1 -0.5077701 1.3979467 -0.75395346 -0.2787043 -0.04471833 0.3422936
#> 2 -0.3239898 -0.2311032 -1.20208195 0.4210638 0.93175992 -0.3699733
#> 3 -1.8553536 -2.4647028 -0.65712133 0.3285993 0.59144748 -0.1307554
#> 4 0.1685455 0.1119932 0.04982723 0.6816996 -0.12956715 -0.8348541
#> 5 0.9070900 0.4994626 -0.99964057 -0.7660710 1.53962908 -0.9308951
#> 6 -1.0357105 0.6359685 1.06954128 -0.2441489 -1.52176072 0.4657214
#> m5 m6 m7 m8
#> 1 -0.66227113 -0.30925865 1.58001664 0.008326522
#> 2 1.09811497 -0.09969085 1.02369272 0.045104531
#> 3 -0.30196963 0.38853526 -0.05841533 -0.436429826
#> 4 0.08936191 -0.69699157 0.41615473 0.973411472
#> 5 1.12107670 1.07603088 0.37449777 -0.289794580
#> 6 -1.55992443 -0.42705075 -0.98761802 -0.639473238
The mediators have an internal correlation structure that is based
off the covariance matrix from the Multi-Ethnic Study of Atherosclerosis
(MESA) data. However, hdbm
does not model internal
correlation between mediators. Instead, 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]
#> [1] "m12" "m65" "m89"
Here, we calculate the posterior inclusion probability
r1 = r3 = 1 | Data
, and classify a mediator as active if
its posterior probability is greater than 0.5.
Yanyi Song, Xiang Zhou et al. Bayesian Shrinkage Estimation of High Dimensional Causal Mediation Effects in Omics Studies. bioRxiv 10.1101/467399
##Contact information If you would like to report a bug, ask
questions, or suggest something, please e-mail Alexander Rix at
[email protected]
.