Bootstrap fix

SEQTaRget/DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SEQTaRget
 Type: Package
 Title: Sequential Trial Emulation
-Version: 1.3.1
+Version: 1.3.2
 Authors@R: c(person(given = "Ryan",
                     family = "O'Dea",
                     role = c("aut", "cre"),

SEQTaRget/R/internal_misc.R

@@ -3,9 +3,11 @@
 #' @keywords internal
 create.risk <- function(data, params) {
   variable <- followup <- V1 <- V2 <- NULL
-  i.value <- i.LCI <- i.UCI <- NULL
-  UCI <- LCI <- NULL
+  i.value <- i.LCI <- i.UCI <- i.SE <- NULL
+  UCI <- LCI <- SE <- NULL
   rd_lb <- rd_ub <- rr_lb <- rr_ub <- NULL
+  rd <- rd_se <- rr <- rr_se <- NULL
+
   var <- if ("inc0" %in% data[["variable"]]) "inc" else "risk"
   table <- data[, .SD[.N], by = "variable"
                 ][variable %like% var, 
@@ -14,28 +16,24 @@ create.risk <- function(data, params) {
   out <- CJ(table$variable, table$variable)[table, on = c("V2" = "variable")
                                             ][table, on = c("V1" = "variable")][V1 != V2, ]
 
-  out[, `:=` (rr = value / i.value, rd = value - i.value)
-      ][, `:=` (value = NULL, i.value = NULL)]
-
+  out[, `:=` (rr = value / i.value, rd = value - i.value)]
 
   table[, `:=` (A = sub(".*_", "", variable), 
                 Method = params@method,
                 variable = NULL)]
 
   if (all(c("LCI", "UCI") %in% names(out))) {
-    out[, `:=` (
-      rd_lb = LCI - i.LCI,
-      rd_ub = UCI - i.UCI,
-      rr_lb = LCI / i.LCI,
-      rr_ub = UCI / i.UCI
-      )][, `:=` (
-      rd_lci = pmin(rd_lb, rd_ub),
-      rd_uci = pmax(rd_lb, rd_ub),
-      rr_lci = pmin(rr_lb, rr_ub),
-      rr_uci = pmax(rr_lb, rr_ub))
-      ][, `:=` (LCI = NULL, UCI = NULL, i.LCI = NULL, i.UCI = NULL,
-                rd_lb = NULL, rd_ub = NULL, rr_lb = NULL, rr_ub = NULL)]
+    z <- qnorm(1 - (1 - params@bootstrap.CI)/2)
 
+    out[, `:=` (rd_se = sqrt(SE^2 + i.SE^2),
+                rr_se = sqrt((SE/value)^2 + (i.SE/i.value)^2))
+        ][, `:=` (rd_lci = rd - z*rd_se,
+                  rd_uci = rd + z*rd_se,
+                  rr_lci = exp(log(rr) - z*rr_se),
+                  rr_uci = exp(log(rr) + z*rr_se))
+          ][, `:=` (value = NULL, i.value = NULL, LCI = NULL, UCI = NULL, 
+                    i.LCI = NULL, i.UCI = NULL, SE = NULL, i.SE = NULL, 
+                    rd_se = NULL, rr_se = NULL)]
     setnames(out, names(out), c("A_x", "A_y", 
                                 "Risk Ratio", "Risk Differerence",
                                 "RD 95% LCI", "RD 95% UCI", "RR 95% LCI", "RR 95% UCI"))
@@ -45,6 +43,7 @@ create.risk <- function(data, params) {
     setnames(table, c("value", "LCI", "UCI"), c("Risk", "95% LCI", "95% UCI"))
     setcolorder(table, c("Method", "A", "Risk", "95% LCI", "95% UCI"))
   } else {
+    out[, `:=` (value = NULL, i.value = NULL)]
     setnames(out, names(out), c("A_x", "A_y", "Risk Ratio", "Risk Difference"))
     setnames(table, "value", "Risk")
     setcolorder(table, c("Method", "A", "Risk"))

SEQTaRget/R/internal_survival.R

@@ -6,6 +6,7 @@
 #'
 #' @keywords internal
 internal.survival <- function(params, outcome) {
+  SE <- NULL
   result <- local({
     on.exit({
       rm(list = setdiff(ls(), "result"))
@@ -115,22 +116,20 @@ 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)
+      DT.se <- rbindlist(data)[, list(SE = sd(value)), by = c("followup", "variable")]
+
       if (params@bootstrap.CI_method == "se") {
         z <- qnorm(1 - (1 - params@bootstrap.CI)/2)
-        DT <- rbindlist(data)[, list(se = sd(value)),
-                               by = c("followup", "variable")]
-
-        surv <- full$data[DT, on = c("followup", "variable")
-                          ][, `:=` (LCI = value - z*se, UCI = value + z*se)
-                            ][, se := NULL]
+        surv <- full$data[DT.se, on = c("followup", "variable")
+                          ][, `:=` (LCI = max(0, value - z*SE), UCI = min(1, value + z*SE)), by = .I]
       } 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]
+        DT.q<- rbindlist(data)[, list(LCI = quantile(value, (1 - params@bootstrap.CI)/2),
+                                   UCI = quantile(value, 1 - (1 - params@bootstrap.CI)/2)),
+                               by = c("followup", "variable")]
+
+        surv <- full$data[DT.se, on = c("followup", "variable")
+                          ][DT.q, on = c("followup", "variable")]
       }
-
     } else  surv <- full$data
     out <- list(data = surv, 
                 ce.model = if (!is.na(params@compevent)) if (params@bootstrap) c(list(full$ce.model), ce.models) else list(full$ce.model) else list())

SEQTaRget/tests/testthat/test_survival.R

@@ -13,3 +13,11 @@ test_that("Bootstrapped Survival", {
   expect_s4_class(model, "SEQoutput")
   expect_s3_class(model@survival.curve[[1]], "ggplot")
 })
+
+test_that("Bootstrapped Survival - Percentile", {
+  data <- copy(SEQdata)
+  model <- SEQuential(data, "ID", "time", "eligible", "tx_init", "outcome", list("N", "L", "P"), list("sex"),
+                      method = "ITT", options = SEQopts(km.curves = TRUE, bootstrap = TRUE, bootstrap.nboot = 2, bootstrap.CI_method = "percentile"))
+  expect_s4_class(model, "SEQoutput")
+  expect_s3_class(model@survival.curve[[1]], "ggplot")
+})