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

# Load the mgus2 dataset: monoclonal-gammopathy patients with competing-risk follow-up
mg = gw.load_dataset("mgus2", backend="polars")

# Event time: progression time if progressed, otherwise overall follow-up time
etime = np.where(mg["pstat"] == 1, mg["ptime"], mg["futime"])

# Cause code: 1 = pcm progression, 2 = death without progression, 0 = censored
cause = np.where(mg["pstat"] == 1, 1, 2 * mg["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.

# Show the first eight event times (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.

# Show the first eight cause codes (0 censored, 1 pcm, 2 death)
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.

# Build the multi-state response, naming the two competing causes
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.

# Rebuild the dataset and competing-risks response
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"))

# Fit the Aalen-Johansen estimator to get a CIF for each competing cause
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.

# Export to a tidy frame and show the first two rows per cause
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	death	1	1384	0.0303468208092	0.00461101747295	0.0213093926302	0.0393842489883
1	death	2	1341	0.0505931213442	0.00589210798256	0.0390448019053	0.062141440783
2	pcm	1	1384	0	0	0	0
3	pcm	2	1341	0.00144616432392	0.00102185289691	0	0.00344895919936

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.

Stratified cumulative incidence

You can compare cumulative incidence across groups by passing by= to fit(), just as with Kaplan-Meier. Each group receives its own independent set of CIF curves, estimated and stored together in the fitted object.

# Fit separate CIFs for each level of sex
aj_sex = gw.AalenJohansen().fit(y, by=mg["sex"])
aj_sex

AalenJohansen (Aalen-Johansen cumulative incidence)

states: pcm, death
strata: 2

The summary shows the final cumulative incidence per cause for each sex group. To drill into a specific cause across strata, filter the tidy frame and keep the last row per stratum, which is the cumulative incidence at the final observed event time in that group.

import polars as pl

# Final pcm cumulative incidence by sex stratum (last event time in each group)
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	F	pcm	394	1	0.198554321073	0.0399976113701	0.12016044332	0.276948198826
1	M	pcm	424	1	0.104460229977	0.0158156196773	0.0734621850159	0.135458274937

Each row gives the probability of pcm progression accumulated by the last observed event time in that sex group. The Aalen-Johansen estimator handles both the unstratified and stratified case through the same interface: construct the estimator, then call fit() with an optional by= argument to produce all group curves in one step.

Confidence intervals

The confidence band around each CIF curve is computed from the Aalen (Marubini-Valsecchi) delta-method variance. Greenwood supports three transforms for the interval, chosen with conf_type=. The default is "plain", which gives a symmetric Wald interval directly on the probability scale, matching R’s survfit default for competing risks. The "log" transform constructs the interval on the log-CIF scale before back-transforming, and "log-log" uses the same log(-log(F)) scale that R applies to the CIF.

You set the transform at construction time.

# Rebuild the dataset and competing-risks response
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"])
y = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))

# Fit with the log-log confidence interval transform
aj_loglog = gw.AalenJohansen(conf_type="log-log").fit(y)

# Inspect the pcm curve with its log-log bounds
aj_loglog.to_frame(format="polars").filter(pl.col("cause") == "pcm").head(5)[
    ["cause", "time", "estimate", "conf_low", "conf_high"]
]

	cause str	time f64	estimate f64	conf_low f64	conf_high f64
PolarsRows5Columns5
0	pcm	1	0	NaN	NaN
1	pcm	2	0.00144616432392	0.000309260554073	0.00503788210186
2	pcm	3	0.00144616432392	0.000309260554073	0.00503788210186
3	pcm	4	0.00216924648589	0.000626978552936	0.00609026601181
4	pcm	5	0.00289232864785	0.000996829775409	0.00712140883971

The point estimates are identical across transforms; only the interval limits differ. At time = 0 (before any events), conf_low and conf_high are NaN for the "log" and "log-log" transforms, since log(0) is undefined.

Comparing confidence interval transforms

Three transforms are available, all built on the same delta-method standard error. The code below fits all three to the same data and assembles a comparison frame from sampled rows of the pcm curve.

# Fit all three CI transforms
aj_plain  = gw.AalenJohansen(conf_type="plain").fit(y)
aj_log    = gw.AalenJohansen(conf_type="log").fit(y)
aj_loglog = gw.AalenJohansen(conf_type="log-log").fit(y)

# PCM curve: take a handful of rows from each to compare limits
rows = []
for aj, name in [(aj_plain, "plain"), (aj_log, "log"), (aj_loglog, "log-log")]:
    df = aj.to_frame(format="polars").filter(pl.col("cause") == "pcm")
    rows.append(
        df[::40]
        .select(["time", "estimate", "conf_low", "conf_high"])
        .with_columns(pl.lit(name).alias("transform"))
    )

pl.concat(rows)

	time f64	estimate f64	conf_low f64	conf_high f64	transform str
PolarsRows21Columns5
0	1	0	0	0	plain
1	41	0.0246078535599	0.0164387541886	0.0327769529311	plain
2	81	0.0472283168356	0.0359254787521	0.0585311549192	plain
3	121	0.0656390406572	0.0520838109942	0.0791942703203	plain
4	161	0.080546864453	0.0651315150291	0.0959622138769	plain
5	203	0.0937689374365	0.0763226592751	0.111215215598	plain
6	260	0.104027675794	0.083315654599	0.12473969699	plain
7	1	0	NaN	NaN	log
18	161	0.080546864453	0.0660226169531	0.0968470522727	log-log
19	203	0.0937689374365	0.0772671987375	0.112141892046	log-log
20	260	0.104027675794	0.0844817169034	0.125870192136	log-log

Transform summary:

plain (default): symmetric interval on the probability scale. Can stray outside [0, 1] at very early or late times, but is clipped.
log: interval on the log-CIF scale (bounds stay positive).
log-log: interval on the log(-log(F)) scale; matches R’s conf.type = "log-log". Produces NaN at F = 0.

Choose the transform before fitting and report which one you used. The "plain" option is easiest to explain to a non-statistical audience. The "log" and "log-log" options give more reliable bounds near the edges of the probability scale.

Adjusting the confidence level

By default, Greenwood computes 95% confidence intervals. You can change this with conf_level=.

# Fit with a 90% confidence level instead of the default 95%
aj_90 = gw.AalenJohansen(conf_level=0.90).fit(y)

# Inspect the narrower interval on the first three rows of the pcm curve
aj_90.to_frame(format="polars").filter(pl.col("cause") == "pcm").head(3)[
    ["time", "estimate", "conf_low", "conf_high"]
]

	time f64	estimate f64	conf_low f64	conf_high f64
PolarsRows3Columns4
0	1	0	0	0
1	2	0.00144616432392	0	0.00312696276761
2	3	0.00144616432392	0	0.00312696276761

Narrowing the confidence level to 0.90 tightens the bounds on both sides of the estimate. The conf_level= parameter applies uniformly to whichever conf_type is in use.

The conf_type= and conf_level= parameters are independent: mix and match any combination to suit your reporting requirements. All choices produce the same point estimate. Only the interval width and its scale transform differ.

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.

# Rebuild the dataset and response
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"))

# Fit the Fine-Gray model for the pcm cause with age and sex as covariates
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.

# Exponentiate to get subdistribution hazard ratios
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.04967555089e-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.

Plot the results with plot_cif() from Visualizing survival. It accepts a fitted AalenJohansen directly and facets by cause, with optional grouping by covariate and an aligned numbers-at-risk table.
Multi-state models generalize competing risks to settings where subjects move between several states over time.
Revisit Survival data for the details of the multi-state response.