Competing Risks Analysis

Competing risks analysis addresses situations where a subject is at risk of experiencing more than one type of event, but only one event can occur first and doing so removes the subject from further observation. A subject “competes” among several possible failure causes.

Classic examples:

A patient may die from cancer, heart disease, or another cause; the first to occur ends their observation period.
A mechanical component may fail by fatigue, corrosion, or overload; which failure mode occurs first determines both the failure time and its cause.
A customer may churn, upgrade, or downgrade; the event that happens first changes the analysis for the remaining outcomes.

Competing risks require special treatment because treating the competing events as ordinary independent censoring gives a quantity that cannot be interpreted as a real-world probability. A deeper subtlety is the identifiability problem: from competing-risks data alone the marginal (net, latent) distribution of each cause — the distribution that would be seen if the other causes were removed — cannot be identified without an untestable assumption about the dependence between causes [Tsiatis1975cr]. This is why the observable, well-defined target is the cumulative incidence function rather than a marginal cause-specific survival.

This page builds the theory from first principles: what the data look like, the two quantities that can be estimated (the cause-specific hazard and the cumulative incidence function), why the tempting shortcut of “one minus Kaplan-Meier with the other causes censored” is wrong, and then the non-parametric, parametric and regression estimators SurPyval provides. Every method described here has a runnable example on the Competing Risks SurPyval Modelling page, and the full API is on the Competing Risks reference page.

Relationship to Univariate Analysis

Standard survival methods (Kaplan-Meier, parametric MLE) applied to a single cause — ignoring others — estimate the cause-specific survival function. Naïvely applying KM while censoring the competing events yields a quantity that cannot be interpreted as the probability of experiencing the event in the real world because the competing events are not truly independent censoring mechanisms. The subsections below make this precise.

The data: a time and a cause

In ordinary (univariate) survival analysis each unit contributes a time $x_i$ and a flag saying whether that time is an observed failure or a censoring time. Competing-risks data adds one more piece of information to each observed failure: its cause (also called the failure mode or event type). Write

$T$ for the time of the first event of any kind,
$K \in \{1, \dots, m\}$ for the cause of that event, where $m$ is the number of distinct causes,
$C$ for the censoring time, so that we observe $x = \min(T, C)$ and, if $T \leq C$, the cause $K$.

A censored unit therefore has no cause: we only know that none of the $m$ events had happened by $x$. In SurPyval the cause is passed as the array e (any hashable labels: integers, strings, …), and a missing cause (None or NaN) marks a censored row.

Throughout, censoring is assumed to be independent: $C$ carries no information about $(T, K)$ (for the regression models, conditionally on the covariates). This is the same assumption the Kaplan-Meier estimator makes; it is about the censoring, not about the causes, and it says nothing about whether the causes are independent of each other.

The cause-specific hazard

The natural building block is the rate at which cause $k$ strikes units that are still event-free. The cause-specific hazard is

\[h_k(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t,\; K = k \mid T \geq t)}{\Delta t}, \qquad k = 1, \dots, m.\]

Intuitively $h_k(t)\,\Delta t$ is the probability that a unit which has survived everything up to $t$ fails from cause $k$ in the next small interval. Because a unit can fail from only one cause at a time, the all-cause hazard is the sum of the cause-specific hazards, and so the all-cause survival function (the probability of being free of every event) is

\[h(t) = \sum_{k=1}^{m} h_k(t), \qquad S(t) = P(T > t) = \exp\!\Big(-\sum_{k=1}^{m} H_k(t)\Big),\]

where $H_k(t) = \int_0^t h_k(u)\,du$ is the cumulative cause-specific hazard. The cause-specific hazards are identifiable: they can be estimated from competing-risks data without any assumption about how the causes depend on one another, because each one only involves units that are observed to be still at risk.

The cumulative incidence function

The correct marginal quantity of interest is the Cumulative Incidence Function (CIF), also called the sub-distribution function. For cause $k$:

\[F_k(t) = P(T \leq t,\; K = k)\]

where $T$ is the event time and $K$ is the event type. In words: the probability that a unit has failed from cause $k$ by time $t$, in the real world where all the other causes are also operating. It is built from the cause-specific hazard by asking, at each instant, “is the unit still event-free, and does cause $k$ strike now?”:

\[F_k(t) = \int_0^t h_k(u)\, S(u^-)\, du .\]

