Introduction

In survival analysis, a common challenge is to stratify subjects based on a continuous predictor. The conventional approach is often to dichotomise this predictor using a single cut-point, typically the median. However, this method can be statistically weak and may fail to capture more complex, non-linear relationships. The fundamental question is often not just where to make a cut, but how many cuts are statistically justified.

This vignette demonstrates how the OptSurvCutR package provides a comprehensive workflow to solve this problem. We will address the following research question:

Based on simulated data, can we find the optimal temperature threshold(s) that best separate fast and slow germination in rapeseed?

To answer this, we will use the germination dataset, which is included with the package. This dataset was simulated based on the findings of Haj Sghaier et al. (2022) to provide a realistic case study. The data was generated based on 7 core temperature groups (5, 10, 15, 20, 25, 30, and 35°C) and 6 intermediate groups (7.5, 12.5, 17.5, 22.5, 27.5, and 32.5°C) to create a rich, continuous predictor.

Analysis Workflow

This guide covers a complete workflow for answering this question using OptSurvCutR:

  1. Installation & Setup: Get the package and load our toolkit.
  2. Data Exploration: Visualise the raw germination patterns.
  3. Step 1: Discovering the Optimal Number of Cut-points: Use find_cutpoint_number() to determine if 1, 2, or more thresholds are best.
  4. Step 2: Discovering the Optimal Value of Cut-points: Use find_cutpoint() to identify the specific temperature thresholds.
  5. Step 3: Validating the Cut-points: Use validate_cutpoint() to assess the stability of our findings.
  6. Visualisation and Interpretation: Create publication-ready plots to understand and present the results.

1. Installation & Setup

Installing OptSurvCutR

First, let us get the necessary tools. OptSurvCutR is available on CRAN, and the development version is on GitHub.

# Option 1: Install the development version from GitHub for the latest features
# You will need the 'remotes' package first: install.packages("remotes")
# remotes::install_github("paytonyau/OptSurvCutR")
# library("OptSurvCutR")
# Option 2: Install the stable version from CRAN (when available)
# install.packages("OptSurvCutR")

Loading Our Toolkit

Now, we load all the R packages we will need for this analysis. We use pacman to automatically install any missing packages.

# Check if pacman is installed, if not, install it
if (!require("pacman")) install.packages("pacman")

# Use pacman to load/install all required packages
pacman::p_load(
  OptSurvCutR, # The core package for this analysis
  survival,    # For core survival analysis functions (Surv, survfit)
  survminer,   # For plotting survival curves (ggsurvplot, ggforest)
  dplyr,       # For data manipulation (pipes %>%, select, mutate)
  ggplot2,     # For creating custom plots
  patchwork,   # For combining multiple ggplots into one figure
  knitr,       # For creating formatted tables (kable)
  tidyr,       # For data tidying
  cli,         # For stylish console messages
  rgenoud      # Required for method = "genetic"
)

The Rapeseed Germination Data

To create a realistic dataset for demonstrating the package’s capabilities, time-to-event data was simulated based on the findings of Haj Sghaier et al. (2022). The simulation was designed in two parts to generate a rich, continuous predictor with a non-linear relationship to the outcome. First, a core dataset was generated based on parameters reported directly in the source manuscript for seven distinct temperature points (5, 10, 15, 20, 25, 30, and 35°C). These parameters included the germination time window (start and end day), coefficients for a linear growth model (slope and intercept), and the overall germination rate for each temperature. Second, to create a more continuous variable suitable for cut-point optimisation, these parameters were linearly interpolated for six intermediate temperature points (7.5, 12.5, 17.5, 22.5, 27.5, and 32.5°C). The final dataset, germination, is included with the package and was formed by combining these two simulated subsets.

# Load the pre-simulated data included with the package
data("germination", package = "OptSurvCutR")
analysis_data <- germination
head(analysis_data)
##   temperature replicate time    growth germinated
## 1           5         1   22 8.1091549          1
## 2           5         1   18 0.0000000          0
## 3           5         1   15 1.2184871          1
## 4           5         1   22 6.0553687          1
## 5           5         1   19 5.1456349          1
## 6           5         1   14 0.4994083          1

