R/gs_power_npe.R

In gsDesign2: Group Sequential Design with Non-Constant Effect

#  Copyright (c) 2025 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 computation with non-constant effect and information.
#'
#' @description The following two functions allow a non-constant treatment effect over time,
#' but also can be applied for the usual homogeneous effect size designs.
#' They require treatment effect and statistical information at each analysis
#' as well as a method of deriving bounds, such as spending.
#' Initial bound types supported are spending bounds and fixed bounds.
#' These routines enables two things not available in the gsDesign package: 1)
#' non-constant effect, 2) more flexibility in boundary selection.
#'
#' \code{gs_power_npe()} derives group sequential bounds and boundary crossing probabilities
#' for a design, given treatment effect and information at each analysis and the
#' method of deriving bounds, such as spending.
#'
#' \code{gs_design_npe()} derives group sequential design size,
#' bounds and boundary crossing probabilities based on proportionate
#' information and effect size at analyses, as well as the
#' method of deriving bounds, such as spending.
#'
#' The only differences in arguments between the two functions are the \code{alpha} and \code{beta}
#' parameters used in the \code{gs_design_npe()}.
#'
#' @param theta Natural parameter for group sequential design representing
#'   expected cumulative drift at all analyses; used for power calculation.
#'   It can be a scalar (constant treatment effect) or a vector (non-constant treatment effect).
#'   The user must provide a value for `theta`.
#' @param theta0 Natural parameter for null hypothesis.
#'   It can be a scalar (constant treatment effect) or a vector (non-constant treatment effect).
#'   The default is 0. If a value other than 0 is provided, it affects upper bound computation.
#' @param theta1 Natural parameter for alternate hypothesis,
#'   if needed for lower bound computation.
#'   It can be a scalar (constant treatment effect) or a vector (non-constant treatment effect).
#'   The default is the same as `theta`, which yields the usual beta-spending.
#'   If set to 0, spending is 2-sided under the null hypothesis.
#' @param info Statistical information at all analyses for input `theta`.
#'   It is a vector of positive numbers in increasing order.
#'   The user must provide values corresponding to `theta`.
#' @param info0 Statistical information under null hypothesis.
#'   It is a vector of all positive numbers with increasing order.
#'   Default is set to be the same as `info`.
#'   If `info0` is different than `info`, it impacts null hypothesis bound calculation.
#' @param info1 Statistical information under hypothesis used for futility bound calculation.
#'   It is a vector of all positive numbers with increasing order.
#'   Default is set to be the same as `info`.
#'   If `info1` is different from `info`, it impacts futility bound calculation.
#' @param info_scale Information scale for calculation. Options are:
#'   - `"h0_h1_info"` (default): variance under both null and alternative hypotheses is used.
#'   - `"h0_info"`: variance under null hypothesis is used.
#'   This is often used for testing methods that use local alternatives, such as the Schoenfeld method.
#'   - `"h1_info"`: variance under alternative hypothesis is used.
#' @param binding Indicator of whether futility bound is binding;
#'   default of `FALSE` is recommended.
#' @param upper Function to compute upper bound.
#'   - \code{gs_spending_bound()}: alpha-spending efficacy bounds.
#'   - \code{gs_b()}: fixed efficacy bounds.
#' @param lower Function to compute lower bound, which can be set up similarly as `upper`.
#' See [this vignette](https://merck.github.io/gsDesign2/articles/story-seven-test-types.html).
#' @param upar Parameters passed to `upper`.
#'   - If `upper = gs_b`, then `upar` is a numerical vector specifying the fixed efficacy bounds per analysis.
#'   - If `upper = gs_spending_bound`, then `upar` is a list including
#'       - `sf` for the spending function family.
#'       - `total_spend` for total alpha spend.
#'       - `param` for the parameter of the spending function.
#'       - `timing` specifies spending time if different from information-based spending; see details.
#' @param lpar Parameters passed to `lower`, which can be set up similarly as `upar.`
#' @param test_upper Indicator of which analyses should include
#'   an upper (efficacy) bound;
#'   single value of `TRUE` (default) indicates all analyses; otherwise,
#'   a logical vector of the same length as `info` should
#'   indicate which analyses will have an efficacy bound.
#' @param test_lower Indicator of which analyses should include a lower bound;
#'   single value of `TRUE` (default) indicates all analyses;
#'   single value of `FALSE` indicated no lower bound; otherwise,
#'   a logical vector of the same length as `info` should
#'   indicate which analyses will have a lower bound.
#' @param r Integer value controlling grid for numerical integration as in
#'   Jennison and Turnbull (2000); default is 18, range is 1 to 80.
#'   Larger values provide larger number of grid points and greater accuracy.
#'   Normally, `r` will not be changed by the user.
#' @param tol Tolerance parameter for boundary convergence (on Z-scale); normally not changed by the user.
#'
#' @return A tibble with columns of
#'   - `analysis`: analysis index.
#'   - `bound`: either of value `"upper"` or `"lower"`, indicating the upper and lower bound.
#'   - `z`: the Z-score bounds.
#'   - `probability`: cumulative probability of crossing the bound at or before the analysis.
#'   - `theta`: same as the input.
#'   - `theta1`: same as the input.
#'   - `info`: statistical information at each analysis.
#'      + If it is returned by `gs_power_npe`, the `info`, `info0`, `info1` are same as the input.
#'      + If it is returned by `gs_design_npe`, the `info`, `info0`, `info1` are changed by a constant scale factor.
#' factor to ensure the design has power `1 - beta`.
#'   - `info0`: statistical information under the null at each analysis.
#'   - `info1`: statistical information under the alternative at each analysis.
#'   - `info_frac`: information fraction at each analysis, i.e., `info / max(info)`.
#'
#' @section Specification:
#' \if{latex}{
#'  \itemize{
#'    \item Extract the length of input info as the number of interim analysis.
#'    \item Validate if input info0 is NULL, so set it equal to info.
#'    \item Validate if the length of inputs info and info0 are the same.
#'    \item Validate if input theta is a scalar, so replicate
#'    the value for all k interim analysis.
#'    \item Validate if input theta1 is NULL and if it is a scalar.
#'    If it is NULL, set it equal to input theta. If it is a scalar,
#'    replicate the value for all k interim analysis.
#'    \item Validate if input test_upper is a scalar,
#'    so replicate the value for all k interim analysis.
#'    \item Validate if input test_lower is a scalar,
#'    so replicate the value for all k interim analysis.
#'    \item Define vector a to be -Inf with
#'    length equal to the number of interim analysis.
#'    \item Define vector b to be Inf with
#'    length equal to the number of interim analysis.
#'    \item Define hgm1_0 and hgm1 to be NULL.
#'    \item Define upper_prob and lower_prob to be
#'    vectors of NA with length of the number of interim analysis.
#'    \item Update lower and upper bounds using \code{gs_b()}.
#'    \item If there are no interim analysis, compute probabilities
#'    of crossing upper and lower bounds
#'    using \code{h1()}.
#'    \item Compute cross upper and lower bound probabilities
#'    using \code{hupdate()} and \code{h1()}.
#'    \item Return a tibble of analysis number, bound, z-values,
#'    probability of crossing bounds,
#'    theta, theta1, info, and info0.
#'   }
#' }
#' \if{html}{The contents of this section are shown in PDF user manual only.}
#'
#' @name gs_power_design_npe
#' @rdname gs_power_design_npe
#' @export
#'
#' @examples
#'
#' # Example 6 ----
#' # Default of gs_power_npe (single analysis; Type I error controlled)
#' gs_power_npe(theta = 0) |> dplyr::filter(bound == "upper")
#'
#' # Example 7 ----
#' # gs_power_npe with fixed bound
#' gs_power_npe(
#'   theta = c(.1, .2, .3),
#'   info = (1:3) * 40,
#'   upper = gs_b,
#'   upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
#'   lower = gs_b,
#'   lpar = c(-1, 0, 0)
#' )
#'
#' # Same fixed efficacy bounds, no futility bound (i.e., non-binding bound), null hypothesis
#' gs_power_npe(
#'   theta = rep(0, 3),
#'   info = (1:3) * 40,
#'   upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
#'   lpar = rep(-Inf, 3)
#' ) |>
#'   dplyr::filter(bound == "upper")
#'
#' # Example 8 ----
#' # gs_power_npe with fixed bound testing futility only at analysis 1; efficacy only at analyses 2, 3
#' gs_power_npe(
#'   theta = c(.1, .2, .3),
#'   info = (1:3) * 40,
#'   upper = gs_b,
#'   upar = c(Inf, 3, 2),
#'   lower = gs_b,
#'   lpar = c(qnorm(.1), -Inf, -Inf)
#' )
#'
#' # Example 9 ----
#' # gs_power_npe with spending function bounds
#' # Lower spending based on non-zero effect
#' gs_power_npe(
#'   theta = c(.1, .2, .3),
#'   info = (1:3) * 40,
#'   upper = gs_spending_bound,
#'   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
#'   lower = gs_spending_bound,
#'   lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -1, timing = NULL)
#' )
#'
#' # Same bounds, but power under different theta
#' gs_power_npe(
#'   theta = c(.15, .25, .35),
#'   info = (1:3) * 40,
#'   upper = gs_spending_bound,
#'   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
#'   lower = gs_spending_bound,
#'   lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -1, timing = NULL)
#' )
#'
#' # Example 10 ----
#' # gs_power_npe with two-sided symmetric spend, O'Brien-Fleming spending
#' # Typically, 2-sided bounds are binding
#' x <- gs_power_npe(
#'   theta = rep(0, 3),
#'   info = (1:3) * 40,
#'   binding = TRUE,
#'   upper = gs_spending_bound,
#'   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
#'   lower = gs_spending_bound,
#'   lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
#' )
#'
#' # Re-use these bounds under alternate hypothesis
#' # Always use binding = TRUE for power calculations
#' gs_power_npe(
#'   theta = c(.1, .2, .3),
#'   info = (1:3) * 40,
#'   binding = TRUE,
#'   upar = (x |> dplyr::filter(bound == "upper"))$z,
#'   lpar = -(x |> dplyr::filter(bound == "upper"))$z
#' )
#'
#' # Example 11 ----
#' # Different values of `r` and `tol` lead to different numerical accuracy
#' # Larger `r` and smaller `tol` give better accuracy, but leads to slow computation
#' n_analysis <- 5
#' gs_power_npe(
#'   theta = 0.1,
#'   info = 1:n_analysis,
#'   info0 = 1:n_analysis,
#'   info1 = NULL,
#'   info_scale = "h0_info",
#'   upper = gs_spending_bound,
#'   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
#'   lower = gs_b,
#'   lpar = -rep(Inf, n_analysis),
#'   test_upper = TRUE,
#'   test_lower = FALSE,
#'   binding = FALSE,
#'   # Try different combinations of (r, tol) with
#'   # r in 6, 18, 24, 30, 35, 40, 50, 60, 70, 80, 90, 100
#'   # tol in 1e-6, 1e-12
#'   r = 6,
#'   tol = 1e-6
#' )
gs_power_npe <- function(theta = .1, theta0 = 0, theta1 = theta, # 3 theta
                         info = 1, info0 = NULL, info1 = NULL, # 3 info
                         info_scale = c("h0_h1_info", "h0_info", "h1_info"),
                         upper = gs_b, upar = qnorm(.975),
                         lower = gs_b, lpar = -Inf,
                         test_upper = TRUE, test_lower = TRUE, binding = FALSE,
                         r = 18, tol = 1e-6) {
  # Check & set up parameters ----
  n_analysis <- length(info)
  theta  <- check_theta(theta,  n_analysis)
  theta0 <- check_theta(theta0, n_analysis)
  theta1 <- check_theta(theta1, n_analysis)
  upper <- match.fun(upper)
  lower <- match.fun(lower)
  test_upper <- check_test_upper(test_upper, n_analysis)
  test_lower <- check_test_lower(test_lower, n_analysis)

  # Set up info ----
  # impute info
  if (is.null(info0)) {
    info0 <- info
  }

  if (is.null(info1)) {
    info1 <- info
  }

  # set up info_scale
  info_scale <- match.arg(info_scale)

  if (info_scale == "h0_info") {
    info <- info0
    info1 <- info0
  }
  if (info_scale == "h1_info") {
    info <- info1
    info0 <- info1
  }

  # check info
  check_info(info)
  check_info(info0)
  check_info(info1)
  if (length(info0) != length(info)) stop("gs_design_npe(): length of info, info0 must be the same!")
  if (length(info1) != length(info)) stop("gs_design_npe(): length of info, info1 must be the same!")

  # Initialization ----
  a <- rep(-Inf, n_analysis)
  b <- rep(Inf, n_analysis)
  hgm1_0 <- NULL
  hgm1_1 <- NULL
  hgm1 <- NULL
  upper_prob <- rep(NA, n_analysis)
  lower_prob <- rep(NA, n_analysis)

  # Calculate crossing prob under H1 ----
  for (k in 1:n_analysis) {
    # compute/update lower/upper bound
    a[k] <- lower(
      k = k, par = lpar, hgm1 = hgm1_1, info = info1, r = r, tol = tol, test_bound = test_lower,
      theta = theta1, efficacy = FALSE
    )
    b[k] <- upper(k = k, par = upar, hgm1 = hgm1_0, info = info0, r = r, tol = tol, test_bound = test_upper)

    # if it is the first analysis
    if (k == 1) {
      # compute the probability to cross upper/lower bound
      upper_prob[1] <- if (b[1] < Inf) {
        pnorm(sqrt(info[1]) * (theta[1] - b[1] / sqrt(info0[1])))
      } else {
        0
      }
      lower_prob[1] <- if (a[1] > -Inf) {
        pnorm(-sqrt(info[1]) * (theta[1] - a[1] / sqrt(info0[1])))
      } else {
        0
      }
      # update the grids
      hgm1_0 <- h1(r = r, theta = theta0[1], info = info0[1], a = if (binding) {
        a[1]
      } else {
        -Inf
      }, b = b[1])
      hgm1_1 <- h1(r = r, theta = theta1[1], info = info1[1], a = a[1], b = b[1])
      hgm1 <- h1(r = r, theta = theta[1], info = info[1], a = a[1], b = b[1])
    } else {
      # compute the probability to cross upper bound
      upper_prob[k] <- if (b[k] < Inf) {
        sum(hupdate(
          theta = theta[k], thetam1 = theta[k - 1],
          info = info[k], im1 = info[k - 1],
          a = b[k], b = Inf, gm1 = hgm1, r = r
        )$h)
      } else {
        0
      }
      # compute the probability to cross lower bound
      lower_prob[k] <- if (a[k] > -Inf) {
        sum(hupdate(
          theta = theta[k], thetam1 = theta[k - 1],
          info = info[k], im1 = info[k - 1],
          a = -Inf, b = a[k], gm1 = hgm1, r = r
        )$h)
      } else {
        0
      }

      # update the grids
      if (k < n_analysis) {
        hgm1_0 <- hupdate(r = r, theta = theta0[k], info = info0[k], a = if (binding) {
          a[k]
        } else {
          -Inf
        }, b = b[k], thetam1 = 0, im1 = info0[k - 1], gm1 = hgm1_0)
        hgm1_1 <- hupdate(
          r = r, theta = theta1[k], info = info1[k],
          a = a[k], b = b[k], thetam1 = theta1[k - 1],
          im1 = info1[k - 1], gm1 = hgm1_1
        )
        hgm1 <- hupdate(
          r = r, theta = theta[k], info = info[k],
          a = a[k], b = b[k], thetam1 = theta[k - 1],
          im1 = info[k - 1], gm1 = hgm1
        )
      }
    }
  }

  ans <- tibble(
    analysis = rep(1:n_analysis, 2),
    bound = c(rep("upper", n_analysis), rep("lower", n_analysis)),
    z = c(b, a),
    probability = c(cumsum(upper_prob), cumsum(lower_prob)),
    theta = rep(theta, 2),
    theta1 = rep(theta1, 2),
    info_frac = rep(info / max(info), 2),
    info = rep(info, 2)
  ) |>
    mutate(
      info0 = rep(info0, 2),
      info1 = rep(info1, 2)
    ) |>
    # filter(abs(Z) < Inf) |>
    arrange(desc(bound), analysis)

  return(ans)
}