R/summary.R

In keaven/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/>.

#' Summary for fixed design objects
#'
#' Summary for [fixed_design()] objects
#'
#' @param object A fixed design object returned by [fixed_design()]
#' @param ... Additional arguments
#'
#' @return A data frame
#'
#' @rdname summary.fixed_design
#' @method summary fixed_design
#' @export
#'
#' @examples
#' library(dplyr)
#'
#' # Enrollment rate
#' enrollRates <- tibble::tibble(
#'   Stratum = "All",
#'   duration = 18,
#'   rate = 20)
#'
#' # Failure rates
#' failRates <- tibble::tibble(
#'   Stratum = "All",
#'   duration = c(4, 100),
#'   failRate = log(2) / 12,
#'   hr = c(1, .6),
#'   dropoutRate = .001)
#'
#' # Study duration in months
#' studyDuration <- 36
#'
#' # Experimental / Control randomization ratio
#' ratio <- 1
#'
#' # 1-sided Type I error
#' alpha <- 0.025
#' # Type II error (1 - power)
#' beta <- 0.1
#'
#' # ------------------------- #
#' #        AHR                #
#' # ------------------------- #
#' # under fixed power
#' fixed_design(
#'   x = "AHR",
#'   alpha = alpha,
#'   power = 1 - beta,
#'   enrollRates = enrollRates,
#'   failRates = failRates,
#'   studyDuration = studyDuration,
#'   ratio = ratio
#'   ) %>% summary()
#'
#' # ------------------------- #
#' #        FH                 #
#' # ------------------------- #
#' # under fixed power
#' fixed_design(
#'   x = "FH",
#'   alpha = alpha,
#'   power = 1 - beta,
#'   enrollRates = enrollRates,
#'   failRates = failRates,
#'   studyDuration = studyDuration,
#'    ratio = ratio
#'   ) %>% summary()
#'
summary.fixed_design <- function(object, ...){
  x <- object
  x_design <- switch(x$design,
                     "AHR" = {"Average hazard ratio"},
                     "LF" = {"Lachin and Foulkes"},
                     "RD" = {"Risk difference"},
                     "Milestone" = {paste0("Milestone: tau = ", x$design_par$tau)},
                     "RMST" = {paste0("RMST: tau = ", x$design_par$tau)},
                     "MB" = {paste0("Modestly weighted LR: tau = ", x$design_par$tau)},
                     "FH" = {
                       if(x$design_par$rho == 0 & x$design_par$gamma == 0){
                         paste0("Fleming-Harrington FH(0, 0) (logrank)")
                       }else{
                         paste0("Fleming-Harrington FH(", x$design_par$rho, ", ", x$design_par$gamma, ")")
                       }
                     },
                     "MaxCombo" = {
                       temp <- paste0("MaxCombo: FH(",
                                      paste(apply(do.call(rbind, x$design_par[c(1:2)]), 2 , paste , collapse = ", " ), collapse = "), FH("),
                                      ")")
                       gsub(pattern = "FH\\(0, 0\\)", replacement = "logrank", x = temp)
                     }
  )

  ans <- x$analysis %>% mutate(Design = x_design)
  class(ans) <- c("fixed_design", x$design, class(ans))
  return(ans)
}


