Cox proportional-hazards regression

The log-rank test in 13.1 answered a yes-or-no question: do two survival curves differ? It does not tell you by how much, and it cannot adjust for age, stage, or anything else that differs between the groups. Cox proportional-hazards regression fills both gaps. It gives the effect of a treatment or risk factor as a single number, the hazard ratio, and it lets you estimate that effect while holding other variables fixed. This lesson builds the model, reads its output, and shows where its one big assumption can break.

We keep the survival notation from 13.1: each patient has a follow-up time and a status that records whether the event happened or the observation was censored. The running example is a cohort of patients with advanced colorectal cancer treated at a Malaysian oncology centre, followed from the start of treatment until death or censoring.

The hazard and the hazard ratio

The survival function from 13.1 looks at the probability of surviving past time . The hazard looks at the same process from a different angle. The hazard is the instantaneous event rate at time among those still at risk: given that a patient has survived to , how fast is the event arriving right now. It is a rate, not a probability, so it can be larger than 1, and it can rise and fall over follow-up.

To compare two groups you take the ratio of their hazards. The hazard ratio is the effect measure of survival analysis, the way the odds ratio is for logistic regression. A hazard ratio of 2 means the event rate in one group is twice the rate in the other at any given moment. A hazard ratio below 1 means a lower rate, so better survival.

The Cox proportional-hazards model

Cox regression models the hazard for an individual as a baseline hazard multiplied by a factor that depends on their covariates:

Here is the baseline hazard, the hazard for a patient whose covariates are all 0, and to are the covariates such as treatment group and age. The model is called semi-parametric because it leaves completely unspecified. You never have to say whether the underlying risk rises, falls, or wobbles over time. Cox regression estimates the coefficients without ever estimating the shape of . That is what makes it the default tool: you get the effect of your covariates without committing to a curve for the baseline risk.

The trick that allows this is the risk set. At each time an event occurs, Cox compares the covariates of the patient who had the event against everyone still being followed at that instant. The baseline hazard is shared by every patient in the risk set, so it cancels out of the comparison and never needs a value.

Reading a coefficient as a log hazard ratio

Take the simplest case: one binary covariate , coded 1 for the new regimen and 0 for standard care. Write the hazard ratio comparing the two groups and the baseline hazard cancels:

So the regression coefficient is the log hazard ratio, and you exponentiate it to get the hazard ratio itself. Software reports the coefficient on the log scale, then as the HR with a 95% confidence interval. The recipe is always the same.

1. Fit the model and read the coefficient for each covariate.
2. Exponentiate it: is the hazard ratio for a one-unit increase in that covariate.
3. Exponentiate the confidence limits of to get the 95% CI for the HR.
4. Check whether the interval excludes 1. If it does, the effect is significant at the 5% level.

Worked example 1: one binary covariate

Thirty-two patients are followed from the start of treatment. The status variable is 1 for death and 0 for censored, exactly as in 13.1. The covariate group is 1 for a new oral regimen and 0 for the standard infusion. We fit a Cox model with this single covariate and read the hazard ratio for the new regimen.

webr::install("survival")
library(survival)
time   <- c(5,8,12,3,22,18,30,7,15,25,10,40,6,28,35,14,9,45,20,33,4,26,11,38,16,50,13,42,19,55,24,17)
status <- c(1,1,1,1,0,1,0,1,1,0,1,0,1,0,0,1,1,0,1,1,1,1,1,0,1,0,1,0,1,0,1,1)
group  <- c(0,1,0,0,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0)
fit <- coxph(Surv(time, status) ~ group)

exp(coef(fit))      # fill in: read the hazard ratio
exp(confint(fit))   # fill in: read the 95% CI

exp(coef(fit))        # hazard ratio
#     group
# 0.2057
exp(confint(fit))     # 95% CI for the HR
#        2.5 %    97.5 %
# group  0.080    0.527
summary(fit)$coefficients   # coef, se, z, p
# group  coef = -1.581,  z = -3.29,  p = 0.001

