Competing Risks SurPyval Modelling

This page shows how to use SurPyval’s competing risks classes. For the theoretical background see Competing Risks Analysis; for the complete list of arguments and methods see the Competing Risks API reference.

Note

Every example on this page is executed when the documentation is built, so the outputs shown are produced by the installed version of surpyval.

SurPyval’s competing-risks tools, and the question each one answers:

Tool	Answers
CompetingRisks	Non-parametric (Aalen-Johansen) cumulative incidence of each cause.
ParametricCompetingRisks	One distribution per cause; smooth CIFs, extrapolation, simulation.
gray_test	Does the cumulative incidence of a cause differ between groups?
FineGray	How do covariates change the cumulative incidence of one cause?
CompetingRisksProportionalHazards	How do covariates change each cause-specific hazard (model="Cox"), or the incidence of every cause (model="Fine-Gray")?

The classes live in surpyval.univariate.competing_risks; gray_test is also available at the top level as surpyval.gray_test. All of them accept right-censored data only. For repeated events with several failure modes see the cause-specific MCF and NHPP models in Recurrent Event Modelling with SurPyval.

Standard imports used throughout this page:

import surpyval as surv
import numpy as np
from matplotlib import pyplot as plt

Fitting a Competing Risks Model

The CompetingRisks class estimates a non-parametric cumulative incidence function (CIF) for each failure cause with the Aalen-Johansen estimator. Pass the observed times, a cause indicator, and optional censoring flags.

Competing-risks data format

Every competing-risks fit takes the same core arrays:

• x – the observed time of each row (failure or censoring time);

• e – the cause of each row. Any hashable labels work (integers, strings, …). A censored row has no cause: use None (or NaN, or a blank cell in a DataFrame);

• c – optional censoring flags, 0 observed and 1 right-censored. When c is omitted it is derived from e: a missing cause means censored. If you do pass c, every row with c == 1 must have a missing cause and every other row must have one, otherwise a ValueError explains the mismatch;

• n – optional counts, when each row stands for several identical units.

Left- and interval-censored rows (c of -1 or 2) are rejected.

lifelines, scikit-survival and R’s cmprsk code the causes as one integer column with 0 for a censored row. Here 0 is a cause like any other (you may number causes from 0), so convert such data with e = np.where(np.asarray(e) == 0, None, e). Because the two codings cannot be told apart, a fit whose labels are numbers including 0, with none missing and no c, warns and says so; passing c (np.zeros(len(e)) when nothing is censored) silences it.

The smallest possible example is the six-unit data set worked by hand on the Competing Risks Analysis page. The fourth unit is censored:

from surpyval.univariate.competing_risks import CompetingRisks

x = [1, 2, 3, 4, 5, 6]
e = ["A", "B", "A", None, "B", "A"]

small = CompetingRisks.fit(x, e)
print("CIF of A at t=6:", small.cif(6, "A"), "(7/12 = %.4f)" % (7 / 12))
print("CIF of B at t=6:", small.cif(6, "B"), "(5/12 = %.4f)" % (5 / 12))

CIF of A at t=6: 0.5833333333333334 (7/12 = 0.5833)
CIF of B at t=6: 0.4166666666666667 (5/12 = 0.4167)

and “one minus Kaplan-Meier with the other cause censored” gives the misleading answers discussed on the theory page – 1 for A and 0.6 for B, a total “probability” of 1.6:

for k in ["A", "B"]:
    c_k = [0 if ei == k else 1 for ei in e]   # other cause -> censored
    naive = surv.KaplanMeier.fit(x, c=c_k)
    print(k, "naive 1 - KM at t=6:", naive.ff(6))

A naive 1 - KM at t=6: 1.0
B naive 1 - KM at t=6: 0.6

When each row stands for several identical units, give the counts in n rather than repeating rows. Doubling every row of the six-unit data gives the same CIFs, because the estimator only uses the proportions \(d_{k,j}/r_j\). The iif method returns the individual increments \(\hat{S}(x_{j-1})\,d_{k,j}/r_j\), here the three terms \(\tfrac{1}{6}, \tfrac{1}{6}, \tfrac{1}{4}\) of the hand calculation on the theory page:

doubled = CompetingRisks.fit(x, e, n=[2] * 6)
print("CIF of A at t=3, 6     :", doubled.cif([3, 6], "A"))
print("increments at t=1, 3, 6:", doubled.iif([1, 3, 6], "A"))

CIF of A at t=3, 6     : [0.33333333 0.58333333]
increments at t=1, 3, 6: [0.16666667 0.16666667 0.25      ]

Non-parametric cumulative incidence

A more realistic example: components fail either by wear-out (a Weibull latent life) or by random shocks (an exponential latent life), whichever comes first, and some are removed from test (censored) before failing. Because we simulate the data we know the true CIFs and can check the estimates:

rng = np.random.default_rng(0)
N = 300
t_wear = 100 * rng.weibull(3.0, N)       # Weibull(alpha=100, beta=3)
t_shock = rng.exponential(150, N)        # Exponential, mean 150
t_cens = rng.uniform(0, 200, N)          # removal from test

x = np.minimum.reduce([t_wear, t_shock, t_cens])
e = np.where(x == t_cens, None,
             np.where(t_wear < t_shock, "wear", "shock")).astype(object)

model = CompetingRisks.fit(x, e)         # c is derived from e
print(model)
print("event index map:", model.event_idx_map)

