Survival Analysis

Survival analysis models the time until an event occurs — and critically, handles observations where the event hasn't happened yet (censoring). Standard regression can't handle censoring correctly; survival models do. Applications include customer churn (time to cancellation), credit default (time to delinquency), employee attrition, clinical trials, and equipment failure.

1. Censoring

Censoring arises when the event of interest has not been observed for some subjects by the end of the study or data collection window. There are three main types:

• Right censoring (most common): the event hasn't occurred by end of the observation period. The customer is still active; the patient is still alive. We know their survival time is at least T.
• Left censoring: the event occurred before observation began. We know it happened, but not exactly when.
• Interval censoring: the event occurred within a known time window, but the exact time is unknown.

Why Naive Approaches Fail

• Dropping censored observations introduces selection bias and underestimates true survival time — the censored observations are not missing at random; they tend to be longer-surviving.
• Treating end-of-observation as the event time artificially shortens survival times for censored subjects, biasing all estimates downward.

Correct Approach

Maximum likelihood estimation explicitly models censored observations as "survived at least this long," contributing the probability $P(T \geq t_i)$ to the likelihood rather than a point mass at $t_i$.

Example DataFrame with right censoring (event_occurred=0 means censored):

2. Survival and Hazard Functions

Two mathematically equivalent representations fully describe a time-to-event distribution.

Survival Function

$$S(t) = P(T > t) = 1 - F(t)$$

Starts at $S(0) = 1$, is non-increasing, and approaches 0 as $t \to \infty$. It gives the probability a subject survives beyond time $t$.

Hazard Function

The instantaneous event rate at time $t$, given survival up to $t$:

$$h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t \mid T \geq t)}{\Delta t}=\frac{f(t)}{S(t)}$$

Cumulative Hazard and Key Relationship

$$H(t) = \int_0^t h(u)\, du = -\log S(t) \qquad \Longrightarrow \qquad S(t) = \exp\!\bigl(-H(t)\bigr)$$

Given any one of $S(t)$, $h(t)$, or $H(t)$, the others follow directly.

Common Parametric Distributions

3. Kaplan-Meier Estimator

The Kaplan-Meier (KM) estimator is the non-parametric maximum likelihood estimator of $S(t)$. It requires no distributional assumptions and is the standard first step in any survival analysis. The estimate is a step function that drops only at observed event times.

At each ordered event time $t_i$:

$$\hat{S}(t) = \prod_{t_i \leq t} \left(1 - \frac{d_i}{n_i}\right)$$

where $d_i$ is the number of events at $t_i$ and $n_i$ is the number at risk just before $t_i$. Censored observations reduce the risk set $n_i$ without contributing to the numerator.

The log-rank test compares two or more KM curves under the null hypothesis that the survival functions are equal. It is non-parametric and does not assume proportional hazards.

from lifelines import KaplanMeierFitter
from lifelines.statistics import logrank_test
import pandas as pd

kmf = KaplanMeierFitter()
kmf.fit(durations=df['duration_months'], event_observed=df['churned'], label='All customers')
kmf.plot_survival_function()
print(f"Median survival: {kmf.median_survival_time_:.1f} months")

# Compare two groups
mask = df['segment'] == 'premium'
kmf_prem = KaplanMeierFitter().fit(df[mask]['duration_months'], df[mask]['churned'], label='Premium')
kmf_std  = KaplanMeierFitter().fit(df[~mask]['duration_months'], df[~mask]['churned'], label='Standard')

results = logrank_test(df[mask]['duration_months'], df[~mask]['duration_months'],
                       event_observed_A=df[mask]['churned'],
                       event_observed_B=df[~mask]['churned'])
print(f"Log-rank p-value: {results.p_value:.4f}")

4. Cox Proportional Hazards Model

The Cox model is the workhorse of survival regression. It is semi-parametric: the baseline hazard $h_0(t)$ is left completely unspecified, while the effect of covariates enters multiplicatively:

$$h(t \mid \mathbf{x}) = h_0(t) \cdot \exp\!\bigl(\boldsymbol{\beta}^T \mathbf{x}\bigr)$$

Key Properties

• $h_0(t)$ is the unspecified baseline hazard — it cancels out of the partial likelihood, so $\boldsymbol{\beta}$ can be estimated without estimating it.
• $\exp(\beta_j)$ is the hazard ratio for feature $j$: HR > 1 means higher hazard (shorter expected survival); HR < 1 means lower hazard (longer expected survival).
• Parameters are fit by partial likelihood — a clever construction that conditions on the observed event ordering.
• Proportional hazards assumption: hazard ratios are constant over time. The effect of a covariate multiplies the baseline hazard by the same fixed factor at every time point.

Testing the PH Assumption

• Schoenfeld residuals: regress Schoenfeld residuals on time; a statistically significant slope indicates the hazard ratio is not constant.
• Log-log plot: plot $\log(-\log \hat{S}(t))$ vs. $\log(t)$ for each group; parallel lines indicate proportional hazards.

from lifelines import CoxPHFitter

cph = CoxPHFitter()
cph.fit(df[['duration_months', 'churned', 'age', 'monthly_spend', 'segment_premium']],
        duration_col='duration_months', event_col='churned')

cph.print_summary()
# Shows: coef, exp(coef)=HR, 95% CI, p-value for each feature

# Check proportional hazards assumption
cph.check_assumptions(df_clean, p_value_threshold=0.05)

Interpretation example: "A hazard ratio of 1.35 for monthly_spend means customers with one unit higher spend have 35% higher hazard of churning at any given time, assuming proportional hazards."

5. When PH Assumption Fails and ML Approaches

Remedies When PH Fails

• Stratification: allow a separate baseline hazard for each stratum of the offending covariate. The stratified Cox model estimates separate $h_0(t)$ per stratum while sharing $\boldsymbol{\beta}$.
• Time-varying covariates: interact the covariate with a function of time (e.g., $x \cdot \log(t)$) to let the hazard ratio evolve.
• Parametric models: choose a flexible parametric family (e.g., generalized gamma) that doesn't impose proportional hazards.

ML-Based Survival Models

• Random Survival Forests: extend random forests to survival outcomes; capture non-linear interactions and non-proportional hazards; output a full estimated survival function per subject.
• XGBoost with Cox objective: fast, scalable, handles millions of rows; optimizes the Cox partial likelihood as its loss function.
• scikit-survival: sklearn-compatible API providing RSF, gradient boosted survival trees, and survival SVMs.

from sksurv.ensemble import RandomSurvivalForest
import numpy as np

# y must be structured array with (event_occurred, duration)
y = np.array([(bool(e), t) for e, t in zip(df['churned'], df['duration_months'])],
             dtype=[('event', bool), ('time', float)])

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

# Concordance index (C-index): primary metric for survival models
# Proportion of comparable pairs where model correctly orders survival times
c_index = rsf.score(X_test, y_test)
print(f"C-index: {c_index:.3f}")  # 0.5 = random, 1.0 = perfect

6. Applications and Business Framing

C-index (Concordance Index)

The primary evaluation metric for survival models. Given all comparable pairs $(i, j)$ where subject $i$ experienced the event before subject $j$ and both were at risk at that time, the C-index is the fraction of pairs the model ranks correctly (assigns higher risk to $i$).

• C-index = 0.5: random (no discriminative ability)
• C-index = 1.0: perfect ranking
• Analogous to AUC-ROC for binary classification — AUC-ROC is actually the C-index for binary outcomes at a fixed time horizon
• Key difference from AUC-ROC: C-index handles censoring by excluding incomparable pairs from the denominator; AUC-ROC requires a fixed prediction horizon

Practice Questions