R/WR_analysis.R

In winratiosim: Simulate Power for Hierarchical Win Ratio Endpoints

#' Perform Hierarchical Win Ratio Analysis
#'
#' Analyzes treatment-control pairwise comparisons across three prioritized
#' outcome layers. The function computes layer-specific win, tie, and loss
#' counts; sample sizes; Finkelstein-Schoenfeld statistics; and win ratio
#' statistics based on permutation and large-sample variance formulas.
#'
#' @param dataset1 Data frame containing pairwise scores for the first,
#'   highest-priority layer.
#' @param dataset2 Data frame containing pairwise scores through the second
#'   layer.
#' @param dataset3 Data frame containing pairwise scores through the third
#'   layer.
#'
#' @return A named list with four elements:
#' \describe{
#'   \item{win.losses.count.summary}{Counts and proportions of treatment wins,
#'   ties, and treatment losses by layer and overall.}
#'   \item{sample.size.summary}{Treatment, control, total, and pairwise
#'   comparison counts.}
#'   \item{FS.analysis.summary}{Finkelstein-Schoenfeld statistic, variance,
#'   z-score, and one-sided p-value.}
#'   \item{WR.analysis.summary}{Win ratio, log win ratio, variance estimates,
#'   confidence limits, and one-sided p-value.}
#' }
#'
#' @examples
#' dataset1 <- data.frame(
#'   usubjid1 = c(1, 1, 2, 2),
#'   treatment1 = c(1, 1, 0, 0),
#'   usubjid2 = c(2, 3, 1, 3),
#'   treatment2 = c(0, 0, 1, 0),
#'   score = c(1, NA, -1, NA)
#' )
#' dataset2 <- dataset1
#' dataset2$score <- c(1, 1, -1, NA)
#' dataset3 <- dataset2
#'
#' WR_analysis(dataset1, dataset2, dataset3)$sample.size.summary
#'
#' @references
#' Finkelstein, D. M., and Schoenfeld, D. A. (1999). Combining mortality and
#' longitudinal measures in clinical trials. \emph{Statistics in Medicine},
#' 18(11), 1341-1354.
#'
#' Pocock, S. J., Ariti, C. A., Collier, T. J., and Wang, D. (2012). The win
#' ratio: a new approach to the analysis of composite endpoints in clinical
#' trials based on clinical priorities. \emph{European Heart Journal}, 33(2),
#' 176-182.
#'
#' Yu, R. X., and Ganju, J. (2022). Sample size formula for a win ratio
#' endpoint. \emph{Statistics in Medicine}, 41(6), 950-963.
#'
#' @export
WR_analysis <- function(dataset1, dataset2, dataset3) {
  required <- c("usubjid1", "treatment1", "usubjid2", "treatment2", "score")
  for (name in c("dataset1", "dataset2", "dataset3")) {
    missing_cols <- setdiff(required, names(get(name)))
    if (length(missing_cols) > 0L) {
      stop(name, " is missing required columns: ",
           paste(missing_cols, collapse = ", "), call. = FALSE)
    }
  }

  N_trt <- length(unique(dataset1$usubjid1[dataset1$treatment1 == 1]))
  N_ctl <- length(unique(dataset1$usubjid1[dataset1$treatment1 == 0]))
  N <- N_trt + N_ctl
  N_comparison_win_ratio <- N_trt * N_ctl

  if (N_trt == 0L || N_ctl == 0L) {
    stop("Both treatment and control subjects are required.", call. = FALSE)
  }

  wr_rows1 <- dataset1$treatment1 == 1 & dataset1$treatment2 == 0
  wr_rows2 <- dataset2$treatment1 == 1 & dataset2$treatment2 == 0
  wr_rows3 <- dataset3$treatment1 == 1 & dataset3$treatment2 == 0

  dataset1.WR <- dataset1[wr_rows1, , drop = FALSE]
  dataset2.WR <- dataset2[wr_rows2, , drop = FALSE]
  dataset3.WR <- dataset3[wr_rows3, , drop = FALSE]

  W_T_layer1 <- sum(dataset1.WR$score == 1, na.rm = TRUE)
  W_C_layer1 <- sum(dataset1.WR$score == -1, na.rm = TRUE)
  W_0_layer1 <- N_comparison_win_ratio - (W_T_layer1 + W_C_layer1)

  W_T_layer2 <- sum(dataset2.WR$score == 1, na.rm = TRUE) - W_T_layer1
  W_C_layer2 <- sum(dataset2.WR$score == -1, na.rm = TRUE) - W_C_layer1
  W_0_layer2 <- W_0_layer1 - (W_T_layer2 + W_C_layer2)

  W_T_layer3 <- sum(dataset3.WR$score == 1, na.rm = TRUE) -
    (W_T_layer1 + W_T_layer2)
  W_C_layer3 <- sum(dataset3.WR$score == -1, na.rm = TRUE) -
    (W_C_layer1 + W_C_layer2)
  W_0_layer3 <- W_0_layer2 - (W_T_layer3 + W_C_layer3)

  W_T <- W_T_layer1 + W_T_layer2 + W_T_layer3
  W_C <- W_C_layer1 + W_C_layer2 + W_C_layer3
  W_0 <- N_comparison_win_ratio - (W_T + W_C)

  T.stat <- W_T - W_C
  Ui <- as.numeric(tapply(dataset3$score, dataset3$usubjid1, sum,
                          na.rm = TRUE))
  V <- ((N_trt * N_ctl) / ((N_trt + N_ctl) * (N_trt + N_ctl - 1))) *
    sum(Ui^2, na.rm = TRUE)
  z <- if (V > 0) T.stat / sqrt(V) else NA_real_
  p_value_FS <- stats::pnorm(z, lower.tail = FALSE)

  R_w <- if (W_C > 0) W_T / W_C else Inf
  log_R_w <- log(R_w)
  variance_log_R_w_permutation <- if (is.finite(log_R_w) &&
                                      is.finite(z) && z != 0) {
    (log_R_w / z)^2
  } else {
    NA_real_
  }
  UB_R_w_permutation <- exp(log_R_w + 1.96 *
                              sqrt(variance_log_R_w_permutation))
  LB_R_w_permutation <- exp(log_R_w - 1.96 *
                              sqrt(variance_log_R_w_permutation))

  P0 <- W_0 / N_comparison_win_ratio
  k <- N_trt / (N_trt + N_ctl)
  sigma <- sqrt(4 * (1 + P0) / (3 * k * (1 - k) * (1 - P0)))
  variance_log_R_w <- (sigma^2) / N
  p_value_R_w <- 1 - stats::pnorm(
    q = log_R_w * sqrt((3 * N_trt * N_ctl * (1 - P0)) /
                        (4 * N * (1 + P0)))
  )
  UB_R_w <- exp(log_R_w + 1.96 * sqrt(variance_log_R_w))
  LB_R_w <- exp(log_R_w - 1.96 * sqrt(variance_log_R_w))

  total_values <- c(W_T, W_0, W_C)
  total_n <- sum(total_values)
  win.losses.count.summary <- data.frame(
    Count = c("Number of wins in Treatment Group",
              "Number of ties",
              "Number of losses in Treatment Group",
              "Sum"),
    First_Layer = c(W_T_layer1, W_0_layer1, W_C_layer1,
                    sum(c(W_T_layer1, W_0_layer1, W_C_layer1))),
    Second_Layer = c(W_T_layer2, W_0_layer2, W_C_layer2,
                     sum(c(W_T_layer2, W_0_layer2, W_C_layer2))),
    Third_Layer = c(W_T_layer3, W_0_layer3, W_C_layer3,
                    sum(c(W_T_layer3, W_0_layer3, W_C_layer3))),
    Total_count = c(total_values, total_n),
    Total_probability = round(c(total_values, total_n) / total_n, 3)
  )

  sample.size.summary <- data.frame(
    N = N,
    N_trt = N_trt,
    N_ctl = N_ctl,
    N_comparison_win_ratio = N_comparison_win_ratio
  )
  FS.analysis.summary <- data.frame(
    T = T.stat,
    V = V,
    z = z,
    p_value_FS = p_value_FS
  )
  WR.analysis.summary <- data.frame(
    R_w = R_w,
    logR_w = log_R_w,
    variance_log_R_w_permutation = variance_log_R_w_permutation,
    LB_R_w_95p_permutation = LB_R_w_permutation,
    UB_R_w_95p_permutation = UB_R_w_permutation,
    Var_logR_w = variance_log_R_w,
    UB_R_w = UB_R_w,
    LB_R_w = LB_R_w,
    p_value_R_w = p_value_R_w
  )

  list(
    win.losses.count.summary = win.losses.count.summary,
    sample.size.summary = sample.size.summary,
    FS.analysis.summary = FS.analysis.summary,
    WR.analysis.summary = WR.analysis.summary
  )
}