This example demonstrates how to use cmest when there are multiple mediators. 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 count mediator \(M_1\), a categorical mediator \(M_2\) and a binary outcome \(Y\). The true regression models for \(A\), \(M_1\), \(M_2\) and \(Y\) are: \[logit(E(A|C_1,C_2))=0.2+0.5C_1+0.1C_2\] \[log(E(M_1|A,C_1,C_2))=1-2A+0.5C_1+0.8C_2\] \[log\frac{E[M_2=1|A,M_1,C_1,C_2]}{E[M_2=0|A,M_1,C_1,C_2]}=0.1+0.1A+0.4M_1-0.5C_1+0.1C_2\] \[log\frac{E[M_2=2|A,M_1,C_1,C_2]}{E[M_2=0|A,M_1,C_1,C_2]}=0.4+0.2A-0.1M_1-C_1+0.5C_2\] \[logit(E(Y|A,M_1,M_2,C_1,C_2)))=-4+0.8A-1.8M_1+0.5(M_2==1)+0.8(M_2==2)+0.5AM_1-0.4A(M_2==1)-1.4A(M_2==2)+0.3*C_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))
M1 <- rpois(n, exp(1 - 2*A + 0.5*C1 + 0.8*C2))
linpred1 <- 0.1 + 0.1*A + 0.4*M1 - 0.5*C1 + 0.1*C2
linpred2 <- 0.4 + 0.2*A - 0.1*M1 - C1 + 0.5*C2
probm0 = 1 / (1 + exp(linpred1) + exp(linpred2))
probm1 = exp(linpred1) / (1 + exp(linpred1) + exp(linpred2))
probm2 = exp(linpred2) / (1 + exp(linpred1) + exp(linpred2))
M2 = factor(sapply(1:n, FUN = function(x) sample(c(0, 1, 2), size = 1, replace = TRUE,
                                                 prob=c(probm0[x],
                                                        probm1[x],
                                                        probm2[x]))))
Y <- rbinom(n, 1, expit(1 + 0.8*A - 1.8*M1 + 0.5*(M2 == 1) + 0.8*(M2 == 2) + 
                          0.5*A*M1 - 0.4*A*(M2 == 1) - 1.4*A*(M2 == 2)  + 0.3*C1 - 0.6*C2))
data <- data.frame(A, M1, M2, Y, C1, C2)

The DAG for this scientific setting is:

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

In this setting, we have a count mediator, so the marginal structural model is not available. We can use the rest five statistical modeling approaches. The results are shown below.

The Regression-based Approach

res_rb <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
                mediator = c("M1", "M2"), basec = c("C1", "C2"), EMint = TRUE,
                mreg = list("poisson", "multinomial"), yreg = "logistic",
                astar = 0, a = 1, mval = list(0, 2), 
                estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_rb)
