A long hospital ward in a vintage painterly style is engulfed in creeping darkness, with many occupied beds on the left; at the center, a patient, a nurse, a doctor, and other figures form a taut human chain across a deep chasm, linking the dying ward to a bright exit marked “Recovery.”

Survival Analysis

Survival Analysis is a statistical framework for studying the time until an event of interest occurs, such as death, relapse, reoperation, device failure, or hospital readmission.

Unlike standard regression methods, survival analysis explicitly accounts for censoring, time-dependent risk, and non-Gaussian distributions, making it essential for clinical research, biomedical engineering, and reliability studies.

This article provides a concise yet technically rigorous overview of the main concepts, mathematical foundations, dataset organization principles, and a practical Python workflow for Kaplan–Meier estimation, Cox regression, and survival curve interpretation.

In brief

  • Survival analysis studies time‑to‑event data with censoring and time‑varying risk.
  • Key functions: survival S(t), hazard h(t), cumulative hazard H(t) with S(t)=e^{-H(t)}.
  • Kaplan–Meier estimates S(t). Nelson–Aalen estimates H(t).
  • Cox PH models covariate effects via hazard ratios, assuming proportional hazards.
  • Python workflow uses lifelines for KM, Cox, diagnostics, and plots with number‑at‑risk. </aside>

✅ Introduction

Survival analysis is used whenever we are interested not only in whether an event happens, but when it happens.

This applies to dozens of medical and scientific questions:

  • How long do patients survive after surgery?
  • How quickly do devices fail after implantation?
  • When does cancer recur after treatment?
  • How long do patients remain event-free during follow-up?

Unlike standard regression or classification problems, survival data have a unique property:

not all events are observed. Some patients reach the end of follow-up without experiencing the event, or they are lost to follow-up.

timeline with events and a censoring case

We do not want to discard these patients, because they still provide valuable information.

To deal with this, survival analysis uses a family of tools that:

  1. correctly handle partially observed outcomes (censoring),
  2. describe how risk evolves over time,
  3. and allow us to study how predictors influence the timing of the event.

Before diving into mathematics, it is essential to build an intuitive mental model of how survival methods work.

Two ideas form the foundation of nearly everything in survival analysis:

  • Kaplan–Meier curves, which describe how many individuals remain event-free over time,
  • Cox proportional hazards models, which describe how predictors affect the risk of experiencing the event.

We start with these two concepts using intuition and visual reasoning, and only later introduce the mathematical objects behind them.

✅ 1. Kaplan–Meier: An Intuitive Explanation

Imagine following 100 patients after a cardiac operation.

Every day, some patients may experience the event (e.g., death), while others remain alive; some may stop attending follow-up.

If we plot the count of patients who are still alive over time, we obtain a stepwise survival curve.

This is essentially what the Kaplan–Meier estimator does.

What problem does Kaplan–Meier solve?

If we simply compute “percentage alive at 1 year”, we would have to exclude censored patients.

This wastes information.

Kaplan–Meier instead uses all available follow-up:

  • when an event occurs, the curve steps down,
  • when a patient is censored, the curve stays flat but that person is no longer considered “at risk” afterwards.

This gives a non-parametric summary of survival over time that correctly accounts for incomplete follow-up.

How to read a Kaplan–Meier curve

A KM curve answers:

  • What proportion of the cohort is still event-free at time ttt?
  • When does the curve drop most steeply? (periods of higher risk)
  • What is the median survival? (the time when survival reaches 50%)

Comparing groups

If we split patients by a factor — for example, diabetics vs non-diabetics — we can draw two KM curves:

  • If the curves separate early → one group experiences events sooner.
  • If the curves remain close → groups have similar survival.
  • The log-rank test provides a statistical comparison of the curves.

Why is KM important?

KM is the entry point to survival analysis because it:

  • visualises the survival experience of the cohort,
  • handles censored data naturally,
  • reveals temporal patterns (early vs late risk),
  • provides clinical intuition before modelling.

It is descriptive, not inferential.

It tells us what happens, not why.

To understand why some patients have higher risk, we need a model.

✅ 2. Cox Proportional Hazards: An Intuitive Explanation

While Kaplan–Meier describes survival without accounting for predictors, the Cox proportional hazards model explains how predictors modify the risk over time.

The key idea

Cox does not try to predict the “exact survival time”.

