Competing risks

So far every method has assumed a single kind of event. Often that is not realistic. A patient with a precancerous condition might progress to cancer, but they might also die of something else first, and death from another cause prevents the progression from ever being observed. These are competing risks: mutually exclusive event types where the occurrence of one removes the possibility of the others. Analyzing them with ordinary survival methods, one event type at a time, gives biased and even nonsensical results. This page explains why, and shows the two standard tools Greenwood provides.

flowchart LR
    A[Event-free] --> B[Progression]
    A --> C[Death]

Figure 1: A competing-risks model: from an initial state, a subject can move to exactly one of several absorbing event states.

We use the bundled mgus2 dataset, which follows patients with monoclonal gammopathy. The competing events are progression to a plasma-cell malignancy and death without progression. The dataset records each endpoint in its own pair of columns: ptime and pstat for progression, and futime and death for death. A competing-risks analysis needs one time and one cause per subject, and gw.Outcome.first_event() describes exactly that by naming the pairs. Like a formula, it holds column names rather than data, so each fit evaluates it against the frame passed as data=.

import greenwood as gw

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

# Name each endpoint's (time, event) columns, and the first observed endpoint wins
pcm_or_death = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)

# Evaluate the description against the frame to see the response it produces
y = pcm_or_death.bind(mg)
y

Surv(type=mright, n=1384, events=975, states=('pcm', 'death'))

For each subject, the endpoint with the earliest observed event becomes the cause, and its time becomes the event time. Because a subject can experience only the first of the competing events, one time per subject is all we need. Ties go to the endpoint listed first, so a progression diagnosed at the same visit as death counts as progression. A subject with no observed event is censored at the latest of the endpoint times, the last time they were known to be event-free. Here are the first several subjects, with the times in months.

# Show each subject's event time and state
y.to_frame(format="polars").head(8)

PolarsRows8Columns2
	time f64	status f64
0	30	2
1	25	2
2	46	2
3	92	2
4	8	2
5	4	2
6	151	2
7	2	2

NoteThe same response, built by hand

first_event() replaces a recipe you may have seen elsewhere, which picks ptime for subjects who progressed and futime for everyone else, then codes the cause as 1, 2, or 0. With those two columns in the frame, R’s factor() turns the codes into states, whose first label means censoring, and the formula response is the one you would write in R:

import polars as pl

# The recipe: one event time and one cause code per subject
mg_coded = mg.with_columns(
    etime=pl.when(pl.col("pstat") == 1).then(pl.col("ptime")).otherwise(pl.col("futime")),
    cause=pl.when(pl.col("pstat") == 1).then(1).otherwise(2 * pl.col("death")),
)
gw.AalenJohansen().fit(
    "Surv(etime, factor(cause, 0:2, c('censor', 'pcm', 'death')))", data=mg_coded
)

Both give identical responses on mgus2. Outcome.first_event() states the rule once and needs no extra columns. When the values are already in hand as arrays, gw.first_event() builds the same response from them directly.

first_event() also checks the data. If a subject’s event is recorded after another endpoint’s follow-up had already ended (a progression at month 80 when survival follow-up stopped at month 60, say), it warns. That usually means the time columns do not share an origin or a unit.

Cumulative incidence, not one-minus-survival

The intuitive but wrong approach is to treat one cause as the event, censor the other, run Kaplan-Meier, and report one minus the survival curve as the probability of that cause. This overestimates the probability, sometimes badly, because it treats subjects who died of the other cause as if they could still progress. The correct quantity is the cumulative incidence function, which gives the probability of experiencing a specific cause by a given time while accounting for the competing cause.

The Aalen-Johansen estimator computes it. You fit it to the multi-state response, and it returns a curve for every cause.

# Reload the dataset and describe the competing-risks response
mg = gw.load_dataset("mgus2", backend="polars")

pcm_or_death = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)

# Fit the Aalen-Johansen estimator to get a CIF for each competing cause
aj = gw.AalenJohansen().fit(pcm_or_death, data=mg)
aj

