---
title: "Competing risks and cumulative incidence"
description: "When animals die from several mutually exclusive causes, one minus the Kaplan-Meier curve overstates each cause. Estimate the cumulative incidence function."
date: "2026-06-09 09:00"
categories: [R, survival analysis, competing risks, ecology tutorial]
image: thumbnail.png
image-alt: "Two rising curves showing that one minus Kaplan-Meier overstates the cumulative incidence of a cause of death relative to the cumulative incidence function."
date-modified: "2026-10-15"
---
*Corrected 15 October 2026: the caption and alt text of the hazard-ratio figure said one covariate group's predation curve sits a little above the other's throughout; they now say the gap opens in the first year and narrows towards the end of follow-up.*
*Corrected 7 October 2026: the post said the Fine-Gray weight given to an animal that died of the competing cause decays over time; it falls only as censoring accumulates, and with every survivor censored at the end of follow-up the weights here never fall, which the text now says.*
*Corrected 30 September 2026: the post said most of the collared animals were removed by other causes before a predator found them; fewer than half died of other causes, so the text now says more animals died of other causes than of predation.*
*Corrected 27 September 2026: the caption of the hazard-ratio figure said the cumulative incidence of predation was nearly identical in the two habitats; it now says the x = 1 curve sits a little above the x = 0 curve throughout, which is what the figure shows.*
A radio-collared animal rarely just "dies". It dies of predation, or starvation, or disease, and the moment one cause strikes, the others can no longer occur. These are *competing risks*, and the habit carried over from single-event survival analysis, taking one minus a Kaplan-Meier curve to get the probability of a given cause, quietly inflates every number it touches.
This post shows the size of that inflation, the estimator that fixes it (the cumulative incidence function, via the Aalen-Johansen estimator), and the trap waiting on the modelling side: the cause-specific hazard and the Fine-Gray subdistribution hazard answer different questions, and a covariate can push one up while leaving the other flat.
```{r}
#| label: setup
#| message: false
#| warning: false
#| code-fold: true
#| code-summary: "Setup: packages, colours and plot theme"
library(ggplot2)
library(survival)
library(grid)
FOREST <- "#275139"; GOLD <- "#c9b458"; RUST <- "#b5534e"; SAGE <- "#93a87f"
INK <- "#16241d"; BODY <- "#2c3a31"; LINE <- "#dad9ca"
theme_te <- function() {
theme_minimal(base_size = 12) +
theme(plot.title = element_text(colour = INK, face = "bold", size = 12),
plot.subtitle = element_text(colour = BODY, size = 10),
axis.title = element_text(colour = BODY),
axis.text = element_text(colour = BODY),
legend.position = "none",
panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = LINE, linewidth = 0.3),
plot.background = element_rect(fill = "#f5f4ee", colour = NA),
panel.background = element_rect(fill = "#f5f4ee", colour = NA))
}
two_panel <- function(pa, pb) {
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
print(pa, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(pb, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
}
```
## The set-up: two causes, one clock
Suppose we follow 300 collared animals for two and a half years. Two causes remove them: predation (cause 1) and everything else pooled as "other" (cause 2). Each animal has a latent time to predation and a latent time to an other-cause death; whichever comes first (or the end of the study) is what we observe.
```{r}
#| label: sim-a
set.seed(4256)
n <- 300; rate_pred <- 0.30; rate_other <- 0.45; tau <- 2.5
t_pred <- rexp(n, rate_pred)
t_other <- rexp(n, rate_other)
obs_t <- pmin(t_pred, t_other, tau)
cause <- ifelse(obs_t == tau, 0L, ifelse(t_pred < t_other, 1L, 2L))
n_pred <- sum(cause == 1)
n_other <- sum(cause == 2)
n_cens <- sum(cause == 0)
```
Of the 300 animals, `r n_pred` were predated, `r n_other` died of other causes, and `r n_cens` were still alive (right-censored) at the end. Predation is the minority cause here: the other-cause hazard is higher, so more animals are removed by something else than by a predator.
## One minus Kaplan-Meier overstates the cause
The tempting move is to treat predation as "the event", censor the other-cause deaths, fit a Kaplan-Meier curve, and read one minus the survival as the probability of dying by predation. Alongside it, the correct estimator: the cumulative incidence function (CIF), which the Aalen-Johansen estimator returns directly from a multi-state fit.
```{r}
#| label: cif-a
km_pred <- survfit(Surv(obs_t, cause == 1) ~ 1)
naive_cif <- 1 - summary(km_pred, times = tau, extend = TRUE)$surv
ev <- factor(cause, levels = c(0, 1, 2), labels = c("censor", "pred", "other"))
aj <- survfit(Surv(obs_t, ev) ~ 1)
pred_col <- which(aj$states == "pred")
other_col <- which(aj$states == "other")
si <- summary(aj, times = tau, extend = TRUE)
aj_cif_pred <- si$pstate[, pred_col]
aj_cif_other <- si$pstate[, other_col]
theo_cif <- (rate_pred / (rate_pred + rate_other)) *
(1 - exp(-(rate_pred + rate_other) * tau))
overstate_pct <- 100 * (naive_cif - aj_cif_pred) / aj_cif_pred
```
The naive curve puts the probability of death by predation at `r sprintf("%.3f", naive_cif)`. The cumulative incidence function puts it at `r sprintf("%.3f", aj_cif_pred)`. That is an overstatement of `r sprintf("%.0f", overstate_pct)` per cent, and it is not a fluke of this sample; the mechanism is structural.
One minus Kaplan-Meier asks "what fraction would eventually be predated if other-cause deaths were merely a temporary interruption?" But an animal killed by disease is not interrupted; it is gone, and can never be predated. Censoring those animals treats them as if they re-enter the predation risk set later, which inflates the estimate. The CIF keeps the accounting honest: the other-cause deaths stay out of the numerator, and the two causes plus the survivors add up exactly to one.
```{r}
#| label: coherence-a
cif_total <- aj_cif_pred + aj_cif_other
km_all <- survfit(Surv(obs_t, cause > 0) ~ 1) # all-cause mortality, estimated
allcause <- 1 - summary(km_all, times = tau, extend = TRUE)$surv
surv_total <- 1 - exp(-(rate_pred + rate_other) * tau) # population value from the true rates
```
Here the two cumulative incidences sum to `r sprintf("%.3f", cif_total)`, which is exactly one minus the all-cause Kaplan-Meier survival, `r sprintf("%.3f", allcause)`; the estimator enforces that identity in every sample. The true rates put the population figure at `r sprintf("%.3f", surv_total)`, and the small gap is sampling error, not a failure of the accounting. The naive curve has no such guarantee: add up one minus Kaplan-Meier across both causes and the total can exceed one.
```{r}
#| label: fig-cif
#| fig-cap: "Left: one minus Kaplan-Meier for predation overstates the risk relative to the cumulative incidence function. Right: the stacked cumulative incidence splits total mortality into the two causes, always summing to one minus overall survival."
#| fig-alt: "Left panel, the naive one minus Kaplan-Meier curve rises well above the cumulative incidence curve for predation. Right panel, a stacked cumulative incidence with a lower band for predation and an upper band for other causes."
#| code-fold: true
#| code-summary: "Show the code behind this output"
naive_df <- data.frame(
time = c(0, km_pred$time),
cif = c(0, 1 - km_pred$surv),
what = "one minus Kaplan-Meier")
aj_df <- data.frame(
time = c(0, aj$time),
cif = c(0, aj$pstate[, pred_col]),
what = "cumulative incidence")
cif_lines <- rbind(naive_df, aj_df)
p1 <- ggplot(cif_lines, aes(time, cif, colour = what)) +
geom_step(linewidth = 0.9) +
scale_colour_manual(values = c("one minus Kaplan-Meier" = RUST,
"cumulative incidence" = FOREST)) +
annotate("text", x = 1.15, y = 0.55, label = "one minus KM", colour = RUST, hjust = 0, size = 3.4) +
annotate("text", x = 1.35, y = 0.24, label = "cumulative incidence", colour = FOREST, hjust = 0, size = 3.4) +
labs(title = "Predation risk, two estimators", x = "years", y = "probability of death by predation") +
coord_cartesian(ylim = c(0, 0.7)) +
theme_te()
stack_df <- data.frame(
time = c(0, aj$time),
lower = c(0, aj$pstate[, pred_col]),
upper = c(0, aj$pstate[, pred_col] + aj$pstate[, other_col]))
p2 <- ggplot(stack_df, aes(time)) +
geom_step(aes(y = upper), colour = GOLD, linewidth = 0.9) +
geom_step(aes(y = lower), colour = FOREST, linewidth = 0.9) +
annotate("text", x = 1.2, y = 0.18, label = "predation", colour = FOREST, hjust = 0, size = 3.4) +
annotate("text", x = 1.5, y = 0.5, label = "other causes", colour = "#8a7420", hjust = 0, size = 3.4) +
annotate("text", x = 1.7, y = 0.92, label = "still alive", colour = SAGE, hjust = 0, size = 3.4) +
labs(title = "Stacked cumulative incidence", x = "years", y = "cumulative probability") +
coord_cartesian(ylim = c(0, 1)) +
theme_te()
two_panel(p1, p2)
```
## Two hazards, two questions
The modelling side hides a subtler point. Add a covariate, say a habitat indicator that raises the risk of both causes, and ask how it affects predation. There are two hazards you could model, and they do not agree.
The *cause-specific hazard* is the instantaneous predation rate among animals still alive and at risk; you fit it with a Cox model that censors the other-cause deaths. The *subdistribution hazard* of Fine and Gray is built so that its coefficient maps directly onto the cumulative incidence; it keeps animals who died of a competing cause in the risk set, with a weight that falls only as censoring accumulates (here every survivor is censored at the end of follow-up, so the weights never fall).
```{r}
#| label: sim-b
set.seed(42562)
nb <- 1200
x <- rbinom(nb, 1, 0.5)
h_pred <- 0.25 * exp(log(1.8) * x)
h_other <- 0.40 * exp(log(1.8) * x)
taub <- 3
tp <- rexp(nb, h_pred); to <- rexp(nb, h_other)
tb <- pmin(tp, to, taub)
cb <- ifelse(tb == taub, 0L, ifelse(tp < to, 1L, 2L))
db <- data.frame(time = tb, cause = cb, x = x)
```
```{r}
#| label: fit-b
#| warning: false
csh <- coxph(Surv(time, cause == 1) ~ x, data = db)
csh_hr <- exp(coef(csh)); csh_ci <- exp(confint(csh))
evb <- factor(cb, levels = c(0, 1, 2), labels = c("censor", "pred", "other"))
fg <- finegray(Surv(time, evb) ~ ., data = data.frame(time = tb, x = x, evb = evb),
etype = "pred")
fgm <- coxph(Surv(fgstart, fgstop, fgstatus) ~ x, weight = fgwt, data = fg)
fg_hr <- exp(coef(fgm)); fg_ci <- exp(confint(fgm))
cif_group <- function(sub) {
e <- factor(sub$cause, levels = c(0, 1, 2), labels = c("c", "pred", "other"))
a <- survfit(Surv(sub$time, e) ~ 1)
s <- summary(a, times = taub, extend = TRUE)
s$pstate[, which(a$states == "pred")]
}
cif_x0 <- cif_group(db[db$x == 0, ])
cif_x1 <- cif_group(db[db$x == 1, ])
```
The habitat covariate raises the cause-specific predation hazard by a clear margin: the cause-specific hazard ratio is `r sprintf("%.2f", csh_hr)` (95 per cent interval `r sprintf("%.2f", csh_ci[1])` to `r sprintf("%.2f", csh_ci[2])`). Read as a rate among survivors, animals in this habitat are predated faster.
Yet the Fine-Gray subdistribution hazard ratio is only `r sprintf("%.2f", fg_hr)` (interval `r sprintf("%.2f", fg_ci[1])` to `r sprintf("%.2f", fg_ci[2])`, spanning one), and the cumulative incidence of predation moves far less than the rate does: `r sprintf("%.3f", cif_x0)` at `x` equal to zero against `r sprintf("%.3f", cif_x1)` at `x` equal to one, a difference the Fine-Gray interval cannot separate from none. The reason is the competing cause. The same habitat also raises other-cause mortality, so animals there tend to die of something else before a predator reaches them. The elevated predation *rate* is mostly cancelled by faster removal through the competing risk, and the *probability of actually dying by predation* moves only a little.
```{r}
#| label: fig-hazards
#| echo: false
#| fig-cap: "Left: the cause-specific hazard ratio for the covariate sits clearly above one, while the Fine-Gray subdistribution ratio is close to one. Right: the cumulative incidence of predation is higher for x = 1 throughout, but the gap that opens in the first year narrows towards the end of follow-up, because the covariate also raises the competing cause."
#| fig-alt: "Left panel, a point and interval for the cause-specific hazard ratio near 1.7 and a Fine-Gray ratio near 1.15 spanning one. Right panel, two cumulative incidence step curves for the two covariate groups, the x = 1 curve above the x = 0 curve throughout, with a gap that opens in the first year and narrows towards the end."
forest_df <- data.frame(
model = factor(c("cause-specific\nhazard", "Fine-Gray\nsubdistribution"),
levels = c("Fine-Gray\nsubdistribution", "cause-specific\nhazard")),
hr = c(csh_hr, fg_hr),
lo = c(csh_ci[1], fg_ci[1]),
hi = c(csh_ci[2], fg_ci[2]))
p3 <- ggplot(forest_df, aes(hr, model)) +
geom_vline(xintercept = 1, colour = LINE, linewidth = 0.5) +
geom_pointrange(aes(xmin = lo, xmax = hi), colour = FOREST, linewidth = 0.8, size = 0.6) +
labs(title = "Two hazard ratios", x = "hazard ratio for predation", y = NULL) +
coord_cartesian(xlim = c(0.8, 2.3)) +
theme_te()
g0 <- db[db$x == 0, ]; g1 <- db[db$x == 1, ]
grid_cif <- function(sub, lab) {
e <- factor(sub$cause, levels = c(0, 1, 2), labels = c("c", "pred", "other"))
a <- survfit(Surv(sub$time, e) ~ 1)
pc <- which(a$states == "pred")
data.frame(time = c(0, a$time), cif = c(0, a$pstate[, pc]), grp = lab)
}
cif_x_df <- rbind(grid_cif(g0, "x = 0"), grid_cif(g1, "x = 1"))
p4 <- ggplot(cif_x_df, aes(time, cif, colour = grp)) +
geom_step(linewidth = 0.9) +
scale_colour_manual(values = c("x = 0" = SAGE, "x = 1" = FOREST)) +
annotate("text", x = 2.0, y = 0.33, label = "x = 1", colour = FOREST, hjust = 0, size = 3.4) +
annotate("text", x = 2.0, y = 0.20, label = "x = 0", colour = "#5d7a5f", hjust = 0, size = 3.4) +
labs(title = "Predation incidence by habitat", x = "years", y = "cumulative incidence of predation") +
coord_cartesian(ylim = c(0, 0.4)) +
theme_te()
two_panel(p3, p4)
```
## How reliable is the overstatement
To confirm that one minus Kaplan-Meier is systematically biased rather than merely unlucky here, repeat the first scenario a thousand times and compare both estimators against the known cumulative incidence.
```{r}
#| label: mc-c
set.seed(425699)
B <- 1000
truth <- (rate_pred / (rate_pred + rate_other)) *
(1 - exp(-(rate_pred + rate_other) * tau))
naive_mc <- aj_mc <- numeric(B)
for (b in 1:B) {
tp <- rexp(300, rate_pred); to <- rexp(300, rate_other)
tm <- pmin(tp, to, tau)
cc <- ifelse(tm == tau, 0L, ifelse(tp < to, 1L, 2L))
naive_mc[b] <- 1 - summary(survfit(Surv(tm, cc == 1) ~ 1), times = tau, extend = TRUE)$surv
e <- factor(cc, levels = c(0, 1, 2), labels = c("c", "p", "o"))
aa <- survfit(Surv(tm, e) ~ 1)
aj_mc[b] <- summary(aa, times = tau, extend = TRUE)$pstate[, which(aa$states == "p")]
}
naive_bias <- mean(naive_mc) - truth
aj_bias <- mean(aj_mc) - truth
```
Against a true cumulative incidence of `r sprintf("%.3f", truth)`, one minus Kaplan-Meier averages `r sprintf("%.3f", mean(naive_mc))`, a bias of `r sprintf("%+.3f", naive_bias)`. The Aalen-Johansen cumulative incidence averages `r sprintf("%.3f", mean(aj_mc))`, a bias of `r sprintf("%+.3f", aj_bias)`. The overstatement is a property of the estimator, not the draw.
## What to take away
With competing causes of death the right target is the cumulative incidence function, estimated by Aalen-Johansen, not one minus Kaplan-Meier. When you model a covariate, keep the two hazards distinct: the cause-specific hazard describes the rate among survivors and answers a mechanistic question ("does this habitat make predation faster?"), while the Fine-Gray subdistribution hazard tracks the cumulative incidence and answers a population question ("are more animals in this habitat actually killed by predators?"). They can disagree, and each is right about its own question. Report the cumulative incidence when the quantity of interest is how many animals a cause removes, and never read one minus Kaplan-Meier as if a competing death were a temporary pause.
## References
Aalen OO, Johansen S 1978 Scandinavian Journal of Statistics 5:141-150
Andersen PK, Geskus RB, de Witte T, Putter H 2012 International Journal of Epidemiology 41:861-870 (10.1093/ije/dyr213)
Austin PC, Lee DS, Fine JP 2016 Circulation 133:601-609 (10.1161/CIRCULATIONAHA.115.017719)
Fine JP, Gray RJ 1999 Journal of the American Statistical Association 94:496-509 (10.1080/01621459.1999.10474144)
Gray RJ 1988 Annals of Statistics 16:1141-1154 (10.1214/aos/1176350951)
Heisey DM, Patterson BR 2006 Journal of Wildlife Management 70:1544-1555 (10.2193/0022-541X(2006)70[1544:AROMTE]2.0.CO;2)
Putter H, Fiocco M, Geskus RB 2007 Statistics in Medicine 26:2389-2430 (10.1002/sim.2712)
## Related tutorials
- [Parametric survival and the AFT model](../parametric-survival-aft/)
- [Lost collar signals and informative censoring](../lost-signals-and-informative-censoring/)
- [Checking a survival model](../checking-a-survival-model/)
- [Cox proportional hazards and the PH assumption](../cox-proportional-hazards/)
