OptSurvCutR: Validated Cut-point Selection for Survival Analysis

R-CMD-check Lifecycle: Stable Codecov License: GPL-3 OptSurvCutR logo

OptSurvCutR finds data-driven, covariate-adjusted thresholds for stratifying patients in survival (time-to-event) data — for biomarkers, gene expression, or any continuous predictor — and tells you whether those thresholds would hold up in a different sample. It is developed against rOpenSci software standards and is currently under rOpenSci review.

Most threshold-selection tools stop once a split is statistically significant. OptSurvCutR’s third step is what happens next: does the same threshold reappear under resampling, or was it a feature of this one dataset?

Why OptSurvCutR?

Feature Benefit
Optimal number of cuts AIC, AICc, or BIC choose between 0 and k cut-points.
Covariate adjustment Thresholds are selected conditional on clinical confounders, not marginally.
Bootstrap stability grade Every threshold gets a 95% CI and an automated Tier 1–4 grade.
Schoenfeld diagnostics Proportional hazards is checked automatically alongside the fit.
Two search engines A deterministic grid, or a multithreaded genetic search (rgenoud).
Publication-ready plots KM curves, distribution splits, forest plots, objective surfaces, 2D stability maps.

Installation

install.packages("OptSurvCutR")

# development version
if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes")
remotes::install_github("paytonyau/OptSurvCutR")

Example: Serum Bilirubin in Primary Biliary Cholangitis

Using the Mayo Clinic PBC trial cohort shipped with survival, reproducing the manuscript’s analysis:

library(OptSurvCutR)
library(survival)

data("pbc", package = "survival")
pbc_clean <- na.omit(pbc[, c("time", "status", "bili", "age", "sex", "edema")])

# Composite endpoint: transplant-free survival (0 = censored, 1 = transplant, 2 = death)
pbc_clean$event <- as.integer(pbc_clean$status %in% c(1, 2))

# 1. How many cut-points does the data support?
num_res <- find_cutpoint_number(
  data = pbc_clean, predictor = "bili",
  outcome_time = "time", outcome_event = "event",
  covariates = c("age", "sex", "edema"),
  method = "genetic", criterion = "AIC",
  max_cuts = 5, nmin = 0.15, seed = 123
)
summary(num_res)   # AIC minimised at 3 cut-points

# 2. Where are the boundaries, adjusting for covariates?
cut_res <- find_cutpoint(
  data = pbc_clean, predictor = "bili",
  outcome_time = "time", outcome_event = "event",
  covariates = c("age", "sex", "edema"),
  method = "genetic", criterion = "logrank",
  num_cuts = num_res$optimal_num_cuts,   # carried forward, not re-chosen
  nmin = 0.15, n_perm = 1000, n_cores = 2, seed = 123
)
summary(cut_res)   # thresholds, adjusted hazard ratios, Schoenfeld test

# 3. Are those boundaries reproducible?
val_res <- validate_cutpoint(
  cutpoint_result = cut_res,
  num_replicates = 500, n_cores = 2, seed = 123
)
summary(val_res)                       # Tier 3 (OVERLAPPING): adjacent intervals overlap
summary(val_res)$stability$tier        # same result, programmatically

# 4. Act on it - a single threshold proves highly reproducible
cut_1 <- find_cutpoint(
  data = pbc_clean, predictor = "bili",
  outcome_time = "time", outcome_event = "event",
  covariates = c("age", "sex", "edema"),
  method = "systematic", criterion = "logrank",
  num_cuts = 1, nmin = 0.15, n_perm = 1000, n_cores = 2, seed = 123
)
summary(validate_cutpoint(cut_1, num_replicates = 500, n_cores = 2, seed = 123))
# Tier 1: threshold 2.3 mg/dL, bootstrap interval 2.0-3.4

# Visualisation
plot(cut_res, type = "distribution")            # predictor density with thresholds
plot(cut_res, type = "outcome")                 # adjusted Kaplan-Meier curves
plot(cut_res, type = "forest")                  # adjusted hazard ratios
plot(cut_res, type = "diagnostic")               # Schoenfeld residual check
plot_validation(val_res, focus_cuts = c(1, 2))   # joint bootstrap distribution

What this shows. The 4-group model has 98.4% of the AIC weight and hazard ratios up to 20.9 - it would pass any significance test. Bootstrap validation showed its boundaries don’t survive resampling. A single threshold at 2.3 mg/dL does. Significance and stability are different properties; step 3 is what tells them apart.

Full analysis: vignette("bilirubin", package = "OptSurvCutR").

Workflow

  1. find_cutpoint_number() - how many thresholds does the data support (AIC/AICc/BIC)?
  2. find_cutpoint() - where are they, adjusted for covariates? Includes Schoenfeld diagnostics.
  3. validate_cutpoint() - do they reproduce under bootstrap resampling? Returns a Tier 1-4 grade.

What’s New

v0.11.1 - two correctness fixes, one label change, one new argument

If you used either combination below under v0.11.0, re-run with v0.11.1.

Full technical detail in NEWS.md.

v0.11.0 - stability metric fix

Recompute any Tier reported from v0.10.1 or earlier. validate_cutpoint() didn’t carry the predictor into its result, so relative CI width was computed against the wrong quantity - widths were inflated and Tier 1 was unreachable for any input.

Contributing

Bug reports and pull requests are welcome - see CONTRIBUTING.md for code standards and local test workflows.

Resources

Citation

@article{yau2025optsurvcutr,
  author    = {Yau, Payton T. O.},
  title     = {OptSurvCutR: Validated Cut-point Selection for Survival Analysis},
  year      = {2025},
  doi       = {10.1101/2025.10.08.681246},
  publisher = {Cold Spring Harbor Laboratory},
  journal   = {bioRxiv},
  url       = {https://www.biorxiv.org/content/10.1101/2025.10.08.681246}
}

License

GPL-3.

Contact

Questions, feature requests, or bug reports: open an issue.