AalenJohansen (Aalen-Johansen cumulative incidence)

states: pcm, death
n = 1384

       final CIF
pcm       0.1613
death     0.8387

The fitted estimator holds a cumulative incidence curve for each cause. To inspect it we call to_frame() and take the first couple of rows within each cause, so both curves are visible side by side.

# Export to a tidy frame and show the first two rows per cause
aj.to_frame(format="polars").group_by("cause").head(2)

PolarsRows4Columns7
	cause str	time f64	n_risk f64	estimate f64	std_error f64	conf_low f64	conf_high f64
0	pcm	1	1384	0	0	0	0
1	pcm	2	1341	0.00144616432392	0.00102185289691	0	0.00344895919936
2	death	1	1384	0.0303468208092	0.00461101747295	0.0213093926302	0.0393842489883
3	death	2	1341	0.0505931213442	0.00589210798256	0.0390448019053	0.062141440783

Each cause has its own cumulative incidence estimate, standard error, and confidence interval at every event time. Because the causes are competing, their cumulative incidences plus the event-free probability sum to one at every time, which is exactly the accounting that one-minus-Kaplan-Meier gets wrong.

WarningDo not use one-minus-Kaplan-Meier for a competing cause

Censoring the competing event and applying Kaplan-Meier treats those subjects as if they remained at risk for the cause of interest, which they do not. Always use the cumulative incidence function for competing-risks probabilities.

Stratified cumulative incidence

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

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

AalenJohansen (Aalen-Johansen cumulative incidence)

states: pcm, death
strata: 2

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

import polars as pl

# Final pcm cumulative incidence by sex stratum (last event time in each group)
aj_sex.to_frame(format="polars").filter(pl.col("cause") == "pcm").group_by("strata").tail(1)

PolarsRows2Columns8
	strata str	cause str	time f64	n_risk f64	estimate f64	std_error f64	conf_low f64	conf_high f64
0	F	pcm	394	1	0.198554321073	0.0399976113701	0.12016044332	0.276948198826
1	M	pcm	424	1	0.104460229977	0.0158156196773	0.0734621850159	0.135458274937

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

Confidence intervals

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

NoteAalenJohansen’s default differs from KaplanMeier’s

AalenJohansen defaults to conf_type="plain", matching R’s survfit for competing risks. This is not the same as KaplanMeier’s default of conf_type="log" (see Kaplan-Meier). If an analysis fits both estimators, set conf_type= explicitly on each rather than assuming the defaults match.

You set the transform at construction time.

# Reload the dataset and describe the competing-risks response
mg = gw.load_dataset("mgus2", backend="polars")

pcm_or_death = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)

# Fit with the log-log confidence interval transform
aj_loglog = gw.AalenJohansen(conf_type="log-log").fit(pcm_or_death, data=mg)

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

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

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

Comparing confidence interval transforms

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

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

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

# Assemble the comparison across all three transforms
pl.concat(rows)

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

Transform summary:

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

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

Adjusting the confidence level

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

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

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

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

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

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

Testing for group differences with Gray’s test

The Aalen-Johansen estimator and stratified CIF show the cumulative incidence curves by group, but do not formally test whether the groups differ. Gray’s test (1988) is the competing-risks analogue of the log-rank test. It compares the cumulative incidence function for a specific cause across groups, properly accounting for competing events.

# Reload the dataset and describe the response
mg = gw.load_dataset("mgus2", backend="polars")

pcm_or_death = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)

# Test whether the CIF of pcm differs between sexes
gw.grays_test(pcm_or_death, group="sex", data=mg, cause="pcm")

TestResult(method="Gray's test (cause='pcm')", statistic=1.1945, df=1, p_value=0.2744)

The test returns a TestResult with the chi-square statistic, degrees of freedom, p-value, and per-group observed vs. expected event counts. You can test any cause in the competing-risks response by changing the cause= argument.