Competing Risk model with events:
['shock', 'wear']

event index map: {'shock': 0, 'wear': 1}

Causes are sorted and mapped to internal indices (event_idx_map); you always refer to a cause by its own label. Once fitted, you can query the CIF for each cause:

t_plot = np.linspace(0, 200, 400)
for k in ["wear", "shock"]:
    plt.step(t_plot, model.cif(t_plot, event=k), where="post", label=k)
plt.xlabel('Time')
plt.ylabel('Cumulative Incidence')
plt.legend()
plt.title('Competing Risks CIF by Cause')

Text(0.5, 1.0, 'Competing Risks CIF by Cause')

The CIFs of all causes add up to the all-cause failure probability. The Aalen-Johansen increments are weighted with the Kaplan-Meier all-cause survival, so the sum equals one minus the ordinary Kaplan-Meier curve of the data with every cause treated as a failure – exactly, not approximately:

t = np.array([25.0, 50.0, 100.0, 150.0])
total = model.cif(t, "wear") + model.cif(t, "shock")
all_cause_km = surv.KaplanMeier.fit(x, c=(e == None).astype(int))
print("sum of CIFs   :", np.round(total, 4))
print("1 - KM (all)  :", np.round(all_cause_km.ff(t), 4))

sum of CIFs   : [0.1485 0.3722 0.7605 0.9567]
1 - KM (all)  : [0.1485 0.3722 0.7605 0.9567]

Now compare the CIF of wear-out with the naive “1 - Kaplan-Meier with shocks censored” curve, and with the truth. For this simulation the true CIF of wear is \(\int_0^t f_{\text{wear}}(u)\, S_{\text{shock}}(u)\, du\):

from scipy.integrate import quad

def true_cif_wear(t):
    dens = lambda u: surv.Weibull.df(u, 100, 3.0) * np.exp(-u / 150)
    return quad(dens, 0, t)[0]

naive_wear = surv.KaplanMeier.fit(x, c=(e != "wear").astype(int))

plt.step(t_plot, model.cif(t_plot, "wear"), where="post",
         label="Aalen-Johansen CIF")
plt.step(t_plot, naive_wear.ff(t_plot), where="post",
         label="naive 1 - KM (shocks censored)")
plt.plot(t_plot, [true_cif_wear(ti) for ti in t_plot], "k--",
         label="true CIF")
plt.xlabel("Time")
plt.ylabel("Probability of wear-out failure")
plt.legend()

<matplotlib.legend.Legend at 0x76f8df902120>

The naive curve climbs towards one because it pretends that a component destroyed by a shock could still have worn out later. The Aalen-Johansen estimate tracks the true incidence, which levels off at the probability that wear-out, rather than a shock, is what ends a component’s life.

What the fitted model returns

All methods take the query times first and a cause label as event. The CIF methods need a cause; the others treat event=None (the default) as “all causes combined”. Every function is a step function that is zero (or one, for sf) before the first observed time.

Method	Returns
cif(x, event)	Aalen-Johansen cumulative incidence \(\hat{F}_k(x)\) – the real-world probability of having failed from event by x.
iif(x, event)	The CIF jump at the most recent observed time at or before x (zero if no failure from event happened at that time).
hf(x, event=None)	The Nelson-Aalen hazard increment \(d_{k,j}/r_j\) at the most recent observed time at or before x (a jump size, not a rate).
Hf(x, event=None)	Cumulative (cause-specific) hazard \(\hat{H}_k(x)\).
sf(x, event=None) / ff(x, event=None)	The survival and 1 - sf, by the method the model was fitted with: \(e^{-\hat{H}_k(x)}\) (Nelson-Aalen, the default) or the product limit \(\prod (1 - d_{k,j}/r_j)\) (Kaplan-Meier). With an event these are the net quantities (the cause acting alone); with event=None they are all-cause.
df(x, event=None)	hf * sf.

print("CIF wear         :", np.round(model.cif(t, "wear"), 4))
print("net ff wear      :", np.round(model.ff(t, event="wear"), 4))
print("cum. hazard wear :", np.round(model.Hf(t, event="wear"), 4))
print("all-cause sf     :", np.round(model.sf(t), 4))

CIF wear         : [0.0179 0.0792 0.3338 0.5018]
net ff wear      : [0.0195 0.1005 0.5334 0.8829]
cum. hazard wear : [0.0197 0.1059 0.7623 2.1449]
all-cause sf     : [0.8517 0.6286 0.2417 0.0497]

Warning

ff(x, event=k) is the net failure probability: “one minus Kaplan-Meier with the other causes censored” (exactly that with method="Kaplan-Meier", its Nelson-Aalen version by default). It answers “what if cause k were the only cause?”, which is only meaningful if the causes act independently. For “how likely is a failure from cause k?” always use cif.

The method argument ("Nelson-Aalen", the default, or "Kaplan-Meier") selects the survival estimator that sf, ff and Hf report (Hf is -log sf, so the three stay consistent); the all-cause estimate is also stored as model.S at the distinct observed times model.x. It does not change cif: the Aalen-Johansen weights are always the product-limit (Kaplan-Meier) survival, the only one for which the CIFs sum to the all-cause failure probability.

