Using bootstrapping with regemedint

Kazuki Yoshida

2021-05-11

Bootstrapping support

The regmedint function itself does not contain a bootstrap standard error option, which may be perferred in some settings. However, it is relatively easy to implmement in R using the regmedint() function and the corresponding coef() point estimate extraction method.

set.seed(138087069)
library(regmedint)
library(tidyverse)
## Prepare dataset
data(vv2015)
## Main fit
regmedint_obj <- regmedint(data = vv2015,
                           ## Variables
                           yvar = "y",
                           avar = "x",
                           mvar = "m",
                           cvar = c("c"),
                           eventvar = "event",
                           ## Values at which effects are evaluated
                           a0 = 0,
                           a1 = 1,
                           m_cde = 1,
                           c_cond = 0.5,
                           ## Model types
                           mreg = "logistic",
                           yreg = "survAFT_weibull",
                           ## Additional specification
                           interaction = TRUE,
                           casecontrol = FALSE)
coef(summary(regmedint_obj))
##              est         se         Z           p       lower      upper
## cde  0.541070807 0.29422958 1.8389409 0.065923882 -0.03560858 1.11775019
## pnde 0.488930417 0.21049248 2.3227928 0.020190284  0.07637274 0.90148809
## tnie 0.018240025 0.03706111 0.4921608 0.622605663 -0.05439841 0.09087846
## tnde 0.498503455 0.21209540 2.3503737 0.018754573  0.08280410 0.91420281
## pnie 0.008666987 0.02730994 0.3173565 0.750973092 -0.04485951 0.06219348
## te   0.507170442 0.21090051 2.4047853 0.016181972  0.09381303 0.92052785
## pm   0.045436278 0.01725694 2.6329276 0.008465238  0.01161329 0.07925926

boot package

The boot package is the classical way to perform bootstrapping in R. It requires defining a wrapper function.

library(boot)
## Define a wrapper function
regmedint_boot <- function(data, ind)  {
    ## Note the change in the data argument.
    regmedint_obj <- regmedint(data = data[ind,],
                               ## Variables
                               yvar = "y",
                               avar = "x",
                               mvar = "m",
                               cvar = c("c"),
                               eventvar = "event",
                               ## Values at which effects are evaluated
                               a0 = 0,
                               a1 = 1,
                               m_cde = 1,
                               c_cond = 0.5,
                               ## Model types
                               mreg = "logistic",
                               yreg = "survAFT_weibull",
                               ## Additional specification
                               interaction = TRUE,
                               casecontrol = FALSE)
    coef(regmedint_obj)
}
## Run bootstrapping
ncpus <- 1
## For parallization, use the following instead.
## ncpus <- parallel::detectCores()
boot_obj <- boot(data = vv2015, statistic = regmedint_boot, R = 1000,
                 ## For palatalization
                 ## See https://cran.r-project.org/package=boot
                 parallel = "multicore",
                 ncpus = ncpus)
## Confidence interval for the pm
boot.ci(boot_obj, type = "basic", index = 7)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_obj, type = "basic", index = 7)
## 
## Intervals : 
## Level      Basic         
## 95%   (-0.4322,  0.4106 )  
## Calculations and Intervals on Original Scale

modelr package

In the tidyverse ecosystem, the modelr package can be used to provide a potentially more flexible workflow in some settings.

library(modelr)

library(future)
future::plan(sequential)
## For parallization, use the following instead.
## future::plan(multiprocess)
library(furrr)

## Bootstrapping
tib_obj <- vv2015 %>%
    modelr::bootstrap(n = 1000) %>%
    ## Resampled data objects are in the list column named strap.
    mutate(boot_fit = future_map(strap, function(strap) {
        ## Note the change in the data argument.
        regmedint_obj <- regmedint(data = as_tibble(strap),
                                   ## Variables
                                   yvar = "y",
                                   avar = "x",
                                   mvar = "m",
                                   cvar = c("c"),
                                   eventvar = "event",
                                   ## Values at which effects are evaluated
                                   a0 = 0,
                                   a1 = 1,
                                   m_cde = 1,
                                   c_cond = 0.5,
                                   ## Model types
                                   mreg = "logistic",
                                   yreg = "survAFT_weibull",
                                   ## Additional specification
                                   interaction = TRUE,
                                   casecontrol = FALSE)
        ## Trick to return a row tibble
        mat <- t(matrix(coef(regmedint_obj)))
        colnames(mat) <- names(coef(regmedint_obj))
        as_tibble(mat)
    })) %>%
    select(-strap) %>%
    unnest(cols = c(boot_fit))

tib_obj2 <- tib_obj %>%
    pivot_longer(-.id) %>%
    mutate(name = factor(name, levels = names(coef(regmedint_obj)))) %>%
    group_by(name) %>%
    summarize(lower_boot = quantile(value, probs = 0.025),
              upper_boot = quantile(value, probs = 0.975))
tib_obj2
## # A tibble: 7 x 3
##   name  lower_boot upper_boot
##   <fct>      <dbl>      <dbl>
## 1 cde      -0.0968     1.17  
## 2 pnde      0.0339     0.962 
## 3 tnie     -0.0800     0.144 
## 4 tnde      0.0425     0.971 
## 5 pnie     -0.0655     0.0936
## 6 te        0.0519     0.984 
## 7 pm       -0.298      0.413

Comparison to the delta method confidence intervals

tib_obj2 %>%
    mutate(lower_delta = confint(regmedint_obj)[,"lower"],
           upper_delta = confint(regmedint_obj)[,"upper"])
## # A tibble: 7 x 5
##   name  lower_boot upper_boot lower_delta upper_delta
##   <fct>      <dbl>      <dbl>       <dbl>       <dbl>
## 1 cde      -0.0968     1.17       -0.0356      1.12  
## 2 pnde      0.0339     0.962       0.0764      0.901 
## 3 tnie     -0.0800     0.144      -0.0544      0.0909
## 4 tnde      0.0425     0.971       0.0828      0.914 
## 5 pnie     -0.0655     0.0936     -0.0449      0.0622
## 6 te        0.0519     0.984       0.0938      0.921 
## 7 pm       -0.298      0.413       0.0116      0.0793