tests/testthat/helper-old-version-gs_design_ahr_.R

In gsDesign2: Group Sequential Design with Non-Constant Effect

#  Copyright (c) 2022 Merck Sharp & Dohme Corp. a subsidiary of Merck & Co., Inc., Rahway, NJ, USA.
#
#  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/>.

#' @importFrom tibble tibble
#' @importFrom gsDesign gsDesign sfLDOF
#' @importFrom stats qnorm
#' @importFrom dplyr mutate full_join select arrange desc
NULL
#' Group sequential design using average hazard ratio under non-proportional hazards
#'
#' @param enrollRates enrollment rates
#' @param failRates failure and dropout rates
#' @param ratio Experimental:Control randomization ratio (not yet implemented)
#' @param alpha One-sided Type I error
#' @param beta Type II error
#' @param IF Targeted information fraction at each analysis
#' @param analysisTimes Minimum time of analysis
#' @param binding indicator of whether futility bound is binding; default of FALSE is recommended
#' @param upper Function to compute upper bound
#' @param upar Parameter passed to \code{upper()}
#' @param lower Function to compute lower bound
#' @param lpar Parameter passed to \code{lower()}
#' @param h1_spending Indicator that lower bound to be set by spending under alternate hypothesis (input \code{failRates})
#' if spending is used for lower bound
#' @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 \code{info} should indicate which analyses will have an efficacy bound
#' @param test_lower indicator of which analyses should include an lower bound; single value of TRUE (default) indicates all analyses;
#' single value FALSE indicated no lower bound; otherwise, a logical vector of the same length as \code{info} should indicate which analyses will have a
#' lower bound
#' @param r  Integer, at least 2; default of 18 recommended by Jennison and Turnbull
#' @param tol Tolerance parameter for boundary convergence (on Z-scale)
#' @section Specification:
#' \if{latex}{
#'  \itemize{
#'    \item Validate if input analysisTimes is a positive number or positive increasing sequence.
#'    \item Validate if input IF is a positive number or positive increasing sequence
#'    on (0, 1] with final value of 1.
#'    \item Validate if input IF and analysisTimes  have the same length if both have length > 1.
#'    \item Get information at input analysisTimes
#'    \itemize{
#'      \item Use \code{gs_info_ahr()} to get the information and effect size based on AHR approximation.
#'      \item Extract the final event.
#'      \item Check if input If needed for (any) interim analysis timing.
#'    }
#'    \item Add the analysis column to the information at input analysisTimes.
#'    \item Add the sample size column to the information at input analysisTimes using \code{eAccrual()}.
#'    \item Get sample size and bounds using \code{gs_design_npe()} and save them to bounds.
#'    \item Add Time, Events, AHR, N that have already been calculated to the bounds.
#'    \item Return a list of design enrollment, failure rates, and bounds.
#'   }
#' }
#' \if{html}{The contents of this section are shown in PDF user manual only.}
#'
#' @return a \code{tibble} with columns Analysis, Bound, Z, Probability, theta, Time, AHR, Events
#' @details Need to be added
#' @noRd
#'
#' @examples
#' library(gsDesign)
#' library(gsDesign2)
#' library(dplyr)
#' # call with defaults
#' gs_design_ahr()
#'
#' # Single analysis
#' gs_design_ahr(analysisTimes = 40)
#'
#' # Multiple analysisTimes
#' gs_design_ahr(analysisTimes = c(12, 24, 36))
#'
#' # Specified information fraction
#' gs_design_ahr(IF = c(.25, .75, 1), analysisTimes = 36)
#'
#' # multiple analysis times & IF
#' # driven by times
#' gs_design_ahr(IF = c(.25, .75, 1), analysisTimes = c(12, 25, 36))
#' # driven by IF
#' gs_design_ahr(IF = c(1 / 3, .8, 1), analysisTimes = c(12, 25, 36))
#'
#' # 2-sided symmetric design with O'Brien-Fleming spending
#' gs_design_ahr(
#'   analysisTimes = c(12, 24, 36),
#'   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
#'   ),
#'   h1_spending = FALSE
#' )
#'
#' # 2-sided asymmetric design with O'Brien-Fleming upper spending
#' # Pocock lower spending under H1 (NPH)
#' gs_design_ahr(
#'   analysisTimes = c(12, 24, 36),
#'   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::sfLDPocock, total_spend = 0.1,
#'     param = NULL, timing = NULL
#'   ),
#'   h1_spending = TRUE
#' )
gs_design_ahr_ <- function(enrollRates = tibble::tibble(
                             Stratum = "All",
                             duration = c(2, 2, 10),
                             rate = c(3, 6, 9)
                           ),
                           failRates = tibble::tibble(
                             Stratum = "All",
                             duration = c(3, 100),
                             failRate = log(2) / c(9, 18),
                             hr = c(.9, .6),
                             dropoutRate = rep(.001, 2)
                           ),
                           ratio = 1, # Experimental:Control randomization ratio
                           alpha = 0.025, # One-sided Type I error
                           beta = 0.1, # NULL if enrollment is not adapted
                           IF = NULL, # relative information fraction timing (vector, if not NULL; increasing to 1)
                           analysisTimes = 36, # Targeted times of analysis or just planned study duration
                           binding = FALSE,
                           upper = gs_spending_bound,
                           # Default is LDOF
                           upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
                           lower = gs_spending_bound,
                           lpar = list(sf = gsDesign::sfLDOF, total_spend = beta),
                           h1_spending = TRUE,
                           test_upper = TRUE,
                           test_lower = TRUE,
                           r = 18,
                           tol = 1e-6) {
  ################################################################################
  # Check input values
  msg <- "analysisTimes must be a positive number or positive increasing sequence"
  if (!is.vector(analysisTimes, mode = "numeric")) stop(msg)
  if (min(analysisTimes - dplyr::lag(analysisTimes, def = 0)) <= 0) stop(msg)
  msg <- "gs_design_ahr(): IF must be a positive number or positive increasing sequence on (0, 1] with final value of 1"
  if (is.null(IF)) {
    IF <- 1
  }
  if (!is.vector(IF, mode = "numeric")) stop(msg)
  if (min(IF - dplyr::lag(IF, def = 0)) <= 0) stop(msg)
  if (max(IF) != 1) stop(msg)
  msg <- "gs_design_ahr() IF and analysisTimes must have the same length if both have length > 1"
  if ((length(analysisTimes) > 1) & (length(IF) > 1) & (length(IF) != length(analysisTimes))) stop(msg)
  # end check input values
  ################################################################################
  # Get information at input analysisTimes
  y <- gs_info_ahr_(enrollRates, failRates, ratio = ratio, events = NULL, analysisTimes = analysisTimes)
  finalEvents <- y$Events[nrow(y)]
  IFalt <- y$Events / finalEvents
  # Check if IF needed for (any) IA timing
  K <- max(length(analysisTimes), length(IF))
  nextTime <- max(analysisTimes)
  if (length(IF) == 1) {
    IF <- IFalt
  } else {
    IFindx <- IF[1:(K - 1)]
    for (i in seq_along(IFindx)) {
      if (length(IFalt) == 1) {
        y <-
          rbind(
            tEvents_(enrollRates, failRates,
              targetEvents = IF[K - i] * finalEvents, ratio = ratio,
              interval = c(.01, nextTime)
            ) %>% dplyr::mutate(theta = -log(AHR), Analysis = K - i),
            y
          )
      } else if (IF[K - i] > IFalt[K - i]) {
        y[K - i, ] <-
          tEvents_(enrollRates, failRates,
            targetEvents = IF[K - i] * finalEvents, ratio = ratio,
            interval = c(.01, nextTime)
          ) %>%
          dplyr::transmute(Analysis = K - i, Time, Events, AHR, theta = -log(AHR), info, info0)
      }
      nextTime <- y$Time[K - i]
    }
  }
  y$Analysis <- 1:K
  y$N <- eAccrual_(x = y$Time, enrollRates = enrollRates)
  if (h1_spending) {
    theta1 <- y$theta
    info1 <- y$info
  } else {
    theta1 <- 0
    info1 <- y$info0
  }

  # Get sample size and bounds using gs_design_npe
  bounds <- gs_design_npe_(
    theta = y$theta,
    theta1 = theta1,
    info = y$info,
    info0 = y$info0,
    info1 = info1,
    alpha = alpha,
    beta = beta,
    binding = binding,
    upper = upper,
    lower = lower,
    upar = upar,
    lpar = lpar,
    test_upper = test_upper,
    test_lower = test_lower,
    r = r,
    tol = tol
  ) %>%
    # Add Time, Events, AHR, N from gs_info_ahr call above
    dplyr::full_join(y %>% dplyr::select(-c(info, info0, theta)), by = "Analysis") %>%
    dplyr::select(c("Analysis", "Bound", "Time", "N", "Events", "Z", "Probability", "AHR", "theta", "info", "info0")) %>%
    dplyr::arrange(dplyr::desc(Bound), Analysis)
  bounds$Events <- bounds$Events * bounds$info[K] / y$info[K]
  bounds$N <- bounds$N * bounds$info[K] / y$info[K]

  # Document design enrollment, failure rates, and bounds
  return(list(
    enrollRates = enrollRates %>%
      dplyr::mutate(rate = rate * bounds$info[K] / y$info[K]),
    failRates = failRates,
    bounds = bounds
  ))
}