BKMRCMA Quick Start Example

Source: vignettes/BKMRCMA_QuickStart.Rmd

In this scenario, we have a continuous mediator MM, a continuous outcome YY, and x2 as the effect modifier on YY. The sample size is 50 and there are 3 covariates.

Sample simulated Dataset 

dat <-  cma_sampledata(N = 50, L=3, P=3, scenario=1, seed=7) 
head(dat$data, n = 3L)
#>           z1         z2         z3           M         x1          x2        x3
#> 1  2.2767428  0.1432685  0.7973699  0.58827130 -0.3556944  2.02334405 -1.477592
#> 2 -1.0893537 -0.0944783  0.4035083  0.01382552  1.0973004  0.86249250 -2.973256
#> 3 -0.6389109  0.2491581 -0.1243989 -0.42344188 -0.9066920 -0.02490949 -1.339762
#>           y
#> 1  1.964451
#> 2 -1.526532
#> 3 -1.305984
dat = dat$data

Fit BKMR Models

let AA be a nn-by-LL matrix containing an exposure mixture comprised of LL elements,E.M and E.Y be effect modifiers of exposure-mediator and exposure-outcome relationship respectively, yy, a vector of outcome data, and mm, a vector of mediator data.

Z.M <- cbind(A,E.M)

Z.Y <- cbind(A,E.Y)

Let Z.M be the exposures and effect modifers E.M and let Z.Y be the exposures and effect modifers E.Y, create one more matrix containing the exposures, effect modifier Z.Y and mediator, precisely in that order.

Zm.Y <- cbind(Z.Y,m)

NOTE: the last column of the Zm.Y matrix must be your mediator in order for the functions to work properly.

A <- cbind(dat$z1, dat$z2, dat$z3)
X <- cbind(dat$x1, dat$x2, dat$x3)
y  <- dat$y
m  <- dat$M 

E.M <- NULL
E.Y <- dat$x2

Z.M <- cbind(A,E.M) 
Z.Y <- cbind(A, E.Y) 
Zm.Y <- cbind(Z.Y, m)

set.seed(1)
fit.y <- kmbayes(y=y, Z=Zm.Y, X=X, iter=10000, verbose=TRUE, varsel=FALSE) 
#save(fit.y,file="bkmr_y.RData")

set.seed(2)
fit.y.TE <- kmbayes(y=y, Z=Z.Y, X=X, iter=10000, verbose=TRUE, varsel=FALSE) 
#save(fit.y.TE,file="bkmr_y_TE.RData")

set.seed(3)
fit.m <- kmbayes(y=m, Z=Z.M, X=X, iter=10000, verbose=TRUE, varsel=FALSE) 
#save(fit.m,file="bkmr_m.RData")

Values at which to predict counterfactuals

Mean level of confounders:

X.predict <- matrix(colMeans(X),nrow=1)

We define the change in exposure for which you want to estimate the mediation effects: in this example, we will consider a change in all exposures from their 25th to 75th percentiles, fixing age (E.Y) at testing to its 10th and 90th percentiles. However, this contrast can be anything.

Note: If modifiers are considered, you should fix the levels of the modifiers

astar <- c(apply(A, 2, quantile, probs=0.25))  
a <- c(apply(A, 2, quantile, probs=0.75))

e.y10 = quantile(E.Y, probs=0.1)
e.y90 = quantile(E.Y, probs=0.9)

The index of the MCMC iterations to be used for inference:

sel<-seq(5000,10000,by=10)

Estimate TE for BKMR

Estimate the TE for a change in the exposures from a*a^* to aa fixing Effect modifier at testing to its 10th percentile or 90th percentile:

e.y10 = quantile(E.Y, probs=0.1)
e.y90 = quantile(E.Y, probs=0.9)

TE.ey10 <- TE.bkmr(a=a, astar=astar, e.y = e.y10, fit.y.TE=fit.y.TE, X.predict=X.predict, alpha=0.05, sel=sel, seed=122)

TE.ey90 <- TE.bkmr(a=a, astar=astar, e.y = e.y90, fit.y.TE=fit.y.TE, X.predict=X.predict, alpha=0.05, sel=sel, seed=122)

Look at the posterior mean, median, and 95% CI for TE:

TE.ey10$est
#>         mean    median     lower     upper         sd
#> TE 0.4414075 0.4427607 0.2848167 0.5929055 0.07907911
TE.ey90$est
#>        mean   median   lower    upper         sd
#> TE 1.339614 1.338395 1.17107 1.494782 0.07832421

plotdf <- as.data.frame(TE.ey10$est)
plotdf["Effect"] <- rownames(plotdf)
ggplot(plotdf, aes(Effect, mean, ymin = lower, ymax = upper))  + 
  geom_pointrange(position = position_dodge(width = 0.75))  +  coord_flip()

Estimate CDE for BKMR

Estimate the CDE for a change in the exposures from a*a^* to aa, fixing the mediator at its 10th, 50th, and 75th percentile and the effect modifier at testing at its 10th percentile:

CDE.ey10 <- CDE.bkmr(a=a, astar=astar, e.y = e.y10, m.quant=c(0.1,.25, 0.5,0.75), fit.y=fit.y, alpha=0.05, sel=sel, seed=777)
#> [1] "Running 4 mediator values for CDE:"
#> [1] "1 out of 4"
#> [1] "2 out of 4"
#> [1] "3 out of 4"
#> [1] "4 out of 4"

CDE.ey90 <- CDE.bkmr(a=a, astar=astar, e.y = e.y90, m.quant=c(0.1,.25, 0.5,0.75), fit.y=fit.y, alpha=0.05, sel=sel, seed=777)
#> [1] "Running 4 mediator values for CDE:"
#> [1] "1 out of 4"
#> [1] "2 out of 4"
#> [1] "3 out of 4"
#> [1] "4 out of 4"

Look at the posterior mean, median, and 95% CI for the CDEs:

CDE.ey10$est
#>             mean    median     lower     upper         sd
#> CDE10% 0.3281490 0.3348591 0.1343525 0.5058129 0.09798879
#> CDE25% 0.3370905 0.3438450 0.1502369 0.5111669 0.09441366
#> CDE50% 0.3505605 0.3589534 0.1608827 0.5232460 0.09380879
#> CDE75% 0.3612455 0.3687880 0.1603481 0.5388214 0.09718429
CDE.ey90$est
#>            mean   median     lower    upper         sd
#> CDE10% 1.220641 1.224680 0.9837519 1.415078 0.10790548
#> CDE25% 1.229732 1.233787 1.0112267 1.413634 0.10027420
#> CDE50% 1.239655 1.244461 1.0452505 1.416089 0.09336095
#> CDE75% 1.244387 1.252115 1.0555353 1.419519 0.09252889

Plotting:

plotdf <- as.data.frame(CDE.ey10$est)
plotdf["Effect"] <- rownames(plotdf)
ggplot(plotdf, aes(Effect, mean, ymin = lower, ymax = upper ))  + 
  geom_pointrange(position = position_dodge(width = 0.5))  +  coord_flip()

Estimate NDE/NIE for BKMR

Estimate the TE, NDE and NIE for a change in the exposures from a*a^* to aa fixing age at testing to its 90th percentile.

Note: this step takes a while to run and will take longer for more complex BKMR fits, longer sel vectors and larger.

mediationeffects.ey10  <- mediation.bkmr(a=a, astar=astar, e.y = e.y10, fit.m=fit.m, fit.y=fit.y, fit.y.TE=fit.y.TE, X.predict.M=X.predict, X.predict.Y=X.predict, alpha=0.05, sel=sel, seed=22, K=10)

mediationeffects.ey90  <- mediation.bkmr(a=a, astar=astar, e.y = e.y90, fit.m=fit.m, fit.y=fit.y, fit.y.TE=fit.y.TE, X.predict.M=X.predict, X.predict.Y=X.predict, alpha=0.05, sel=sel, seed=22, K=10)

Look at the posterior mean, median, and 95% CI for the TE, NDE, and NIE