## Causal Mediation Analysis
## 
## # Outcome regression:
## 
## Call:
## glm(formula = Y ~ A + M1 + M2 + A * M1 + A * M2 + 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  
## -2.2066  -0.5055  -0.0007   0.6383   3.3715  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.125104   0.478893   2.349 0.018804 *  
## A            0.460265   0.400463   1.149 0.250419    
## M1          -1.809307   0.180027 -10.050  < 2e-16 ***
## M21         -0.007101   0.364067  -0.020 0.984438    
## M22          0.841167   0.405154   2.076 0.037879 *  
## C1           0.504243   0.284671   1.771 0.076508 .  
## C2          -0.643001   0.062025 -10.367  < 2e-16 ***
## A:M1         0.535812   0.183808   2.915 0.003556 ** 
## A:M21        0.129896   0.370967   0.350 0.726222    
## A:M22       -1.447096   0.412172  -3.511 0.000447 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13336.6  on 9999  degrees of freedom
## Residual deviance:  7347.9  on 9990  degrees of freedom
## AIC: 7367.9
## 
## Number of Fisher Scoring iterations: 10
## 
## 
## # Mediator regressions: 
## 
## Call:
## glm(formula = M1 ~ A + C1 + C2, family = poisson(), data = getCall(x$reg.output$mreg[[1L]])$data, 
##     weights = getCall(x$reg.output$mreg[[1L]])$weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4516  -1.0946  -0.2769   0.5079   3.4944  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)  1.06240    0.05631   18.868  < 2e-16 ***
## A           -2.00634    0.01331 -150.709  < 2e-16 ***
## C1           0.45102    0.05496    8.206 2.29e-16 ***
## C2           0.80084    0.01314   60.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 43205  on 9999  degrees of freedom
## Residual deviance: 10846  on 9996  degrees of freedom
## AIC: 33081
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## Call:
## nnet::multinom(formula = M2 ~ A + C1 + C2, data = getCall(x$reg.output$mreg[[2]])$data, 
##     trace = FALSE, weights = getCall(x$reg.output$mreg[[2]])$weights)
## 
## Coefficients:
##   (Intercept)         A         C1        C2
## 1    1.946698 -1.785773 -0.3220412 0.7256062
## 2    0.654218  0.688984 -1.7259629 0.4295399
## 
## Std. Errors:
##   (Intercept)          A        C1         C2
## 1   0.2620667 0.06393051 0.2550181 0.05235024
## 2   0.3129261 0.10149606 0.2996053 0.06106820
## 
## Residual Deviance: 17966.15 
## AIC: 17982.15 
## 
## # 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            0.840859  0.046685  0.788612   0.851 <2e-16 ***
## Rpnde           2.540010  0.249404  2.203176   2.538 <2e-16 ***
## Rtnde           1.198066  0.147841  1.024302   1.223 <2e-16 ***
## Rpnie          22.289482  0.880313 22.434122  23.617 <2e-16 ***
## Rtnie          10.513453  0.127523 10.808614  10.980 <2e-16 ***
## Rte            26.704271  2.414681 24.190746  27.435 <2e-16 ***
## ERcde          -6.683714  1.717727 -8.708034  -6.400 <2e-16 ***
## ERintref        8.223724  1.468323  7.939139   9.912 <2e-16 ***
## ERintmed        2.874780  3.045590 -0.625298   3.466      1    
## ERpnie         21.289482  0.880313 21.434936  22.618 <2e-16 ***
## ERcde(prop)    -0.260023  0.099297 -0.375913  -0.243 <2e-16 ***
## ERintref(prop)  0.319936  0.094601  0.300688   0.428 <2e-16 ***
## ERintmed(prop)  0.111841  0.117704 -0.027554   0.131      1    
## ERpnie(prop)    0.828247  0.122399  0.811238   0.976 <2e-16 ***
## pm              0.940087  0.004695  0.941819   0.948 <2e-16 ***
## int             0.431777  0.023102  0.400231   0.431 <2e-16 ***
## pe              1.260023  0.099297  1.242507   1.376 <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] 0
## 
## $mval[[2]]
## [1] 2

The Weighting-based Approach

res_wb <- cmest(data = data, model = "wb", outcome = "Y", exposure = "A",
                mediator = c("M1", "M2"), basec = c("C1", "C2"), EMint = TRUE,
                ereg = "logistic", yreg = "logistic",
                astar = 0, a = 1, mval = list(0, 2), 
                estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_wb)
