---
title: "Analysing Clinical Survival Data with OptSurvCutR"
subtitle: "Discovering biomarker thresholds and tier systems in liver disease"
author:
  - "Payton Yau"
  - "Suhirthakumar Puvanendran"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    fig_caption: yes
vignette: >
  %\VignetteIndexEntry{Analysing Clinical Survival Data with OptSurvCutR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, echo = TRUE, comment = "+",
  fig.width = 8, fig.height = 4.5, fig.align = "center",
  warning = TRUE, message = TRUE
)
options(cli.num_colors = 1)
# Suppress cli progress bars (validate_cutpoint()'s "Bootstrapping" bar):
# a live progress bar re-renders on every step during a vignette knit,
# producing hundreds of lines of output in the rendered document rather
# than a live terminal display. Interim measure until every
# validate_cutpoint() call below is updated to pass quiet = TRUE directly.
options(cli.progress_show_after = Inf)
```

# Introduction

Continuous biomarkers are routinely split into risk groups using a median
or a threshold borrowed from another study. Both are arbitrary, and a
single split cannot represent a relationship where risk rises in stages.

`OptSurvCutR` selects thresholds by optimising a survival criterion, and
then assesses how far those thresholds move under resampling.

| Function | Question | Output |
|:---|:---|:---|
| `find_cutpoint_number()` | How many groups? | A number, *k* |
| `find_cutpoint()` | Where are the boundaries? | *k* thresholds |
| `validate_cutpoint()` | How much do they move? | Intervals + stability tier |

## How this differs from a single optimal cut-point

Established tools solve part of this problem well. `maxstat` computes a
maximally selected rank statistic with a corrected *p*-value, which
addresses the inflation that comes from testing many candidate
thresholds; `survminer::surv_cutpoint()` wraps it in a convenient
interface. If a single unadjusted threshold is what the question needs,
they are a sound choice.

`OptSurvCutR` covers three things they do not:

| | `maxstat` / `survminer` | `OptSurvCutR` |
|:---|:---|:---|
| Number of cut-points | One | Chosen from the data, one or more |
| Covariate adjustment | No | Yes, inside the search |
| Threshold stability | Not assessed | Bootstrap intervals and a tier |

The third is the one that changes conclusions. A corrected *p*-value
tells you a threshold is unlikely to have arisen by chance; it says
nothing about whether the same threshold would be found again in a
comparable sample. Those are different questions, and in this vignette
they give different answers: the four-group model is significant by any
test, yet its boundaries are not reproducible.

## The data

The Mayo Clinic trial of D-penicillamine in **primary biliary
cholangitis** (Dickson et al., 1989) followed 418 patients with this
chronic autoimmune liver disease, in which progressive destruction of the
bile ducts leads to cholestasis, cirrhosis and eventually liver failure.

**Serum bilirubin** is the classic marker of that progression: as bile
drainage fails, bilirubin accumulates. It is the strongest single
component of the Mayo risk score, and clinicians have long used
thresholds of it to time referral for transplantation. Values above
roughly 1.2 mg/dL are considered abnormal.

Here we ask what threshold the data themselves support, adjusting for
age, sex and the presence of edema.

------------------------------------------------------------------------

# 1. Data preparation

```{r load-libraries, message=FALSE, warning=FALSE}
library(survival)
library(dplyr)
library(ggplot2)
library(knitr)
library(OptSurvCutR)
```

```{r load-data}
data(pbc, package = "survival")

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

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

nrow(pbc_clean)
```

All six variables are complete across `pbc`, so the analysis retains all
418 patients. The figure of 312 often quoted for this dataset refers to
records with complete `trt` and laboratory values, which this model does
not use.

Transplant is treated as an event rather than a censoring. Transplants
are allocated preferentially to the sickest patients, so censoring them
would remove the highest-risk patients from the high-bilirubin groups
just as they were about to fail.

Adjusting for `age`, `sex` and `edema` means thresholds are selected
inside a Cox model containing those variables, so a threshold cannot
appear useful merely by tracking one of them.

------------------------------------------------------------------------

# 2. How many cut-points?

`find_cutpoint_number()` compares models of increasing complexity using
an information criterion. AIC penalises complexity lightly and suits
exploratory work; BIC is stricter and is the better choice when a
threshold is intended for clinical use. AICc corrects AIC for small
samples ($N < 200$).

```{r find-number}
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,          # each group holds at least 15% of patients
  max.generations = NULL, pop.size = NULL,
  boundary.enforcement = 2, seed = 123
)

