inst/doc/vig_06_delta_boot.R

In regmedint: Regression-Based Causal Mediation Analysis with Interaction and Effect Modification Terms

## ----message = FALSE, tidy = FALSE, echo = F----------------------------------
## knitr configuration: http://yihui.name/knitr/options#chunk_options
library(knitr)
showMessage <- FALSE
showWarning <- TRUE
set_alias(w = "fig.width", h = "fig.height", res = "results")
opts_chunk$set(comment = "##", error= TRUE, warning = showWarning, message = showMessage,
               tidy = FALSE, cache = FALSE, echo = TRUE,
               fig.width = 7, fig.height = 7,
               fig.path = "man/figures")

## -----------------------------------------------------------------------------
library(regmedint)
library(tidyverse)

expit <- function(x){exp(x)/(1 + exp(x))}

## -----------------------------------------------------------------------------
library(parallel)
library(MASS)
numCores <- detectCores()
numCores

## -----------------------------------------------------------------------------
# Number of bootstrap
trials <- 1:5000
seed <-  3104

## -----------------------------------------------------------------------------
# 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))

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----eval = FALSE-------------------------------------------------------------
#  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)

## ----message = FALSE, tidy = FALSE, echo = F----------------------------------
m_lin_y_lin <- cbind.data.frame(c(0.54832158, 0.37202753, 0.28120386, 0.58513575,
                                  0.06809564, 0.65323139, 0.43048124),
                                c(0.02201021, 0.02273628, 0.01480052, 0.02334807,
                                  0.01196713, 0.02075177, 0.02380022),
                                c(0.5476528, 0.3719104, 0.2807417, 0.5843535,
                                  0.0682987, 0.6526522, 0.4304126),
                                c(0.02170323, 0.02241993, 0.01490996, 0.02312960,
                                  0.01233101, 0.02123155, 0.02344534))
row.names(m_lin_y_lin) <- c("CDE", "PNDE", "TNIE", "TNDE", "PNIE", "TE", "PM")
colnames(m_lin_y_lin) <- c(rep(c("Point Estimate", "Standard Error"), 2))


m_log_y_lin <- cbind.data.frame(c(0.2674768, 0.5532311, 0.0869584, 0.6227790,
                                  0.0174105, 0.6401895, 0.1358323),
                                c(0.02753958, 0.02144037, 0.01342216, 0.02032353,
                                  0.00459631, 0.02035298, 0.02033647),
                                c(0.27144584, 0.55465628, 0.08809824, 0.62500131,
                                  0.01775321, 0.64275452, 0.13703290),
                                c(0.02890304, 0.02129513, 0.01480665, 0.01874403,
                                  0.00467464, 0.01879988, 0.02272391))

m_lin_y_log <- cbind.data.frame(c(0.9731831, 1.0001259, 0.6608177, 0.6089993,
                                  1.0519444, 1.6609436, 0.5969717),
                                c(0.2696701, 0.2838602, 0.2557117, 0.3558141,
                                  0.3044409, 0.1614044, 0.1616112),
                                c(0.9779244, 1.0044763, 0.6606327, 0.6078226,
                                  1.0572864, 1.6651091, 0.5776469),
                                c(0.2712117, 0.2857491, 0.2528156, 0.3542694,
                                  0.3059525, 0.1599455, 0.1685690))

m_log_y_log <- cbind.data.frame(c(0.20341749, 0.76137841, 0.06535229, 0.82960632,
                                  -0.00287562, 0.82673070, 0.11246227),
                                c(0.11831758, 0.08613598, 0.01561403, 0.08759429,
                                  0.01138042, 0.08688896, 0.02606395),
                                c(0.20089667, 0.76132524, 0.06497380, 0.82951442,
                                  -0.00321537, 0.82629904, 0.11228721),
                                c(0.12168979, 0.08675318, 0.01562340, 0.08783860,
                                  0.01164672, 0.08757983, 0.02634418))

row_name <- c("CDE", "PNDE", "TNIE", "TNDE", "PNIE", "TE", "PM")
col_name <- c(rep(c("Point Estimate", "Standard Error"), 2))
row.names(m_lin_y_lin) <- row.names(m_log_y_lin) <-
  row.names(m_lin_y_log) <- row.names(m_log_y_log) <- row_name
colnames(m_lin_y_lin) <- colnames(m_log_y_lin) <-
  colnames(m_lin_y_log) <- colnames(m_log_y_log) <- col_name

## ----message = FALSE, tidy = FALSE, echo = F----------------------------------
# output table
library(kableExtra)
library(formattable)

kbl1 <- knitr::kable(m_lin_y_lin, align = "cccc", col.names = col_name,  digits = 8)
kbl2 <- knitr::kable(m_log_y_lin, align = "cccc", col.names = col_name,  digits = 8)
kbl3 <- knitr::kable(m_lin_y_log, align = "cccc", col.names = col_name,  digits = 8)
kbl4 <- knitr::kable(m_log_y_log, align = "cccc", col.names = col_name,  digits = 8)

# formatting:
# https://haozhu233.github.io/kableExtra/awesome_table_in_html.html

## ----message = FALSE, tidy = FALSE, echo = F----------------------------------
kbl1 %>%
  kable_classic(html_font = "Garamond") %>%
  add_header_above(c(" " = 1, "Non-bootstrap" = 2, "Bootstrap" = 2))

## ----message = FALSE, tidy = FALSE, echo = F----------------------------------
kbl2 %>%
  kable_classic(html_font = "Garamond") %>%
  add_header_above(c(" " = 1, "Non-bootstrap" = 2, "Bootstrap" = 2))

## ----message = FALSE, tidy = FALSE, echo = F----------------------------------
kbl3 %>%
  kable_classic(html_font = "Garamond") %>%
  add_header_above(c(" " = 1, "Non-bootstrap" = 2, "Bootstrap" = 2))

## ----message = FALSE, tidy = FALSE, echo = F----------------------------------
kbl4 %>%
  kable_classic(html_font = "Garamond") %>%
  add_header_above(c(" " = 1, "Non-bootstrap" = 2, "Bootstrap" = 2))