import numpy as np
from statsmodels.duration.hazard_regression import PHReg
time   = np.array([5,8,12,3,22,18,30,7,15,25,10,40,6,28,35,14,9,45,20,33,4,26,11,38,16,50,13,42,19,55,24,17])
status = np.array([1,1,1,1,0,1,0,1,1,0,1,0,1,0,0,1,1,0,1,1,1,1,1,0,1,0,1,0,1,0,1,1])
group  = np.array([0,1,0,0,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0])
res = PHReg(time, group, status=status).fit()

hr = # fill in: exponentiate the coefficient
ci = # fill in: exponentiate the confidence limits

hr = np.exp(res.params[0])        # hazard ratio
ci = np.exp(res.conf_int()[0])    # 95% CI
print(round(hr, 3), np.round(ci, 3))
# 0.206 [0.08  0.527]
print(round(res.pvalues[0], 4))   # p-value
# 0.001

The hazard ratio is 0.21, with a 95% confidence interval of 0.08 to 0.53. At any moment during follow-up, patients on the new regimen die at about one-fifth the rate of patients on standard care. The interval sits well below 1, so the benefit is statistically clear in this cohort. Notice the link back to 13.1: a Cox model with a single binary covariate is the regression form of the log-rank test. The Cox score test and the log-rank test ask the same null hypothesis, that the two hazards are equal, and they give almost identical p-values. The difference is that Cox also hands you the size of the effect, not just the verdict.

Adjusting for several covariates

The point of regression is adjustment. Patients on the new regimen might be younger, and younger patients survive longer for reasons that have nothing to do with the drug. Put both covariates in one model and each coefficient is the effect of that variable with the other held fixed. The result is an adjusted hazard ratio: the treatment effect among patients of the same age. The model assumes the covariate effects multiply, so a coefficient is read after controlling for the rest.

Worked example 2: an adjusted model with two covariates

We add age in years to the model. Now the treatment hazard ratio is adjusted for age, and the age hazard ratio is adjusted for treatment.

webr::install("survival")
library(survival)
time   <- c(5,8,12,3,22,18,30,7,15,25,10,40,6,28,35,14,9,45,20,33,4,26,11,38,16,50,13,42,19,55,24,17)
status <- c(1,1,1,1,0,1,0,1,1,0,1,0,1,0,0,1,1,0,1,1,1,1,1,0,1,0,1,0,1,0,1,1)
group  <- c(0,1,0,0,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0)
age    <- c(71,65,58,71,49,73,53,73,61,70,55,58,77,65,57,72,62,52,75,50,65,61,65,57,57,48,77,63,74,59,75,74)
fit2 <- coxph(Surv(time, status) ~ group + age)

summary(fit2)$conf.int   # fill in: adjusted HRs and their CIs

summary(fit2)$conf.int
#        exp(coef)  exp(-coef)  lower .95  upper .95
# group     0.251       3.98       0.095      0.660
# age       1.069       0.935      1.012      1.130
cox.zph(fit2)            # check proportional hazards (see below)

import numpy as np, pandas as pd
from statsmodels.duration.hazard_regression import PHReg
time   = np.array([5,8,12,3,22,18,30,7,15,25,10,40,6,28,35,14,9,45,20,33,4,26,11,38,16,50,13,42,19,55,24,17])
status = np.array([1,1,1,1,0,1,0,1,1,0,1,0,1,0,0,1,1,0,1,1,1,1,1,0,1,0,1,0,1,0,1,1])
group  = np.array([0,1,0,0,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0])
age    = np.array([71,65,58,71,49,73,53,73,61,70,55,58,77,65,57,72,62,52,75,50,65,61,65,57,57,48,77,63,74,59,75,74])
exog = pd.DataFrame({"group": group, "age": age})
res2 = PHReg(time, exog, status=status).fit()

hr = # fill in: exponentiate the parameters
ci = # fill in: exponentiate the confidence intervals

hr = np.exp(res2.params)        # adjusted hazard ratios
ci = np.exp(res2.conf_int())    # 95% CIs
print(np.round(hr, 3).tolist())
print(np.round(ci, 3).tolist())
# group: HR 0.251  (0.095, 0.660)
# age:   HR 1.069  (1.012, 1.130)

