Scikit Survival Analysis

pip install scikit-survival scikit-learn pandas numpy matplotlib

decisions:
  - id: D1
    param: timeAndEventColumns
    kind: required
    source: data
    ask: "Which column holds follow-up time, which holds the event indicator, and does the indicator mark the event or censoring?"
    default: null

  - id: D2
    param: eventDefinition
    kind: required
    source: user
    depends_on: [D1]
    ask: "Which event is being modelled, and are competing events treated as censoring?"
    default: null

  - id: D3
    param: covariates
    kind: required
    source: data
    ask: "Which variables enter the model, and are any measured after baseline?"
    default: null

  - id: D4
    param: modelFamily
    kind: required
    source: user
    ask: "A proportional-hazards model whose coefficients are interpretable, or an ensemble that predicts better but does not yield hazard ratios?"
    default: "Cox proportional hazards"

  - id: D5
    param: regularization
    kind: required
    source: user
    depends_on: [D3, D4]
    ask: "With more covariates than events, should coefficients be penalized - and toward selection or toward shrinkage?"
    default: "none"

  - id: D6
    param: validationScheme
    kind: required
    source: user
    ask: "How should performance be estimated - cross-validation, a held-out split, or an external cohort?"
    default: "cross-validation"

  - id: D7
    param: ensembleHyperparameters
    kind: optional_conditional
    source: user
    depends_on: [D4]
    ask: "How many trees, how deep, and at what learning rate?"
    default: "100 estimators, depth 3, learning rate 0.1"
    skip_if: "a proportional-hazards model was chosen"

from sksurv.datasets import load_breast_cancer
from sksurv.ensemble import RandomSurvivalForest
from sksurv.metrics import concordance_index_ipcw
from sklearn.model_selection import train_test_split

X, y = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

rsf = RandomSurvivalForest(n_estimators=100, random_state=42)
rsf.fit(X_train, y_train)

risk_scores = rsf.predict(X_test)
c_index = concordance_index_ipcw(y_train, y_test, risk_scores)[0]
print(f"C-index: {c_index:.3f}")  # e.g., 0.68

# Individual survival curves
surv_fns = rsf.predict_survival_function(X_test[:2])
for fn in surv_fns:
    print(f"5-year survival: {fn(365 * 5):.3f}")

import numpy as np
import pandas as pd
from sksurv.util import Surv
from sksurv.preprocessing import OneHotEncoder, encode_categorical
from sksurv.datasets import load_gbsg2, load_breast_cancer
from sklearn.preprocessing import StandardScaler

# Create survival outcome from arrays
event = np.array([True, False, True, True, False])
time = np.array([120.0, 365.0, 200.0, 90.0, 400.0])
y = Surv.from_arrays(event=event, time=time)
print(y.dtype)  # [('event', '?'), ('time', '<f8')]

# From DataFrame columns
# y = Surv.from_dataframe("event_col", "time_col", df)

# Load built-in datasets
# Available: load_gbsg2, load_breast_cancer, load_veterans_lung_cancer,
#            load_whas500, load_aids, load_flchain
X, y = load_gbsg2()
print(f"Shape: {X.shape}, Events: {y['event'].sum()}, "
      f"Censoring rate: {1 - y['event'].mean():.1%}")

# Encode categoricals (survival-aware one-hot)
X_encoded = encode_categorical(X)  # auto-detect and encode all categorical cols

# Standardize (critical for Cox and SVM models)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_encoded)

from sksurv.io import loadarff

# Load ARFF format (Weka format)
data = loadarff("survival_data.arff")
X_arff, y_arff = data[0], data[1]  # DataFrame, structured array

from sksurv.linear_model import CoxPHSurvivalAnalysis, CoxnetSurvivalAnalysis, IPCRidge

# Standard Cox PH model
cox = CoxPHSurvivalAnalysis(alpha=0.0, ties="breslow")
cox.fit(X_train, y_train)
print(f"Coefficients: {cox.coef_}")  # log hazard ratios
# Hazard ratio interpretation: exp(coef) = HR for 1-unit increase
risk_scores = cox.predict(X_test)  # Higher = higher risk

# Survival function for individual patients
surv_funcs = cox.predict_survival_function(X_test[:3])
for fn in surv_funcs:
    print(f"5-year survival: {fn(365 * 5):.3f}")

# Penalized Cox (elastic net) -- for high-dimensional data (p > n)
coxnet = CoxnetSurvivalAnalysis(
    l1_ratio=0.9,           # 0=Ridge, 1=Lasso, between=Elastic Net
    alpha_min_ratio=0.01,   # smallest alpha / largest alpha ratio
    n_alphas=100,           # steps in regularization path
)
coxnet.fit(X_train, y_train)

# Feature selection: non-zero coefficients
selected = np.where(coxnet.coef_ != 0)[0]
print(f"Selected {len(selected)} / {X_train.shape[1]} features")