Instead, it compares individuals who are at risk at the same moment, and quantifies how much more (or less) likely one patient is to experience the event compared to another.

This comparison produces the hazard ratio (HR).

What is the hazard ratio? (Intuition)

If two patients are identical except for one characteristic:

  • HR = 2 → one patient has twice the instantaneous risk of experiencing the event at any time
  • HR = 0.5 → one patient has half the instantaneous risk
  • HR = 1 → no difference

Importantly:

  • The Cox model does not assume anything about the shape of the baseline risk over time. Early risk, late risk, U-shaped risk — all allowed.
  • It only assumes that the ratio between risks of two individuals stays constant over time. This is the proportional hazards assumption.

Why Cox is powerful

Cox allows us to:

  • quantify the effect of continuous predictors (age, creatinine, EuroSCORE II),
  • adjust for multiple risk factors simultaneously,
  • estimate adjusted survival curves for selected profiles,
  • test whether predictors violate proportionality (Schoenfeld residuals).

What Cox does not do

  • It does not estimate absolute survival probability without additional steps.
  • It does not assume a distribution (e.g., exponential, Weibull).
  • It does not tell you the baseline hazard explicitly (though it can be estimated).

KM vs Cox: intuition

GoalMethod
How many patients survive over time?Kaplan–Meier
Do groups differ in unadjusted survival?Kaplan–Meier + log-rank
How do multiple predictors affect risk?Cox model
How does risk accumulate?Nelson–Aalen

Having built an intuitive understanding of how survival curves and hazard ratios work, we can now formalize these ideas using their underlying mathematical functions.

✅ Mathematical Foundations of Survival Analysis

Survival analysis relies on mathematical functions that describe time-to-event data. Understanding these interconnected functions is essential for interpreting survival curves, hazards, and Cox models.

Time-to-Event Variable

Let:

    \[T \ge 0\]

be a non-negative random variable representing the time until the event of interest (death, relapse, reoperation, etc.).

  • T = 0 represents the time origin (baseline: surgery, diagnosis, randomization).
  • T = t represents the time the event occurs.

Distribution Function and Density

Cumulative Distribution Function (CDF)

    \[F(t) = \Pr(T \le t)\]

  • Represents the probability that the event has occurred by time t.
  • It is a non-decreasing function, starting at 0 and converging to 1 (if all eventually experience the event).

Probability Density Function (PDF)

If T is continuous, it admits a density:

    \[f(t) = \frac{dF(t)}{dt}\]

  • Represents the instantaneous accumulation rate of event probability around time t.
  • It is not itself a probability, but integrates to 1.

Survival Function

The survival function is the central object of survival analysis:

    \[S(t) = \Pr(T > t)\]

This means:

  • The probability that an individual survives beyond time t without experiencing the event.
  • It is a non-increasing function:
    • S(0) = 1 (everyone is alive at time zero),
    • S(t) to 0 as t to infinity (if the event is certain),
    • flatter curves → better survival.

Relationship with the CDF:

    \[S(t) = 1 - F(t)\]

Hazard Function (Instantaneous Risk)

The hazard function describes the instantaneous risk of experiencing the event at time t, given survival up to that time.

    \[h(t) = \lim_{\Delta t\to 0} \frac{\Pr(t \le T < t + \Delta t \mid T \ge t)}{\Delta t}\]

Interpretation:

  • It is not a probability: hazards can exceed 1.
  • It measures how dangerous the moment t is for someone who has survived up to t.
  • High hazard at t: steep drop in survival curve at t.

Relationship with the PDF and survival function:

    \[h(t) = \frac{f(t)}{S(t)}\]

Meaning:

  • The instantaneous density of events divided by the probability of still being at risk.

Cumulative Hazard Function

Defined as:

    \[H(t) = \int_0^t h(u)\,du\]

  • Aggregates the total hazard accumulated up to time t.
  • It is a non-decreasing function.
  • Has a key link to survival:

    \[S(t) = \exp[-H(t)]\]

Meaning:

  • Survival decreases exponentially with cumulative hazard.

Relationship Between Survival and Hazard

Because the hazard, survival, density, and cumulative hazard all derive from one another, we can express them interchangeably. For instance:

    \[H(t) = \int_0^t \frac{f(u)}{S(u)}\,du\]

This expresses the cumulative hazard in terms of the PDF and survival.

Median Survival Time