After adjusting for age, the treatment hazard ratio is 0.25 (95% CI 0.10 to 0.66). The benefit holds and even sharpens a little once age is accounted for. The age hazard ratio is 1.069 (95% CI 1.01 to 1.13), meaning each extra year of age multiplies the death rate by about 1.07, roughly a 7% higher hazard per year, with treatment held fixed. Both intervals exclude 1. The two models tell a consistent story, which is reassuring, but the adjusted estimate is the one to report because it separates the drug effect from the age difference between the arms.

The proportional-hazards assumption

The whole model rests on one assumption, and it is in the name. Proportional hazards means the hazard ratio between two groups is constant over time. The new regimen having one-fifth the death rate must hold at month 3 and at month 40 alike. The two hazards can rise and fall in any pattern, but their ratio stays fixed. That is why a single number can describe the whole follow-up.

The assumption fails when the effect of a covariate changes over time. The clearest failure is hazards that cross. Suppose a surgery carries a high early risk but a long-term benefit, while a drug is gentle early and weaker later. Their survival curves cross, the early hazard ratio is above 1 and the late one below 1, and no single HR can be right for both periods. Reporting one number would average them into something that describes neither.

When hazards are non-proportional, a single Cox hazard ratio is a time-weighted average of effects that genuinely differ across follow-up. It is not wrong arithmetic, it is the wrong summary. A drug that helps for a year then loses effect can show an unremarkable overall HR near 1 that hides both the early benefit and the late catch-up. Always look at whether the assumption holds before you trust one hazard ratio.

Checking proportional hazards

Two standard checks tell you whether the assumption is reasonable.

• Log-log survival plot. Plot against for each group. Under proportional hazards the curves are roughly parallel, separated by a constant gap. Lines that converge, diverge, or cross signal trouble.
• Schoenfeld residuals. These residuals should show no trend against time if the hazard ratio is constant. A clear slope means the effect is drifting over follow-up. In R, cox.zph(fit) tests each covariate and the model as a whole, returning a small p-value when proportionality is in doubt.

If a covariate fails the check, you have three options. You can split follow-up into periods and estimate a separate HR in each. You can add an interaction between the covariate and time. Or you can stratify on the offending variable when it is a nuisance confounder rather than the exposure of interest. Stratification lets the baseline hazard differ across strata without forcing a constant ratio.

Common mistakes

• Reading the coefficient as the hazard ratio. The raw coefficient is the log hazard ratio. A coefficient of -1.58 is not an HR of -1.58; exponentiate it to get 0.21. Always report , never the bare coefficient.
• Treating the hazard ratio as a risk ratio or a cure rate. An HR of 0.21 does not mean 79% of patients are cured or that absolute risk drops by 79%. It compares instantaneous event rates. The effect on five-year survival depends on the baseline hazard too.
• Never checking proportional hazards. A single HR is only meaningful if the assumption holds. Crossing or converging curves make it misleading. Run cox.zph or a log-log plot before you trust the number.
• Adjusting for a variable on the causal pathway. Putting a mediator in the model removes part of the effect you want to measure. Adjust for confounders, not for steps between the exposure and the outcome.
• Forgetting that censoring must be non-informative. Cox inherits the same assumption as the Kaplan-Meier estimator from 13.1. If patients drop out because they are doing badly, the hazard ratios are biased.

Tips

• Report the hazard ratio, its 95% confidence interval, and the p-value together. The interval carries more information than the p-value alone.
• Centre or scale continuous covariates for readable output. A hazard ratio per 10-year increase in age often reads better than the per-year figure, and you get it by dividing age by 10 before fitting.
• Fit the unadjusted model first, then the adjusted one. A large shift in a hazard ratio after adjustment points to real confounding worth explaining.
• With no tied event times, R's coxph and Python's PHReg give the same estimates. When ties are common, check which tie-handling method each uses, since the defaults differ.
• A Cox model with one binary covariate is the regression twin of the log-rank test from 13.1. Use the test for a quick comparison, the model when you need the effect size or adjustment.