Cox Proportional Hazards Model Tutorial

Introduction

The Cox proportional hazards model, introduced by Sir David Cox in a single 1972 paper, is the dominant regression method for time-to-event data in clinical research — the model behind most reported hazard ratios in oncology, cardiovascular, and virtually every other therapeutic area with a survival or time-to-progression endpoint. Its central innovation, the partial likelihood, made it possible to estimate covariate effects on the hazard without ever specifying the shape of the underlying baseline hazard function — a genuinely unusual and powerful idea in statistics.

This tutorial goes deep on the full mechanics: the model and its partial likelihood derived from first principles, a real Newton-Raphson fit on simulated trial data, the exact relationship between the log-rank test and the Cox score test, how ties are handled (Breslow vs. Efron, with a real numeric comparison), how to test the proportional hazards assumption with Schoenfeld residuals, and — deliberately given equal weight — the model's real shortcomings, illustrated with a genuine non-proportional-hazards example.

The Hazard Function and Model Specification

The hazard function $$h(t)$$ is the instantaneous event rate at time $$t$$ , conditional on survival up to $$t$$ :

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

The Cox model specifies how covariates $X = (X_1,\dots,X_p)$ multiplicatively scale a common baseline hazard:

h(t \mid X) = h_0(t)\,\exp(\beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p)

Two things make this specification distinctive. First, $$h_0(t)$$ — the hazard for a subject with all covariates equal to zero — is left completely unspecified; it can take any shape over time. Second, and this is the defining proportional hazards assumption: the ratio of hazards between any two subjects is constant over time, since it depends only on the difference in their linear predictors, not on $$t$$ :

\frac{h(t\mid X_A)}{h(t\mid X_B)} = \exp\big(\beta^\top(X_A - X_B)\big) \quad \text{— constant in } t

The Partial Likelihood — Cox's Key Insight

Because $$h_0(t)$$ is unspecified, the ordinary full likelihood can't be written down or maximized directly. Cox's solution: condition, at each observed event time $t_{(i)}$ , on which subject among those still at risk had the event — not on when it happened. This conditional probability doesn't involve $$h_0(t)$$ at all, because it cancels out of the ratio:

P(\text{subject } i \text{ has the event} \mid \text{one event occurs at } t_{(i)}, \text{ risk set } R(t_{(i)})) = \frac{\exp(\beta^\top X_i)}{\sum_{j \in R(t_{(i)})} \exp(\beta^\top X_j)}

Multiplying this conditional probability across every observed event time gives the partial likelihood:

L(\beta) = \prod_{i:\,\delta_i=1} \frac{\exp(\beta^\top X_i)}{\sum_{j \in R(t_i)} \exp(\beta^\top X_j)}

