AdvancedIndependently reviewed, not yet spot-checked by a human

The Cox proportional hazards model

What ratio a hazard ratio actually is (not a ratio of survival times, and not an odds ratio), why the coefficients are estimable without ever specifying the baseline hazard, what multivariable adjustment is really doing, how ties are handled, and how to write up a confidence interval that crosses 1.

What problem this model solves

Kaplan-Meier and the log-rank test can tell you that two curves differ, but they cannot tell you two things: by how much, and by how much once everything else is accounted for.

Real clinical questions are almost always the second kind. Women with lung cancer survive longer than men — is that the effect of age? Of performance status? Or of sex itself? And if the grouping variable is continuous (age, tumour size, eGFR), the log-rank test has nothing to work with at all; you would have to chop it into categories, at cut points you chose yourself.

The Cox proportional hazards model handles exactly this: put several variables in at once, continuous or categorical, get an effect estimate for each, and do it without assuming survival times follow any particular distribution.

The core idea: the hazard function and the HR

Start with the hazard function $h(t)$ : for someone who has not yet had the event by time $t$ , the rate at which the event occurs at that instant.

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

Two things in that expression matter: the condition $T \ge t$ in the numerator (only people still in the risk set count), and the fact that it is a rate, not a probability, so it can exceed 1.

The Cox model is written:

h(t \mid x) = h_0(t) \cdot \exp(\beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p)

which is two pieces multiplied together:

$h_0(t)$ is the baseline hazard — how risk moves over time for a person whose covariates are all zero. It can take any shape at all: high then low, peaked, wavy, anything.
$\exp(\beta_1 x_1 + \cdots)$ is a multiplier that does not change with time, pulling that entire baseline curve up or down.

So the ratio of the hazards of two people is:

\frac{h(t \mid x = 1)}{h(t \mid x = 0)} = \exp(\beta)

$h_0(t)$ cancels. That quantity is the hazard ratio (HR), and the cancellation explains two things at once: why $h_0(t)$ never has to be specified (it disappears from the ratio — this is what semi-parametric means: the effect part is parameterised, the time part is not), and why the model is called proportional hazards — that ratio is assumed to be the same number throughout follow-up.

Run it yourself

library(survival)
data(cancer, package = "survival")   # lung, colon and mgus2 all live here

lung$sex_f <- factor(lung$sex, levels = c(1, 2),
                     labels = c("Male", "Female"))

fit <- coxph(Surv(time, status) ~ age + sex_f + ph.ecog + wt.loss,
             data = lung)
summary(fit)        # coef, exp(coef) = HR, 95% CI, p, concordance

# Age rescaled to "per 10 years": multiply on the coefficient scale,
# not by multiplying the HR by 10
b  <- coef(fit)["age"]; se <- sqrt(vcov(fit)["age", "age"])
exp(10 * c(b, b - 1.96 * se, b + 1.96 * se))

# A larger real trial: colon adjuvant chemotherapy, etype 2 = death
colon_d <- subset(colon, etype == 2)
colon_d$rx <- factor(colon_d$rx, levels = c("Obs", "Lev", "Lev+5FU"))
fit2 <- coxph(Surv(time, status) ~ rx + sex + age + nodes + extent + differ,
              data = colon_d)
summary(fit2)

# Adjusted survival curves: all three arms evaluated at the same covariates
ref <- data.frame(rx = factor(levels(colon_d$rx), levels = levels(colon_d$rx)),
                  sex = 1, age = 60, nodes = 2, extent = 3, differ = 2)
plot(survfit(fit2, newdata = ref), lwd = 2,
     col = c("#8c6d4a", "#4d6a8c", "#c44d4d"),
     xlab = "Days", ylab = "Adjusted survival")

Verified with R 4.6.0 and survival 3.8.6

import statsmodels.api as sm
from lifelines import CoxPHFitter

# Naming trap: on Rdatasets, lung sits under the "cancer" help page
lung = sm.datasets.get_rdataset("cancer", "survival").data
d = lung[["time", "status", "age", "sex", "ph.ecog", "wt.loss"]].dropna().copy()
d["event"]  = (d["status"] == 2).astype(int)
d["female"] = (d["sex"] == 2).astype(int)
d = d.drop(columns=["status", "sex"])

cph = CoxPHFitter().fit(d, duration_col="time", event_col="event")
cph.print_summary()          # coef, exp(coef), CI, p, concordance
cph.plot()                   # this is a forest plot

