---
title: "Using winratiosim"
author: "Se Yoon Lee"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using winratiosim}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```

## Overview

`winratiosim` simulates operating characteristics for two-arm clinical trials
with a hierarchical win ratio endpoint. A simulated trial includes three
prioritized outcome layers:

1. time to death,
2. annualized recurrent event count, and
3. a continuous quality-of-life change score.

The package was created for the simulation workflow used in Lee (2025), which
compares Finkelstein-Schoenfeld permutation variance calculations with the
large-sample variance formula discussed by Yu and Ganju.

## Quick start

The main function is `winratiosim()`. The example below uses only two
simulated trials and a small sample size so that the vignette builds quickly.
Increase `nsim` and `N` for a real operating-characteristic study.

```{r quick-start}
library(winratiosim)

quick_res <- winratiosim(
  nsim = 2,
  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 = 20250518
)

quick_res$df_WR.analysis.summary
quick_res$df_sample.size.summary
```

The returned object is a named list:

```{r result-names}
names(quick_res)
```

The most commonly used elements are:

* `df_FS.analysis.summary`: Finkelstein-Schoenfeld statistic, variance,
  z-score, and p-value for each simulated trial.
* `df_WR.analysis.summary`: win ratio estimates, confidence limits, variance
  estimates, and p-values for each simulated trial.
* `df_Total_probability`: treatment win, tie, and control win probabilities.
* `df_sample.size.summary`: treatment and control sample sizes generated under
  the requested randomization ratio.

## Estimating power

For a one-sided superiority analysis at level 0.025, one common summary is the
proportion of simulated trials with a significant result. Exact binomial
confidence intervals can be calculated with `binom.conf.exact()`.

```{r power-summary}
fs_success <- quick_res$df_FS.analysis.summary$p_value_FS < 0.025
wr_success <- quick_res$df_WR.analysis.summary$LB_R_w > 1

data.frame(
  Method = c("FS test", "YG win ratio test"),
  Estimated_power = c(
    mean(fs_success, na.rm = TRUE),
    mean(wr_success, na.rm = TRUE)
  )
)

binom.conf.exact(
  x = sum(wr_success, na.rm = TRUE),
  n = sum(!is.na(wr_success))
)
```

This small example is intended only to show the workflow. Power estimates from
two simulations are not scientifically meaningful.

## Paper-style simulation workflow

The following code mirrors the larger simulation workflow used for the paper.
It is not evaluated when this vignette is built because `nsim = 10000` can take
substantial time.

```{r paper-workflow, eval = FALSE}
library(winratiosim)

power.design_parameters <- list(
  nsim = 10000,
  N = 400,
  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 = 45 + 7.5,
  sd.delta.x_trt = 20,
  xbase_ctl = 45,
  xfinal_ctl = 45,
  sd.delta.x_ctl = 20,
  censorrate_trt = 0.2,
  censorrate_ctl = 0.2,
  nc = 10,
  seed = 20250518
)

power.sim_res <- do.call(winratiosim, power.design_parameters)

Power_binom_CI_one_sided_FS_Permutation <- binom.conf.exact(
  x = sum(power.sim_res$df_FS.analysis.summary$p_value_FS < 0.025,
          na.rm = TRUE),
  n = sum(!is.na(power.sim_res$df_FS.analysis.summary$p_value_FS))
)

Power_binom_CI_one_sided_WR_Ron_Yu <- binom.conf.exact(
  x = sum(power.sim_res$df_WR.analysis.summary$LB_R_w > 1,
          na.rm = TRUE),
  n = sum(!is.na(power.sim_res$df_WR.analysis.summary$LB_R_w))
)

t1e.design_parameters <- list(
  nsim = power.design_parameters$nsim,
  N = power.design_parameters$N,
  Randomization.ratio = power.design_parameters$Randomization.ratio,
  alpha.JFM = power.design_parameters$alpha.JFM,
  theta.JFM = power.design_parameters$theta.JFM,
  lambda_trt = power.design_parameters$lambda_ctl,
  lambda_ctl = power.design_parameters$lambda_ctl,
  ann.icr_trt = power.design_parameters$ann.icr_ctl,
  ann.icr_ctl = power.design_parameters$ann.icr_ctl,
  xbase_trt = power.design_parameters$xbase_ctl,
  xfinal_trt = power.design_parameters$xfinal_ctl,
  sd.delta.x_trt = power.design_parameters$sd.delta.x_trt,
  xbase_ctl = power.design_parameters$xbase_ctl,
  xfinal_ctl = power.design_parameters$xfinal_ctl,
  sd.delta.x_ctl = power.design_parameters$sd.delta.x_ctl,
  censorrate_trt = power.design_parameters$censorrate_trt,
  censorrate_ctl = power.design_parameters$censorrate_ctl,
  nc = power.design_parameters$nc,
  seed = 20250518
)

