Competing risks

So far every method has assumed a single kind of event. Often that is not realistic. A patient with a precancerous condition might progress to cancer, but they might also die of something else first, and death from another cause prevents the progression from ever being observed. These are competing risks: mutually exclusive event types where the occurrence of one removes the possibility of the others. Analyzing them with ordinary survival methods, one event type at a time, gives biased and even nonsensical results. This page explains why, and shows the two standard tools Greenwood provides.

flowchart LR
    A[Event-free] --> B[Progression]
    A --> C[Death]

Figure 1: A competing-risks model: from an initial state, a subject can move to exactly one of several absorbing event states.

We use the bundled mgus2 dataset, which follows patients with monoclonal gammopathy. The competing events are progression to a plasma-cell malignancy and death without progression. We first build the competing-risks endpoint from two intermediate arrays: a single event time and a cause code that is 0 for censored, 1 for progression, and 2 for death.

import numpy as np
import greenwood as gw

mg = gw.load_dataset("mgus2", backend="polars")

etime = np.where(mg["pstat"] == 1, mg["ptime"], mg["futime"])
cause = np.where(mg["pstat"] == 1, 1, 2 * mg["death"])   # 0 censor, 1 pcm, 2 death

The etime array picks each subject’s event time: the progression time ptime when the subject progressed (pstat == 1), and otherwise the follow-up time futime. Because a subject can experience only the first of the competing events, one time per subject is all we need. Here are the first several values, in months.

etime[:8]

array([ 30,  25,  46,  92,   8,   4, 151,   2])

The cause array records which event that time refers to, using the coding noted in the comment above. Where a subject progressed it is 1 (pcm); otherwise it is 2 * death, which is 2 when the subject died without progressing and 0 when the subject was censored. Reading etime and cause position by position gives each subject’s (time, cause) pair.

cause[:8]

array([2, 2, 2, 2, 2, 2, 2, 2])

We pass both arrays to gw.Surv.multistate, mapping cause code 1 to pcm and code 2 to death, with code 0 left as censoring. The resulting response object shows the event times alongside their decoded states.

y = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))
y

Surv(type=right, n=1384, events=975, states=('pcm', 'death'))

Cumulative incidence, not one-minus-survival

The intuitive but wrong approach is to treat one cause as the event, censor the other, run Kaplan-Meier, and report one minus the survival curve as the probability of that cause. This overestimates the probability, sometimes badly, because it treats subjects who died of the other cause as if they could still progress. The correct quantity is the cumulative incidence function, which gives the probability of experiencing a specific cause by a given time while accounting for the competing cause.

The Aalen-Johansen estimator computes it. You fit it to the multi-state response, and it returns a curve for every cause.

mg = gw.load_dataset("mgus2", backend="polars")

etime = np.where(mg["pstat"] == 1, mg["ptime"], mg["futime"])
cause = np.where(mg["pstat"] == 1, 1, 2 * mg["death"])   # 0 censor, 1 pcm, 2 death
y = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))

aj = gw.AalenJohansen().fit(y)
aj

AalenJohansen (Aalen-Johansen cumulative incidence)

states: pcm, death
n = 1384

       final CIF
pcm       0.1613
death     0.8387

The fitted estimator holds a cumulative incidence curve for each cause. To inspect it we call to_frame() and take the first couple of rows within each cause, so both curves are visible side by side.

aj.to_frame(format="polars").group_by("cause").head(2)

	cause str	time f64	n_risk f64	estimate f64	std_error f64	conf_low f64	conf_high f64
PolarsRows4Columns7
0	pcm	1	1384	0	0	0	0
1	pcm	2	1341	0.00144616432392	0.00102185289691	0	0.00344895919936
2	death	1	1384	0.0303468208092	0.00461101747295	0.0213093926302	0.0393842489883
3	death	2	1341	0.0505931213442	0.00589210798256	0.0390448019053	0.062141440783

Each cause has its own cumulative incidence estimate, standard error, and confidence interval at every event time. Because the causes are competing, their cumulative incidences plus the event-free probability sum to one at every time, which is exactly the accounting that one-minus-Kaplan-Meier gets wrong.

Do not use one-minus-Kaplan-Meier for a competing cause

Censoring the competing event and applying Kaplan-Meier treats those subjects as if they remained at risk for the cause of interest, which they do not. Always use the cumulative incidence function for competing-risks probabilities.

You can stratify by a grouping variable, just as with Kaplan-Meier, to compare cumulative incidence across groups.

aj_sex = gw.AalenJohansen().fit(y, by=mg["sex"])
aj_sex

