Survival Analysis

Survival Analysis

Time-to-event analysis for cancer clinical data. Covers Kaplan-Meier, Cox proportional hazards, competing risks, restricted mean survival time, and optimal cutpoint selection using the survival, ggsurvfit, tidycmprsk, and survRM2 packages.

When to Use This Skill

Activate when the user requests:

• Kaplan-Meier survival curves with risk tables
• Log-rank or stratified log-rank tests between groups
• Cox proportional hazards modeling (univariable or multivariable)
• Proportional hazards assumption checking
• Competing risks analysis (cause-specific or Fine-Gray)
• Restricted mean survival time (RMST) comparisons
• Optimal biomarker cutpoint selection for survival
• Forest plots for multivariate Cox models

Inputs

| Data Type | Format | Source | |-----------|--------|--------| | Clinical | Tabular (time, event, covariates) | TCGA GDC via TCGAbiolinks, cBioPortal | | Time variable | Days/months to event or last follow-up | days_to_death, days_to_last_follow_up | | Event indicator | Binary (0 = censored, 1 = event) | Derived from vital_status | | Covariates | Categorical or continuous | Age, stage, gene expression, mutations |

---

Preparing TCGA Survival Data

library(TCGAbiolinks)
library(SummarizedExperiment)

query  0, ]

Kaplan-Meier Estimation

library(survival)    # v3.8+
library(ggsurvfit)   # v1.2+, pharmaverse — modern replacement for survminer

# survfit2() tracks the calling environment for cleaner legend labels
# Use survfit2() with ggsurvfit; use survfit() with base R plotting
km_fit  log-rank is most powerful
  Curves cross or separate late (common in immunotherapy)
    -> log-rank has low power. Use RMST or weighted log-rank.
  > 2 groups
    -> log-rank gives a global test. Follow up with pairwise comparisons
       using p.adjust() for multiple testing correction.

Cox Proportional Hazards

Univariable

# Continuous covariate
cox_age = 10 events per predictor in the model.
  With 50 events, fit at most 5 covariates.
  More recent work (Riley 2019) recommends >= 20 EPV for adequate precision.
  Too many covariates relative to events = overfitting and unstable HRs.

Continuous vs categorical:
  Keep continuous variables continuous in Cox models. Dichotomizing age at
  median discards information and reduces power.
  If a cutpoint is needed for clinical communication, derive it separately
  (see Optimal Cutpoints below) and validate externally.

Proportional Hazards Assumption

# Schoenfeld residual test
zph  0.05: PH assumption not rejected (but doesn't prove PH holds)
# p  PH holds
# Systematic trend (slope) -> PH violated, effect changes over time
plot(zph)
# Or with survminer: survminer::ggcoxzph(zph)

Interpreting Schoenfeld residuals:
  Positive slope: effect increases over time (HR drifts upward)
  Negative slope: effect weakens over time (early benefit that fades)
  U-shape or crossing zero: complex time-varying effect

  The p-value is sample-size dependent: large datasets reject PH
  for clinically negligible departures. Always look at the plot.

Handling PH Violations

# Option 1: Time-varying coefficient with tt()
# Preferred when a specific covariate violates PH
cox_tv  HR is more efficient. Report both HR and RMST.
  Curves cross or separate late (immunotherapy, delayed effect)?
    -> RMST captures the difference that HR misses.
  Need a clinically interpretable number?
    -> RMST: "Treatment gives 3.2 extra months of survival on average"
       HR: "Treatment reduces the hazard by 30%" (less intuitive)

Optimal Cutpoint Selection

library(survminer)  # v0.5+, provides surv_cutpoint

# maxstat approach: maximally selected log-rank statistic
# Finds the cutpoint that maximally separates high/low groups
cp  14%, PSA > 4 ng/mL)
     - Communication to non-statistical audiences

  Use continuous modeling for primary analysis, cutpoints for clinical translation.

Forest Plots

library(forestmodel)  # v0.6+

# Directly from a coxph object — one function call
forest_model(cox_multi)
# Produces a ggplot2-based forest plot with HRs, 95% CIs, p-values

