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)## Error in library(furrr): there is no package called '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))## Error in `mutate()`:
## ℹ In argument: `boot_fit = future_map(...)`.
## Caused by error in `future_map()`:
## ! could not find function "future_map"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))## Error in eval(expr, envir, enclos): object 'tib_obj' not found## Error in eval(expr, envir, enclos): object 'tib_obj2' not found