This example demonstrates how to use cmest when there are mediator-outcome confounders affected by the exposure. 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 continuous mediator-outcome confounder affected by the exposure \(L\), a binary mediator \(M\) and a binary outcome \(Y\). The true regression models for \(A\), \(L\), \(M\) and \(Y\) are: \[logit(E(A|C_1,C_2))=0.2+0.5C_1+0.1C_2\] \[E(L|A,C_1,C_2)=1+A-C_1-0.5C_2\] \[logit(E(M|A,L,C_1,C_2))=1+2A-L+1.5C_1+0.8C_2\] \[logit(E(Y|A,L,M,C_1,C_2)))=-3-0.4A-1.2M+0.5AM-0.5L+0.3C_1-0.6C_2\]

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))
L <- rnorm(n, mean = 1 + A - C1 - 0.5*C2, sd = 0.5)
M <- rbinom(n, 1, expit(1 + 2*A - L + 1.5*C1 + 0.8*C2))
Y <- rbinom(n, 1, expit(-3 - 0.4*A - 1.2*M + 0.5*A*M - 0.5*L + 0.3*C1 - 0.6*C2))
data <- data.frame(A, M, Y, C1, C2, L)

The DAG for this scientific setting is:

library(CMAverse)
cmdag(outcome = "Y", exposure = "A", mediator = "M",
      basec = c("C1", "C2"), postc = "L", node = TRUE, text_col = "white")

In this setting, we can use the marginal structural model and the \(g\)-formula approach. The results are shown below.

The Marginal Structural Model