The integrand $f_k^{\text{sub}}(t) = h_k(t)S(t^-)$ is the sub-distribution density (SurPyval calls it the instantaneous incidence function, iif). Two properties follow immediately:

The CIFs sum to the overall failure probability:

\[\sum_{k=1}^{m} F_k(t) = F(t) = 1 - S(t).\]
Each CIF is improper: it plateaus at the eventual probability of that cause, $F_k(\infty) = P(K = k) < 1$, rather than rising to one. A unit that fails from cause 2 can never go on to fail from cause 1.

Notice that $F_k$ depends on all the cause-specific hazards through $S$. Raising the hazard of cause 2 lowers the incidence of cause 1, even if $h_1$ is untouched, simply because fewer units survive long enough to fail from cause 1. This is the single most important idea in competing risks, and it is why the regression section below distinguishes models for the hazard from models for the incidence.

Why one minus Kaplan-Meier of one cause is wrong

The tempting shortcut is to analyse cause $k$ on its own: treat failures from every other cause as censored and compute a Kaplan-Meier curve $\hat{S}_k^{\text{KM}}$. This estimates

\[1 - S_k(t) = 1 - \exp\{-H_k(t)\},\]

which is the probability of failing from cause $k$ by $t$ in a hypothetical world where cause $k$ is the only cause in operation (the “net” probability). Two things are wrong with reading it as a real-world probability:

It overstates the incidence. Since $S(u) \leq S_k(u)$, $F_k(t) = \int_0^t h_k S \, du \leq \int_0^t h_k S_k\, du = 1 - S_k(t)$. Treating competing failures as censored pretends those units could still go on to fail from cause $k$, when in reality they never can. Summed over the causes, the “one minus KM” curves can exceed one.
It is only meaningful under an untestable assumption. The hypothetical world “with the other causes removed” corresponds to the distribution of a latent failure time $T_k$, and it equals $S_k$ only if the latent times of the different causes are independent. Competing-risks data cannot confirm or refute that independence [Tsiatis1975cr].

A small example makes the first point concrete. Six units are followed; the fourth is censored and the rest fail from cause A or B:

time $x$	1	2	3	4	5	6
cause	A	B	A	censored	B	A
at risk $r$	6	5	4	3	2	1
all-cause KM $\hat{S}(x)$	5/6	4/6	3/6	3/6	1/4	0

By time 6 every unit has either failed or been censored. The Aalen-Johansen estimator described next gives $\hat{F}_A(6) = 7/12$ and $\hat{F}_B(6) = 5/12$, which sum to $1 - \hat{S}(6) = 1$. The naive “1 - KM with B censored” for cause A instead reaches 1 at time 6 (the last unit at risk fails from A), and the naive curve for B reaches 0.6: together they claim a total failure probability of 1.6. The naive curve for A says that every unit eventually fails from A, when in fact 5/12 of the incidence was due to B. The same numbers are reproduced with code on the Competing Risks SurPyval Modelling page.

Net and crude probabilities

The two quantities above have traditional names. The CIF $F_k$ is the crude probability: it is what actually happens, and it is always identifiable. $1 - S_k$ is the net probability: what would happen if cause $k$ acted alone. In reliability engineering the net quantity is sometimes exactly what is wanted — “if we eliminate the corrosion failure mode by a design change, what will the fatigue life look like?” — but that is a counterfactual question and its answer relies on the causes acting independently. SurPyval exposes both; the CIF (cif) should be the default answer to “how likely is failure from this cause?”, and the net quantities (ff/sf with an event) should be used knowingly.

Non-Parametric CIF Estimation

The Aalen-Johansen estimator

The non-parametric estimate replaces every quantity in $F_k(t) = \int_0^t h_k(u) S(u^-) du$ with its empirical counterpart, in exactly the way the Nelson-Aalen and Kaplan-Meier estimators do for single-cause data (see Non-Parametric Estimation). Order the distinct observed times $x_1 < x_2 < \dots < x_J$ and at each $x_j$ let

$r_j$ be the number at risk (units with an observed or censored time $\geq x_j$),
$d_{k,j}$ be the number of cause-$k$ events at $x_j$, and $d_j = \sum_k d_{k,j}$ the number of events of any cause.

The empirical CIF for cause $k$ is then estimated from the cause-specific hazard rates:

\[\hat{F}_k(t) = \sum_{x_j \leq t} \hat{h}_k(x_j)\, \hat{S}(x_{j}^-) = \sum_{x_j \leq t} \frac{d_{k,j}}{r_j}\,\hat{S}(x_{j-1})\]

where $\hat{h}_k(x_j) = d_{k,j} / r_j$ is the cause-specific hazard increment at event time $x_j$, and $\hat{S}(x_j^-) = \hat{S}(x_{j-1})$ is the all-cause Kaplan-Meier survival just before $x_j$,

\[\hat{S}(t) = \prod_{x_j \leq t}\left(1 - \frac{d_j}{r_j}\right), \qquad \hat{S}(x_0) = 1 .\]

This is the Aalen-Johansen estimator [AalenJohansen1978cr]. Read each term as “the fraction of the population still event-free just before $x_j$” times “the fraction of those at risk that fail from cause $k$ at $x_j$”. With a single cause and no censoring it reduces to the empirical CDF, and with a single cause and censoring to $1 - \text{KM}$.

Two details of the formula are easy to get wrong, and both matter:

Weight by $\hat{S}(x_j^-)$, not $\hat{S}(x_j)$. The hazard at $x_j$ acts on the units alive just before $x_j$. Using the survival after the jump removes the units failing at $x_j$ before they have been counted, and makes the CIFs systematically too small (with one cause and no censoring the total incidence stops well short of one).
Use the product-limit (Kaplan-Meier) survival as the weight. Only it satisfies the telescoping identity $\hat{S}(x_{j-1}) - \hat{S}(x_j) = \hat{S}(x_{j-1})\,d_j/r_j$, which is what makes the estimated CIFs sum exactly to $1 - \hat{S}(t)$. Pairing the discrete increments $d_{k,j}/r_j$ with the exponential (Nelson-Aalen) survival $e^{-\hat{H}}$ instead inflates the CIFs and can push the total incidence above one in small samples.

The worked example above is this formula applied by hand: for cause A, $\hat{F}_A(6) = 1 \cdot \tfrac{1}{6} + \tfrac{4}{6} \cdot \tfrac{1}{4} + \tfrac{1}{4} \cdot \tfrac{1}{1} = \tfrac{7}{12}$.

What SurPyval computes

SurPyval estimates the non-parametric CIF for each cause with the CompetingRisks class, which uses this cause-specific-hazard construction directly. The incidence increments are always weighted by the Kaplan-Meier $\hat{S}(x_j^-)$, whichever method ("Nelson-Aalen", the default, or "Kaplan-Meier") is requested; the shared helper aalen_johansen_iif() implements the weighting once for the non-parametric CIF, the cause-specific Cox CIF and the per-group CIFs inside Gray’s test. So cif and iif do not depend on method.

Alongside the CIFs, the fitted model also carries the cause-specific hazard quantities. hf returns the hazard increment $d_{k,j}/r_j$ at the most recent observed time (a jump size, not a rate; zero if no cause-$k$ failure occurred there). The survival functions follow method:

`method`	`sf(t, event=k)`	`Hf(t, event=k)`
`"Nelson-Aalen"` (default)	$\exp\{-\hat{H}_k(t)\}$, with $\hat{H}_k(t) = \sum_{x_j \leq t} d_{k,j}/r_j$	$\hat{H}_k(t)$
`"Kaplan-Meier"`	$\prod_{x_j \leq t} (1 - d_{k,j}/r_j)$	$-\log$ of that product

so that sf == exp(-Hf) and ff == 1 - sf for either method. With an event these are the net quantities discussed above (cause $k$ acting alone, the other causes treated as censoring); without one they refer to all causes combined, using $d_j$ in place of $d_{k,j}$. The two methods differ little while the risk sets are large: since $e^{-h} \geq 1 - h$, the Nelson-Aalen survival is never below the product limit, and the gap grows in the tail, where the risk sets are small and each increment is large. All of these are step functions, equal to zero (or one, for survival) before the first observed time. The fitted all-cause survival at the distinct times is also stored as the attribute S.

Censoring and truncation

Censoring enters the Aalen-Johansen estimator only through the risk sets $r_j$: a unit censored at $x_i$ counts as at risk at every event time up to and including $x_i$ and then leaves, exactly as in the Kaplan-Meier estimator. The competing-risks classes in SurPyval support right censoring only. Left- and interval-censored rows (c of -1 or 2) are rejected with a ValueError: with an interval-censored failure of known cause the unit’s contribution to $F_k$ would have to be spread across the interval using all of the cause-specific hazards, which these estimators do not do.

