R/winratiosim.R

In winratiosim: Simulate Power for Hierarchical Win Ratio Endpoints

#' Simulate Hierarchical Win Ratio Trials
#'
#' Simulates replicated two-arm clinical trials and analyzes each trial using a
#' three-layer hierarchical win ratio framework: time to death, annualized
#' recurrent event count, and a continuous quality-of-life score.
#'
#' @param nsim Integer. Number of simulated trials.
#' @param N Integer. Total number of subjects in each simulated trial.
#' @param Randomization.ratio Numeric vector of length 2 giving the treatment
#'   and control allocation ratio, for example \code{c(1, 1)}.
#' @param alpha.JFM Numeric. Alpha parameter for the joint frailty model.
#' @param theta.JFM Numeric. Frailty variance parameter for the joint frailty
#'   model. Must be positive.
#' @param lambda_trt,lambda_ctl Numeric. Annual mortality probabilities for the
#'   treatment and control arms.
#' @param ann.icr_trt,ann.icr_ctl Numeric. Annual recurrent event incidence
#'   rates for the treatment and control arms.
#' @param xbase_trt,xfinal_trt Numeric. Baseline and expected final continuous
#'   outcome values in the treatment arm.
#' @param xbase_ctl,xfinal_ctl Numeric. Baseline and expected final continuous
#'   outcome values in the control arm.
#' @param sd.delta.x_trt,sd.delta.x_ctl Numeric. Standard deviations for the
#'   continuous outcome change in the treatment and control arms.
#' @param censorrate_trt,censorrate_ctl Numeric. Annual censoring probabilities
#'   for the treatment and control arms.
#' @param nc Integer. Number of worker processes to use. The default is 1.
#' @param seed Optional integer seed. If supplied, results are reproducible
#'   across different values of \code{nc}.
#'
#' @return A named list with the following elements:
#' \describe{
#'   \item{df_FS.analysis.summary}{Finkelstein-Schoenfeld analysis summary for
#'   each simulation.}
#'   \item{df_WR.analysis.summary}{Win ratio analysis summary for each
#'   simulation.}
#'   \item{df_sample.size.summary}{Sample sizes used in each simulated trial.}
#'   \item{df_Total_probability}{Win, tie, loss, and total probabilities for
#'   each simulation.}
#'   \item{df_Total_count}{Win, tie, loss, and total counts for each
#'   simulation.}
#' }
#'
#' @examples
#' result <- winratiosim(
#'   nsim = 1,
#'   N = 20,
#'   Randomization.ratio = c(1, 1),
#'   alpha.JFM = 0,
#'   theta.JFM = 1,
#'   lambda_trt = 0.13,
#'   lambda_ctl = 0.15,
#'   ann.icr_trt = 0.32,
#'   ann.icr_ctl = 0.55,
#'   xbase_trt = 45,
#'   xfinal_trt = 52.5,
#'   xbase_ctl = 45,
#'   xfinal_ctl = 45,
#'   sd.delta.x_trt = 20,
#'   sd.delta.x_ctl = 20,
#'   censorrate_trt = 0.2,
#'   censorrate_ctl = 0.2,
#'   nc = 1,
#'   seed = 2025
#' )
#' result$df_WR.analysis.summary
#'
#' @references
#' Lee, S. Y. (2025). A note on the sample size formula for a win ratio
#' endpoint. \emph{Statistics in Medicine}, 44, e70165.
#' \doi{10.1002/sim.70165}
#'
#' @export
winratiosim <- function(nsim, N, Randomization.ratio, alpha.JFM, theta.JFM,
                        lambda_trt, lambda_ctl, ann.icr_trt, ann.icr_ctl,
                        xbase_trt, xfinal_trt, xbase_ctl, xfinal_ctl,
                        sd.delta.x_trt, sd.delta.x_ctl,
                        censorrate_trt, censorrate_ctl, nc = 1,
                        seed = NULL) {
  validate_integer <- function(x, name, lower) {
    if (length(x) != 1L || !is.finite(x) || x < lower || x != as.integer(x)) {
      stop(name, " must be an integer greater than or equal to ", lower, ".",
           call. = FALSE)
    }
    as.integer(x)
  }
  validate_probability <- function(x, name) {
    if (length(x) != 1L || !is.finite(x) || x < 0 || x >= 1) {
      stop(name, " must be a probability in [0, 1).", call. = FALSE)
    }
  }
  validate_nonnegative <- function(x, name) {
    if (length(x) != 1L || !is.finite(x) || x < 0) {
      stop(name, " must be a non-negative number.", call. = FALSE)
    }
  }

  nsim <- validate_integer(nsim, "nsim", 1L)
  N <- validate_integer(N, "N", 2L)
  nc <- validate_integer(nc, "nc", 1L)
  nc <- min(nc, nsim)

  if (length(Randomization.ratio) != 2L ||
      any(!is.finite(Randomization.ratio)) ||
      any(Randomization.ratio <= 0)) {
    stop("Randomization.ratio must contain two positive finite values.",
         call. = FALSE)
  }
  if (length(theta.JFM) != 1L || !is.finite(theta.JFM) || theta.JFM <= 0) {
    stop("theta.JFM must be positive.", call. = FALSE)
  }

  validate_probability(lambda_trt, "lambda_trt")
  validate_probability(lambda_ctl, "lambda_ctl")
  validate_probability(censorrate_trt, "censorrate_trt")
  validate_probability(censorrate_ctl, "censorrate_ctl")
  validate_nonnegative(ann.icr_trt, "ann.icr_trt")
  validate_nonnegative(ann.icr_ctl, "ann.icr_ctl")
  validate_nonnegative(sd.delta.x_trt, "sd.delta.x_trt")
  validate_nonnegative(sd.delta.x_ctl, "sd.delta.x_ctl")

  if (!is.null(seed)) {
    seed <- validate_integer(seed, "seed", 1L)
    set.seed(seed)
  }
  trial_seeds <- sample.int(.Machine$integer.max, nsim)

  sim_data <- SimData_per_group
  score_tte <- Scoring_TTE
  score_conti <- Scoring_Conti
  analyze_wr <- WR_analysis

  simulate_one_trial <- function(trial_seed) {
    set.seed(trial_seed)

    allocation_probability <- Randomization.ratio[1] / sum(Randomization.ratio)
    n_treatment <- stats::rbinom(n = 1, size = N,
                                 prob = allocation_probability)
    n_treatment <- max(1L, min(N - 1L, n_treatment))
    n_control <- N - n_treatment

    surv_1 <- sim_data(
      treatment = 1,
      ngroup = n_treatment,
      alpha.JFM = alpha.JFM,
      theta.JFM = theta.JFM,
      ann.icr = ann.icr_trt,
      lambda = lambda_trt,
      censorrate = censorrate_trt,
      xbase = xbase_trt,
      xfinal = xfinal_trt,
      sd.delta.x = sd.delta.x_trt
    )

    surv_0 <- sim_data(
      treatment = 0,
      ngroup = n_control,
      alpha.JFM = alpha.JFM,
      theta.JFM = theta.JFM,
      ann.icr = ann.icr_ctl,
      lambda = lambda_ctl,
      censorrate = censorrate_ctl,
      xbase = xbase_ctl,
      xfinal = xfinal_ctl,
      sd.delta.x = sd.delta.x_ctl
    )

    df_trial <- rbind(surv_1$surv_1, surv_0$surv_0)
    df_trial$subjid[df_trial$treatment == 0] <-
      df_trial$subjid[df_trial$treatment == 0] + 1000
    df_trial$HFH_Annual <- (df_trial$HFH / df_trial$censortime) * 360
    names(df_trial)[names(df_trial) == "subjid"] <- "usubjid"
    df_trial <- df_trial[order(df_trial$usubjid), , drop = FALSE]

    df_base <- df_trial[rep(seq_len(nrow(df_trial)), each = nrow(df_trial)),
                        , drop = FALSE]
    names(df_base) <- paste0(names(df_base), "1")

    df_compare <- df_trial[rep(seq_len(nrow(df_trial)), times = nrow(df_trial)),
                           , drop = FALSE]
    names(df_compare) <- paste0(names(df_compare), "2")

    df_FS_input <- cbind(df_base, df_compare)
    df_FS_input$score <- NA_real_
    df_FS_input$WR_cat <- ""

    df_outcome1 <- score_tte(
      dataset = df_FS_input,
      var1 = "deathdays1",
      var2 = "deathdays2",
      censor1 = "death1",
      censor2 = "death2"
    )

    df_outcome2 <- score_conti(
      dataset = df_outcome1,
      higher_better = "No",
      var1 = "HFH_Annual1",
      var2 = "HFH_Annual2"
    )

    df_outcome3 <- score_conti(
      dataset = df_outcome2,
      higher_better = "Yes",
      var1 = "kccq1",
      var2 = "kccq2"
    )

    layer_cols <- c("usubjid1", "treatment1", "usubjid2", "treatment2",
                    "score")
    df_layer1 <- df_outcome1[, layer_cols, drop = FALSE]
    df_layer2 <- df_outcome2[, layer_cols, drop = FALSE]
    df_layer3 <- df_outcome3[, layer_cols, drop = FALSE]

    analyze_wr(
      dataset1 = df_layer1,
      dataset2 = df_layer2,
      dataset3 = df_layer3
    )
  }

  if (nc > 1L) {
    cl <- parallel::makeCluster(nc)
    on.exit(parallel::stopCluster(cl), add = TRUE)
    sim_results <- parallel::parLapply(cl, trial_seeds, simulate_one_trial)
  } else {
    sim_results <- lapply(trial_seeds, simulate_one_trial)
  }

  df_FS_summary <- do.call(
    rbind,
    lapply(sim_results, function(x) x$FS.analysis.summary)
  )
  df_WR_summary <- do.call(
    rbind,
    lapply(sim_results, function(x) x$WR.analysis.summary)
  )
  df_sample_size <- do.call(
    rbind,
    lapply(sim_results, function(x) x$sample.size.summary)
  )
  df_prob_summary <- do.call(
    rbind,
    lapply(sim_results, function(x) x$win.losses.count.summary$Total_probability)
  )
  df_count_summary <- do.call(
    rbind,
    lapply(sim_results, function(x) x$win.losses.count.summary$Total_count)
  )

  colnames(df_prob_summary) <- c(
    "Prob_of_Win_trt",
    "Prob_of_tie",
    "Prob_of_Win_ctl",
    "ALL"
  )
  colnames(df_count_summary) <- c(
    "Count_Win_trt",
    "Count_tie",
    "Count_Win_ctl",
    "ALL"
  )

  list(
    df_FS.analysis.summary = df_FS_summary,
    df_WR.analysis.summary = df_WR_summary,
    df_sample.size.summary = df_sample_size,
    df_Total_probability = as.data.frame(df_prob_summary),
    df_Total_count = as.data.frame(df_count_summary)
  )
}