20 add bootstrap in maic unanchored tte

R/bootstrapping.R

@@ -0,0 +1,56 @@
+#' Bootstrapping for MAIC weighted hazard ratios
+#'
+#' @param intervention_data A data frame containing individual patient data from the internal IPD study including columns:
+#' treatment, time, status, centered baseline characteristics.
+#' @param centered_colnames A character vector giving the names of the covariates to use
+#'   in matching. These names must match the column names in intervention_data.
+#' @param i Index used to select a sample within \code{\link{boot}}.
+#' @param model A model formula in the form 'Surv(Time, Event==1) ~ ARM'.
+#'   Variable names need to match the corresponding columns in intervention_data.
+#' @param comparator_data A data frame containing pseudo individual patient data
+#'   from the comparator study needed to derive the relative treatment effect
+#'   including columns: time, event, weight=1, scaled_weight=1.
+#'   The outcome variables names must match intervention_data.
+#' @param min_weight A numeric value that defines the minimum weight allowed.
+#'   This value (default 0.0001) will replace weights estimated at 0 in a sample.
+#' @param trt_ext A character giving the name of the comparator treatment to be used as reference
+#'
+#' @details This function is intended to be used in conjunction with the
+#'   \code{\link{boot}} function to return the statistic to be
+#'   bootstrapped. In this case by performing MAIC weighting using
+#'   {\link{estimate_weights}} and returning a weighted hazard ratio (HR) from a
+#'   Cox proportional hazards model. This is used as the 'statistic' argument in
+#'   the boot function.
+#'
+#' @return The HR as a numeric value.
+#' @seealso \code{\link{estimate_weights}}, \code{\link{boot}}
+#' @example inst/examples/MAIC_example_analysis.R
+#' @export
+
+
+
+bootstrap_HR <- function(intervention_data, centered_colnames, i, model, comparator_data, min_weight = 0.0001, trt_ext) {
+  # create a visible binding for R CMD check
+  wt <- NULL
+
+  # Samples the data
+  bootstrap_data <- intervention_data[i, ]
+
+  # Estimates weights
+  perform_wt <- estimate_weights(data = bootstrap_data, centered_colnames = centered_colnames)
+
+
+  # Give comparator data weights of 1
+  comparator_data_wts <- comparator_data %>% dplyr::mutate(weights = 1, scaled_weights = 1)
+
+  # Add the comparator data
+  combined_data <- dplyr::bind_rows(perform_wt$data, comparator_data_wts)
+  combined_data$treatment <- stats::relevel(as.factor(combined_data$treatment), ref = trt_ext)
+
+  # set weights that are below min_weight to min_weight to avoid issues with 0 values
+  combined_data$wt <- ifelse(combined_data$weights < min_weight, min_weight, combined_data$weights)
+
+  # survival data stat
+  cox_model <- survival::coxph(model, data = combined_data, weights = wt)
+  HR <- exp(cox_model$coefficients)
+}

R/maic_unanchored_tte.R

@@ -2,6 +2,8 @@
 #'
 #' @param useWt a numeric vector of individual MAIC weights, length should the same as \code{nrow(dat)}
 #' @param dat a data frame that meet format requirements in 'Details', individual patient data (IPD) of internal trial
+#' @param dat_centered a data frame containing individual patient data from the internal IPD study including columns:
+#' treatment, time, status, centered baseline characteristics.
 #' @param dat_ext a data frame, pseudo IPD from digitized KM curve of external trial
 #' @param trt  a character string, name of the interested treatment in internal trial (real IPD)
 #' @param trt_ext character string, name of the interested comparator in external trial used to subset \code{dat_ext} (pseudo IPD)
@@ -22,7 +24,7 @@
 #' @export
 #'
 #' @examples
-maic_tte_unanchor <- function(useWt, dat, dat_ext, trt, trt_ext,
+maic_tte_unanchor <- function(useWt, dat, dat_centered, dat_ext, trt, trt_ext,
                               time_scale = "month", endpoint_name = "OS",
                               transform = "log") {
   timeUnit <- list("year" = 365.24, "month" = 30.4367, "week" = 7, "day" = 1)
@@ -32,6 +34,17 @@ maic_tte_unanchor <- function(useWt, dat, dat_ext, trt, trt_ext,
 
   res <- list()
 
+  # Bootstrap the confidence interval of the weighted HR
+
+  HR_bootstraps <- boot(
+    data = dat_centered, # intervention data
+    statistic = bootstrap_HR, # bootstrap the HR (defined in the maicplus package)
+    R = 1000, # number of bootstrap samples
+    comparator_data = dat_ext, # comparator pseudo data
+    centered_colnames = grep("_CENTERED$", names(dat_centered)), # matching variables
+    model = Surv(time, status == 1) ~ treatment # model to fit
+  )
+
   # set up IPD 'dat' with maic weights
   dat$weight <- useWt
   dat <- dat[, c("treatment", "time", "status", "weight")]
@@ -129,6 +142,30 @@ maic_tte_unanchor <- function(useWt, dat, dat_ext, trt, trt_ext,
   )
   res[["plot_resid_after"]] <- grDevices::recordPlot()
 
+
+
+  # ==> Report 4: Bootstrapping
+
+  ## Bootstrapping diagnostics
+  # Summarize bootstrap estimates in a histogram
+  # Vertical lines indicate the median and upper and lower CIs
+  hist(HR_bootstraps$t, main = "", xlab = "Boostrapped HR")
+  abline(v = quantile(HR_bootstraps$t, probs = c(0.025, 0.5, 0.975)), lty = 2)
+
+  # Median of the bootstrap samples
+  HR_median <- median(HR_bootstraps$t)
+
+  # Bootstrap CI - Percentile CI
+  boot_ci_HR <- boot.ci(boot.out = HR_bootstraps, index = 1, type = "perc")
+
+  # Bootstrap CI - BCa CI
+  boot_ci_HR_BCA <- boot.ci(boot.out = HR_bootstraps, index = 1, type = "bca")
+
+  res[["boot_hist"]] <- grDevices::recordPlot()
+  res[["HR_median"]] <- HR_median
+  res[["boot_ci_HR"]] <- boot_ci_HR
+  res[["boot_ci_HR_BCA"]] <- boot_ci_HR_BCA
+
   # output
   res
 }