summary(num_res)
```

```{r criterion-curve}
plot(num_res)
```

AIC is minimised at **three cut-points**, carrying 98.4% of the AIC
weight against 1.6% for the two-cut model.

The four- and five-cut models return no valid solution: with
`nmin = 0.15` each group needs 62 patients, and the genetic search found
no partition that satisfied this. A search returning nothing does not
prove that nothing exists — raising `max.generations` and `pop.size`, or
lowering `nmin`, may find one.

------------------------------------------------------------------------

# 3. Where are the boundaries?

```{r find-cutpoint}
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 from step 1
  nmin = 0.15,
  n_perm = 20,          # low for build speed; use >= 1000 when reporting
  max.generations = NULL, pop.size = NULL,
  boundary.enforcement = 2, seed = 123, n_cores = 1
)

summary(cut_res)
```

`num_cuts` is taken from the Step 1 object rather than specified
directly. Selecting the number of groups after inspecting the survival
curves would reintroduce the selection problem that Step 1 is intended to
control.

The thresholds are **0.7, 2.3 and 5.945 mg/dL**, giving groups of 100,
175, 80 and 63 patients. Hazard ratios rise across them (3.27, 12.20,
20.93 relative to the lowest group) and median transplant-free survival
falls from not reached to 3445, 1504 and 930 days. Concordance is 0.811.

The reported permutation *p* of 0.0476 is exactly 1/(20+1), the smallest
value `n_perm = 20` can return. It means no permutation exceeded the
observed statistic, not that *p* equals 0.0476.

The Schoenfeld test gives *p* = 0.103, so there is no evidence against
proportional hazards and the hazard ratios can be read as approximately
constant over follow-up.

```{r plot-distribution}
plot(cut_res, type = "distribution") +
  geom_rug(alpha = 0.5) +
  labs(caption = "Bilirubin thresholds on the marker distribution")
```

------------------------------------------------------------------------

# 4. Are the boundaries stable?

`validate_cutpoint()` re-runs the search on resampled cohorts and records
where each boundary lands.

```{r validate-cutpoint}
val_res <- validate_cutpoint(
  cutpoint_result = cut_res,
  num_replicates = 30,      # reduced for build speed; use >= 500 when reporting
  n_cores = 1,      # single core: vignette builds cannot use parallel workers
  max.generations = NULL, pop.size = NULL,
  boundary.enforcement = 2, seed = 123
)

summary(val_res)
```

The bootstrap `nmin` is relaxed to 55 automatically. Resampled cohorts
contain duplicates and fewer distinct values, so holding the original
constraint would cause replicates to fail.

## Reading the tier

Two quantities are assessed: whether adjacent intervals **overlap**, and
each interval's **width** relative to the predictor's 10th–90th
percentile range.

| Tier | Overlap | Width | Meaning |
|:---|:---|:---|:---|
| 1 — OPTIMAL | None | < 30% | Precise boundaries, distinct groups |
| 2 — DISTINCT | None | Any | Groups distinct, boundaries move |
| 3 — CAUTION | Present | 30–60% | Adjacent groups not cleanly separated |
| 4 — UNSTABLE | Present | > 60% | Boundaries not reproducible |

These are two independent diagnostics rather than an ordinal scale:
overlap matters most for a multi-tier rule, width for a single reported
boundary. The percentages are practical conventions.

The number of replicates affects the width of the intervals. With few
replicates the 2.5th and 97.5th percentiles fall close to the extremes of
the resampled values, and the intervals are narrower than they should be.
Use at least 500 replicates for any reported analysis. The chunk above
uses 30 to keep the vignette build short; the section below reports the
500-replicate results.

------------------------------------------------------------------------

# 5. Results at 500 replicates

Running the same validation with `num_replicates = 500` gives:

| Model | Thresholds | Widest interval | Tier |
|:---|:---|---:|:---|
| Three cut-points | 0.7, 2.3, 5.945 | 54.2% | 3 — CAUTION |
| Two cut-points | 2.3, 5.945 | 59.9% | 3 — CAUTION |
| **One cut-point** | **2.3** | **18.8%** | **1 — OPTIMAL** |

For the three-cut model the intervals are 0.6–1.6, 1.4–3.3 and
3.0–7.03 mg/dL. Both adjacent pairs intersect, so a patient with a
bilirubin of 1.5 mg/dL is assigned to a different group depending on the
resample. The component widths are uneven: the two lower boundaries are
well resolved at 13.5% and 25.6%, both within Tier 1 bounds on width
alone, whereas the upper boundary, situated in the sparse right tail, is
not reproducible at 54.2%.

The two-cut model performs less well, with a widest interval of 59.9% and
the lowest boundary widening from 0.6–1.6 to 0.7–3.0 mg/dL. With three
groups the minimum stratum constraint admits a wider feasible region, so
the optimum varies more across resamples. Reducing the number of
cut-points does not necessarily improve stability, and is worth testing
rather than assuming.

The threshold is **2.3 mg/dL** with a bootstrap interval of 2.0–3.4. The
median across 500 resamples is 2.3, identical to the value found in the
full data, and the interquartile range is 2.3–3.0. This is also the
middle boundary of the three-cut model: the search returned to the same
point regardless of how many groups it was asked for.

The reduction from three cut-points to one followed the stability
assessment and constitutes a post-hoc simplification, which should be
reported as such. Two considerations support it: the sequence is
prescribed by the workflow in advance rather than selected for this
dataset, and `validate_cutpoint()` issues the recommendation as part of
its diagnostic output. The adjusted hazard ratio for the
single-threshold model nonetheless remains optimistic, since the
threshold was selected from the same data used to estimate it.

------------------------------------------------------------------------

# 6. Reporting

## Group composition

```{r group-composition}
final_dataset <- plot(cut_res, return_data = TRUE)