# colon and mgus2 are not on Rdatasets (the same help-page naming problem),
# so to use them from Python, export them from R first:
#   write.csv(subset(colon, etype == 2), "colon_death.csv", row.names = FALSE)

lifelines' CoxPHFitter also defaults to Efron for ties, and its coefficients agree with R to several decimal places.

Forest plot of a multivariable Cox model fitted to survival::lung, showing hazard ratios and 95% confidence intervals for four variables. The intervals for female sex and ECOG performance status exclude 1, while those for age and weight loss cross 1. — The four-variable Cox model on survival::lung. The dashed line is HR = 1. An interval crossing that line means this study did not detect an association between that variable and the hazard of death.Plotting script `figures/scripts/B3-03-cox.R`

How to read the report

What summary(coxph(...)) prints maps onto every column of a paper’s table. Walking through the lung model above cell by cell — but first, how many patients the model was actually fitted on: the dataset has 228 patients, 15 of whom were dropped by complete-case analysis because a covariate was missing, so the model was fitted on 213 patients with 151 events.

(Say so when patients drop out like this. coxph() issues no warning; it quietly fits on fewer rows. And missingness is often related to prognosis, so this is not a harmless rounding detail. The cohort study chapter has a whole section on it.)

Variable	HR	95% CI	p
Age, per year	1.01	0.99–1.03	0.165
Female vs male	0.55	0.39–0.78	< 0.001
ECOG PS, per point	1.67	1.31–2.14	< 0.001
Weight loss, per kg	0.99	0.98–1.00	0.176

Column by column:

The direction of the HR. Above 1 is higher hazard, below 1 is lower. Women versus men gives HR = 0.55, which reads as “among patients of the same age, the same ECOG performance status and the same weight loss, a woman’s hazard of death at any instant is about 0.55 times a man’s”. “Holding the other variables constant” is not a politeness — it is part of the definition of the HR; remove it and you are quoting a different number.
The unit of a continuous variable. The HR for age is 1.013, which looks too small to matter — because it is the effect of one additional year. Per 10 years it becomes HR 1.14 (0.95–1.38). The rescaling is done on the coefficient scale, $\exp(10\beta)$ ; you cannot multiply the HR by 10. An HR for a continuous variable whose unit is not stated cannot be interpreted at all.
The confidence interval is more useful than the p-value. The interval for age is 0.995–1.033, which crosses 1; for weight loss it is 0.978–1.004, which also crosses 1.
Concordance (Harrell’s C) = 0.647. Pick two patients at random: this is the probability that the model gives the higher risk to the one who has the event first. 0.5 is a coin flip, 1 is perfect. 0.647 means the model carries signal but discriminates only moderately — which is entirely ordinary for clinical prognostic models, and is why “the p-value is highly significant” and “the model can predict an individual patient” are two separate questions.
The overall model test (likelihood ratio test) = 31.0, df = 4, p < 0.001. It asks whether all the coefficients are simultaneously zero — the survival analogue of the overall F test in linear regression.

How to write up a confidence interval that crosses 1. The HR for age is 1.013 (0.995–1.033), p = 0.165. That should be written as “this study did not detect a significant association between age and the hazard of death”, and must not be written as any of the following:

❌ “Age is unrelated to survival” — failing to detect is not the same as absent, and an upper bound of 1.033 leaves ample room for a clinically meaningful effect
❌ “Each additional year of age raises the hazard of death by 1.3%” — reporting a point estimate whose interval crosses 1 as though it were an effect is the most common form of over-interpretation
❌ “The two are equivalent” — claiming equivalence requires a non-inferiority or equivalence design with a margin fixed in advance

The correct sentence carries the interval with it: “not statistically significant, with a confidence interval of 0.995–1.033, which is not narrow enough to rule out a moderate effect.”

What multivariable adjustment is actually doing

“Adjusting for age” sounds like subtracting the influence of age, but what happens inside the model is more concrete: at every event time, each person in the risk set is weighted by the risk their own covariates predict, and the model asks whether the person who actually had the event was the one predicted to be at higher risk. The coefficients are the set of numbers that makes this as consistent as possible across all event times — technically, the partial likelihood.

The consequence is that the reported HR for sex is a comparison between two people of the same age, the same ECOG score and the same weight loss. Three practical implications follow:

