# 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.
*\[Rich HTML output -- view on the documentation site\]*
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.
``` python
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.
``` python
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.
``` python
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.
``` python
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.
``` python
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()](../reference/AFT.md#greenwood.AFT.to_frame) and take the first couple of rows within each cause, so both curves are visible side by side.
``` python
aj.to_frame(format="polars").group_by("cause").head(2)
```
PolarsRows4Columns7
|
|
cause
str
|
time
f64
|
n_risk
f64
|
estimate
f64
|
std_error
f64
|
conf_low
f64
|
conf_high
f64
|
| 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.
> **Warning: 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.
``` python
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.
``` python
import polars as pl
aj_sex.to_frame(format="polars").filter(pl.col("cause") == "pcm").group_by("strata").tail(1)
```
PolarsRows2Columns8
|
|
strata
str
|
cause
str
|
time
f64
|
n_risk
f64
|
estimate
f64
|
std_error
f64
|
conf_low
f64
|
conf_high
f64
|
| 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.
``` python
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.
``` python
gw.tidy(fg, exponentiate=True, format="polars")
```
PolarsRows2Columns7
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
| 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.
``` python
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])
```
PandasRows4Columns8
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
cause
str
|
| 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.
> **Note: 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](multi-state.md) generalize competing risks to settings where subjects move between several states over time.
- [Visualizing survival](visualization.md) shows the grammar-of-graphics approach you can reuse to plot cumulative incidence curves.
- Revisit [Survival data](survival-data.md) for the details of the multi-state response.