# IPCRidge: accelerated failure time model (predicts log survival time)
ipcridge = IPCRidge(alpha=1.0)
ipcridge.fit(X_train, y_train)
log_survival_time = ipcridge.predict(X_test)

from sksurv.ensemble import (
    RandomSurvivalForest,
    GradientBoostingSurvivalAnalysis,
    ComponentwiseGradientBoostingSurvivalAnalysis,
    ExtraSurvivalTrees,
)

# Random Survival Forest -- robust, minimal tuning
rsf = RandomSurvivalForest(
    n_estimators=200, min_samples_split=10, min_samples_leaf=15,
    max_features="sqrt", random_state=42, n_jobs=-1,
)
rsf.fit(X_train, y_train)
risk = rsf.predict(X_test)

# Gradient Boosting -- best performance, needs tuning
gbs = GradientBoostingSurvivalAnalysis(
    loss="coxph",           # "coxph" or "ipcwls" (AFT)
    n_estimators=300, learning_rate=0.05, max_depth=3,
    subsample=0.8, dropout_rate=0.1, random_state=42,
)
gbs.fit(X_train, y_train)

# ComponentwiseGB -- linear model with automatic feature selection
cgbs = ComponentwiseGradientBoostingSurvivalAnalysis(
    n_estimators=100, learning_rate=0.1,
)
cgbs.fit(X_train, y_train)
print(f"Non-zero coefficients: {np.sum(cgbs.coef_ != 0)}")

# ExtraSurvivalTrees -- more regularized than RSF, faster training
est = ExtraSurvivalTrees(n_estimators=100, random_state=42)
est.fit(X_train, y_train)

# Survival curves from any ensemble model
surv_funcs = rsf.predict_survival_function(X_test[:1])
chf_funcs = rsf.predict_cumulative_hazard_function(X_test[:1])

from sksurv.svm import FastSurvivalSVM, FastKernelSurvivalSVM, HingeLossSurvivalSVM

# Linear SVM -- fast, for linear relationships
lsvm = FastSurvivalSVM(alpha=1.0, rank_ratio=1.0, max_iter=100, random_state=42)
lsvm.fit(X_train_scaled, y_train)
risk = lsvm.predict(X_test_scaled)

# Kernel SVM -- for non-linear relationships (rbf, poly, sigmoid)
ksvm = FastKernelSurvivalSVM(
    alpha=1.0, kernel="rbf", gamma="scale",
    max_iter=50, random_state=42,
)
ksvm.fit(X_train_scaled, y_train)

# Hinge loss variant
hsvm = HingeLossSurvivalSVM(alpha=1.0, random_state=42)
hsvm.fit(X_train_scaled, y_train)
# NaiveSurvivalSVM also available but slower (O(n^3))

from sksurv.kernels import ClinicalKernelTransform

# Clinical kernel: combines clinical + molecular features
# Weighs clinical variables separately from high-dimensional molecular data
transform = ClinicalKernelTransform(fit_once=True)
transform.prepare(X_train)  # auto-detect clinical features
X_kern = transform.fit_transform(X_train)

from sksurv.nonparametric import kaplan_meier_estimator, nelson_aalen_estimator
import matplotlib.pyplot as plt

# Kaplan-Meier survival curve
time_km, surv_prob = kaplan_meier_estimator(y["event"], y["time"])
plt.step(time_km, surv_prob, where="post")
plt.xlabel("Time (days)")
plt.ylabel("Survival probability")
plt.title("Kaplan-Meier Estimate")
plt.savefig("km_curve.png", dpi=150)

# With confidence intervals
time_km, surv_prob, conf_int = kaplan_meier_estimator(
    y["event"], y["time"], conf_type="log-log",
)

# Nelson-Aalen cumulative hazard
time_na, cum_hazard = nelson_aalen_estimator(y["event"], y["time"])

# Stratified KM by group
for group_name, mask in [("Treated", treated_mask), ("Control", control_mask)]:
    t, s = kaplan_meier_estimator(y[mask]["event"], y[mask]["time"])
    plt.step(t, s, where="post", label=group_name)
plt.legend()

from sksurv.metrics import (
    concordance_index_censored,
    concordance_index_ipcw,
    cumulative_dynamic_auc,
    integrated_brier_score,
    brier_score,
    as_concordance_index_ipcw_scorer,
    as_integrated_brier_score_scorer,
)
import numpy as np

risk_scores = model.predict(X_test)

# Harrell's C-index (simple, biased with high censoring)
c_harrell = concordance_index_censored(
    y_test["event"], y_test["time"], risk_scores
)[0]

# Uno's C-index (recommended -- robust to censoring)
c_uno = concordance_index_ipcw(y_train, y_test, risk_scores)[0]
print(f"C-index (Harrell): {c_harrell:.3f}, (Uno): {c_uno:.3f}")

# Time-dependent AUC at clinically relevant timepoints
times = np.array([365, 730, 1095])  # 1, 2, 3 years
auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, risk_scores, times)
print(f"AUC at 1/2/3yr: {auc}, Mean: {mean_auc:.3f}")