Change the adjustment set and the HR changes. Two papers reporting different HRs may simply have fitted different models, not contradicted each other. Never read an HR without reading what was in the model with it.
Do not put a mediator in. If the treatment works by lowering an inflammatory marker, adjusting for that marker subtracts part of the treatment effect you were trying to estimate. What belongs in the model is decided by a causal diagram (DAG), not by p-values.
The number of events limits how many variables you can afford. The traditional rule of thumb is at least 10 events per variable (EPV); recent simulation work argues it can be relaxed depending on circumstances, but the direction is unchanged: with 151 events, four variables is comfortable and fifteen would be fitting noise.

For a look at adjustment on a real trial, survival::colon is an adjuvant chemotherapy trial after colon cancer surgery. The file ships as a two-endpoint long format — one recurrence row and one death row per patient, 1858 rows in all. This page first takes the 929 death-endpoint rows with etype == 2, and then, after dropping covariate missingness, the model uses 888 patients and 430 deaths:

Adjusted survival curves for the three arms of the colon trial. All three curves are evaluated at the same covariate values (same age, same node count, same extent), so the differences between them come only from treatment. This is not the same as plotting KM curves — KM shows the patients each arm actually received.Plotting script figures/scripts/B3-03-cox.R

Variable	HR	95% CI	p
Levamisole vs observation	0.913	0.731–1.141	0.423
Levamisole+5FU vs observation	0.673	0.530–0.854	0.001
Male sex	0.967	0.799–1.169	0.728
Age, per year	1.006	0.998–1.014	0.131
Positive nodes, per node	1.091	1.071–1.111	< 0.001
Local extent, per grade	1.617	1.293–2.023	< 0.001
Differentiation, per grade	1.157	0.949–1.409	0.149

This table turns the warning above into a concrete case: levamisole alone has an HR of 0.913 with an interval of 0.731–1.141 that crosses 1 — this trial did not detect an association between levamisole alone and the hazard of death. Levamisole plus 5-FU has an HR of 0.673 (0.530–0.854), an interval lying entirely to the left of 1. Historically this trial is one of the pieces of evidence that made 5-FU-based combinations the standard adjuvant therapy for colon cancer.

Ties: two patients with the event on the same day

The derivation of the Cox model assumes event times are continuous — that no two people have the event at exactly the same instant. Real data of course do: record time in days, months or years and ties are guaranteed. R offers three ways to handle them:

Method	What it does	When to use it
Efron	A weighted approximation over the tied events	R’s default; choose it essentially always
Breslow	A cruder approximation that treats tied events as independent	SAS’s default; with many ties it pulls coefficients toward zero
exact	Enumerates every possible ordering of the tied events	Theoretically the most correct, and very slow when ties are abundant

This section switches back to the full colon death-endpoint data (no complete-case exclusion, because comparing tie handling involves only time and treatment arm and does not need the covariates with missingness), so the event count is 452 rather than the 430 in the model table above — two different numbers on one page because the analysis sets differ, not because one of them is wrong.

With time recorded in days, the 452 deaths fall on 409 distinct days, and the three methods barely differ:

Method	HR for Lev+5FU	HR for node count
efron	0.6709	1.0957
breslow	0.6710	1.0956
exact	0.6709	1.0958

Coarsen the same data so that time is measured in years (452 events crammed into 8 time points) and the gap opens up:

Method	HR for Lev+5FU	HR for node count
efron	0.6747	1.0926
breslow	0.6955	1.0853
exact	0.6570	1.1210

Coarsening to years was done deliberately as a demonstration, not as analytical advice. But the situation it demonstrates is entirely real: national health insurance claims databases, cancer registries and long-term follow-up studies reporting in months or years produce tie densities exactly like this. Under those conditions Breslow pulls the HR toward 1 (0.695 against Efron’s 0.675) — that is, it conservatively understates the effect.

Both variables are shown so you can confirm this is not a cherry-picked column: the node count behaves the same way (Breslow 1.085 is closer to 1 than Efron’s 1.093). A pattern that holds for one variable only becomes a general claim once you check a second one — with only one column reported, the reader has no way to tell.

The practical conclusion is short: use Efron, R’s default, and do not switch to Breslow yourself. When reading older papers produced in SAS, just be aware that its default is Breslow.

Assumptions to check before you rely on it

The Cox model has three assumptions, and the first is both the easiest to skip and the most damaging to get wrong:

