R/summary.R

In Merck/simtrial: Clinical Trial Simulation

#  Copyright (c) 2024 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
#  All rights reserved.
#
#  This file is part of the simtrial program.
#
#  simtrial 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 of group sequential simulations.
#'
#' @param object Simulation results generated by [sim_gs_n()]
#' @param design Asymptotic design generated by [gsDesign2::gs_design_ahr()],
#' [gsDesign2::gs_power_ahr()], [gsDesign2::gs_design_wlr()], or [gsDesign2::gs_power_wlr].
#' @param bound The boundaries.
#' @param ... Additional parameters (not used).
#'
#' @rdname summary
#' @return A data frame
#'
#' @importFrom data.table ":=" .N .SD as.data.table dcast merge.data.table
#'
#' @export
#'
#' @examples
#' library(gsDesign2)
#'
#' # Parameters for enrollment
#' enroll_rampup_duration <- 4 # Duration for enrollment ramp up
#' enroll_duration <- 16 # Total enrollment duration
#' enroll_rate <- define_enroll_rate(
#'   duration = c(
#'     enroll_rampup_duration, enroll_duration - enroll_rampup_duration),
#'  rate = c(10, 30))
#'
#' # Parameters for treatment effect
#' delay_effect_duration <- 3 # Delay treatment effect in months
#' median_ctrl <- 9 # Survival median of the control arm
#' median_exp <- c(9, 14) # Survival median of the experimental arm
#' dropout_rate <- 0.001
#' fail_rate <- define_fail_rate(
#'   duration = c(delay_effect_duration, 100),
#'   fail_rate = log(2) / median_ctrl,
#'   hr = median_ctrl / median_exp,
#'   dropout_rate = dropout_rate)
#'
#' # Other related parameters
#' alpha <- 0.025 # Type I error
#' beta <- 0.1 # Type II error
#' ratio <- 1 # Randomization ratio (experimental:control)
#'
#' # Build a one-sided group sequential design
#' design <- gs_design_ahr(
#'   enroll_rate = enroll_rate, fail_rate = fail_rate,
#'   ratio = ratio, alpha = alpha, beta = beta,
#'   analysis_time = c(12, 24, 36),
#'   upper = gs_spending_bound,
#'   upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
#'   lower = gs_b,
#'   lpar = rep(-Inf, 3))
#'
#' # Define cuttings of 2 IAs and 1 FA
#' ia1_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[1]))
#' ia2_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[2]))
#' fa_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[3]))
#'
#' # Run simulations
#' simulation <- sim_gs_n(
#'   n_sim = 3,
#'   sample_size = ceiling(design$analysis$n[3]),
#'   enroll_rate = design$enroll_rate,
#'   fail_rate = design$fail_rate,
#'   test = wlr,
#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
#'   weight = fh(rho = 0, gamma = 0.5))
#'
#' # Summarize simulations
#' bound <- gsDesign::gsDesign(k = 3, test.type = 1, sfu = gsDesign::sfLDOF)$upper$bound
#' simulation |> summary(bound = bound)
#'
#' # Summarize simulation and compare with the planned design
#' simulation |> summary(design = design)
summary.simtrial_gs_wlr <- function(object,
                                    design = NULL,
                                    bound = NULL,
                                    ...) {
  # get the total number of analysis and simulations
  n_analysis <- nrow(object[object$sim_id == 1, ])
  n_sim <- nrow(object) / n_analysis

  object <- as.data.table(object)
  # if the design input is NULL
  # then simply output the simulated n, event, power
  if (is.null(design)) {
    ans1 <- object[,
                   .(
                     sim_n = mean(n),
                     sim_event = mean(event),
                     sim_time = mean(cut_date)
                   ),
                   by = "analysis"]

    bound_dt <- data.table(analysis = 1:n_analysis, upper_bound = bound)
    ans2 <- merge.data.table(object, bound_dt, all.x = TRUE, sort = FALSE)
    ans2[, cross_upper := z >= upper_bound]
    ans2 <- ans2[cross_upper == TRUE, ]
    ans2 <- ans2[, .SD[1], by = "sim_id"]
    ans2 <- ans2[, .(n_cross_upper = .N), by = "analysis"]
    ans2 <- ans2[order(analysis),
                 .(analysis, sim_upper_prob = cumsum(n_cross_upper) / n_sim)]

    ans <- merge.data.table(ans1, ans2, all.x = TRUE)

    attr(ans, "compare_with_design") <- "no"
  } else {
    # get the design type, 1-sided or 2-sided
    design_type <- if(length(unique(design$bound$bound)) == 1) "one-sided" else "two-sided"

    # add the futility and efficacy bounds to the simulation results
    if (design_type == "one-sided") {
      bound_dt <- as.data.table(design$bound)
      bound_dt <- bound_dt[, .(analysis, upper_bound = z, bound)]
      sim_tbl <- merge.data.table(object, bound_dt, all.x = TRUE, sort = FALSE)
      sim_tbl[, cross_upper := z >= upper_bound]
    } else {
      bound_dt <- as.data.table(design$bound)
      bound_dt <- dcast(bound_dt, analysis ~ bound, value.var = "z")
      bound_dt <- bound_dt[, .(analysis, lower_bound = lower, upper_bound = upper)]

      sim_tbl <- merge.data.table(object, bound_dt, all.x = TRUE, sort = FALSE)
      sim_tbl[, cross_lower := z <= lower_bound]
      sim_tbl[, cross_upper := z >= upper_bound]
    }

    # calculate the prob of crossing efficacy bounds
    tbl_upper <- sim_tbl[cross_upper == TRUE, ]
    tbl_upper <- tbl_upper[, .SD[1], by = "sim_id"]
    tbl_upper <- tbl_upper[, .(n_cross_upper = .N), by = "analysis"]
    tbl_upper <- tbl_upper[order(analysis),
                           .(analysis, sim_upper_prob = cumsum(n_cross_upper) / n_sim)]

    # calculate the prob of crossing futility bounds
    if (design_type == "two-sided") {
      tbl_lower <- sim_tbl[cross_lower == TRUE, ]
      tbl_lower <- tbl_lower[, .SD[1], by = "sim_id"]
      tbl_lower <- tbl_lower[, .(n_cross_lower = .N), by = "analysis"]
      tbl_lower <- tbl_lower[order(analysis),
                             .(analysis, sim_lower_prob = cumsum(n_cross_lower) / n_sim)]
    }

    # combining prob of crossing efficacy and futility bounds under H1
    if (design_type == "one-sided") {
      tbl_asy_prob <- as.data.table(design$bound)
      tbl_asy_prob <- tbl_asy_prob[, .(analysis, asy_upper_prob = probability)]
    } else {
      tbl_asy_prob <- as.data.table(design$bound)
      tbl_asy_prob <- dcast(tbl_asy_prob, analysis ~ bound, value.var = "probability")
      tbl_asy_prob <- tbl_asy_prob[, .(analysis, asy_upper_prob = upper, asy_lower_prob = lower)]
    }

    # calculate the number of analysis time, events and sample size
    tbl_event <- object[,
                        .(
                          sim_event = mean(event),
                          sim_n = mean(n),
                          sim_time = mean(cut_date)
                        ),
                        by = "analysis"]
    analysis_dt <- as.data.table(design$analysis)
    analysis_dt <- analysis_dt[,
                               .(
                                 analysis,
                                 asy_time = time,
                                 asy_n = n,
                                 asy_event = event
                               )]
    tbl_event <- merge.data.table(tbl_event, analysis_dt, all.y = TRUE, sort = FALSE)
    # combine all the information together
    if (design_type == "one-sided") {
      ans <- tbl_asy_prob |>
        merge.data.table(tbl_upper, all.x = TRUE, sort = FALSE) |>
        merge.data.table(tbl_event, all.x = TRUE, sort = FALSE)
    } else {
      ans <- tbl_asy_prob |>
        merge.data.table(tbl_upper, all.x = TRUE, sort = FALSE) |>
        merge.data.table(tbl_lower, all.x = TRUE, sort = FALSE) |>
        merge.data.table(tbl_event, all.x = TRUE, sort = FALSE)
    }

    attr(ans, "compare_with_design") <- "yes"
    attr(ans, "design_type") <- design_type
  }

  ans <- as.data.frame(ans)
  class(ans) <- c("simtrial_gs_wlr", class(ans))
  attr(ans, "method") <- attributes(object)$method

  return(ans)
}