where $\delta_i=1$ indicates subject $$i$$ had an observed event (not censored), and $$R(t_i)$$ is the risk set — everyone still under observation and event-free just before $$t_i$$ . This is not a true likelihood in the usual sense (it isn't derived from a full joint density), but Cox showed it can be maximized and yields consistent, asymptotically normal estimates exactly as if it were one — a result that took several years and additional theoretical work (notably by Tsiatis, and Andersen & Gill using counting-process theory) to fully justify.

Estimation via Newton-Raphson: A Real Worked Fit

The log partial likelihood, its score (first derivative), and observed information (negative second derivative) for a single covariate are:

\ell(\beta) = \sum_{i:\delta_i=1}\left[\beta X_i - \log\!\left(\sum_{j\in R(t_i)} e^{\beta X_j}\right)\right]$$ $$U(\beta) = \sum_{i:\delta_i=1}\left[X_i - \frac{S_1(\beta,t_i)}{S_0(\beta,t_i)}\right], \qquad S_0=\sum_{j\in R} e^{\beta X_j},\ \ S_1=\sum_{j \in R} X_j e^{\beta X_j}

Newton-Raphson updates $\beta \leftarrow \beta + U(\beta)/I(\beta)$ until the score is (numerically) zero. Below is the actual iteration history fitting a univariate Cox model (treatment indicator only) on a simulated 200-subject oncology trial (124 events, 38% censored) — not a textbook illustration, the real numeric path from a from-scratch implementation of the formulas above:

Iteration	β	exp(β)	Score U(β)	Information I(β)	Log-Likelihood
0	0.000000	1.000000	−9.138358	29.578707	−561.9678
1	−0.308951	0.734217	0.124237	30.164512	−560.5707
2	−0.304832	0.737247	−0.000002	30.165392	−560.5704
3 (converged)	−0.304832	0.737247	0.000000	30.165391	−560.5704

Convergence in three iterations is typical — the partial log-likelihood is well-behaved (globally concave for a single covariate), so Newton-Raphson converges quadratically once near the optimum. From the converged Information, $SE(\hat\beta) = 1/\sqrt{I(\hat\beta)} = 0.18207$ , giving a Wald 95% CI for the hazard ratio of (0.516, 1.053) and $$p=0.094$$ . This matches, to four decimal places, what a production implementation (the lifelines Python package) returns on the identical data — confirming the from-scratch derivation above is correct, not just illustrative.

The Multivariate Model: A Real Trial Example

Extending to three covariates — treatment assignment, age, and ECOG performance status (poor vs. good) — on the same simulated trial:

Covariate	β̂	HR = exp(β̂)	SE(β̂)	95% CI for HR	p-value
Treatment (vs. control)	−0.2983	0.742	0.1823	(0.519, 1.061)	0.102
Age (per year)	0.0125	1.013	0.0103	(0.992, 1.033)	0.227
ECOG PS ≥1 (vs. 0)	0.5866	1.798	0.1879	(1.244, 2.598)	0.0018

Hazard ratios and 95% confidence intervals from the multivariate model above. The vertical line at HR=1 marks no effect; ECOG performance status is the only covariate whose interval excludes 1 in this particular sample.

Model fit: concordance index 0.600 (a rank-based measure similar to AUC, adapted for censored survival data), overall likelihood ratio test $\chi^2=14.75$ on 3 df, $p \lt 0.005$ — strong evidence the model as a whole predicts survival better than an intercept-only model, even though the treatment effect alone falls short of conventional significance in this particular (deliberately modest) sample size.

Interpreting the Hazard Ratio — Carefully

HR = 0.742 for treatment means: at any instant, among subjects still event-free and matched on the other covariates, the instantaneous event rate in the treatment arm is 74.2% of the control arm's rate. Several things this does not mean, that are worth stating explicitly because they're common sources of misinterpretation:

It is not a relative risk. Relative risk compares cumulative probabilities over a fixed period; the hazard ratio compares instantaneous rates and is mathematically distinct, though the two are sometimes numerically close when events are rare.
It does not, by itself, tell you how long a "typical" treated patient survives — that requires converting to the survival scale (e.g., via the baseline hazard, below), since the same HR can correspond to very different absolute survival gains depending on the baseline event rate.
It assumes the ratio truly is constant over time — if it isn't, the fitted HR is a peculiar (event-count-weighted) average of a time-varying effect, and can be misleading. This is explored fully in the Shortcomings section below.

Three Ways to Test a Cox Coefficient

Every coefficient in a Cox model can be tested three (asymptotically equivalent, but numerically distinct in finite samples) ways:

Test	Statistic	Requires refitting?
Wald test	$Z = \hat\beta/SE(\hat\beta)$ , compared to standard normal	No — uses only the fitted model
Score test	$$U(0)^2/I(0)$$ , evaluated at the null $\beta=0$	No — needs only the null model's score/information
Likelihood Ratio test	$2[\ell(\hat\beta) - \ell(0)]$ , compared to $\chi^2_1$	Yes — needs both fitted and null log-likelihoods

On the univariate treatment model above: Wald $$Z=-1.674$$ ( $$p=0.094$$ ); the likelihood ratio test gives $\chi^2 = 2.80$ ( $$p=0.09$$ ). In small-to-moderate samples the three tests can disagree modestly; the likelihood ratio test is generally considered the most reliable of the three in smaller samples, since the Wald test's accuracy depends on the log-likelihood being well-approximated by a quadratic near $\hat\beta$ , an approximation that degrades for large effect sizes or sparse data.

The Log-Rank Test — and Its Exact Connection to the Cox Score Test

The log-rank test predates the Cox model but turns out to be a special case of it: the log-rank test is algebraically identical to the Cox score test for a single binary covariate with no tied event times. This isn't a loose analogy — it's the same formula. The log-rank statistic compares, at every observed event time, the number of events actually observed in each group to the number expected under the null of equal hazards, using a hypergeometric variance:

\chi^2 = \frac{(O_1 - E_1)^2}{V}, \qquad E_{1,i} = d_i \frac{n_{1,i}}{n_i}, \qquad V = \sum_i d_i\frac{n_{1,i}n_{0,i}(n_i-d_i)}{n_i^2(n_i-1)}

Hand-Worked Example

A 14-subject subsample, sorted by time, illustrates the full calculation (n = at risk, n1/n0 = at risk by arm, d = events at that time, d1 = events in treatment arm, E1 = expected events in treatment arm):

t	n	n1	n0	d	d1	E1	Var
0.06	14	10	4	1	0	0.7143	0.20408
0.78	12	9	3	1	1	0.7500	0.18750
4.59	11	8	3	1	1	0.7273	0.19835
7.14	9	7	2	1	1	0.7778	0.17284
10.26	8	6	2	1	1	0.7500	0.18750
14.32	7	5	2	1	0	0.7143	0.20408
15.04	6	5	1	1	0	0.8333	0.13889
15.52	5	5	0	1	1	1.0000	0.00000
16.66	4	4	0	1	1	1.0000	0.00000
21.42	2	2	0	1	1	1.0000	0.00000

Summing: $$O_1 = 7$$ , $$E_1 = 8.267$$ , $$V = 1.293$$ , giving $\chi^2 = (7-8.267)^2/1.293 = 1.241$ ( $$p=0.265$$ , 1 df) — small subsample, so a wide, inconclusive interval is expected.

Validating the Equivalence on the Full Trial

On the complete 200-subject dataset, the log-rank test gives $\chi^2 = 2.826$ ( $$p=0.093$$ ). The univariate Cox score/likelihood-ratio test on the identical data gave $\chi^2=2.80$ ( $$p=0.09$$ ) — the same result to within rounding, confirming the theoretical equivalence directly rather than merely asserting it.

Handling Tied Event Times: Breslow vs. Efron

The partial likelihood as written assumes no two events happen at exactly the same time — a reasonable assumption for continuously measured time, but real trial data are often recorded to the nearest day or month, producing genuine ties. Two standard approximations extend the partial likelihood to handle them:

Breslow's approximation treats tied events as occurring in an arbitrary order and uses the full risk-set sum for each tied event, an approach that is simple but biased toward the null when there are many ties:

L_{\text{Breslow}}(\beta) = \prod_i \frac{\exp(\beta^\top s_i)}{\left(\sum_{j\in R(t_i)} e^{\beta^\top X_j}\right)^{d_i}}

Efron's approximation (the standard default in modern software, including lifelines and R's survival package) is more accurate: it progressively removes a fraction of the tied cases' contribution to the risk set as it works through them, better approximating the true discrete partial likelihood.