km_model = CompetingRisks.fit(x, e, how="Kaplan-Meier")
print("same CIFs          :", np.allclose(km_model.cif(t, "wear"), model.cif(t, "wear")))
print("sf (Kaplan-Meier)  :", np.round(km_model.sf(t), 4))
print("sf (Nelson-Aalen)  :", np.round(model.sf(t), 4))
print("CIFs sum to 1 - KM :", np.allclose(
    km_model.cif(t, "wear") + km_model.cif(t, "shock"), km_model.ff(t)))

same CIFs          : True
sf (Kaplan-Meier)  : [0.8515 0.6278 0.2395 0.0433]
sf (Nelson-Aalen)  : [0.8517 0.6286 0.2417 0.0497]
CIFs sum to 1 - KM : True

Like the single-cause step estimates, every function starts at its initial value before the first observed time (sf 1, the others 0) and holds its last value after the last one, however far from the data. To give the estimate an explicit range instead, call set_support(lower, upper): the functions keep that convention only within [lower, upper] and are nan outside it. The bounds must contain the observed times; lower may be negative and either bound infinite. set_support returns the model, and to_dict saves the bounds with it.

bounded = CompetingRisks.fit(x, e).set_support(0, 2 * x.max())
q = [-1, 0, 1.5 * x.max(), 3 * x.max()]
print("bounded CIF wear :", np.round(bounded.cif(q, "wear"), 4))
print("CIF wear         :", np.round(model.cif(q, "wear"), 4))

bounded CIF wear : [   nan 0.     0.5451    nan]
CIF wear         : [0.     0.     0.5451 0.5451]

Data held in a pandas DataFrame can be passed with fit_from_df, naming the time and cause columns (and optionally c_col and n_col). The frame is kept on the model as source_df:

import pandas as pd

df = pd.DataFrame({"time": x, "mode": e})
model_df = CompetingRisks.fit_from_df(df, x_col="time", e_col="mode")
np.allclose(model_df.cif(t, "shock"), model.cif(t, "shock"))

True

Parametric cumulative incidence

ParametricCompetingRisks fits one parametric distribution to each cause (the fit separates by cause, as explained on the theory page). The default is a Weibull for every cause; pass a {cause: distribution} dict to choose per cause. Here we use the distributions the data were simulated from:

from surpyval.univariate.competing_risks import ParametricCompetingRisks

pmodel = ParametricCompetingRisks.fit(
    x, e, dist={"wear": surv.Weibull, "shock": surv.Exponential}
)
print(pmodel)
for k in pmodel.causes:
    print(k, pmodel.models[k].params)

Parametric Competing Risks SurPyval Model
=========================================
Causes              : ['shock', 'wear']
Cause distributions : shock: Exponential, wear: Weibull
shock [0.00652956]
wear [112.40369541   2.89653776]

The fitted per-cause models are ordinary SurPyval models in pmodel.models (true values: Weibull \(\alpha = 100, \beta = 3\) and an exponential rate of \(1/150 \approx 0.0067\); the sample has only 83 wear-out failures, most of them below 100, so the Weibull scale carries the most sampling error). The smooth parametric CIFs sit on top of the non-parametric steps:

for k, colour in [("wear", "C0"), ("shock", "C1")]:
    plt.step(t_plot, model.cif(t_plot, k), where="post", color=colour,
             alpha=0.5, label=k + " (Aalen-Johansen)")
    plt.plot(t_plot, pmodel.cif(t_plot, k), color=colour,
             label=k + " (parametric)")
plt.xlabel("Time")
plt.ylabel("Cumulative incidence")
plt.legend()

<matplotlib.legend.Legend at 0x76f8dd1d4a70>

Unlike the step function, the parametric model can be evaluated beyond the last observation. probability_of_cause gives \(F_k(\infty)\), the long-run share of failures from each cause (these sum to one unless a cause has a cure fraction):

for k in pmodel.causes:
    print(k, "eventual probability: %.3f" % pmodel.probability_of_cause(k))

shock eventual probability: 0.465

wear eventual probability: 0.535

The other methods follow the same conventions as CompetingRisks: hf/Hf take an optional event (cause-specific) and otherwise sum over causes; sf and ff are the all-cause survival and failure probability; iif(x, event) is the sub-distribution density; and cif(x) without an event is the all-cause incidence 1 - sf. cif returns a float for a scalar time:

print("CIF wear at 100    :", pmodel.cif(100.0, "wear"))
print("all-cause ff at 100:", pmodel.ff(100.0))
print("hazard of shock    :", pmodel.hf(100.0, event="shock"))

CIF wear at 100    : 0.32488468151932903
all-cause ff at 100: 0.7447889974470009
hazard of shock    : 0.006529558412992116

Because the joint likelihood is the product of the per-cause likelihoods, neg_ll and aic are sums over causes. bic is the joint criterion 2 neg_ll + K ln(n), with K the parameters of every cause and n the observed failures of any cause (here 189: 83 wear and 106 shock), not the sum of the causes’ own BICs. Either can be used to compare candidate distributions – for example, whether the shock mode needs a Weibull rather than an exponential:

all_weibull = ParametricCompetingRisks.fit(x, e)   # Weibull for both causes
print("Weibull + Exponential AIC: %.1f, BIC: %.1f"
      % (pmodel.aic(), pmodel.bic()))
print("Weibull + Weibull     AIC: %.1f, BIC: %.1f"
      % (all_weibull.aic(), all_weibull.bic()))
print("fitted shock shape:", all_weibull.models["shock"].params[1])