t1e.sim_res <- do.call(winratiosim, t1e.design_parameters)

t1e_binom_CI_one_sided_FS_Permutation <- binom.conf.exact(
  x = sum(t1e.sim_res$df_FS.analysis.summary$p_value_FS < 0.025,
          na.rm = TRUE),
  n = sum(!is.na(t1e.sim_res$df_FS.analysis.summary$p_value_FS))
)

t1e_binom_CI_one_sided_WR_Ron_Yu <- binom.conf.exact(
  x = sum(t1e.sim_res$df_WR.analysis.summary$LB_R_w > 1,
          na.rm = TRUE),
  n = sum(!is.na(t1e.sim_res$df_WR.analysis.summary$LB_R_w))
)

df.power.type1 <- data.frame(
  Method = c("FS test", "YG test"),
  Power = paste(
    round(c(Power_binom_CI_one_sided_FS_Permutation[1],
            Power_binom_CI_one_sided_WR_Ron_Yu[1]), 3),
    "(",
    round(c(Power_binom_CI_one_sided_FS_Permutation[2],
            Power_binom_CI_one_sided_WR_Ron_Yu[2]), 3),
    ", ",
    round(c(Power_binom_CI_one_sided_FS_Permutation[3],
            Power_binom_CI_one_sided_WR_Ron_Yu[3]), 3),
    ")",
    sep = ""
  ),
  Type_I_Error = paste(
    round(c(t1e_binom_CI_one_sided_FS_Permutation[1],
            t1e_binom_CI_one_sided_WR_Ron_Yu[1]), 3),
    "(",
    round(c(t1e_binom_CI_one_sided_FS_Permutation[2],
            t1e_binom_CI_one_sided_WR_Ron_Yu[2]), 3),
    ", ",
    round(c(t1e_binom_CI_one_sided_FS_Permutation[3],
            t1e_binom_CI_one_sided_WR_Ron_Yu[3]), 3),
    ")",
    sep = ""
  )
)

df.variance <- data.frame(
  Median_Variance_under_Power = c(
    median(power.sim_res$df_WR.analysis.summary$variance_log_R_w_permutation,
           na.rm = TRUE),
    median(power.sim_res$df_WR.analysis.summary$Var_logR_w,
           na.rm = TRUE)
  ),
  Median_Variance_under_Type_I_Error = c(
    median(t1e.sim_res$df_WR.analysis.summary$variance_log_R_w_permutation,
           na.rm = TRUE),
    median(t1e.sim_res$df_WR.analysis.summary$Var_logR_w,
           na.rm = TRUE)
  )
)

df.combined <- cbind(df.power.type1, round(df.variance, 4))
df.combined

median(power.sim_res$df_WR.analysis.summary$R_w, na.rm = TRUE)
median(power.sim_res$df_Total_probability[, "Prob_of_tie"], na.rm = TRUE)
```

## Output from the 10,000-simulation run

When the paper-style workflow above is run with `nsim = 10000`, `N = 400`,
`nc = 10`, and `seed = 20250518`, the final summary commands produce the
following output. The full simulation is not rerun during vignette building.

```r
df.power.type1
#>    Method               Power        Type_I_Error
#> 1 FS test  0.86(0.853, 0.866) 0.024(0.021, 0.028)
#> 2 YG test 0.811(0.803, 0.819) 0.018(0.015, 0.021)

df.variance
#>   Median_Variance_under_Power Median_Variance_under_Type_I_Error
#> 1                  0.01677787                         0.01607165
#> 2                  0.01969007                         0.01882680

median(power.sim_res$df_WR.analysis.summary$R_w, na.rm = TRUE)
#> [1] 1.474161

median(power.sim_res$df_Total_probability[, "Prob_of_tie"], na.rm = TRUE)
#> [1] 0.191
```

## Parameter notes

* `lambda_trt` and `lambda_ctl` are annual mortality probabilities.
* `ann.icr_trt` and `ann.icr_ctl` are annual recurrent event incidence rates.
* `xbase_*` and `xfinal_*` define the mean continuous outcome change in each
  arm.
* `censorrate_*` gives the annual censoring probability.
* `nc` controls the number of worker processes. Use `nc = 1` when debugging.
* `seed` makes the simulation reproducible.

## References

Lee, S. Y. (2025). A note on the sample size formula for a win ratio endpoint.
*Statistics in Medicine*, 44, e70165. <https://doi.org/10.1002/sim.70165>

Finkelstein, D. M., and Schoenfeld, D. A. (1999). Combining mortality and
longitudinal measures in clinical trials. *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. *European Heart Journal*, 33(2), 176-182.

Yu, R. X., and Ganju, J. (2022). Sample size formula for a win ratio endpoint.
*Statistics in Medicine*, 41(6), 950-963.