AalenJohansen (Aalen-Johansen cumulative incidence)

states: pcm, death
strata: 2

The stratified fit carries a separate set of curves for each level of sex. Here we pull the pcm cause and keep the last row within each stratum, which is the cumulative incidence of progression at the final event time in that group.

import polars as pl

aj_sex.to_frame(format="polars").filter(pl.col("cause") == "pcm").group_by("strata").tail(1)

	strata str	cause str	time f64	n_risk f64	estimate f64	std_error f64	conf_low f64	conf_high f64
PolarsRows2Columns8
0	M	pcm	424	1	0.104460229977	0.0158156196773	0.0734621850159	0.135458274937
1	F	pcm	394	1	0.198554321073	0.0399976113701	0.12016044332	0.276948198826

Modeling with the Fine-Gray model

The cumulative incidence function describes groups. To model how covariates affect the incidence of a specific cause, use the Fine-Gray model. It is a regression on the subdistribution hazard, constructed so that its coefficients describe the effect of covariates on the cumulative incidence itself, which is usually what you want to report.

mg = gw.load_dataset("mgus2", backend="polars")

etime = np.where(mg["pstat"] == 1, mg["ptime"], mg["futime"])
cause = np.where(mg["pstat"] == 1, 1, 2 * mg["death"])   # 0 censor, 1 pcm, 2 death
y = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))

fg = gw.FineGray("pcm").fit(y, mg[["age", "sex"]])
fg

FineGray (Fine-Gray subdistribution hazard model, cause='pcm')

         coef  exp(coef)  se(coef)       z         p
age   -0.0173     0.9828  0.005701  -3.035  0.002408
sexM  -0.2597     0.7713    0.1856  -1.399    0.1617

n = 1384, events = 115
Standard errors: robust (clustered)

We name pcm as the target cause and fit against age and sex. To read the fitted model we tidy it into a coefficient table, asking for exponentiated estimates so the numbers are on the hazard-ratio scale rather than the log scale.

gw.tidy(fg, exponentiate=True, format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows2Columns7
0	age	0.982848064921	0.00570095802298	-3.03470633443	0.00240769998927	0.9719271695	0.993891671139
1	sexM	0.771282230528	0.185579308546	-1.39940662934	0.161691080051	0.536102622861	1.1096313537

The exponentiated coefficients are subdistribution hazard ratios. A value above 1 means the covariate increases the cumulative incidence of the target cause, pcm here, over time. The model reports clustered robust standard errors by default, which are appropriate given the inverse-probability weighting it uses internally.

Modeling multiple competing causes

In a competing-risks study with more than one cause of interest, you can fit separate Fine-Gray models, one for each cause, to compare the covariate effects across causes.

import pandas as pd

# Fit FineGray models for both competing causes
fg_pcm = gw.FineGray("pcm").fit(y, mg[["age", "sex"]])
fg_death = gw.FineGray("death").fit(y, mg[["age", "sex"]])

# Compare the hazard ratios side-by-side
pcm_tidy = gw.tidy(fg_pcm, exponentiate=True, format="pandas")
death_tidy = gw.tidy(fg_death, exponentiate=True, format="pandas")

pcm_tidy["cause"] = "PCM"
death_tidy["cause"] = "Death"

pd.concat([pcm_tidy, death_tidy])

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64	cause str
PandasRows4Columns8
0	age	0.982848064921	0.00570095802298	-3.03470633443	0.00240769998927	0.9719271695	0.993891671139	PCM
1	sexM	0.771282230528	0.185579308546	-1.39940662934	0.161691080051	0.536102622861	1.1096313537	PCM
2	age	1.06035605364	0.0036793490261	15.9280217679	4.04967555101e-57	1.05273691623	1.06803033421	Death
3	sexM	1.44900400068	0.0667809170816	5.55362879898	2.79799326688e-08	1.27123419034	1.65163320019	Death

This comparison reveals how covariates affect each cause differently. A covariate might increase the incidence of progression to PCM while decreasing the incidence of death from other causes (or vice versa) which both tell important clinical stories.

Two kinds of hazard, two kinds of question

A cause-specific Cox model answers “what drives the rate of this event among those still at risk?”, while the Fine-Gray model answers “what drives the cumulative probability of this event?”. They can point in different directions for the same covariate, and both can be correct. Decide which question you are asking before choosing the model.

Next steps

You can now estimate and model competing-risks data correctly.

Multi-state models generalize competing risks to settings where subjects move between several states over time.
Visualizing survival shows the grammar-of-graphics approach you can reuse to plot cumulative incidence curves.
Revisit Survival data for the details of the multi-state response.