Left truncation (delayed entry) is not an argument of the non-parametric or regression competing-risks classes. For parametric models it can be handled exactly by fitting each cause separately and assembling the result, as the next section explains.

Parametric Competing Risks

The latent-failure-time picture

A parametric competing risks model specifies a separate parametric distribution for each cause. The most intuitive way to picture it is the latent failure time model: each cause $k$ has its own clock $T_k$ with distribution $F_k^{\text{net}} = 1 - S_k$, the unit fails at $T = \min_k T_k$, and the cause is whichever clock ran out first. The overall survival function is the product of the cause-specific survival functions (assuming independent latent failure times):

\[S(t) = \prod_{k=1}^{m} S_k(t)\]

The overall density is:

\[f(t) = \sum_{k=1}^{m} f_k(t) \prod_{j \neq k} S_j(t)\]

The term inside the sum is the sub-distribution density of cause $k$: the density of a cause-$k$ failure at $t$ multiplied by the probability that no other cause has struck first. Integrating it gives the model’s cumulative incidence,

\[F_k(t) = \int_0^t f_k(u) \prod_{j \neq k} S_j(u)\, du = \int_0^t h_k(u)\, S(u)\, du ,\]

and the eventual probability of each cause is $F_k(\infty)$.

A more careful — and assumption-free — reading of the same model is that it specifies the cause-specific hazards parametrically: $h_k(t)$ is the hazard of the chosen distribution for cause $k$, and $S_k = e^{-H_k}$ is just a convenient way of writing $\exp(-\int h_k)$. Under that reading the fitted CIFs, the all-cause survival and the simulated (time, cause) pairs are all valid whether or not the latent times are independent. Only the net interpretation of each $S_k$ (“the life if the other causes were removed”) needs independence.

Why the fit separates by cause

SurPyval provides the ParametricCompetingRisks class for this model. Write $\delta_{ik} = 1$ if unit $i$ was observed to fail from cause $k$ and 0 otherwise (a censored unit has $\delta_{ik} = 0$ for every $k$). The likelihood of a unit is its cause-specific hazard at the failure time (if it failed) times the probability of surviving all causes until $x_i$:

\[L = \prod_i \Big[\prod_k h_k(x_i)^{\delta_{ik}}\Big]\, S(x_i) = \prod_{k=1}^{m} \Big[\prod_i h_k(x_i)^{\delta_{ik}}\, S_k(x_i)\Big].\]

Each bracket on the right is an ordinary right-censored likelihood for cause $k$’s distribution in which failures from cause $k$ are the observed events and everything else — failures from the other causes and genuinely censored units — is right-censored. Under the independent-latent-times assumption (or, equivalently, with the cause-specific hazards parametrised separately) the joint likelihood therefore separates, so each cause’s distribution is fitted independently by MLE with the other causes’ events treated as right-censored. This is not an approximation: it is the exact maximum-likelihood estimate of the joint model, and its log-likelihood, AIC and BIC are the sums of the per-cause values.

It is worth being clear about how this squares with the warning against “1 - KM of one cause”. Treating the other causes as censored is the correct way to estimate each cause-specific hazard; it is only wrong to read $1 - S_k$ as the probability of failing from cause $k$. The parametric model gets its CIFs right by recombining the fitted hazards through the formula for $F_k$ above.

The same factorisation extends to left truncation. A unit that enters observation at age $\tau_i$ contributes $\prod_k h_k(x_i)^{\delta_{ik}}\, S(x_i)/S(\tau_i)$, and because $S(x)/S(\tau) = \prod_k S_k(x)/S_k(\tau)$ this also splits into one left-truncated, right-censored likelihood per cause. Delayed-entry data can therefore be handled by fitting each cause’s distribution with its truncation bounds and combining the fits with ParametricCompetingRisks.from_fitted; the how-to page shows this. (Interval censoring does not factorise in this way, because the probability of a cause-$k$ failure inside an interval involves every cause’s survival across the interval.)

Computing the model’s quantities

Most model quantities are closed-form combinations of the per-cause models: the all-cause survival is $\prod_k S_k(t)$, the all-cause hazard is $\sum_k h_k(t)$, and the instantaneous incidence of cause $k$ is $f_k(t)\prod_{j\neq k} S_j(t)$. The CIF integral generally has no closed form, so SurPyval evaluates it numerically, on the cause’s own probability scale: substituting $p = F_k(u)$,