Rounding the same trial's event times to the nearest month manufactures substantial ties (up to 15 tied events at a single month). Refitting the univariate treatment model with each method, from scratch:

Method	β̂	HR	SE
Breslow	−0.2885	0.749	0.1820
Efron	−0.2963	0.744	0.1820

The difference is modest here but grows with the number and size of ties — with heavy tying (e.g., grouped/interval data with many events per interval), the two methods can diverge meaningfully, which is why Efron (or, for heavily grouped data, a fully discrete logistic model) is generally preferred over Breslow in modern practice.

Checking the Proportional Hazards Assumption: Schoenfeld Residuals

The proportional hazards assumption is checkable, not just assumable. Schoenfeld residuals — one per covariate per observed event — measure, at each event time, the difference between the covariate value of the subject who had the event and the risk-set-weighted average covariate value at that time. Under proportional hazards, these residuals should fluctuate randomly around zero with no trend over time; a systematic trend indicates the covariate's effect is changing over time, i.e., a PH violation. The standard test regresses the scaled residuals against (a transformation of) event time and tests the slope against zero.

A Well-Behaved Example

Kaplan-Meier curves by treatment arm from the multivariate trial example above — the two curves separate early and the gap stays roughly proportional throughout, visually consistent with proportional hazards.

Schoenfeld residuals for the treatment covariate, plotted against event time. No systematic trend is visible — consistent with the formal test result: $\chi^2=1.184$ , $$p=0.277$$ (age: $$p=0.744$$ ; ECOG: $$p=0.837$$ ) — no evidence against proportional hazards for any covariate in this model.

