Statistical Modeling with a Single Mediator
2022-03-30
Source:vignettes/single_mediator.Rmd
single_mediator.RmdThis example demonstrates how to use cmest when there is a single mediator. For this purpose, we simulate some data containing a continuous baseline confounder \(C_1\), a binary baseline confounder \(C_2\), a binary exposure \(A\), a binary mediator \(M\) and a binary outcome \(Y\). The true regression models for \(A\), \(M\) and \(Y\) are: \[logit(E(A|C_1,C_2))=0.2+0.5C_1+0.1C_2\] \[logit(E(M|A,C_1,C_2))=1+2A+1.5C_1+0.8C_2\] \[logit(E(Y|A,M,C_1,C_2)))=-3-0.4A-1.2M+0.5AM+0.3C_1-0.6C_2\]
## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: mstate## Loading required package: survival## Loading required package: tidyr
set.seed(1)
expit <- function(x) exp(x)/(1+exp(x))
n <- 10000
C1 <- rnorm(n, mean = 1, sd = 0.1)
C2 <- rbinom(n, 1, 0.6)
A <- rbinom(n, 1, expit(0.2 + 0.5*C1 + 0.1*C2))
M <- rbinom(n, 1, expit(1 + 2*A + 1.5*C1 + 0.8*C2))
Y <- rbinom(n, 1, expit(-3 - 0.4*A - 1.2*M + 0.5*A*M + 0.3*C1 - 0.6*C2))
data <- data.frame(A, M, Y, C1, C2)The DAG for this scientific setting is:
cmdag(outcome = "Y", exposure = "A", mediator = "M",
basec = c("C1", "C2"), postc = NULL, node = TRUE, text_col = "white")
In this setting, we can use all of the six statistical modeling approaches. The results are shown below:
The Regression-based Approach
Closed-form Parameter Function Estimation and Delta Method Inference
res_rb_param_delta <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
mreg = list("logistic"), yreg = "logistic",
astar = 0, a = 1, mval = list(1),
estimation = "paramfunc", inference = "delta")
summary(res_rb_param_delta)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ A + M + A * M + C1 + C2, family = binomial(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5681 -0.2414 -0.1741 -0.1604 3.0521
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5731 0.7713 -4.633 3.61e-06 ***
## A 0.7084 0.5283 1.341 0.17990
## M -1.0471 0.3736 -2.803 0.00507 **
## C1 0.9337 0.6973 1.339 0.18054
## C2 -0.8415 0.1443 -5.829 5.56e-09 ***
## A:M -0.4222 0.5541 -0.762 0.44607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2022.7 on 9999 degrees of freedom
## Residual deviance: 1965.0 on 9994 degrees of freedom
## AIC: 1977
##
## Number of Fisher Scoring iterations: 7
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ A + C1 + C2, family = binomial(), data = getCall(x$reg.output$mreg[[1L]])$data,
## weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2578 0.1145 0.1647 0.2519 0.5563
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4353 0.6410 0.679 0.49712
## A 1.7076 0.1443 11.837 < 2e-16 ***
## C1 1.9982 0.6474 3.086 0.00203 **
## C2 0.9024 0.1345 6.709 1.96e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2294.0 on 9999 degrees of freedom
## Residual deviance: 2072.5 on 9996 degrees of freedom
## AIC: 2080.5
##
## Number of Fisher Scoring iterations: 7
##
##
## # Effect decomposition on the risk ratio 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
## Rcde 1.33138 0.22278 0.95912 1.848 0.0872 .
## Rpnde 1.41988 0.23708 1.02358 1.970 0.0358 *
## Rtnde 1.34928 0.22056 0.97940 1.859 0.0669 .
## Rpnie 0.93338 0.03576 0.86587 1.006 0.0719 .
## Rtnie 0.88697 0.05295 0.78902 0.997 0.0445 *
## Rte 1.25939 0.19804 0.92535 1.714 0.1425
## ERcde 0.30416 0.20019 -0.08820 0.697 0.1287
## ERintref 0.11571 0.11472 -0.10914 0.341 0.3131
## ERintmed -0.09387 0.09323 -0.27659 0.089 0.3140
## ERpnie -0.06662 0.03576 -0.13670 0.003 0.0624 .
## ERcde(prop) 1.17262 0.23540 0.71124 1.634 6.32e-07 ***
## ERintref(prop) 0.44610 0.49051 -0.51528 1.407 0.3631
## ERintmed(prop) -0.36189 0.39881 -1.14354 0.420 0.3642
## ERpnie(prop) -0.25684 0.23474 -0.71691 0.203 0.2739
## pm -0.61872 0.50289 -1.60437 0.367 0.2186
## int 0.08421 0.09264 -0.09736 0.266 0.3633
## pe -0.17262 0.23540 -0.63400 0.289 0.4634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rcde: controlled direct effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnde: total natural direct effect risk ratio; Rpnie: pure natural indirect effect risk ratio; Rtnie: total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; ERintref: excess relative risk due to reference interaction; ERintmed: excess relative risk due to mediated interaction; ERpnie: excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; ERintref(prop): proportion ERintref; ERintmed(prop): proportion ERintmed; ERpnie(prop): proportion ERpnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"
##
## $mval
## $mval[[1]]
## [1] 1
##
##
## $basecval
## $basecval[[1]]
## [1] 0.9993463
##
## $basecval[[2]]
## [1] 0.6061Closed-form Parameter Function Estimation and Bootstrap Inference
res_rb_param_bootstrap <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
mreg = list("logistic"), yreg = "logistic",
astar = 0, a = 1, mval = list(1),
estimation = "paramfunc", inference = "bootstrap", nboot = 2)
summary(res_rb_param_bootstrap)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ A + M + A * M + C1 + C2, family = binomial(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5681 -0.2414 -0.1741 -0.1604 3.0521
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5731 0.7713 -4.633 3.61e-06 ***
## A 0.7084 0.5283 1.341 0.17990
## M -1.0471 0.3736 -2.803 0.00507 **
## C1 0.9337 0.6973 1.339 0.18054
## C2 -0.8415 0.1443 -5.829 5.56e-09 ***
## A:M -0.4222 0.5541 -0.762 0.44607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2022.7 on 9999 degrees of freedom
## Residual deviance: 1965.0 on 9994 degrees of freedom
## AIC: 1977
##
## Number of Fisher Scoring iterations: 7
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ A + C1 + C2, family = binomial(), data = getCall(x$reg.output$mreg[[1L]])$data,
## weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2578 0.1145 0.1647 0.2519 0.5563
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4353 0.6410 0.679 0.49712
## A 1.7076 0.1443 11.837 < 2e-16 ***
## C1 1.9982 0.6474 3.086 0.00203 **
## C2 0.9024 0.1345 6.709 1.96e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2294.0 on 9999 degrees of freedom
## Residual deviance: 2072.5 on 9996 degrees of freedom
## AIC: 2080.5
##
## Number of Fisher Scoring iterations: 7
##
##
## # Effect decomposition on the risk ratio scale via the regression-based approach
##
## Closed-form parameter function estimation with
## bootstrap standard errors, percentile confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## Rcde 1.331376 0.025783 1.168386 1.203 <2e-16 ***
## Rpnde 1.419875 0.122168 1.268237 1.432 <2e-16 ***
## Rtnde 1.349275 0.042976 1.188801 1.247 <2e-16 ***
## Rpnie 0.933380 0.002995 0.894425 0.898 <2e-16 ***
## Rtnie 0.886970 0.042062 0.781893 0.838 <2e-16 ***
## Rte 1.259387 0.042173 1.063293 1.120 <2e-16 ***
## ERcde 0.304162 0.023479 0.146643 0.178 <2e-16 ***
## ERintref 0.115713 0.098689 0.121859 0.254 <2e-16 ***
## ERintmed -0.093869 0.082990 -0.211091 -0.100 <2e-16 ***
## ERpnie -0.066620 0.002995 -0.105575 -0.102 <2e-16 ***
## ERcde(prop) 1.172621 0.625254 1.495704 2.336 <2e-16 ***
## ERintref(prop) 0.446102 0.147895 1.919311 2.118 <2e-16 ***
## ERintmed(prop) -0.361887 0.140542 -1.756808 -1.568 <2e-16 ***
## ERpnie(prop) -0.256836 0.617901 -1.687057 -0.857 <2e-16 ***
## pm -0.618723 0.477360 -3.255047 -2.614 <2e-16 ***
## int 0.084215 0.007353 0.351321 0.361 <2e-16 ***
## pe -0.172621 0.625254 -1.335736 -0.496 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rcde: controlled direct effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnde: total natural direct effect risk ratio; Rpnie: pure natural indirect effect risk ratio; Rtnie: total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; ERintref: excess relative risk due to reference interaction; ERintmed: excess relative risk due to mediated interaction; ERpnie: excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; ERintref(prop): proportion ERintref; ERintmed(prop): proportion ERintmed; ERpnie(prop): proportion ERpnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"
##
## $mval
## $mval[[1]]
## [1] 1
##
##
## $basecval
## $basecval[[1]]
## [1] 0.9993463
##
## $basecval[[2]]
## [1] 0.6061Direct Counterfactual Imputation Estimation and Bootstrap Inference
res_rb_impu_bootstrap <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
mreg = list("logistic"), yreg = "logistic",
astar = 0, a = 1, mval = list(1),
estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_rb_impu_bootstrap)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ A + M + A * M + C1 + C2, family = binomial(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5681 -0.2414 -0.1741 -0.1604 3.0521
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5731 0.7713 -4.633 3.61e-06 ***
## A 0.7084 0.5283 1.341 0.17990
## M -1.0471 0.3736 -2.803 0.00507 **
## C1 0.9337 0.6973 1.339 0.18054
## C2 -0.8415 0.1443 -5.829 5.56e-09 ***
## A:M -0.4222 0.5541 -0.762 0.44607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2022.7 on 9999 degrees of freedom
## Residual deviance: 1965.0 on 9994 degrees of freedom
## AIC: 1977
##
## Number of Fisher Scoring iterations: 7
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ A + C1 + C2, family = binomial(), data = getCall(x$reg.output$mreg[[1L]])$data,
## weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2578 0.1145 0.1647 0.2519 0.5563
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4353 0.6410 0.679 0.49712
## A 1.7076 0.1443 11.837 < 2e-16 ***
## C1 1.9982 0.6474 3.086 0.00203 **
## C2 0.9024 0.1345 6.709 1.96e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2294.0 on 9999 degrees of freedom
## Residual deviance: 2072.5 on 9996 degrees of freedom
## AIC: 2080.5
##
## Number of Fisher Scoring iterations: 7
##
##
## # Effect decomposition on the risk ratio scale via the regression-based approach
##
## Direct counterfactual imputation estimation with
## bootstrap standard errors, percentile confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## Rcde 1.32298 0.11254 1.15262 1.304 <2e-16 ***
## Rpnde 1.40305 0.18350 1.30134 1.548 <2e-16 ***
## Rtnde 1.34137 0.11990 1.18225 1.343 <2e-16 ***
## Rpnie 0.93320 0.03611 0.90219 0.951 <2e-16 ***
## Rtnie 0.89218 0.06009 0.78299 0.864 <2e-16 ***
## Rte 1.25177 0.06547 1.12398 1.212 <2e-16 ***
## ERcde 0.29546 0.09356 0.14337 0.269 <2e-16 ***
## ERintref 0.10759 0.08995 0.15853 0.279 <2e-16 ***
## ERintmed -0.08448 0.08192 -0.23863 -0.129 <2e-16 ***
## ERpnie -0.06680 0.03611 -0.09777 -0.049 <2e-16 ***
## ERcde(prop) 1.17354 0.08514 1.15331 1.268 <2e-16 ***
## ERintref(prop) 0.42734 0.02993 1.27696 1.317 <2e-16 ***
## ERintmed(prop) -0.33556 0.06693 -1.12447 -1.035 <2e-16 ***
## ERpnie(prop) -0.26532 0.04814 -0.46040 -0.396 <2e-16 ***
## pm -0.60088 0.11507 -1.58487 -1.430 <2e-16 ***
## int 0.09178 0.03700 0.19270 0.242 <2e-16 ***
## pe -0.17354 0.08514 -0.26770 -0.153 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rcde: controlled direct effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnde: total natural direct effect risk ratio; Rpnie: pure natural indirect effect risk ratio; Rtnie: total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; ERintref: excess relative risk due to reference interaction; ERintmed: excess relative risk due to mediated interaction; ERpnie: excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; ERintref(prop): proportion ERintref; ERintmed(prop): proportion ERintmed; ERpnie(prop): proportion ERpnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"
##
## $mval
## $mval[[1]]
## [1] 1The Weighting-based Approach
res_wb <- cmest(data = data, model = "wb", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
ereg = "logistic", yreg = "logistic",
astar = 0, a = 1, mval = list(1),
estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_wb)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ A + M + A * M + C1 + C2, family = binomial(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5681 -0.2414 -0.1741 -0.1604 3.0521
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5731 0.7713 -4.633 3.61e-06 ***
## A 0.7084 0.5283 1.341 0.17990
## M -1.0471 0.3736 -2.803 0.00507 **
## C1 0.9337 0.6973 1.339 0.18054
## C2 -0.8415 0.1443 -5.829 5.56e-09 ***
## A:M -0.4222 0.5541 -0.762 0.44607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2022.7 on 9999 degrees of freedom
## Residual deviance: 1965.0 on 9994 degrees of freedom
## AIC: 1977
##
## Number of Fisher Scoring iterations: 7
##
##
## # Exposure regression for weighting:
##
## Call:
## glm(formula = A ~ C1 + C2, family = binomial(), data = getCall(x$reg.output$ereg)$data,
## weights = getCall(x$reg.output$ereg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6126 -1.4791 0.8573 0.8867 0.9750
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.08342 0.21440 0.389 0.69723
## C1 0.60899 0.21208 2.872 0.00409 **
## C2 0.10532 0.04375 2.407 0.01606 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 12534 on 9999 degrees of freedom
## Residual deviance: 12520 on 9997 degrees of freedom
## AIC: 12526
##
## Number of Fisher Scoring iterations: 4
##
##
## # Effect decomposition on the risk ratio scale via the weighting-based approach
##
## Direct counterfactual imputation estimation with
## bootstrap standard errors, percentile confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## Rcde 1.322984 0.110877 0.994574 1.144 1
## Rpnde 1.415242 0.006241 1.203978 1.212 <2e-16 ***
## Rtnde 1.342561 0.085848 1.040273 1.156 <2e-16 ***
## Rpnie 0.922806 0.031552 0.940719 0.983 <2e-16 ***
## Rtnie 0.875415 0.044185 0.843561 0.903 <2e-16 ***
## Rte 1.238924 0.047933 1.022702 1.087 <2e-16 ***
## ERcde 0.291843 0.102790 -0.005249 0.133 1
## ERintref 0.123399 0.109031 0.071129 0.218 <2e-16 ***
## ERintmed -0.099124 0.085726 -0.172743 -0.058 <2e-16 ***
## ERpnie -0.077194 0.031552 -0.059256 -0.017 <2e-16 ***
## ERcde(prop) 1.221488 1.384502 -0.372872 1.487 1
## ERintref(prop) 0.516477 6.900553 1.001163 10.272 <2e-16 ***
## ERintmed(prop) -0.414878 5.467707 -8.153047 -0.807 <2e-16 ***
## ERpnie(prop) -0.323088 0.048344 -0.746158 -0.681 <2e-16 ***
## pm -0.737965 5.516051 -8.899204 -1.488 <2e-16 ***
## int 0.101600 1.432846 0.193996 2.119 <2e-16 ***
## pe -0.221488 1.384502 -0.487211 1.373 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rcde: controlled direct effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnde: total natural direct effect risk ratio; Rpnie: pure natural indirect effect risk ratio; Rtnie: total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; ERintref: excess relative risk due to reference interaction; ERintmed: excess relative risk due to mediated interaction; ERpnie: excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; ERintref(prop): proportion ERintref; ERintmed(prop): proportion ERintmed; ERpnie(prop): proportion ERpnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"
##
## $mval
## $mval[[1]]
## [1] 1The Inverse Odds-ratio Weighting Approach
res_iorw <- cmest(data = data, model = "iorw", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"),
ereg = "logistic", yreg = "logistic",
astar = 0, a = 1, mval = list(1),
estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_iorw)## Causal Mediation Analysis
##
## # Outcome regression for the total effect:
##
## Call:
## glm(formula = Y ~ A + C1 + C2, family = binomial(), data = getCall(x$reg.output$yregTot)$data,
## weights = getCall(x$reg.output$yregTot)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.3031 -0.2480 -0.1742 -0.1622 3.0225
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.4182 0.7137 -6.191 5.98e-10 ***
## A 0.2183 0.1565 1.395 0.163
## C1 0.8540 0.6956 1.228 0.220
## C2 -0.8831 0.1435 -6.156 7.47e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2022.7 on 9999 degrees of freedom
## Residual deviance: 1980.4 on 9996 degrees of freedom
## AIC: 1988.4
##
## Number of Fisher Scoring iterations: 7
##
##
## # Outcome regression for the direct effect:
##
## Call:
## glm(formula = Y ~ A + C1 + C2, family = binomial(), data = getCall(x$reg.output$yregDir)$data,
## weights = getCall(x$reg.output$yregDir)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4839 -0.2000 -0.1433 -0.1146 4.2345
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.7823 0.8563 -4.417 1.00e-05 ***
## A 0.3675 0.1735 2.118 0.0342 *
## C1 0.2839 0.8440 0.336 0.7366
## C2 -1.0638 0.1803 -5.901 3.61e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.8 on 9999 degrees of freedom
## Residual deviance: 1312.7 on 9996 degrees of freedom
## AIC: 798.25
##
## Number of Fisher Scoring iterations: 7
##
##
## # Exposure regression for weighting:
##
## Call:
## glm(formula = A ~ M + C1 + C2, family = binomial(), data = getCall(x$reg.output$ereg)$data,
## weights = getCall(x$reg.output$ereg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6138 -1.5035 0.8465 0.8684 1.6561
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4666 0.2542 -5.770 7.91e-09 ***
## M 1.7067 0.1442 11.832 < 2e-16 ***
## C1 0.5218 0.2140 2.438 0.0148 *
## C2 0.0648 0.0443 1.463 0.1436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 12534 on 9999 degrees of freedom
## Residual deviance: 12360 on 9996 degrees of freedom
## AIC: 12368
##
## Number of Fisher Scoring iterations: 4
##
##
## # Effect decomposition on the risk ratio scale via the inverse odds ratio weighting approach
##
## Direct counterfactual imputation estimation with
## bootstrap standard errors, percentile confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## Rte 1.23749 0.02336 1.43410 1.465 <2e-16 ***
## Rpnde 1.42948 0.04037 1.55304 1.607 <2e-16 ***
## Rtnie 0.86569 0.03823 0.89226 0.944 <2e-16 ***
## pm -0.80841 0.15694 -0.39937 -0.189 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rte: total effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnie: total natural indirect effect risk ratio; pm: proportion mediated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"The Natural Effect Model
res_ne <- cmest(data = data, model = "ne", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
yreg = "logistic",
astar = 0, a = 1, mval = list(1),
estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_ne)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ A + M + A * M + C1 + C2, family = binomial(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5681 -0.2414 -0.1741 -0.1604 3.0521
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5731 0.7713 -4.633 3.61e-06 ***
## A 0.7084 0.5283 1.341 0.17990
## M -1.0471 0.3736 -2.803 0.00507 **
## C1 0.9337 0.6973 1.339 0.18054
## C2 -0.8415 0.1443 -5.829 5.56e-09 ***
## A:M -0.4222 0.5541 -0.762 0.44607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2022.7 on 9999 degrees of freedom
## Residual deviance: 1965.0 on 9994 degrees of freedom
## AIC: 1977
##
## Number of Fisher Scoring iterations: 7
##
##
##
## # Effect decomposition on the risk ratio scale via the natural effect model
##
## Direct counterfactual imputation estimation with
## bootstrap standard errors, percentile confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## Rcde 1.32298 0.16150 1.45870 1.676 <2e-16 ***
## Rpnde 1.41685 0.08099 1.51129 1.620 <2e-16 ***
## Rtnde 1.34451 0.14151 1.47225 1.662 <2e-16 ***
## Rpnie 0.92215 0.01006 0.84766 0.861 <2e-16 ***
## Rtnie 0.87507 0.04308 0.82576 0.884 <2e-16 ***
## Rte 1.23984 0.13668 1.24796 1.432 <2e-16 ***
## ERcde 0.29180 0.13669 0.37185 0.555 <2e-16 ***
## ERintref 0.12505 0.05570 0.06471 0.140 <2e-16 ***
## ERintmed -0.09915 0.04562 -0.11076 -0.049 <2e-16 ***
## ERpnie -0.07785 0.01006 -0.15234 -0.139 <2e-16 ***
## ERcde(prop) 1.21661 0.15878 1.28856 1.502 <2e-16 ***
## ERintref(prop) 0.52138 0.30942 0.15459 0.570 <2e-16 ***
## ERintmed(prop) -0.41341 0.24890 -0.45276 -0.118 <2e-16 ***
## ERpnie(prop) -0.32459 0.21930 -0.61942 -0.325 <2e-16 ***
## pm -0.73799 0.46820 -1.07218 -0.443 <2e-16 ***
## int 0.10798 0.06052 0.03623 0.118 <2e-16 ***
## pe -0.21661 0.15878 -0.50188 -0.289 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rcde: controlled direct effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnde: total natural direct effect risk ratio; Rpnie: pure natural indirect effect risk ratio; Rtnie: total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; ERintref: excess relative risk due to reference interaction; ERintmed: excess relative risk due to mediated interaction; ERpnie: excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; ERintref(prop): proportion ERintref; ERintmed(prop): proportion ERintmed; ERpnie(prop): proportion ERpnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"
##
## $mval
## $mval[[1]]
## [1] 1The Marginal Structural Model
res_msm <- cmest(data = data, model = "msm", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
ereg = "logistic", yreg = "logistic", mreg = list("logistic"),
wmnomreg = list("logistic"), wmdenomreg = list("logistic"),
astar = 0, a = 1, mval = list(1),
estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_msm)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ A + M + A * M, family = binomial(), data = getCall(x$reg.output$yreg)$data,
## weights = getCall(x$reg.output$yreg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5182 -0.2088 -0.2057 -0.1827 3.4604
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.1163 0.3762 -8.284 <2e-16 ***
## A 0.4632 0.6200 0.747 0.4550
## M -0.9908 0.4028 -2.460 0.0139 *
## A:M -0.1803 0.6421 -0.281 0.7789
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1996.3 on 9999 degrees of freedom
## Residual deviance: 1985.5 on 9996 degrees of freedom
## AIC: 2014.4
##
## Number of Fisher Scoring iterations: 6
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ A, family = binomial(), data = getCall(x$reg.output$mreg[[1L]])$data,
## weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1423 0.1427 0.1446 0.3227 0.3593
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.87310 0.07858 36.56 <2e-16 ***
## A 1.69598 0.14370 11.80 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2271.0 on 9999 degrees of freedom
## Residual deviance: 2112.9 on 9998 degrees of freedom
## AIC: 2132.1
##
## Number of Fisher Scoring iterations: 7
##
##
## # Mediator regressions for weighting (denominator):
##
## Call:
## glm(formula = M ~ A + C1 + C2, family = binomial(), data = getCall(x$reg.output$wmdenomreg[[1L]])$data,
## weights = getCall(x$reg.output$wmdenomreg[[1L]])$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2578 0.1145 0.1647 0.2519 0.5563
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4353 0.6410 0.679 0.49712
## A 1.7076 0.1443 11.837 < 2e-16 ***
## C1 1.9982 0.6474 3.086 0.00203 **
## C2 0.9024 0.1345 6.709 1.96e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2294.0 on 9999 degrees of freedom
## Residual deviance: 2072.5 on 9996 degrees of freedom
## AIC: 2080.5
##
## Number of Fisher Scoring iterations: 7
##
##
## # Mediator regressions for weighting (nominator):
##
## Call:
## glm(formula = M ~ A, family = binomial(), data = getCall(x$reg.output$wmnomreg[[1L]])$data,
## weights = getCall(x$reg.output$wmnomreg[[1L]])$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0301 0.1428 0.1428 0.3355 0.3355
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.84922 0.07775 36.65 <2e-16 ***
## A 1.73145 0.14383 12.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2294 on 9999 degrees of freedom
## Residual deviance: 2128 on 9998 degrees of freedom
## AIC: 2132
##
## Number of Fisher Scoring iterations: 7
##
##
## # Exposure regression for weighting:
##
## Call:
## glm(formula = A ~ C1 + C2, family = binomial(), data = getCall(x$reg.output$ereg)$data,
## weights = getCall(x$reg.output$ereg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6126 -1.4791 0.8573 0.8867 0.9750
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.08342 0.21440 0.389 0.69723
## C1 0.60899 0.21208 2.872 0.00409 **
## C2 0.10532 0.04375 2.407 0.01606 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 12534 on 9999 degrees of freedom
## Residual deviance: 12520 on 9997 degrees of freedom
## AIC: 12526
##
## Number of Fisher Scoring iterations: 4
##
##
## # Effect decomposition on the risk ratio scale via the marginal structural model
##
## Direct counterfactual imputation estimation with
## bootstrap standard errors, percentile confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## Rcde 1.31995 0.31567 1.25898 1.683 <2e-16 ***
## Rpnde 1.34854 0.11471 1.49641 1.651 <2e-16 ***
## Rtnde 1.32626 0.28303 1.29556 1.676 <2e-16 ***
## Rpnie 0.93913 0.01814 0.92371 0.948 <2e-16 ***
## Rtnie 0.92361 0.08707 0.82083 0.938 <2e-16 ***
## Rte 1.24553 0.23792 1.22830 1.548 <2e-16 ***
## ERcde 0.29539 0.27915 0.24430 0.619 <2e-16 ***
## ERintref 0.05315 0.16444 0.03139 0.252 <2e-16 ***
## ERintmed -0.04214 0.14134 -0.21544 -0.026 <2e-16 ***
## ERpnie -0.06087 0.01814 -0.07628 -0.052 <2e-16 ***
## ERcde(prop) 1.20309 0.04775 1.06316 1.127 <2e-16 ***
## ERintref(prop) 0.21646 0.79400 0.07396 1.141 <2e-16 ***
## ERintmed(prop) -0.17165 0.67972 -0.97411 -0.061 <2e-16 ***
## ERpnie(prop) -0.24790 0.06654 -0.22976 -0.140 <2e-16 ***
## pm -0.41955 0.74625 -1.20387 -0.201 <2e-16 ***
## int 0.04481 0.11428 0.01306 0.167 <2e-16 ***
## pe -0.20309 0.04775 -0.12731 -0.063 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rcde: controlled direct effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnde: total natural direct effect risk ratio; Rpnie: pure natural indirect effect risk ratio; Rtnie: total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; ERintref: excess relative risk due to reference interaction; ERintmed: excess relative risk due to mediated interaction; ERpnie: excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; ERintref(prop): proportion ERintref; ERintmed(prop): proportion ERintmed; ERpnie(prop): proportion ERpnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"
##
## $mval
## $mval[[1]]
## [1] 1The g-formula Approach
res_gformula <- cmest(data = data, model = "gformula", outcome = "Y", exposure = "A",
mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
mreg = list("logistic"), yreg = "logistic",
astar = 0, a = 1, mval = list(1),
estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_gformula)## Causal Mediation Analysis
##
## # Outcome regression:
##
## Call:
## glm(formula = Y ~ A + M + A * M + C1 + C2, family = binomial(),
## data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5681 -0.2414 -0.1741 -0.1604 3.0521
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5731 0.7713 -4.633 3.61e-06 ***
## A 0.7084 0.5283 1.341 0.17990
## M -1.0471 0.3736 -2.803 0.00507 **
## C1 0.9337 0.6973 1.339 0.18054
## C2 -0.8415 0.1443 -5.829 5.56e-09 ***
## A:M -0.4222 0.5541 -0.762 0.44607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2022.7 on 9999 degrees of freedom
## Residual deviance: 1965.0 on 9994 degrees of freedom
## AIC: 1977
##
## Number of Fisher Scoring iterations: 7
##
##
## # Mediator regressions:
##
## Call:
## glm(formula = M ~ A + C1 + C2, family = binomial(), data = getCall(x$reg.output$mreg[[1L]])$data,
## weights = getCall(x$reg.output$mreg[[1L]])$weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2578 0.1145 0.1647 0.2519 0.5563
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4353 0.6410 0.679 0.49712
## A 1.7076 0.1443 11.837 < 2e-16 ***
## C1 1.9982 0.6474 3.086 0.00203 **
## C2 0.9024 0.1345 6.709 1.96e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2294.0 on 9999 degrees of freedom
## Residual deviance: 2072.5 on 9996 degrees of freedom
## AIC: 2080.5
##
## Number of Fisher Scoring iterations: 7
##
##
## # Effect decomposition on the risk ratio scale via the g-formula approach
##
## Direct counterfactual imputation estimation with
## bootstrap standard errors, percentile confidence intervals and p-values
##
## Estimate Std.error 95% CIL 95% CIU P.val
## Rcde 1.322984 0.178003 1.157574 1.397 <2e-16 ***
## Rpnde 1.408302 0.168338 1.324566 1.551 <2e-16 ***
## Rtnde 1.343111 0.179689 1.187007 1.428 <2e-16 ***
## Rpnie 0.929004 0.018661 0.934447 0.960 <2e-16 ***
## Rtnie 0.886001 0.034558 0.837403 0.884 <2e-16 ***
## Rte 1.247756 0.194578 1.109195 1.371 <2e-16 ***
## ERcde 0.293559 0.172309 0.145948 0.377 <2e-16 ***
## ERintref 0.114743 0.003972 0.173750 0.179 <2e-16 ***
## ERintmed -0.089550 0.007579 -0.149575 -0.139 <2e-16 ***
## ERpnie -0.070996 0.018661 -0.065544 -0.040 <2e-16 ***
## ERcde(prop) 1.184870 0.242721 1.022668 1.349 <2e-16 ***
## ERintref(prop) 0.463128 0.905395 0.491039 1.707 <2e-16 ***
## ERintmed(prop) -0.361443 0.768216 -1.427107 -0.395 <2e-16 ***
## ERpnie(prop) -0.286555 0.379901 -0.629098 -0.119 <2e-16 ***
## pm -0.647998 1.148117 -2.056205 -0.514 <2e-16 ***
## int 0.101685 0.137180 0.096032 0.280 <2e-16 ***
## pe -0.184870 0.242721 -0.348764 -0.023 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Rcde: controlled direct effect risk ratio; Rpnde: pure natural direct effect risk ratio; Rtnde: total natural direct effect risk ratio; Rpnie: pure natural indirect effect risk ratio; Rtnie: total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; ERintref: excess relative risk due to reference interaction; ERintmed: excess relative risk due to mediated interaction; ERpnie: excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; ERintref(prop): proportion ERintref; ERintmed(prop): proportion ERintmed; ERpnie(prop): proportion ERpnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
##
## Relevant variable values:
## $a
## [1] 1
##
## $astar
## [1] 0
##
## $yval
## [1] "1"
##
## $mval
## $mval[[1]]
## [1] 1