In LindaValeri/CMAverse: Causal Mediation Analysis

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)

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)

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)

The Natural Effect Model

```r

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)

```

```r

summary(res_ne)

```

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)

LindaValeri/CMAverse documentation built on July 16, 2024, 11:58 p.m.