\[\mathrm{CIF}_k(t) = \int_0^{F_k(t)} \prod_{j\neq k} S_j\big(F_k^{-1}(p)\big)\,dp ,\]

an integral over a finite interval of a bounded, monotone function, which adaptive quadrature evaluates to about ten significant digits for each requested time. It stays accurate when the requested times span many orders of magnitude, when a density is infinite at zero (a Weibull shape below one) and for very heavy-tailed causes, and $F_k(\infty)$ is the same integral up to $F_k(\infty)$ (1, or the cure fraction), with no finite horizon to choose. The CIFs of the causes sum to the all-cause ff to that accuracy.

Because the parametric model is a full generative model, it can also be simulated: draw a latent time from every cause’s distribution and keep the earliest, together with its cause. The per-cause models need not be from the same family, and a cause may carry a cure (limited-failure-population) fraction, in which case some units never fail from that cause; if every cause has a cure fraction some units never fail at all and the probabilities of the causes sum to less than one.

When to prefer a parametric model: when you need smooth CIFs, extrapolation beyond the last observed failure, a compact description of each failure mode (e.g. a Weibull shape telling wear-out from random failures), or simulation for Monte-Carlo studies. The non-parametric estimator makes no shape assumption and is the right first look — and a good check of the parametric fit.

Regression: Fine-Gray and Cause-Specific PH

When covariates $Z$ (a row vector per unit: treatment, load, material, …) are available there are two fundamentally different things one can model, mirroring the two quantities defined above: the cause-specific hazards, or the cumulative incidence directly. Two main regression approaches are used in competing risks:

Cause-specific proportional hazards — fits a separate Cox or parametric PH model for each cause, with all other cause events treated as censored:

\[h_k(t \mid Z) = h_{k,0}(t)\, e^{Z \beta_k}\]

This estimates the effect of covariates on the hazard of each cause independently.

Fine-Gray sub-distribution hazards — models the effect of covariates directly on the CIF via a proportional hazards model on the sub-distribution hazard:

\[h_k^*(t \mid Z) = h_{k,0}^*(t)\, e^{Z \gamma_k}\]

This is the natural choice when the scientific question is about the probability of a cause occurring in the presence of competing risks (e.g. clinical risk scores).

SurPyval provides the FineGray and CompetingRisksProportionalHazards classes. The subsections below explain each model, how it is estimated and how to read its coefficients.

Cause-specific proportional hazards in detail

Here $h_{k,0}(t)$ is an unspecified baseline hazard for cause $k$ and $\beta_k$ a vector of coefficients for that cause (one per column of $Z$). By the same factorisation as in the parametric case, the partial likelihood separates by cause, so each $\beta_k$ is estimated by an ordinary Cox model in which cause-$k$ failures are events and all other rows are censored. $e^{\beta_{k,p}}$ is a cause-specific hazard ratio: the multiplicative change in the rate of cause $k$ among units still event-free, per unit increase of covariate $p$.

CompetingRisksProportionalHazards with model="Cox" fits one CoxPH model per cause (see Cox Proportional Hazards) and keeps each cause’s baseline cumulative hazard $\hat{\Lambda}_{k,0}$. Its tie_method argument is passed on as the Cox tie-handling method, "efron" by default as for CoxPH, and the baseline follows it (Efron’s tie correction after an Efron fit, Breslow’s estimator otherwise; the two agree when no failure times are tied). The CIF at a covariate vector $Z$ is then assembled step by step from the covariate-specific hazard increments $\Delta\hat{\Lambda}_l(x_j \mid Z) = \Delta\hat{\Lambda}_{l,0}(x_j) e^{Z\hat{\beta}_l}$, with $\Delta\hat{\Lambda}(x_j \mid Z)$ their sum over the causes:

\[\hat{F}_k(t \mid Z) = \sum_{x_j \leq t} \hat{S}(x_{j-1} \mid Z)\, \frac{\Delta\hat{\Lambda}_k(x_j \mid Z)}{\Delta\hat{\Lambda}(x_j \mid Z)} \Big(1 - e^{-\Delta\hat{\Lambda}(x_j \mid Z)}\Big), \qquad \hat{S}(t \mid Z) = \exp\Big(-\sum_{l=1}^{m} \hat{\Lambda}_{l,0}(t)\, e^{Z\hat{\beta}_l}\Big).\]