2. Initial Data Exploration

Before modelling, we explore the raw data. The summary table below shows key germination metrics for each of the thirteen temperature groups.

# Create a summary table of the data
summary_table <- analysis_data %>%
  group_by(temperature) %>%
  summarise(
    N_Seeds = n(),
    N_Germinated = sum(germinated),
    Germination_Rate_Pct = mean(germinated) * 100,
    Avg_Time_to_Germinate_Days = mean(time[germinated == 1], na.rm = TRUE)
  ) %>%
  rename(`Temperature (°C)` = temperature)

# Display the table using kable for nice formatting
knitr::kable(summary_table,
      digits = 2,
      caption = "Summary of Germination Outcomes by Temperature Group")
Summary of Germination Outcomes by Temperature Group
Temperature (°C) N_Seeds N_Germinated Germination_Rate_Pct Avg_Time_to_Germinate_Days
5.0 80 70 87.50 18.73
7.5 80 71 88.75 16.08
10.0 80 72 90.00 13.93
12.5 80 73 91.25 11.85
15.0 80 74 92.50 10.64
17.5 80 76 95.00 9.72
20.0 80 78 97.50 9.59
22.5 80 77 96.25 10.18
25.0 80 75 93.75 8.51
27.5 80 74 92.50 9.36
30.0 80 73 91.25 10.75
32.5 80 58 72.50 12.78
35.0 80 44 55.00 14.59

Interpretation: This table confirms our data simulates the expected biological response. The highest germination rates and fastest germination times occur in the 15-25°C range, with performance dropping off at colder and hotter temperatures.

3. The Three-Step Analysis Workflow

Step 1: Determine the Optimal Number of Cut-points

Our first step is to determine how many cut-points the data supports. Forcing a single cut-point might be too simple, while too many might overfit the data. find_cutpoint_number() uses information criteria to provide statistical evidence for this decision.

Choosing Function Arguments

The function’s behaviour is controlled by two key arguments: method and criterion.

Method How It Works Recommendation Rating (Accuracy & Performance)
"systematic" Exhaustive Search: Tests every single possible cut-point. Best for 1 cut-point; guarantees the optimal result. Accuracy: ★★★★★ Performance: ★★★☆☆
"genetic" Evolutionary Search: Uses a smart algorithm to efficiently find a near-perfect solution. Highly recommended for 2+ cut-points. Much faster than the systematic search. Accuracy: ★★★★☆ Performance: ★★★★★
Criterion What It Is Recommendation Rating (Accuracy & Performance)
"AIC" Akaike Information Criterion Balances model fit and complexity. A good general-purpose choice. Accuracy: ★★★★☆ Performance: ★★★★★
"AICc" Corrected AIC A version of AIC with a greater penalty for extra parameters. It is specifically recommended for smaller sample sizes. Accuracy: ★★★★★ Performance: ★★★★★
"BIC" Bayesian Information Criterion Similar to AIC, but applies a stronger penalty for complexity, especially in larger datasets. It tends to favour simpler models. Accuracy: ★★★★★ Performance: ★★★★★

We will use the "genetic" search method and the "BIC" criterion, which is excellent for balancing model fit and complexity.

# --- Step 1: Find the number of cut(s) using BIC ---
cli_h2("Step 1: Find the number of cut(s) using BIC")
number_result_bic <- find_cutpoint_number(
  data = analysis_data,
  predictor = "temperature",
  outcome_time = "time",
  outcome_event = "germinated",
  method = "genetic",        # Use the fast genetic algorithm
  criterion = "BIC",         # Use BIC for model selection
  use_parallel = TRUE,       # Enable parallel processing
  n_cores = 2,               # Specify number of cores (set to your preference)
  max_cuts = 5,              # Test models with 0, 1, 2, 3, 4, 5 cuts
  nmin = 0.1,                # Each group must contain at least 10% of data
  maxiter = 500,             # Iterations for the genetic algorithm
  seed = 42                  # Set seed for reproducible results
)

Interpretation of the BIC Table: This table above is our guide for choosing the right number of cut-points.

The model with the lowest BIC score is considered the best. Here, the model with 3 cut-points has the lowest BIC (10457.27).