The median survival time is the earliest time at which the survival function drops to 0.5:

    \[\text{Median survival} = \inf \{ t : S(t) \le 0.5 \}\]

This measure is robust and often preferred to the mean in right-skewed distributions.

Cox Proportional Hazards Model

The Cox model relates survival time to covariates while leaving the baseline hazard unspecified:

    \[h(t \mid X_i) = h_0(t)\exp(\beta_1 X_{i1} + \beta_2 X_{i2} + \dots + \beta_p X_{ip})\]

Here:

  • h0(t) is the unspecified baseline hazard.
  • β coefficients measure the effect of covariates on the hazard.
  • Covariates may be binary (e.g., sex), continuous (e.g., age), or categorical.

The central quantity is the hazard ratio (HR):

    \[\mathrm{HR} = \exp(\beta)\]

Interpretation:

  • HR > 1: increased hazard (worse prognosis).
  • HR < 1: decreased hazard (protective).
  • HR = 2: hazard doubled relative to baseline.

Proportional Hazards Assumption

The defining assumption of the Cox model is that hazard ratios are constant over time. Formally:

    \[\frac{h(t \mid X_1)}{h(t \mid X_2)} = \exp\big[\beta^\top (X_1 - X_2)\big]\]

Meaning:

  • The effect of covariates multiplies the hazard but does not depend on time t.
  • If hazards cross visibly in KM curves, the assumption is likely violated.

Complete Cheat Sheet (All Formulas Together)

For convenience, here is a single block containing the full set of fundamental equations:

    \[\begin{aligned} T &\ge 0 \\ F(t) &= \Pr(T \le t) \\ f(t) &= \frac{dF(t)}{dt} \\ S(t) &= \Pr(T > t) = 1 - F(t) \\ h(t) &= \lim_{\Delta t\to 0} \frac{\Pr(t \le T < t + \Delta t \mid T \ge t)}{\Delta t} \\ h(t) &= \frac{f(t)}{S(t)} \\ H(t) &= \int_0^t h(u)\,du \\ S(t) &= \exp[-H(t)] \\ \text{Median survival} &= \inf \{ t : S(t) \le 0.5 \} \\ h(t \mid X) &= h_0(t)\exp(\beta^\top X) \\ \mathrm{HR} &= \exp(\beta) \end{aligned}\]

✅ Nelson–Aalen: A Simple View of How Risk Accumulates Over Time

While the Kaplan–Meier estimator describes how many individuals remain event-free at each time point, the Nelson–Aalen estimator describes how much risk has accumulated up to that moment.

It is the natural companion to Kaplan–Meier, but it looks at survival from the opposite direction.

Instead of asking:

What is the probability of still being alive at time t?

the Nelson–Aalen estimator answers:

How much hazard has the cohort accumulated by time ttt?

Every time an event occurs, the cumulative hazard increases by a small step.

If many patients experience the event while few remain at risk, the step is larger;

if the risk is low or no events occur, the curve stays flat.

This produces a monotonic, stepwise curve that grows whenever events happen and remains unchanged during event-free periods.

In practice:

  • a steep Nelson–Aalen curve = a period of high event intensity,
  • a flat curve = a low-risk interval.

Although clinicians tend to report Kaplan–Meier survival curves, Nelson–Aalen is essential behind the scenes: it forms the basis for estimating the baseline hazard in the Cox model and is useful for understanding when the cohort faces elevated risk.

Kaplan-meier vs nelson-aalen line graph

This cumulative-hazard perspective is also the foundation for estimating the Cox model’s baseline hazard, which explains why Nelson–Aalen and Cox are often presented together

✅Assumptions of Survival Analysis

A correct survival analysis depends critically on a small set of fundamental assumptions. If these assumptions are satisfied, survival estimates, hazard ratios, and model diagnostics have a clear interpretation. If they are violated, the results may become biased or misleading.

Time Origin, Time Scale, and Follow-Up

Every survival model begins with a clear definition of time zero, because all subsequent assumptions depend on it.

Time zero may be:

  • the date of surgery,
  • the date of diagnosis,
  • the date of randomization,
  • the date of discharge.

Once the origin is fixed, every subject must have their follow-up measured consistently from that same moment. Likewise, the time scale (days, months, years) must be the same for all individuals.

If different patients effectively start counting from different types of events (e.g., some from diagnosis, others from treatment), even perfect statistical techniques will give distorted results, because the underlying timeline is inconsistent.