Key Assumptions, Stated Explicitly

Proportional hazards — the covariate effect is constant over the entire follow-up period (checkable via Schoenfeld residuals, above).
Non-informative censoring — the reason a subject is censored is unrelated to their underlying risk of the event, conditional on the covariates in the model.
Log-linearity for continuous covariates — the model as specified assumes a linear relationship between each continuous covariate and the log-hazard; a nonlinear true relationship (e.g., a U-shaped age effect) will be poorly approximated unless splines or categorization are used.
Independence between subjects — violated by clustered data (e.g., multiple observations per patient, or patients clustered within sites with shared unmeasured risk), which requires a frailty (random-effects) extension or robust/cluster-adjusted standard errors.

Shortcomings and Limitations

The Proportional Hazards Assumption Is Often Wrong in Practice

This is the model's best-known limitation, and it's worth demonstrating rather than just stating. Consider a genuinely different simulated trial where treatment causes early harm followed by a later benefit — a realistic pattern for, for example, some immunotherapies (early immune-related toxicity, later durable response):

Kaplan-Meier curves for a simulated trial with genuine crossing hazards: the treatment arm has worse survival than control for roughly the first 15 months, then better survival afterward — the curves visibly cross.

Fitting a standard (single, time-fixed coefficient) Cox model to this data gives HR = 1.146 ( $$p=0.43$$ ) — suggesting no meaningful treatment effect at all. This is wrong, and dangerously so: it doesn't mean there's no effect, it means the model's core assumption doesn't hold, and the single fitted HR is a nearly meaningless average of two opposing effects that happen to roughly cancel out.

Schoenfeld residuals for the same crossing-hazards data show a clear downward trend over time — formally confirmed by the test: $\chi^2=9.373$ , $$p=0.0022$$ , strong evidence against proportional hazards. This is exactly the pattern that a single naive HR conceals and the Schoenfeld test is designed to catch.

The practical fix is not to abandon the Cox framework but to extend it: split follow-up time into pre- and post-crossover intervals and fit a time-varying treatment coefficient (or an explicit treatment×time interaction term), which recovers the true early-harm/late-benefit story that the single-coefficient model hides entirely.

The Hazard Ratio Is Not Collapsible

Unlike a risk difference, the hazard ratio is non-collapsible: adjusting for a covariate that is prognostic (predicts the outcome) but not a confounder (unrelated to treatment assignment) can still change the estimated hazard ratio for treatment, even in a perfectly randomized trial with no confounding whatsoever. This is a mathematical property of the hazard ratio itself, not a sign of model misspecification, and it means "the adjusted HR differs from the unadjusted HR" cannot automatically be interpreted as evidence of confounding the way it typically can for a risk difference or (to a lesser degree) an odds ratio.