# Test the death cause instead
gw.grays_test(pcm_or_death, group="sex", cause="death", data=mg)

TestResult(method="Gray's test (cause='death')", statistic=11.6513, df=1, p_value=0.0006416)

Gray’s test uses the subdistribution risk set, the same framework behind the Fine-Gray model. Subjects who experience a competing event remain in the risk set with inverse-probability-of- censoring weights. This makes the test sensitive to differences in the cumulative incidence function itself, not just the cause-specific hazard rate. For comparing more than two groups, the test produces a multi-degree-of-freedom chi-square statistic.

NoteGray’s test vs. the log-rank test

The log-rank test applied to cause-specific hazards asks “do the cause-specific hazard rates differ?”. Gray’s test asks “do the cumulative incidence functions differ?”. These questions can give different answers because the CIF depends on all cause-specific hazards, not just the target one.

Modeling with the Fine-Gray model

The cumulative incidence function describes groups. To model how covariates affect the incidence of a specific cause, use the Fine-Gray model. It is a regression on the subdistribution hazard, constructed so that its coefficients describe the effect of covariates on the cumulative incidence itself, which is usually what you want to report.

# Reload the dataset and describe the response
mg = gw.load_dataset("mgus2", backend="polars")

pcm_or_death = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)

# Fit the Fine-Gray model for the pcm cause with age and sex as covariates
fg = gw.FineGray(cause="pcm").fit(pcm_or_death, covariates=["age", "sex"], data=mg)
fg

FineGray (Fine-Gray subdistribution hazard model, cause='pcm')

         coef  exp(coef)  se(coef)       z         p
age   -0.0173     0.9828  0.005701  -3.035  0.002408
sexM  -0.2597     0.7713    0.1856  -1.399    0.1617

n = 1384, events = 115
Standard errors: robust (clustered)

We name pcm as the target cause and fit against age and sex. To read the fitted model we tidy it into a coefficient table, asking for exponentiated estimates so the numbers are on the hazard-ratio scale rather than the log scale.

# Exponentiate to get subdistribution hazard ratios
gw.tidy(fg, exponentiate=True, format="polars")

PolarsRows2Columns7
	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
0	age	0.982848064921	0.00570095802298	-3.03470633443	0.00240769998927	0.9719271695	0.993891671139
1	sexM	0.771282230528	0.185579308546	-1.39940662934	0.161691080051	0.536102622861	1.1096313537

The exponentiated coefficients are subdistribution hazard ratios. A value above 1 means the covariate increases the cumulative incidence of the target cause, pcm here, over time. The model reports clustered robust standard errors by default, which are appropriate given the inverse-probability weighting it uses internally.

Modeling multiple competing causes

In a competing-risks study with more than one cause of interest, you can fit separate Fine-Gray models, one for each cause, to compare the covariate effects across causes.

import pandas as pd

# Fit FineGray models for both competing causes
fg_pcm = gw.FineGray(cause="pcm").fit(pcm_or_death, covariates=["age", "sex"], data=mg)
fg_death = gw.FineGray(cause="death").fit(pcm_or_death, covariates=["age", "sex"], data=mg)

# Compare the hazard ratios side-by-side
pcm_tidy = gw.tidy(fg_pcm, exponentiate=True, format="pandas")
death_tidy = gw.tidy(fg_death, exponentiate=True, format="pandas")

# Label each table by cause and stack them for comparison
pcm_tidy["cause"] = "PCM"
death_tidy["cause"] = "Death"

pd.concat([pcm_tidy, death_tidy])

PandasRows4Columns8
	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64	cause str
0	age	0.982848064921	0.00570095802298	-3.03470633443	0.00240769998927	0.9719271695	0.993891671139	PCM
1	sexM	0.771282230528	0.185579308546	-1.39940662934	0.161691080051	0.536102622861	1.1096313537	PCM
2	age	1.06035605364	0.0036793490261	15.9280217679	4.04967555089e-57	1.05273691623	1.06803033421	Death
3	sexM	1.44900400068	0.0667809170816	5.55362879898	2.79799326688e-08	1.27123419034	1.65163320019	Death