In short:

  • A single, shared time origin
  • A single, coherent time scale
  • A consistent definition of follow-up

These form the foundation on which survival models can operate meaningfully.

Event Observation and Censoring: What We Know About Each Patient

In survival data, two quantities are recorded for every subject:

  • How long they were observed
  • Whether they experienced the event or not
  • If the event does not occur within the observation window, the observation is right-censored.

This leads to one of the most important assumptions of all survival models:

Survival analysis assumes that the reason a subject becomes censored does not carry additional information about their underlying risk—beyond what is already measured in the dataset.

In other words, a censored patient should be statistically similar to a non-censored patient with the same covariates.

This assumption is satisfied when:

  • the study ends for everyone at the same administrative date;
  • patients miss follow-up for reasons unrelated to prognosis;
  • the time of censoring is independent of the time of the event.

It is violated when:

  • frail or critically ill patients systematically drop out earlier;
  • healthier patients preferentially attend follow-up visits;
  • loss to follow-up is correlated with the progression of disease.

If censoring is informative, Kaplan–Meier and Cox estimates tend to be biased, often in optimistic directions.

Independence of Observations

Most classical survival methods assume that observations are independent:

  • each row corresponds to a single subject,
  • subjects do not influence each other’s risks.

In real data this assumption may fail:

  • patients treated within the same hospital or by the same surgeon may share unmeasured characteristics;
  • individuals in family or cluster settings may have correlated risks;
  • repeated events within the same subject require specialized modeling.

When independence is doubtful, robust standard errors, frailty terms, or mixed-effects models may be appropriate.

In standard clinical datasets (e.g., single event per patient), independence is usually a reasonable working assumption.

Assumptions for Kaplan–Meier Estimation

Kaplan–Meier is mathematically simple and robust. Its assumptions are:

  1. Censoring is non-informative (as explained above).
  2. Event times are precisely known (not interval-censored).
  3. Observations are independent.
  4. Competing risks are either absent or treated as censoring.

Note: when competing risks are substantial, Kaplan–Meier can overestimate the incidence of the event of interest. Use the Cumulative Incidence Function and, for regression, Fine–Gray models.

Quando queste condizioni tengono, la curva di Kaplan–Meier fornisce un quadro pulito e non distorto della sopravvivenza nel tempo.

Assumption of the Cox Model

The Cox proportional hazards model introduces a major simplification:

it assumes that risk factors multiply the hazard by a constant factor that does not vary with time.

In other words:

  • If a certain covariate doubles the hazard at day 30, it also doubles it at day 300 and at day 1000.
  • The baseline hazard may change with time in any arbitrary way, but the effect of the covariates remains proportional across the entire timeline.

This is the Proportional Hazards (PH) assumption.

When PH holds, hazard ratios have a simple and intuitive interpretation.

When PH is violated (e.g., hazards cross visibly), Cox estimates may remain “statistically significant” but become scientifically misleading, because they represent an average of effects that change over time.

The Cox model also assumes that continuous covariates influence the log-hazard linearly.

For example:

  • age is assumed to increase the log-hazard at a constant rate per year;
  • creatinine is assumed to have a smooth, linear effect.

But biological processes are rarely linear.

When continuous covariates have nonlinear effects, analysts may need:

  • splines,
  • polynomial terms,
  • piecewise transformations.

Correct functional form is critical; otherwise hazard ratios may misrepresent the underlying relationship.

In many clinical settings, important variables change over time (e.g., postoperative complications, treatment changes, biomarkers).

Standard Cox regression assumes:

  • either these time-dependent influences are small,
  • or they are unrelated to the predictors in the model.

When important time-varying confounders exist, the Cox model is misspecified and advanced methods are needed:

  • extended Cox models with time-dependent covariates,
  • joint models,
  • marginal structural models.

Summary of Assumptions

Instead of memorizing lists, the assumptions of survival analysis form a coherent logic:

  • Define time correctly
    • same time zero, same time scale.
  • Observe follow-up consistently
    • no informative censoring.
  • Ensure independence or adjust for clustering
    • independence is the default assumption.
  • Use appropriate model assumptions
    • Kaplan–Meier requires minimal assumptions.
    • Cox additionally requires:
      • proportional hazards,
      • linearity of log-hazard,
      • no major time-dependent confounding.
  • Check assumptions diagnostically
    • graphs, residuals, and tests help confirm model validity.

