This 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\]

library(CMAverse)
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.6061

Closed-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.6061

Direct 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] 1

The 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] 1

The 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] 1

The 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] 1

The 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