res_msm <- cmest(data = data, model = "msm", outcome = "Y", exposure = "A",
                 mediator = "M", basec = c("C1", "C2"), postc = "L", 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.6465  -0.1964  -0.1428  -0.1399   3.3658  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.49700    0.47457  -7.369 1.72e-13 ***
## A            0.09466    0.67862   0.139    0.889    
## M           -0.40808    0.49218  -0.829    0.407    
## A:M         -0.79067    0.70198  -1.126    0.260    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1432.5  on 9999  degrees of freedom
## Residual deviance: 1412.9  on 9996  degrees of freedom
## AIC: 1432.5
## 
## Number of Fisher Scoring iterations: 7
## 
## 
## # 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  
## -2.9114   0.1972   0.1999   0.3089   0.3437  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.96393    0.08185  36.211  < 2e-16 ***
## A            0.95266    0.11992   7.944 1.95e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2623.2  on 9999  degrees of freedom
## Residual deviance: 2560.6  on 9998  degrees of freedom
## AIC: 2580.5
## 
## Number of Fisher Scoring iterations: 6
## 
## 
## # Mediator regressions for weighting (denominator): 
## 
## Call:
## glm(formula = M ~ A + C1 + C2 + L, 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.4241   0.1334   0.1885   0.2668   0.8260  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.9946     0.6067   1.639   0.1012    
## A             1.8884     0.1726  10.944  < 2e-16 ***
## C1            1.4955     0.6090   2.456   0.0141 *  
## C2            0.8166     0.1410   5.793 6.92e-09 ***
## L            -0.9171     0.1215  -7.546 4.50e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2646.0  on 9999  degrees of freedom
## Residual deviance: 2397.3  on 9995  degrees of freedom
## AIC: 2407.3
## 
## 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  
## -2.8106   0.1972   0.1972   0.3224   0.3224  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.93070    0.08064  36.345   <2e-16 ***
## A            0.99963    0.11952   8.364   <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: 2646.0  on 9999  degrees of freedom
## Residual deviance: 2576.3  on 9998  degrees of freedom
## AIC: 2580.3
## 
## Number of Fisher Scoring iterations: 6
## 
## 
## # 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             0.5035510  0.1696963  0.3942646   0.622 <2e-16 ***
## rRpnde           0.5444178  0.1812002  0.4051440   0.648 <2e-16 ***
## rRtnde           0.5208059  0.1743300  0.3991930   0.633 <2e-16 ***
## rRpnie           0.9867847  0.0087895  1.0023348   1.014 <2e-16 ***
## rRtnie           0.9439869  0.0151504  0.9788925   0.999 <2e-16 ***
## Rte              0.5139232  0.1712168  0.4048390   0.635 <2e-16 ***
## ERcde           -0.4852025  0.1803464 -0.6204261  -0.378 <2e-16 ***
## rERintref        0.0296203  0.0008538  0.0269944   0.028 <2e-16 ***
## rERintmed       -0.0172792  0.0011939 -0.0161869  -0.015 <2e-16 ***
## rERpnie         -0.0132153  0.0087895  0.0023366   0.014 <2e-16 ***
## ERcde(prop)      0.9982014  0.0040755  1.0391856   1.045 <2e-16 ***
## rERintref(prop) -0.0609374  0.0238952 -0.0778943  -0.046 <2e-16 ***
## rERintmed(prop)  0.0355483  0.0149357  0.0247657   0.045 <2e-16 ***
## rERpnie(prop)    0.0271878  0.0130350 -0.0236356  -0.006 <2e-16 ***
## rpm              0.0627361  0.0279707  0.0011300   0.039 <2e-16 ***
## rint            -0.0253891  0.0089595 -0.0330625  -0.021 <2e-16 ***
## rpe              0.0017986  0.0040755 -0.0446609  -0.039 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Rcde: controlled direct effect risk ratio; rRpnde: randomized analogue of pure natural direct effect risk ratio; rRtnde: randomized analogue of total natural direct effect risk ratio; rRpnie: randomized analogue of pure natural indirect effect risk ratio; rRtnie: randomized analogue of total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; rERintref: randomized analogue of excess relative risk due to reference interaction; rERintmed: randomized analogue of excess relative risk due to mediated interaction; rERpnie: randomized analogue of excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; rERintref(prop): proportion rERintref; rERintmed(prop): proportion rERintmed; rERpnie(prop): proportion rERpnie; rpm: randomized analogue of overall proportion mediated; rint: randomized analogue of overall proportion attributable to interaction; rpe: randomized analogue of 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"), postc = "L", EMint = TRUE,
                      mreg = list("logistic"), yreg = "logistic", postcreg = list("linear"),
                      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 + L, family = binomial(), 
##     data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.3617  -0.1788  -0.1545  -0.1289   3.1833  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -2.7488     0.9365  -2.935  0.00333 **
## A            -0.1267     0.6617  -0.192  0.84812   
## M            -0.7274     0.4136  -1.759  0.07863 . 
## C1           -0.1906     0.8705  -0.219  0.82664   
## C2           -0.5566     0.1933  -2.880  0.00397 **
## L            -0.2242     0.1709  -1.312  0.18964   
## A:M          -0.3425     0.6635  -0.516  0.60572   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1447.7  on 9999  degrees of freedom
## Residual deviance: 1415.5  on 9993  degrees of freedom
## AIC: 1429.5
## 
## Number of Fisher Scoring iterations: 7
## 
## 
## # Mediator regressions: 
## 
## Call:
## glm(formula = M ~ A + C1 + C2 + L, 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.4241   0.1334   0.1885   0.2668   0.8260  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.9946     0.6067   1.639   0.1012    
## A             1.8884     0.1726  10.944  < 2e-16 ***
## C1            1.4955     0.6090   2.456   0.0141 *  
## C2            0.8166     0.1410   5.793 6.92e-09 ***
## L            -0.9171     0.1215  -7.546 4.50e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2646.0  on 9999  degrees of freedom
## Residual deviance: 2397.3  on 9995  degrees of freedom
## AIC: 2407.3
## 
## Number of Fisher Scoring iterations: 7
## 
## 
## # Regressions for mediator-outcome confounders affected by the exposure: 
## 
## Call:
## glm(formula = L ~ A + C1 + C2, family = gaussian(), data = getCall(x$reg.output$postcreg[[1L]])$data, 
##     weights = getCall(x$reg.output$postcreg[[1L]])$weights)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.82456  -0.34194   0.00712   0.33751   1.82462  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.00268    0.05077   19.75   <2e-16 ***
## A            1.00202    0.01081   92.68   <2e-16 ***
## C1          -1.00369    0.04980  -20.15   <2e-16 ***
## C2          -0.49437    0.01032  -47.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2539334)
## 
##     Null deviance: 5324.2  on 9999  degrees of freedom
## Residual deviance: 2538.3  on 9996  degrees of freedom
## AIC: 14678
## 
## Number of Fisher Scoring iterations: 2
## 
## 
## # 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             0.504938  0.003275  0.631908   0.636 <2e-16 ***
## rRpnde           0.524665  0.062186  0.545245   0.629 <2e-16 ***
## rRtnde           0.513482  0.026241  0.597452   0.633 <2e-16 ***
## rRpnie           0.972198  0.029936  0.941736   0.982 <2e-16 ***
## rRtnie           0.951478  0.032623  0.988079   1.032      1    
## Rte              0.499207  0.043654  0.562642   0.621 <2e-16 ***
## ERcde           -0.471286  0.012399 -0.350420  -0.334 <2e-16 ***
## rERintref       -0.004050  0.074586 -0.120836  -0.021 <2e-16 ***
## rERintmed        0.002344  0.048468  0.010442   0.076 <2e-16 ***
## rERpnie         -0.027802  0.029936 -0.058241  -0.018 <2e-16 ***
## ERcde(prop)      0.941078  0.120817  0.763855   0.926 <2e-16 ***
## rERintref(prop)  0.008087  0.165221  0.053557   0.276 <2e-16 ***
## rERintmed(prop) -0.004680  0.108147 -0.172266  -0.027 <2e-16 ***
## rERpnie(prop)    0.055515  0.063743  0.047240   0.133 <2e-16 ***
## rpm              0.050835  0.044404 -0.039387   0.020      1    
## rint             0.003407  0.057073  0.026588   0.103 <2e-16 ***
## rpe              0.058922  0.120817  0.073828   0.236 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Rcde: controlled direct effect risk ratio; rRpnde: randomized analogue of pure natural direct effect risk ratio; rRtnde: randomized analogue of total natural direct effect risk ratio; rRpnie: randomized analogue of pure natural indirect effect risk ratio; rRtnie: randomized analogue of total natural indirect effect risk ratio; Rte: total effect risk ratio; ERcde: excess relative risk due to controlled direct effect; rERintref: randomized analogue of excess relative risk due to reference interaction; rERintmed: randomized analogue of excess relative risk due to mediated interaction; rERpnie: randomized analogue of excess relative risk due to pure natural indirect effect; ERcde(prop): proportion ERcde; rERintref(prop): proportion rERintref; rERintmed(prop): proportion rERintmed; rERpnie(prop): proportion rERpnie; rpm: randomized analogue of overall proportion mediated; rint: randomized analogue of overall proportion attributable to interaction; rpe: randomized analogue of overall proportion eliminated)
## 
## Relevant variable values: 
## $a
## [1] 1
## 
## $astar
## [1] 0
## 
## $yval
## [1] "1"
## 
## $mval
## $mval[[1]]
## [1] 1