---
title: "Speed comparison"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Speed comparison}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  eval     = FALSE
)
```

# Purpose

FastSurvival is designed for repeated evaluation inside large simulation
loops. This vignette shows how to benchmark each estimation and testing
function against an established reference and reports representative results.
The benchmark code is shown but not executed when the vignette is built,
because timing many microbenchmark replicates would exceed the build-time
limits. To reproduce the numbers, run the code blocks interactively. The same
code is collected in the `tools/benchmark_speed.R` script of the package's
GitHub repository.

The reported figures are median times of 1,000 microbenchmark replicates,
measured on 2026-09-29 with R 4.6.0 on Windows 11 (x86_64) for the data set of
500 subjects below, all of whom have an event. The FastSurvival functions are
timed on presorted input, so the single sort of the data is excluded from
their times, whereas the reference functions sort internally; in a simulation
loop the sort is paid once per data set. `coxph_fast()` computes a closed-form
approximation of the Cox estimate (the Pike-Halley Estimator) rather than the
iterative maximum partial likelihood estimate, so its row compares two
estimators of the same quantity. Absolute timings depend on hardware, sample
size, and event rate, so the ratios matter more than the raw values.

```{r load}
library(FastSurvival)
library(survival)
library(microbenchmark)
```

# Setup

The key to the speed gain is that the analysis functions accept pre-sorted
vectors. Inside a simulation loop the data are sorted once and reused, so the
sort cost is paid a single time rather than on every call. We build a single
two-group dataset of 500 subjects with `simdata_fast()` and prepare the
sorted vectors, the binary arm indicator, and the restriction horizon used by
the time-restricted methods.

```{r data}
dataset <- simdata_fast(
  nsim     = 1,
  n        = 500,
  a.time   = c(0, 12.5),
  a.rate   = 40,
  e.median = list(5.811, 4.3),
  seed     = 1
)

# Sort once and reuse, the intended pattern for the pre-sorted fast path.
ord <- order(dataset$tte)
t_s <- dataset$tte[ord]
e_s <- dataset$event[ord]
g_s <- dataset$group[ord]

# Control is group 1, treatment is group 2.
arm <- as.integer(dataset$group == 2)

# Restriction horizon within both arms' follow-up.
tau <- floor(min(tapply(t_s, g_s, max)))

