Chapter 12 · Regression models
Which factors change the time to an event?
Hazard ratios, the proportional hazards assumption, adjusted survival and competing risks.
Read first: Risk and survival
The problem first
You are a final-year Surgical Oncology resident. Your thesis follows 400 patients for up to 5 years after surgery for colorectal cancer. About half of them (203) received adjuvant chemotherapy. The Kaplan-Meier curves separate well and the log-rank test gives p < 0.001, so you write that chemotherapy reduced relapse and death by more than half.
At the pre-submission seminar, your guide looks at the table of baseline characteristics and asks one question.
“The patients who got chemotherapy were twelve years younger. How do you know it was the chemotherapy?”
Start here
What exactly is a hazard?
Survival tells you how many patients are still event-free by a given time. The hazard tells you how fast events are happening now, among the patients who are still event-free.
Think of a ward of 200 patients who are alive at the start of month 6. During that month 10 of them die. The death rate in that month is about 10 ÷ 200 = 0.05 per patient per month. Make the interval shorter and shorter, and this rate becomes the hazard at that moment. The hazard relates to the event rate now. Survival reflects everything that has happened up to now (Clark et al., 2003).
h(t) ≈ events in a short interval after t ÷ (patients at risk at t × length of the interval)
Cumulative hazard H(t) = −ln S(t)
- h(t)
- the hazard at time t: events per patient per unit of time among patients still at risk
- S(t)
- the probability of being event-free at time t, as on a Kaplan-Meier curve
- H(t)
- the hazard added up from time 0 to t
The difference between risk, rate and hazard, and the Kaplan-Meier curve itself, are covered in the risk and survival chapter.
The Cox model
How can one model compare patients who differ in age, stage and treatment at the same time?
The Cox model says that each factor multiplies the hazard by a fixed amount. It never needs to know the shape of the hazard over time.
h(t) = h₀(t) × exp(b₁x₁ + b₂x₂ + … + bₚxₚ)
- h₀(t)
- the baseline hazard: the hazard for a patient whose x values are all 0. Its shape over time is left unspecified.
- x₁ … xₚ
- the explanatory variables (covariates): chemotherapy 1 or 0, age, stage III 1 or 0
- b
- the regression coefficient: the change in the log of the hazard for a one-unit change in x
- exp(b)
- the hazard ratio for a one-unit change in x, with the other variables held constant
How are the b values found without knowing h₀(t)? At each moment an event happens, the model looks only at the patients still at risk and asks: given that one event happened now, how likely is it that it happened to this patient rather than to another patient at risk?
This is why the Cox model is called semi-parametric. The effect of the covariates is a formula with parameters (the b values), while the baseline hazard is estimated without assuming any distribution. Bradburn et al. (2003) put it this way: in the Cox model the survival times are not assumed to follow a particular statistical distribution. Fully parametric models, such as the exponential or Weibull model, assume a shape for the hazard as well.
Reading an adjusted hazard ratio
- Binary variable (chemotherapy yes = 1, no = 0): HR 0.66 means the hazard with chemotherapy is 0.66 times the hazard without it, at any time, for patients of the same age and stage.
- Continuous variable (age per 10 years): HR 1.54 is the hazard ratio for a 10-year difference. For a 20-year difference it is 1.54² = exp(2 × b) = 2.36.
- Adjusted means holding the other variables in the model constant. It does not remove confounding by variables that were not measured.
Interactive
Change the age gap between the groups and watch the unadjusted and adjusted hazard ratios
Older patients relapse and die sooner. If the chemotherapy group is younger, part of its better survival is the effect of age. The Cox model separates the two.
Try this: move the age gap from 12 years to 0. Watch the unadjusted hazard ratio for chemotherapy come close to the adjusted one: when the two groups are the same age, age can no longer distort the comparison.
Teaching example, simulated cohort of 400 patients after colorectal cancer surgery. The true hazard ratios are set by you (chemotherapy) and fixed at 1.5 per 10 years of age and 2.0 for stage III.
Oncologists often avoid chemotherapy in older patients
1.00 means chemotherapy has no effect
172 of 400 patients relapsed or died. Mean age 56.5 years with chemotherapy and 68.6 years without. Unadjusted, the hazard looks 58% lower with chemotherapy. Adjusted for age and stage, it is 34% lower (95% CI 0.46 to 0.96). The true value you set is 0.70.
| Model | Variable | b | SE | HR (95% CI) | p |
|---|---|---|---|---|---|
| Unadjusted | Chemotherapy | -0.871 | 0.161 | 0.42 (0.31 to 0.57) | < 0.001 |
| Adjusted | Chemotherapy | -0.410 | 0.187 | 0.66 (0.46 to 0.96) | 0.029 |
| Adjusted | Age, per 10 years | 0.430 | 0.086 | 1.54 (1.30 to 1.82) | < 0.001 |
| Adjusted | Stage III | 0.628 | 0.154 | 1.87 (1.39 to 2.54) | < 0.001 |
The formula
Worked example: reading a Cox table line by line
Made-up teaching data: the first cohort in the interactive element above. 400 patients, 172 relapsed or died, median follow-up 40.4 months. Mean age 56.5 years with chemotherapy and 68.6 years without.
| Model | Variable | b | SE | z | HR (95% CI) | p |
|---|---|---|---|---|---|---|
| Unadjusted | Chemotherapy | −0.871 | 0.161 | −5.41 | 0.42 (0.31 to 0.57) | < 0.001 |
| Adjusted | Chemotherapy | −0.410 | 0.187 | −2.19 | 0.66 (0.46 to 0.96) | 0.029 |
| Adjusted | Age, per 10 years | 0.430 | 0.086 | 5.01 | 1.54 (1.30 to 1.82) | < 0.001 |
| Adjusted | Stage III vs II | 0.628 | 0.154 | 4.08 | 1.87 (1.39 to 2.54) | < 0.001 |
- Hazard ratio = exp(b) = exp(−0.410) = 0.664.
- 95% CI = exp(b ± 1.96 × SE) = exp(−0.410 − 0.367) to exp(−0.410 + 0.367) = 0.460 to 0.958.
- Wald z = b ÷ SE = −0.410 ÷ 0.187 = −2.19, giving p = 0.029.
- Unadjusted, chemotherapy appeared to lower the hazard by 58%. After adjustment for age and stage, by 34%. The rest of the unadjusted effect was the age difference.
- The simulation set the true values at 0.70 for chemotherapy, 1.5 per 10 years of age and 2.0 for stage III. Each true value lies inside its 95% CI.
Fitted in Python with lifelines 0.30.3 and checked against the interactive element's own Newton-Raphson fit. There were no tied event times, so the Breslow and Efron methods give the same answer. The concordance (C) statistic of the adjusted model was 0.68. Tests on scaled Schoenfeld residuals gave p = 0.16 for chemotherapy, 0.87 for age and 0.78 for stage. Events per variable: 172 ÷ 3 = 57.
The key assumption
What does "proportional hazards" assume, and how do you check it?
The model gives one hazard ratio for the whole follow-up. That is only a fair summary if the hazard ratio stays roughly the same over time.
Proportional hazards means that the hazard in one group is a fixed multiple of the hazard in the other at every time: the hazard curves are proportional and cannot cross (Bradburn et al., 2003). IBM's SPSS documentation states the same assumption: the hazard ratio should be constant across time.
Three checks
- Look at the Kaplan-Meier curves. Curves that cross, or that separate only after a long time, suggest that the hazard ratio changes.
- Log-minus-log plot. Plot log(−log S(t)) against log(time) for each group. If proportional hazards hold, the lines are parallel (Bradburn et al., 2003).
- Scaled Schoenfeld residuals. For each patient with an event, the residual compares that patient's covariate value with the average among patients still at risk. Plotted against time, the residuals of a variable with a constant effect scatter around a flat line. A trend means the effect changes over time. The test based on these residuals (Grambsch and Therneau, 1994) is cox.zph() in R.
Interactive
Make the hazard ratio change after 12 months and watch the curves, the log-minus-log lines and the single hazard ratio
When the hazard ratio changes over time, a Cox model still reports one number. That number is an average over follow-up, weighted by when the events happen, and it can hide an early harm and a late benefit.
Try this: press "Early harm, late benefit". Watch the survival curves cross at 36 months, the log-minus-log lines cross, and the single hazard ratio move to 1.09, as if the treatment did almost nothing.
Teaching example, made-up hazards: control patients die at 3% per month, follow-up ends at 48 months, equal groups. The single hazard ratio is the value a Cox model would report in a very large trial with these hazards.
Below 1: fewer deaths on the new treatment
Set it equal to the first to make the hazards proportional
The hazards are proportional, so this one number describes the whole follow-up. The log-minus-log lines are parallel.
| Teaching trial | Control | New treatment |
|---|---|---|
| Alive at 12 months | 69.8% | 80.6% |
| Alive at 48 months | 23.7% | 42.1% |
| Mean months alive in the first 48 (restricted mean) | 25.4 | 32.1 |
difference in restricted mean survival to 48 months: +6.7 months
When the assumption fails
What do you do when the hazards are not proportional?
Choose the fix according to which variable breaks the assumption.
| Situation | What to do | What you lose or gain |
|---|---|---|
| A variable you only need to adjust for (centre, tumour grade) breaks the assumption | Stratified Cox model: a separate baseline hazard for each category of that variable | The variable must be categorical, and no hazard ratio is estimated for it, so use this only for variables that are not of primary interest (Bradburn et al., 2003) |
| The exposure you are studying has an effect that changes over time | Time-varying effect: split follow-up (for example before and after 12 months) and estimate a hazard ratio for each period, or include an interaction of the variable with time | You report two or more hazard ratios, which describe the effect honestly |
| The survival curves cross, and a single ratio is not meaningful | Restricted mean survival time (RMST): the area under each survival curve up to a chosen time, compared between groups | An absolute difference in months, which needs no proportional hazards assumption (Royston and Parmar, 2013) |
In the interactive element, "early harm, late benefit" gives a single hazard ratio of 1.09. Split at 12 months, the hazard ratios are 2.0 and 0.5. The restricted mean survival over 48 months is 25.4 months with control and 22.1 months with the new treatment.
Time-varying covariates
What if a patient's exposure starts after follow-up begins?
A patient who starts a drug at month 5 has, by definition, survived the first 5 months. If that patient is called treated from month 0, those survived months are credited to the drug.
Immortal time is a span of cohort follow-up during which, because of the exposure definition, the outcome could not occur (Suissa, 2008). Suissa traces the bias to cohort studies of the survival benefit of heart transplantation in the 1970s, and to later observational studies that reported large benefits of medications.
How to set up the data
Give each patient one row for each period with a constant exposure: start time, stop time, the exposure in that period, and whether the event happened at the stop time. Patient 1 in the figure has two rows: months 0 to 5 untreated with no event, and months 5 to 30 treated with the event. The R survival package fits time-dependent variables with this counting process formulation of Andersen and Gill. In SPSS, time-dependent covariates are defined with the time variable T_.
Laboratory values measured repeatedly (HbA1c, eGFR, CD4 count) are also time-varying covariates. Use only values known at each moment. A value measured later must not be used to predict an earlier event.
Competing risks
Why does 1 − Kaplan-Meier overestimate the risk when patients can die of something else?
A competing risk is an event that prevents the event of interest from happening. A patient who dies cannot later relapse or start dialysis.
Kaplan-Meier treats a patient who dies as censored, which assumes the patient could still have had the event later. So 1 − Kaplan-Meier estimates the risk in an imaginary world where nobody could die first. Austin, Lee and Fine (2016) state that in the presence of competing risks this approach gives estimates of incidence that are biased upward, and that the cumulative incidence function should be used instead.
Cumulative incidence of cause k by time t = probability that the first event happens by t and is of type k
Two kinds of hazard ratio
- Cause-specific hazard: the rate of the event among patients who are currently free of all events. Fit an ordinary Cox model and censor patients at the competing event. Austin et al. (2016) suggest this family of models for questions about causes (aetiology).
- Subdistribution hazard (Fine and Gray): its risk set also keeps patients who have already had the competing event. It links covariates directly to the cumulative incidence, and Austin et al. (2016) suggest it for estimating a patient's prognosis.
Interactive
Raise the death rate and compare 1 − Kaplan-Meier with the cumulative incidence
1 − Kaplan-Meier depends only on the relapse rate. The true proportion of patients who relapse also depends on how many die first.
Try this: raise the death rate from 15 to 40 per 100 patient-years. The 1 − Kaplan-Meier bar for relapse does not move (39.3%), while the cumulative incidence of relapse falls to 18.4%, because more patients die before they can relapse.
Teaching example, made-up constant rates: patients after cancer surgery can relapse, or can die without relapse (the competing event). Both rates are per 100 patient-years.
The event you want to count
The competing event: higher in older, frailer patients
The gap grows with longer follow-up
At 5 years, 1 − Kaplan-Meier overstates the risk of relapse by 10.8 percentage points. It answers a different question: the risk of relapse if no patient could die first. The two 1 − Kaplan-Meier estimates add up to 92.1%; the cumulative incidences add up to 71.3%.
Sample size
How many patients does a Cox analysis need?
Precision and power in survival analysis come from the number of events, not the number of patients (Bradburn et al., 2003).
Events needed (two equal groups) = 4 × (Zα/2 + Zβ)² ÷ (ln HR)²
4 × (1.960 + 0.842)² = 4 × 7.849 = 31.40
31.40 ÷ 0.1272 = 246.8 → 247 events
- Detect HR 0.7 with 80% power at a two-sided α of 5%: 247 events. With 90% power: 331 events.
- If 30% of patients are expected to have the event during follow-up: 247 ÷ 0.30 = 823.3 → 824 patients.
- Allow 10% loss to follow-up: 824 ÷ 0.9 = 915.6 → 916 patients.
- Larger effects need far fewer events: HR 0.5 needs 66 events at 80% power.
Events per variable
In simulations of Cox models, Peduzzi et al. (1995) found biased coefficients and unreliable confidence intervals when there were few events per variable, and judged 10 events per variable the most prudent boundary. Vittinghoff and McCulloch (2007) found circumstances in which fewer than 10 was acceptable, and concluded the rule can be relaxed. With 60 events, a model with 10 candidate predictors is already at that boundary.
Schoenfeld (1983) derived a sample size formula for the proportional hazards model when other covariates are included. Sample size in general is covered in the sample size chapter.
Decide
Which survival model fits your question?
Answer the questions in order.
How many factors do you want to study at once?
Which survival model: all possible answers
- Kaplan-Meier curves and the log-rank test: Describe survival in each group and compare the curves. The hazard ratio from a Cox model with only this factor gives the size of the difference.
- Cause-specific Cox model: Fit an ordinary Cox model for your event and censor patients at the competing event. Report cumulative incidence curves as well.
- Cumulative incidence and the Fine-Gray model: Estimate cumulative incidence (Aalen-Johansen), not 1 − Kaplan-Meier. The Fine-Gray model gives subdistribution hazard ratios.
- Cox model with a time-varying covariate: Split each patient’s follow-up at the moment the exposure changes, so time before the exposure counts as unexposed.
- Standard Cox proportional hazards model: Report adjusted hazard ratios with 95% confidence intervals and say how you checked proportional hazards.
- Stratified Cox model: Each stratum gets its own baseline hazard. No hazard ratio is estimated for the stratifying variable.
- Time-varying effect, or a different summary: Estimate separate hazard ratios for periods of follow-up, or report the difference in restricted mean survival time.
Comparison
Pairs that residents confuse
Log-rank test or Cox regression?
| Feature | Log-rank test | Cox regression |
|---|---|---|
| Answers | Are the survival curves different? | By how much does each factor change the hazard? |
| Result | A p-value only | Hazard ratio with 95% CI for each variable |
| Adjustment | None (a stratified log-rank test allows one stratifying factor) | Several variables, categorical and continuous |
| Continuous predictors | Must be grouped | Used as they are |
| Assumes | Works best when hazards are proportional | Proportional hazards |
Cause-specific or subdistribution hazard?
| Feature | Cause-specific hazard | Subdistribution hazard (Fine-Gray) |
|---|---|---|
| Risk set | Patients free of every event | Patients free of the event of interest, including those who had the competing event |
| How to fit | Ordinary Cox model, censoring at the competing event | Fine-Gray model (crr, or finegray with a weighted Cox model in R) |
| Suited to | Questions about causes (aetiology) | Predicting the proportion who will have the event |
| Link to cumulative incidence | Indirect: depends on the hazards of all causes | Direct |
Stratify or adjust?
| Feature | Adjust for the variable | Stratify on the variable |
|---|---|---|
| Baseline hazard | One, shared | One for each stratum |
| Hazard ratio for that variable | Estimated | Not estimated |
| Assumes proportional hazards for that variable | Yes | No |
| Variable type | Any | Categorical, or grouped into categories |
In a real paper
Cox regression in published studies
Abraham JE, Hiller L, Dorling L, Vallier AL, Dunn J, Bowden S, et al. A nested cohort study of 6,248 early breast cancer patients treated in neoadjuvant and adjuvant chemotherapy trials investigating the prognostic value of chemotherapy-related toxicities. BMC Med. 2015;13:306.
doi:10.1186/s12916-015-0547-5 · PMID 26715442 · Free full text (PMC4693418)
The study: A cohort of 6,248 women with early breast cancer, nested in four UK randomised trials. The authors asked whether patients who had severe side effects of chemotherapy, such as neutropenia, had different relapse-free survival. The patients also differed in nodal status, tumour size, grade and receptor status, so the association had to be adjusted in a Cox model.
“Multivariable analysis identified severe neutropenia (grades ≥3) as an independent predictor of relapse-free survival (hazard ratio (HR) = 0.86; 95 % confidence interval (CI), 0.76–0.97; P = 0.02).”
Multivariable analysis: a Cox model with neutropenia and the other prognostic factors together. Independent predictor: the association remained after adjustment for those factors. HR 0.86: at any time, women with grade 3 or 4 neutropenia had a hazard of relapse or death 0.86 times that of women without it, a 14% lower hazard. 95% CI 0.76 to 0.97: the interval excludes 1, which matches P = 0.02. It is an association, not proof that neutropenia itself protects.
| Reported | Value | What it means |
|---|---|---|
| Relapse-free survival | HR 0.86 (0.76–0.97) | Adjusted hazard ratio for severe neutropenia |
| Breast cancer-specific survival | HR 0.87 (0.75–1.01) | Similar size, but the CI includes 1 (P = 0.06) |
| Events for relapse-free survival | 1,888 (30 %) | Plenty of events for the handful of variables in the model |
| Censored observations | 4,360 (70 %) | Women alive and relapse-free at their last follow-up |
| Median follow-up | 6.2 years | Long enough to see relapses after treatment |
- How proportional hazards were checked: “The proportional hazards assumption was checked using the Schoenfeld residuals method”.
- What was done when it failed: “Subsequent models were stratified by variables that violated the proportional hazards assumption.” The base model was stratified by tumour size, tumour grade and estrogen receptor status, so no hazard ratios are reported for those three variables.
- How variables were chosen: factors significant at P < 0.05 in univariable analysis entered the multivariable model. Selecting variables by p-value can leave out confounders that were not significant on their own.
Al-Wahsh H, Tangri N, Quinn R, Liu P, Ferguson T, Fiocco M, et al. Accounting for the competing risk of death to predict kidney failure in adults with stage 4 chronic kidney disease. JAMA Netw Open. 2021;4(5):e219225.
doi:10.1001/jamanetworkopen.2021.9225 · PMID 33944922 · Free full text (PMC8097501)
The study: A population-based prognostic study of adults with stage 4 chronic kidney disease in two Canadian provinces: 14,619 in the development cohort and 2,295 in the validation cohort. Many of these patients die before they reach kidney failure, so the authors compared a standard Cox model with a cause-specific Cox model and a Fine-Gray model.
“The sum of the 2 Kaplan-Meier failure functions corresponding to the 2 competing events (kidney failure and death) was more than 1 before year 8 in the development cohort and year 6 in the validation cohort (which is impossible because probabilities or risks range between 0 and 1).”
Kaplan-Meier failure function: 1 − Kaplan-Meier, calculated once for kidney failure (censoring deaths) and once for death (censoring kidney failure). More than 1: added together, the two estimates said that more than 100% of patients had one of the events, which cannot be true. This happens because each estimate pretends the other event cannot occur. The authors used the Aalen-Johansen method, the cumulative incidence, for the observed risks.
| Reported | Value | What it means |
|---|---|---|
| Kidney failure rate | 6.1 vs 10.3 per 100 person-years | Development vs validation cohort |
| Mortality rate | 12.3 vs 12.4 per 100 person-years | Death was about twice as common as kidney failure in the development cohort |
| Models compared | Standard Cox, cause-specific Cox, Fine-Gray | The competing-risks models were compared with the model that ignores death |
| Beyond 2 years | Standard Cox over-predicted | Its predicted risks exceeded the observed risks |
- Whether crude risks were estimated with 1 − Kaplan-Meier or with cumulative incidence (here, both were shown).
- Whether the question is about causes (cause-specific hazard) or about the proportion of patients who will have the event (cumulative incidence and Fine-Gray).
- Whether the competing event is common in the patients studied: most patients with stage 4 CKD are 70 years or older and are 3 to 10 times more likely to die than to develop kidney failure.
How to run it
Running a Cox model in SPSS, R and jamovi
Menu paths and functions checked against the IBM SPSS Statistics 24 manuals, the CRAN reference manuals and the jSurvival module documentation.
SPSS (Advanced Statistics option)
- Analyze > Survival > Cox Regression...
- Move the time variable to Time. Move the status variable to Status and click Define Event to give the code for the event.
- Move the explanatory variables to Covariates. Click Categorical to declare categorical variables. A variable placed in Strata gives a stratified analysis (IBM's manual: separate models for different groups).
- Plots: choose Survival, Hazard, Log minus log or One minus survival. Log minus log plots help check proportional hazards.
- Save: Partial residuals. IBM's manual: you can plot partial residuals against survival time to test the proportional hazards assumption.
- Options: tick CI for exp(B) to get 95% confidence intervals for the hazard ratios.
- Time-dependent covariates: Analyze > Survival > Cox w/ Time-Dep Cov..., enter an expression using the time variable T_ (for example T_ multiplied by a covariate), then click Model. The new variable T_COV_ must be included as a covariate.
R (survival and cmprsk packages)
library(survival)
fit <- coxph(Surv(months, event) ~ chemo + age10 + stage3, data = d)
summary(fit) # coef, exp(coef) = hazard ratio, 95% CI
cox.zph(fit) # test of proportional hazards (transform = "km" by default)
plot(cox.zph(fit)) # scaled Schoenfeld residuals against time
concordance(fit) # concordance (C) statistic
# stratified model: separate baseline hazard for each stage
coxph(Surv(months, event) ~ chemo + age10 + strata(stage3), data = d)
# competing risks (cmprsk): status 0 = censored, 1 = relapse, 2 = death
library(cmprsk)
cuminc(ftime = d$months, fstatus = d$status, group = d$chemo, cencode = 0)
crr(ftime = d$months, fstatus = d$status, cov1 = cbind(d$chemo, d$age10),
failcode = 1, cencode = 0)- coxph() uses the Efron method for tied event times by default. Time-dependent variables use the counting process data layout, Surv(start, stop, event).
- survSplit() splits each patient's record at the cut times you give, for estimating separate effects before and after a time point.
- cuminc() estimates cumulative incidence functions and tests equality across groups. crr() fits the Fine-Gray regression. survival::finegray() creates a data set that is then fitted with a weighted Cox model.
jamovi
Cox regression is available through the jSurvival module of the ClinicoPath suite: install it from the jamovi library, then use Survival → ClinicoPath Survival → Multivariable Survival Analysis for adjusted hazard ratios with 95% CIs.
Name the software, the version and the method for tied event times in your methods section.
How to report it
What to write in your thesis or paper
- Define time zero, the event, and when patients were censored.
- Give the number of patients, the number of events and the median follow-up.
- Say which variables were in the model and why (clinical reasoning, not only univariable p-values).
- Report unadjusted and adjusted hazard ratios with 95% CIs, and state the unit of each continuous variable (per 10 years).
- Say how the proportional hazards assumption was checked and what you did if it failed.
- If a competing event is common, say how it was handled and report cumulative incidence.
- Name the software, version and ties method.
Worked sentence, with the teaching numbers:
“Of 400 patients, 172 relapsed or died during a median follow-up of 40.4 months. In a Cox proportional hazards model adjusted for age and stage, adjuvant chemotherapy was associated with a lower hazard of relapse or death (hazard ratio 0.66, 95% CI 0.46 to 0.96; p = 0.029). The unadjusted hazard ratio was 0.42 (95% CI 0.31 to 0.57). Proportional hazards were checked with log-minus-log plots and scaled Schoenfeld residuals, with no evidence of a time-varying effect (p = 0.16 for chemotherapy).”
By specialty
Cox regression problems by specialty
A retrospective cohort of patients with colorectal cancer. Chemotherapy was given more often to younger, fitter patients, and survival is better with chemotherapy on Kaplan-Meier curves.
The unadjusted hazard ratio mixes the effect of chemotherapy with the effect of age. Adjust for age, stage and performance status in a Cox model, and still expect some confounding you cannot measure.
Operation versus medical treatment. Early deaths after surgery are more common, but patients who survive the operation do better later.
The hazard ratio changes over time and the survival curves may cross. Report hazard ratios for early and late periods, or the restricted mean survival time, not one averaged hazard ratio.
A registry study of a drug started at a clinic visit some weeks after myocardial infarction. Patients who died before that visit could never be "users".
Counting users as exposed from the day of the infarction creates immortal time bias. Treat drug use as a time-varying covariate.
Time to kidney failure in older adults with stage 4 chronic kidney disease. Many die before they reach dialysis.
Death is a competing event. 1 − Kaplan-Meier overestimates the risk of kidney failure. Use the cumulative incidence function.
Time to full enteral feeds in very preterm babies. Some babies die before reaching full feeds.
Censoring the babies who die treats them as if they could still reach full feeds. Report the cumulative incidence of full feeds with death as a competing event.
Time to loss to follow-up in a tuberculosis treatment cohort, with deaths during treatment.
A patient who dies cannot be lost to follow-up later. Count deaths as a competing event, not as censored.
Time to extubation in the ICU, with sedation depth recorded every day.
A daily sedation score is a time-varying covariate. Using only the admission value misclassifies patients whose sedation changed.
Time to revision surgery after hip replacement in patients aged over 80 years.
Many patients die with a working implant. 1 − Kaplan-Meier overstates the risk of revision. Report the cumulative incidence of revision with death as a competing event.
Common myths
Common mistaken beliefs
Myth: A semi-parametric model makes no assumptions.
The Cox model leaves the baseline hazard unspecified, but it assumes proportional hazards, independent observations and censoring unrelated to prognosis. For continuous variables it assumes a straight-line effect on the log hazard unless you model it otherwise.
Myth: An adjusted hazard ratio shows that the treatment caused the difference.
Adjustment holds the measured variables constant. In an observational study, unmeasured differences between the groups can still explain part of the result.
Myth: A non-significant Schoenfeld test proves the hazards are proportional.
With few events the test can miss a change in the hazard ratio. Look at the residual plot and the log-minus-log plot as well.
Myth: 1 − Kaplan-Meier gives the risk of relapse even when patients die of other causes.
It treats deaths as censored and overestimates the risk. Use the cumulative incidence function (Austin et al., 2016).
Myth: A patient who ever received the drug can be counted as treated from the start of follow-up.
That credits the drug with the time the patient had to survive before starting it: immortal time bias. In the simulation, a drug with no effect appeared to lower the hazard by 38%.
Myth: More patients always means more power in a survival study.
Power depends on the number of events. 1,000 patients with 20 events tell you less than 300 patients with 150 events.
Why this matters to you
Five consequences of getting Cox regression wrong
Protocol sent back
A sample size based on patients instead of expected events cannot show that the study will observe enough events, and an analysis plan with no check of proportional hazards cannot meet the SAMPL item "Confirm that the assumptions of survival analysis were met".
A wrong conclusion in your thesis
An unadjusted log-rank comparison credits chemotherapy with the effect of youth. In the teaching cohort the hazard ratio moved from 0.42 to 0.66 after adjustment.
Questions you cannot answer at the viva
"Did you check proportional hazards?" and "What did you do with patients who died before relapse?" need specific answers.
Misreading papers
Observational studies of drugs that define users by later prescriptions can report large benefits that are entirely immortal time.
Harm to patients
Ignoring the competing risk of death overestimates the risk of kidney failure. Al-Wahsh et al. (2021) suggest competing-risks predictions to reduce unnecessary interventions or overtreatment.
History
History of Cox regression and its extensions
- 1972
D. R. Cox
Published "Regression models and life-tables": the hazard is a function of explanatory variables multiplied by an unspecified function of time.
- 1974
Nathan Mantel and David Byar
Published "Evaluation of response-time data involving transient states: an illustration using heart-transplant data". Suissa (2008) dates the recognition of immortal time bias to such heart transplantation studies of the 1970s.
- 1975
D. R. Cox
Published "Partial likelihood" in Biometrika.
- 1982
David Schoenfeld
Published "Partial residuals for the proportional hazards regression model", the residuals now used to check proportional hazards.
- 1982
Per Kragh Andersen and Richard Gill
Published "Cox's regression model for counting processes: a large sample study". The R survival package handles time-dependent variables with this counting process formulation.
- 1983
David Schoenfeld
Derived a sample size formula for the proportional hazards regression model when other covariates are included.
- 1994
Patricia Grambsch and Terry Therneau
Published "Proportional hazards tests and diagnostics based on weighted residuals", the basis of the scaled Schoenfeld residuals test.
- 1995
Peduzzi, Concato, Feinstein and Holford
Showed in simulations that Cox models with few events per variable give biased coefficients, and judged 10 events per variable the most prudent boundary.
- 1999
Jason Fine and Robert Gray
Published "A proportional hazards model for the subdistribution of a competing risk", introducing the subdistribution hazard.
- 2008
Samy Suissa
Described immortal time bias in pharmaco-epidemiology and showed that the bias grows with the duration of immortal time.
Guidelines
Survival analysis in reporting guidelines
“Specify the circumstances under which data were censored.” “Confirm that the assumptions of survival analysis were met.” “Report the regression model used to assess the associations between the explanatory variables and survival or time-to-event.” “Report a measure of risk (e.g., a hazard ratio) for each explanatory variable, with a confidence interval.”
Give “unadjusted estimates and, if applicable, confounder-adjusted estimates and their precision (e.g., 95% confidence interval)”.
Cohort study: “if applicable, explain how loss to follow-up was addressed”.
Clarity about the treatment effect “requires a thoughtful envisioning of ‘intercurrent events’”, and the addendum lists “terminal events such as death” among them.
Reading a paper
Five checks on any paper that reports a Cox model
- Are the number of events and the number of variables both given, with roughly 10 or more events per variable?
- Are both unadjusted and adjusted hazard ratios reported with 95% CIs, and is the unit of each continuous variable stated?
- Did the authors say how they checked proportional hazards, and do the Kaplan-Meier curves cross?
- Could any exposure have started after time zero, creating immortal time?
- Is there a common competing event, and was 1 − Kaplan-Meier used where cumulative incidence was needed?
Why is it called that?
Word surgery
“Hazard”
A danger, such as a fire hazard.
The event rate at a given moment among patients who have not yet had the event.
“Proportional”
In fair proportion, a fair share.
The hazard in one group is a fixed multiple of the hazard in the other at every time.
“Competing”
Taking part in a contest.
Describes an event, such as death, that prevents the event of interest from ever happening.
Hazard ratio = exp(b). CI = exp(b ± 1.96 SE). Check proportional hazards. Count time before the exposure as unexposed. Count deaths as a competing event, not as censored.
Hazard
Old French hasard: a game of chance played with dice
In survival analysis, the event rate at a moment among patients still at risk.
Proportional
Latin pro portione: "according to the relation" of parts to each other
Proportional hazards: the relation between the two hazards stays the same over time.
Parameter
Greek para- "beside" + metron "measure"
Semi-parametric: the covariate effects are parameters, while the baseline hazard has no assumed distribution.
Stratum
Latin stratum: "thing spread out"; in English, "horizontal layer" from the 1590s
A stratified Cox model fits a separate baseline hazard in each layer of the stratifying variable.
Compete
Late Latin competere: "strive in common", from com "with, together" + petere "to strive, seek"
Competing risks: two or more events striving for the same patient, only the first of which can be observed.
Censor
Latin censere: "to appraise, value, judge"
A censored patient's follow-up ended before the event was seen.
Exam answers
Short-note outlines
Cox proportional hazards regression
- Purpose: to study how several factors together change the time to an event, with censored follow-up, and to give adjusted hazard ratios.
- Model: h(t) = h₀(t) × exp(b₁x₁ + … + bₚxₚ). h₀(t) is the baseline hazard, left unspecified. exp(b) is the hazard ratio for a one-unit change in x.
- Why semi-parametric: the covariate effects are parameters, but no distribution is assumed for survival times. Estimation uses the partial likelihood, in which h₀(t) cancels.
- Assumptions: proportional hazards, independent observations, non-informative censoring, log-linear effect of continuous variables.
- Interpretation: HR above 1, higher hazard. Below 1, lower hazard. Teaching example: chemotherapy HR 0.66 (95% CI 0.46 to 0.96) adjusted for age and stage, against 0.42 unadjusted.
- Checks: log-minus-log plots, scaled Schoenfeld residuals (cox.zph), the concordance (C) statistic for discrimination.
- Sample size: driven by the number of events. About 10 events per variable (Peduzzi et al., 1995).
- Extensions: stratified Cox, time-varying covariates, time-varying effects, cause-specific and Fine-Gray models for competing risks.
- Reporting: events, follow-up, unadjusted and adjusted HRs with 95% CIs, assumption checks, software.
Proportional hazards assumption: how to check it, and what to do if it fails
- Meaning: the hazard ratio between groups is constant over follow-up. The hazard curves are proportional and do not cross.
- Graphical checks: crossing Kaplan-Meier curves suggest a problem. Log(−log survival) against log time: parallel lines under proportional hazards.
- Residual check: scaled Schoenfeld residuals against time, with a smoothed line: flat if the effect is constant. Grambsch-Therneau test (1994).
- Time interaction: add the covariate multiplied by a function of time (in SPSS, a time-dependent covariate using T_). A significant term means non-proportional hazards.
- If it fails for an adjustment variable: stratified Cox model, which gives no hazard ratio for that variable.
- If it fails for the main exposure: hazard ratios for separate periods of follow-up, a time-varying effect, or restricted mean survival time.
- Example: hazard ratio 2.0 for 12 months and 0.5 afterwards gives a single Cox hazard ratio of 1.09 in a large trial, and the survival curves cross at 36 months.
Competing risks and immortal time bias in survival analysis
Competing risks
- A competing event prevents the event of interest: death before relapse, death before dialysis.
- 1 − Kaplan-Meier censors the competing event and overestimates incidence (Austin et al., 2016). Teaching data: 42.9% against a true 30%.
- Use the cumulative incidence function (Aalen-Johansen estimator) and Gray's test to compare groups.
- Regression: cause-specific Cox model (aetiology) or Fine-Gray subdistribution hazard model (prognosis).
Immortal time bias
- Definition: follow-up during which, because of the exposure definition, the outcome could not occur (Suissa, 2008).
- Arises when patients are classed as exposed from time zero on the basis of something that happened later, such as a later prescription or transplant.
- Effect: makes the exposure look protective. Simulation of a drug with no effect: HR 0.62 wrongly, 1.01 correctly.
- Prevention: treat the exposure as a time-varying covariate, using (start, stop] records for each patient.
The one thing to remember
A Cox model gives one hazard ratio for the whole follow-up, so before you report it, check that the hazard ratio stays the same over time, that no exposure starts after time zero, and that deaths from other causes are not simply censored.
Worked examples and interactive elements use made-up or simulated data for teaching. Numbers were calculated in Python (lifelines 0.30.3, statsmodels 0.14.5, SciPy 1.15.3) and checked against the interactive elements' own code.
Sources: Cox DR. J R Stat Soc Ser B 1972;34(2):187–202 (doi:10.1111/j.2517-6161.1972.tb00899.x). Mantel N, Byar DP. J Am Stat Assoc 1974;69(345):81–6 (doi:10.1080/01621459.1974.10480131). Cox DR. Biometrika 1975;62(2):269–76 (doi:10.1093/biomet/62.2.269). Schoenfeld D. Biometrika 1982;69(1):239–41 (doi:10.1093/biomet/69.1.239). Andersen PK, Gill RD. Ann Stat 1982;10(4) (doi:10.1214/aos/1176345976). Schoenfeld DA. Biometrics 1983;39(2):499–503 (PMID 6354290). Grambsch PM, Therneau TM. Biometrika 1994;81(3):515–26 (doi:10.1093/biomet/81.3.515). Peduzzi P et al. J Clin Epidemiol 1995;48(12):1503–10 (doi:10.1016/0895-4356(95)00048-8). Fine JP, Gray RJ. J Am Stat Assoc 1999;94(446):496–509 (doi:10.1080/01621459.1999.10474144). Clark TG et al. Br J Cancer 2003;89(2):232–8 (doi:10.1038/sj.bjc.6601118). Bradburn MJ et al. Br J Cancer 2003;89(3):431–6 (doi:10.1038/sj.bjc.6601119). Bradburn MJ et al. Br J Cancer 2003;89(4):605–11 (doi:10.1038/sj.bjc.6601120). Vittinghoff E, McCulloch CE. Am J Epidemiol 2007;165(6):710–8 (doi:10.1093/aje/kwk052). von Elm E et al. STROBE statement. Lancet 2007;370:1453–7 (PMID 18064739). Suissa S. Am J Epidemiol 2008;167(4):492–9 (doi:10.1093/aje/kwm324). Royston P, Parmar MKB. BMC Med Res Methodol 2013;13:152 (doi:10.1186/1471-2288-13-152). Lang TA, Altman DG. SAMPL guidelines. Int J Nurs Stud 2015;52(1):5–9 (doi:10.1016/j.ijnurstu.2014.09.006). Abraham JE et al. BMC Med 2015;13:306 (doi:10.1186/s12916-015-0547-5). Austin PC, Lee DS, Fine JP. Circulation 2016;133(6):601–9 (doi:10.1161/CIRCULATIONAHA.115.017719). ICH E9(R1) addendum on estimands and sensitivity analysis, adopted 20 November 2019. Al-Wahsh H et al. JAMA Netw Open 2021;4(5):e219225 (doi:10.1001/jamanetworkopen.2021.9225). IBM SPSS Advanced Statistics 24 manual. R packages survival 3.8-12 and cmprsk 2.2-12, CRAN reference manuals. jSurvival (ClinicoPath) module documentation.
Further reading: Indrayan A. Medical Biostatistics, chapters on survival analysis. Jawlik A. Statistics from A to Z: Confusing Concepts Clarified. Wiley; 2016 (Exponential Distribution).