The formula shows the catch in interpreting cause-specific coefficients: the incidence of cause $k$ depends on every cause’s coefficients through $\hat{S}(t \mid Z)$. A covariate can raise the hazard of cause $k$ ($\beta_k > 0$) and yet lower its incidence, if it raises a competing cause’s hazard even more.

Each step is the Aalen-Johansen step with the transition probabilities of the matrix exponential of that step’s hazards, which is how R’s survival computes the state probabilities of a multi-state coxph; SurPyval matches it to rounding. A unit still event-free before the step fails over it with probability $1 - e^{-\Delta\hat{\Lambda}}$, and each cause takes its share of that. The increments telescope, so the cause-specific CIFs sum to exactly $1 - \hat{S}(t \mid Z)$, the model’s ff, and never exceed one, even at a covariate value far from the data, where a step’s total hazard increment can exceed one (a small risk set times a large multiplier). The plain product limit $\prod (1 - \Delta\hat{\Lambda})$ would need a clip there, and would not be the survival sf reports.

The Fine-Gray model in detail

Fine and Gray [FineGray1999cr] asked a different question: is there a hazard whose proportional-hazards model acts directly on the CIF, the way the ordinary hazard acts on the survival function? The answer is the sub-distribution hazard

\[h_k^*(t) = \lim_{\Delta t \to 0} \frac{P\big(t \leq T < t + \Delta t,\, K = k \;\big|\; T \geq t \;\text{or}\; (T < t,\, K \neq k)\big)}{\Delta t} = \frac{f_k^{\text{sub}}(t)}{1 - F_k(t)} .\]

The conditioning event is the peculiar part: units that have already failed from a competing cause are kept “at risk” for cause $k$ — they will never fail from $k$, which is exactly what makes $1 - F_k(t)$ the right denominator. Because $h_k^*$ is the hazard of the improper distribution $F_k$, we get the familiar relation $F_k(t) = 1 - \exp\{-\Lambda_k^*(t)\}$ with $\Lambda_k^* = \int h_k^*$. The proportional sub-distribution hazards model $h_k^*(t \mid Z) = h_{k,0}^*(t) e^{Z\gamma_k}$ therefore gives

\[F_k(t \mid Z) = 1 - \exp\big\{-\Lambda_{k,0}^*(t)\, e^{Z\gamma_k}\big\},\]

so a positive coefficient always raises the incidence of cause $k$, and $e^{\gamma_{k,p}}$ is a sub-distribution hazard ratio. (In the code and on the how-to page the Fine-Gray coefficients are called beta.)

Estimation with censoring. Without censoring the sub-distribution risk set at time $t$ is simply “everyone who has not failed from cause $k$ by $t$”. Under right censoring we do not know whether a unit that failed from a competing cause at $x_i < t$ would still have been under observation at $t$. Fine and Gray keep it in the risk set with the inverse-probability-of-censoring weight

\[\begin{split}w_i(t) = \begin{cases} 1 & x_i \geq t \quad \text{(still under observation)},\\[2pt] \hat{G}(t^-)/\hat{G}(x_i^-) & x_i < t \text{ and failed from a competing cause},\\[2pt] 0 & \text{otherwise (censored, or failed from cause } k \text{, before } t), \end{cases}\end{split}\]

where $\hat{G}(t)$ is the Kaplan-Meier estimate of the censoring survival function $P(C > t)$ (the censored rows play the role of “events”; everyone with $x_j \geq s$ is at risk of censoring at $s$) and $\hat{G}(t^-)$ its value just before $t$. $\hat{G}(t^-)/\hat{G}(x_i^-)$ is the estimated probability that the unit would have remained uncensored from $x_i$ to $t$. One $\hat{G}$ is estimated from the whole sample, which assumes the censoring does not depend on the covariates.