# Factor arm for the nphRCT reference used in the rmw_fast benchmark.
df_rmw <- data.frame(
  tte   = dataset$tte,
  event = dataset$event,
  arm   = factor(ifelse(dataset$group == 1, "control", "treatment"),
                 levels = c("control", "treatment"))
)
```

# survfit_fast vs survfit + summary

```{r bench-survfit}
microbenchmark(
  fast = survfit_fast(t_s, e_s, t_eval = tau, presorted = TRUE),
  ref  = summary(survfit(Surv(tte, event) ~ 1, data = dataset), times = tau),
  times = 1000
)
```

# survdiff_fast vs survdiff

```{r bench-survdiff}
microbenchmark(
  fast = survdiff_fast(t_s, e_s, g_s, control = 1, side = 1, presorted = TRUE),
  ref  = survdiff(Surv(tte, event) ~ group, data = dataset),
  times = 1000
)
```

# coxph_fast vs coxph

```{r bench-coxph}
microbenchmark(
  fast = coxph_fast(t_s, e_s, g_s, control = 1, side = 1, presorted = TRUE),
  ref  = coxph(Surv(tte, event) ~ I(group == 2), data = dataset),
  times = 1000
)
```

# rmst_fast vs survRM2::rmst2

```{r bench-rmst}
microbenchmark(
  fast = rmst_fast(t_s, e_s, g_s, control = 1, tau = tau, side = 1,
                   presorted = TRUE),
  ref  = survRM2::rmst2(time = dataset$tte, status = dataset$event,
                        arm = arm, tau = tau),
  times = 1000
)
```

# survdiff_fast(weight = "fh") vs nph::logrank.test

```{r bench-wlr}
microbenchmark(
  fast = survdiff_fast(t_s, e_s, g_s, control = 1, side = 1,
                       weight = "fh", rho = 0, gamma = 1, presorted = TRUE),
  ref  = nph::logrank.test(dataset$tte, dataset$event, dataset$group,
                           rho = 0, gamma = 1),
  times = 1000
)
```

# wmst_fast vs survWMST::wmst

The window mean survival time is benchmarked against `wmst()` from the
survWMST package. survWMST is distributed on GitHub (pauknemj/survWMST), not
CRAN, so this benchmark is shown as a static block rather than a live chunk,
and the vignette carries no undeclared dependency. Install survWMST with
`remotes::install_github("pauknemj/survWMST")` and run the block to reproduce
it.

```r
microbenchmark(
  fast = wmst_fast(t_s, e_s, g_s, control = 1, tau1 = 0, tau2 = tau,
                   side = 1, presorted = TRUE),
  ref  = survWMST::wmst(time = dataset$tte, status = dataset$event,
                        arm = arm, tau0 = 0, tau1 = tau),
  times = 1000
)
```

# milestone_fast vs survfit + summary

```{r bench-milestone}
microbenchmark(
  fast = milestone_fast(t_s, e_s, g_s, control = 1, tau = tau, side = 1,
                        presorted = TRUE),
  ref  = summary(survfit(Surv(tte, event) ~ group, data = dataset),
                 times = tau),
  times = 1000
)
```

# medsurv_fast vs nph::nphparams

```{r bench-medsurv}
microbenchmark(
  fast = medsurv_fast(t_s, e_s, g_s, control = 1, side = 1,
                      method = "nph", presorted = TRUE),
  ref  = nph::nphparams(time = dataset$tte, event = dataset$event,
                        group = as.integer(dataset$group == 2),
                        param_type = "Q", param_par = 0.5),
  times = 1000
)
```

# maxcombo_fast vs nph::logrank.maxtest

```{r bench-maxcombo}
microbenchmark(
  fast = maxcombo_fast(t_s, e_s, g_s, control = 1, side = 1,
                       rho = c(0, 0, 1), gamma = c(0, 1, 0), presorted = TRUE),
  ref  = nph::logrank.maxtest(dataset$tte, dataset$event,
                              as.integer(dataset$group == 2)),
  times = 1000
)
```

# rmw_fast vs nphRCT::wlrt

`rmw_fast()` combines a standard and a modestly-weighted log-rank statistic,
so the reference computes both weighted log-rank components with nphRCT.

```{r bench-rmw}
microbenchmark(
  fast = rmw_fast(t_s, e_s, g_s, control = 1, side = 1, s_star = 0.5,
                  presorted = TRUE),
  ref  = {
    nphRCT::wlrt(Surv(tte, event) ~ arm, data = df_rmw,
                 method = "mw", s_star = 1)
    nphRCT::wlrt(Surv(tte, event) ~ arm, data = df_rmw,
                 method = "mw", s_star = 0.5)
  },
  times = 1000
)
```

# wkm_fast vs nphsim::wkm.Stat

The weighted Kaplan-Meier (Pepe-Fleming) test is benchmarked against
`wkm.Stat()` from the nphsim package. nphsim is distributed on GitHub
(keaven/nphsim), not CRAN, so this benchmark is shown as a static block.
Install nphsim with `remotes::install_github("keaven/nphsim")` and run the
block to reproduce it.

```r
microbenchmark(
  fast = wkm_fast(t_s, e_s, g_s, control = 1, side = 1, weight = "PF",
                  presorted = TRUE),
  ref  = nphsim::wkm.Stat(survival = dataset$tte, cnsr = 1 - dataset$event,
                          trt = ifelse(dataset$group == 1,
                                       "control", "experimental")),
  times = 1000
)
```

# ahsw_fast vs survAH::ah2

```{r bench-ahsw}
microbenchmark(
  fast = ahsw_fast(t_s, e_s, g_s, control = 1, tau = tau, side = 1,
                   presorted = TRUE),
  ref  = survAH::ah2(time = dataset$tte, status = dataset$event,
                     arm = arm, tau = tau),
  times = 1000
)
```

# ahr_fast vs AHR::ahrKM

The Kalbfleisch-Prentice average hazard ratio is benchmarked against `ahrKM()`
from the AHR package (Dormuth et al. (2024) computed the average hazard ratio
with the AHR package). Because AHR has been archived on CRAN, this benchmark is shown as a
static block rather than a live chunk. Install AHR with
`remotes::install_github("cran/AHR")` and run the block to reproduce it.

```r
microbenchmark(
  fast = ahr_fast(t_s, e_s, g_s, control = 1, tau = tau, side = 1,
                  presorted = TRUE),
  ref  = AHR::ahrKM(tau, Surv(tte, event) ~ group, dataset),
  times = 1000
)
```

# Representative results

The table below summarizes representative median timings on the n = 500
two-group dataset generated above, with `presorted = TRUE` and one-sided tests
where applicable. The exact values will differ on your machine, but the order
of magnitude of the speedup is stable. The `wkm_fast()` row is missing because
nphsim was not installed when the table was produced.

| Function | Replaces | Approximate speed gain |
|----------|----------|------------------------|
| `survfit_fast()` | `survfit()` + `summary()` at one time point | ~40x |
| `survdiff_fast()` | `survdiff()` | ~25x |
| `coxph_fast()` | `coxph()` (point estimate + Wald CI) | ~35x |
| `rmst_fast()` | `survRM2::rmst2()` | ~35x |
| `survdiff_fast(weight = "fh")` | `nph::logrank.test()` | ~300x |
| `wmst_fast()` | `survWMST::wmst()` | ~900x |
| `milestone_fast()` | `survfit()` + `summary()` at a milestone | ~20x |
| `medsurv_fast()` | `nph::nphparams()` | ~35x |
| `maxcombo_fast()` | `nph::logrank.maxtest()` | ~350x |
| `rmw_fast()` | `nphRCT::wlrt()` (two components) | ~75x |
| `ahsw_fast()` | `survAH::ah2()` | ~450x |
| `ahr_fast()` | `AHR::ahrKM()` | ~200x |

# Why it is faster

Each function avoids the overhead that the standard implementations incur on
every call. The standard functions parse a formula, build an S3 model object,
and construct intermediate vectors before producing the result, which is
appropriate for interactive use but wasteful when the same operation is
repeated thousands of times. The FastSurvival functions take plain vectors,
do the core computation in a single C++ pass over the data, and return a
lightweight numeric vector. When the input is already sorted the sort cost is
avoided entirely. In a simulation loop these savings accumulate across every
iteration.

# Whole simulation studies

The per-call gains carry over to complete simulation studies. The scripts in
`tools/paper/` of the package's GitHub repository run the same designs with
FastSurvival and with other simulation packages and record the operating
characteristics and the elapsed time. For a two-arm group-sequential design
with 600 subjects and two event-driven log-rank analyses, the time per
simulated trial was about two orders of magnitude longer with simtrial and
with TrialSimulator than with FastSurvival, and for a crossover after a
positive progression-free survival analysis in an illness-death model it was
more than an order of magnitude longer with TrialSimulator. The ratios varied
between runs on the same computer. The power and the analysis times agreed
within Monte Carlo error. In a separate benchmark with the same sample size and
looks and a delayed treatment effect, generating 10,000 trials and computing
the log-rank and RMST statistics at the two looks with FastSurvival each took
on the order of one second. The max-combo p-values, one multivariate normal integral
per trial and look, take most of the computing time; when only the decisions at
given nominal levels are needed, the `mc.alpha` argument of `analysis_fast()`
restricts the integration to the p-values near those levels.

# References

Dormuth, I., Pauly, M., Rauch, G., & Herrmann, C. (2024). Sample size
calculation under nonproportional hazards using average hazard ratios.
*Biometrical Journal*, 66(6), e202300271.