Calculate sequential p-values for AHR objects (#612)

* Development: 
   - Add a new function to calculate the sequential p-values for AHR objects
   - Add developer tests
   - Update the Rd file of `sequential_pval`

* Fix CI/CD checking errors by (1) adding `sequential_pval` to _pkgdown.yml, (2) reducing the long row in the example, (3) correcting the code by removing the call of `x` to `gs_design`, and (4) adding namespace to `gsBound1`

* Review, Co-authored-by: John Blischak <jdblischak@gmail.com>
   - Update R/sequential_pval.R , improve the decision of binding and non-binding
   - Update R/sequential_pval.R, apply `get_sf` 
   - Update R/sequential_pval.R, remove `@keywords internal` 
   - Update R/sequential_pval.R, apply `get_sf`
   - remove namespace of `gsDesign2::`

Author: LittleBeannie

NAMESPACE

@@ -52,6 +52,7 @@ export(gs_update_ahr)
 export(ppwe)
 export(pw_info)
 export(s2pwe)
+export(sequential_pval)
 export(text_summary)
 export(to_integer)
 export(wlr_weight_1)

R/sequential_pval.R

@@ -0,0 +1,154 @@
+#' @title Sequential p-value computation for gsDesign2
+#' @description \code{sequential_pval} computes a sequential p-value for a group sequential design
+#' created with \code{gs_design_ahr()} using a spending function approach.
+#' It is the minimum of repeated p-values computed at each analysis (Jennison and Turnbull, 2000).
+#' This is particularly useful for multiplicity methods such as the graphical method for group sequential designs
+#' where sequential p-values for multiple hypotheses can be used as nominal p-values to plug into a multiplicity graph.
+#' A sequential p-value is described as the minimum alpha level at which a one-sided group sequential bound would
+#' be rejected given interim and final observed results.
+#' Mild restrictions are required on spending functions used, but these are satisfied for commonly used spending functions
+#' such as the Lan-DeMets spending function approximating an O'Brien-Fleming bound or a Hwang-Shih-DeCani spending function.
+#'
+#' @param gs_design Group sequential design generated by \code{gs_design_ahr()}.
+#' @param event Event counts; numeric vector with increasing, positive values with at most one
+#' value greater than or equal to largest value in \code{gs_design$analysis$event};
+#' NOTE: if NULL, planned \code{event} will be used (\code{gs_design$analysis$event})
+#' @param z z-value tests corresponding to analyses in \code{info}; positive values indicate a positive finding; must have the same length as \code{info}.
+#' @param ustime Spending time for upper bound at specified analyses; specify default: \code{NULL} if this is to be based on information fraction;
+#' if not \code{NULL}, must have the same length as \code{info}; increasing positive values with at most 1 greater than or equal to 1.
+#' @param interval Interval for search to derive p-value; Default: \code{c(1e-05, 0.9999)}. Lower end of interval must be >0 and upper end must be < 1.
+#' The primary reason to not use the defaults would likely be if a test were vs a Type I error <0.0001.
+#' @return Sequential p-value (single numeric one-sided p-value between 0 and 1).
+#' Note that if the sequential p-value is less than the lower end of the input interval,
+#' the lower of interval will be returned.
+#' Similarly, if the sequential p-value is greater than the upper end of the input interval,
+#' then the upper end of interval is returned.
+#' @references
+#' Jennison C and Turnbull BW (2000), \emph{Group Sequential
+#' Methods with Applications to Clinical Trials}. Boca Raton: Chapman and Hall.
+#'
+#' Liu, Qing, and Keaven M. Anderson. "On adaptive extensions of group sequential trials for clinical investigations." \emph{Journal of the American Statistical Association} 103.484 (2008): 1621-1630.
+#'
+#' Maurer, Willi, and Frank Bretz. "Multiple testing in group sequential trials using graphical approaches." \emph{Statistics in Biopharmaceutical Research} 5.4 (2013): 311-320.
+#' @details
+#' Solution is found with a search using \code{uniroot}.
+#' This finds the maximum alpha-level for which an efficacy bound is crossed,
+#' completely ignoring any futility bound.
+#' @export
+#'
+#' @examples
+#' # Derive Group Sequential Design using gsDesign2
+#' x <- gs_design_ahr(
+#'   enroll_rate = define_enroll_rate(duration = c(2, 2, 2, 6),
+#'                                    rate = c(2.5, 5, 7.5, 10)),
+#'   fail_rate = define_fail_rate(duration = Inf, fail_rate = log(2) / 6,
+#'                                hr = 0.6, dropout_rate = .001),
+#'   info_frac = c(.5, .65, .8, 1),
+#'   analysis_time = 30,
+#'   upper = gs_spending_bound,
+#'   upar = list(sf = "sfLDOF", total_spend = 0.025),
+#'   lower = gs_spending_bound,
+#'   lpar = list(sf = "sfHSD", total_spend = 0.1, param = 2),
+#'   binding = FALSE
+#' )
+#'
+#' # Analysis at interim analysis 2
+#' sequential_pval(gs_design = x, event = c(100, 160), z = c(1.5, 2))
+#'
+#' # Use planned spending for interim and final analyses
+#' sequential_pval(gs_design = x,
+#'                 event = c(100, 160, 190),
+#'                 z = c(1.5, 2, 2.5))
+sequential_pval <- function(gs_design,
+                            event = NULL,
+                            z = NULL,
+                            ustime = NULL,
+                            interval = c(.00001, .9999)) {
+  q_e <- gs_design$input$ratio / (1 + gs_design$input$ratio)
+
+  # gs_design should be object returned from gs_design_ahr or gs_power_ahr with >=2 analyses
+  if (gs_design$design != "ahr") stop("gs_design must be an object created with gs_design_ahr() or gs_power_ahr()")
+  if (nrow(gs_design$analysis) == 1) stop("gs_design must be an GSD object with more than 1 analysis")
+
+  # check if gs_design is 1-sided, or has non-binding futility bounds
+  one_sided <- all(gs_design$bound$bound == "upper")
+  nonbinding <- !gs_design$input$binding
+  if (!(one_sided || nonbinding)) stop("gs_design must be one-sided or have non-binding futility bounds")
+
+  # if event is specified, check if it is an increasing numerical vector
+  if (!is.null(event)) check_increasing(event)
+
+  # if ustime is specified, check if it is an increasing numerical vector ranging within (0, 1]
+  if (!is.null(ustime)) check_info_frac(ustime)
+
+  # if event not specified, assume planned values used
+  if (is.null(event)) event <- gs_design$analysis$event
+
+  # if ustime not specified, use information fraction from design
+  if (is.null(ustime)) ustime <- event/max(gs_design$analysis$event)
+
+  # if z not specified, stop
+  if (is.null(z)) stop("z must be provided")
+
+  # check that z is numeric vector
+  if (!is.vector(z, mode = "numeric")) stop("z must be numeric vector")
+
+  # check lengths match
+  if (length(event) != length(z)) stop("event and z must have the same length")
+  if (length(event) != length(ustime)) stop("event and ustime must have the same length")
+
+  # check that ustime has at most 1 value >= 1
+  if (sum(ustime >= 1) > 1) stop("No more than one value of ustime >= 1")
+
+  # check the interval is 2 values in (0,1), exclusive of endpoints
+  if (!is.vector(interval, mode = "numeric") || length(interval) != 2 || min(interval) <= 0 || max(interval) >= 1) {
+    stop("interval must be 2 distinct values strictly between 0 and 1")
+  }
+
+  # Get spending function and parameters from gs_design
+  # sf is a character string, so we need to get the actual function from gsDesign namespace
+  sf_upper <- get_sf(gs_design$input$upar$sf)
+  sf_param <- gs_design$input$upar$param
+
+  # check upper end of p-value interval input
+  probhi <- sf_upper(alpha = max(interval), t = ustime, param = sf_param)$spend
+  if (length(z) > 1) probhi <- probhi - c(0, probhi[1:(length(z) - 1)])
+  if (min(gsDesign::gsBound1(I = event * (1 - q_e) * q_e, theta = 0, a = rep(-20, length(z)),
+                             probhi = probhi)$b - z) > 0) return(max(interval))
+
+  # check lower end of p-value interval input
+  probhi <- sf_upper(alpha = min(interval), t = ustime, param = sf_param)$spend
+  if (length(z) > 1) probhi <- probhi - c(0, probhi[1:(length(z) - 1)])
+  if (min(gsDesign::gsBound1(I = event * (1 - q_e) * q_e, theta = 0, a = rep(-20, length(z)),
+                             probhi = probhi)$b - z) < 0) return(min(interval))
+
+  # if answer is between interval bounds, find it with root-finding
+  x <- try(uniroot(sequential_zdiff, interval = -qnorm(interval), gs_design = gs_design,
+                   event = event, z = z, ustime = ustime))
+
+  if (inherits(x, "try-error")) stop("Failed to find root for sequential p-value")
+
+  pnorm(-x$root)
+}
+
+#' This function calculates the difference between:
+#' Computed efficacy bound at a given alpha level (using the spending function)
+#' Observed Z-statistic
+#' @noRd
+sequential_zdiff <- function(x,
+                             gs_design,
+                             event = (1:3) * 50,
+                             z = 1:3,
+                             ustime = NULL) {
+  q_e <- gs_design$input$ratio / (1 + gs_design$input$ratio)
+  alpha <- pnorm(-x)
+
+  sf_upper <- get_sf(gs_design$input$upar$sf)
+  sf_param <- gs_design$input$upar$param
+
+  probhi <- sf_upper(alpha = alpha, t = ustime, param = sf_param)$spend
+  if (length(z) > 1) probhi <- probhi - c(0, probhi[1:(length(z) - 1)])
+
+  return(min(gsDesign::gsBound1(I = event * (1 - q_e) * q_e, theta = 0, a = rep(-20, length(z)),
+                                probhi = probhi)$b - z))
+}

_pkgdown.yml

@@ -49,6 +49,7 @@ reference:
     - gs_design_ahr
     - gs_update_ahr
     - ahr_blinded
+    - sequential_pval
 - title: "Weighted logrank"
   desc: >
     Functions for the weighted logrank test (WLR) method.

man/sequential_pval.Rd

@@ -0,0 +1,56 @@
+test_that("Comparision with gsDesign::sequentialPValue", {
+  alpha <- 0.025
+  beta <- 0.1
+  ratio <- 1
+
+  # Enrollment
+  enroll_rate <- define_enroll_rate(
+    duration = c(2, 2, 10),
+    rate = (1:3) / 3)
+
+  # Failure and dropout
+  fail_rate <- define_fail_rate(
+    duration = Inf, fail_rate = log(2) / 9,
+    hr = 0.6, dropout_rate = .0001)
+
+  # IA and FA analysis time
+  study_duration <- 36
+
+  # Randomization ratio
+  ratio <- 1
+
+  # Spending
+  upper <- gs_spending_bound
+  lower <- gs_b
+  upar <- list(sf = "sfLDOF", total_spend = alpha)
+  lpar <- rep(-Inf, 3)
+
+  # ------------------------------ #
+  # original design of gsDesign    #
+  # ------------------------------ #
+  x_gsd <- gsSurv(k = 3, test.type = 1, alpha = alpha, beta = beta,
+                  astar = 0, timing = 1:3/3,
+                  sfu = sfLDOF, sfupar = 0,
+                  sfl = sfLDOF, sflpar = 0,
+                  lambdaC = log(2) / 9, hr = 0.6, hr0 = 1,
+                  eta = fail_rate$dropout_rate |> unique(),
+                  gamma = enroll_rate$rate,
+                  R = enroll_rate$duration,
+                  S = NULL, T = study_duration,
+                  minfup = study_duration - sum(enroll_rate$duration),
+                  ratio = ratio)
+
+  # ------------------------------ #
+  # original design of gsDesign2  #
+  # ------------------------------ #
+  x_gsd2 <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate,
+                          alpha = alpha, beta = beta, info_frac = 1:3/3, ratio = ratio,
+                          upper = upper, upar = upar,
+                          lower = lower, lpar = lpar,
+                          analysis_time = study_duration)
+
+  seq_pva_gsd <- sequentialPValue(x_gsd, n.I = c(40, 120), Z = c(1.5,2))
+  seq_pva_gsd2 <- sequential_pval(x_gsd2, event = c(40, 120) * (x_gsd2$analysis$event[3] / x_gsd$n.I[3]), z = c(1.5, 2))
+
+  expect_equal(seq_pva_gsd, seq_pva_gsd2)
+})