mediationeffects.ey10$est
#>               mean     median       lower     upper         sd
#> TE      0.43344519  0.4373605  0.28608048 0.5719413 0.07642777
#> NDE     0.46376164  0.4662025 -0.02268424 0.9612434 0.24502403
#> NIE    -0.03031645 -0.0360115 -0.54497247 0.4376174 0.24559061
#> CDE10%  0.31886211  0.3202047  0.12609272 0.5173728 0.09825198
#> CDE50%  0.34154761  0.3432965  0.15536080 0.5271944 0.09374263
#> CDE75%  0.35179029  0.3559092  0.14178624 0.5365690 0.09660131
mediationeffects.ey90$est
#>              mean     median      lower     upper         sd
#> TE     1.33206315 1.33310908  1.1860711 1.4739473 0.07412271
#> NDE    1.30327228 1.29575422  0.7502507 1.8822087 0.26921873
#> NIE    0.02879087 0.02299106 -0.5259615 0.5433754 0.26874946
#> CDE10% 1.21024791 1.21009273  1.0026922 1.4179148 0.10804486
#> CDE50% 1.23179114 1.23258318  1.0495918 1.4110840 0.09297286
#> CDE75% 1.23695452 1.23410534  1.0502652 1.4193225 0.09199348

Plotting

plotdf <- as.data.frame(mediationeffects.ey10$est)
plotdf["Effect"] <- rownames(plotdf)
ggplot(plotdf, aes(Effect, mean, ymin = lower, ymax = upper ))  + 
  geom_pointrange(position = position_dodge(width = 0.5))  + theme(axis.text = element_text(size = 5))

Summary statistics of the predictor-response function

Risk Summary For Totel Effect 

riskSummary10 = TERiskSummaries.CMA(fit.TE = fit.y.TE, e.y=e.y10, e.y.name = "E.Y", sel=sel)

p_riskSummary10 <-ggplot(riskSummary10,
       aes(quantile,
           est,
           ymin = est - 1.96 * sd,
           ymax = est + 1.96 * sd)) +
  geom_pointrange() + theme(axis.text = element_text(size = 5))

riskSummary90 = TERiskSummaries.CMA(fit.TE = fit.y.TE, e.y=e.y90, e.y.name = "E.Y", sel=sel)

p_riskSummary90 <- ggplot(riskSummary90,
       aes(quantile,
           est,
           ymin = est - 1.96 * sd,
           ymax = est + 1.96 * sd)) +
  geom_pointrange() + theme(axis.text = element_text(size = 5))

ggarrange(p_riskSummary10, p_riskSummary90, ncol=2, nrow =1)

Risk Summary For Controlled Direct Effect 

# CDE 
CDEriskSummary10 = CDERiskSummaries.CMA(fit.y = fit.y, e.y = e.y10, e.y.name = "E.Y", m.name = "m", sel = sel)
 
plot_CDErisk10 <- ggplot(CDEriskSummary10, aes(quantile, est, ymin = est - 1.96*sd,  ymax = est + 1.96*sd, col = m)) + theme(legend.position="none") + geom_pointrange(size = 0.5, position = position_dodge(width = 0.2))   

CDEriskSummary90 = CDERiskSummaries.CMA(fit.y = fit.y, e.y = e.y90, e.y.name = "E.Y", m.name = "m", sel = sel)

plot_CDErisk90 <- ggplot(CDEriskSummary90, aes(quantile, est, ymin = est - 1.96*sd,  ymax = est + 1.96*sd, col = m)) + 
  geom_pointrange(size = 0.5, position = position_dodge(width = 0.2))
ggarrange(plot_CDErisk10, plot_CDErisk90, ncol=2, nrow =1)

Single-predictor health risks

Total Effect 

risks.singvar10 = SingVarRiskSummaries.CMA(BKMRfits = fit.y.TE, which.z = 1:3,
                                           e.y=e.y10, e.y.names="E.Y",
                                           sel=sel)

plot_singvar10 <- ggplot(risks.singvar10, aes(variable, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd, col = q.fixed)) + theme(legend.position="none")+
  geom_pointrange(size= .3, position = position_dodge(width = 0.2)) +
  coord_flip()
 
plot_singvar10

Controlled Direct Effect 

# single variable controlled direct effects
CDErisks.singvar10 = CDESingVarRiskSummaries.CMA(BKMRfits = fit.y,
                                           e.y=e.y10, e.y.names="E.Y", m.name = "m",
                                           sel=sel) 
plotCDE10 <- ggplot(CDErisks.singvar10, aes(variable, est, ymin = est - 1.96*sd,
                            ymax = est + 1.96*sd, col = q.fixed, linetype = m.fixed)) +
  geom_pointrange(size = 0.5, position = position_dodge(width = 0.2)) +theme(legend.position="none")+
  coord_flip() 

plotCDE10