406 conditional power

R/gs_cp.R

@@ -0,0 +1,97 @@
+#  Copyright (c) 2024 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the gsDesign2 program.
+#
+#  gsDesign2 is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+#' @title Group sequential design - Conditional Probability under non-proportional hazards
+#' @description  \code{gs_cp()} computes conditional boundary crossing probabilities at future planned analyses
+#' for a given group sequential design assuming an interim z-statistic at a specified interim analysis.
+#' @param x An object of type \code{gsDesign2}
+#' @param i analysis at which interim z-value is given; must be from 1 to max(x$analysis$analysis)
+#' @param zi interim z-value at analysis i (scalar)
+#' @return A list with: theta -- a list containing theta0 (theta under H0, i.e., 0,) and theta1 (a vector of theta under H1)
+#'                      upper_bound -- the upper bound value return by any gs_design_ahr function
+#'                      upper_prob -- a list contains the conditional probability given zi under H0 and H1, respectively
+#'
+#' @examples
+#' library(gsDesign2)
+#' library(dplyr)
+#'
+#' # Example 1
+#' x <- gs_design_ahr(info_frac = c(.25, .75, 1), analysis_time = 36)
+#' gs_cp(x = x, i = 1, zi = 0, j = 2)
+#'
+gs_cp <- function(x, i, zi, j){
+
+  # input check
+  # Check if 'x' has an 'analysis' element and it's a matrix or data frame
+  # !! `x` should be an output from `gs_update_ahr`
+  if (!is.list(x) || !"analysis" %in% names(x) || (!is.matrix(x$analysis) && !is.data.frame(x$analysis))) {
+    stop("'x' should be a list containing an 'analysis' matrix or data frame.")
+  }
+  # Check if 'i' and 'j' are numeric and within the appropriate range
+  if (!is.numeric(i) || !is.numeric(j)) {
+    stop("'i' and 'j' should be numeric.")
+  }
+
+  if (!(1 <= i && i < j && j <= dim(x$analysis)[1])) {
+    stop("Please ensure that 1 <= i < j and j <= dim(x$analysis)[1].")
+  }
+
+  # Check the # of upper bound is equal to # of analysis
+  # ! what is efficacy is only at IA2 and FA?
+  if(length(which(x$bound$bound == "upper")) != dim(x$analysis)[1]){
+    stop("'x' should contains the same number of upper bounds as the number of analysis")
+  }
+
+  # obtain necessary information from x
+  #!! please change the variable name `t_frac` to `info_frac` for consistency of the entire pkg
+  t_frac <- x$analysis$info_frac
+  # !! the variable name `h0` is confusing, suggest to call it `info0` or `info_h0`
+  h0 <- x$analysis$info0 # under local asymptotic, assume H0 ~= H1
+
+  # default theta: under H0 and under H1
+  theta0 <- c(0, 0, 0)
+  theta <- x$analysis$theta
+
+  # compute the conditional probability under H0
+  ## ! x$bound is not required to switch from tibble to data.frame
+  eff_bound <- as.data.frame(x$bound) %>%
+    filter(bound == "upper") %>%
+    select(z)
+
+  # compute conditional power under H0
+  # prob0 <- 1 - pnorm((eff_bound[i+1, ] - sqrt(t_frac[i]/t_frac[i+1]) * zi)/sqrt((t_frac[i+1]-t_frac[i])/t_frac[i+1]))
+  prob0 <- 1 - pnorm((sqrt(t_frac[j])*eff_bound[j, ] - sqrt(t_frac[i])*zi)/sqrt(t_frac[j] - t_frac[i]))
+
+
+  # compute conditional power under H1
+  #mu_star <- sqrt(H0[i+1])*theta[i+1] + sqrt(t_frac[i]/t_frac[i+1])*(zi - sqrt(H0[i])*theta[i])
+  #sigma2_star <- H0[i+1]/H1[i+1] - (t_frac[i]/t_frac[i+1])*(H0[i]/H1[i])
+  mu_star <- sqrt(t_frac[j])*sqrt(h0[j])*theta[j] - sqrt(t_frac[i])*sqrt(h0[i])*theta[i]
+  sigma2_star <- t_frac[j] - t_frac[i]
+
+  prob1 <- 1 - pnorm((sqrt(t_frac[j])*eff_bound[j, ] - sqrt(t_frac[i])*zi - mu_star)/sqrt(sigma2_star))
+
+  # return list of results
+  output <- list(theta = list(theta0 = 0, theta1 = theta),
+                 upper_bound = eff_bound,
+                 upper_prob = list(prob0 = prob0, prob1 = prob1)
+  )
+
+  return(output)
+}

tests/testthat/test-developer-gs_cp.R