Proportional hazards — the HR is the same number throughout follow-up. When it is violated, the single HR is an average of effects that pointed in different directions at different times and has no clear meaning. How to check it and what to do instead are on the proportional hazards assumption and Schoenfeld residuals.
Continuous variables are linear in the log hazard — age need not act on the hazard of death in a straight line. Check by plotting martingale residuals against the variable, or sidestep the issue with restricted cubic splines.
Non-informative censoring — the foundation shared by every survival method; see censoring and the structure of survival data.

One more: when competing risks are present, a cause-specific Cox model answers a different question from Fine-Gray, which is the subject of competing risks: CIF and Fine-Gray.

Common misuses

Misuse	Why it is wrong
Reading an HR as a ratio of survival times	An HR compares instantaneous hazards; a ratio of survival times is what an AFT model estimates
Reading an HR as an odds ratio or a risk ratio	All three have different denominators and different relationships to time; close in value when events are rare, different in meaning always
Reporting an HR with no baseline risk or absolute risk difference	HR 0.5 on a rare event may buy a very small absolute benefit
Reporting an HR for a continuous variable without stating the unit	Per year, per decade and per standard deviation are different numbers, and the reader cannot interpret any of them
Converting an HR to “per 10 units” by multiplying by 10	The rescaling happens on the coefficient scale: $\exp(10\beta)$ is not $10 \exp(\beta)$
Describing the size of an effect whose confidence interval crosses 1	Write “did not detect a significant association”, and give the interval to show the range of uncertainty
Claiming “the two are equivalent” from a non-significant result	Equivalence needs a non-inferiority or equivalence design with a pre-specified margin
Fitting many covariates on few events	Overfitting; neither the coefficients nor the confidence intervals can be trusted
Adjusting away a mediator of the treatment effect	It subtracts part of the effect you set out to estimate; what to include is decided by a causal diagram, not by p-values
Fitting a Cox model without checking proportional hazards	If the assumption fails, that HR has no clear interpretation
Grouping by a status only knowable after baseline	Immortal time bias; see B3-01
Comparing HRs between two papers directly	Different adjustment sets mean the two HRs define different contrasts

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B3-03-cox.R

Read the figure

The answer comes from the same statistical output that produced this page's figures, not from a number typed in beside them.

This Cox model has four covariates, and the summary prints a hazard ratio on each row. Which row is ECOG PS?

Show the answer and why

Correct answer: Hazard ratio 1.67 - each further ECOG point raises the hazard of death by about two thirds

A higher ECOG score means worse performance status, so its hazard ratio has to exceed one - which rules out the other two options before reading any output at all. 0.55 is in fact female versus male, and 0.99 is per kilogram of weight loss. The three sit next to each other in the summary, and the order R prints them in follows the formula rather than your sense of which matters. Checking the label on the row before reading the hazard ratio off it is the cheapest step in reading regression output, and the one most often skipped.

The n= line in the output is smaller than the number of rows in the dataset. Which statement is right?

Show the answer and why

Correct answer: The model used 213 patients - any row with a missing value on any covariate was dropped whole

coxph does complete-case analysis by default: one missing covariate removes the entire row, so 213 of the 228 rows went in. The output says so only in the single n= / number of events= line tucked under the coefficient table, which is easy to slide past; 151 is the event count, not a patient count. Those fifteen patients did not disappear at random - they tend to be the sickest and the least completely measured, so this is not just a smaller sample but a different one.

For the same model's likelihood ratio test, what determines the degrees of freedom?

Show the answer and why

Correct answer: Degrees of freedom are 4, the number of covariates the model estimates

The likelihood ratio test compares the model with these covariates against one with none, and the degrees of freedom are the difference in parameters - four covariates here. 15 is how many rows were dropped for missingness and 151 is the event count; those integers sit within a few lines of each other in the output and mean nothing like the same thing. The event count does govern how much power the test has, but it is not its degrees of freedom - and the test asks whether the four terms together explain anything, not whether any one of them is significant.

Chapters that use this method

Watch next

ENThe Pharmacist Academy· 3 minExplains what a hazard ratio is in under three minutes. Watch it before the second section and the algebra will go down easily.

ENnumiqo· 6 minA six-minute tour of the model, with animation that makes "the hazard ratio stays constant" visible.

ENBiostatsquid· 11 minLines up the coefficients, the hazard ratios and the output columns cell by cell — exactly the "how to read the report" section below.

繁中Ming-Chieh Shih· 18 minIn Mandarin, continuing the first part recommended on B3-01. This episode covers the Cox model and the hazard ratio.

Sources and licences

This page is original writing