Delta_BIC shows how much worse each other model is compared to the best one.

BIC_Weight tells us that the 3-cut-point model is overwhelmingly the most plausible, with a 95% probability of being the best model in this set.

Now, call summary() on this new result object for a full report.

summary(number_result_bic)
##  num_cuts      BIC Delta_BIC BIC_Weight    Evidence                             cuts
##         0 10878.79    421.52         0%     Minimal                               NA
##         1 10687.53    230.25         0%     Minimal                             7.52
##         2 10532.11     74.83         0%     Minimal                     10.71, 30.69
##         3 10457.27      0.00      88.2% Substantial                  9, 13.16, 30.62
##         4 10461.76      4.49       9.3%    Moderate            8.5, 13, 18.21, 30.38
##         5 10464.44      7.16       2.5%     Minimal 8.56, 14.34, 18.55, 25.14, 30.45
##   Group   N Events
## 1    G1 160    141
## 2    G2 160    145
## 3    G3 560    527
## 4    G4 160    102
## Call: survfit(formula = survival::Surv(time, event) ~ group, data = data)
## 
##            n events median 0.95LCL 0.95UCL
## group=G1 160    141     18      17      19
## group=G2 160    145     13      12      14
## group=G3 560    527     10      10      11
## group=G4 160    102     15      14      16
## Call:
## survival::coxph(formula = as.formula(formula_str), data = model_data)
## 
##   n= 1040, number of events= 915 
## 
##            coef exp(coef) se(coef)      z Pr(>|z|)    
## groupG2  1.3185    3.7378   0.1372  9.613  < 2e-16 ***
## groupG3  2.3464   10.4480   0.1306 17.965  < 2e-16 ***
## groupG4  0.8722    2.3922   0.1463  5.963 2.48e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## groupG2     3.738    0.26754     2.857     4.891
## groupG3    10.448    0.09571     8.088    13.496
## groupG4     2.392    0.41802     1.796     3.187
## 
## Concordance= 0.706  (se = 0.008 )
## Likelihood ratio test= 497.3  on 3 df,   p=<2e-16
## Wald test            = 393.6  on 3 df,   p=<2e-16
## Score (logrank) test = 475.4  on 3 df,   p=<2e-16

A plot makes the choice clear. The lowest point on the curve indicates the optimal number of cuts.

# The S3 plot method shows the BIC for each model.
# The lowest point (in orange) is the best.
plot(number_result_bic)
Information criterion (BIC) scores by number of cut-points. The lowest score indicates the optimal number.

Information criterion (BIC) scores by number of cut-points. The lowest score indicates the optimal number.

Interpretation: The BIC analysis strongly suggests that a model with 3 cut-points is optimal, as it has the lowest information criterion score. This provides statistical justification for creating four distinct temperature groups (e.g., Cool, Sub-Optimal, Optimal, Warm).

Step 2: Find the Optimal Temperature Thresholds

Now that we know we need three cut-points, we use find_cutpoint() to discover their specific values. We will optimise for the "logrank" statistic to find the thresholds that create the most statistically significant separation between the survival curves of the resulting groups.

Choosing an Optimisation criterion

Criterion What It Optimises Recommendation Rating (Accuracy & Performance)
"logrank" The statistical significance of the separation between survival curves (maximises the chi-squared statistic). The most common and standard method. Best when the primary goal is to prove a significant difference. Accuracy: ★★★★★ Performance: ★★★★★
"hazard_ratio" The effect size (maximises the Hazard Ratio). Best for clinical interpretability. Finds the cut-point that creates the largest practical difference between groups. Accuracy: ★★★★★ Performance: ★★★★☆
"p_value" The p-value from a Cox model (minimises the p-value). A powerful way to find the most significant split, but the p-value itself should be interpreted with caution due to multiple testing. Accuracy: ★★★★☆ Performance: ★★★★☆ 
cli_alert_info("Step 2: Finding optimal cut-point values...")

