Merge pull request #499 from Merck/498-allow-event-table-as-input-of-gs_update_ahr

Allow `event_tbl` as input of `gs_update_ahr`

Author: keaven

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: gsDesign2
 Title: Group Sequential Design with Non-Constant Effect
-Version: 1.1.3
+Version: 1.1.3.1
 Authors@R: c(
     person("Keaven", "Anderson", email = "keaven_anderson@merck.com", role = c("aut")),
     person("Yilong", "Zhang", email = "elong0527@gmail.com", role = c("aut")),

R/gs_update_ahr.R

@@ -25,10 +25,18 @@
 #' with the spending time at each analysis.
 #' @param lstime Default is NULL in which case lower bound spending time is determined by timing.
 #' Otherwise, this should be a vector of length k (total number of analyses)
-#' with the spending time at each analysis
-#' @param observed_data a list of observed datasets by analyses.
+#' with the spending time at each analysis.
+#' @param event_tbl A data frame with two columns: (1) analysis and (2) event,
+#' which represents the events observed at each analysis per piecewise interval.
+#' This can be defined via the `pw_observed_event()` function or manually entered.
+#' For example, consider a scenario with two intervals in the piecewise model:
+#' the first interval lasts 6 months with a hazard ratio (HR) of 1,
+#' and the second interval follows with an HR of 0.6.
+#' The data frame `event_tbl = data.frame(analysis = c(1, 1, 2, 2), event = c(30, 100, 30, 200))`
+#' indicates that 30 events were observed during the delayed effect period,
+#' 130 events were observed at the IA, and 230 events were observed at the FA.
 #'
-#' @return A list with input parameters, enrollment rate, analysis, and bound.
+#' @return A list with input parameters, enrollment rate, failure rate, analysis, and bound.
 #'
 #' @export
 #'
@@ -58,85 +66,7 @@
 #' ratio <- 1
 #'
 #' # ------------------------------------------------- #
-#' # Example A: one-sided design (efficacy only)
-#' # ------------------------------------------------- #
-#' # Original design
-#' upper <- gs_spending_bound
-#' upar <- list(sf = sfLDOF, total_spend = alpha)
-#' x <- gs_design_ahr(
-#'   enroll_rate = enroll_rate, fail_rate = fail_rate,
-#'   alpha = alpha, beta = beta, ratio = ratio,
-#'   info_scale = "h0_info",
-#'   info_frac = NULL,
-#'   analysis_time = c(20, 36),
-#'   upper = gs_spending_bound, upar = upar,
-#'   lower = gs_b, lpar = rep(-Inf, 2),
-#'   test_upper = TRUE, test_lower = FALSE) |> to_integer()
-#'
-#' # Observed dataset at IA and FA
-#' set.seed(123)
-#'
-#' observed_data <- simtrial::sim_pw_surv(
-#'   n = x$analysis$n[x$analysis$analysis == 2],
-#'   stratum = data.frame(stratum = "All", p = 1),
-#'   block = c(rep("control", 2), rep("experimental", 2)),
-#'   enroll_rate = x$enroll_rate,
-#'   fail_rate = (fail_rate |> simtrial::to_sim_pw_surv())$fail_rate,
-#'   dropout_rate = (fail_rate |> simtrial::to_sim_pw_surv())$dropout_rate)
-#'
-#' observed_data_ia <- observed_data |> simtrial::cut_data_by_date(x$analysis$time[1])
-#' observed_data_fa <- observed_data |> simtrial::cut_data_by_date(x$analysis$time[2])
-#'
-#' observed_event_ia <- sum(observed_data_ia$event)
-#' observed_event_fa <- sum(observed_data_fa$event)
-#'
-#' planned_event_ia <- x$analysis$event[1]
-#' planned_event_fa <- x$analysis$event[2]
-#'
-#' # Example A1 ----
-#' # IA spending = observed events / final planned events
-#' # the remaining alpha will be allocated to FA.
-#' ustime <- c(observed_event_ia / planned_event_fa, 1)
-#' gs_update_ahr(
-#'   x = x,
-#'   ustime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
-#'
-#' # Example A2 ----
-#' # IA, FA spending = observed events / final planned events
-#' ustime <- c(observed_event_ia, observed_event_fa) / planned_event_fa
-#' gs_update_ahr(
-#'   x = x,
-#'   ustime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
-#'
-#' # Example A3 ----
-#' # IA spending = min(observed events, planned events) / final planned events
-#  # the remaining alpha will be allocated to FA.
-#' ustime <- c(min(observed_event_ia, planned_event_ia) / planned_event_fa, 1)
-#' gs_update_ahr(
-#'   x = x,
-#'   ustime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
-#'
-#' # Example A4 ----
-#' # IA spending = min(observed events, planned events) / final planned events
-#' ustime <- c(min(observed_event_ia, planned_event_ia),
-#'             min(observed_event_fa, planned_event_fa)) / planned_event_fa
-#' gs_update_ahr(
-#'   x = x,
-#'   ustime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
-#'
-#' # alpha is upadted to 0.05
-#' gs_update_ahr(
-#'   x = x,
-#'   alpha = 0.05,
-#'   ustime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
-#'
-#' # ------------------------------------------------- #
-#' # Example B: Two-sided asymmetric design,
+#' # Two-sided asymmetric design,
 #' # beta-spending with non-binding lower bound
 #' # ------------------------------------------------- #
 #' # Original design
@@ -153,83 +83,48 @@
 #'   test_lower = c(TRUE, FALSE),
 #'   binding = FALSE) |> to_integer()
 #'
-#' # Example B1 ----
-#' # IA spending = observed events / final planned events
-#' # the remaining alpha will be allocated to FA.
-#' ustime <- c(observed_event_ia / planned_event_fa, 1)
-#' gs_update_ahr(
-#'   x = x,
-#'   ustime = ustime,
-#'   lstime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
+#' planned_event_ia <- x$analysis$event[1]
+#' planned_event_fa <- x$analysis$event[2]
 #'
-#' # Example B2 ----
-#' # IA, FA spending = observed events / final planned events
-#' ustime <- c(observed_event_ia, observed_event_fa) / planned_event_fa
-#' gs_update_ahr(
-#'   x = x,
-#'   ustime = ustime,
-#'   lstime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
 #'
-#' # Example B3 ----
-#' ustime <- c(min(observed_event_ia, planned_event_ia) / planned_event_fa, 1)
+#' # Updated design with 190 events observed at IA,
+#' # where 50 events observed during the delayed effect.
+#' # IA spending = observed events / final planned events, the remaining alpha will be allocated to FA.
 #' gs_update_ahr(
 #'   x = x,
-#'   ustime = ustime,
-#'   lstime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
+#'   ustime = c(190 / planned_event_fa, 1),
+#'   lstime = c(190 / planned_event_fa, 1),
+#'   event_tbl = data.frame(analysis = c(1, 1),
+#'                          event = c(50, 140)))
 #'
-#' # Example B4 ----
-#' # IA spending = min(observed events, planned events) / final planned events
-#' ustime <- c(min(observed_event_ia, planned_event_ia),
-#'             min(observed_event_fa, planned_event_fa)) / planned_event_fa
+#' # Updated design with 190 events observed at IA, and 300 events observed at FA,
+#' # where 50 events observed during the delayed effect.
+#' # IA spending = observed events / final planned events, the remaining alpha will be allocated to FA.
 #' gs_update_ahr(
 #'   x = x,
-#'   ustime = ustime,
-#'   lstime = ustime,
-#'   observed_data = list(observed_data_ia, observed_data_fa))
+#'   ustime = c(190 / planned_event_fa, 1),
+#'   lstime = c(190 / planned_event_fa, 1),
+#'   event_tbl = data.frame(analysis = c(1, 1, 2, 2),
+#'                          event = c(50, 140, 50, 250)))
 #'
-#' # Example B5 ----
-#' # alpha is updated to 0.05 ----
-#' gs_update_ahr(x = x, alpha = 0.05)
-#'
-#' # Example B6 ----
-#' # updated boundaries only when IA data is observed
-#' ustime <- c(observed_event_ia / planned_event_fa, 1)
+#' # Updated design with 190 events observed at IA, and 300 events observed at FA,
+#' # where 50 events observed during the delayed effect.
+#' # IA spending = minimal of planned and actual information fraction spending
 #' gs_update_ahr(
 #'   x = x,
-#'   ustime = ustime,
-#'   lstime = ustime,
-#'   observed_data = list(observed_data_ia, NULL))
-#'
-#' # ------------------------------------------------- #
-#' # Example C: Two-sided asymmetric design,
-#' # with calendar spending for efficacy and futility bounds
-#' # beta-spending with non-binding lower bound
-#' # ------------------------------------------------- #
-#' # Original design
-#' x <- gs_design_ahr(
-#'   enroll_rate = enroll_rate, fail_rate = fail_rate,
-#'   alpha = alpha, beta = beta, ratio = ratio,
-#'   info_scale = "h0_info",
-#'   info_frac = NULL, analysis_time = c(20, 36),
-#'   upper = gs_spending_bound,
-#'   upar = list(sf = sfLDOF, total_spend = alpha, timing = c(20, 36) / 36),
-#'   test_upper = TRUE,
-#'   lower = gs_spending_bound,
-#'   lpar = list(sf = sfLDOF, total_spend = beta, timing = c(20, 36) / 36),
-#'   test_lower = c(TRUE, FALSE),
-#'   binding = FALSE) |> to_integer()
+#'   ustime = c(min(190, planned_event_ia) / planned_event_fa, 1),
+#'   lstime = c(min(190, planned_event_ia) / planned_event_fa, 1),
+#'   event_tbl = data.frame(analysis = c(1, 1, 2, 2),
+#'                          event = c(50, 140, 50, 250)))
 #'
-#' # Updated design due to potential change of multiplicity graph
+#' # Alpha is updated to 0.05
 #' gs_update_ahr(x = x, alpha = 0.05)
 gs_update_ahr <- function(
     x = NULL,
     alpha = NULL,
     ustime = NULL,
     lstime = NULL,
-    observed_data = NULL) {
+    event_tbl = NULL) {
 
   # ----------------------------------- #
   #         Check inputs                #
@@ -277,10 +172,9 @@ gs_update_ahr <- function(
   #       At design stage,              #
   #       with different alpha          #
   # ----------------------------------- #
-  # If users do not input observed data
-  # which means they are still at the design stage
-  # but with different alpha
-  if (is.null(observed_data)) {
+  # If users do not input observed data, nor event_tbl
+  # which means they update design with a different value of alpha
+  if (is.null(event_tbl)) {
 
     # Check if ustime and lstime matches the spending time of the original design
     if (!is.null(ustime) && any(ustime != x$input$upar$timing)) {
@@ -337,24 +231,11 @@ gs_update_ahr <- function(
     #       At analysis stage,            #
     #       with different alpha          #
     # ----------------------------------- #
-    # Get the piecewise exp model for the failure rates
-    fr_duration <- x$input$fail_rate$duration
-    fr_hr <- x$input$fail_rate$hr
-    all_t <- sort(c(fr_duration, x$analysis$time))
-
-    if (is.infinite(max(x$input$fail_rate$duration))) {
-      hr_interval <- cumsum(c(fr_duration[-length(fr_duration)], max(x$analysis$time) + 50))
-    } else {
-      hr_interval <- cumsum(fr_duration)
-    }
-
-    pw_hr <- stepfun(x = hr_interval, y = c(fr_hr, last(fr_hr)), right = TRUE)
-
     # Calculate the blinded estimation of AHR
     blinded_est <- NULL
     observed_event <- NULL
     for (i in 1:n_analysis) {
-      if (is.null(observed_data[[i]])) {
+      if (!(i %in% event_tbl$analysis)) {
         # if there is no observed data at analysis i,
         # for example, we only observed IA data and FA data is unavailable yet
         blinded_est_new <- data.frame(event = x$analysis$event[i],
@@ -363,14 +244,15 @@ gs_update_ahr <- function(
                                       info0 = x$analysis$info0[i])
         event_new <- x$analysis$event[i]
       } else {
-        # if there is observed data at analysis i,
-        # we calculate the blinded estimation
-        blinded_est_new <- ahr_blinded(surv = survival::Surv(time = observed_data[[i]]$tte,
-                                                             event = observed_data[[i]]$event),
-                                       intervals = all_t[all_t <= x$analysis$time[i]],
-                                       hr = pw_hr(all_t[all_t <= x$analysis$time[i]]),
-                                       ratio = x$input$ratio)
-        event_new <- sum(observed_data[[i]]$event)
+        q_e <- x$input$ratio / (1 + x$input$ratio)
+        event_i <- event_tbl$event[event_tbl$analysis == i]
+        hr_i <- x$fail_rate$hr
+        event_new <- sum(event_i)
+
+        blinded_est_new <- data.frame(event = sum(event_i),
+                                      theta = -sum(log(hr_i) * event_i) / sum(event_i),
+                                      info0 = sum(event_i) * (1 - q_e) * q_e)
+        blinded_est_new$ahr <- exp(-blinded_est_new$theta)
       }
 
       blinded_est <- rbind(blinded_est, blinded_est_new)
@@ -441,37 +323,37 @@ gs_update_ahr <- function(
     analysis = 1:n_analysis,
     time = x$analysis$time,
     n = x$analysis$n,
-    event = if (is.null(observed_data)) {
+    event = if (is.null(event_tbl)) {
       x$analysis$event
     } else {
       observed_event
     },
-    ahr = if (is.null(observed_data)) {
+    ahr = if (is.null(event_tbl)) {
       x$analysis$ahr
     } else {
       exp(-blinded_est$theta)
     },
-    theta = if (is.null(observed_data)) {
+    theta = if (is.null(event_tbl)) {
       x$analysis$theta
     } else {
       blinded_est$theta
     },
-    info = if (is.null(observed_data)) {
+    info = if (is.null(event_tbl)) {
       x$analysis$info
     } else {
       blinded_est$info0
     },
-    info0 = if (is.null(observed_data)) {
+    info0 = if (is.null(event_tbl)) {
       x$analysis$info0
     } else {
       blinded_est$info0
     },
-    info_frac = if (is.null(observed_data)) {
+    info_frac = if (is.null(event_tbl)) {
       x$analysis$info_frac
     } else {
       upar_update$timing
     },
-    info_frac0 = if (is.null(observed_data)) {
+    info_frac0 = if (is.null(event_tbl)) {
       x$analysis$info_frac0
     } else {
       observed_event / max(observed_event)