Weibull + Exponential AIC: 2213.3, BIC: 2223.1
Weibull + Weibull     AIC: 2215.3, BIC: 2228.3
fitted shock shape: 1.012735551113601

The shock mode’s fitted Weibull shape is close to one (an exponential), and the Weibull + Exponential model has the lower AIC and BIC: the extra shape parameter is not worth its cost. ParametricCompetingRisks also has a fit_from_df(df, x_col, e_col, c_col=None, n_col=None, dist=Weibull, how="MLE"); how is passed to each cause’s distribution fit.

Assembling a model from separately fitted causes

ParametricCompetingRisks.from_fitted combines already-fitted single-cause models – a {cause: model} dict, or a list whose causes become 0, 1, ... – into one competing-risks model. Each model can use any SurPyval option (a different family, an offset, a limited-failure or zero-inflated fit, …). The one rule: fit each model to the cause-specific view of the data, with that cause’s failures observed and every other row right-censored.

This is also how to handle delayed entry (left truncation), which the competing-risks fit methods do not take as an argument. Suppose the components were only put under observation at a random age entry, so that units which failed before entering were never seen. Ignoring the entry ages biases the fit; fitting each cause with tl=entry and assembling the result is the exact maximum-likelihood fit (the truncated likelihood factorises by cause too):

truth = ParametricCompetingRisks.from_fitted({
    "wear": surv.Weibull.from_params([100.0, 3.0]),
    "shock": surv.Exponential.from_params([1 / 150]),
})
sim = truth.random(2000, random_state=4)
entry = np.random.default_rng(5).uniform(0, 80, 2000)
seen = sim["x"] > entry                      # only survivors to entry
x_t, e_t, tl = sim["x"][seen], sim["e"][seen], entry[seen]

naive = ParametricCompetingRisks.fit(
    x_t, e_t, dist={"wear": surv.Weibull, "shock": surv.Exponential}
)

per_cause = {}
for k, dist in [("wear", surv.Weibull), ("shock", surv.Exponential)]:
    c_k = np.where(e_t == k, 0, 1)           # cause-specific view
    per_cause[k] = dist.fit(x_t, c=c_k, tl=tl)
adjusted = ParametricCompetingRisks.from_fitted(per_cause)

print("true shock rate    : %.5f" % (1 / 150))
print("ignoring entry     : %.5f" % naive.models["shock"].params[0])
print("with truncation    : %.5f" % adjusted.models["shock"].params[0])
print("P(wear) true / naive / adjusted: %.3f / %.3f / %.3f" % (
    truth.probability_of_cause("wear"),
    naive.probability_of_cause("wear"),
    adjusted.probability_of_cause("wear"),
))

true shock rate    : 0.00667
ignoring entry     : 0.00401
with truncation    : 0.00678

P(wear) true / naive / adjusted: 0.564 / 0.689 / 0.556

Early failures – which are disproportionately shocks – were never observed, so the naive fit underestimates the shock rate and overstates the share of wear-out. The truncation-aware fit recovers both.

Simulating competing-risks data

The previous example already used random. It draws a latent time from each cause’s distribution (by inverse-transform sampling through its qf) and keeps the earliest, returning a structured array with fields x (time) and e (cause). A unit whose latent times are all infinite (possible with cure models) never fails and is returned with x = inf and cause None. Censoring is not part of the model, so add it yourself when simulating a study:

draws = truth.random(5, random_state=1)
print(draws)
print(draws["x"], draws["e"])

[( 81.95590243, 'wear') (120.70266429, 'wear') ( 23.35070755, 'shock')
 ( 92.73905296, 'wear') ( 30.34635267, 'wear')]
[ 81.95590243 120.70266429  23.35070755  92.73905296  30.34635267] ['wear' 'wear' 'shock' 'wear' 'wear']

A small Monte-Carlo study: how well does a 150-unit, censored test estimate the long-run share of wear-out failures?

rng = np.random.default_rng(2)
estimates = []
for _ in range(30):
    s = truth.random(150, random_state=rng.integers(1_000_000))
    cens = rng.uniform(0, 200, 150)
    x_sim = np.minimum(s["x"], cens)
    e_sim = np.where(s["x"] <= cens, s["e"], None)
    fit = ParametricCompetingRisks.fit(
        x_sim, e_sim, dist={"wear": surv.Weibull, "shock": surv.Exponential}
    )
    estimates.append(fit.probability_of_cause("wear"))
print("true P(wear) = %.3f;  estimates: mean %.3f, sd %.3f" % (
    truth.probability_of_cause("wear"), np.mean(estimates), np.std(estimates)))

true P(wear) = 0.564;  estimates: mean 0.561, sd 0.043

Saving and loading competing-risks models

Every competing-risks model – CompetingRisks, ParametricCompetingRisks, the Fine-Gray model returned by FineGray.fit and CompetingRisksProportionalHazards – can be serialised with to_dict/to_json and restored with the class’s from_dict/from_json or with the package-level surpyval.from_dict/surpyval.from_json, which work out the class from the dictionary. The restored model reproduces every prediction exactly:

import json

restored = surv.from_dict(json.loads(json.dumps(model.to_dict())))
print(type(restored).__name__,
      np.allclose(restored.cif(t, "wear"), model.cif(t, "wear")))

restored_p = surv.from_dict(pmodel.to_dict())
print(type(restored_p).__name__, restored_p.cif(100.0, "wear"))