Taking $\hat{G}$ just before each time settles ties: when a censoring and an event happen at the same recorded time, the event is taken to come first, so the censorings at $t$ do not yet reduce the weights at $t$, and the censorings at $x_i$ do not count against a unit that failed from a competing cause at $x_i$. This is the weight of R’s cmprsk::crr, $\hat{G}(t^-)/\hat{G}(x_i^-)$ with $\hat{G}$ the reverse Kaplan-Meier (crr reads its censoring curve at ftime*(1 - 100*.Machine$double.eps)), and SurPyval’s coefficients and likelihood reproduce a transcription of crr’s code on data with such ties. When no censoring time equals an event time the left limits are simply $\hat{G}(t)$ and $\hat{G}(x_i)$. (SurPyval 0.20 and earlier used $\hat{G}(t)/\hat{G}(x_i)$ and so differed slightly from crr on tied data. R’s survival::finegray also evaluates $\hat{G}$ just before each time, but it removes the failures at a tied time from the censoring risk set there, so on tied data its weights can differ slightly from crr’s and SurPyval’s.) The coefficients maximise the weighted partial log-likelihood

\[\ell(\gamma) = \sum_{i:\,K_i = k} n_i \Big[ Z_i\gamma - \log \sum_j w_j(x_i)\, n_j\, e^{Z_j\gamma} \Big],\]

with $n_i$ the count of each row. (Tied event times are handled by this Breslow form: each tied event sees the same weighted risk set.) SurPyval’s FineGray maximises this with BFGS (exact gradients by automatic differentiation), reports standard errors from the inverse of the Hessian of $\ell$ at the optimum, and estimates the baseline by the Breslow-type increments $\Delta\hat{\Lambda}_{k,0}^*(t) = d_k(t) / \sum_j w_j(t)\, n_j\, e^{Z_j\hat{\gamma}}$. The fitted CIF is a step function, flat before the first and after the last cause-$k$ event time. In time order the weighted risk set at $t$ is a suffix of the rows ($x_j \geq t$) plus $\hat{G}(t^-)$ times a prefix of the competing failures, each weighted $1/\hat{G}(x_j^-)$, so every sum is a cumulative sum and a fit takes time and memory linear in the number of rows (after a sort): about a second at $10^5$ rows and three covariates. Those standard errors are the model-based (inverse information) ones; the robust sandwich variance that Fine and Gray derived to account for the estimated weights is not implemented, so treat the reported se and p_values as approximate.

Choosing between them

The two models answer different questions and are complementary rather than competing [Latouche2013cr]:

“What drives the rate of this failure mode?” — aetiology, physics of failure, mechanism. Use cause-specific hazards. The coefficients have the usual hazard-ratio meaning among units still at risk.
“Who is most likely to end up failing from this mode?” — prognosis, risk scores, warranty or spares forecasting. Use Fine-Gray. The coefficients map monotonically onto the incidence.

The two coincide only in special cases (for instance, when the covariate does not affect the competing causes at all). Two practical caveats of Fine-Gray: the model is fitted one cause at a time, and separate Fine-Gray fits for different causes are not constrained to be mutually consistent (their CIFs can sum to more than one); and the sub-distribution hazard itself has no direct physical meaning, so interpret the model through its CIFs. In SurPyval, CompetingRisksProportionalHazards with model="Fine-Gray" fits the Fine-Gray model for every cause at once.

Comparing incidence across groups: Gray’s test

The log-rank test compares survival curves; its competing-risks analogue is Gray’s test [Gray1988cr], which compares the cumulative incidence functions of a chosen cause across groups. The distinction matters. A cause-specific log-rank compares the cause-specific hazards — the instantaneous rate of the cause among those still at risk — whereas Gray’s test compares the CIFs themselves, i.e. the actual incidence of the cause in a population that is also being depleted by the competing causes.

Gray’s test achieves this by modifying the risk set. Instead of removing subjects who fail from a competing cause (as a cause-specific analysis would), it keeps them in the subdistribution risk set, which Gray estimates within each group $g$ as

\[R_g(t) = \frac{Y_g(t)\,\{1 - \hat{F}_g(t^-)\}}{\hat{S}_g(t^-)},\]

with $Y_g$ the number at risk in the group, $\hat{F}_g$ its Aalen-Johansen CIF of the cause and $\hat{S}_g$ its all-cause Kaplan-Meier survival. Since $Y_g(t)/\hat{S}_g(t^-)$ estimates the group’s size times its probability of still being uncensored, this is the inverse-probability-of-censoring weighted risk set of the Fine-Gray model (below), with a censoring distribution estimated separately in every group: subjects who have already failed from a competing cause continue to count — with a weight that decays as the group’s censoring accrues — which is precisely what makes the comparison one of incidence rather than of instantaneous rate. Under the null hypothesis of equal CIFs the resulting statistic is approximately $\chi^2$ distributed with $G - 1$ degrees of freedom for $G$ groups. Reach for it when the question is “how many fail of this cause”, and for the cause-specific log-rank when the question is “how fast”.

