This document applies the code snippets from McGrath & Wang et al. (2026) to a simulated dataset.
library(CMAverse) # Used for cmest()
library(factoextra) # Used for fviz_eig()
library(gcdnet) # Used for elastic net for ERS
library(bkmr) # Used for kmbayes()
library(causalbkmr) # Used for mediation.bkmr()
# causalbkmr may need to be installed from GitHub:
# devtools::install_github("zc2326/causalbkmr")set.seed(123)
n <- 300
# Confounders
C1 <- rnorm(n); C2 <- rnorm(n)
# Correlated exposures
X1 <- 0.3 * C1 + 0.3 * C2 + rnorm(n)
X2 <- 0.3 * C1 + 0.3 * C2 + 0.7 * X1 + rnorm(n, 0, 0.5)
X3 <- 0.3 * C1 + 0.3 * C2 + rnorm(n)
X4 <- 0.3 * C1 + 0.3 * C2 + 0.7 * X3 + rnorm(n, 0, 0.5)
X5 <- 0.3 * C1 + 0.3 * C2 + rnorm(n)
# Mediator
M <- 0.4 * X1 + 0.3 * X3 + 0.2 * C1 + rnorm(n)
# Outcome
Y <- 0.4 * M + 0.3 * X1 + 0.2 * X3 + 0.2 * C2 + rnorm(n)
data <- data.frame(X1, X2, X3, X4, X5, M, Y, C1, C2)
exposure_names <- c("X1","X2","X3","X4","X5")
confounder_names <- c("C1","C2")
X <- as.matrix(data[, exposure_names])
C <- as.matrix(data[, confounder_names])est <- CMAverse::cmest(
data = data,
model = "rb",
outcome = "Y",
exposure = "X1",
mediator = "M",
basec = c("C1", "C2"),
EMint = FALSE,
mreg = list("linear"),
yreg = "linear",
estimation = "paramfunc",
inference = "delta",
mval = list(median(data$M))
)
summary(est)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ X1 + M + C1 + C2, family = gaussian(), data = getCall(x$reg.output$yreg)$data,
## weights = getCall(x$reg.output$yreg)$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.008349 0.059700 -0.140 0.889
## X1 0.260941 0.063948 4.080 5.79e-05 ***
## M 0.486537 0.062020 7.845 8.02e-14 ***
## C1 0.009204 0.066296 0.139 0.890
## C2 0.272679 0.062269 4.379 1.66e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.064442)
##
## Null deviance: 512.43 on 299 degrees of freedom
## Residual deviance: 314.01 on 295 degrees of freedom
## AIC: 877.06
##
## Number of Fisher Scoring iterations: 2
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ X1 + C1 + C2, family = gaussian(), data = getCall(x$reg.output$mreg[[1L]])$data,
## weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05000 0.05587 -0.895 0.371541
## X1 0.43956 0.05421 8.108 1.37e-14 ***
## C1 0.22522 0.06074 3.708 0.000249 ***
## C2 0.03651 0.05832 0.626 0.531758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.9348966)
##
## Null deviance: 372.84 on 299 degrees of freedom
## Residual deviance: 276.73 on 296 degrees of freedom
## AIC: 837.14
##
## Number of Fisher Scoring iterations: 2
##
##
## # Effect decomposition on the mean difference scale via the regression-based approach
##
## Closed-form parameter function estimation with
## delta method standard errors, confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## cde 0.26094 0.06395 0.13560 0.386 4.49e-05 ***
## pnde 0.26094 0.06395 0.13560 0.386 4.49e-05 ***
## tnde 0.26094 0.06395 0.13560 0.386 4.49e-05 ***
## pnie 0.21386 0.03793 0.13952 0.288 1.72e-08 ***
## tnie 0.21386 0.03793 0.13952 0.288 1.72e-08 ***
## te 0.47480 0.06358 0.35020 0.599 8.13e-14 ***
## pm 0.45042 0.08509 0.28365 0.617 1.20e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (cde: controlled direct effect; pnde: pure natural direct effect; tnde: total natural direct effect; pnie: pure natural indirect effect; tnie: total natural indirect effect; te: total effect; pm: overall proportion mediated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $mval
## $mval[[1]]
## [1] -0.05955207
##
##
## $basecval
## $basecval[[1]]
## [1] 0.03444141
##
## $basecval[[2]]
## [1] 0.009109353
est <- CMAverse::cmest(
data = data,
model = "rb",
outcome = "Y",
exposure = "X1",
mediator = "M",
basec = c("C1", "C2", "X2", "X3", "X4", "X5"),
EMint = FALSE,
mreg = list("linear"),
yreg = "linear",
estimation = "paramfunc",
inference = "delta",
mval = list(median(data$M))
)
summary(est)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ X1 + M + C1 + C2 + X2 + X3 + X4 + X5, family = gaussian(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.013752 0.058685 -0.234 0.81489
## X1 0.331062 0.103396 3.202 0.00152 **
## M 0.437357 0.063084 6.933 2.66e-11 ***
## C1 -0.021969 0.082070 -0.268 0.78913
## C2 0.216610 0.083429 2.596 0.00990 **
## X2 -0.071390 0.117034 -0.610 0.54234
## X3 0.226462 0.101222 2.237 0.02603 *
## X4 0.003979 0.118265 0.034 0.97318
## X5 0.071976 0.058284 1.235 0.21786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.024344)
##
## Null deviance: 512.43 on 299 degrees of freedom
## Residual deviance: 298.08 on 291 degrees of freedom
## AIC: 869.44
##
## Number of Fisher Scoring iterations: 2
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ X1 + C1 + C2 + X2 + X3 + X4 + X5, family = gaussian(),
## data = getCall(x$reg.output$mreg[[1L]])$data, weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05140 0.05436 -0.946 0.345135
## X1 0.31738 0.09410 3.373 0.000845 ***
## C1 0.18318 0.07537 2.430 0.015690 *
## C2 -0.03393 0.07737 -0.439 0.661339
## X2 0.16296 0.10815 1.507 0.132937
## X3 0.32053 0.09201 3.484 0.000570 ***
## X4 -0.13186 0.10944 -1.205 0.229224
## X5 -0.03875 0.05402 -0.717 0.473732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.8815097)
##
## Null deviance: 372.84 on 299 degrees of freedom
## Residual deviance: 257.40 on 292 degrees of freedom
## AIC: 823.42
##
## Number of Fisher Scoring iterations: 2
##
##
## # Effect decomposition on the mean difference scale via the regression-based approach
##
## Closed-form parameter function estimation with
## delta method standard errors, confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## cde 0.33106 0.10340 0.12841 0.534 0.00137 **
## pnde 0.33106 0.10340 0.12841 0.534 0.00137 **
## tnde 0.33106 0.10340 0.12841 0.534 0.00137 **
## pnie 0.13881 0.04577 0.04911 0.229 0.00242 **
## tnie 0.13881 0.04577 0.04911 0.229 0.00242 **
## te 0.46987 0.10947 0.25532 0.684 1.77e-05 ***
## pm 0.29542 0.09845 0.10247 0.488 0.00269 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (cde: controlled direct effect; pnde: pure natural direct effect; tnde: total natural direct effect; pnie: pure natural indirect effect; tnie: total natural indirect effect; te: total effect; pm: overall proportion mediated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $mval
## $mval[[1]]
## [1] -0.05955207
##
##
## $basecval
## $basecval[[1]]
## [1] 0.03444141
##
## $basecval[[2]]
## [1] 0.009109353
##
## $basecval[[3]]
## [1] 0.04743073
##
## $basecval[[4]]
## [1] 0.02314481
##
## $basecval[[5]]
## [1] 0.05519588
##
## $basecval[[6]]
## [1] 0.01961338
# Step 1: Perform unsupervised PCA on the standardized exposure matrix
res_pca <- prcomp(X, center = TRUE, scale. = TRUE)
# Step 2: Visualize explained variance to decide how many PCs to retain
suppressWarnings(print(factoextra::fviz_eig(res_pca)))
# Step 3: Extract PC scores that explain 80% of the total variance
explained_variance <- summary(res_pca)$importance[3, ]
id_80var <- which.max(explained_variance >= 0.8)
PC_scores <- res_pca$x[, 1:id_80var]
# Step 4: Append PC scores to the dataset
colnames(PC_scores) <- paste0("PC", 1:id_80var)
data <- cbind(data, PC_scores)
# Step 5: Perform mediation analysis using the first PC as the exposure
est <- CMAverse::cmest(
data = data,
model = "rb",
outcome = "Y",
exposure = "PC1",
mediator = "M",
basec = c("C1", "C2", paste0("PC", 2:id_80var)),
EMint = FALSE,
mreg = list("linear"),
yreg = "linear",
estimation = "paramfunc",
inference = "delta",
mval = list(median(data$M))
)
summary(est)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ PC1 + M + C1 + C2 + PC2 + PC3, family = gaussian(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.004063 0.058892 0.069 0.945
## PC1 -0.285505 0.055139 -5.178 4.18e-07 ***
## M 0.455116 0.062684 7.260 3.49e-12 ***
## C1 -0.111994 0.072713 -1.540 0.125
## C2 0.118654 0.072552 1.635 0.103
## PC2 0.035094 0.051780 0.678 0.498
## PC3 -0.019956 0.063023 -0.317 0.752
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.037004)
##
## Null deviance: 512.43 on 299 degrees of freedom
## Residual deviance: 303.84 on 293 degrees of freedom
## AIC: 871.18
##
## Number of Fisher Scoring iterations: 2
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ PC1 + C1 + C2 + PC2 + PC3, family = gaussian(),
## data = getCall(x$reg.output$mreg[[1L]])$data, weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03118 0.05476 -0.569 0.569579
## PC1 -0.36800 0.04660 -7.898 5.71e-14 ***
## C1 0.10044 0.06740 1.490 0.137240
## C2 -0.12679 0.06710 -1.890 0.059790 .
## PC2 0.15740 0.04729 3.328 0.000985 ***
## PC3 0.11025 0.05828 1.892 0.059516 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.8976786)
##
## Null deviance: 372.84 on 299 degrees of freedom
## Residual deviance: 263.92 on 294 degrees of freedom
## AIC: 826.92
##
## Number of Fisher Scoring iterations: 2
##
##
## # Effect decomposition on the mean difference scale via the regression-based approach
##
## Closed-form parameter function estimation with
## delta method standard errors, confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## cde -0.28550 0.05514 -0.39357 -0.177 2.24e-07 ***
## pnde -0.28550 0.05514 -0.39357 -0.177 2.24e-07 ***
## tnde -0.28550 0.05514 -0.39357 -0.177 2.24e-07 ***
## pnie -0.16748 0.03133 -0.22890 -0.106 9.04e-08 ***
## tnie -0.16748 0.03133 -0.22890 -0.106 9.04e-08 ***
## te -0.45299 0.05439 -0.55958 -0.346 < 2e-16 ***
## pm 0.36973 0.07166 0.22928 0.510 2.47e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (cde: controlled direct effect; pnde: pure natural direct effect; tnde: total natural direct effect; pnie: pure natural indirect effect; tnie: total natural indirect effect; te: total effect; pm: overall proportion mediated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $mval
## $mval[[1]]
## [1] -0.05955207
##
##
## $basecval
## $basecval[[1]]
## [1] 0.03444141
##
## $basecval[[2]]
## [1] 0.009109353
##
## $basecval[[3]]
## [1] -1.438519e-17
##
## $basecval[[4]]
## [1] 2.031397e-17
# Split data into training and analysis sets
n_train <- floor(nrow(data) / 2)
train_ids <- sample(seq_len(nrow(data)), n_train)
analysis_ids <- setdiff(seq_len(nrow(data)), train_ids)
X_train <- as.matrix(data[train_ids, exposure_names])
C_train <- as.matrix(data[train_ids, confounder_names])
Y_train <- data[train_ids, "Y"]
# Combine exposures and confounders from the training set into a matrix
Z_train <- cbind(X_train, C_train)
# Penalty factors for L1 and L2 regularization: Exposures are penalized (1),
# confounders are unpenalized (0)
pf <- c(rep(1, ncol(X_train)), rep(0, ncol(C_train)))
# Range of lambda2 values for L2 regularization
lambda2_values <- exp(seq(log(1e-4), log(1e2), length.out = 100))
# Select optimal lambda2 using cross-validation
cv_lambda2 <- sapply(lambda2_values, function(lambda2) {
min(cv.gcdnet(x = Z_train, y = Y_train, lambda2 = lambda2, method = "ls",
nfolds = 5, pf = pf, pf2 = pf)$cvm)
})
optimal_lambda2 <- lambda2_values[which.min(cv_lambda2)]
# Select optimal lambda1 for the chosen lambda2
cv_result <- cv.gcdnet(x = Z_train, y = Y_train, lambda2 = optimal_lambda2,
method = "ls", nfolds = 5, pf = pf, pf2 = pf)
optimal_lambda1 <- cv_result$lambda.min
# Fit final elastic net model
best_model <- gcdnet(x = Z_train, y = Y_train,
lambda = optimal_lambda1, lambda2 = optimal_lambda2,
method = "ls", pf = pf, pf2 = pf)
# Extract exposure coefficients (exclude intercept and confounders)
beta_hat <- coef(best_model)
exposure_coef <- beta_hat[exposure_names, , drop = FALSE]
# Compute ERS in the analysis set
ERS <- as.matrix(data[analysis_ids, exposure_names]) %*% exposure_coef
data_analysis <- data[analysis_ids, ]
data_analysis$ERS <- as.numeric(ERS)
# Perform mediation analysis using ERS as the exposure
est <- CMAverse::cmest(
data = data_analysis,
model = "rb",
outcome = "Y",
exposure = "ERS",
mediator = "M",
basec = confounder_names,
EMint = FALSE,
mreg = list("linear"),
yreg = "linear",
estimation = "paramfunc",
inference = "delta",
mval = list(median(data_analysis$M))
)
summary(est)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ ERS + M + C1 + C2, family = gaussian(), data = getCall(x$reg.output$yreg)$data,
## weights = getCall(x$reg.output$yreg)$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.12996 0.08288 -1.568 0.119030
## ERS 0.47703 0.12917 3.693 0.000313 ***
## M 0.36558 0.08372 4.367 2.38e-05 ***
## C1 0.10153 0.09415 1.078 0.282663
## C2 0.26268 0.09051 2.902 0.004282 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.015232)
##
## Null deviance: 230.59 on 149 degrees of freedom
## Residual deviance: 147.21 on 145 degrees of freedom
## AIC: 434.86
##
## Number of Fisher Scoring iterations: 2
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ ERS + C1 + C2, family = gaussian(), data = getCall(x$reg.output$mreg[[1L]])$data,
## weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.06544 0.08175 -0.800 0.425
## ERS 0.55824 0.11905 4.689 6.24e-06 ***
## C1 0.11055 0.09262 1.194 0.235
## C2 0.00241 0.08947 0.027 0.979
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.9921601)
##
## Null deviance: 173.29 on 149 degrees of freedom
## Residual deviance: 144.86 on 146 degrees of freedom
## AIC: 430.45
##
## Number of Fisher Scoring iterations: 2
##
##
## # Effect decomposition on the mean difference scale via the regression-based approach
##
## Closed-form parameter function estimation with
## delta method standard errors, confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## cde 0.47703 0.12917 0.22385 0.730 0.000222 ***
## pnde 0.47703 0.12917 0.22385 0.730 0.000222 ***
## tnde 0.47703 0.12917 0.22385 0.730 0.000222 ***
## pnie 0.20408 0.06386 0.07892 0.329 0.001395 **
## tnie 0.20408 0.06386 0.07892 0.329 0.001395 **
## te 0.68112 0.12805 0.43015 0.932 1.04e-07 ***
## pm 0.29963 0.09756 0.10843 0.491 0.002131 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (cde: controlled direct effect; pnde: pure natural direct effect; tnde: total natural direct effect; pnie: pure natural indirect effect; tnie: total natural indirect effect; te: total effect; pm: overall proportion mediated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $mval
## $mval[[1]]
## [1] -0.1539442
##
##
## $basecval
## $basecval[[1]]
## [1] -0.03930162
##
## $basecval[[2]]
## [1] -0.02632793
Note: The analyses below use 100 MCMC iterations (iter = 100) for
computational simplicity. A larger number of MCMC iterations is needed
in practice.
fit.y <- bkmr::kmbayes(y = Y, Z = cbind(X, M), X = C,
verbose = FALSE, varsel = TRUE,
iter = 100)
fit.m <- bkmr::kmbayes(y = M, Z = X, X = C,
verbose = FALSE, varsel = TRUE,
iter = 100)
fit.y.TE <- bkmr::kmbayes(y = Y, Z = X, X = C,
verbose = FALSE, varsel = TRUE,
iter = 100)cor_mat <- cor(X, method = "pearson")
hc <- hclust(as.dist(1 - cor_mat))
groups <- cutree(hc, k = 3)
fit.y.hier <- bkmr::kmbayes(y = Y, Z = cbind(X, M), X = C,
verbose = FALSE, varsel = TRUE,
iter = 100,
groups = c(groups, max(groups) + 1))
fit.m.hier <- bkmr::kmbayes(y = M, Z = X, X = C,
verbose = FALSE, varsel = TRUE,
iter = 100,
groups = groups)
fit.y.TE.hier <- bkmr::kmbayes(y = Y, Z = X, X = C,
verbose = FALSE, varsel = TRUE,
iter = 100,
groups = groups)# Define the MCMC iterations to be used for inference
sel <- seq(51, 100)
# Specify exposure levels to compare (e.g., 25th vs. 75th percentile)
a <- apply(X, 2, quantile, probs = 0.75)
astar <- apply(X, 2, quantile, probs = 0.25)
# Specify covariates used for predicting M and Y (e.g., mean of confounders)
X.predict <- matrix(colMeans(C), nrow = 1)
# Estimate TE, NDE, and NIE based on component-wise selection
mediation.effects <- causalbkmr::mediation.bkmr(
a = a,
astar = astar,
e.y = NULL,
e.m = NULL,
fit.m = fit.m,
fit.y = fit.y,
fit.y.TE = fit.y.TE,
X.predict.M = X.predict,
X.predict.Y = X.predict,
m.quant = c(0.1, 0.25, 0.5, 0.75),
alpha = 0.05,
sel = sel,
K = 50,
seed = 123
)## [1] "iter 50 time: 0.01 min"
## [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"
mediation.effects## $est
## mean median lower upper sd
## TE 1.0827234 1.1851573 -0.2899952 2.144708 0.6760127
## NDE 0.8552039 0.8789151 -0.5139456 1.808211 0.6718745
## NIE 0.2275195 0.2300885 -0.7829272 1.174658 0.5127300
## CDE10% 1.1412329 1.3025811 -0.9697934 3.063295 1.1232726
## CDE25% 0.4773654 0.5980800 -1.4208777 1.994660 0.9247246
## CDE50% 0.5945804 0.7181701 -1.4828886 2.047765 0.9633614
## CDE75% 0.9629236 1.0097753 -1.0860203 2.486692 0.9565665
##
## $TE.samp
## [1] 0.84317368 1.75330013 0.10936702 2.03212961 0.91259431 1.37766924
## [7] 1.12616418 -0.31697174 2.36359285 1.56076278 0.66430471 1.00072450
## [13] 1.70564450 1.40390460 -0.19707600 1.45610621 1.01601105 1.03463970
## [19] 1.31770434 1.16893555 0.84158522 -0.71285242 2.16397234 1.15485192
## [25] 1.48492419 1.20137902 0.27311937 0.12469439 0.04653138 1.04462392
## [31] 1.48197395 1.41815850 0.42450169 1.25087405 0.38146593 2.07835370
## [37] 1.86099932 0.18930571 1.57206676 1.28186940 0.96602932 0.58740117
## [43] 0.94522437 1.34761432 0.45676358 1.31362867 1.62665588 2.07217706
## [49] 1.33299089 1.59260457
##
## $NDE.samp
## [1] 1.7979905 0.8713656 -0.8083674 1.5603686 1.1932170 0.6528135
## [7] 1.1009084 -0.2896902 2.7524729 0.9441043 0.9520047 0.6381791
## [13] 1.6345224 1.2427659 0.2061380 1.2417740 0.2836434 0.3585799
## [19] 1.0485565 1.1060023 1.6444189 -0.5790519 1.7111085 0.8864645
## [25] 0.2582880 1.1804535 0.1555098 0.5882798 0.7608917 0.5174578
## [31] 1.0047197 1.3402640 0.5453075 0.6365460 -0.1010308 1.8111781
## [37] 1.2365112 0.4819407 1.3262220 1.6625145 0.4701675 0.8019439
## [43] 0.7962546 0.3519915 0.2039394 0.6741085 1.4615554 1.1682927
## [49] -0.1987476 1.4753483
##
## $NIE.samp
## [1] -0.95481679 0.88193450 0.91773447 0.47176102 -0.28062267 0.72485577
## [7] 0.02525579 -0.02728149 -0.38888008 0.61665847 -0.28769995 0.36254543
## [13] 0.07112211 0.16113870 -0.40321403 0.21433220 0.73236769 0.67605984
## [19] 0.26914780 0.06293328 -0.80283366 -0.13380049 0.45286389 0.26838743
## [25] 1.22663619 0.02092552 0.11760958 -0.46358544 -0.71436028 0.52716613
## [31] 0.47725427 0.07789448 -0.12080578 0.61432806 0.48249675 0.26717560
## [37] 0.62448812 -0.29263497 0.24584479 -0.38064509 0.49586184 -0.21454274
## [43] 0.14896979 0.99562279 0.25282416 0.63952016 0.16510043 0.90388440
## [49] 1.53173850 0.11725632
##
## $`CDE10%.samp`
## [1] 0.77109171 2.74145404 0.08355497 2.64224536 0.51795718 2.07345418
## [7] 1.28488925 -0.57857053 2.92712238 2.12824338 0.22371219 0.50161859
## [13] 1.51453941 1.90531943 -0.99584239 1.90944128 1.40927553 1.60317126
## [19] 1.77547804 1.20524020 0.59818118 -2.08155483 3.13997234 1.17184043
## [25] 1.75843249 1.22860592 0.08093740 -0.09862687 -0.88006925 0.96753335
## [31] 1.79278556 1.32027295 -0.28095892 1.49279451 0.43914853 1.96745380
## [37] 2.50212037 -0.51218785 1.51047915 1.21203114 0.85204135 0.28852802
## [43] 0.45512148 1.97986841 0.06213734 1.75257650 1.57734460 3.10282940
## [49] 2.25091124 1.76770003
##
## $`CDE25%.samp`
## [1] 0.20036158 1.58520652 -0.68947264 1.81132086 -0.05816402 0.96121846
## [7] 0.38965183 -1.24092528 2.02918222 1.15205050 -0.27449157 -0.03972426
## [13] 0.87953141 0.91017280 -1.47312199 1.02901873 0.38550452 0.60154787
## [19] 0.87190825 0.54929039 0.05043120 -2.46228316 2.06113116 0.48845227
## [25] 0.99095111 0.67313021 -0.33380146 -0.41875288 -0.62770700 0.48245797
## [31] 1.21076424 1.04261071 -0.18813357 0.90609686 -0.41152923 1.72398410
## [37] 1.61773303 -0.70280166 1.14328247 0.59461217 0.33884630 -0.08780490
## [43] 0.01542559 0.92449505 -0.33404445 0.78993465 0.90260928 1.87574825
## [49] 0.75961547 1.26275114
##
## $`CDE50%.samp`
## [1] 0.16908166 1.32652105 -0.89258189 1.70908795 0.03502249 0.77186283
## [7] 0.27082471 -1.31831858 2.06139704 1.04052604 -0.22875437 0.12768548
## [13] 1.00275108 0.69882368 -1.53066693 0.91568262 0.27897260 0.42695830
## [19] 0.73751653 0.52602336 0.07226569 -2.54139828 1.99925923 0.51551345
## [25] 1.16628192 0.84870976 -0.22789956 -0.25176519 -0.49033721 0.68249761
## [31] 1.41008596 1.28001515 0.03390347 1.17198691 -0.21465840 2.00081038
## [37] 2.00065141 -0.46325770 1.46098785 1.05749660 0.79473994 0.24078643
## [43] 0.48330624 1.31381538 -0.02386153 1.15206171 1.32618209 2.15092750
## [49] 1.16241015 1.48908867
##
## $`CDE75%.samp`
## [1] 0.6810355 1.4758735 -0.8007080 2.1643009 0.7439076 0.9923627
## [7] 0.7806869 -1.0896949 2.6695565 1.5105997 0.3802439 0.8381000
## [13] 1.8884904 1.0119390 -1.0733636 1.2263517 0.3828786 0.5144784
## [19] 1.0722636 0.9603551 0.5972606 -2.1721101 2.3752704 1.0076116
## [25] 1.6639247 1.3996299 0.3124645 0.1809480 0.1087142 1.1585348
## [31] 1.7630139 1.7749787 0.4573660 1.2837258 0.1179423 2.5190406
## [37] 2.1243993 0.0827436 1.8393345 1.3666018 0.8818894 0.4850066
## [43] 0.7573544 1.4230763 0.1870290 1.2610505 1.6524208 2.3252145
## [49] 1.1311038 1.7509801
# Estimate TE, NDE, and NIE based on hierarchical selection
mediation.effects.hier <- causalbkmr::mediation.bkmr(
a = a,
astar = astar,
e.y = NULL,
e.m = NULL,
fit.m = fit.m.hier,
fit.y = fit.y.hier,
fit.y.TE = fit.y.TE.hier,
X.predict.M = X.predict,
X.predict.Y = X.predict,
m.quant = c(0.1, 0.25, 0.5, 0.75),
alpha = 0.05,
sel = sel,
K = 50,
seed = 123
)## [1] "iter 50 time: 0.01 min"
## [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"
mediation.effects.hier## $est
## mean median lower upper sd
## TE 1.1006171 1.1198304 0.6139017 1.5359266 0.2484828
## NDE 0.8879738 0.9184710 0.2362151 1.3135324 0.2879491
## NIE 0.2126433 0.2372203 -0.2059980 0.5487169 0.2363560
## CDE10% 0.9337695 0.9231778 -0.2924228 2.0962866 0.6735259
## CDE25% 0.7169421 0.7012659 -0.3697995 1.8263375 0.5767113
## CDE50% 0.4027754 0.3891099 -0.6094594 1.5735345 0.5776849
## CDE75% 0.7285577 0.8001787 -0.1351033 2.0100127 0.6445774
##
## $TE.samp
## [1] 1.1229820 1.3753847 0.6268104 1.3450677 0.9409682 1.2213496 1.0683931
## [8] 0.6101540 1.4933899 1.2670633 0.9030157 1.0219524 1.2386694 1.1919661
## [15] 0.6754900 1.1653084 1.0204238 1.0310707 1.1351700 0.9599263 0.9680148
## [22] 0.5033966 1.5482759 1.1134749 1.1566209 1.2789354 0.9239718 0.7975352
## [29] 0.8148052 1.1166788 1.3524915 1.3258421 1.0177915 1.1416460 1.0512334
## [36] 1.6117038 1.4895120 0.6884614 1.3274472 1.2795564 0.9884702 0.9643126
## [43] 0.9988417 1.2587459 0.7779854 1.1634378 1.2678289 1.4171335 1.1156241
## [50] 1.1565247
##
## $NDE.samp
## [1] 1.3309913 0.8251250 0.2055898 1.0176944 1.1151874 0.8669919 1.0397569
## [8] 0.4005002 1.6813586 0.9314539 1.0714752 0.7972334 1.1001672 1.1525718
## [15] 0.5998251 1.1710317 0.7177078 0.5745398 1.1295274 0.8673662 1.1827790
## [22] 0.3080485 1.2011453 0.7329119 0.2153603 1.2390642 0.6742503 0.9966054
## [29] 0.8963270 0.7312620 1.0333445 1.2533960 0.9842715 0.7062055 0.7140805
## [36] 1.2180488 1.1172894 0.6862273 0.9249049 0.9645722 0.5651836 0.8361146
## [43] 0.7835598 0.9391096 0.6048349 0.6200353 0.9274693 0.9120372 0.8335835
## [50] 1.0005721
##
## $NIE.samp
## [1] -0.208009334 0.550259724 0.421220589 0.327373296 -0.174219255
## [6] 0.354357706 0.028636204 0.209653817 -0.187968604 0.335609396
## [11] -0.168459420 0.224719034 0.138502241 0.039394362 0.075664916
## [16] -0.005723292 0.302716070 0.456530919 0.005642583 0.092560113
## [21] -0.214764229 0.195348074 0.347130639 0.380562965 0.941260597
## [26] 0.039871218 0.249721488 -0.199070132 -0.081521881 0.385416792
## [31] 0.319147087 0.072446097 0.033519939 0.435440438 0.337152925
## [36] 0.393655073 0.372222618 0.002234086 0.402542295 0.314984252
## [41] 0.423286543 0.128197962 0.215281955 0.319636237 0.173150544
## [46] 0.543402534 0.340359667 0.505096333 0.282040591 0.155952603
##
## $`CDE10%.samp`
## [1] 0.6069784 1.8329932 0.2374303 1.7566095 0.5306964 1.5030103
## [7] 1.1238946 -0.2945796 2.1190956 1.5136034 0.3941415 0.6753610
## [13] 1.2923502 1.4638451 -0.2849936 1.3410020 1.1443105 1.2455426
## [19] 1.4782791 0.9812801 0.5451865 -0.7267107 2.2918871 0.9417423
## [25] 1.4344747 1.1546371 0.3869453 0.1805910 -0.1888383 0.9159801
## [31] 2.0177219 1.8439008 0.9303756 1.6214892 1.1022269 1.7252844
## [37] 1.9749615 0.6518904 0.5083504 0.3028902 0.2404177 0.4393557
## [43] 0.4903429 0.6421231 0.3706082 0.6320940 0.8099560 1.0841147
## [49] 0.8519078 0.8517181
##
## $`CDE25%.samp`
## [1] 0.59608541 1.52102566 0.09382949 1.56810105 0.57510049 1.18780377
## [7] 0.91342185 -0.32452431 1.89132234 1.31696905 0.37229280 0.62152409
## [13] 1.12098693 1.09147440 -0.38294396 1.01173894 0.72610815 0.83834751
## [19] 1.14889165 0.79693281 0.43215038 -0.77070359 1.87175720 0.79155654
## [25] 1.08429806 0.93536571 0.07798124 -0.01816479 -0.22135127 0.68452247
## [31] 1.66989207 1.55881752 0.85722286 1.09112530 0.71800940 1.31400897
## [37] 1.38325697 0.46143183 0.38520988 0.30075311 0.24359533 0.32256910
## [43] 0.35734404 0.50566411 0.24157243 0.49336724 0.51042504 0.72315803
## [49] 0.55489051 0.60288920
##
## $`CDE50%.samp`
## [1] 0.51946887 1.37063519 -0.13567802 1.43270459 0.57012230 0.94532136
## [7] 0.67279625 -0.50742001 1.70719510 1.20169470 0.33260170 0.56606144
## [13] 0.92912216 0.79873430 -0.63908369 0.70099086 0.31594094 0.44754919
## [19] 0.89524060 0.62289696 0.28957972 -0.96665691 1.61442061 0.66833817
## [25] 0.80776790 0.79153677 -0.37018236 -0.28730898 -0.43710716 0.44561800
## [31] 0.85112005 0.82397834 0.16927262 0.68281599 0.23825179 0.80907992
## [37] 0.96858247 0.06813117 0.03950647 0.16561911 0.11149113 -0.10459644
## [43] -0.08713172 0.04211127 -0.16476644 0.06463505 -0.03667983 0.14596664
## [49] -0.03156179 0.08004367
##
## $`CDE75%.samp`
## [1] 0.9647135008 1.6618021694 0.2466463346 1.9036314896 1.0192641986
## [6] 1.3439262193 1.1365741522 -0.0814464581 2.2509611650 1.5853343098
## [11] 0.7914734756 0.9861034588 1.4875842106 1.2311935310 -0.1386332876
## [16] 1.2082289271 0.8088838522 0.9196820654 1.2917714869 1.0668149166
## [21] 0.7666063189 -0.5579018828 2.0408975652 1.0874065348 1.2708745633
## [26] 1.1577456436 0.1718032335 0.1247003147 0.0787241170 0.8930388033
## [31] 1.1784238689 1.2492968511 0.4610935560 0.7857976463 0.0991333227
## [36] 0.9713019391 1.0779004589 -0.0156267430 0.0904223770 0.1031780485
## [41] 0.0383218566 -0.0381028972 -0.0004428131 0.1204488440 -0.1229444213
## [46] 0.1071940150 0.3489740915 0.5258944934 0.2888792521 0.4403371451