CompetingRisks True
ParametricCompetingRisks 0.32488468151932903

Comparing incidence across groups: Gray’s test

Having estimated a cumulative incidence function for each group, surpyval.gray_test tests whether the CIFs differ for a chosen cause. It is the competing-risks analogue of the log-rank test, but with an important distinction: where a cause-specific log-rank compares the instantaneous hazards of a cause, Gray’s test compares the incidence — the CIFs directly. It does this by keeping competing-cause failures in the subdistribution risk set, rather than removing them, with each group’s risk set estimated from that group’s own data (and so its own censoring distribution). Reach for it when the clinical or engineering question is “how many fail of this cause”, not “how fast”.

Pass the observed times x, the per-observation cause label e, the group label, and the cause of interest; optional n gives row counts. Censored rows follow the same rules as for the model classes (see Competing-risks data format): a missing cause (None or NaN) marks a censored row, and an optional c must agree with it.

Here two groups have genuinely different cause-1 incidence:

from surpyval import gray_test

rng = np.random.default_rng(7)

def simulate(n, p_cause1):
    is1 = rng.random(n) < p_cause1
    t = rng.exponential(6.0, n)
    return t, np.where(is1, 1, 2)         # causes labelled 1 and 2

x_a, e_a = simulate(300, 0.35)
x_b, e_b = simulate(300, 0.60)            # higher cause-1 incidence
x = np.concatenate([x_a, x_b])
e = np.concatenate([e_a, e_b])
group = np.array([0] * 300 + [1] * 300)

result = gray_test(x, e, group, event=1)
print('statistic = %.2f   df = %d   p = %.3g'
      % (result.statistic, result.df, result.p_value))

statistic = 47.93   df = 1   p = 4.41e-12

The tiny p-value correctly flags the difference in cause-1 incidence. The result is a named tuple (statistic, df, p_value, cause, groups); df is the number of groups minus one, so more than two groups are compared in one test. The statistic is Gray’s (1988): an observed-minus-expected count of cause-1 failures on the subdistribution risk sets, with Gray’s asymptotic variance (the Competing Risks Analysis page gives the formulas).

Calibration

On data where the groups share the same incidence the test is calibrated, returning p-values spread over [0, 1] — including under censoring, and when the groups are censored differently. A quick check: simulate many pairs of groups with identical cause-specific hazards and independent censoring, and count how often p < 0.05. The second pair of groups is censored very differently, exponentially with means 2 and 50:

rng = np.random.default_rng(3)

def simulate_cr(n, h1, h2, cens_mean):
    t1 = rng.exponential(1 / h1, n)       # latent cause-1 time
    t2 = rng.exponential(1 / h2, n)       # latent cause-2 time
    cz = rng.exponential(cens_mean, n)    # censoring time
    x = np.minimum.reduce([t1, t2, cz])
    e = np.where(x == cz, None, np.where(t1 < t2, 1, 2)).astype(object)
    return x, e

def rejection_rate(cens_mean_0, cens_mean_1, n=100, reps=200):
    p_values = []
    for _ in range(reps):
        x0, e0 = simulate_cr(n, 0.1, 0.2, cens_mean_0)
        x1, e1 = simulate_cr(n, 0.1, 0.2, cens_mean_1)
        res = gray_test(np.concatenate([x0, x1]),
                        np.concatenate([e0, e1]),
                        np.repeat([0, 1], n), event=1)
        p_values.append(res.p_value)
    return np.mean(np.array(p_values) < 0.05)

same = rejection_rate(10.0, 10.0)
different = rejection_rate(2.0, 50.0)
print("same censoring:      rejection rate at 5%%: %.3f" % same)
print("different censoring: rejection rate at 5%%: %.3f" % different)

same censoring:      rejection rate at 5%: 0.065
different censoring: rejection rate at 5%: 0.035

Both rejection rates are close to the nominal 5%, as they should be when the null hypothesis is true.

Gray’s test versus a cause-specific log-rank

The difference between “how fast” and “how many” is easiest to see in an example. Give two groups the same cause-1 hazard, but a much higher cause-2 hazard in group B. Cause 1 strikes at the same rate among the survivors in both groups, but in group B most units are removed by cause 2 before cause 1 has a chance, so far fewer of them ever fail from cause 1. A cause-specific log-rank (surpyval.logrank with cause 2 treated as censored) compares the rates and finds nothing; Gray’s test compares the incidence and finds the difference:

from surpyval import logrank

x_a, e_a = simulate_cr(300, 0.1, 0.05, 20.0)    # group A
x_b, e_b = simulate_cr(300, 0.1, 0.30, 20.0)    # group B: more cause 2
x = np.concatenate([x_a, x_b])
e = np.concatenate([e_a, e_b])
group = np.repeat(["A", "B"], 300)

cs = logrank(x, group, c=np.where(e == 1, 0, 1))
gray = gray_test(x, e, group, event=1)
print("cause-specific log-rank p = %.3f" % cs.p_value)
print("Gray's test             p = %.2g" % gray.p_value)

t_grid = np.linspace(0, 40, 400)
for g, xs, es in [("A", x_a, e_a), ("B", x_b, e_b)]:
    cr = CompetingRisks.fit(xs, es)
    plt.step(t_grid, cr.cif(t_grid, 1), where="post", label="group " + g)
plt.xlabel("Time")
plt.ylabel("Cumulative incidence of cause 1")
plt.legend()