#' Generate a table summarizing the bounds in the group sequential design
#'
#' Generate a table summarizing the bounds in the group sequential design
#' generated by [gs_design_ahr()], [gs_design_wlr()], or [gs_design_combo()].
#'
#' @param object An object returned by [gs_design_ahr()], [gs_design_wlr()], or [gs_design_combo()]
#' @param analysis_vars The variables to be put at the summary header of each analysis
#' @param analysis_decimals The displayed number of digits of `analysis_vars`
#' @param col_vars The variables to be displayed
#' @param col_decimals The decimals to be displayed for the displayed variables in `col_vars`
#' @param bound_names Names for bounds; default is `c("Efficacy", "Futility")`.
#' @param ... Additional arguments
#'
#' @return A summary table
#'
#' @rdname summary.gs_design
#' @method summary gs_design
#' @export
#'
#' @examples
#' # ---------------------------- #
#' #     design parameters        #
#' # ---------------------------- #
#' library(tibble)
#' library(gsDesign)
#' library(gsDesign2)
#' library(dplyr)
#'
#' # enrollment/failure rates
#' enrollRates <- tibble(Stratum = "All",
#'                       duration = 12,
#'                       rate = 1)
#' failRates <- tibble(Stratum = "All", duration = c(4, 100),
#'                     failRate = log(2) / 12,
#'                     hr = c(1, .6),
#'                      dropoutRate = .001)
#'
#' # Information fraction
#' IF <- (1:3)/3
#'
#' # Analysis times in months; first 2 will be ignored as IF will not be achieved
#' analysisTimes <- c(.01, .02, 36)
#'
#' # Experimental / Control randomization ratio
#' ratio <- 1
#'
#' # 1-sided Type I error
#' alpha <- 0.025
#'
#' # Type II error (1 - power)
#' beta <- .1
#'
#' # Upper bound
#' upper <- gs_spending_bound
#' upar <- list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
#'
#' # Lower bound
#' lower <- gs_spending_bound
#' lpar <- list(sf = gsDesign::sfHSD, total_spend = 0.1, param = 0, timing = NULL)
#'
#' # weight function in WLR
#' wgt00 <- function(x, arm0, arm1){wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = 0)}
#' wgt05 <- function(x, arm0, arm1){wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = .5)}
#'
#' # test in COMBO
#' fh_test <- rbind(
#'   data.frame(rho = 0, gamma = 0, tau = -1, test = 1, Analysis = 1:3,analysisTimes = c(12, 24, 36)),
#'   data.frame(rho = c(0, 0.5), gamma = 0.5, tau = -1, test = 2:3, Analysis = 3, analysisTimes = 36)
#' )
#'
#' # ---------------------------- #
#' #          ahr                 #
#' # ---------------------------- #
#' x_ahr <- gs_design_ahr(
#'   enrollRates = enrollRates,
#'   failRates = failRates,
#'   IF = IF, # Information fraction
#'   analysisTimes = analysisTimes,
#'   ratio = ratio,
#'   alpha = alpha,
#'   beta = beta,
#'   upper = upper,
#'   upar = upar,
#'   lower = lower,
#'   lpar = lpar)
#'
#' x_ahr %>% summary()
#' x_ahr %>% summary(analysis_vars = c("Time", "Events", "IF"), analysis_decimals = c(1, 0, 2))
#' x_ahr %>% summary(bound_names = c("A is better", "B is better"))
#'
#' # ---------------------------- #
#' #         wlr                  #
#' # ---------------------------- #
#' x_wlr <- gs_design_wlr(
#'   enrollRates = enrollRates,
#'   failRates = failRates,
#'   weight = wgt05,
#'   IF = NULL,
#'   analysisTimes = sort(unique(x_ahr$analysis$Time)),
#'   ratio = ratio,
#'   alpha = alpha,
#'   beta = beta,
#'   upper = upper,
#'   upar = upar,
#'   lower = lower,
#'   lpar = lpar
#' )
#' x_wlr %>% summary()
#'
#' # ---------------------------- #
#' #         max combo            #
#' # ---------------------------- #
#' x_combo <- gs_design_combo(
#'   ratio = 1,
#'   alpha = 0.025,
#'   beta = 0.2,
#'   enrollRates = tibble::tibble(Stratum = "All", duration = 12, rate = 500/12),
#'   failRates = tibble::tibble(Stratum = "All", duration = c(4, 100),
#'                              failRate = log(2) / 15, hr = c(1, .6), dropoutRate = .001),
#'   fh_test = fh_test,
#'   upper = gs_spending_combo,
#'   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
#'   lower = gs_spending_combo,
#'   lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.2))
#' x_combo %>% summary()
#'
#' # ---------------------------- #
#' #      risk difference         #
#' # ---------------------------- #
#' gs_design_rd(
#'   p_c = tibble(Stratum = "All", Rate = .2),
#'   p_e = tibble(Stratum = "All", Rate = .15),
#'   IF = c(0.7, 1),
#'   rd0 = 0,
#'   alpha = .025,
#'   beta = .1,
#'   ratio = 1,
#'   stratum_prev = NULL,
#'   weight = "un-stratified",
#'   upper = gs_b,
#'   lower = gs_b,
#'   upar = gsDesign::gsDesign(k = 3, test.type = 1, sfu = gsDesign::sfLDOF, sfupar = NULL)$upper$bound,
#'   lpar = c(qnorm(.1), rep(-Inf, 2))
#' ) %>% summary()
#'
summary.gs_design <- function(
  object,
  analysis_vars = NULL,
  analysis_decimals = NULL,
  col_vars = NULL,
  col_decimals = NULL,
  bound_names = c("Efficacy", "Futility"),
  ...
){
  x <- object
  method <- class(x)[class(x) %in% c("ahr", "wlr", "combo", "rd")]
  x_bounds <- x$bounds
  x_analysis <- x$analysis
  K <- max(x_analysis$Analysis)

  # --------------------------------------------- #
  #     prepare the columns decimals              #
  # --------------------------------------------- #
  if(method == "ahr"){
    if(is.null(col_vars) & is.null(col_decimals)){
      x_decimals <- tibble::tibble(
        col_vars = c("Analysis", "Bound", "Z", "~HR at bound", "Nominal p", "Alternate hypothesis", "Null hypothesis"),
        col_decimals = c(NA, NA, 2, 4, 4, 4, 4))
    }else{
      x_decimals <- tibble::tibble(col_vars = col_vars, col_decimals = col_decimals)
    }
  }
  if(method == "wlr"){
    if(is.null(col_vars) & is.null(col_decimals)){
      x_decimals <- tibble::tibble(
        col_vars = c("Analysis", "Bound", "Z", "~wHR at bound", "Nominal p", "Alternate hypothesis", "Null hypothesis"),
        col_decimals = c(NA, NA, 2, 4, 4, 4, 4))
    }else{
      x_decimals <- tibble::tibble(col_vars = col_vars, col_decimals = col_decimals)
    }
  }
  if(method == "combo"){
    if(is.null(col_vars) & is.null(col_decimals)){
      x_decimals <- tibble::tibble(
        col_vars = c("Analysis", "Bound", "Z", "Nominal p", "Alternate hypothesis", "Null hypothesis"),
        col_decimals = c(NA, NA, 2, 4, 4, 4))
    }else{
      x_decimals <- tibble::tibble(col_vars = col_vars, col_decimals = col_decimals)
    }
  }

  if(method == "rd"){
    if(is.null(col_vars) & is.null(col_decimals)){
      x_decimals <- tibble::tibble(
        col_vars = c("Analysis", "Bound", "Z", "~Risk difference at bound", "Nominal p", "Alternate hypothesis", "Null hypothesis"),
        col_decimals = c(NA, NA, 2, 4, 4, 4, 4))
    }else{
      x_decimals <- tibble::tibble(col_vars = col_vars, col_decimals = col_decimals)
    }
  }

  # --------------------------------------------- #
  #     prepare the analysis summary row          #
  # --------------------------------------------- #
  # get the
  # (1) analysis variables to be displayed on the header
  # (2) decimals to be displayed for the analysis variables in (3)
  if(is.null(analysis_vars) & is.null(analysis_decimals)){
    if(method %in% c("ahr", "wlr")){
      analysis_vars <- c("Time", "N", "Events", "AHR", "IF")
      analysis_decimals <- c(1, 1, 1, 2, 2)
    }
    if(method == "combo"){
      analysis_vars <- c("Time", "N", "Events", "AHR", "EF")
      analysis_decimals <- c(1, 1, 1, 2, 2)
    }
    if(method == "rd"){
      analysis_vars <- c("N", "rd", "IF")
      analysis_decimals <- c(1, 4, 2)
    }
  }else if(is.null(analysis_vars) & !is.null(analysis_decimals)){
    stop("summary: please input analysis_vars and analysis_decimals in pairs!")
  }else if(!is.null(analysis_vars) & is.null(analysis_decimals)){
    stop("summary: please input analysis_vars and analysis_decimals in pairs!")
  }
  # set the analysis summary header
  analyses <- x_analysis %>%
    dplyr::group_by(Analysis) %>%
    dplyr::filter(dplyr::row_number() == 1) %>%
    dplyr::select(all_of(c("Analysis", analysis_vars))) %>%
    dplyr::arrange(Analysis)

  # --------------------------------------------- #
  #             merge 2 tables:                   #
  #         (1) alternate hypothesis table        #
  #         (2) null hypothesis table             #
  # --------------------------------------------- #
  # table A: a table under alternative hypothesis
  xy <- x_bounds %>%
    dplyr::rename("Alternate hypothesis" = Probability) %>%
    dplyr::rename("Null hypothesis" = Probability0) %>%
    # change Upper -> bound_names[1], e.g., Efficacy
    # change Lower -> bound_names[2], e.g., Futility
    dplyr::mutate(Bound = dplyr::recode(Bound, "Upper" = bound_names[1], "Lower" = bound_names[2]))

  if("Probability0" %in% colnames(x_bounds)){
    xy <- x_bounds %>%
      dplyr::rename("Alternate hypothesis" = Probability) %>%
      dplyr::rename("Null hypothesis" = Probability0)
  }else{
    xy <- x_bounds %>%
      dplyr::rename("Alternate hypothesis" = Probability) %>%
      tibble::add_column("Null hypothesis" = "-")
  }
  # change Upper -> bound_names[1], e.g., Efficacy
  # change Lower -> bound_names[2], e.g., Futility
  xy <- xy %>%
    dplyr::mutate(Bound = dplyr::recode(Bound, "Upper" = bound_names[1], "Lower" = bound_names[2]))  %>%
    dplyr::arrange(Analysis,desc(Bound))

  # tbl_a <- x_bounds %>%
  #   dplyr::filter(hypothesis == "H1") %>%
  #   dplyr::rename("Alternate hypothesis" = Probability) %>%
  #   # change Upper -> bound_names[1], e.g., Efficacy
  #   # change Lower -> bound_names[2], e.g., Futility
  #   dplyr::mutate(Bound = dplyr::recode(Bound, "Upper" = bound_names[1], "Lower" = bound_names[2]))
  #
  # # table B: a table under null hypothesis
  # tbl_b <- x_bounds %>%
  #   dplyr::filter(hypothesis == "H0") %>%
  #   dplyr::rename("Null hypothesis" = Probability) %>%
  #   dplyr::mutate(Bound = dplyr::recode(Bound, "Upper" = bound_names[1], "Lower" = bound_names[2])) %>%
  #   dplyr::select(all_of(c("Analysis", "Bound", "Null hypothesis")))
  #
  # xy <- full_join(tbl_a, tbl_b, by = c("Analysis", "Bound"))

  # --------------------------------------------- #
  #             merge 2 tables:                   #
  #         (1) analysis summary table            #
  #         (2) xy: bound_summary_detail table    #
  # --------------------------------------------- #
  # Merge 3 tables: 1 line per analysis, alternate hypothesis table, null hypothesis table
  # if the method is AHR
  if(method == "ahr"){
    # header
    analysis_summary_header <- analyses %>% dplyr::select(all_of(c("Analysis", analysis_vars)))
    # bound details
    bound_summary_detail <- xy
  }

  # if the method is WLR, change AHR to wAHR
  if(method == "wlr"){
    # header
    analysis_summary_header <- analyses %>% dplyr::select(all_of(c("Analysis", analysis_vars)))
    if("AHR" %in% analysis_vars){
      analysis_summary_header <- analysis_summary_header %>% dplyr::rename(wAHR = AHR)
    }
    # bound details
    if("~HR at bound" %in% names(xy)){
      bound_summary_detail <- xy %>% dplyr::rename("~wHR at bound" = "~HR at bound")
    }else{
      bound_summary_detail <- xy
    }
  }

  # if the method is COMBO, remove the column of "~HR at bound", and remove AHR from header
  if(method == "combo"){
    # header
    analysis_summary_header <- analyses %>% dplyr::select(all_of(c("Analysis", analysis_vars)))
    # bound details
    if("~HR at bound" %in% names(xy)){
      stop("summary: ~HR at bound can't be display!")
    }else{
      bound_summary_detail <- xy
    }
  }

  # if the method is RD
  if(method == "rd"){
    # header
    analysis_summary_header <- analyses %>%
      dplyr::select(all_of(c("Analysis", analysis_vars))) %>%
      dplyr::rename("risk difference" = rd)
    # bound details
    bound_summary_detail <- xy
  }

  output <- table_ab(
    # A data frame to be show as the summary header
    # It has only ONE record for each value of `byvar`
    table_a = analysis_summary_header,
    # A data frame to be shown as the listing details
    # It has >= 1 records for each value of `byvar`
    table_b = bound_summary_detail,
    decimals = c(0, analysis_decimals),
    byvar = "Analysis"
  ) %>%
    dplyr::group_by(Analysis)


  if(method == "ahr"){
    output <- output %>% select(Analysis, Bound, Z, `~HR at bound`, `Nominal p`, `Alternate hypothesis`, `Null hypothesis`)
  }else if(method == "wlr"){
    output <- output %>% select(Analysis, Bound, Z, `~wHR at bound`, `Nominal p`, `Alternate hypothesis`, `Null hypothesis`)
  }else if(method == "combo"){
    output <- output %>% select(Analysis, Bound, Z, `Nominal p`, `Alternate hypothesis`, `Null hypothesis`)
  }else if(method == "rd"){
    output <- output %>% select(Analysis, Bound, Z, `~Risk difference at bound`, `Nominal p`, `Alternate hypothesis`, `Null hypothesis`)
  }

  # --------------------------------------------- #
  #     set the decimals to display               #
  # --------------------------------------------- #
  output <- output %>% select(x_decimals$col_vars)
  if("Z" %in% colnames(output)){
    output <- output %>% dplyr::mutate_at("Z", round, (x_decimals %>% filter(col_vars == "Z"))$col_decimals)
  }
  if("~HR at bound" %in% colnames(output)){
    output <- output %>% dplyr::mutate_at("~HR at bound", round, (x_decimals %>% filter(col_vars == "~HR at bound"))$col_decimals)
  }
  if("~Risk difference at bound" %in% colnames(output)){
    output <- output %>% dplyr::mutate_at("~Risk difference at bound", round, (x_decimals %>% filter(col_vars == "~Risk difference at bound"))$col_decimals)
  }
  if("Nominal p" %in% colnames(output)){
    output <- output %>% dplyr::mutate_at("Nominal p", round, (x_decimals %>% filter(col_vars == "Nominal p"))$col_decimals)
  }
  if("Alternate hypothesis" %in% colnames(output)){
    output <- output %>% dplyr::mutate_at("Alternate hypothesis", round, (x_decimals %>% filter(col_vars == "Alternate hypothesis"))$col_decimals)
  }
  if("Null hypothesis" %in% colnames(output) & is.vector(output[["Null hypothesis"]], mode = "numeric")){
    output <- output %>% dplyr::mutate_at("Null hypothesis", round, (x_decimals %>% filter(col_vars == "Null hypothesis"))$col_decimals)
  }

  class(output) <- c(method, "gs_design", class(output))
  return(output)
}