# This function searches for the 3 cut-points that best
# separate the data, according to the log-rank test.
multi_cut_result <- find_cutpoint(
  data = analysis_data,
  predictor = "temperature",
  outcome_time = "time",
  outcome_event = "germinated",
  method = "genetic",      # Use genetic algorithm
  criterion = "logrank",   # Optimise for the log-rank statistic
  num_cuts = 3,            # We are looking for 3 cuts (from Step 1)
  nmin = 0.1,              # Each group must have >= 10% of data
  maxiter = 500,           # Iterations for a precise search
  seed = 123               # Use a different seed
)

Now, let us call summary(multi_cut_result) on this result to see the comparison table.

summary(multi_cut_result)
##   Group   N Events
## 1    G1 160    141
## 2    G2 160    145
## 3    G3 560    527
## 4    G4 160    102
## Call: survfit(formula = survival::Surv(time, event) ~ group, data = data)
## 
##            n events median 0.95LCL 0.95UCL
## group=G1 160    141     18      17      19
## group=G2 160    145     13      12      14
## group=G3 560    527     10      10      11
## group=G4 160    102     15      14      16
## Call:
## survival::coxph(formula = as.formula(formula_str), data = data)
## 
##   n= 1040, number of events= 915 
## 
##            coef exp(coef) se(coef)      z Pr(>|z|)    
## groupG2  1.3185    3.7378   0.1372  9.613  < 2e-16 ***
## groupG3  2.3464   10.4480   0.1306 17.965  < 2e-16 ***
## groupG4  0.8722    2.3922   0.1463  5.963 2.48e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## groupG2     3.738    0.26754     2.857     4.891
## groupG3    10.448    0.09571     8.088    13.496
## groupG4     2.392    0.41802     1.796     3.187
## 
## Concordance= 0.706  (se = 0.008 )
## Likelihood ratio test= 497.3  on 3 df,   p=<2e-16
## Wald test            = 393.6  on 3 df,   p=<2e-16
## Score (logrank) test = 475.4  on 3 df,   p=<2e-16
##        chisq df    p
## group   1.79  3 0.62
## GLOBAL  1.79  3 0.62

A Note on Final Outputs The summary() output is optimised for the R console. To create a more polished final document, it’s better to build your own table. The following chunk shows how to pull the key results from the output object and use knitr::kable() to generate a clean, professional-looking table, giving you complete control over your report.

# Manually create a results data frame from the find_cutpoint object
results_summary_df <- data.frame(
  Parameter = c("Predictor", "Criterion", "Optimal Log-Rank Statistic", "Recommended Cut-point(s)"),
  Value = c(
    multi_cut_result$parameters$predictor,
    multi_cut_result$parameters$criterion,
    round(multi_cut_result$optimal_stat, 4),
    paste(round(multi_cut_result$optimal_cuts, 3), collapse = ', ')
  )
)

# Use kable() to create a beautiful and robust table
knitr::kable(
  results_summary_df,
  col.names = c("Metric", "Result"), # Nicer column names
  caption = "Summary of Optimal Temperature Threshold Analysis"
)
Summary of Optimal Temperature Threshold Analysis
Metric Result
Predictor temperature
Criterion logrank
Optimal Log-Rank Statistic 454.4678
Recommended Cut-point(s) 8.703, 12.804, 31.826

A plot of the temperature distribution with the discovered cut-points overlaid gives us a clear visual confirmation.

# The S3 plot method (type = "distribution") shows a histogram
# of the predictor with the optimal cut-points as vertical lines.
plot(multi_cut_result, 
     type = "distribution", 
     palette = "jco", 
     title = "Optimized Biomarker Stratification")
Distribution of the predictor variable with optimal cut-points overlaid.

Distribution of the predictor variable with optimal cut-points overlaid.

Interpretation: The algorithm has identified three temperature thresholds: 8.7°C, 12.8°C, and 31.8°C. These values are not arbitrary; they represent the points that best separate the rapeseed seeds into four groups with distinct germination behaviours.

Step 3: Validate the Cut-point Stability

A key question is whether these cut-points are stable or just an artefact of our specific sample. We run a bootstrap analysis with validate_cutpoint() to generate 95% confidence intervals.

cli_alert_info("Step 3: Validating cut-point stability with bootstrapping...")