cause-specific log-rank p = 0.330
Gray's test             p = 6e-10

<matplotlib.legend.Legend at 0x76f8dd2ad220>

Neither answer is wrong. If the question is whether the groups differ in the mechanism behind cause 1, the log-rank is the relevant test; if it is whether they differ in how many units end up failing from cause 1, it is Gray’s.

The rho argument weights each event time by \(\{1 - \hat{F}^0(t^-)\}^{\rho}\), where \(\hat{F}^0\) is Gray’s pooled CIF of the cause. The default rho=0 is the standard test; rho > 0 emphasises differences in early incidence:

print("rho = 1: p = %.2g" % gray_test(x, e, group, event=1, rho=1.0).p_value)

rho = 1: p = 4e-09

Fine-Gray Sub-distribution Hazards

The Fine-Gray model estimates the effect of covariates directly on the cumulative incidence function of a chosen cause, using an inverse-probability-of-censoring-weighted subdistribution risk set. e is the per-observation cause label (None for a censored row) and cause selects the cause of interest. Z is the covariate matrix, one row per observation.

Here we simulate two-cause data whose cause-1 incidence follows a Fine-Gray model with coefficients \((0.7, -0.4)\), apply right-censoring, and recover the coefficients. (This is the simulation design of Fine and Gray’s paper: the true CIF of cause 1 is \(F_1(t \mid Z) = 1 - \{1 - p(1 - e^{-t})\}^{\exp(Z\beta)}\) with \(p = 0.5\).)

from surpyval.univariate.competing_risks import FineGray

rng = np.random.default_rng(1)
N, beta, p = 800, np.array([0.7, -0.4]), 0.5
Z = rng.uniform(-1, 1, size=(N, 2))
phi = np.exp(Z @ beta)
p1 = 1 - (1 - p) ** phi                    # P(cause = 1 | Z)
is1 = rng.uniform(size=N) < p1

x = np.empty(N)
e = np.empty(N, dtype=object)
v = rng.uniform(size=N)
w = 1 - (1 - v * p1) ** (1 / phi)          # invert the cause-1 CIF
x[is1] = (-np.log(np.clip(1 - w / p, 1e-12, 1.0)))[is1]
e[is1] = 1
x[~is1] = rng.exponential(1.0, size=N)[~is1]   # cause 2 mops up the rest
e[~is1] = 2

cens = rng.exponential(3.0, size=N)        # independent right-censoring
c = (x > cens).astype(int)
x = np.minimum(x, cens)
e[c == 1] = None

model = FineGray.fit(x, Z, e, c=c, event=1)
model

