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.06592388 -0.03560858 1.11775019
## pnde 0.488930417 0.21049248 2.3227928 0.02019028 0.07637274 0.90148809
## tnie 0.018240025 0.03706111 0.4921608 0.62260566 -0.05439841 0.09087846
## tnde 0.498503455 0.21209540 2.3503737 0.01875457 0.08280410 0.91420281
## pnie 0.008666987 0.02730994 0.3173565 0.75097309 -0.04485951 0.06219348
## te 0.507170442 0.21090051 2.4047853 0.01618197 0.09381303 0.92052785
## pm 0.045436278 0.09119614 0.4982259 0.61832484 -0.13330488 0.22417743boot packageThe 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 Scalemodelr packageIn 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 × 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.413tib_obj2 %>%
mutate(lower_delta = confint(regmedint_obj)[,"lower"],
upper_delta = confint(regmedint_obj)[,"upper"])## # A tibble: 7 × 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.133 0.224