## Causal Mediation Analysis
## 
## # Outcome regression:
## 
## Call:
## glm(formula = Y ~ A + M1 + M2 + A * M1 + A * M2 + 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  
## -2.2066  -0.5055  -0.0007   0.6383   3.3715  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.125104   0.478893   2.349 0.018804 *  
## A            0.460265   0.400463   1.149 0.250419    
## M1          -1.809307   0.180027 -10.050  < 2e-16 ***
## M21         -0.007101   0.364067  -0.020 0.984438    
## M22          0.841167   0.405154   2.076 0.037879 *  
## C1           0.504243   0.284671   1.771 0.076508 .  
## C2          -0.643001   0.062025 -10.367  < 2e-16 ***
## A:M1         0.535812   0.183808   2.915 0.003556 ** 
## A:M21        0.129896   0.370967   0.350 0.726222    
## A:M22       -1.447096   0.412172  -3.511 0.000447 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13336.6  on 9999  degrees of freedom
## Residual deviance:  7347.9  on 9990  degrees of freedom
## AIC: 7367.9
## 
## Number of Fisher Scoring iterations: 10
## 
## 
## # 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            0.840859  0.001153  0.843270   0.845 <2e-16 ***
## Rpnde           2.184496  0.124689  2.056035   2.224 <2e-16 ***
## Rtnde           1.197960  0.017578  1.135918   1.160 <2e-16 ***
## Rpnie          20.345349  1.376126 18.331408  20.180 <2e-16 ***
## Rtnie          11.157227  0.294594 10.127728  10.524 <2e-16 ***
## Rte            24.372911  1.917943 20.822968  23.400 <2e-16 ***
## ERcde          -6.150317  0.544402 -5.782384  -5.051 <2e-16 ***
## ERintref        7.334812  0.669092  6.107188   7.006 <2e-16 ***
## ERintmed        2.843066  0.417128  1.436970   1.997 <2e-16 ***
## ERpnie         19.345349  1.376126 17.333768  19.183 <2e-16 ***
## ERcde(prop)    -0.263139  0.002492 -0.258089  -0.255 <2e-16 ***
## ERintref(prop)  0.313817  0.003497  0.308008   0.313 <2e-16 ***
## ERintmed(prop)  0.121639  0.012419  0.072417   0.089 <2e-16 ***
## ERpnie(prop)    0.827682  0.013425  0.856280   0.874 <2e-16 ***
## pm              0.949322  0.001005  0.945382   0.947 <2e-16 ***
## int             0.435456  0.015916  0.380425   0.402 <2e-16 ***
## pe              1.263139  0.002492  1.254741   1.258 <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] 0
## 
## $mval[[2]]
## [1] 2

The Inverse Odds-ratio Weighting Approach

res_iorw <- cmest(data = data, model = "iorw", outcome = "Y", exposure = "A",
                  mediator = c("M1", "M2"), basec = c("C1", "C2"), EMint = TRUE,
                  ereg = "logistic", yreg = "logistic",
                  astar = 0, a = 1, mval = list(0, 2), 
                  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  
## -1.6731  -1.0673  -0.1524   0.7688   2.5278  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.9946     0.2726 -10.987   <2e-16 ***
## A             4.1997     0.1205  34.855   <2e-16 ***
## C1           -0.1306     0.2473  -0.528    0.597    
## C2           -1.3288     0.0534 -24.885   <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: 13336.6  on 9999  degrees of freedom
## Residual deviance:  9397.2  on 9996  degrees of freedom
## AIC: 9405.2
## 
## Number of Fisher Scoring iterations: 6
## 
## 
## # 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  
## -4.7789  -0.0549  -0.0069   0.0818   8.4799  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.5600     0.7405  -2.107   0.0351 *  
## A             1.7294     0.1460  11.848   <2e-16 ***
## C1           -1.2645     0.7413  -1.706   0.0880 .  
## C2           -3.7661     0.4333  -8.691   <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: 1926.6  on 9999  degrees of freedom
## Residual deviance: 1423.0  on 9996  degrees of freedom
## AIC: 1214.1
## 
## Number of Fisher Scoring iterations: 8
## 
## 
## # Exposure regression for weighting: 
## 
## Call:
## glm(formula = A ~ M1 + M2 + 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  
## -3.6520  -0.0024   0.0343   0.1447   3.1329  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.14163    0.58304   1.958   0.0502 .  
## M1          -2.00321    0.05823 -34.402  < 2e-16 ***
## M21          0.08020    0.13700   0.585   0.5583    
## M22         -0.01835    0.17781  -0.103   0.9178    
## C1           3.43446    0.58316   5.889 3.88e-09 ***
## C2           4.84059    0.18684  25.907  < 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: 12534.4  on 9999  degrees of freedom
## Residual deviance:  2139.2  on 9994  degrees of freedom
## AIC: 2151.2
## 
## Number of Fisher Scoring iterations: 8
## 
## 
## # 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   23.713831  1.745561 24.962887  27.308 <2e-16 ***
## Rpnde  4.483896  0.448867  4.168979   4.772 <2e-16 ***
## Rtnie  5.288667  0.197403  5.722561   5.988 <2e-16 ***
## pm     0.846618  0.008286  0.856622   0.868 <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 = c("M1", "M2"), basec = c("C1", "C2"), EMint = TRUE,
                yreg = "logistic",
                astar = 0, a = 1, mval = list(0, 2), 
                estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_ne)