When these principles are respected, survival analysis provides reliable, interpretable, and clinically meaningful estimates of risk and prognosis.

✅ Dataset Organization for Survival Analysis

A robust survival analysis starts long before writing code. It starts with precise design of the dataset. Here we detail how to define variables, structure the data, and prepare them for Python.

Core Variables: Time and Event

The two fundamental columns in any survival dataset are:

  • Duration variable (time, follow_up, etc.)
    • Represents the time elapsed between a clearly defined start and end:
      • from surgery to death
      • from diagnosis to relapse
      • from discharge to readmission
    • Must be non-negative and in a consistent unit (days, months, years).
  • Event indicator (event, status, etc.)
    • Binary:
      • 1 = event occurred
      • 0 = censored
    • Must clearly reflect the event definition used in your study.
time: 365
event: 1   # death within 1 year

time: 540
event: 0   # alive at last follow-up (censored)

Dates vs. Durations

In clinical databases, you often start with dates, not durations:

  • date_surgery
  • date_last_follow_up
  • date_death (if applicable)

You then derive:

  • time = difference between two dates (e.g., days)
  • event = 1 if death occurred before or at date_last_follow_up, 0 otherwise

Python example:

import pandas as pd

df = pd.read_csv("raw_clinical_data.csv", parse_dates=["date_surgery", "date_last_fu", "date_death"]) 

# duration in days from surgery to last observed time (death or last FU)
df["time"] = (df["date_last_fu"] - df["date_surgery"]).dt.days

# event = 1 if death date is not null and <= last follow-up, else 0
df["event"] = ((df["date_death"].notna()) &
               (df["date_death"] <= df["date_last_fu"])).astype(int)

This step is critical: any mistake in date handling will directly affect survival times and bias the analysis.

Covariates: Types and Coding

Beyond time and event, your dataset typically includes covariates that may influence risk:

1. Continuous variables

Examples:

  • age (years)
  • creatinine (mg/dL)
  • left_ventricular_ejection_fraction

Recommendations:

  • Keep units clinically meaningful.
  • Consider transformations when distributions are highly skewed (e.g., log(creatinine)).

2. Categorical variables

Examples:

  • sex (male/female)
  • diabetes (yes/no)
  • NYHA_class (I, II, III, IV)

In Python and statistical modeling, these are typically converted into:

  • Binary indicators (dummy variables) for each category (or all but one reference category). Example: sex → sex_male (1 = male, 0 = female as reference).

With pandas and lifelines, you can often pass categorical variables directly and let the library handle encoding, but explicit control is often preferable in clinical research.

3. Ordinal variables

Variables that have a natural order but not a precise numerical distance:

  • NYHA
  • Killip class
  • tumor_stage (I–IV)

Options:

  • Treat as continuous (if the scale is approximately linear in risk),
  • Or convert to dummy variables (e.g., NYHA III–IV vs I–II),
  • Or model with splines if risk increases non-linearly.

Explicitly justify your choice in a scientific context.

Baseline vs Time-Dependent Covariates

In basic survival analysis, we often use baseline covariates:

  • measured once at time zero (e.g., age, sex, preoperative risk scores).

However, some variables change over time:

  • hemodynamic parameters,
  • occurrence of complications,
  • treatment changes,
  • lab values.

If these time-varying variables are important, you need a time-dependent survival dataset, often in long format:

  • One row per subject per time interval
  • With updated covariate values

In this introductory lesson, we focus on baseline-only models.

Wide vs Long Format Datasets

Two common structures:

  • Wide format (most common in simple survival models)
    • One row per patient
    • Columns: id, time, event, age, sex, bmi, euroscore, etc.
    • Suitable for:
      • Kaplan–Meier
      • Cox PH with baseline covariates
  • Long format
    • Multiple rows per patient
    • Columns may include start_time, end_time, event, time-varying covariates
    • Used for:
      • time-dependent covariates,
      • recurrent events,
      • counting-process notation models.

For a first implementation with lifelines and standard Cox, wide format is usually sufficient.

Handling Missing Data

As in any regression model, missing data can cause bias and loss of power.

Options:

  • Complete-case analysis (exclude rows with missing covariates)
    • Simple but may bias results if missingness is not completely at random.
  • Imputation (e.g., multiple imputation)
    • More advanced but statistically preferable in serious clinical research.