Fine-Gray Subdistribution Hazard Model
======================================
Cause of interest   : 1
Coefficients (beta'Z acts on the subdistribution hazard):
   beta_0  :   0.770136  (se 0.104097, p 0.0000)
   beta_1  :  -0.416207  (se 0.096274, p 0.0000)

The IPCW correction is what lets the coefficients come back near their true values under censoring — a naive unweighted subdistribution risk set would be biased. The fitted model stores the estimates as arrays, one entry per column of Z; np.exp(model.beta) gives the sub-distribution hazard ratios:

print("beta     :", np.round(model.beta, 3))     # also model.coefficients
print("se       :", np.round(model.se, 3))
print("p-values :", model.p_values)
print("SHR      :", np.round(np.exp(model.beta), 3))
print("cov      :\n", np.round(model.cov, 4))

beta     : [ 0.77  -0.416]
se       : [0.104 0.096]
p-values : [1.37889700e-13 1.53830841e-05]
SHR      : [2.16 0.66]
cov      :
 [[ 0.0108 -0.001 ]
 [-0.001   0.0093]]

The standard errors come from the inverse Hessian of the weighted partial likelihood (the robust variance of Fine and Gray is not implemented), so treat them as approximate. Because the model targets the incidence directly, cif reads off the cumulative incidence of the cause at any covariate value (one covariate vector per call). The dashed lines are the true CIFs:

t = np.linspace(0, 3, 200)
for z1, label, colour in [(-1.0, 'Z1 = -1', 'C0'), (1.0, 'Z1 = +1', 'C1')]:
    z = np.array([z1, 0.0])
    plt.plot(t, model.cif(t, Z=z), color=colour, label=label)
    plt.plot(t, 1 - (1 - p * (1 - np.exp(-t))) ** np.exp(z @ beta),
             '--', color=colour)
plt.legend()
plt.xlabel('Time')
plt.ylabel('Cumulative incidence of cause 1')

Text(0, 0.5, 'Cumulative incidence of cause 1')

A positive \(\beta_0\) raises the cause-1 incidence, so the Z1 = +1 curve sits above Z1 = -1. The fitted CIF is a step function built on the observed cause-1 event times, so it is flat after the last of them (about \(t = 5.6\) in this sample) rather than extrapolating. sf(x, Z) returns 1 - cif(x, Z) and phi(Z) the multiplier \(e^{Z\beta}\). As the Cox fit does, the Fine-Gray fit centres the covariates on their means, which leaves the coefficients and every prediction unchanged but keeps \(e^{Z\beta}\) from overflowing on a covariate far from zero; the baseline is then reported at \(Z = 0\), or, with center=True, at the means (model.center, and phi is relative to them).

Things to watch:

• cause must be given when the data contain more than one cause, and it must be one that is observed; with a single cause it may be omitted.

• The Fine-Gray model is fitted for one cause at a time. To model every cause, use CompetingRisksProportionalHazards with model="Fine-Gray" (below); the separate fits are not constrained to be mutually consistent.

• Censoring times tied with event times (common when times are recorded in whole days or months) follow R’s cmprsk::crr: the event is taken to come first. A unit that failed from a competing cause at \(x_i\) keeps the weight \(\hat{G}(t^-)/\hat{G}(x_i^-)\) at a later event time \(t\), with \(\hat{G}\) the Kaplan-Meier estimate of the censoring survival taken just before each time, so censorings at \(t\) or at \(x_i\) do not count against the events there. These are crr’s weights; without such ties the left limits are just \(\hat{G}(t)\) and \(\hat{G}(x_i)\). See Competing Risks Analysis for the details.

The fitted model serialises like the other competing-risks models; the optimiser result (model.res) is not stored:

reloaded = surv.from_dict(model.to_dict())
np.allclose(reloaded.cif([0.5, 1.0], Z=[1.0, 0.0]),
            model.cif([0.5, 1.0], Z=[1.0, 0.0]))

True

Cause-Specific Proportional Hazards

CompetingRisksProportionalHazards fits a proportional-hazards model per cause. With model="Cox" (the default) each cause is a Cox model with the other causes treated as censored; model="Fine-Gray" fits the subdistribution model above for every cause. It reuses the simulated data from the previous section:

from surpyval.univariate.competing_risks import (
    CompetingRisksProportionalHazards,
)

csph = CompetingRisksProportionalHazards.fit(x, Z, e, c=c, model="Cox")
# cumulative incidence of cause 1 at a covariate vector
csph.cif(np.array([0.5, 1.0, 2.0]), Z=[0.5, -0.5], event=1)

array([0.32677816, 0.4955077 , 0.65502479])

The cause-specific model answers “what drives the rate of this cause among those still at risk”, while Fine-Gray answers “what drives the eventual incidence of this cause”; the two coincide only when the competing causes are unaffected by the covariates.

Coefficients and predictions

The fitted coefficients are in betas, one row per cause. Look the row up through event_idx_map rather than assuming an order:

for cause, row in csph.event_idx_map.items():
    print("cause", cause, "cause-specific log hazard ratios:",
          np.round(csph.betas[row], 3))

cause 1 cause-specific log hazard ratios: [ 0.705 -0.323]
cause 2 cause-specific log hazard ratios: [-0.397  0.379]

In this simulation \(Z_1\) raises the cause-1 incidence (the Fine-Gray coefficient is +0.7), and it does so by raising the cause-1 hazard and lowering the cause-2 hazard: the two sets of cause-specific coefficients together produce the incidence effect.

The causes are sorted, so the row order of betas is reproducible. phi_e(Z, row) is a cause’s hazard multiplier \(e^{Z\hat\beta_k}\) (relative to a unit at the covariate means, center, for a fit with center=True, which keeps the baselines there), and results holds each cause’s optimiser result. The model also has beta and phi. These are kept for backward compatibility: beta is the sum of the rows of betas, which is not a quantity of the model, and no prediction uses it. Read the coefficients from betas.

row = csph.event_idx_map[1]
print("cause-1 hazard multiplier at z = [0.5, -0.5]:",
      np.round(csph.phi_e(np.array([0.5, -0.5]), row), 3))

cause-1 hazard multiplier at z = [0.5, -0.5]: 1.672

tie_method chooses how the Cox fits handle tied failure times (see Cox Proportional Hazards): "efron" (the default, as for CoxPH), "breslow", "exact" or "kalbfleisch-prentice" ("kp"). The simulated times are continuous, so there are no ties and every method gives the same fit. Rounding the times up to the next 0.25 creates heavy ties (the largest is 107 cause-1 failures at a single time):

x_tied = np.ceil(x * 4) / 4
tied = {}
for tm in ["breslow", "efron", "exact", "kp"]:
    fit_t = CompetingRisksProportionalHazards.fit(x_tied, Z, e, c=c,
                                                  tie_method=tm)
    tied[tm] = fit_t.betas[row]
    print("%-8s cause 1: %s" % (tm, np.round(tied[tm], 3)))

breslow  cause 1: [ 0.674 -0.306]
efron    cause 1: [ 0.72  -0.335]

exact    cause 1: [ 0.722 -0.336]

kp       cause 1: [ 0.77  -0.353]

"exact" averages the partial likelihood over every order in which the tied failures could have happened, which is the right treatment when the ties come from rounding a continuous time, as here. Efron’s approximation is very close to it; Breslow’s pulls the coefficients towards zero. "kp" fits a different model, one in which time is genuinely discrete, so its coefficients are log odds ratios rather than log hazard ratios and come out larger. Every method is fast, even with ties this heavy.

With model="Cox" the model has the usual functions, each taking the times, one covariate vector Z and an optional event:

• cif(x, Z, event) – the cumulative incidence of event at Z: over each step \(x_j\) the unit, event-free with probability \(\hat{S}(x_{j-1} \mid Z) = e^{-H}\), fails with probability \(1 - e^{-\Delta H}\), and cause \(k\) takes the share \(\Delta H_k / \Delta H\) of it (the matrix-exponential form R’s survival uses for a multi-state coxph);

• Hf/hf – the cause-specific cumulative hazard and its increment at the most recent event time; with event=None they are summed over causes, each cause with its own coefficients;

• sf/ff/df – \(e^{-H}\), \(1 - e^{-H}\) and hf * sf of those hazards. sf(x, Z) with no event is the all-cause survival; with an event it is the net (cause-alone) survival, not a probability of the real-world outcome.

times = np.array([0.5, 1.0, 2.0])
z = [0.5, -0.5]
cif1 = csph.cif(times, Z=z, event=1)
cif2 = csph.cif(times, Z=z, event=2)
print("CIF cause 1          :", np.round(cif1, 4))
print("CIF cause 2          :", np.round(cif2, 4))
print("sum of CIFs          :", np.round(cif1 + cif2, 4))
print("1 - all-cause sf     :", np.round(1 - csph.sf(times, Z=z), 4))

CIF cause 1          : [0.3268 0.4955 0.655 ]
CIF cause 2          : [0.1151 0.182  0.2434]
sum of CIFs          : [0.4418 0.6775 0.8984]
1 - all-cause sf     : [0.4418 0.6775 0.8984]

The CIFs add up to the all-cause failure probability 1 - sf, exactly, so their total never exceeds one.

With model="Fine-Gray", cif, sf (1 - cif), ff and Hf need an event and come from each cause’s Fine-Gray model; hf and df raise a ValueError because the step baseline has no pointwise density. Comparing the two fits on the same data:

fg_all = CompetingRisksProportionalHazards.fit(x, Z, e, c=c, model="Fine-Gray")
for cause, row in fg_all.event_idx_map.items():
    print("cause", cause, "Fine-Gray coefficients:",
          np.round(fg_all.betas[row], 3))
print("Fine-Gray CIF of cause 1:",
      np.round(fg_all.cif(times, Z=z, event=1), 4))

cause 1 Fine-Gray coefficients: [ 0.77  -0.416]
cause 2 Fine-Gray coefficients: [-0.64   0.462]
Fine-Gray CIF of cause 1: [0.3347 0.4975 0.649 ]

The Fine-Gray coefficients for cause 1 match FineGray.fit above. For cause 2 they have the opposite sign to cause 1, even though the covariates were only built into the cause-1 incidence: whatever raises the incidence of one cause necessarily lowers the incidence of the other.

Both kinds of fit can be saved and restored like the other competing-risks models. The dictionary holds the per-cause coefficients and baselines (and, for model="Fine-Gray", each cause’s Fine-Gray model); the optimiser results in results are not stored, so a restored model has results = None:

for fitted in [csph, fg_all]:
    reloaded = surv.from_dict(json.loads(json.dumps(fitted.to_dict())))
    print(fitted.model, type(reloaded).__name__,
          np.array_equal(reloaded.cif(times, Z=z, event=1),
                         fitted.cif(times, Z=z, event=1)))

Cox CompetingRisksProportionalHazards True
Fine-Gray CompetingRisksProportionalHazards True

Fitting from a DataFrame

fit_from_df takes the time and cause column names, and the covariates either as Z_cols (a column name or list of names) or as a formula; c_col, n_col, model and tie_method (default "efron", passed to each cause’s Cox fit) are optional. A blank/NaN cause marks a censored row, and rows with a missing (or infinite) covariate are dropped, with a warning giving the count – as are rows of Z containing NaN in fit, for both the Cox and Fine-Gray fits and for FineGray.fit. The fitted model predicts from a DataFrame of the covariate columns, read by name (their order and any other columns do not matter), or still from an array Z in the fitted column order:

frame = pd.DataFrame({"time": x, "cause": e, "z1": Z[:, 0], "z2": Z[:, 1]})
csph_df = CompetingRisksProportionalHazards.fit_from_df(
    frame, x_col="time", e_col="cause", Z_cols=["z1", "z2"]
)
print(csph_df.feature_names)
new = pd.DataFrame({"z2": [z[1]], "z1": [z[0]]})
print(np.allclose(csph_df.cif(times, Z=new, event=1), cif1),
      np.allclose(csph_df.cif(times, Z=z, event=1), cif1))

['z1', 'z2']
True True

With a formula the prediction DataFrame holds the raw covariates and the fitted formula expands them – categorical levels coded against the fitted reference level, transforms with their fitted statistics – exactly as CoxPH does, before and after to_dict / from_dict. A categorical level the model was not fitted with raises a ValueError naming the column and the level (there is no coefficient for it), and a row with a missing covariate – or a missing (NaN) time – predicts nan, leaving the other rows as they are:

import textwrap

frame["band"] = np.where(Z[:, 1] > 0, "high", "low")
csph_f = CompetingRisksProportionalHazards.fit_from_df(
    frame, x_col="time", e_col="cause", formula="z1 + band"
)
print(csph_f.feature_names)
rows = pd.DataFrame({"z1": [0.5, 0.5], "band": ["low", "high"]})
print(csph_f.cif([2.0, 2.0], rows, event=1).round(4))
try:
    csph_f.cif([2.0], pd.DataFrame({"z1": [0.5], "band": ["mid"]}), 1)
except ValueError as err:
    print(textwrap.fill(str(err), 79))

['z1', 'band[T.low]']
[0.6487 0.5012]
Unknown categorical level(s): column 'band' has the level(s) ['mid'], which are
not among the levels ['high', 'low'] of the formula term 'band'. The model has
no coefficient for a level it was not fitted with; only the levels with rows in
the fitted data can be used.