# This function re-runs the find_cutpoint algorithm many times on
# resampled data to see how much the cut-points vary.
validation_result <- validate_cutpoint(
  cutpoint_result = multi_cut_result, # Pass in the results from Step 2
  num_replicates = 100, # Use a lower number for a fast demo (>= 500 for publication)
  use_parallel = TRUE,  # Enable parallel processing
  n_cores = 2,          # Specify number of cores (set to your preference)
  maxiter = 100,        # Use fewer iterations *within each replicate* to speed up the demo
  seed = 456            # Use a different seed
)

summary(validation_result)
## Cut-point Stability Analysis (Bootstrap)
## ----------------------------------------
## Original Optimal Cut-point(s): 8.703, 12.804, 31.826 
## 
## Bootstrap Distribution Summary
## -----------------------------
##       Cut   Mean    SD Median     Q1     Q3
## 25%  Cut1  8.994 0.886  8.907  8.274  9.627
## 25%1 Cut2 14.707 1.472 14.583 13.858 15.168
## 25%2 Cut3 30.986 1.096 31.027 30.209 31.956
## 
## 95% Confidence Intervals
## ------------------------
##        Lower  Upper
## Cut 1  7.617 11.136
## Cut 2 12.696 17.145
## Cut 3 28.405 32.462
## 
## Validation Parameters
## ---------------------
## Replicates Requested: 100 
## Successful Replicates: 100 / 100 ( 100 %)
## Failed Replicates: 0 
## Cores Used: 2 
## Seed: 456 
## Minimum Group Size (nmin): 93 
## Method: genetic 
## Criterion: logrank 
## Covaricates: None

A density plot of the bootstrap results provides a powerful visual for assessing stability. Narrow, sharp peaks indicate a highly stable cut-point.

# The S3 plot method for a validation object shows the density
# plots of the cut-points found during bootstrapping.
plot(validation_result)
Bootstrap distribution of the three optimal cut-points. The solid red line represents the original cut-point.

Bootstrap distribution of the three optimal cut-points. The solid red line represents the original cut-point.

Interpretation: The bootstrap results give us confidence in our findings. The confidence intervals for the cut-points are relatively narrow, especially for the lower and middle thresholds. This suggests that these are robust findings and not just random chance.

4. Advanced Visualisation & Reporting

OptSurvCutR provides a unified, highly extensible S3 plotting engine. Rather than forcing you to extract data, manually format groups, and fit complex Cox models from scratch, you can simply call plot() on your result object and use the type argument to generate a suite of publication-ready clinical graphics.

By default, the package mathematically assigns your data into groups labelled G1, G2, G3, etc., ordered sequentially from the lowest to highest predictor value.

4.1 The “Big Picture” Dashboard (type = "all")

The most powerful way to interpret a cut-point is to see exactly where the cuts were made on the underlying data, stacked directly above the resulting survival curves. This composite “cause and effect” dashboard is the perfect starting point for new analyses.

# Generate the composite cause-and-effect dashboard
plot(multi_cut_result, type = "all")

Interpretation: The top graph (Predictor Distribution) shows the natural distribution of your temperatures with dashed lines indicating exactly where the algorithm split the data. The bottom graph (Kaplan-Meier) instantly shows how those splits impacted germination rates.

4.2 Verifying Group Composition (The Escape Hatch)

Before diving into deeper diagnostics, you may want to verify exactly which temperatures fell into G1, G2, G3, and G4. You can use return_data = TRUE to extract the raw data with the groups perfectly assigned, acting as an “escape hatch” for custom tables.

library(dplyr)

# Extract the raw data frame with the optimal groups natively assigned
final_dataset <- plot(multi_cut_result, return_data = TRUE)

# Example logic for viewing which temperatures map to which group:
composition_table <- final_dataset %>%
  group_by(group) %>%
  summarise(
    Temperatures_in_Group = paste(sort(unique(factor)), collapse = ", ")
  )

knitr::kable(composition_table, caption = "Composition of Discovered Temperature Groups")
Composition of Discovered Temperature Groups
group Temperatures_in_Group
G1 5, 7.5
G2 10, 12.5
G3 15, 17.5, 20, 22.5, 25, 27.5, 30
G4 32.5, 35