The statistic

SurPyval computes Gray’s statistic [Gray1988cr]. At each distinct time $\tau$ at which the cause of interest occurs, let $R = \sum_g R_g$, $d_g(\tau)$ be the number of cause-of-interest events in group $g$ and $d = \sum_g d_g$. Under the null hypothesis that every group has the same CIF, the events should be shared out in proportion to the subdistribution risk sets, so the test accumulates observed-minus-expected counts

\[U_g = \sum_\tau W(\tau)\Big[d_g(\tau) - d(\tau)\frac{R_g(\tau)}{R(\tau)}\Big].\]

Its variance is not the hypergeometric variance of an ordinary log-rank test: the risk sets $R_g$ are themselves estimates, built from each group’s incidence and survival, so failures from the competing causes add variability too. Gray linearises $U$ in each group’s counting-process martingales, of the cause and of the competing causes, under the null hypothesis. With $n_g(\tau) = Y_g(\tau)/\hat{S}_g(\tau^-)$ (the group size times its censoring survival), $n = \sum_g n_g$, and $\hat{F}^0$ the pooled CIF of the cause, which steps by $d(\tau)/n(\tau)$, this gives

\[V_{gh} = \sum_r \sum_u \Big[A_{gr}(u)\,A_{hr}(u)\,w_r(u) + B_{gr}(u)\,B_{hr}(u)\,w^{c}_r(u)\Big],\]

\[\begin{split}A_{gr}(u) &= a_{gr}(u) + \Big(1 - \frac{1 - \hat{F}^0(u)}{\hat{S}_r(u)}\Big) Q_{gr}(u), \qquad B_{gr}(u) = \frac{1 - \hat{F}^0(u)}{\hat{S}_r(u)}\,Q_{gr}(u), \\ Q_{gr}(u) &= \sum_{\tau > u} a_{gr}(\tau)\, \frac{d(\tau)/n(\tau)}{1 - \hat{F}^0(\tau^-)}, \qquad a_{gr}(\tau) = W(\tau)\,n_g(\tau)\Big[\delta_{gr} - \frac{n_r(\tau)}{n(\tau)}\Big],\end{split}\]

where the martingale variances are the failures from the cause that group $r$ is expected to have under the null hypothesis, $w_r(u) = d(u)/\{n(u)\,n_r(u)\}$, and those it had from the competing causes, $w^c_r(u) = d^c_r(u)\,\{\hat{S}_r(u^-)/Y_r(u)\}^2$, each times a correction for tied failures ($1 - (d - 1)/(n \hat{S}_r(u^-) - 1)$ and $1 - (d^c_r - 1)/(Y_r - 1)$). Dropping one group to make $V$ invertible, the statistic $U^{\top} V^{-1} U$ is referred to a $\chi^2$ distribution with (number of groups $- 1$) degrees of freedom. The weight $W(\tau) = \{1 - \hat{F}^0(\tau^-)\}^{\rho}$; the default $\rho = 0$ (every event time weighted equally) is the standard test, and $\rho > 0$ down-weights late event times, making the test more sensitive to differences in early incidence.

This is the test R’s cmprsk::cuminc computes, step for step: SurPyval’s statistic agrees with cuminc’s to rounding, with ties and for any $\rho$. The subdistribution risk set $R_g(t)$ is formed from the number at risk at $t$ and the survival just before $t$, so a censoring tied with a failure counts after it, the same ordering as the Fine-Gray weights above. Because every group’s risk set carries its own censoring estimate, the test does not need the groups to be censored alike; the how-to page checks its size by simulation, with groups censored very differently.

A useful way to see the difference from a cause-specific log-rank: imagine two groups with identical cause-1 hazards but a much larger cause-2 hazard in the second group. The cause-specific log-rank for cause 1 sees no difference (the rates are equal), while Gray’s test correctly reports that far fewer units in the second group ever fail from cause 1. Both answers are right; they answer different questions. The how-to page runs exactly this experiment.

For worked examples of estimating cumulative incidence functions, fitting the Fine-Gray and cause-specific proportional hazards models, and comparing groups with Gray’s test, see the Competing Risks SurPyval Modelling page. For repeated events with several failure modes (the recurrent-events analogue of competing risks) see Recurrent Event Analysis.