final_dataset %>%
  group_by(group) %>%
  summarise(
    Bilirubin = paste0(round(min(factor), 2), " – ", round(max(factor), 2)),
    N = n(),
    Events = sum(event),
    Mean_Age = round(mean(age), 1),
    Edema = round(mean(edema), 2)
  ) %>%
  kable(caption = "Composition of the four-group model")
```

Mean age is flat across groups, so the survival gradient is not age in
disguise. Edema is seven times more common in the highest group than the
lowest and is strongly associated with the outcome, which is why
it is included as an adjustment variable. The highest group contains 63
patients against a floor of 62, so its lower boundary is determined
partly by the `nmin` constraint rather than by the data alone.

## Hazard ratios and diagnostics

```{r plot-forest, fig.width=7, fig.height=3.5}
plot(cut_res, type = "forest",
     main = "Adjusted hazard ratios relative to group 1")
```

These estimates come from the data that selected the thresholds and are
therefore optimistic.

```{r plot-diagnostic}
plot(cut_res, type = "diagnostic")
```

```{r plot-km}
plot(cut_res, type = "outcome",
     title = "Transplant-free survival by bilirubin group",
     xlab = "Follow-up (days)", ylab = "Transplant-free survival",
     legend.title = "Bilirubin group")
```

## Joint stability

An interval shows how one threshold moves; it cannot show whether two
move together.

```{r plot-validation-2d-12}
plot_validation(val_res, focus_cuts = c(1, 2),
                main = "Cuts 1 and 2 across resamples")
```

```{r plot-validation-2d-23}
plot_validation(val_res, focus_cuts = c(2, 3),
                main = "Cuts 2 and 3 across resamples")
```
A tight cloud indicates both boundaries are well determined. Elongation
along one axis identifies the imprecise boundary. A diagonal spread means
the two are trading off, with several partitions scoring similarly.

------------------------------------------------------------------------

# 7. Conclusion

AIC selected four bilirubin groups with 98.4% of the model weight,
monotonic hazard ratios and concordance of 0.81. The bootstrap showed
that the boundaries overlapped and would not survive a different sample
of the same size. A single threshold at **2.3 mg/dL** — close to twice
the upper limit of normal — proved highly reproducible.

Significance and stability are different properties. A tool reporting
only thresholds and a *p*-value would have returned four groups here, and
every check it ran would have supported them. The third step is one
function call, and in this analysis it is the difference between
reporting "0.7, 2.3 and 5.945 mg/dL" and reporting "2.3 mg/dL".

When reporting an analysis of this kind: use BIC for confirmatory work,
`n_perm` of at least 1000, at least 500 bootstrap replicates, and check
the composition table for covariate imbalance. Bootstrap stability
describes sampling variability within one cohort; it is not a prediction
about an independent one.

## Reference

Dickson ER, Grambsch PM, Fleming TR, Fisher LD, Langworthy A (1989).
Prognosis in primary biliary cirrhosis: model for decision making.
*Hepatology* 10:1–7.

```{r session_info}
sessionInfo()
```