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.
library(regmedint)
library(tidyverse)
<- function(x){exp(x)/(1 + exp(x))} expit
library(parallel)
library(MASS)
<- detectCores()
numCores numCores
## [1] 8
# Number of bootstrap
<- 1:5000
trials <- 3104 seed
# Model 1: M linear, Y linear
<- function(n, k){
datamaker.s4.m1 <- matrix(rnorm(n*1, 0, 2), ncol = 1)
C <- rbinom(n, 1, expit(C + C^2))
A <- 0.2 + 0.4*A + 0.5*C + 0.2*A*C + rnorm(n, 0, 0.5)
M <- 0.5 + 0.3*A + 0.2*M + 0.5*A*M + 0.2*A*C + k*M*C + rnorm(n, 0, 0.5)
Y list(C = C, A = A, M = M, Y = Y)
}
# Model 2: M logistic, Y linear
<- function(n, k){
datamaker.s4.m2 <- matrix(rnorm(n*1, 0, 2), ncol = 1)
C <- rbinom(n, 1, expit(C + C^2))
A <- rbinom(n, 1, expit(0.2 + 0.4*A + 0.5*C + 0.2*A*C))
M <- 0.5 + 0.3*A + 0.2*M + 0.5*A*M + 0.2*A*C + k*M*C + rnorm(n, 0, 0.5)
Y list(C = C, A = A, M = M, Y = Y)
}
# Model 3: M linear, Y logistic
<- function(n, k){
datamaker.s4.m3 <- matrix(rnorm(n*1, 0, 2), ncol = 1)
C <- rbinom(n, 1, expit(C + C^2))
A <- (0.2 + 0.4*A + 0.5*C + 0.2*A*C + rnorm(n, 0, 0.5)/5)
M <- 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)))
Y list(C = C, A = A, M = M, Y = Y)
}
# Model 4: M logistic, Y logistic
<- function(n, k){
datamaker.s4.m4 <- matrix(rnorm(n*1, 0, 2), ncol = 1)
C <- rbinom(n, 1, expit(C + C^2))
A <- rbinom(n, 1, expit(0.2 + 0.4*A + 0.5*C + 0.2*A*C))
M <- 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))
Y list(C = C, A = A, M = M, Y = Y)
}
set.seed(seed)
<- as.data.frame(datamaker.s4.m1(n = 5000, k = 0.3))
dat_linear_M_linear_Y <- as.data.frame(datamaker.s4.m2(n = 5000, k = 0.3))
dat_logistic_M_linear_Y <- as.data.frame(datamaker.s4.m3(n = 5000, k = 0.7))
dat_linear_M_logistic_Y <- as.data.frame(datamaker.s4.m4(n = 5000, k = 0.3)) dat_logistic_M_logistic_Y
<- regmedint(data = dat_linear_M_linear_Y,
regmedint1 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)
<- dat_linear_M_linear_Y
data1 <- function(trials){
boot1 <- sample(5000, 5000, replace = TRUE)
ind <- data1[ind,]
dat
<- regmedint(data = dat,
regmedint1 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)
out <- out$summary_myreg[1,1]
cde.est.boot <- out$summary_myreg[2,1]
pnde.est.boot <- out$summary_myreg[3,1]
tnie.est.boot <- out$summary_myreg[4,1]
tnde.est.boot <- out$summary_myreg[5,1]
pnie.est.boot <- out$summary_myreg[6,1]
te.est.boot <- out$summary_myreg[7,1]
pm.est.boot 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({
<- mclapply(trials, boot1, mc.cores = numCores)
results1
})
<- as.data.frame(do.call(rbind, results1))
results1.df apply(results1.df, 2, mean)
apply(results1.df, 2, sd)
<- regmedint(data = dat_logistic_M_linear_Y,
regmedint2 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)
<- dat_logistic_M_linear_Y
data2 <- function(trials){
boot2 <- sample(5000, 5000, replace = TRUE)
ind <- data2[ind,]
dat
<- regmedint(data = dat,
regmedint2 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)
out <- out$summary_myreg[1,1]
cde.est.boot <- out$summary_myreg[2,1]
pnde.est.boot <- out$summary_myreg[3,1]
tnie.est.boot <- out$summary_myreg[4,1]
tnde.est.boot <- out$summary_myreg[5,1]
pnie.est.boot <- out$summary_myreg[6,1]
te.est.boot <- out$summary_myreg[7,1]
pm.est.boot 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({
<- mclapply(1:100, boot2, mc.cores = numCores)
results2
})
<- as.data.frame(do.call(rbind, results2))
results2.df apply(results2.df, 2, mean)
apply(results2.df, 2, sd)
<- regmedint(data = dat_linear_M_logistic_Y,
regmedint3 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)
<- dat_linear_M_logistic_Y
data3 <- function(trials){
boot3 <- sample(5000, 5000, replace = TRUE)
ind <- data3[ind,]
dat
<- regmedint(data = dat,
regmedint3 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)
out <- out$summary_myreg[1,1]
cde.est.boot <- out$summary_myreg[2,1]
pnde.est.boot <- out$summary_myreg[3,1]
tnie.est.boot <- out$summary_myreg[4,1]
tnde.est.boot <- out$summary_myreg[5,1]
pnie.est.boot <- out$summary_myreg[6,1]
te.est.boot <- out$summary_myreg[7,1]
pm.est.boot 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({
<- mclapply(trials, boot3, mc.cores = numCores)
results3
})
<- as.data.frame(do.call(rbind, results3))
results3.df apply(results3.df, 2, mean)
apply(results3.df, 2, sd)
<- regmedint(data = dat_logistic_M_logistic_Y,
regmedint4 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)
<- dat_logistic_M_logistic_Y
data4 <- function(trials){
boot4 <- sample(5000, 5000, replace = TRUE)
ind <- data4[ind,]
dat
<- regmedint(data = dat,
regmedint4 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)
out <- out$summary_myreg[1,1]
cde.est.boot <- out$summary_myreg[2,1]
pnde.est.boot <- out$summary_myreg[3,1]
tnie.est.boot <- out$summary_myreg[4,1]
tnde.est.boot <- out$summary_myreg[5,1]
pnie.est.boot <- out$summary_myreg[6,1]
te.est.boot <- out$summary_myreg[7,1]
pm.est.boot 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({
<- mclapply(trials, boot4, mc.cores = numCores)
results4
})
<- as.data.frame(do.call(rbind, results4))
results4.df 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.
Non-bootstrap
|
Bootstrap
|
|||
---|---|---|---|---|
Point Estimate | Standard Error | Point Estimate | Standard Error | |
CDE | 0.54832158 | 0.02201021 | 0.5476528 | 0.02170323 |
PNDE | 0.37202753 | 0.02273628 | 0.3719104 | 0.02241993 |
TNIE | 0.28120386 | 0.01480052 | 0.2807417 | 0.01490996 |
TNDE | 0.58513575 | 0.02334807 | 0.5843535 | 0.02312960 |
PNIE | 0.06809564 | 0.01196713 | 0.0682987 | 0.01233101 |
TE | 0.65323139 | 0.02075177 | 0.6526522 | 0.02123155 |
PM | 0.43048124 | 0.02380022 | 0.4304126 | 0.02344534 |
Non-bootstrap
|
Bootstrap
|
|||
---|---|---|---|---|
Point Estimate | Standard Error | Point Estimate | Standard Error | |
CDE | 0.2674768 | 0.02753958 | 0.27144584 | 0.02890304 |
PNDE | 0.5532311 | 0.02144037 | 0.55465628 | 0.02129513 |
TNIE | 0.0869584 | 0.01342216 | 0.08809824 | 0.01480665 |
TNDE | 0.6227790 | 0.02032353 | 0.62500131 | 0.01874403 |
PNIE | 0.0174105 | 0.00459631 | 0.01775321 | 0.00467464 |
TE | 0.6401895 | 0.02035298 | 0.64275452 | 0.01879988 |
PM | 0.1358323 | 0.02033647 | 0.13703290 | 0.02272391 |
Non-bootstrap
|
Bootstrap
|
|||
---|---|---|---|---|
Point Estimate | Standard Error | Point Estimate | Standard Error | |
CDE | 0.9731831 | 0.2696701 | 0.9779244 | 0.2712117 |
PNDE | 1.0001259 | 0.2838602 | 1.0044763 | 0.2857491 |
TNIE | 0.6608177 | 0.2557117 | 0.6606327 | 0.2528156 |
TNDE | 0.6089993 | 0.3558141 | 0.6078226 | 0.3542694 |
PNIE | 1.0519444 | 0.3044409 | 1.0572864 | 0.3059525 |
TE | 1.6609436 | 0.1614044 | 1.6651091 | 0.1599455 |
PM | 0.5969717 | 0.1616112 | 0.5776469 | 0.1685690 |
Non-bootstrap
|
Bootstrap
|
|||
---|---|---|---|---|
Point Estimate | Standard Error | Point Estimate | Standard Error | |
CDE | 0.20341749 | 0.11831758 | 0.20089667 | 0.12168979 |
PNDE | 0.76137841 | 0.08613598 | 0.76132524 | 0.08675318 |
TNIE | 0.06535229 | 0.01561403 | 0.06497380 | 0.01562340 |
TNDE | 0.82960632 | 0.08759429 | 0.82951442 | 0.08783860 |
PNIE | -0.00287562 | 0.01138042 | -0.00321537 | 0.01164672 |
TE | 0.82673070 | 0.08688896 | 0.82629904 | 0.08757983 |
PM | 0.11246227 | 0.02606395 | 0.11228721 | 0.02634418 |