Bootstrap CI

SEQTaRget/DESCRIPTION

@@ -1,11 +1,11 @@
 Package: SEQTaRget
 Type: Package
 Title: Sequential Trial Emulation
-Version: 1.2.0
+Version: 1.3.0
 Authors@R: c(person(given = "Ryan",
                     family = "O'Dea",
                     role = c("aut", "cre"),
-                    email = "ryanodea@hsph.harvard.edu",
+                    email = "ryan.odea@psi.ch",
                     comment = c(ORCID="0009-0000-0103-9546")),
              person(given = "Alejandro",
                     family = "Szmulewicz",
@@ -35,7 +35,7 @@ Description: Implementation of sequential trial emulation for the analysis of ob
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Suggests: 
     rmarkdown,
     testthat (>= 3.0.0)

SEQTaRget/R/SEQopts.R

@@ -1,6 +1,8 @@
 #' Parameter Builder for SEQuential Model and Estimates
 #'
 #' @param bootstrap Logical: defines if SEQuential should run bootstrapping, default is FALSE
+#' @param bootstrap.CI Numeric: defines the confidence interval after bootstrapping, default is 0.95 (95\% CI)
+#' @param bootstrap.CI_method Character: selects which way to calculate bootstraps confidence intervals ("se", "percentile")
 #' @param bootstrap.nboot Integer: number of bootstraps
 #' @param bootstrap.sample Numeric: percentage of data to use when bootstrapping, should in [0, 1], default is 0.8
 #' @param cense String: column name for additional censoring variable, e.g. loss-to-follow-up
@@ -53,13 +55,12 @@
 #' @param weight.preexpansion Logical: whether weighting should be done on pre-expanded data
 #' @param weight.upper Numeric: weights truncated at upper end at this weight
 #' @param weighted Logical: whether or not to preform weighted analysis, default is FALSE
-#'
+#' @returns An object of class 'SEQopts'
 #' @export
 #' @importFrom stats runif
 #' @importFrom parallelly availableCores
 #' @import data.table
-#' @returns An object of class 'SEQopts'
-SEQopts <- function(bootstrap = FALSE, bootstrap.nboot = 100, bootstrap.sample = 0.8,
+SEQopts <- function(bootstrap = FALSE, bootstrap.nboot = 100, bootstrap.sample = 0.8, bootstrap.CI = 0.95, bootstrap.CI_method = "se",
                     cense = NA, cense.denominator = NA, cense.eligible = NA, cense.numerator = NA,
                     compevent = NA, covariates = NA, data.return = FALSE, denominator = NA,
                     deviation = FALSE, deviation.col = NA, deviation.conditions = c(NA, NA), deviation.excused = FALSE, deviation.excused_cols = c(NA, NA),
@@ -144,6 +145,8 @@ SEQopts <- function(bootstrap = FALSE, bootstrap.nboot = 100, bootstrap.sample =
       bootstrap = bootstrap,
       bootstrap.nboot = bootstrap.nboot,
       bootstrap.sample = bootstrap.sample,
+      bootstrap.CI = bootstrap.CI,
+      bootstrap.CI_method = bootstrap.CI_method,
       seed = seed,
       followup.min = followup.min,
       followup.max = followup.max,

SEQTaRget/R/class_definitions.R

@@ -14,6 +14,8 @@ setClass("SEQopts",
     bootstrap.nboot = "integer",
     bootstrap = "logical",
     bootstrap.sample = "numeric",
+    bootstrap.CI_method = "character",
+    bootstrap.CI = "numeric",
     seed = "integer",
     followup.min = "numeric",
     followup.max = "numeric",
@@ -71,6 +73,8 @@ setClass("SEQopts",
     ncores = parallelly::availableCores(),
     bootstrap = FALSE,
     bootstrap.sample = 0.8,
+    bootstrap.CI_method = "se",
+    bootstrap.CI = 0.95,
     seed = 1636L,
     followup.min = -Inf,
     followup.max = Inf,

SEQTaRget/R/class_setters.R

@@ -24,6 +24,8 @@ parameter.setter <- function(data, DT,
     bootstrap.nboot = opts@bootstrap.nboot,
     bootstrap = opts@bootstrap,
     bootstrap.sample = opts@bootstrap.sample,
+    bootstrap.CI = opts@bootstrap.CI,
+    bootstrap.CI_method = opts@bootstrap.CI_method,
     seed = opts@seed,
     followup.include = opts@followup.include,
     trial.include = opts@trial.include,
@@ -124,6 +126,8 @@ parameter.simplifier <- function(params) {
   if (params@multinomial & params@method == "dose-response") stop("Multinomial dose-response is not supported")
 
   if (params@excused & params@deviation.excused) stop("Must select either excused from deviation or excused from treatment swap")
+
+  if (!params@bootstrap.CI_method %in% c("se", "percentile")) stop("Invalid confidence interval method defined")
 
   return(params)
 }

SEQTaRget/R/internal_hazard.R

@@ -87,8 +87,13 @@ internal.hazard <- function(model, params) {
     bootstrap <- unlist(bootstrap)
     if (all(is.na(bootstrap))) return(c(Hazard = NA_real_, LCI = NA_real_, UCI = NA_real_))
 
-    se <- sd(bootstrap, na.rm = TRUE) / sqrt(sum(!is.na(bootstrap)))
-    ci <- sort(c(full + se, full - se), decreasing = FALSE)
+    if (params@bootstrap.CI_method == "se") {
+      z <- qnorm(1 - (1 - params@bootstrap.CI)/2)
+      se <- sd(bootstrap, na.rm = TRUE) / sqrt(sum(!is.na(bootstrap)))
+      ci <- sort(c(full + z*se, full - z*se), decreasing = FALSE) 
+    } else ci <- quantile(bootstrap, 
+                          probs = c((1 - params@bootstrap.CI)/2, 
+                                    1 - (1 - params@bootstrap.CI)/2))
   } else {
     ci <- c(NA_real_, NA_real_)
   }

SEQTaRget/R/internal_survival.R

@@ -115,13 +115,21 @@ internal.survival <- function(params, outcome) {
       }
       data <- lapply(seq_along(result), function(x) result[[x]]$data)
       ce.models <- lapply(seq_along(result), function(x) result[[x]]$ce.model)
+      if (params@bootstrap.CI_method == "se") {
+        z <- qnorm(1 - (1 - params@bootstrap.CI)/2)
+        DT <- rbindlist(data)[, list(se = sd(value) / sqrt(params@bootstrap.nboot)),
+                               by = c("followup", "variable")]
 
-      DT <- rbindlist(data)[, list(se = sd(value) / sqrt(params@bootstrap.nboot)),
-                             by = c("followup", "variable")]
-
-      surv <- full$data[DT, on = c("followup", "variable")
-                        ][, `:=` (LCI = value - se, UCI = value + se)
-                          ][, se := NULL]
+        surv <- full$data[DT, on = c("followup", "variable")
+                          ][, `:=` (LCI = value - z*se, UCI = value + z*se)
+                            ][, se := NULL]
+      } else {
+        DT <- rbindlist(data)
+        surv <- full$data[DT, on = c("followup", "variable")
+                          ][, `:=` (LCI = quantile(value, (1 - params@bootstrap.CI)/2),
+                                    UCI = quantile(value, 1 - (1 - params@bootstrap.CI)/2))
+                            ][, se := NULL]
+      }
 
     } else  surv <- full$data
     out <- list(data = surv,