Since there is no gold standard to verify our extended formulas
incorporating effect modifications, we use bootstrap to check the point
estimates and standard errors given in one single run of
regmedint(). The original estimation for standard errors
uses delta method, which agrees with bootstrap asymptotically.
We use simulated data, and present the estimates from 5000 times of boostrap. For the purpose of demonstration, we include all three effect modification terms (i.e. \(A\times C\) in mediator model, \(A\times C\) in outcome model, and \(M\times C\) in outcome model). Due to the long computational time, the code chunks are commented out and only the summary tables are shown. Readers are free to run the code on their side to replicate the results.
## [1] 12# Model 1: M linear, Y linear
datamaker.s4.m1 <- function(n, k){
C <- matrix(rnorm(n*1, 0, 2), ncol = 1)
A <- rbinom(n, 1, expit(C + C^2))
M <- 0.2 + 0.4*A + 0.5*C + 0.2*A*C + rnorm(n, 0, 0.5)
Y <- 0.5 + 0.3*A + 0.2*M + 0.5*A*M + 0.2*A*C + k*M*C + rnorm(n, 0, 0.5)
list(C = C, A = A, M = M, Y = Y)
}
# Model 2: M logistic, Y linear
datamaker.s4.m2 <- function(n, k){
C <- matrix(rnorm(n*1, 0, 2), ncol = 1)
A <- rbinom(n, 1, expit(C + C^2))
M <- rbinom(n, 1, expit(0.2 + 0.4*A + 0.5*C + 0.2*A*C))
Y <- 0.5 + 0.3*A + 0.2*M + 0.5*A*M + 0.2*A*C + k*M*C + rnorm(n, 0, 0.5)
list(C = C, A = A, M = M, Y = Y)
}
# Model 3: M linear, Y logistic
datamaker.s4.m3 <- function(n, k){
C <- matrix(rnorm(n*1, 0, 2), ncol = 1)
A <- rbinom(n, 1, expit(C + C^2))
M <- (0.2 + 0.4*A + 0.5*C + 0.2*A*C + rnorm(n, 0, 0.5)/5)
Y <- rbinom(n, 1, expit((0.5 + 0.3*A + 0.6*M + 0.4*C + 0.5*A*M + 0.2*A*C + k*M*C)))
list(C = C, A = A, M = M, Y = Y)
}
# Model 4: M logistic, Y logistic
datamaker.s4.m4 <- function(n, k){
C <- matrix(rnorm(n*1, 0, 2), ncol = 1)
A <- rbinom(n, 1, expit(C + C^2))
M <- rbinom(n, 1, expit(0.2 + 0.4*A + 0.5*C + 0.2*A*C))
Y <- rbinom(n, 1, expit(0.5 + 0.3*A + 0.2*M + 0.1*C + 0.5*A*M + 0.2*A*C + k*M*C))
list(C = C, A = A, M = M, Y = Y)
}set.seed(seed)
dat_linear_M_linear_Y <- as.data.frame(datamaker.s4.m1(n = 5000, k = 0.3))
dat_logistic_M_linear_Y <- as.data.frame(datamaker.s4.m2(n = 5000, k = 0.3))
dat_linear_M_logistic_Y <- as.data.frame(datamaker.s4.m3(n = 5000, k = 0.7))
dat_logistic_M_logistic_Y <- as.data.frame(datamaker.s4.m4(n = 5000, k = 0.3))regmedint1 <- regmedint(data = dat_linear_M_linear_Y,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0.5012509,
c_cond = -0.0434094,
mreg = "linear",
yreg = "linear",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
summary(regmedint1)data1 <- dat_linear_M_linear_Y
boot1 <- function(trials){
ind <- sample(5000, 5000, replace = TRUE)
dat <- data1[ind,]
regmedint1 <- regmedint(data = dat,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0.5012509,
c_cond = -0.0434094,
mreg = "linear",
yreg = "linear",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
out <- summary(regmedint1)
cde.est.boot <- out$summary_myreg[1,1]
pnde.est.boot <- out$summary_myreg[2,1]
tnie.est.boot <- out$summary_myreg[3,1]
tnde.est.boot <- out$summary_myreg[4,1]
pnie.est.boot <- out$summary_myreg[5,1]
te.est.boot <- out$summary_myreg[6,1]
pm.est.boot <- out$summary_myreg[7,1]
return(c(cde.est.boot,
pnde.est.boot, tnie.est.boot,
tnde.est.boot, pnie.est.boot,
te.est.boot, pm.est.boot))
}
set.seed(seed)
system.time({
results1 <- mclapply(trials, boot1, mc.cores = numCores)
})
results1.df <- as.data.frame(do.call(rbind, results1))
apply(results1.df, 2, mean)
apply(results1.df, 2, sd)regmedint2 <- regmedint(data = dat_logistic_M_linear_Y,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0,
c_cond = -0.0434094,
mreg = "logistic",
yreg = "linear",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
summary(regmedint2)data2 <- dat_logistic_M_linear_Y
boot2 <- function(trials){
ind <- sample(5000, 5000, replace = TRUE)
dat <- data2[ind,]
regmedint2 <- regmedint(data = dat,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0,
c_cond = -0.0434094,
mreg = "logistic",
yreg = "linear",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
out <- summary(regmedint2)
cde.est.boot <- out$summary_myreg[1,1]
pnde.est.boot <- out$summary_myreg[2,1]
tnie.est.boot <- out$summary_myreg[3,1]
tnde.est.boot <- out$summary_myreg[4,1]
pnie.est.boot <- out$summary_myreg[5,1]
te.est.boot <- out$summary_myreg[6,1]
pm.est.boot <- out$summary_myreg[7,1]
return(c(cde.est.boot,
pnde.est.boot, tnie.est.boot,
tnde.est.boot, pnie.est.boot,
te.est.boot, pm.est.boot))
}
set.seed(seed)
system.time({
results2 <- mclapply(1:100, boot2, mc.cores = numCores)
})
results2.df <- as.data.frame(do.call(rbind, results2))
apply(results2.df, 2, mean)
apply(results2.df, 2, sd)regmedint3 <- regmedint(data = dat_linear_M_logistic_Y,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0.5012509,
c_cond = 0.5,
mreg = "linear",
yreg = "logistic",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
summary(regmedint3)data3 <- dat_linear_M_logistic_Y
boot3 <- function(trials){
ind <- sample(5000, 5000, replace = TRUE)
dat <- data3[ind,]
regmedint3 <- regmedint(data = dat,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0.5012509,
c_cond = 0.5,
mreg = "linear",
yreg = "logistic",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
out <- summary(regmedint3)
cde.est.boot <- out$summary_myreg[1,1]
pnde.est.boot <- out$summary_myreg[2,1]
tnie.est.boot <- out$summary_myreg[3,1]
tnde.est.boot <- out$summary_myreg[4,1]
pnie.est.boot <- out$summary_myreg[5,1]
te.est.boot <- out$summary_myreg[6,1]
pm.est.boot <- out$summary_myreg[7,1]
return(c(cde.est.boot,
pnde.est.boot, tnie.est.boot,
tnde.est.boot, pnie.est.boot,
te.est.boot, pm.est.boot))
}
set.seed(seed)
system.time({
results3 <- mclapply(trials, boot3, mc.cores = numCores)
})
results3.df <- as.data.frame(do.call(rbind, results3))
apply(results3.df, 2, mean)
apply(results3.df, 2, sd)regmedint4 <- regmedint(data = dat_logistic_M_logistic_Y,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0,
c_cond = -0.0434094,
mreg = "logistic",
yreg = "logistic",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
summary(regmedint4)data4 <- dat_logistic_M_logistic_Y
boot4 <- function(trials){
ind <- sample(5000, 5000, replace = TRUE)
dat <- data4[ind,]
regmedint4 <- regmedint(data = dat,
yvar = "Y",
avar = "A",
mvar = "M",
cvar = c("C"),
emm_ac_mreg = c("C"),
emm_ac_yreg = c("C"),
emm_mc_yreg = c("C"),
eventvar = NULL,
a0 = 0,
a1 = 1,
m_cde = 0,
c_cond = -0.0434094,
mreg = "logistic",
yreg = "logistic",
interaction = TRUE,
casecontrol = FALSE,
na_omit = FALSE)
out <- summary(regmedint4)
cde.est.boot <- out$summary_myreg[1,1]
pnde.est.boot <- out$summary_myreg[2,1]
tnie.est.boot <- out$summary_myreg[3,1]
tnde.est.boot <- out$summary_myreg[4,1]
pnie.est.boot <- out$summary_myreg[5,1]
te.est.boot <- out$summary_myreg[6,1]
pm.est.boot <- out$summary_myreg[7,1]
return(c(cde.est.boot,
pnde.est.boot, tnie.est.boot,
tnde.est.boot, pnie.est.boot,
te.est.boot, pm.est.boot))
}
set.seed(seed)
system.time({
results4 <- mclapply(trials, boot4, mc.cores = numCores)
})
results4.df <- as.data.frame(do.call(rbind, results4))
apply(results4.df, 2, mean)
apply(results4.df, 2, sd)The following tables shows the point estimates and standard errors
from one single run of regmedint() and bootstrap.
## Error in library(kableExtra): there is no package called 'kableExtra'## Error in library(formattable): there is no package called 'formattable'## Error in add_header_above(., c(` ` = 1, `Non-bootstrap` = 2, Bootstrap = 2)): could not find function "add_header_above"## Error in add_header_above(., c(` ` = 1, `Non-bootstrap` = 2, Bootstrap = 2)): could not find function "add_header_above"## Error in add_header_above(., c(` ` = 1, `Non-bootstrap` = 2, Bootstrap = 2)): could not find function "add_header_above"## Error in add_header_above(., c(` ` = 1, `Non-bootstrap` = 2, Bootstrap = 2)): could not find function "add_header_above"