This comparison reveals how covariates affect each cause differently. A covariate might increase the incidence of progression to PCM while decreasing the incidence of death from other causes (or vice versa) which both tell important clinical stories.

NoteTwo kinds of hazard, two kinds of question

A cause-specific Cox model answers “what drives the rate of this event among those still at risk?”, while the Fine-Gray model answers “what drives the cumulative probability of this event?”. They can point in different directions for the same covariate, and both can be correct. Decide which question you are asking before choosing the model.

Cause-specific Cox regression

The Fine-Gray model targets the cumulative incidence directly, which is useful for prediction and policy questions. But sometimes the scientific question is about the rate at which an event occurs among those still at risk, not the marginal probability of the event. A cause-specific Cox model answers that by treating competing events as censoring: subjects who experience a different cause leave the risk set, and a standard Cox model is fitted on the remaining data.

CauseSpecificCox is a convenience wrapper that automates this re-coding. You pass a multi-state response and name the target cause. Internally it creates a right-censored response where only the target cause counts as an event.

# Fit a cause-specific Cox model for the pcm cause
csc_pcm = gw.CauseSpecificCox(cause="pcm").fit(pcm_or_death, covariates=["age", "sex"], data=mg)
csc_pcm

CauseSpecificCox (cause-specific Cox model, cause='pcm', ties='efron')

          coef  exp(coef)  se(coef)       z       p
age    0.01304      1.013  0.008259   1.579  0.1144
sexM  -0.02514     0.9752    0.1885  -0.133  0.8939

n = 1384, events = 115
Likelihood ratio test = 2.682 on 2 df, p = 0.2616

The fitted object delegates to a standard CoxPH, so all its features are available: hazard ratios, global tests, predictions, Schoenfeld residuals, and cox_zph() diagnostics.

# Hazard ratios for the pcm-specific model
gw.tidy(csc_pcm, exponentiate=True, format="polars")

PolarsRows2Columns7
	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
0	age	1.01312394262	0.00825868591171	1.57877051681	0.114388703517	0.996856793129	1.02965654664
1	sexM	0.975175534022	0.188455848193	-0.133388215518	0.893886364781	0.674013954399	1.41090153393

You can fit a second model for the competing cause and compare the two side by side.

# Fit a cause-specific Cox model for death
csc_death = gw.CauseSpecificCox(cause="death").fit(pcm_or_death, covariates=["age", "sex"], data=mg)
gw.tidy(csc_death, exponentiate=True, format="polars")

PolarsRows2Columns7
	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
0	age	1.06697086195	0.00362027556808	17.9057263304	1.06398034763e-71	1.05942684946	1.07456859418
1	sexM	1.4817530258	0.0696981045797	5.641844439	1.68238072112e-08	1.29255397795	1.6986462979

This is the standard analysis layout in competing-risks studies: one Fine-Gray model (for cumulative incidence and prediction) and one or two cause-specific Cox models (for etiological insight). When the two model types disagree on the direction of an effect, both are usually reported because they answer different questions.

Predicting cumulative incidence from a Fine-Gray model

Fitting a Fine-Gray model tells you which covariates matter, but the clinical question is often “what is this patient’s predicted probability of experiencing the event by time t?” The predict_cumulative_incidence() method answers that. It uses the baseline subdistribution cumulative hazard from the training data together with the fitted coefficients to compute F(t \mid x) = 1 - \exp\bigl(-H_0(t)\,\exp(x\beta)\bigr) for each subject.

# Rebuild the dataset and fit the Fine-Gray model
mg = gw.load_dataset("mgus2", backend="polars")

pcm_or_death = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)

fg = gw.FineGray(cause="pcm").fit(pcm_or_death, covariates=["age", "sex"], data=mg)

# Predict cumulative incidence at 10, 20, and 30 years for three subjects
fg.predict_cumulative_incidence(mg[:3], times=[120, 240, 360], format="polars")