Based on this table, we can interpret the biological meaning of our groups: * G1: Cool * G2: Sub-Optimal * G3: Optimal * G4: Warm

4.3 Survival Outcomes (type = "outcome")

To generate a standalone Kaplan-Meier curve, use the outcome type. Because OptSurvCutR passes additional arguments directly to the survminer engine, you can customize titles, palettes, and labels in a single line of code.

# Generate a standalone Kaplan-Meier plot with custom aesthetics
plot(multi_cut_result, type = "outcome",
     title = "Germination by Optimal Temperature Strata",  
     xlab = "Time (Days)", ylab = "Proportion Ungerminated",  
     legend.title = "Temp Group")
Kaplan-Meier survival curve for germination by temperature group.

Kaplan-Meier survival curve for germination by temperature group.

Interpretation: The “Optimal” group (G3) curve drops most steeply, indicating the fastest germination. The “Cool” (G1) and “Warm” (G4) groups show much flatter curves, confirming inhibited germination. The global p-value confirms these differences are highly statistically significant.

4.4 Hazard Ratios & Forest Plots (type = "forest")

The forest plot visualizes the relative risk (Hazard Ratios) between your groups. By default, standard Cox models use G1 as the baseline. However, OptSurvCutR allows you to instantly change the baseline using the reference_group argument—perfect for comparing everything against our “Optimal” G3 group.

# Generate a Forest Plot comparing all groups against G3 (Optimal)
plot(multi_cut_result, type = "forest", reference_group = "G3",
     main = "Hazard Ratios Relative to Optimal Temperature (G3)")
Forest plot of Hazard Ratios for each temperature group relative to the optimal group.

Forest plot of Hazard Ratios for each temperature group relative to the optimal group.

Interpretation: This plot shows the rate of germination for each group compared to the optimal group (G3). An HR less than 1 (e.g., G1 and G4) means a significantly lower rate of germination compared to the baseline. Because none of the confidence interval lines cross the 1.0 vertical threshold, we can conclude that every other temperature zone is significantly worse than the optimal zone.

4.5 Advanced Model Diagnostics (type = "diagnostic")

For major publications, reviewers often require proof that your Cox models do not violate the Proportional Hazards assumption. OptSurvCutR automates this complex check natively by calculating and plotting the Schoenfeld residuals.

# Automatically run survival::cox.zph() and plot the residuals
plot(multi_cut_result, type = "diagnostic")
## $`1`
Schoenfeld residuals diagnostic plot to check Proportional Hazards assumptions.

Schoenfeld residuals diagnostic plot to check Proportional Hazards assumptions.

Interpretation: For the Cox model to be statistically valid, the hazard ratios between the groups must remain roughly constant over time. If the trend lines in this plot are relatively flat (and the global p-value is > 0.05), the Proportional Hazards assumption is met, giving you confidence in your Forest Plot results.

4.6 Predictive Power Over Time (type = "auc")

Does the predictive power of our cut-points hold up over the entire timeline of the experiment? We can evaluate this by calculating the Time-Dependent Area Under the Curve (AUC).

# Requires the 'timeROC' package
plot(multi_cut_result, type = "auc")
Time-dependent AUC demonstrating the predictive power of the cut-points across the duration of the trial.

Time-dependent AUC demonstrating the predictive power of the cut-points across the duration of the trial.

Interpretation: This line graph tracks model accuracy over time. An AUC of 0.5 is no better than random guessing, while 1.0 is perfect prediction. A high, stable curve indicates that the temperature cut-points remain a strong predictor of germination throughout the entire duration of the trial.

4.7 Interactive Web Mode

If you are rendering HTML documents, RMarkdown reports, or Vignettes, static plots can be limiting. You can wrap any OptSurvCutR graph in optsurv_interactive() to instantly convert it into a dynamic web widget.

# 1. Generate a standard distribution plot
p_dist <- plot(multi_cut_result, type = "distribution")

# 2. Make it interactive! (Hover your mouse over the distribution peaks)
optsurv_interactive(p_dist)

4.8 The “Escape Hatch” (Raw Data Extraction)