# Integrated Brier Score (measures discrimination + calibration)
surv_funcs = model.predict_survival_function(X_test)
preds = np.row_stack([fn(times) for fn in surv_funcs])
ibs = integrated_brier_score(y_train, y_test, preds, times)
print(f"IBS: {ibs:.3f}")  # Lower is better; compare vs KM baseline

from sksurv.nonparametric import cumulative_incidence_competing_risks
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.util import Surv
import numpy as np

# Cumulative Incidence Function (CIF) estimation
# y_competing: structured array with integer event codes
# (0=censored, 1=relapse, 2=death_in_remission)
time_pts, cif_relapse, cif_death = cumulative_incidence_competing_risks(y_competing)

import matplotlib.pyplot as plt
plt.step(time_pts, cif_relapse, where="post", label="Relapse")
plt.step(time_pts, cif_death, where="post", label="Death in remission")
plt.xlabel("Time")
plt.ylabel("Cumulative Incidence")
plt.legend()

# Cause-specific hazard models: fit separate Cox per event type
event_types = np.array([0, 1, 2, 1, 0, 2, 1])
times = np.array([10.2, 5.3, 8.1, 3.7, 12.5, 6.8, 4.2])

# Model for relapse (event type 1); treat type 2 as censored
y_relapse = Surv.from_arrays(event=(event_types == 1), time=times)
cox_relapse = CoxPHSurvivalAnalysis()
cox_relapse.fit(X, y_relapse)

# Model for death (event type 2); treat type 1 as censored
y_death = Surv.from_arrays(event=(event_types == 2), time=times)
cox_death = CoxPHSurvivalAnalysis()
cox_death.fit(X, y_death)
print(f"Relapse HR (age): {np.exp(cox_relapse.coef_[0]):.3f}")
print(f"Death HR (age): {np.exp(cox_death.coef_[0]):.3f}")

High-dimensional (p > n)?
├── Yes → CoxnetSurvivalAnalysis (elastic net)
└── No
    ├── Need interpretable coefficients?
    │   ├── Yes → CoxPHSurvivalAnalysis or ComponentwiseGB
    │   └── No
    │       ├── Large data (n > 1000) → GradientBoostingSurvivalAnalysis
    │       ├── Medium data → RandomSurvivalForest or FastKernelSurvivalSVM
    │       └── Small data → RandomSurvivalForest
    └── For max performance → compare multiple models

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV, train_test_split
from sksurv.datasets import load_breast_cancer
from sksurv.ensemble import GradientBoostingSurvivalAnalysis
from sksurv.metrics import (
    concordance_index_ipcw, cumulative_dynamic_auc,
    integrated_brier_score, as_concordance_index_ipcw_scorer,
)
import numpy as np

X, y = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Build pipeline
pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("model", GradientBoostingSurvivalAnalysis(random_state=42)),
])

# Hyperparameter search
param_grid = {
    "model__learning_rate": [0.01, 0.05, 0.1],
    "model__n_estimators": [100, 200],
    "model__max_depth": [3, 5],
}
cv = GridSearchCV(pipe, param_grid, cv=5,
                  scoring=as_concordance_index_ipcw_scorer(), n_jobs=-1)
cv.fit(X_train, y_train)
print(f"Best CV C-index: {cv.best_score_:.3f}, Params: {cv.best_params_}")

# Final evaluation on test set
best = cv.best_estimator_
risk = best.predict(X_test)
c_uno = concordance_index_ipcw(y_train, y_test, risk)[0]
times = np.percentile(y_test["time"][y_test["event"]], [25, 50, 75])
auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, risk, times)
print(f"Test C-index: {c_uno:.3f}, Mean AUC: {mean_auc:.3f}")

from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.ensemble import RandomSurvivalForest, GradientBoostingSurvivalAnalysis
from sksurv.svm import FastSurvivalSVM
from sksurv.metrics import concordance_index_ipcw, integrated_brier_score
from sklearn.preprocessing import StandardScaler
import numpy as np

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s = scaler.transform(X_test)

models = {
    "CoxPH": CoxPHSurvivalAnalysis(),
    "RSF": RandomSurvivalForest(n_estimators=200, random_state=42),
    "GBS": GradientBoostingSurvivalAnalysis(n_estimators=200, random_state=42),
    "LinearSVM": FastSurvivalSVM(alpha=1.0, random_state=42),
}

times = np.percentile(y_test["time"][y_test["event"]], [25, 50, 75])
for name, model in models.items():
    X_tr = X_train_s if name in ("CoxPH", "LinearSVM") else X_train
    X_te = X_test_s if name in ("CoxPH", "LinearSVM") else X_test
    model.fit(X_tr, y_train)
    risk = model.predict(X_te)
    c = concordance_index_ipcw(y_train, y_test, risk)[0]
    print(f"{name:12s}: C-index = {c:.3f}")