In the dataset used for survival analysis, you should:

  1. Clearly document which variables have missing values.
  2. Decide a strategy (e.g., complete-case for a tutorial, multiple imputation in publications).
  3. Apply the same treatment across all models.

Minimal Dataset for a First Python Survival Analysis

A minimal but realistic dataset for your Python scripts might include:

  • time: follow-up in days from surgery to death/last contact
  • id: patient identifier
  • event: 1 if death occurred, 0 if censored
  • age: at surgery
  • sex: 0 = female, 1 = male
  • euroscore: predicted risk at baseline (continuous)
  • bmi: body mass index
  • diabetes: 0/1
  • copd: 0/1
  • creatinine: mg/dL

In CSV:

id,time,event,age,sex,euroscore,bmi,diabetes,copd,creatinine
1,365,1,72,1,3.4,27,1,0,1.10
2,540,0,81,0,7.1,25,0,1,1.80
3,120,1,65,1,2.1,28,1,0,0.90
4,900,0,59,1,1.2,31,0,0,1.00
...

This structure is immediately compatible with lifelines:

import pandas as pd
from lifelines import CoxPHFitter

df = pd.read_csv("survival_data.csv")

# Fit Cox model with multiple covariates
cph = CoxPHFitter()
cph.fit(df, duration_col="time", event_col="event",
        formula="age + sex + euroscore + bmi + diabetes + copd + creatinine")

cph.print_summary()

Reproducible Workflow: Raw → Clean → Analysis

  • Raw data
    • Direct export from the hospital information system, registry, or database
    • Contains dates, codes, free-text, missing or inconsistent entries
  • Clean analytic dataset for survival
    • Derived columns: time, event
    • Well-defined covariates
    • Clear handling of missing values
    • Saved as survival_analysis_dataset.csv
  • Model-specific subsets
    • e.g., a version with only complete cases for Cox
    • or a dataset restricted to a specific population (e.g., isolated CABG)

✅ Python Workflow: Scripts and Analysis

Importing libraries

import pandas as pd
from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.plotting import add_at_risk_counts
import matplotlib.pyplot as plt

Loading and preparing the dataset

df = pd.read_csv("survival_data.csv")

# Ensure correct dtypes
df["event"] = df["event"].astype(int)
df["time"] = df["time"].astype(float)

Kaplan–Meier Estimator

km = KaplanMeierFitter()
km.fit(durations=df["time"], event_observed=df["event"])

fig, ax = plt.subplots(figsize=(8, 6))
km.plot_survival_function(ax=ax)
add_at_risk_counts(km, ax=ax)
plt.title("Kaplan–Meier Survival Curve")
plt.xlabel("Time (days)")
plt.ylabel("Survival Probability")
plt.tight_layout()
plt.show()

Cox Proportional Hazards Model

cph = CoxPHFitter()
cph.fit(df, duration_col="time", event_col="event")

cph.print_summary()

Checking proportional hazards assumption

cph.check_assumptions(df, p_value_threshold=0.05)

Interpretation of Results

Kaplan–Meier

Look for:

  • median survival time
  • separation between strata (if groups are plotted)
  • width of confidence intervals
  • presence of early vs late risk

Interpretation example:

A steep early drop indicates a high early hazard (perioperative mortality), followed by plateau.

Cox Model

Key parameters:

  • Hazard Ratio (HR):
    • HR = exp(β) HR > 1 → increased risk; HR < 1 → protective effect.
  • p-value: significance of covariate
  • 95% CI: uncertainty
  • Concordance (C-index): predictive accuracy
  • Schoenfeld test: PH assumption check

Example interpretation:

Age (HR 1.04, p < 0.01) is a significant risk factor with a 4% increase in instantaneous mortality risk per year. EuroSCORE II demonstrates an expected positive association with hazard. BMI shows no significant effect.

Graphical Outputs

Kaplan–Meier Curves

  • Survival function
  • Confidence bands
  • Number-at-risk table (optional)

Cox Model Visualizations

  • Forest plot of hazard ratios
  • Schoenfeld residuals

✅ Complete Python Script – Survival Simulation + KM + Cox + Full Plots