# Customize
forest_model(cox_multi,
  factor_separate_line = TRUE,       # categorical levels on separate lines
  format_options = forest_model_format_options(
    text_size = 3.5, point_size = 3
  )
)

# Alternative: survminer::ggforest() for quick plots
# Limitation: does not handle interactions or spline terms properly
survminer::ggforest(cox_multi, data = clinical)

Output Specification

| Output | Format | Description | |--------|--------|-------------| | km_curves.pdf | PDF | Kaplan-Meier curves with risk table and p-value | | cox_summary.csv | CSV | HR, 95% CI, p-value for each covariate | | cox_zph.pdf | PDF | Schoenfeld residual plots for PH assumption | | forest_plot.pdf | PDF | Multivariate Cox forest plot | | cif_curves.pdf | PDF | Cumulative incidence curves for competing risks | | rmst_results.csv | CSV | RMST difference, ratio, CI, p-value | | cutpoint_analysis.csv | CSV | Optimal cutpoint, corrected p-value, HR at cutpoint |

Validation Checks

After running survival analysis on TCGA-GBM, verify:

Overall survival:
  Median OS for all GBM: ~14-16 months
  If median OS > 24 months: likely includes IDH-mutant cases (not true GBM under WHO 2021)
  2-year survival: ~20-25%

IDH1 mutation (strongest single prognostic factor in glioma):
  IDH-mutant median OS: ~31 months
  IDH-wildtype median OS: ~14-16 months
  Log-rank p-value: should be  80)

Proportional hazards:
  Age typically satisfies PH in GBM (chronic risk factor)
  Treatment effect may violate PH (delayed separation with TMZ)
  If global cox.zph p = 10 rule (preferably >= 20).
5. **Stepwise selection for final models**: Forward/backward stepwise selection inflates Type I error and produces non-reproducible models. Pre-specify covariates based on clinical knowledge.
6. **Ignoring collinearity**: Stage and tumor size are correlated. Including both in a Cox model inflates standard errors and produces unstable HRs. Check VIF or pairwise correlations before fitting.

### Competing Risks
7. **Standard KM for competing risks endpoints**: KM treats competing events as censored, which inflates cumulative incidence of the event of interest. Use cumulative incidence functions (CIF) via `tidycmprsk::cuminc()` when competing events are present.
8. **Interpreting Fine-Gray as causal**: Fine-Gray SDH measures the effect on cumulative incidence, not on the rate among those at risk. A covariate that increases competing event risk will artificially appear to decrease the subdistribution hazard for the primary event.

### RMST and Cutpoints
9. **Post-hoc tau selection**: Choosing tau after seeing the data invalidates the RMST test. Pre-specify tau based on clinical context (standard follow-up duration, planned analysis timepoint).
10. **Cutpoint without validation**: An optimal cutpoint found by maxstat will appear highly significant in the discovery cohort even for a non-prognostic biomarker. Without independent validation, the cutpoint is a hypothesis, not a result.

## Related Skills

- [`cancer-multiomics`](../cancer-multiomics/SKILL.md): TCGA data retrieval, mutation status and gene expression for Cox covariates
- [`immune-deconvolution`](../immune-deconvolution/SKILL.md): Immune cell fraction estimates as continuous survival predictors

## Public Datasets for Testing

| Dataset | Samples | Use Case |
|---------|---------|----------|
| TCGA-GBM | 611 | Glioblastoma, strong IDH1/MGMT prognostic markers, poor prognosis |
| TCGA-BRCA | 1098 | Breast cancer, molecular subtype-specific survival differences |
| TCGA-LUAD | 585 | Lung adenocarcinoma, EGFR/KRAS mutation-driven survival |
| TCGA-OV | 608 | Ovarian cancer, classic for competing risks (progression vs death) |

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [zamushwani2](https://github.com/zamushwani2)
- **Source:** [zamushwani2/biomedical-ai-skills](https://github.com/zamushwani2/biomedical-ai-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.