PolarsRows3Columns4
	time f64	subject_1 f64	subject_2 f64	subject_3 f64
0	120	0.190249092054	0.221889387455	0.136459495566
1	240	0.285958812753	0.329968422973	0.208769010861
2	360	0.368295936626	0.420790435227	0.273375369165

Each column is one subject and each row is a time point. The values are the predicted probability of PCM progression by that time, accounting for the competing risk of death. A subject with higher risk covariates accumulates incidence faster.

You can also call predict() for the raw linear predictor or relative subdistribution hazard, which are useful for model evaluation with concordance_index_incidence() and brier_score_incidence().

# Linear predictor (centered) for the first five subjects
fg.predict(mg[:5], type="lp")

array([-0.16279108,  0.01021626, -0.52629639, -0.07647732, -0.19739255])

Penalized Fine-Gray for high-dimensional covariates

When the number of covariates is large relative to the number of target-cause events, the unpenalized Fine-Gray model may overfit or fail to converge. PenalizedFineGray adds an elastic-net penalty to the subdistribution-hazard partial likelihood, enabling variable selection and coefficient shrinkage. It is the competing-risks analogue of CoxNet.

The penalty is \lambda\bigl(\alpha\|\beta\|_1 + \tfrac{1-\alpha}{2}\|\beta\|_2^2\bigr), where penalizer is \lambda and l1_ratio is \alpha. Setting l1_ratio=1 gives lasso (sparse), l1_ratio=0 gives ridge (smooth shrinkage), and values in between blend the two.

# Reload the dataset and describe the response
mg = gw.load_dataset("mgus2", backend="polars")

pcm_or_death = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)

# Fit a lasso-penalized Fine-Gray model
pfg = gw.PenalizedFineGray(cause="pcm", penalizer=0.01, l1_ratio=1.0).fit(
    pcm_or_death, covariates=["age", "sex"], data=mg
)
pfg

PenalizedFineGray (elastic-net Fine-Gray, cause='pcm', penalizer=0.01, l1_ratio=1.0)

           coef
age   -0.008028
sexM  -0.002464

n = 1384, events = 115, nonzero coefficients = 2

The fitted object shows the penalized coefficients and how many survived shrinkage. Just like the unpenalized FineGray, you can predict cumulative incidence for new subjects.

# Predict cumulative incidence at landmark times
pfg.predict_cumulative_incidence(mg[:3], times=[120, 240, 360], format="polars")

PolarsRows3Columns4
	time f64	subject_1 f64	subject_2 f64	subject_3 f64
0	120	0.0553047921496	0.0597870862744	0.0526473972235
1	240	0.0868374562654	0.0937456817667	0.0827325611344
2	360	0.116744848103	0.125863551583	0.111314476378

To examine the penalized coefficients as a tidy frame, call to_frame(). Because penalized coefficients are biased by design, the table reports point estimates and hazard ratios but not standard errors or p-values.

# Export penalized coefficients
pfg.to_frame(format="polars")

PolarsRows2Columns3
	term str	estimate f64	hazard_ratio f64
0	age	-0.00802846325941	0.992003678777
1	sexM	-0.00246350050574	0.997539531421

Choosing the penalty strength

A stronger penalty shrinks more coefficients toward (or exactly to) zero, improving generalizability at the cost of bias. To find the right balance, fit models at several penalty values and evaluate with concordance_index_incidence() or brier_score_incidence() on held-out data.

import numpy as np

# Compare penalized models at different strengths
for lam in [0.001, 0.01, 0.1]:
    m = gw.PenalizedFineGray(cause="pcm", penalizer=lam, l1_ratio=1.0).fit(
        pcm_or_death, covariates=["age", "sex"], data=mg
    )
    n_nz = int(np.count_nonzero(m.coef_))
    print(f"  penalizer={lam:.3f}  nonzero={n_nz}  coef={np.round(m.coef_, 4)}")

  penalizer=0.001  nonzero=2  coef=[-0.0164 -0.2341]
  penalizer=0.010  nonzero=2  coef=[-0.008  -0.0025]
  penalizer=0.100  nonzero=0  coef=[-0. -0.]

