R/bootstrap.R
In CausalInference/gfoRmula: Parametric G-Formula
Defines functions
bootstrap_helper_with_trycatch
bootstrap_helper
Documented in
bootstrap_helper
#' Bootstrap Observed Data and Simulate Under All Interventions
#'
#' This internal function bootstraps the observed data (i.e., resamples the observed data set with replacement to
#' construct bootstrap confidence intervals and standard errors). Then, the function simulates data
#' using the resampled dataset to estimate the survival outcome, binary end-of-follow-up outcome, or
#' continuous end-of-follow-up outcome.
#'
#' @param r Integer specifying the index of the current iteration of the bootstrap.
#' @param time_points Number of time points to simulate.
#' @param obs_data Data table containing the observed data.
#' @param bootseeds Vector of integers specifying the seeds. One seed is used to initialize each bootstrap iteration.
#' @param outcome_type Character string specifying the "type" of the outcome. The possible "types" are: \code{"survival"}, \code{"continuous_eof"}, and \code{"binary_eof"}.
#' @param intvars List, whose elements are vectors of character strings. The kth vector in \code{intvars} specifies the name(s) of the variable(s) to be intervened
#' on in each round of the simulation under the kth intervention in \code{interventions}.
#' @param interventions List, whose elements are lists of vectors. Each list in \code{interventions} specifies a unique intervention on the relevant variable(s) in \code{intvars}. Each vector contains a function
#' implementing a particular intervention on a single variable, optionally
#' followed by one or more "intervention values" (i.e.,
#' integers used to specify the treatment regime).
#' @param int_times List, whose elements are lists of vectors. The kth list in \code{int_times} corresponds to the kth intervention in \code{interventions}. Each vector specifies the time points in which the relevant intervention is applied on the corresponding variable in \code{intvars}.
#' When an intervention is not applied, the simulated natural course value is used. By default, this argument is set so that all interventions are applied in all time points.
#' @param ref_int Integer denoting the intervention to be used as the
#' reference for calculating the risk ratio and risk difference. 0 denotes the
#' natural course, while subsequent integers denote user-specified
#' interventions in the order that they are
#' named in \code{interventions}.
#' @param covparams List of vectors, where each vector contains information for
#' one parameter used in the modeling of the time-varying covariates (e.g.,
#' model statement, family, link function, etc.). Each vector
#' must be the same length as \code{covnames} and in the same order.
#' If a parameter is not required for a certain covariate, it
#' should be set to \code{NA} at that index.
#' @param covnames Vector of character strings specifying the names of the time-varying covariates in \code{obs_data}.
#' @param covtypes Vector of character strings specifying the "type" of each time-varying covariate included in \code{covnames}. The possible "types" are: \code{"binary"}, \code{"normal"}, \code{"categorical"}, \code{"bounded normal"}, \code{"zero-inflated normal"}, \code{"truncated normal"}, \code{"absorbing"}, \code{"categorical time"}, and \code{"custom"}.
#' @param covfits_custom Vector containing custom fit functions for time-varying covariates that
#' do not fall within the pre-defined covariate types. It should be in
#' the same order \code{covnames}. If a custom fit function is not
#' required for a particular covariate (e.g., if the first
#' covariate is of type \code{"binary"} but the second is of type \code{"custom"}), then that
#' index should be set to \code{NA}.
#' @param covpredict_custom Vector containing custom prediction functions for time-varying
#' covariates that do not fall within the pre-defined covariate types.
#' It should be in the same order as \code{covnames}. If a custom
#' prediction function is not required for a particular
#' covariate, then that index should be set to \code{NA}.
#' @param basecovs Vector of character strings specifying the names of baseline covariates in \code{obs_data}.
#' @param histvars List of vectors. The kth vector specifies the names of the variables for which the kth history function
#' in \code{histories} is to be applied.
#' @param histvals List of length two. The first element is a numeric vector specifying the lags used in the model statements (e.g., if \code{lag1_varname} and \code{lag2_varname} were included in the model statements, this vector would be \code{c(1,2)}). The second element is a numeric vector specifying the lag averages used in the model statements.
#' @param histories Vector of history functions to apply to the variables specified in \code{histvars}.
#' @param ymodel Model statement for the outcome variable.
#' @param ymodel_fit_custom Function specifying a custom outcome model. See the vignette "Using Custom Outcome Models in gfoRmula" for details.
#' @param ymodel_predict_custom Function obtaining predictions from the custom outcome model specified in \code{ymodel_fit_custom}. See the vignette "Using Custom Outcome Models in gfoRmula" for details.
#' @param yrestrictions List of vectors. Each vector containins as its first entry
#' a condition and its second entry an integer. When the
#' condition is \code{TRUE}, the outcome variable is simulated
#' according to the fitted model; when the condition is \code{FALSE},
#' the outcome variable takes on the value in the second entry.
#' @param compevent_restrictions List of vectors. Each vector containins as its first entry
#' a condition and its second entry an integer. When the
#' condition is \code{TRUE}, the competing event variable is simulated
#' according to the fitted model; when the condition is \code{FALSE},
#' the competing event variable takes on the value in the
#' second entry.
#' @param restrictions List of vectors. Each vector contains as its first entry a covariate for which
#' \emph{a priori} knowledge of its distribution is available; its second entry a condition
#' under which no knowledge of its distribution is available and that must be \code{TRUE}
#' for the distribution of that covariate given that condition to be estimated via a parametric
#' model or other fitting procedure; its third entry a function for estimating the distribution
#' of that covariate given the condition in the second entry is false such that \emph{a priori} knowledge
#' of the covariate distribution is available; and its fourth entry a value used by the function in the
#' third entry. The default is \code{NA}.
#' @param comprisk Logical scalar indicating the presence of a competing event.
#' @param compevent_model Model statement for the competing event variable.
#' @param time_name Character string specifying the name of the time variable in \code{obs_data}.
#' @param outcome_name Character string specifying the name of the outcome variable in \code{obs_data}.
#' @param compevent_name Character string specifying the name of the competing event variable in \code{obs_data}.
#' @param ranges List of vectors. Each vector contains the minimum and
#' maximum values of one of the covariates in \code{covnames}.
#' @param parallel Logical scalar indicating whether to parallelize simulations of
#' different interventions to multiple cores.
#' @param ncores Integer specifying the number of cores to use in parallel
#' simulation.
#' @param max_visits A vector of one or more values denoting the maximum number of times
#' a binary covariate representing a visit process may be missed before
#' the individual is censored from the data (in the observed data) or
#' a visit is forced (in the simulated data). Multiple values exist in the
#' vector when the modeling of more than covariate is attached to a visit
#' process. A value of \code{NA} should be provided when there is no visit process.
#' @param hazardratio Logical scalar indicating whether the hazard ratio should be computed between two interventions.
#' @param intcomp List of two numbers indicating a pair of interventions to be compared by a hazard ratio.
#' The default is \code{NA}, resulting in no hazard ratio calculation.
#' @param boot_diag Logical scalar indicating whether to return the coefficients, standard errors, and variance-covariance matrices of the parameters of the fitted models in the bootstrap samples. The default is \code{FALSE}.
#' @param nsimul Number of subjects for whom to simulate data. By default, this argument is set
#' equal to the number of subjects in \code{obs_data}.
#' @param baselags Logical scalar for specifying the convention used for lagi and lag_cumavgi terms in the model statements when pre-baseline times are not
#' included in \code{obs_data} and when the current time index, \eqn{t}, is such that \eqn{t < i}. If this argument is set to \code{FALSE}, the value
#' of all lagi and lag_cumavgi terms in this context are set to 0 (for non-categorical covariates) or the reference
#' level (for categorical covariates). If this argument is set to \code{TRUE}, the value of lagi and lag_cumavgi terms
#' are set to their values at time 0. The default is \code{FALSE}.
#' @param below_zero_indicator Logical scalar indicating whether the observed data set contains rows for time \eqn{t < 0}.
#' @param min_time Numeric scalar specifying lowest value of time \eqn{t} in the observed data set.
#' @param show_progress Logical scalar indicating whether to print a progress bar for the number of bootstrap samples completed in the R console. This argument is only applicable when \code{parallel} is set to \code{FALSE} and bootstrap samples are used (i.e., \code{nsamples} is set to a value greater than 0). The default is \code{TRUE}.
#' @param pb Progress bar R6 object. See \code{\link[progress]{progress_bar}} for further details.
#' @param int_visit_type Vector of logicals. The kth element is a logical specifying whether to carry forward the intervened value (rather than the natural value) of the treatment variables(s) when performing a carry forward restriction type for the kth intervention in \code{interventions}.
#' When the kth element is set to \code{FALSE}, the natural value of the treatment variable(s) in the kth intervention in \code{interventions} will be carried forward.
#' By default, this argument is set so that the intervened value of the treatment variable(s) is carried forward for all interventions.
#' @param sim_trunc Logical scalar indicating whether to truncate simulated covariates to their range in the observed data set. This argument is only applicable for covariates of type \code{"normal"}, \code{"bounded normal"}, \code{"truncated normal"}, and \code{"zero-inflated normal"}.
#' @param ... Other arguments
#' @return A list with the following components:
#' \item{Result}{Matrix containing risks over time under the natural course and under each user-specific intervention.}
#' \item{ResultRatio}{Matrix containing risk ratios over time under the natural course and under each user-specific intervention.}
#' \item{ResultDiff}{Matrix containing risk differences over time under the natural course and under each user-specific intervention.}
#' \item{bootcoeffs}{List of the coefficients of the fitted models. If the argument \code{boot_diag} is set to \code{FALSE}, a value of \code{NA} is given.}
#' \item{bootstderrs}{List of the standard errors of the coefficients of the fitted models. If the argument \code{boot_diag} is set to \code{FALSE}, a value of \code{NA} is given.}
#'
#' @keywords internal
#' @import data.table
bootstrap_helper <- function(r, time_points, obs_data, bootseeds, outcome_type,
intvars, interventions, int_times, ref_int,
covparams, covnames, covtypes, covfits_custom, covpredict_custom, basecovs, histvars, histvals, histories,
ymodel, ymodel_fit_custom, ymodel_predict_custom,
yrestrictions, compevent_restrictions, restrictions,
comprisk, compevent_model,
time_name, outcome_name, compevent_name,
ranges, parallel, ncores, max_visits,
hazardratio, intcomp, boot_diag, nsimul, baselags,
below_zero_indicator, min_time, show_progress, pb,
int_visit_type, sim_trunc, ...){
set.seed(bootseeds[r])
data_len <- length(unique(obs_data$newid))
ids <- as.data.table(sample(1:data_len, data_len, replace = TRUE))
ids[, 'bid' := 1:data_len]
colnames(ids) <- c("newid", "bid")
resample_data <- copy(obs_data)
setkey(resample_data, "newid")
resample_data <- resample_data[J(ids), allow.cartesian = TRUE] # create the new data set names "sample"
resample_data[, 'newid' := resample_data$bid]
resample_data[, 'bid' := NULL]
resample_data_geq_0 <- resample_data[resample_data[[time_name]] >= 0]
# Fit models for covariates, outcome, and competing event (if any)
fitcov <- pred_fun_cov(covparams = covparams, covnames = covnames, covtypes = covtypes,
covfits_custom = covfits_custom, restrictions = restrictions,
time_name = time_name, obs_data = resample_data_geq_0,
model_fits = FALSE)
fitY <- pred_fun_Y(ymodel, yrestrictions, outcome_type, outcome_name, time_name, resample_data_geq_0,
model_fits = FALSE, ymodel_fit_custom = ymodel_fit_custom)
if (comprisk){
fitD <- pred_fun_D(compevent_model, compevent_restrictions, resample_data_geq_0, model_fits = FALSE)
} else {
fitD <- NA
}
len <- length(unique(resample_data$newid))
# If the number of user desired simulations differs from the number of individuals in
# the observed dataset, sample the desired number of observed IDs with replacement
if (nsimul < len){
ids <- as.data.table(sort(sample(unique(resample_data$newid), nsimul, replace = TRUE)))
colnames(ids) <- "newid"
ids[, 'bid' := seq_len(.N)]
resample_data <- merge(ids, resample_data, all.x = TRUE, by = "newid")
resample_data[, 'newid' := resample_data$bid]
resample_data[, 'bid' := NULL]
} else if (nsimul > len){
ids <- as.data.table(sample(unique(resample_data$newid), nsimul, replace = TRUE))
ids[, 'newid' := 1:nsimul]
colnames(ids) <- c("newid", "bid")
setkeyv(resample_data, "newid")
resample_data <- resample_data[J(ids), allow.cartesian = TRUE]
resample_data[, 'newid' := resample_data$bid]
resample_data[, 'bid' := NULL]
}
comb_interventions <- interventions
comb_intvars <- intvars
comb_int_times <- int_times
# Simulate data under different interventions
pools <- lapply(seq_along(comb_interventions), FUN = function(i){
simulate(fitcov = fitcov, fitY = fitY, fitD = fitD,
ymodel_predict_custom = ymodel_predict_custom,
yrestrictions = yrestrictions,
compevent_restrictions = compevent_restrictions,
restrictions = restrictions,
outcome_name = outcome_name, compevent_name = compevent_name,
time_name = time_name,
intvars = comb_intvars[[i]], interventions = comb_interventions[[i]], int_times = comb_int_times[[i]],
histvars = histvars, histvals = histvals, histories = histories,
covparams = covparams, covnames = covnames, covtypes = covtypes,
covpredict_custom = covpredict_custom,
basecovs = basecovs,
comprisk = comprisk, ranges = ranges,
outcome_type = outcome_type,
subseed = bootseeds[r], time_points = time_points,
obs_data = resample_data, parallel = FALSE, max_visits = max_visits,
baselags = baselags, below_zero_indicator = below_zero_indicator,
min_time = min_time, show_progress = show_progress, pb = pb,
int_visit_type = int_visit_type[i], sim_trunc = sim_trunc, ...)
})
nat_pool <- pools[[1]] # Simulated data under natural course
pools <- pools[-1] # Simulated data under various interventions
if (outcome_type == 'survival'){
param <- 'poprisk'
} else if (outcome_type == 'continuous_eof'){
param <- 'Ey'
} else if (outcome_type == 'binary_eof'){
param <- 'Py'
}
# Initialize result matrix
if (grepl('eof', outcome_type)){
result_ratio <- result_diff <- int_result <- rep(NA, length(pools) + 1)
} else {
result_ratio <- result_diff <- int_result <-
matrix(NA, nrow = length(pools) + 1, ncol = time_points)
}
# Calculate mean risk at each time point under natural course
if (grepl('eof', outcome_type)){
nat_result <- mean(nat_pool[[param]], na.rm = TRUE)
} else {
nat_result <- tapply(nat_pool[[param]], nat_pool[[time_name]], FUN = mean)
}
if (ref_int == 0){
ref_result <- nat_result
} else {
if (grepl('eof', outcome_type)){
ref_result <- mean(pools[[ref_int]][[param]], na.rm = TRUE)
} else {
# Calculate mean risk at each time point for specified reference intervention
ref_result <- tapply(pools[[ref_int]][[param]], pools[[ref_int]][[time_name]], FUN = mean)
}
}
# Calculate risk ratio under natural course
# Calculate risk ratio for remaining interventions
if (grepl('eof', outcome_type)){
int_result[1] <- nat_result
int_result[-1] <- sapply(pools, FUN = function(pool){
mean(pool[[param]], na.rm = TRUE)
})
result_ratio <- int_result / ref_result
result_diff <- int_result - ref_result
} else {
int_result[1, ] <- nat_result
result_ratio[1, ] <- int_result[1, ] / ref_result
result_diff[1, ] <- int_result[1, ] - ref_result
if (length(comb_interventions) > 1){
for (i in 2:(length(pools) + 1)){
int_result[i, ] <- tapply(pools[[i - 1]][[param]], pools[[i - 1]][[time_name]], FUN = mean)
result_ratio[i, ] <- int_result[i, ] / ref_result
result_diff[i, ] <- int_result[i, ] - ref_result
}
}
}
# Calculate hazard ratio for specified interventions
if (hazardratio){
# Generate dataset containing failure/censor time information for each subject
# under each intervention
pools_hr <- lapply(seq_along(intcomp), FUN = hr_helper, intcomp = intcomp,
time_name = time_name, pools = pools)
data_hr <- rbindlist(pools_hr)
names(data_hr)[names(data_hr) == time_name] <- "t0"
names(data_hr)[names(data_hr) == outcome_name] <- "Y"
# Factor event variable
data_hr$event <- factor(data_hr$Ycomp, 0:2, labels=c("censor", "Y", "D"))
if (comprisk){
# Calculate subdistribution hazard ratio
hr_data <- survival::finegray(survival::Surv(t0, event) ~ ., data = data_hr, etype = "Y")
hr_res <- survival::coxph(survival::Surv(fgstart, fgstop, fgstatus) ~ regime, data = hr_data)
hr_res <- exp(hr_res$coefficients)
}
else {
# Calculate cause-specific hazard ratio
hr_res <- survival::coxph(formula = survival::Surv(t0, Y == "1") ~ regime, data = data_hr)
hr_res <- exp(hr_res$coefficients)
}
if (is.na(hr_res)){
warning(paste0('NA value for the hazard ratio found in bootstrap replicate ', r,
'. The analysis should likely be repeated with more parsimonious models.'))
}
} else {
hr_res <- NA
}
if (time_points > 1){
fits <- fitcov
fits[[length(fits) + 1]] <- fitY
} else {
fits <- list(fitY)
}
if (!is.na(fitD)[[1]]){
fits[[length(fits) + 1]] <- fitD
}
if (boot_diag){
bootcoeffs <- get_coeffs(fits = fits, fitD = fitD, time_points = time_points,
outcome_name = outcome_name, compevent_name = compevent_name,
covnames = covnames)
bootstderrs <- get_stderrs(fits = fits, fitD = fitD, time_points = time_points,
outcome_name = outcome_name, compevent_name = compevent_name,
covnames = covnames)
bootvcovs <- get_vcovs(fits = fits, fitD = fitD, time_points = time_points,
outcome_name = outcome_name, compevent_name = compevent_name,
covnames = covnames)
} else {
bootcoeffs <- NA
bootstderrs <- NA
bootvcovs <- NA
}
# Check NA values
if (any(is.na(int_result))){
warning(paste0('NA values for the g-form risk found in bootstrap replicate ', r,
'. The analysis should likely be repeated with more parsimonious models.'))
}
final <- list(Result = int_result, ResultRatio = result_ratio, ResultDiff = result_diff, ResultHR = hr_res,
bootcoeffs = bootcoeffs, bootstderrs = bootstderrs, bootvcovs = bootvcovs)
return (final)
}
bootstrap_helper_with_trycatch <- function(r, time_points, obs_data, bootseeds, outcome_type,
intvars, interventions, int_times, ref_int,
covparams, covnames, covtypes, covfits_custom, covpredict_custom, basecovs, histvars, histvals, histories,
ymodel, ymodel_fit_custom, ymodel_predict_custom,
yrestrictions, compevent_restrictions, restrictions,
comprisk, compevent_model,
time_name, outcome_name, compevent_name,
ranges, parallel, ncores, max_visits,
hazardratio, intcomp, boot_diag, nsimul, baselags,
below_zero_indicator, min_time, show_progress, pb,
int_visit_type, sim_trunc, ...){
tryCatch({
bootstrap_helper(
r = r, time_points = time_points, obs_data = obs_data, bootseeds = bootseeds, outcome_type = outcome_type,
intvars = intvars, interventions = interventions, int_times = int_times, ref_int = ref_int,
covparams = covparams, covnames = covnames, covtypes = covtypes, covfits_custom = covfits_custom, covpredict_custom = covpredict_custom, basecovs = basecovs, histvars = histvars, histvals = histvals, histories = histories,
ymodel = ymodel, ymodel_fit_custom = ymodel_fit_custom, ymodel_predict_custom = ymodel_predict_custom,
yrestrictions = yrestrictions, compevent_restrictions = compevent_restrictions, restrictions = restrictions,
comprisk = comprisk, compevent_model = compevent_model,
time_name = time_name, outcome_name = outcome_name, compevent_name = compevent_name,
ranges = ranges, parallel = parallel, ncores = ncores, max_visits = max_visits,
hazardratio = hazardratio, intcomp = intcomp, boot_diag = boot_diag, nsimul = nsimul, baselags = baselags,
below_zero_indicator = below_zero_indicator, min_time = min_time, show_progress = show_progress, pb = pb,
int_visit_type = int_visit_type, sim_trunc = sim_trunc, ...)
},
error = function(e){
warning(paste0('An error occured in bootstrap replicate ', r,
'. The analysis should likely be repeated with more parsimonious models. The error message encountered is:\n', e))
if (grepl('eof', outcome_type)){
result_ratio <- result_diff <- int_result <- rep(NA, length(interventions))
} else {
result_ratio <- result_diff <- int_result <-
matrix(NA, nrow = length(interventions), ncol = time_points)
}
final <- list(Result = int_result, ResultRatio = result_ratio, ResultDiff = result_diff, ResultHR = NA,
bootcoeffs = NULL, bootstderrs = NULL, bootvcovs = NULL)
final
})
}
CausalInference/gfoRmula documentation built on Oct. 1, 2024, 8:36 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions bootstrap_helper_with_trycatch bootstrap_helper
Documented in bootstrap_helper
#' Bootstrap Observed Data and Simulate Under All Interventions
#'
#' This internal function bootstraps the observed data (i.e., resamples the observed data set with replacement to
#' construct bootstrap confidence intervals and standard errors). Then, the function simulates data
#' using the resampled dataset to estimate the survival outcome, binary end-of-follow-up outcome, or
#' continuous end-of-follow-up outcome.
#'
#' @param r Integer specifying the index of the current iteration of the bootstrap.
#' @param time_points Number of time points to simulate.
#' @param obs_data Data table containing the observed data.
#' @param bootseeds Vector of integers specifying the seeds. One seed is used to initialize each bootstrap iteration.
#' @param outcome_type Character string specifying the "type" of the outcome. The possible "types" are: \code{"survival"}, \code{"continuous_eof"}, and \code{"binary_eof"}.
#' @param intvars List, whose elements are vectors of character strings. The kth vector in \code{intvars} specifies the name(s) of the variable(s) to be intervened
#' on in each round of the simulation under the kth intervention in \code{interventions}.
#' @param interventions List, whose elements are lists of vectors. Each list in \code{interventions} specifies a unique intervention on the relevant variable(s) in \code{intvars}. Each vector contains a function
#' implementing a particular intervention on a single variable, optionally
#' followed by one or more "intervention values" (i.e.,
#' integers used to specify the treatment regime).
#' @param int_times List, whose elements are lists of vectors. The kth list in \code{int_times} corresponds to the kth intervention in \code{interventions}. Each vector specifies the time points in which the relevant intervention is applied on the corresponding variable in \code{intvars}.
#' When an intervention is not applied, the simulated natural course value is used. By default, this argument is set so that all interventions are applied in all time points.
#' @param ref_int Integer denoting the intervention to be used as the
#' reference for calculating the risk ratio and risk difference. 0 denotes the
#' natural course, while subsequent integers denote user-specified
#' interventions in the order that they are
#' named in \code{interventions}.
#' @param covparams List of vectors, where each vector contains information for
#' one parameter used in the modeling of the time-varying covariates (e.g.,
#' model statement, family, link function, etc.). Each vector
#' must be the same length as \code{covnames} and in the same order.
#' If a parameter is not required for a certain covariate, it
#' should be set to \code{NA} at that index.
#' @param covnames Vector of character strings specifying the names of the time-varying covariates in \code{obs_data}.
#' @param covtypes Vector of character strings specifying the "type" of each time-varying covariate included in \code{covnames}. The possible "types" are: \code{"binary"}, \code{"normal"}, \code{"categorical"}, \code{"bounded normal"}, \code{"zero-inflated normal"}, \code{"truncated normal"}, \code{"absorbing"}, \code{"categorical time"}, and \code{"custom"}.
#' @param covfits_custom Vector containing custom fit functions for time-varying covariates that
#' do not fall within the pre-defined covariate types. It should be in
#' the same order \code{covnames}. If a custom fit function is not
#' required for a particular covariate (e.g., if the first
#' covariate is of type \code{"binary"} but the second is of type \code{"custom"}), then that
#' index should be set to \code{NA}.
#' @param covpredict_custom Vector containing custom prediction functions for time-varying
#' covariates that do not fall within the pre-defined covariate types.
#' It should be in the same order as \code{covnames}. If a custom
#' prediction function is not required for a particular
#' covariate, then that index should be set to \code{NA}.
#' @param basecovs Vector of character strings specifying the names of baseline covariates in \code{obs_data}.
#' @param histvars List of vectors. The kth vector specifies the names of the variables for which the kth history function
#' in \code{histories} is to be applied.
#' @param histvals List of length two. The first element is a numeric vector specifying the lags used in the model statements (e.g., if \code{lag1_varname} and \code{lag2_varname} were included in the model statements, this vector would be \code{c(1,2)}). The second element is a numeric vector specifying the lag averages used in the model statements.
#' @param histories Vector of history functions to apply to the variables specified in \code{histvars}.
#' @param ymodel Model statement for the outcome variable.
#' @param ymodel_fit_custom Function specifying a custom outcome model. See the vignette "Using Custom Outcome Models in gfoRmula" for details.
#' @param ymodel_predict_custom Function obtaining predictions from the custom outcome model specified in \code{ymodel_fit_custom}. See the vignette "Using Custom Outcome Models in gfoRmula" for details.
#' @param yrestrictions List of vectors. Each vector containins as its first entry
#' a condition and its second entry an integer. When the
#' condition is \code{TRUE}, the outcome variable is simulated
#' according to the fitted model; when the condition is \code{FALSE},
#' the outcome variable takes on the value in the second entry.
#' @param compevent_restrictions List of vectors. Each vector containins as its first entry
#' a condition and its second entry an integer. When the
#' condition is \code{TRUE}, the competing event variable is simulated
#' according to the fitted model; when the condition is \code{FALSE},
#' the competing event variable takes on the value in the
#' second entry.
#' @param restrictions List of vectors. Each vector contains as its first entry a covariate for which
#' \emph{a priori} knowledge of its distribution is available; its second entry a condition
#' under which no knowledge of its distribution is available and that must be \code{TRUE}
#' for the distribution of that covariate given that condition to be estimated via a parametric
#' model or other fitting procedure; its third entry a function for estimating the distribution
#' of that covariate given the condition in the second entry is false such that \emph{a priori} knowledge
#' of the covariate distribution is available; and its fourth entry a value used by the function in the
#' third entry. The default is \code{NA}.
#' @param comprisk Logical scalar indicating the presence of a competing event.
#' @param compevent_model Model statement for the competing event variable.
#' @param time_name Character string specifying the name of the time variable in \code{obs_data}.
#' @param outcome_name Character string specifying the name of the outcome variable in \code{obs_data}.
#' @param compevent_name Character string specifying the name of the competing event variable in \code{obs_data}.
#' @param ranges List of vectors. Each vector contains the minimum and
#' maximum values of one of the covariates in \code{covnames}.
#' @param parallel Logical scalar indicating whether to parallelize simulations of
#' different interventions to multiple cores.
#' @param ncores Integer specifying the number of cores to use in parallel
#' simulation.
#' @param max_visits A vector of one or more values denoting the maximum number of times
#' a binary covariate representing a visit process may be missed before
#' the individual is censored from the data (in the observed data) or
#' a visit is forced (in the simulated data). Multiple values exist in the
#' vector when the modeling of more than covariate is attached to a visit
#' process. A value of \code{NA} should be provided when there is no visit process.
#' @param hazardratio Logical scalar indicating whether the hazard ratio should be computed between two interventions.
#' @param intcomp List of two numbers indicating a pair of interventions to be compared by a hazard ratio.
#' The default is \code{NA}, resulting in no hazard ratio calculation.
#' @param boot_diag Logical scalar indicating whether to return the coefficients, standard errors, and variance-covariance matrices of the parameters of the fitted models in the bootstrap samples. The default is \code{FALSE}.
#' @param nsimul Number of subjects for whom to simulate data. By default, this argument is set
#' equal to the number of subjects in \code{obs_data}.
#' @param baselags Logical scalar for specifying the convention used for lagi and lag_cumavgi terms in the model statements when pre-baseline times are not
#' included in \code{obs_data} and when the current time index, \eqn{t}, is such that \eqn{t < i}. If this argument is set to \code{FALSE}, the value
#' of all lagi and lag_cumavgi terms in this context are set to 0 (for non-categorical covariates) or the reference
#' level (for categorical covariates). If this argument is set to \code{TRUE}, the value of lagi and lag_cumavgi terms
#' are set to their values at time 0. The default is \code{FALSE}.
#' @param below_zero_indicator Logical scalar indicating whether the observed data set contains rows for time \eqn{t < 0}.
#' @param min_time Numeric scalar specifying lowest value of time \eqn{t} in the observed data set.
#' @param show_progress Logical scalar indicating whether to print a progress bar for the number of bootstrap samples completed in the R console. This argument is only applicable when \code{parallel} is set to \code{FALSE} and bootstrap samples are used (i.e., \code{nsamples} is set to a value greater than 0). The default is \code{TRUE}.
#' @param pb Progress bar R6 object. See \code{\link[progress]{progress_bar}} for further details.
#' @param int_visit_type Vector of logicals. The kth element is a logical specifying whether to carry forward the intervened value (rather than the natural value) of the treatment variables(s) when performing a carry forward restriction type for the kth intervention in \code{interventions}.
#' When the kth element is set to \code{FALSE}, the natural value of the treatment variable(s) in the kth intervention in \code{interventions} will be carried forward.
#' By default, this argument is set so that the intervened value of the treatment variable(s) is carried forward for all interventions.
#' @param sim_trunc Logical scalar indicating whether to truncate simulated covariates to their range in the observed data set. This argument is only applicable for covariates of type \code{"normal"}, \code{"bounded normal"}, \code{"truncated normal"}, and \code{"zero-inflated normal"}.
#' @param ... Other arguments
#' @return A list with the following components:
#' \item{Result}{Matrix containing risks over time under the natural course and under each user-specific intervention.}
#' \item{ResultRatio}{Matrix containing risk ratios over time under the natural course and under each user-specific intervention.}
#' \item{ResultDiff}{Matrix containing risk differences over time under the natural course and under each user-specific intervention.}
#' \item{bootcoeffs}{List of the coefficients of the fitted models. If the argument \code{boot_diag} is set to \code{FALSE}, a value of \code{NA} is given.}
#' \item{bootstderrs}{List of the standard errors of the coefficients of the fitted models. If the argument \code{boot_diag} is set to \code{FALSE}, a value of \code{NA} is given.}
#'
#' @keywords internal
#' @import data.table
bootstrap_helper <- function(r, time_points, obs_data, bootseeds, outcome_type,
intvars, interventions, int_times, ref_int,
covparams, covnames, covtypes, covfits_custom, covpredict_custom, basecovs, histvars, histvals, histories,
ymodel, ymodel_fit_custom, ymodel_predict_custom,
yrestrictions, compevent_restrictions, restrictions,
comprisk, compevent_model,
time_name, outcome_name, compevent_name,
ranges, parallel, ncores, max_visits,
hazardratio, intcomp, boot_diag, nsimul, baselags,
below_zero_indicator, min_time, show_progress, pb,
int_visit_type, sim_trunc, ...){
set.seed(bootseeds[r])
data_len <- length(unique(obs_data$newid))
ids <- as.data.table(sample(1:data_len, data_len, replace = TRUE))
ids[, 'bid' := 1:data_len]
colnames(ids) <- c("newid", "bid")
resample_data <- copy(obs_data)
setkey(resample_data, "newid")
resample_data <- resample_data[J(ids), allow.cartesian = TRUE] # create the new data set names "sample"
resample_data[, 'newid' := resample_data$bid]
resample_data[, 'bid' := NULL]
resample_data_geq_0 <- resample_data[resample_data[[time_name]] >= 0]
# Fit models for covariates, outcome, and competing event (if any)
fitcov <- pred_fun_cov(covparams = covparams, covnames = covnames, covtypes = covtypes,
covfits_custom = covfits_custom, restrictions = restrictions,
time_name = time_name, obs_data = resample_data_geq_0,
model_fits = FALSE)
fitY <- pred_fun_Y(ymodel, yrestrictions, outcome_type, outcome_name, time_name, resample_data_geq_0,
model_fits = FALSE, ymodel_fit_custom = ymodel_fit_custom)
if (comprisk){
fitD <- pred_fun_D(compevent_model, compevent_restrictions, resample_data_geq_0, model_fits = FALSE)
} else {
fitD <- NA
}
len <- length(unique(resample_data$newid))
# If the number of user desired simulations differs from the number of individuals in
# the observed dataset, sample the desired number of observed IDs with replacement
if (nsimul < len){
ids <- as.data.table(sort(sample(unique(resample_data$newid), nsimul, replace = TRUE)))
colnames(ids) <- "newid"
ids[, 'bid' := seq_len(.N)]
resample_data <- merge(ids, resample_data, all.x = TRUE, by = "newid")
resample_data[, 'newid' := resample_data$bid]
resample_data[, 'bid' := NULL]
} else if (nsimul > len){
ids <- as.data.table(sample(unique(resample_data$newid), nsimul, replace = TRUE))
ids[, 'newid' := 1:nsimul]
colnames(ids) <- c("newid", "bid")
setkeyv(resample_data, "newid")
resample_data <- resample_data[J(ids), allow.cartesian = TRUE]
resample_data[, 'newid' := resample_data$bid]
resample_data[, 'bid' := NULL]
}
comb_interventions <- interventions
comb_intvars <- intvars
comb_int_times <- int_times
# Simulate data under different interventions
pools <- lapply(seq_along(comb_interventions), FUN = function(i){
simulate(fitcov = fitcov, fitY = fitY, fitD = fitD,
ymodel_predict_custom = ymodel_predict_custom,
yrestrictions = yrestrictions,
compevent_restrictions = compevent_restrictions,
restrictions = restrictions,
outcome_name = outcome_name, compevent_name = compevent_name,
time_name = time_name,
intvars = comb_intvars[[i]], interventions = comb_interventions[[i]], int_times = comb_int_times[[i]],
histvars = histvars, histvals = histvals, histories = histories,
covparams = covparams, covnames = covnames, covtypes = covtypes,
covpredict_custom = covpredict_custom,
basecovs = basecovs,
comprisk = comprisk, ranges = ranges,
outcome_type = outcome_type,
subseed = bootseeds[r], time_points = time_points,
obs_data = resample_data, parallel = FALSE, max_visits = max_visits,
baselags = baselags, below_zero_indicator = below_zero_indicator,
min_time = min_time, show_progress = show_progress, pb = pb,
int_visit_type = int_visit_type[i], sim_trunc = sim_trunc, ...)
})
nat_pool <- pools[[1]] # Simulated data under natural course
pools <- pools[-1] # Simulated data under various interventions
if (outcome_type == 'survival'){
param <- 'poprisk'
} else if (outcome_type == 'continuous_eof'){
param <- 'Ey'
} else if (outcome_type == 'binary_eof'){
param <- 'Py'
}
# Initialize result matrix
if (grepl('eof', outcome_type)){
result_ratio <- result_diff <- int_result <- rep(NA, length(pools) + 1)
} else {
result_ratio <- result_diff <- int_result <-
matrix(NA, nrow = length(pools) + 1, ncol = time_points)
}
# Calculate mean risk at each time point under natural course
if (grepl('eof', outcome_type)){
nat_result <- mean(nat_pool[[param]], na.rm = TRUE)
} else {
nat_result <- tapply(nat_pool[[param]], nat_pool[[time_name]], FUN = mean)
}
if (ref_int == 0){
ref_result <- nat_result
} else {
if (grepl('eof', outcome_type)){
ref_result <- mean(pools[[ref_int]][[param]], na.rm = TRUE)
} else {
# Calculate mean risk at each time point for specified reference intervention
ref_result <- tapply(pools[[ref_int]][[param]], pools[[ref_int]][[time_name]], FUN = mean)
}
}
# Calculate risk ratio under natural course
# Calculate risk ratio for remaining interventions
if (grepl('eof', outcome_type)){
int_result[1] <- nat_result
int_result[-1] <- sapply(pools, FUN = function(pool){
mean(pool[[param]], na.rm = TRUE)
})
result_ratio <- int_result / ref_result
result_diff <- int_result - ref_result
} else {
int_result[1, ] <- nat_result
result_ratio[1, ] <- int_result[1, ] / ref_result
result_diff[1, ] <- int_result[1, ] - ref_result
if (length(comb_interventions) > 1){
for (i in 2:(length(pools) + 1)){
int_result[i, ] <- tapply(pools[[i - 1]][[param]], pools[[i - 1]][[time_name]], FUN = mean)
result_ratio[i, ] <- int_result[i, ] / ref_result
result_diff[i, ] <- int_result[i, ] - ref_result
}
}
}
# Calculate hazard ratio for specified interventions
if (hazardratio){
# Generate dataset containing failure/censor time information for each subject
# under each intervention
pools_hr <- lapply(seq_along(intcomp), FUN = hr_helper, intcomp = intcomp,
time_name = time_name, pools = pools)
data_hr <- rbindlist(pools_hr)
names(data_hr)[names(data_hr) == time_name] <- "t0"
names(data_hr)[names(data_hr) == outcome_name] <- "Y"
# Factor event variable
data_hr$event <- factor(data_hr$Ycomp, 0:2, labels=c("censor", "Y", "D"))
if (comprisk){
# Calculate subdistribution hazard ratio
hr_data <- survival::finegray(survival::Surv(t0, event) ~ ., data = data_hr, etype = "Y")
hr_res <- survival::coxph(survival::Surv(fgstart, fgstop, fgstatus) ~ regime, data = hr_data)
hr_res <- exp(hr_res$coefficients)
}
else {
# Calculate cause-specific hazard ratio
hr_res <- survival::coxph(formula = survival::Surv(t0, Y == "1") ~ regime, data = data_hr)
hr_res <- exp(hr_res$coefficients)
}
if (is.na(hr_res)){
warning(paste0('NA value for the hazard ratio found in bootstrap replicate ', r,
'. The analysis should likely be repeated with more parsimonious models.'))
}
} else {
hr_res <- NA
}
if (time_points > 1){
fits <- fitcov
fits[[length(fits) + 1]] <- fitY
} else {
fits <- list(fitY)
}
if (!is.na(fitD)[[1]]){
fits[[length(fits) + 1]] <- fitD
}
if (boot_diag){
bootcoeffs <- get_coeffs(fits = fits, fitD = fitD, time_points = time_points,
outcome_name = outcome_name, compevent_name = compevent_name,
covnames = covnames)
bootstderrs <- get_stderrs(fits = fits, fitD = fitD, time_points = time_points,
outcome_name = outcome_name, compevent_name = compevent_name,
covnames = covnames)
bootvcovs <- get_vcovs(fits = fits, fitD = fitD, time_points = time_points,
outcome_name = outcome_name, compevent_name = compevent_name,
covnames = covnames)
} else {
bootcoeffs <- NA
bootstderrs <- NA
bootvcovs <- NA
}
# Check NA values
if (any(is.na(int_result))){
warning(paste0('NA values for the g-form risk found in bootstrap replicate ', r,
'. The analysis should likely be repeated with more parsimonious models.'))
}
final <- list(Result = int_result, ResultRatio = result_ratio, ResultDiff = result_diff, ResultHR = hr_res,
bootcoeffs = bootcoeffs, bootstderrs = bootstderrs, bootvcovs = bootvcovs)
return (final)
}
bootstrap_helper_with_trycatch <- function(r, time_points, obs_data, bootseeds, outcome_type,
intvars, interventions, int_times, ref_int,
covparams, covnames, covtypes, covfits_custom, covpredict_custom, basecovs, histvars, histvals, histories,
ymodel, ymodel_fit_custom, ymodel_predict_custom,
yrestrictions, compevent_restrictions, restrictions,
comprisk, compevent_model,
time_name, outcome_name, compevent_name,
ranges, parallel, ncores, max_visits,
hazardratio, intcomp, boot_diag, nsimul, baselags,
below_zero_indicator, min_time, show_progress, pb,
int_visit_type, sim_trunc, ...){
tryCatch({
bootstrap_helper(
r = r, time_points = time_points, obs_data = obs_data, bootseeds = bootseeds, outcome_type = outcome_type,
intvars = intvars, interventions = interventions, int_times = int_times, ref_int = ref_int,
covparams = covparams, covnames = covnames, covtypes = covtypes, covfits_custom = covfits_custom, covpredict_custom = covpredict_custom, basecovs = basecovs, histvars = histvars, histvals = histvals, histories = histories,
ymodel = ymodel, ymodel_fit_custom = ymodel_fit_custom, ymodel_predict_custom = ymodel_predict_custom,
yrestrictions = yrestrictions, compevent_restrictions = compevent_restrictions, restrictions = restrictions,
comprisk = comprisk, compevent_model = compevent_model,
time_name = time_name, outcome_name = outcome_name, compevent_name = compevent_name,
ranges = ranges, parallel = parallel, ncores = ncores, max_visits = max_visits,
hazardratio = hazardratio, intcomp = intcomp, boot_diag = boot_diag, nsimul = nsimul, baselags = baselags,
below_zero_indicator = below_zero_indicator, min_time = min_time, show_progress = show_progress, pb = pb,
int_visit_type = int_visit_type, sim_trunc = sim_trunc, ...)
},
error = function(e){
warning(paste0('An error occured in bootstrap replicate ', r,
'. The analysis should likely be repeated with more parsimonious models. The error message encountered is:\n', e))
if (grepl('eof', outcome_type)){
result_ratio <- result_diff <- int_result <- rep(NA, length(interventions))
} else {
result_ratio <- result_diff <- int_result <-
matrix(NA, nrow = length(interventions), ncol = time_points)
}
final <- list(Result = int_result, ResultRatio = result_ratio, ResultDiff = result_diff, ResultHR = NA,
bootcoeffs = NULL, bootstderrs = NULL, bootvcovs = NULL)
final
})
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.