## Causal Mediation Analysis
## 
## # Outcome regression:
## 
## Call:
## glm(formula = Y ~ A + M1 + M2 + A * M1 + A * M2 + 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  
## -2.2066  -0.5055  -0.0007   0.6383   3.3715  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.125104   0.478893   2.349 0.018804 *  
## A            0.460265   0.400463   1.149 0.250419    
## M1          -1.809307   0.180027 -10.050  < 2e-16 ***
## M21         -0.007101   0.364067  -0.020 0.984438    
## M22          0.841167   0.405154   2.076 0.037879 *  
## C1           0.504243   0.284671   1.771 0.076508 .  
## C2          -0.643001   0.062025 -10.367  < 2e-16 ***
## A:M1         0.535812   0.183808   2.915 0.003556 ** 
## A:M21        0.129896   0.370967   0.350 0.726222    
## A:M22       -1.447096   0.412172  -3.511 0.000447 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13336.6  on 9999  degrees of freedom
## Residual deviance:  7347.9  on 9990  degrees of freedom
## AIC: 7367.9
## 
## Number of Fisher Scoring iterations: 10
## 
## 
## 
## # 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            0.840859  0.086137  0.895442   1.011      1    
## Rpnde           2.397002  0.104623  2.429955   2.571 <2e-16 ***
## Rtnde           1.212344  0.101779  1.330542   1.467 <2e-16 ***
## Rpnie          20.806505  1.917377 18.176285  20.752 <2e-16 ***
## Rtnie          10.523413  0.173974 10.741630  10.975 <2e-16 ***
## Rte            25.224641  0.701059 26.669640  27.612 <2e-16 ***
## ERcde          -6.332043  3.477253 -4.288821   0.383      1    
## ERintref        7.729045  3.581876  1.047182   5.859 <2e-16 ***
## ERintmed        3.021133  1.320940  5.284538   7.059 <2e-16 ***
## ERpnie         19.806505  1.917377 17.180793  19.757 <2e-16 ***
## ERcde(prop)    -0.261389  0.131063 -0.160993   0.015      1    
## ERintref(prop)  0.319057  0.133527  0.040616   0.220 <2e-16 ***
## ERintmed(prop)  0.124713  0.056884  0.198651   0.275 <2e-16 ***
## ERpnie(prop)    0.817618  0.054420  0.669220   0.742 <2e-16 ***
## pm              0.942331  0.002464  0.940984   0.944 <2e-16 ***
## int             0.443770  0.076643  0.315690   0.419 <2e-16 ***
## pe              1.261389  0.131063  0.984910   1.161 <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] 0
## 
## $mval[[2]]
## [1] 2

The g-formula Approach