No Direct Estimate of Absolute Risk

The partial likelihood's entire trick is to cancel out $$h_0(t)$$ — which means the fitted model, by itself, gives no absolute survival probabilities. Getting those requires a second step: the Breslow estimator for the cumulative baseline hazard,

\hat H_0(t) = \sum_{t_i \le t} \frac{d_i}{\sum_{j \in R(t_i)} e^{\hat\beta^\top X_j}}

From the multivariate model above, the estimated baseline cumulative hazard at selected times: $\hat H_0(3)=0.222$ , $\hat H_0(6)=0.451$ , $\hat H_0(12)=0.866$ , $\hat H_0(18)=1.229$ , $\hat H_0(24)=1.736$ , from which subject-specific survival curves $\hat S(t\mid X) = \exp(-\hat H_0(t)e^{\hat\beta^\top X})$ can finally be constructed — but this is an extra, separate estimation step, not something the Cox coefficients alone provide.

Competing Risks Are Treated as Ordinary Censoring

A standard Cox model treats death from any competing cause (e.g., cardiovascular death when the endpoint is cancer progression) as if it were ordinary, non-informative censoring — the subject simply leaves the risk set. But a competing death is not "missing data" the way a lost-to-follow-up subject is: it's a real, informative outcome that removes the subject from ever experiencing the event of interest. This can bias interpretation, especially in older or sicker populations with meaningful competing mortality; the Fine-Gray subdistribution hazard model is the standard alternative when competing risks are a genuine concern (a topic for a dedicated tutorial).

Sensitivity to Omitted Covariates and Informative Censoring

Because the hazard ratio is non-collapsible (above), omitting a prognostic covariate changes the estimand being estimated, not just its precision — unlike ordinary linear regression, where an omitted, non-confounding covariate mainly costs precision. And, as with any survival method, if censoring is informative (related to unmeasured risk) even after conditioning on the model's covariates, all downstream estimates — hazard ratios, baseline hazard, survival curves — are biased, with no fully general diagnostic that detects this from the observed data alone.

Extensions Worth Knowing About

Stratified Cox model — allows a different, unspecified baseline hazard per stratum (e.g., per study site) while sharing covariate effects, useful when a stratifying variable violates PH but isn't of direct interest.
Time-varying covariates and coefficients — the direct fix for the crossing-hazards example above; the covariate (or its coefficient) is allowed to change value over follow-up.
Frailty models — add a random effect to handle clustered or recurrent-event data, relaxing the independence assumption.
Fine-Gray subdistribution hazard model — the standard competing-risks alternative referenced above.

Summary

The Cox model's partial likelihood is a genuinely elegant piece of statistical machinery — it estimates covariate effects on the hazard without ever specifying the baseline hazard's shape, and its score test for a single binary covariate is exactly the familiar log-rank test, not merely similar to it. But every part of that elegance rests on the proportional hazards assumption, which is checkable (Schoenfeld residuals) and, in a meaningful share of real trials, false — and when it's false, a single fitted hazard ratio can actively conceal a clinically important, time-varying treatment effect rather than merely losing some precision. Treat the model's convenience and its assumption's fragility as equally important parts of the same tool.

References

Cox, D.R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society, Series B.
Breslow, N.E. (1974). Covariance analysis of censored survival data. Biometrics.
Efron, B. (1977). The efficiency of Cox's likelihood function for censored data. Journal of the American Statistical Association.
Schoenfeld, D. (1982). Partial residuals for the proportional hazards regression model. Biometrika.
Andersen, P.K. & Gill, R.D. (1982). Cox's regression model for counting processes: a large sample study. Annals of Statistics.
Fine, J.P. & Gray, R.J. (1999). A proportional hazards model for the subdistribution of a competing risk. Journal of the American Statistical Association.

Survival Analysis: Cox Proportional Hazards