Adjust the FA futility bound with same value as the efficacy bound for asymmetric beta-spending non-binding/binding designs

R/gs_power_npe.R

@@ -350,6 +350,9 @@ gs_power_npe <- function(theta = .1, theta0 = 0, theta1 = theta, # 3 theta
   upper_prob <- rep(NA, n_analysis)
   lower_prob <- rep(NA, n_analysis)
   harm_prob <- rep(NA, n_analysis)
+  # Close the continuation region at the final analysis for beta-spending designs
+  close_final_bound <- test_upper[n_analysis] && test_lower[n_analysis] &&
+    identical(lower, gs_spending_bound) && theta1[n_analysis] != theta0[n_analysis]
 
   # Calculate crossing prob under H1 ----
   for (k in 1:n_analysis) {
@@ -359,6 +362,10 @@ gs_power_npe <- function(theta = .1, theta0 = 0, theta1 = theta, # 3 theta
       theta = theta1, efficacy = FALSE
     )
     b[k] <- upper(k = k, par = upar, hgm1 = hgm1_0, info = info0, r = r, tol = tol, test_bound = test_upper)
+    # Set the final futility bound equal to the efficacy bound
+    if (k == n_analysis && close_final_bound) {
+      a[k] <- b[k]
+    }
     harm_z[k] <- harm(
       k = k, par = hpar, hgm1 = hgm1_harm0, info = info0, r = r, tol = tol, test_bound = test_harm,
       theta = theta0, efficacy = FALSE
@@ -460,6 +467,11 @@ gs_power_npe <- function(theta = .1, theta0 = 0, theta1 = theta, # 3 theta
         )
       }
     }
+    # Assign the remaining probability mass to futility when no final continuation region remains
+    if (k == n_analysis && close_final_bound) {
+      previous_lower_prob <- if (k > 1) sum(lower_prob[seq_len(k - 1)]) else 0
+      lower_prob[k] <- 1 - sum(upper_prob) - previous_lower_prob
+    }
   }
 
   if (all(!test_harm)) {

tests/testit/test-developer-gs_design_ahr.R

@@ -43,8 +43,7 @@ assert("Multiple analysisTimes", {
   (x1$analysis$theta %==% u2$theta)
   (x1$analysis$info %==% u2$info)
   (x1$analysis$info0 %==% u2$info0)
-  (l1$z %==% l2$Z)
-  (l1$probability %==% l2$Probability)
+  (head(l1$z, -1) %==% head(l2$Z, -1))
 })
 
 assert("Specified information fraction", {
@@ -62,8 +61,7 @@ assert("Specified information fraction", {
   (all.equal(x1$analysis$theta, u2$theta, tolerance = 1e-7))
   (all.equal(x1$analysis$info, u2$info, tolerance = 4e-8))
   (all.equal(x1$analysis$info0, u2$info0, tolerance = 4e-8))
-  (all.equal(l1$z, l2$Z, tolerance = 2e-7))
-  (all.equal(l1$probability, l2$Probability))
+  (all.equal(head(l1$z, -1), head(l2$Z, -1), tolerance = 2e-7))
 })
 
 assert("Multiple analysis times & IF and driven by times", {
@@ -81,8 +79,7 @@ assert("Multiple analysis times & IF and driven by times", {
   (x1$analysis$theta %==% u2$theta)
   (x1$analysis$info %==% u2$info)
   (x1$analysis$info0 %==% u2$info0)
-  (l1$z %==% l2$Z)
-  (l1$probability %==% l2$Probability)
+  (head(l1$z, -1) %==% head(l2$Z, -1))
 })
 
 assert("Multiple analysis times & IF and driven by IF", {
@@ -100,8 +97,7 @@ assert("Multiple analysis times & IF and driven by IF", {
   (x1$analysis$theta %==% u2$theta)
   (x1$analysis$info %==% u2$info)
   (x1$analysis$info0 %==% u2$info0)
-  (l1$z %==% l2$Z)
-  (l1$probability %==% l2$Probability)
+  (head(l1$z, -1) %==% head(l2$Z, -1))
 })
 
 assert("2-sided symmetric design with O'Brien-Fleming spending", {
@@ -166,8 +162,99 @@ assert("Pocock lower spending under H1 (NPH)", {
   (x1$analysis$theta %==% u2$theta)
   (x1$analysis$info %==% u2$info)
   (x1$analysis$info0 %==% u2$info0)
-  (l1$z %==% l2$Z)
-  (l1$probability %==% l2$Probability)
+  (head(l1$z, -1) %==% head(l2$Z, -1))
+})
+
+assert("Final futility and efficacy bounds match with lower spending under beta-spending, binding/non-binding when it is PH", {
+  enroll_rate <- define_enroll_rate(
+    duration = c(2, 2, 10),
+    rate = (1:3) / 3
+  )
+  fail_rate <- define_fail_rate(
+    duration = Inf,
+    fail_rate = log(2) / 9,
+    hr = 0.6,
+    dropout_rate = 0.0001
+  )
+
+  # binding
+  x <- gs_design_ahr(
+      enroll_rate = enroll_rate, fail_rate = fail_rate,
+      alpha = 0.025, beta = 0.1,
+      info_scale = "h0_h1_info",
+      info_frac = 1:3 / 3, analysis_time = 36,
+      upper = gs_spending_bound,
+      upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+      lower = gs_spending_bound,
+      lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = 3),
+      binding = TRUE, h1_spending = TRUE
+    )
+  fa_efficacy_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "upper", ]
+  fa_futility_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "lower", ]
+  (fa_efficacy_bound$z == fa_futility_bound$z)
+
+  # non-binding
+  x <- gs_design_ahr(
+      enroll_rate = enroll_rate, fail_rate = fail_rate,
+      alpha = 0.025, beta = 0.1,
+      info_scale = "h0_h1_info",
+      info_frac = 1:3 / 3, analysis_time = 36,
+      upper = gs_spending_bound,
+      upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+      lower = gs_spending_bound,
+      lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = 3),
+      binding = FALSE, h1_spending = TRUE
+    )
+  fa_efficacy_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "upper", ]
+  fa_futility_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "lower", ]
+  (fa_efficacy_bound$z == fa_futility_bound$z)
+
+})
+
+assert("Final futility and efficacy bounds match with lower spending under beta-spending, binding/non-binding when it is NPH", {
+  enroll_rate <- define_enroll_rate(
+    duration = c(2, 2, 10),
+    rate = (1:3) / 3
+  )
+  fail_rate <- define_fail_rate(
+    duration = c(3, Inf),
+    fail_rate = log(2) / 9,
+    hr = c(1, 0.6),
+    dropout_rate = 0.0001
+  )
+
+  # binding
+  x <- gs_design_ahr(
+      enroll_rate = enroll_rate, fail_rate = fail_rate,
+      alpha = 0.025, beta = 0.1,
+      info_scale = "h0_h1_info",
+      info_frac = 1:3 / 3, analysis_time = 36,
+      upper = gs_spending_bound,
+      upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+      lower = gs_spending_bound,
+      lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = 3),
+      binding = TRUE, h1_spending = TRUE
+    )
+  fa_efficacy_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "upper", ]
+  fa_futility_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "lower", ]
+  (fa_efficacy_bound$z == fa_futility_bound$z)
+
+  # non-binding
+  x <- gs_design_ahr(
+      enroll_rate = enroll_rate, fail_rate = fail_rate,
+      alpha = 0.025, beta = 0.1,
+      info_scale = "h0_h1_info",
+      info_frac = 1:3 / 3, analysis_time = 36,
+      upper = gs_spending_bound,
+      upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+      lower = gs_spending_bound,
+      lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = 3),
+      binding = FALSE, h1_spending = TRUE
+    )
+  fa_efficacy_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "upper", ]
+  fa_futility_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "lower", ]
+  (fa_efficacy_bound$z == fa_futility_bound$z)
+
 })
 
 assert("Spending time when both efficacy and futility bound are fixed", {

tests/testit/test-developer-gs_design_npe.R

@@ -216,7 +216,10 @@ assert("spending bounds", {
     dplyr::mutate(bound = tolower(bound)) |>
     dplyr::select(analysis, bound, z, probability, theta, info, info0, info1) |>
     dplyr::arrange(analysis, bound)
-  (as.data.frame(x1_c) %==% as.data.frame(x2))
+  legacy_rows <- x1_c$analysis < 3 | x1_c$bound != "lower"
+  (as.data.frame(x1_c[legacy_rows, ]) %==% as.data.frame(x2[legacy_rows, ]))
+  (x1_c$z[x1_c$analysis == 3 & x1_c$bound == "lower"] ==
+    x1_c$z[x1_c$analysis == 3 & x1_c$bound == "upper"])
 })
 
 assert("2-sided symmetric spend", {
@@ -266,7 +269,10 @@ assert("2-sided symmetric spend", {
     dplyr::mutate(bound = tolower(bound)) |>
     dplyr::select(analysis, bound, z, probability, theta, info, info0, info1) |>
     dplyr::arrange(analysis, bound)
-  (as.data.frame(x1_c) %==% as.data.frame(x2))
+  legacy_rows <- x1_c$analysis < 3 | x1_c$bound != "lower"
+  (as.data.frame(x1_c[legacy_rows, ]) %==% as.data.frame(x2[legacy_rows, ]))
+  (x1_c$z[x1_c$analysis == 3 & x1_c$bound == "lower"] ==
+    x1_c$z[x1_c$analysis == 3 & x1_c$bound == "upper"])
 })
 
 assert("Harm bound is not provided for fixed designs", {
@@ -277,4 +283,110 @@ assert("Harm bound is not provided for fixed designs", {
       lower = gs_b, lpar = -Inf, test_lower = FALSE,
       harm = gs_b, hpar = -2, test_harm = TRUE)
   ))
+})
+
+assert("Comparison with gsDesign when test.type = 4", {
+  timing <- 1:3/3
+  alpha <- 0.025
+  beta <- 0.1
+  effect <- 0.5
+  standard_deviation <- 1
+
+  fixed_normal <- nNormal(
+    delta1 = effect, sd = standard_deviation,
+    alpha = alpha, beta = beta, ratio = 1, outtype = 3)
+
+  gsdesign_normal <- gsDesign(
+    k = 3, test.type = 4, alpha = alpha, beta = beta,
+    n.fix = fixed_normal$n, timing = timing,
+    sfu = sfLDOF, sfl = sfLDOF, delta0 = 0, delta1 = effect)
+
+  fixed_information <- fixed_normal$n / (4 * standard_deviation^2)
+
+  gsdesign2_normal <- gs_design_npe(
+    theta = effect, theta0 = 0, theta1 = effect,
+    info = fixed_information * timing, info_scale = "h0_info",
+    alpha = alpha, beta = beta, binding = FALSE,
+    upper = gs_spending_bound,
+    upar = list(sf = sfLDOF, total_spend = alpha),
+    lower = gs_spending_bound,
+    lpar = list(sf = sfLDOF, total_spend = beta))
+
+    # the efficacy bounds from gsDesign match gsDesign2
+    (all.equal(
+      gsdesign2_normal$z[gsdesign2_normal$bound == "upper"],
+      gsdesign_normal$upper$bound,
+      tolerance = 7e-6,
+      scale = 1
+    ))
+
+    # the futility bounds from gsDesign match gsDesign2
+    (all.equal(
+      gsdesign2_normal$z[gsdesign2_normal$bound == "lower"],
+      gsdesign_normal$lower$bound,
+      tolerance = 1e-5,
+      scale = 1
+    ))
+
+    # the FA efficacy bound match with futility bound in gsDesign2
+    (all.equal(
+      gsdesign2_normal$z[gsdesign2_normal$bound == "upper" & gsdesign2_normal$analysis == 3],
+      gsdesign2_normal$z[gsdesign2_normal$bound == "lower" & gsdesign2_normal$analysis == 3],
+      tolerance = 1e-8,
+      scale = 1
+    ))
+
+})
+
+assert("Comparison with gsDesign when test.type = 3", {
+  timing <- 1:3/3
+  alpha <- 0.025
+  beta <- 0.1
+  effect <- 0.5
+  standard_deviation <- 1
+
+  fixed_normal <- nNormal(
+    delta1 = effect, sd = standard_deviation,
+    alpha = alpha, beta = beta, ratio = 1, outtype = 3)
+
+  gsdesign_normal <- gsDesign(
+    k = 3, test.type = 3, alpha = alpha, beta = beta,
+    n.fix = fixed_normal$n, timing = timing,
+    sfu = sfLDOF, sfl = sfLDOF, delta0 = 0, delta1 = effect)
+
+  fixed_information <- fixed_normal$n / (4 * standard_deviation^2)
+
+  gsdesign2_normal <- gs_design_npe(
+    theta = effect, theta0 = 0, theta1 = effect,
+    info = fixed_information * timing, info_scale = "h0_info",
+    alpha = alpha, beta = beta, binding = TRUE,
+    upper = gs_spending_bound,
+    upar = list(sf = sfLDOF, total_spend = alpha),
+    lower = gs_spending_bound,
+    lpar = list(sf = sfLDOF, total_spend = beta))
+
+    # the efficacy bounds from gsDesign match gsDesign2
+    (all.equal(
+      gsdesign2_normal$z[gsdesign2_normal$bound == "upper"],
+      gsdesign_normal$upper$bound,
+      tolerance = 7e-6,
+      scale = 1
+    ))
+
+    # the futility bounds from gsDesign match gsDesign2
+    (all.equal(
+      gsdesign2_normal$z[gsdesign2_normal$bound == "lower"],
+      gsdesign_normal$lower$bound,
+      tolerance = 1e-5,
+      scale = 1
+    ))
+
+    # the FA efficacy bound match with futility bound in gsDesign2
+    (all.equal(
+      gsdesign2_normal$z[gsdesign2_normal$bound == "upper" & gsdesign2_normal$analysis == 3],
+      gsdesign2_normal$z[gsdesign2_normal$bound == "lower" & gsdesign2_normal$analysis == 3],
+      tolerance = 1e-8,
+      scale = 1
+    ))
+
 })

tests/testit/test-developer-gs_power_ahr.R

@@ -44,8 +44,8 @@ assert("calendar based cut", {
   (x1$analysis$info0 %==% u2$info0)
   (x1$analysis$time %==% l2$Time)
   (x1$analysis$event %==% l2$Events)
-  (l1$z %==% l2$Z)
-  (l1$probability %==% l2$Probability)
+  (head(l1$z, -1) %==% head(l2$Z, -1))
+  (head(l1$probability, -1) %==% head(l2$Probability, -1))
   (x1$analysis$ahr %==% l2$AHR)
   (x1$analysis$theta %==% l2$theta)
   (x1$analysis$info %==% l2$info)
@@ -81,8 +81,8 @@ assert("event based cut", {
   (all.equal(x1$analysis$info0, x2$info0[x2$Bound == "Upper"], tolerance = 2e-7))
   (all.equal(x1$analysis$time, x2$Time[x2$Bound == "Lower"], tolerance = 2e-7))
   (all.equal(x1$analysis$event, x2$Events[x2$Bound == "Lower"], tolerance = 2e-7))
-  (all.equal(x1$bound$z[x1$bound$bound == "lower"], x2$Z[x2$Bound == "Lower"], tolerance = 2e-7))
-  (all.equal(x1$bound$probability[x1$bound$bound == "lower"], x2$Probability[x2$Bound == "Lower"], tolerance = 3e-7))
+  (all.equal(head(x1$bound$z[x1$bound$bound == "lower"], -1), head(x2$Z[x2$Bound == "Lower"], -1), tolerance = 1e-6, scale = 1))
+  (all.equal(head(x1$bound$probability[x1$bound$bound == "lower"], -1), head(x2$Probability[x2$Bound == "Lower"], -1), tolerance = 3e-7, scale = 1))
   (all.equal(x1$analysis$ahr, x2$AHR[x2$Bound == "Lower"], tolerance = 3e-7))
   (all.equal(x1$analysis$theta, x2$theta[x2$Bound == "Lower"], tolerance = 1e-7))
   (all.equal(x1$analysis$info, x2$info[x2$Bound == "Lower"], tolerance = 3e-6))
@@ -118,8 +118,8 @@ assert("calendar + event based cut", {
   (all.equal(x1$analysis$info0, x2$info0[x2$Bound == "Upper"], tolerance = 3e-6))
   (all.equal(x1$analysis$time, x2$Time[x2$Bound == "Lower"], tolerance = 2e-6))
   (all.equal(x1$analysis$event, x2$Events[x2$Bound == "Lower"], tolerance = 3e-6))
-  (all.equal(x1$bound$z[x1$bound$bound == "lower"], x2$Z[x2$Bound == "Lower"], tolerance = 2e-6))
-  (all.equal(x1$bound$probability[x1$bound$bound == "lower"], x2$Probability[x2$Bound == "Lower"], tolerance = 5e-7))
+  (all.equal(head(x1$bound$z[x1$bound$bound == "lower"], -1), head(x2$Z[x2$Bound == "Lower"], -1), tolerance = 2e-6))
+  (all.equal(head(x1$bound$probability[x1$bound$bound == "lower"], -1), head(x2$Probability[x2$Bound == "Lower"], -1), tolerance = 5e-7, scale = 1))
   (all.equal(x1$analysis$ahr, x2$AHR[x2$Bound == "Lower"], tolerance = 3e-7))
   (all.equal(x1$analysis$theta, x2$theta[x2$Bound == "Lower"], tolerance = 3e-7, scale = 1))
   (all.equal(x1$analysis$info, x2$info[x2$Bound == "Lower"], tolerance = 3e-6))
@@ -156,3 +156,47 @@ assert("Validate the boundary is symmetric in symmetric designs.", {
   lower_z <- x$bound$z[x$bound$bound == "lower"]
   (all.equal(upper_z, -lower_z))
 })
+
+assert("Final futility and efficacy bounds match with lower spending under beta-spending, binding/non-binding when it is NPH", {
+  enroll_rate <- define_enroll_rate(
+    duration = 12,
+    rate = 600/12
+  )
+  fail_rate <- define_fail_rate(
+    duration = c(3, Inf),
+    fail_rate = log(2) / 9,
+    hr = c(1, 0.6),
+    dropout_rate = 0.0001
+  )
+
+  # binding
+  x <- gs_power_ahr(
+      enroll_rate = enroll_rate, fail_rate = fail_rate,
+      info_scale = "h0_h1_info",
+      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::sfHSD, total_spend = 0.1, param = -3),
+      binding = TRUE, h1_spending = TRUE
+    )
+  fa_efficacy_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "upper", ]
+  fa_futility_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "lower", ]
+  (fa_efficacy_bound$z == fa_futility_bound$z)
+
+  # non-binding
+  x <- gs_power_ahr(
+      enroll_rate = enroll_rate, fail_rate = fail_rate,
+      info_scale = "h0_h1_info",
+      analysis_time = c(24, 30, 36),
+      upper = gs_spending_bound,
+      upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+      lower = gs_spending_bound,
+      lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = 3),
+      binding = FALSE, h1_spending = TRUE
+    )
+  fa_efficacy_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "upper", ]
+  fa_futility_bound <- x$bound[x$bound$analysis == 3 & x$bound$bound == "lower", ]
+  (fa_efficacy_bound$z == fa_futility_bound$z)
+
+})