TipRidge vs. lasso for competing risks

Use lasso (l1_ratio=1) when you want automatic variable selection, driving irrelevant coefficients to exactly zero. Use ridge (l1_ratio=0) when all covariates are expected to contribute but you want to stabilize estimates against multicollinearity. Elastic net (0 < l1_ratio < 1) blends both behaviors.

Competing-risk random forest

The parametric and semiparametric models above all assume log-linear covariate effects. When the relationship between covariates and cumulative incidence is non-linear or involves interactions, a tree-based model can capture these patterns automatically.

CompetingRiskForest grows an ensemble of survival trees that split nodes using the log-rank CR statistic. Unlike a standard random survival forest, the splitting criterion uses a Lau-inclusive risk set: subjects who experienced a competing event remain at risk under the subdistribution framework, so each split directly targets separation of the cumulative incidence functions. At each leaf the Aalen-Johansen CIF is estimated for every cause, and the forest prediction averages these across all trees.

# Fit a competing-risk forest
crf = gw.CompetingRiskForest(
    n_estimators=100, max_depth=5, oob_score=True, random_state=0
).fit(pcm_or_death, covariates=["age", "sex"], data=mg)
crf

CompetingRiskForest (100 logrankCR trees, max_features='sqrt')
n = 1384, events = 975, features = 2
causes: pcm, death
out-of-bag concordance = 0.4511

Predict the cumulative incidence for a specific cause.

# Predict CIF for the pcm cause at landmark times
crf.predict(mg[:5], cause="pcm", times=[365, 730, 1095], format="polars")

PolarsRows3Columns6
	time f64	subject_1 f64	subject_2 f64	subject_3 f64	subject_4 f64	subject_5 f64
0	365	0.0913585599972	0.10547895823	0.0473698968742	0.128019373571	0.0905998149289
1	730	0.101280698164	0.115401096397	0.0473698968742	0.128019373571	0.100521953096
2	1095	0.101280698164	0.115401096397	0.0473698968742	0.128019373571	0.100521953096

Because the forest uses bootstrap aggregation, you can inspect out-of-bag variable importance to understand which covariates drive the CIF.

# Permutation variable importance for the death cause
crf.variable_importance(cause="death", n_repeats=3, format="polars")

PolarsRows2Columns2
	term str	importance f64
0	age	0.144957144767
1	sexM	0.0138373553059

The predict_risk() method returns a scalar summary (summed CIF over the event-time grid) that can be used for ranking or as input to concordance_index_incidence().

# Scalar risk scores
risk = crf.predict_risk(mg[:5], cause="pcm")
risk

array([ 9.96696612, 11.7327757 ,  5.99817941, 12.83037389, 10.15096961])

Next steps

You can now estimate and model competing-risks data correctly.

• Plot the results with plot_cif() from Visualizing survival. It accepts a fitted AalenJohansen directly and facets by cause, with optional grouping by covariate and an aligned numbers-at-risk table.
• Multi-state models generalize competing risks to settings where subjects move between several states over time.
• Evaluate CIF predictions with brier_score_incidence(), concordance_index_incidence(), calibration_incidence(), and accuracy_in_time() from Prediction performance, which provide discrimination, calibration, and accuracy metrics for cause-specific cumulative incidence.
• Machine-learning predictions with SurvivalBoost from Machine learning. Where Fine-Gray assumes log-linear effects and models one cause at a time, SurvivalBoost estimates all CIFs jointly using gradient-boosted trees, capturing non-linear effects and interactions automatically.
• Generate synthetic data with simulate_competing_risks() for testing, benchmarking, or teaching. It produces latent-Weibull competing risks with optional covariates and censoring.
• Revisit Survival data for the details of the multi-state response.