Merge pull request #563 from Merck/562-gs_design_wlr-does-not-approximate-gs_design_ahr-when-it-is-logrank-test

Fix the issue where WLR design is different from AHR design even under logrank tests

Author: keaven

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: gsDesign2
 Title: Group Sequential Design with Non-Constant Effect
-Version: 1.1.5
+Version: 1.1.5.1
 Authors@R: c(
     person("Keaven", "Anderson", email = "[email protected]", role = c("aut")),
     person("Yujie", "Zhao", email = "[email protected]", role = c("aut", "cre")),

R/gs_power_wlr.R

@@ -169,6 +169,7 @@ gs_power_wlr <- function(enroll_rate = define_enroll_rate(duration = c(2, 2, 10)
                          event = c(30, 40, 50),
                          analysis_time = NULL,
                          binding = FALSE,
+                         h1_spending = TRUE,
                          upper = gs_spending_bound,
                          lower = gs_spending_bound,
                          upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
@@ -208,12 +209,23 @@ gs_power_wlr <- function(enroll_rate = define_enroll_rate(duration = c(2, 2, 10)
     interval = interval
   )
 
+  # set up H1 spending
+  if (h1_spending) {
+    theta1 <- x$theta
+    info1 <- x$info
+  } else {
+    theta1 <- 0
+    info1 <- x$info0
+  }
+
   # Given the above statistical information calculate the power ----
   y_h1 <- gs_power_npe(
     theta = x$theta,
+    theta0 = 0,
+    theta1 = theta1,
     info = x$info,
     info0 = x$info0,
-    info1 = x$info,
+    info1 = info1,
     info_scale = info_scale,
     binding = binding,
     upper = upper,
@@ -229,10 +241,10 @@ gs_power_wlr <- function(enroll_rate = define_enroll_rate(duration = c(2, 2, 10)
   y_h0 <- gs_power_npe(
     theta = 0,
     theta0 = 0,
-    theta1 = x$theta,
+    theta1 = theta1,
     info = x$info0,
     info0 = x$info0,
-    info1 = x$info,
+    info1 = info1,
     info_scale = info_scale,
     binding = binding,
     upper = upper,

R/utility_wlr.R

@@ -299,6 +299,20 @@ gs_sigma2_wlr <- function(arm0,
   p1 <- arm1$size / (arm0$size + arm1$size)
   p0 <- 1 - p1
 
+  # ----------------------- #
+  # get enroll duration and relative rate
+  # ----------------------- #
+  enroll_duration <- arm0$accr_interval |> diff()
+  enroll_total_duration <- arm0$accr_interval |> max()
+  n_enroll_piece <- length(enroll_duration)
+  enroll_relative_rate <- rep(-10, n_enroll_piece)
+  for (s in 1:n_enroll_piece) {
+    enroll_relative_rate[s] <- arm0$accr_param[s] / arm0$accr_param[n_enroll_piece] * enroll_duration[n_enroll_piece] / enroll_duration[s]
+  }
+
+  # ----------------------- #
+  #    get the weights
+  # ----------------------- #
   if (identical(weight, "logrank")) {
     weight_fun <- wlr_weight_1
   } else {
@@ -311,10 +325,34 @@ gs_sigma2_wlr <- function(arm0,
     )
   }
 
+  # ----------------------- #
+  #     derive sigma2
+  # ----------------------- #
+  # from Section 2.3.3 of Yung and Liu 2020, the WLR test statistics is U/sqrt(V)
+  # for the denominator V, by low of large numbers, we have V/n -> sigma2 in probability, where n is the number of subjects
   if (approx == "event_driven") {
+    # “event_driven" is only for logrank test of PH
     nu <- p0 * prob_event(arm0, tmax = tmax) + p1 * prob_event(arm1, tmax = tmax)
     sigma2 <- p0 * p1 * nu
   } else if (approx %in% c("asymptotic", "generalized_schoenfeld")) {
+    # if IA time is early than enroll finish, we adjust the accrual to
+    # ensure the info0 is approximately proportional to events when it is logrnak test
+    if (tmax < enroll_total_duration) {
+      arm0 <- arm0
+      arm1 <- arm1
+
+      # adjust the accrual time to the time of analysis
+      arm0$accr_time <- tmax
+      arm1$accr_time <- tmax
+      # truncate the accrual interval to be stop at the time of analysis
+      arm0$accr_interval <- c(tmax, arm0$accr_interval)[which(c(tmax, arm0$accr_interval) <= tmax)] |> sort()
+      arm1$accr_interval <- arm0$accr_interval
+      truncated_enroll_duration <- diff(arm0$accr_interval)
+      # adjust the accrual probability per interval by the truncated interval
+      arm0$accr_param <- truncated_enroll_duration * enroll_relative_rate[1:length(truncated_enroll_duration)] / sum(truncated_enroll_duration * enroll_relative_rate[1:length(truncated_enroll_duration)])
+      arm1$accr_param <- arm0$accr_param
+    }
+
     sigma2 <- stats::integrate(
       function(x) {
         weight_fun(x, arm0, arm1)^2 *

man/gs_power_wlr.Rd

@@ -217,3 +217,93 @@ test_that("Validate if the output info-frac match the planned info-frac, when th
   expect_equal(x2$analysis$info_frac0[1], 0.75, tolerance = 5e-6)
   expect_equal(x3$analysis$info_frac[1], 0.75, tolerance = 5e-6)
 })
+
+test_that("Validate if WLR design under logrank test generates similar design as in AHR", {
+  x1 <- gs_design_ahr(
+    enroll_rate = enroll_rate,
+    fail_rate = fail_rate,
+    alpha = 0.025,
+    beta = 0.1,
+    analysis_time = c(9, 27, 36),
+    info_frac = NULL,
+    ratio = ratio,
+    binding = FALSE,
+    upper = "gs_spending_bound",
+    upar = list(sf = "sfLDOF", total_spend = 0.025, param = NULL),
+    lower = "gs_spending_bound",
+    lpar = list(sf = "sfLDOF", total_spend = 0.1, param = NULL),
+    h1_spending = TRUE,
+    info_scale = "h0_h1_info")
+
+  x2 <- gs_design_wlr(
+    enroll_rate = enroll_rate,
+    fail_rate = fail_rate,
+    alpha = 0.025,
+    beta = 0.1,
+    analysis_time = c(9, 27, 36),
+    info_frac = NULL,
+    ratio = ratio,
+    binding = FALSE,
+    weight = "logrank",
+    upper = "gs_spending_bound",
+    upar = list(sf = "sfLDOF", total_spend = 0.025, param = NULL),
+    lower = "gs_spending_bound",
+    lpar = list(sf = "sfLDOF", total_spend = 0.1, param = NULL),
+    h1_spending = TRUE,
+    info_scale = "h0_h1_info")
+
+  # sample size is approximately close
+  expect_equal(x1$analysis$n, x2$analysis$n, tolerance = 1e-2)
+
+  # events is approximately close
+  expect_equal(x1$analysis$event, x2$analysis$event, tolerance = 1e-2)
+
+  # ahr is approximately close
+  expect_equal(x1$analysis$ahr, x2$analysis$ahr, tolerance = 1e-2)
+
+  # theta is approximately close
+  expect_equal(x1$analysis$theta, x2$analysis$theta, tolerance = 5e-2)
+
+  # info is approximately close
+  expect_equal(x1$analysis$info, x2$analysis$info, tolerance = 1e-2)
+
+  # info0 is approximately close
+  expect_equal(x1$analysis$info0, x2$analysis$info0, tolerance = 1e-2)
+
+  # info_frac is approximately close
+  expect_equal(x1$analysis$info_frac, x2$analysis$info_frac, tolerance = 1e-2)
+
+  # info_frac0 is approximately close
+  expect_equal(x1$analysis$info_frac0, x2$analysis$info_frac0, tolerance = 1e-2)
+
+  # upper bound is approximately close
+  expect_equal(x1$bound$z[x1$bound$bound == "upper"], x2$bound$z[x2$bound$bound == "upper"], tolerance = 1e-2)
+
+  # lower bound is approximately close
+  expect_equal(x1$bound$z[x1$bound$bound == "lower"], x2$bound$z[x2$bound$bound == "lower"], tolerance = 5e-2)
+
+  # probability of crossing upper bound under H0 is approximately close
+  expect_equal(x1$bound$probability0[x1$bound$bound == "upper"], x2$bound$probability0[x2$bound$bound == "upper"], tolerance = 1e-2)
+
+  # probability of crossing lower bound under H0 is approximately close
+  expect_equal(x1$bound$probability0[x1$bound$bound == "lower"], x2$bound$probability0[x2$bound$bound == "lower"], tolerance = 1e-2)
+
+  # probability of crossing upper bound under H1 is approximately close
+  expect_equal(x1$bound$probability[x1$bound$bound == "upper"], x2$bound$probability[x2$bound$bound == "upper"], tolerance = 1e-2)
+
+  # probability of crossing lower bound under H1 is approximately close
+  expect_equal(x1$bound$probability[x1$bound$bound == "lower"], x2$bound$probability[x2$bound$bound == "lower"], tolerance = 1e-2)
+
+  # hr at upper bound is approximately close
+  expect_equal(x1$bound$`~hr at bound`[x1$bound$bound == "upper"], x2$bound$`~hr at bound`[x2$bound$bound == "upper"], tolerance = 1e-2)
+
+  # hr at lower bound is approximately close
+  expect_equal(x1$bound$`~hr at bound`[x1$bound$bound == "lower"], x2$bound$`~hr at bound`[x2$bound$bound == "lower"], tolerance = 5e-2)
+
+  # nominal p at upper bound is approximately close
+  expect_equal(x1$bound$`nominal p`[x1$bound$bound == "upper"], x2$bound$`nominal p`[x2$bound$bound == "upper"], tolerance = 1e-2)
+
+  # nominal p at lower bound is approximately close
+  expect_equal(x1$bound$`nominal p`[x1$bound$bound == "upper"], x2$bound$`nominal p`[x2$bound$bound == "upper"], tolerance = 1e-2)
+
+})

tests/testthat/test-developer-gs_design_wlr.R

@@ -0,0 +1,27 @@
+test_that("statistcial information of WLR under logrank test is approximately same as that from AHR", {
+  enroll_rate <- define_enroll_rate(duration = c(2, 2, 2, 18), rate = c(1, 2, 3, 4), stratum = "All")
+
+  duration <- diff(c(0, 4, 6, 44))
+  control_rate <- log(2)/c(6, 6, 6)
+  fail_rate <- define_fail_rate(duration, fail_rate = control_rate,
+                                dropout_rate = c(0.001, 0.001, 0.001),
+                                hr = c(1, 0.8, 0.6))
+
+  x1 <- gs_info_ahr(enroll_rate = enroll_rate,
+                    fail_rate = fail_rate,
+                    ratio = 1,
+                    event = NULL,
+                    analysis_time = c(9, 27, 36))
+
+  x2 <- gs_info_wlr(enroll_rate = enroll_rate,
+                    fail_rate = fail_rate,
+                    ratio = 1,
+                    event = NULL,
+                    analysis_time = c(9, 27, 36),
+                    weight = "logrank")
+
+  # info0 of AHR and WLR is approximately close
+  expect_equal(x1$info0, x2$info0, tolerance = 1e-2)
+  # info0 of WLR is approximately proportional to events
+  expect_equal(x2$event / x2$info0, c(4, 4, 4), tolerance = 1e-2)
+})

tests/testthat/test-developer-gs_info_wlr.R

@@ -0,0 +1,13 @@
+test_that("Validate 2-sided symetric design", {
+  x <- gs_power_wlr(enroll_rate = define_enroll_rate(duration = 12, rate = 50),
+                    analysis_time = NULL, event = c(100, 200, 300),
+                    upper = gs_spending_bound,
+                    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+                    lower = gs_spending_bound,
+                    lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+                    binding = TRUE, h1_spending = FALSE)
+
+  expect_equal(x$bound$z[x$bound$bound == "upper"],
+               -x$bound$z[x$bound$bound == "lower"])
+
+})