@@ -0,0 +1,106 @@
+test_that("Compare the conditional power of gsDesign and gsDesign2 under PH with efficacy only", {
+  # ------------------------------ #
+  #         parameters             #
+  # ------------------------------ #
+  alpha <- 0.025
+  beta <- 0.1
+  k <- 3
+  test.type <- 1
+  astar <- 0
+  timing <- c(0.5, 0.7)
+  sfu <- sfLDOF
+  sfupar <- c(0)
+  sfl <- sfLDOF
+  sflpar <- c(0)
+  hr <- 0.6
+  hr0 <- 1
+  eta <- 0.01
+  gamma <- c(2.5, 5, 7.5, 10)
+  R <- c(2, 2, 2, 6)
+  S <- NULL
+  T <- 18
+  minfup <- 6
+  ratio <- 1
+  lambdaC <- log(2) / 6
+
+  enroll_rate <- define_enroll_rate(duration = R, rate = gamma)
+  fail_rate <- define_fail_rate(duration = Inf, fail_rate = lambdaC, hr = hr, dropout_rate = eta)
+
+  # ------------------------------ #
+  # conditional power of gsDesign  #
+  # ------------------------------ #
+  # reference: https://keaven.github.io/gsDesign/articles/ConditionalPowerPlot.html
+  # original design
+  x0 <- gsSurv(k = k, test.type = test.type, alpha = alpha, beta = beta,
+               astar = astar, timing = timing, sfu = sfu, sfupar = sfupar,
+               sfl = sfl, sflpar = sflpar, lambdaC = lambdaC, hr = hr, hr0 = hr0,
+               eta = eta, gamma = gamma, R = R, S = S, T = T, minfup = minfup, ratio = ratio)
+
+  # update the design at IA1
+  x1 <- gsDesign(k = x0$k, test.type = x0$test.type,
+                 alpha = x0$alpha, beta = x0$beta,
+                 delta = x0$delta, delta1 = x0$delta1, delta0 = x0$delta0,
+                 n.I = c(85, x0$n.I[2:3]),
+                 maxn.IPlan = x0$n.I[3],
+                 sfu = sfu, sfupar = sfupar,
+                 sfl = sfl, sflpar = sflpar)
+  # we assume an interim p-value of 0.04, one-sided.
+  # this does not come close to the first efficacy or futility bound above. However, it is a trend in the right direction.
+  p <- 0.04
+  x2 <- gsCP(x1, i = 1, zi = -qnorm(p))
+
+  # ------------------------------ #
+  # conditional power of gsDesign2 #
+  # ------------------------------ #
+  # original design
+  y0 <- gs_power_ahr(enroll_rate = enroll_rate |> mutate(rate = rate * (x0$eNC[3] + x0$eNE[3]) / sum(R * gamma)),
+                     fail_rate = fail_rate,
+                     upper = gs_spending_bound, upar = list(sf = sfu, total_spend = alpha),
+                     lower = gs_b, lpar = rep(-Inf, 3),
+                     event = x0$n.I, analysis_time = x0$T)
+
+  # update design at IA1
+  y1 <- gs_power_ahr(enroll_rate = y0$enroll_rate,
+                     fail_rate = y0$fail_rate,
+                     upper = gs_spending_bound, upar = list(sf = sfu, total_spend = alpha,
+                                                            timing = c(85, y0$analysis$event[2:3]) / y0$analysis$event[3]
+                                                            ),
+                     lower = gs_b, lpar = rep(-Inf, 3),
+                     event = c(85, y0$analysis$event[2:3]), analysis_time = y0$analysis$time)
+
+  y2 <- gs_cp(x = y1, i = 1, zi = -qnorm(p), j = 2)
+  y3 <- gs_cp(x = y1, i = 1, zi = -qnorm(p), j = 3)
+  # ------------------------------ #
+  #       comparison              #
+  # ------------------------------ #
+  # comparison of sample size of the original design
+  expect_equal((x0$eNC + x0$eNE) |> as.vector(), y0$analysis$n, tolerance = 1e-5)
+
+  # comparison of events of the original design and updated design
+  expect_equal(x0$n.I, y0$analysis$event, tolerance = 1e-5)
+  expect_equal(x1$n.I, y1$analysis$event, tolerance = 1e-5)
+
+  # comparison of timing of the original design and updated design
+  expect_equal((x0$T) |> as.vector(), y0$analysis$time, tolerance = 1e-5)
+  expect_equal((x1$timing) |> as.vector(), y1$analysis$info_frac0, tolerance = 1e-5)
+
+  # comparison of crossing probability under H1
+  expect_equal(x0$upper$prob[, 2] |> cumsum(), y0$bound$probability, tolerance = 1e-2)
+  expect_equal(x1$upper$prob[, 2] |> cumsum(), y1$bound$probability, tolerance = 1e-2)
+
+  # comparison of crossing probability under H0
+  expect_equal(x0$upper$prob[, 1] |> cumsum(), y0$bound$probability0, tolerance = 1e-2)
+  expect_equal(x1$upper$prob[, 1] |> cumsum(), y1$bound$probability0, tolerance = 1e-2)
+
+  # comparison of conditional power
+  # theta = H0
+  expect_equal(x2$upper$prob[1, 2], y2$upper_prob$prob0, tolerance = 1e-3)
+  # ??? expect_equal(x2$upper$prob[2, 2], y3$upper_prob$prob0, tolerance = 1e-3)
+  expect_equal(x2$upper$prob[2, 2], y3$upper_prob$prob0 - y2$upper_prob$prob0, tolerance = 1e-1)
+  # theta = H1
+  expect_equal(x2$upper$prob[1, 3], y2$upper_prob$prob1, tolerance = 1e-3)
+  # ??? expect_equal(x2$upper$prob[2, 3], y3$upper_prob$prob1, tolerance = 1e-3)
+  expect_equal(x2$upper$prob[2, 3], y3$upper_prob$prob1 - y2$upper_prob$prob1, tolerance = 1e-3)
+
+  # ??? theta = IA estimated theta
+})