res_gformula <- cmest(data = data, model = "gformula", outcome = "Y", exposure = "A",
                      mediator = c("M1", "M2"), basec = c("C1", "C2"), EMint = TRUE,
                      mreg = list("poisson", "multinomial"), yreg = "logistic",
                      astar = 0, a = 1, mval = list(0, 2), 
                      estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(res_gformula)
## Causal Mediation Analysis
## 
## # Outcome regression:
## 
## Call:
## glm(formula = Y ~ A + M1 + M2 + A * M1 + A * M2 + 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  
## -2.2066  -0.5055  -0.0007   0.6383   3.3715  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.125104   0.478893   2.349 0.018804 *  
## A            0.460265   0.400463   1.149 0.250419    
## M1          -1.809307   0.180027 -10.050  < 2e-16 ***
## M21         -0.007101   0.364067  -0.020 0.984438    
## M22          0.841167   0.405154   2.076 0.037879 *  
## C1           0.504243   0.284671   1.771 0.076508 .  
## C2          -0.643001   0.062025 -10.367  < 2e-16 ***
## A:M1         0.535812   0.183808   2.915 0.003556 ** 
## A:M21        0.129896   0.370967   0.350 0.726222    
## A:M22       -1.447096   0.412172  -3.511 0.000447 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13336.6  on 9999  degrees of freedom
## Residual deviance:  7347.9  on 9990  degrees of freedom
## AIC: 7367.9
## 
## Number of Fisher Scoring iterations: 10
## 
## 
## # Mediator regressions: 
## 
## Call:
## glm(formula = M1 ~ A + C1 + C2, family = poisson(), data = getCall(x$reg.output$mreg[[1L]])$data, 
##     weights = getCall(x$reg.output$mreg[[1L]])$weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4516  -1.0946  -0.2769   0.5079   3.4944  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)  1.06240    0.05631   18.868  < 2e-16 ***
## A           -2.00634    0.01331 -150.709  < 2e-16 ***
## C1           0.45102    0.05496    8.206 2.29e-16 ***
## C2           0.80084    0.01314   60.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 43205  on 9999  degrees of freedom
## Residual deviance: 10846  on 9996  degrees of freedom
## AIC: 33081
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## Call:
## nnet::multinom(formula = M2 ~ A + C1 + C2, data = getCall(x$reg.output$mreg[[2]])$data, 
##     trace = FALSE, weights = getCall(x$reg.output$mreg[[2]])$weights)
## 
## Coefficients:
##   (Intercept)         A         C1        C2
## 1    1.946698 -1.785773 -0.3220412 0.7256062
## 2    0.654218  0.688984 -1.7259629 0.4295399
## 
## Std. Errors:
##   (Intercept)          A        C1         C2
## 1   0.2620667 0.06393051 0.2550181 0.05235024
## 2   0.3129261 0.10149606 0.2996053 0.06106820
## 
## Residual Deviance: 17966.15 
## AIC: 17982.15 
## 
## # 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.840859  0.065592  0.771998   0.860 <2e-16 ***
## Rpnde           2.628791  0.301244  2.110514   2.515 <2e-16 ***
## Rtnde           1.206834  0.190209  1.003274   1.259      1    
## Rpnie          23.580908  1.786324 22.156678  24.556 <2e-16 ***
## Rtnie          10.825600  0.435251 11.088659  11.673 <2e-16 ***
## Rte            28.458238  2.421510 24.636897  27.890 <2e-16 ***
## ERcde          -7.113482  2.684539 -9.594347  -5.988 <2e-16 ***
## ERintref        8.742273  2.383294  7.503828  10.706 <2e-16 ***
## ERintmed        3.248539  3.906589 -1.029096   4.219      1    
## ERpnie         22.580908  1.786324 21.159954  23.560 <2e-16 ***
## ERcde(prop)    -0.259065  0.136412 -0.406495  -0.223 <2e-16 ***
## ERintref(prop)  0.318384  0.129441  0.279567   0.453 <2e-16 ***
## ERintmed(prop)  0.118308  0.149237 -0.044273   0.156      1    
## ERpnie(prop)    0.822373  0.156208  0.787431   0.997 <2e-16 ***
## pm              0.940681  0.006971  0.943658   0.953 <2e-16 ***
## int             0.436693  0.019796  0.409198   0.436 <2e-16 ***
## pe              1.259065  0.136412  1.223225   1.406 <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] 0
## 
## $mval[[2]]
## [1] 2