"""
Survival Analysis Demo: Simulation + Kaplan–Meier + Cox PH Model
Author: Michele D. Pierri

This script creates a fully synthetic survival dataset, performs:
    - Kaplan–Meier curves
    - Stratified Kaplan–Meier curves
    - Nelson–Aalen cumulative hazard
    - Cox Proportional Hazards regression
    - Model diagnostics (Schoenfeld residuals)
    - Martingale & Deviance residual plots
with extremely verbose commentary for teaching purposes.

Libraries required:
    pandas, numpy, lifelines, matplotlib, seaborn
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from lifelines import KaplanMeierFitter, CoxPHFitter, NelsonAalenFitter
from lifelines.statistics import logrank_test

sns.set(style="whitegrid")


# ============================================================
# 1. Simulate a Survival Dataset
# ============================================================

# For reproducibility
np.random.seed(42)

N = 400  # number of subjects

# Continuous covariates
age = np.random.normal(loc=65, scale=10, size=N)         # age in years
bmi = np.random.normal(loc=27, scale=4, size=N)           # BMI
euroscore = np.clip(np.random.normal(3.5, 2.0, N), 0, 15) # simplified numeric risk score

# Binary covariates
sex = np.random.binomial(1, 0.60, N)  # 1=male, 0=female
diabetes = np.random.binomial(1, 0.25, N)

# True hazard coefficients (for simulation only)
beta_age = 0.03       # hazard increases with age
beta_sex = 0.40       # males higher hazard
beta_diabetes = 0.60  # diabetes increases hazard
beta_euro = 0.12      # risk score effect

# Linear predictor
lp = (
    beta_age * age +
    beta_sex * sex +
    beta_diabetes * diabetes +
    beta_euro * euroscore
)

# Convert linear predictor to an exponential hazard multiplier
hazard_multiplier = np.exp(lp / 20)

# Simulate event times using exponential distribution
base_time = np.random.exponential(scale=600, size=N)
time = base_time / hazard_multiplier  # adjusted by hazard

# Add random censoring
censoring_time = np.random.uniform(200, 1200, size=N)

# Observed time and event indicator
observed_time = np.minimum(time, censoring_time)
event = (time <= censoring_time).astype(int)

df = pd.DataFrame({
    "time": observed_time,
    "event": event,
    "age": age,
    "sex": sex,
    "bmi": bmi,
    "euroscore": euroscore,
    "diabetes": diabetes
})

print("\n=== HEAD OF SIMULATED DATASET ===")
print(df.head())
print("\nDataset size:", df.shape)


# ============================================================
# 2. Kaplan–Meier Estimation
# ============================================================

kmf = KaplanMeierFitter()

# Fit the KM estimator
kmf.fit(durations=df["time"], event_observed=df["event"], label="Overall Survival")

# Print numerical summary
print("\n=== KAPLAN–MEIER ESTIMATOR SUMMARY ===")
print(kmf.survival_function_.head())

# KM Plot
plt.figure(figsize=(8, 6))
kmf.plot(ci_show=True)
plt.title("Kaplan–Meier Survival Curve (Overall)")
plt.xlabel("Time (days)")
plt.ylabel("Survival Probability")
plt.savefig("km_overall.png", dpi=200)
plt.close()


# ============================================================
# 3. Stratified Kaplan–Meier (e.g., by sex)
# ============================================================

plt.figure(figsize=(8, 6))

for group in [0, 1]:
    mask = df["sex"] == group
    kmf.fit(df.loc[mask, "time"], df.loc[mask, "event"], label=f"Sex {group}")
    kmf.plot(ci_show=True)

plt.title("Kaplan–Meier Curves Stratified by Sex")
plt.xlabel("Time (days)")
plt.ylabel("Survival Probability")
plt.savefig("km_by_sex.png", dpi=200)
plt.close()

# Log-rank test
result = logrank_test(
    df.loc[df.sex == 0, "time"], df.loc[df.sex == 1, "time"],
    df.loc[df.sex == 0, "event"], df.loc[df.sex == 1, "event"]
)
print("\n=== LOG-RANK TEST (SEX) ===")
print(result.summary)


# ============================================================
# 4. Nelson–Aalen Cumulative Hazard
# ============================================================

naf = NelsonAalenFitter()
a f i t = naf.fit(df["time"], df["event"])  # keeps API explicit, avoid shadowing

plt.figure(figsize=(8, 6))
afit.plot()
plt.title("Nelson–Aalen Cumulative Hazard Estimate")
plt.xlabel("Time (days)")
plt.ylabel("Cumulative Hazard")
plt.savefig("nelson_aalen.png", dpi=200)
plt.close()


# ============================================================
# 5. Cox Proportional Hazards Model
# ============================================================

cph = CoxPHFitter()

# Fit Cox model with selected covariates
cph.fit(df, duration_col="time", event_col="event",
        formula="age + sex + diabetes + euroscore + bmi")

print("\n=== COX MODEL SUMMARY ===")
print(cph.summary)

# Forest plot of hazard ratios
plt.figure(figsize=(8, 6))
cph.plot()
plt.title("Cox Proportional Hazards – Hazard Ratios")
plt.savefig("cox_forest.png", dpi=200)
plt.close()


# ============================================================
# 6. Cox Diagnostics (Proportional Hazards Assumption)
# ============================================================

print("\n=== TESTING PH ASSUMPTION ===")
cph.check_assumptions(df, p_value_threshold=0.05, show_plots=True)

# All residual plots will be displayed or saved depending on the environment.


# ============================================================
# 7. Martingale and Deviance Residual Plots
# ============================================================

# Extract residuals
martingale_res = cph.compute_residuals(df, kind="martingale")
deviance_res = cph.compute_residuals(df, kind="deviance")

plt.figure(figsize=(8, 6))
plt.scatter(df["age"], martingale_res, alpha=0.7)
plt.title("Martingale Residuals vs Age")
plt.xlabel("Age")
plt.ylabel("Martingale Residual")
plt.savefig("martingale_age.png", dpi=200)
plt.close()

plt.figure(figsize=(8, 6))
plt.scatter(df["age"], deviance_res, alpha=0.7)
plt.title("Deviance Residuals vs Age")
plt.xlabel("Age")
plt.ylabel("Deviance Residual")
plt.savefig("deviance_age.png", dpi=200)
plt.close()


print("\n=== ALL ANALYSES COMPLETED SUCCESSFULLY ===\n")
print("Generated files:")
print("- km_overall.png")
print("- km_by_sex.png")
print("- nelson_aalen.png")
print("- cox_forest.png")
print("- martingale_age.png")
print("- deviance_age.png")
print("\nAll PH diagnostic plots printed by lifelines.")

Kaplan- Meier survival curve (Overall)

Kaplan-Meier Curves by Sex

Cox Model - Hazard Ratios

Scaled Schoenfeld residual of the "age" variable

✅ Conclusion

Survival analysis provides a rigorous framework for understanding not only whether an event occurs, but when it occurs—an aspect that is central to most clinical questions. By integrating Kaplan–Meier curves, cumulative hazard estimators such as Nelson–Aalen, and the Cox proportional hazards model, clinicians and researchers can quantify risk over time, adjust for multiple predictors, and evaluate prognosis with far greater nuance than traditional statistical tools allow.

In practice, these methods help identify periods of heightened vulnerability, clarify the role of specific risk factors, and support individualized decision-making. Whether estimating postoperative survival, assessing treatment effects, or stratifying patients by risk, survival analysis transforms longitudinal data into actionable clinical insight.

As with any method, careful attention to assumptions, data structure, and proper diagnostics is essential. But when applied correctly, survival analysis remains one of the most powerful and clinically meaningful tools for interpreting time-to-event data and informing evidence-based care.

✅ See also

✅ FAQ

  • What is the difference between survival S(t) and hazard h(t)? S(t) is the probability of surviving beyond time t. h(t) is the instantaneous event rate at time t among those still at risk. They are linked by S(t) = exp[−H(t)].
  • When does Kaplan–Meier overestimate risk? When competing risks are substantial. Prefer the Cumulative Incidence Function and Fine–Gray regression for the event of interest.
  • How do I check the proportional hazards assumption? Inspect Schoenfeld residuals, use cph.check_assumptions in lifelines, and test time‑interaction terms. Consider stratification or time‑varying effects if violated.
  • What is a good C‑index? Depends on context. Around 0.5 indicates no discrimination. Values 0.7–0.8 are often considered useful, but calibration and clinical utility also matter.
  • How do I handle time‑dependent covariates in Cox models? Use long‑format data with start/stop times or extended Cox models; consider joint models or marginal structural models when time‑varying confounding exists.

✅Recommended Resources

Documentation & Python

Clinical Sources

  • Altman DG – Practical Statistics for Medical Research
  • Harrell FE – Regression Modeling Strategies