While the built-in graphs are highly capable, sometimes you need the exact grouped data to build a custom visualization in another platform or run a bespoke statistical test. By using return_data = TRUE, the package acts as a data-prepper, returning your original dataset with a new, perfectly assigned group column.

# Extract the raw data frame with the optimal groups assigned
final_dataset <- plot(multi_cut_result, return_data = TRUE)

# View the assigned groups
head(final_dataset[, c("factor", "time", "event", "group")])
##   factor time event group
## 1      5   22     1    G1
## 2      5   18     0    G1
## 3      5   15     1    G1
## 4      5   22     1    G1
## 5      5   19     1    G1
## 6      5   14     1    G1

5. Conclusion & Next Steps

This vignette has demonstrated the three-step workflow for cut-point analysis using OptSurvCutR. By following this workflow, users can confidently identify and validate robust, statistically-optimal thresholds in their own survival data, moving beyond simple median splits to uncover more nuanced relationships.

We encourage you to try OptSurvCutR with your own data.

  • Install the package from GitHub: remotes::install_github("paytonyau/OptSurvCutR")
  • Report issues or suggest features on our GitHub page.
  • Star the repository if you find it useful.
  • Read the preprint: For a detailed methodological discussion, please see our paper on bioRxiv.

6. Session Information

For reproducibility, the session information below lists the R version and all attached packages used to run this analysis.

sessionInfo()
## R version 4.5.3 (2026-03-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8  LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8 LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/London
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rgenoud_5.9-0.11  cli_3.6.5         tidyr_1.3.1       knitr_1.50        patchwork_1.3.2   dplyr_1.1.4      
##  [7] survminer_0.5.1   ggpubr_0.6.2      ggplot2_4.0.1     survival_3.8-6    OptSurvCutR_0.2.1 pacman_0.5.1     
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1    viridisLite_0.4.2   farver_2.1.2        S7_0.2.1            lazyeval_0.2.2     
##  [6] fastmap_1.2.0       digest_0.6.39       lifecycle_1.0.5     magrittr_2.0.4      compiler_4.5.3     
## [11] rlang_1.1.7         sass_0.4.10         rngtools_1.5.2      tools_4.5.3         yaml_2.3.10        
## [16] data.table_1.17.8   ggsignif_0.6.4      htmlwidgets_1.6.4   timereg_2.0.7       labeling_0.4.3     
## [21] doRNG_1.8.6.2       xml2_1.5.0          RColorBrewer_1.1-3  abind_1.4-8         numDeriv_2016.8-1.1
## [26] withr_3.0.2         purrr_1.2.0         pec_2025.06.24      grid_4.5.3          xtable_1.8-4       
## [31] future_1.68.0       globals_0.18.0      scales_1.4.0        iterators_1.0.14    mvtnorm_1.3-3      
## [36] rmarkdown_2.30      generics_0.1.4      rstudioapi_0.17.1   future.apply_1.20.0 km.ci_0.5-6        
## [41] httr_1.4.7          cachem_1.1.0        stringr_1.6.0       splines_4.5.3       parallel_4.5.3     
## [46] survMisc_0.5.6      vctrs_0.6.5         Matrix_1.7-4        jsonlite_2.0.0      carData_3.0-5      
## [51] car_3.1-3           rstatix_0.7.3       Formula_1.2-5       listenv_0.10.0      crosstalk_1.2.2    
## [56] foreach_1.5.2       plotly_4.12.0       jquerylib_0.1.4     parallelly_1.45.1   glue_1.8.0         
## [61] codetools_0.2-20    ggtext_0.1.2        cowplot_1.2.0       stringi_1.8.7       gtable_0.3.6       
## [66] tibble_3.3.0        pillar_1.11.1       htmltools_0.5.8.1   lava_1.8.2          R6_2.6.1           
## [71] KMsurv_0.1-6        doParallel_1.0.17   evaluate_1.0.5      lattice_0.22-9      backports_1.5.0    
## [76] gridtext_0.1.5      broom_1.0.10        ggsci_4.1.0         bslib_0.9.0         timeROC_0.4.1      
## [81] Rcpp_1.1.0          gridExtra_2.3       prodlim_2026.03.11  xfun_0.54           zoo_1.8-14         
## [86] pkgconfig_2.0.3