Competing Risks

Models for units that can fail from one of several distinct causes, where the occurrence of one cause removes the unit from risk of the others. Import them from surpyval.univariate.competing_risks.

For a narrative introduction with worked examples, see Competing Risks SurPyval Modelling; for the statistical background, see Competing Risks Analysis. The Gray test for comparing cumulative incidence between groups is documented with the other hypothesis tests in Comparison Tests and Validation Metrics.

Non-Parametric (Aalen-Johansen)

The cumulative incidence of each cause estimated without a model, and the incidence-increment helper it (and Gray’s test) is built on.

class surpyval.univariate.competing_risks.nonparametric.competing_risks.CompetingRisks

Bases: SerialisableMixin

Non-parametric competing-risks estimate: the Aalen-Johansen cumulative incidence function (CIF) of each cause, and the all-cause and cause-specific (net) hazards and survival.

Each unit fails from one of several causes e (or is right-censored, with no cause). The cumulative incidence of cause \(j\) is the probability of failing from that cause by time \(t\) while the others still act,

\[F_j(t) = \sum_{t_i \le t} S(t_{i-1}) \frac{d_{ij}}{r_i},\]

with \(S\) the all-cause Kaplan-Meier survival, \(d_{ij}\) the failures from cause \(j\) at \(t_i\) and \(r_i\) the number at risk. The CIFs of all causes add up to the all-cause failure probability.

Call the class method CompetingRisks.fit (or fit_from_df); it returns a fitted instance.

Hf(x: ArrayLike, event: Any = None) → NDArray

Cumulative hazard, all causes (event=None) or one cause. With the Nelson-Aalen method it is the sum of the hazard increments d / r; with Kaplan-Meier it is -log of the product-limit survival, so that sf == exp(-Hf) either way.

cif(x: ArrayLike, event: Any) → NDArray

Cumulative incidence of cause event at x: the probability of having failed from that cause by x, with the other causes still acting. event is required.

df(x: ArrayLike, event: Any = None) → NDArray

hf * sf: the probability mass at each event time, all causes (event=None) or from one cause’s net survival.

ff(x: ArrayLike, event: Any = None) → NDArray

1 - sf: all causes, or the net failure probability from one cause with the others removed (treated as censoring). For the probability of failing from a cause while the others still act, use cif().

classmethod fit(x: ArrayLike, e: ArrayLike, c: ArrayLike | None = None, n: ArrayLike | None = None, how: str = 'Nelson-Aalen') → CompetingRisks

Fit the non-parametric competing-risks model.

Parameters:

  • x (array_like) – Failure or censoring times.

  • e (array_like) – The cause of each failure: any hashable labels (integers, strings, tuples, or a mix); they are sorted to fix their order in event_idx_map (labels of different types by type name, then text). A missing value (None, NaN) marks a right-censored row.

  • c (array_like, optional) – Censoring flags: 0 a failure (with a cause in e), 1 right-censored (with e missing). Derived from e if not given. Left and interval censoring are not supported.

  • n (array_like, optional) – Counts. Defaults to 1.

  • how (str, optional) – The all-cause survival estimator that sf, ff and Hf report: "Nelson-Aalen" (the default, exp(-H)) or "Kaplan-Meier". The cumulative incidence always uses the Kaplan-Meier survival, so the CIFs add up to the Kaplan-Meier failure probability either way.

Returns:

The fitted model.

Return type:

CompetingRisks

Examples

>>> from surpyval.univariate.competing_risks import CompetingRisks
>>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> e = ['a', 'b', 'a', None, 'a', 'b', 'a', None, 'b', 'a']
>>> model = CompetingRisks.fit(x, e)
>>> model.cif([5, 10], 'a').round(4)
array([0.3167, 0.6083])
>>> model.cif([5, 10], 'b').round(4)
array([0.1   , 0.3917])
classmethod fit_from_df(df: Any, x_col: str, e_col: str, c_col: str | None = None, n_col: str | None = None, how: str = 'Nelson-Aalen') → CompetingRisks

Fit from the columns of a pandas.DataFrame.

Parameters:

  • df (DataFrame) – The data.

  • x_col (str) – The column of failure / censoring times.

  • e_col (str) – The column of causes (missing for a censored row).

  • c_col (str, optional) – The column of censoring flags; derived from e_col if not given.

  • n_col (str, optional) – The column of counts.

  • how (str, optional) – As for fit().

Returns:

The fitted model; the frame is kept as source_df.

Return type:

CompetingRisks

classmethod from_dict(model_dict: dict) → CompetingRisks

Rebuild a nonparametric competing-risks model from a dict.

classmethod from_json(fp: str | PathLike) → Any

Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

hf(x: ArrayLike, event: Any = None) → NDArray

Hazard (the Nelson-Aalen increment d / r at each event time, 0 between them), all causes (event=None) or one cause.

how: str = 'Nelson-Aalen'

The survival estimator sf/ff/Hf report: "Nelson-Aalen" (exp(-H)) or "Kaplan-Meier" (product limit).

iif(x: ArrayLike, event: Any) → NDArray

Instantaneous incidence of cause event: the step of the cumulative incidence at each event time, \(S(t_{i-1}) d_{ij} / r_i\) (0 between event times).

plot(stacked: bool = True, ax: Any = None) → Any

Plot the cumulative incidence of every cause (#485).

Parameters:

  • stacked (bool, optional) – Stack the causes’ cumulative incidences (the default), so the top of the stack is the all-cause failure probability \(1 - S\); False draws each as its own step curve.

  • ax (matplotlib.axes.Axes, optional) – The axes to draw on; the current axes by default.

Returns:

The axes drawn on.

Return type:

matplotlib.axes.Axes

Examples

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> import matplotlib.pyplot as plt
>>> from surpyval.univariate.competing_risks import CompetingRisks
>>> x = [1, 2, 3, 4, 5, 6, 7, 8]
>>> e = ["a", "b", "a", "b", "a", None, "a", "b"]
>>> c = [0, 0, 0, 0, 0, 1, 0, 0]
>>> model = CompetingRisks.fit(x, e, c=c)
>>> fig, ax = plt.subplots()
>>> model.plot(ax=ax).get_ylabel()
'Cumulative incidence'
>>> plt.close(fig)
set_support(lower: float, upper: float) → CompetingRisks

Give the estimate an explicit support, [lower, upper].

Without one, every function starts at its initial value before the first time and holds its last value after the last, however far from the data. With a support set, they do so only within them:

  • in [lower, x[0]): sf 1, and ff, Hf, hf, df, iif and cif 0;

  • in (x[-1], upper]: the value at the last time, carried (for hf, df and iif, as without bounds, that of the step containing the last time);

  • outside [lower, upper]: NaN.

x[0] and x[-1] are the first and last observed times, failures or censorings. The variable need not be time: lower may be negative, and either bound infinite. The bounds are kept by to_dict.

Parameters:

  • lower (float) – The lower end of the support; at most the first time, x[0].

  • upper (float) – The upper end; at least the last time, x[-1], and above lower.

Returns:

The model itself, so the call can be chained.

Return type:

CompetingRisks

Raises:

ValueError – If a bound is NaN or not a number, lower is not below upper, or the bounds do not contain the observed times.

Examples

>>> from surpyval.univariate.competing_risks import CompetingRisks
>>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> e = ['a', 'b', 'a', None, 'a', 'b', 'a', None, 'b', 'a']
>>> model = CompetingRisks.fit(x, e).set_support(0, 20)
>>> model.cif([-1, 0.5, 5, 15, 25], 'a').round(4)
array([   nan, 0.    , 0.3167, 0.6083,    nan])
sf(x: ArrayLike, event: Any = None) → NDArray

Survival, all causes (event=None) or the net survival from one cause, by the estimator the model was fitted with: exp(-H) (Nelson-Aalen, the default) or the product limit (Kaplan-Meier). The one-cause survival treats the other causes as censoring, so it is not the probability of escaping that cause in the presence of the others – use cif() for that.

source_df: Any

The source DataFrame when fitted via fit_from_df (named so it does not shadow the density method df, #253).

support: tuple[float, float] | None = None

The (lower, upper) interval the estimate is defined on, set by set_support(); None (the default) when it has not been set.

to_dict() → dict

Serialise this fitted nonparametric competing-risks model to a plain, JSON-serialisable dict: the event index map and the fitted step arrays (the shared baseline plus the per-event incidence and cumulative- incidence functions). The reloaded model reproduces every prediction exactly.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

  • fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.

  • with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

surpyval.univariate.competing_risks.aalen_johansen.aalen_johansen_iif(S: NDArray, hazard_increments: NDArray) → NDArray

Instantaneous incidence increments.

The cause-specific hazard increment at t_i acts on the population still alive just before t_i, so each increment is weighted by S(t_i-) — the survival after the previous event time — not S(t_i) (#253). S and the last axis of hazard_increments must be aligned on the same event-time grid; the cumulative incidence is the cumulative sum of the result.

Parameters:

  • S (array_like) – The all-cause (Kaplan-Meier) survival at each event time.

  • hazard_increments (array_like) – The cause-specific hazard increments d_j / r at the same times, one row per cause (or a single 1-D row).

Returns:

hazard_increments * S(t-), the same shape as hazard_increments.

Return type:

numpy array

Examples

>>> import numpy as np
>>> from surpyval.univariate.competing_risks.aalen_johansen import (
...     aalen_johansen_iif,
... )
>>> S = np.array([0.9, 0.8, 0.6])
>>> aalen_johansen_iif(S, np.array([0.1, 0.0, 0.25]))
array([0.1, 0. , 0.2])

Parametric

One distribution per cause, combined into cumulative incidence functions.

class surpyval.univariate.competing_risks.parametric.parametric_competing_risks.ParametricCompetingRisks

Bases: SerialisableMixin

A parametric competing-risks model: one distribution per cause, combined into cumulative-incidence functions. Build it in one step from data with fit() / fit_from_df(), or assemble it from already-fitted per-cause models – each of any distribution family – with from_fitted().

Hf(x: ArrayLike, event: Any = None) → NDArray

Cumulative hazard. event=None gives the all-cause cumulative hazard \(\sum_k H_k\); event=k gives cause k’s.

aic() → float

Akaike information criterion of the joint model, 2 K + 2 neg_ll with K the number of parameters estimated over all causes: the sum of the causes’ AICs, since both terms add over causes.

bic() → float

Bayesian information criterion of the joint model.

2 neg_ll + K ln(n), with K the number of parameters estimated over all causes and n the sample size every SurPyval BIC uses, counted on the whole data: the observed failures of any cause, weighted by their counts (right-censored units add nothing). That is the sum of the causes’ own sample sizes, since each counts the failures of its cause; for a model assembled with from_fitted() it is that sum, which is the whole data’s count when each model was fitted to the cause-specific view of the same data.

It is not the sum of the causes’ BICs, which charged each cause’s parameters ln of its own cause’s failures (and was what this method returned before v0.21), a smaller penalty.

Examples

>>> import numpy as np
>>> from surpyval import Exponential
>>> from surpyval.univariate.competing_risks import (
...     ParametricCompetingRisks,
... )
>>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> e = ['a', 'b', 'a', None, 'a', 'b', 'a', None, 'b', 'a']
>>> model = ParametricCompetingRisks.fit(x, e, dist=Exponential)
>>> round(model.bic(), 4)
61.5902
>>> round(float(2 * model.neg_ll() + 2 * np.log(8)), 4)
61.5902
cif(x: ArrayLike, event: Any = None) → NDArray | float

Cumulative incidence function. event=k returns \(\int_0^t f_k^{\mathrm{sub}}(u)\,du\), the probability of having failed from cause k by t; event=None gives the all-cause incidence \(1 - S(t) = \sum_k \mathrm{CIF}_k(t)\).

The integral is taken over the cause’s own probability scale,

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

(the substitution \(p = F_k(u)\)), by adaptive quadrature to a relative accuracy of about \(10^{-10}\). The integrand is bounded by 1 and monotone, so the result is accurate for every requested time independently of the others, over any span of times, and also where a cause’s density is infinite (a Weibull shape below 1) or very heavy-tailed. The causes’ CIFs sum to the all-cause ff to that accuracy. Each cause’s model needs a quantile function qf.

Parameters:

  • x (array_like or float) – Times at which to evaluate the incidence; inf gives the eventual probability of the cause (see probability_of_cause()).

  • event (optional) – The cause; None for all causes combined.

Returns:

The cumulative incidence at each time, in the shape of x (a float for a scalar x).

Return type:

numpy array or float

ff(x: ArrayLike) → NDArray

All-cause failure probability \(1 - S(t)\) (the total cumulative incidence over all causes).

classmethod fit(x: ArrayLike, e: ArrayLike, c: ArrayLike | None = None, n: ArrayLike | None = None, dist: ~typing.Any = <surpyval.univariate.parametric.distributions.weibull.Weibull_ object>, how: str = 'MLE') → ParametricCompetingRisks

Fit a parametric distribution to each cause’s cause-specific hazard.

Parameters:

  • x (array_like) – Observed times.

  • e (array_like) – The cause of each observation: any hashable labels (integers, strings, tuples, or a mix), kept in sorted order in causes. A missing value (None / NaN) marks a censored observation with no attributed cause.

  • c (array_like, optional) – Censoring flag (0 observed, 1 right-censored). If omitted it is derived from e – a missing event is censored, an event present is observed – so data can be passed as (x, e) alone. Left/interval censoring is not supported.

  • n (array_like, optional) – Counts per observation.

  • dist (ParametricFitter or dict, optional) – The distribution fitted to each cause (default Weibull). Pass a {cause: distribution} mapping, with an entry for every cause, to use a different distribution per cause.

  • how (str, optional) – Estimation method passed to each distribution’s fit (default "MLE").

Returns:

The fitted model.

Return type:

ParametricCompetingRisks

Examples

>>> from surpyval import Exponential
>>> from surpyval.univariate.competing_risks import (
...     ParametricCompetingRisks,
... )
>>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> e = ['a', 'b', 'a', None, 'a', 'b', 'a', None, 'b', 'a']
>>> model = ParametricCompetingRisks.fit(x, e, dist=Exponential)
>>> model.cif([5, 10], 'a').round(4)
array([0.323 , 0.4791])
>>> round(model.probability_of_cause('a'), 4)
0.625
classmethod fit_from_df(df: ~typing.Any, x_col: str, e_col: str, c_col: str | None = None, n_col: str | None = None, dist: ~typing.Any = <surpyval.univariate.parametric.distributions.weibull.Weibull_ object>, how: str = 'MLE') → ParametricCompetingRisks

Fit from the columns of a pandas.DataFrame; see fit().

Parameters:

  • df (DataFrame) – The data.

  • x_col (str) – The time and cause columns.

  • e_col (str) – The time and cause columns.

  • c_col (str, optional) – The censoring-flag and count columns.

  • n_col (str, optional) – The censoring-flag and count columns.

  • dist (optional) – As for fit().

  • how (optional) – As for fit().

Returns:

The fitted model.

Return type:

ParametricCompetingRisks

classmethod from_dict(model_dict: dict) → ParametricCompetingRisks

Rebuild a parametric competing-risks model from a dict.

classmethod from_fitted(models: dict) → ParametricCompetingRisks

Assemble a competing-risks model from already-fitted single-cause models – one per cause – instead of fitting them here.

Each cause’s model may be of a completely different family: a Weibull with a limited-failure (cure) fraction for one cause, a LogNormal for another, a discrete distribution for a third, and so on. The only requirement is that every model exposes the standard surpyval model interface (sf / ff / df / hf / Hf, and the quantile function qf, over which the cumulative incidence is integrated and the samples are drawn); the cumulative incidence, all-cause survival and sampling are then assembled from them exactly as for a fit() model.

This is the right entry point when each cause has been modelled separately – for example fitted with its own distribution, offset, limited-failure or zero-inflated options – and you want to combine them into one competing-risks object.

Parameters:

models (dict or sequence) – Either a {cause: model} mapping, or a sequence of models whose causes are taken to be their positions 0, 1, 2, ....

Returns:

The assembled model.

Return type:

ParametricCompetingRisks

Notes

Each per-cause model should be fitted to the cause-specific view of the data (that cause’s events observed, every other cause’s events and every censored unit treated as right-censored) for the assembled CIFs to be the competing-risks quantities. If the causes carry a cure fraction the all-cause survival need not fall to zero, so the cause probabilities need not sum to one – some units never fail.

classmethod from_json(fp: str | PathLike) → Any

Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

hf(x: ArrayLike, event: Any = None) → NDArray

Hazard rate. event=None is the all-cause hazard \(\sum_k h_k\); event=k is the cause-specific hazard.

iif(x: ArrayLike, event: Any) → NDArray

Instantaneous incidence function (the subdistribution density) of a cause: \(f_k^{\mathrm{sub}}(t) = h_k(t) S(t) = f_k(t) \prod_{j\neq k} S_j(t)\).

property maximum: str

What the causes’ fits reached, one of MAXIMUM_STATES (surpyval.utils.no_maximum): the worst of the causes’ own maximum, since the joint likelihood is maximised a cause at a time (principles 12 and 13).

Examples

>>> from surpyval import Exponential
>>> from surpyval.univariate.competing_risks import (
...     ParametricCompetingRisks,
... )
>>> x = [1, 2, 3, 4, 5, 6, 7, 8]
>>> e = ["a", "b", "a", "b", "a", "b", "a", "b"]
>>> ParametricCompetingRisks.fit(x, e, dist=Exponential).maximum
'verified'
neg_ll() → float

Total negative log-likelihood: the sum over the per-cause fits.

probability_of_cause(event: Any) → Any

The eventual probability that a unit fails from event, \(\mathrm{CIF}_k(\infty)\). These sum to one over all causes unless a cause has a cure (limited-failure) fraction, in which case they sum to the all-cause probability of ever failing.

It is cif() at inf: the integral runs over the whole of the cause’s probability scale, so no finite horizon is chosen and a very heavy-tailed cause (a LogNormal with a large \(\sigma\)) is as accurate as any other.

random(size: int, random_state: int | None = None) → NDArray

Draw size samples of (time, cause) from the model, using the latent-failure-time representation: draw a latent time from each cause and take the earliest, recording its cause.

Each latent time is drawn by inverse-transform sampling through the cause’s quantile function, so it works for any per-cause model, including limited-failure (cure) models – there the quantile is infinite above the cure ceiling, so a draw in the cure region yields an infinite latent time. A unit whose every latent time is infinite never fails; it is returned with x = inf and cause None.

Returns a structured array with fields x and e.

sf(x: ArrayLike) → NDArray

All-cause survival \(S(t) = \prod_k S_k(t)\).

to_dict() → dict

Serialise this fitted parametric competing-risks model to a plain, JSON-serialisable dict: the list of causes and each cause’s fitted distribution (via its own to_dict). The reloaded model reproduces every cumulative-incidence / hazard function exactly.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

  • fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.

  • with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

Regression

Covariate models for competing risks: the Fine-Gray subdistribution hazards model (FineGray is an instance of FineGray_ below; its fit returns a FineGrayModel), and CompetingRisksProportionalHazards, which fits either a cause-specific Cox model per cause or a Fine-Gray model per cause.

class surpyval.univariate.competing_risks.regression.fine_gray.FineGray_

Bases: object

The Fine-Gray subdistribution-hazards regression for one cause of a competing-risks problem: the covariates act proportionally on the subdistribution hazard of the cause of interest, so a coefficient describes its effect on that cause’s cumulative incidence directly,

\[F_k(t \mid Z) = 1 - \exp\left(-\Lambda_{k0}(t)\, e^{\beta' Z}\right).\]

Estimated by inverse-probability-of-censoring weighting (IPCW), with one Kaplan-Meier censoring distribution for the whole sample, so censoring is assumed not to depend on the covariates. A subject 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\): the censoring survival is taken just before each time, so a censoring tied with an event counts after it, as in R’s cmprsk::crr. FineGray (from surpyval.univariate.competing_risks) is an instance of this class; its fit returns a FineGrayModel.

fit(x: ArrayLike, Z: ArrayLike, e: ArrayLike, c: ArrayLike | None = None, n: ArrayLike | None = None, event: Any = None, center: bool = False) → FineGrayModel

Fit the Fine-Gray model for a cause of interest.

Parameters:

  • x (array_like) – Observed times.

  • Z (ndarray) – Covariate matrix, one row per observation. Rows with a missing (NaN) or infinite covariate are dropped, with a warning.

  • e (array_like) – Event-type (cause) labels; None for a censored observation.

  • c (array_like, optional) – Censoring flags (0 observed, 1 right-censored). Defaults to deriving them from e: a missing event (None/NaN) is right-censored, any other is observed. Left/interval censoring is not supported.

  • n (array_like, optional) – Counts per observation. Defaults to 1.

  • event (optional) – The cause of interest (a label in e). May be omitted only when the data contains a single event type. The fitted model keeps it as cause.

  • center (bool, optional) – False (the default) reports the baseline cumulative subdistribution hazard at Z = 0; True reports it at the covariate means, stored as model.center, and phi is then relative to them. The fit runs on centred covariates either way, so the coefficients and every prediction are the same; the default refuses, with a ValueError, covariates so far from 0 that the baseline there over- or underflows.

Returns:

The fitted model, with cif() prediction.

Return type:

FineGrayModel

Examples

>>> from surpyval.univariate.competing_risks import FineGray
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> Z = rng.binomial(1, 0.5, (200, 1)).astype(float)
>>> t_a = rng.exponential(1 / (0.1 * np.exp(0.7 * Z[:, 0])))
>>> t_b = rng.exponential(1 / 0.05, 200)
>>> t_c = rng.uniform(0, 20, 200)  # censoring times
>>> x = np.minimum(np.minimum(t_a, t_b), t_c).round(3)
>>> first = np.where(t_a < t_b, "a", "b")
>>> e = np.where(t_c < np.minimum(t_a, t_b), None, first)
>>> model = FineGray.fit(x, Z, e, event="a")
>>> model.beta.round(3)
array([0.908])
>>> model.cif([5, 10], [[1]]).round(4)
array([0.5808, 0.7395])
fit_from_df(df: Any, x_col: str, e_col: str, Z_cols: str | list[str], c_col: str | None = None, n_col: str | None = None, **fit_options: Any) → FineGrayModel

Fit the Fine-Gray model from the columns of a pandas.DataFrame.

The column names are passed in place of the arrays fit() takes, with the names of every competing-risks fit_from_df (CompetingRisksProportionalHazards.fit_from_df too); event and center are passed to fit() unchanged.

Parameters:

  • df (pandas.DataFrame) – The data.

  • x_col (str) – Column of observed times.

  • e_col (str) – Column of event-type (cause) labels. Use None (or a blank/NaN cell) for a censored observation.

  • Z_cols (str or list of str) – Covariate column(s), in the order of beta.

  • c_col (str, optional) – Column of censoring flags (0 observed, 1 right-censored).

  • n_col (str, optional) – Column of counts per row.

  • **fit_options – event (the cause of interest) and center, as for fit().

Returns:

The model fit() returns for the same arrays.

Return type:

FineGrayModel

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from surpyval.univariate.competing_risks import FineGray
>>> rng = np.random.default_rng(0)
>>> z = rng.binomial(1, 0.5, 200).astype(float)
>>> t_a = rng.exponential(1 / (0.1 * np.exp(0.7 * z)))
>>> t_b = rng.exponential(1 / 0.05, 200)
>>> df = pd.DataFrame({
...     "time": np.minimum(t_a, t_b).round(3),
...     "cause": np.where(t_a < t_b, "a", "b"),
...     "z": z,
... })
>>> model = FineGray.fit_from_df(
...     df, x_col="time", e_col="cause", Z_cols="z", event="a"
... )
>>> model.beta.round(3)
array([0.663])
class surpyval.univariate.competing_risks.regression.fine_gray.FineGrayModel(fit: dict)

Bases: LinearPredictorMixin, SerialisableMixin

A fitted Fine-Gray subdistribution-hazard model for one cause of interest.

The natural prediction is the cumulative incidence function cif(); coefficients/se/p_values describe the (log) subdistribution hazard ratios.

property aliased: NDArray

The columns of Z whose coefficients the data cannot determine (#476): a constant column, which the baseline subdistribution hazard absorbs, or a linear combination of the others. Their beta is nan (R’s NA), and predictions take it as 0.

cause: Any

The cause of interest the subdistribution hazard is of.

center: Any

The covariate point the baseline is at, and phi relative to: zeros (Z = 0) by default, the covariate means for a fit with center=True (#463).

cif(x: ArrayLike, Z: ArrayLike) → NDArray

Cumulative incidence of the cause of interest at times x: 1 - exp(-Lambda0(x) * exp(beta'(Z - center))), the baseline Lambda0 that of a unit at center. Z is one covariate vector (a 1-D array or a single row), used at every time, or one row per time in x (row i with x[i]). The CIF is flat before the first event time and after the last (the baseline is a step function estimated only on the observed range). A missing (NaN) time or covariate gives nan in its place.

coefficients: NDArray

The coefficients (beta, the log subdistribution hazard ratios), their standard errors, Wald p-values and covariance.

classmethod from_dict(model_dict: dict) → FineGrayModel

Rebuild a Fine-Gray model from a to_dict() dictionary.

classmethod from_json(fp: str | PathLike) → Any

Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

maximum: str

What the fit reached, one of MAXIMUM_STATES (surpyval.utils.no_maximum), as its warnings say; "unknown" for a model restored from a dict saved without it.

phi(Z: ArrayLike) → NDArray

The subdistribution hazard multiplier \(e^{\beta' (Z - \text{center})}\), relative to a unit at center, where the baseline is (Z = 0 unless fitted with center=True), one value per row of Z (a scalar for a single covariate vector). It can overflow to inf on covariates far from center; cif() does not, as it combines it with the baseline on the log scale.

res: Any

The optimiser’s result (None on a restored model).

sf(x: ArrayLike, Z: ArrayLike) → NDArray

One minus the cumulative incidence (the cause-of-interest-free probability under the subdistribution).

to_dict() → dict

Serialise this fitted Fine-Gray model to a plain, JSON-serialisable dict.

Stores the coefficients and their covariance, plus the fitted subdistribution baseline cumulative-hazard step arrays (and, for a fit with center=True, the covariate center they are at, which makes the dict schema 2: a schema-1 reader would take the baseline for that at Z = 0), so the reloaded model reproduces cif/sf exactly and can still report the coefficient summary. The optimiser objects are not stored.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

  • fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.

  • with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

class surpyval.univariate.competing_risks.regression.competing_risks_proportional_hazard.CompetingRisksProportionalHazards

Bases: LinearPredictorMixin, SerialisableMixin

Competing-risks proportional-hazards regression.

Fits either a cause-specific proportional-hazards model (model="Cox", one Cox model per cause with the other causes treated as censored) or a Fine-Gray subdistribution-hazards model (model="Fine-Gray"). The naming follows the package convention (compare CompetingRisks and ProportionalHazards).

Call the class method CompetingRisksProportionalHazards.fit (or fit_from_df); it returns a fitted instance. Every prediction takes the covariates Z and, for one cause, its label event: an array in the fitted column order or, for a model fitted with fit_from_df, a DataFrame of the raw covariate columns (a formula is applied to it, as for CoxPH). The baselines are step functions, so the interp of sf, ff, Hf, hf and df takes only "step"; another value raises a ValueError. A fitted model can be saved with to_dict/to_json and restored with from_dict/from_json (or surpyval.from_dict).

Hf(x: ArrayLike, Z: ArrayLike, event: Any = None, interp: str = 'step') → NDArray

Cumulative hazard at x for covariates Z: one cause’s cause-specific cumulative hazard (event) or the all-cause sum (event=None). For a Fine-Gray model, the cumulative subdistribution hazard of event.

property aliased: NDArray

The columns of Z whose coefficients the data cannot determine in some cause’s fit (#476): a constant column, which the cause’s baseline hazard absorbs, or a linear combination of the others. Their coefficients are nan in that cause’s row of betas (R’s NA), and predictions take them as 0.

center: NDArray

The covariate point the baselines h0_e are at, and phi_e relative to: zeros (Z = 0) by default, the covariate means for a fit with center=True (#459, #463). The per-cause fits centre on the means either way.

cif(x: ArrayLike, Z: ArrayLike, event: Any) → NDArray

Cumulative incidence of cause event at x for covariates Z: the probability of failing from that cause by x with the other causes acting. The cause-specific (model="Cox") model builds it step by step from the causes’ hazard increments, as R’s survfit does for a multi-state coxph (the Aalen-Johansen estimate with each step’s matrix exponential), so the causes’ incidences sum to ff = 1 - exp(-H); the Fine-Gray model evaluates the subdistribution directly.

Z is one covariate vector (a 1-D array or a single row), used at every time, or one row per time in x (row i with x[i]), as for sf() and Hf(); a DataFrame of raw covariates for a model fitted with fit_from_df. event must be one of the fitted causes.

df(x: ArrayLike, Z: ArrayLike, event: Any = None, interp: str = 'step') → NDArray

hf * sf at x for covariates Z. Not available for a Fine-Gray model.

ff(x: ArrayLike, Z: ArrayLike, event: Any = None, interp: str = 'step') → NDArray

1 - sf at x for covariates Z. For a Fine-Gray model, the cumulative incidence of event.

classmethod fit(x: ArrayLike, Z: ArrayLike, e: ArrayLike, c: ArrayLike | None = None, n: ArrayLike | None = None, model: str = 'Cox', tie_method: str = 'efron', center: bool = False) → CompetingRisksProportionalHazards

Fit the competing-risks proportional-hazards model.

Parameters:

  • x (array like) – Failure or censoring times.

  • Z (ndarray like) – Covariate matrix, one row per observation. Rows with a missing (NaN) or infinite covariate are dropped, with a warning.

  • e (array like) – The cause of each failure; None (or NaN) for a right-censored observation.

  • c (array like, optional) – Censoring flags: 0 a failure, 1 right-censored. Derived from e if not given (a missing cause is censored). Left and interval censoring are not supported.

  • n (array like, optional) – Array of counts for each x. If data is provided as counts, then this can be provided. If None will assume each observation is 1.

  • model ({'Cox', 'Fine-Gray'}, optional) – 'Cox' (default) fits cause-specific proportional hazards – one Cox model per cause, the other causes treated as censored; 'Fine-Gray' fits one subdistribution-hazards model per cause.

  • tie_method (str, optional) – Tie handling for the 'Cox' path, passed to CoxPH.fit(). Default 'efron'.

  • center (bool, optional) – False (the default) reports each cause’s baseline at Z = 0; True at the covariate means, stored as model.center, with phi_e then relative to them. Passed to each cause’s fit (CoxPH.fit(), FineGray.fit), which centres either way; the default refuses covariates so far from 0 that the baseline there over- or underflows.

Returns:

model – A competing-risks proportional-hazards model. betas holds one row of coefficients per cause, in the order of event_idx_map (causes sorted); phi_e(Z, i) is cause i’s hazard multiplier, relative to a unit at center (where its baseline is: Z = 0 unless center=True). beta and phi (the sum of the per-cause coefficients and its multiplier) are kept for backward compatibility but are not a model quantity: every prediction uses the per-cause coefficients.

Return type:

CompetingRisksProportionalHazards

Examples

Two causes; the covariate doubles cause a’s hazard and leaves cause b’s alone:

>>> from surpyval.univariate.competing_risks import (
...     CompetingRisksProportionalHazards,
... )
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> Z = rng.binomial(1, 0.5, (200, 1)).astype(float)
>>> t_a = rng.exponential(1 / (0.1 * np.exp(0.7 * Z[:, 0])))
>>> t_b = rng.exponential(1 / 0.05, 200)
>>> t_c = rng.uniform(0, 20, 200)  # censoring times
>>> x = np.minimum(np.minimum(t_a, t_b), t_c).round(3)
>>> first = np.where(t_a < t_b, "a", "b")
>>> e = np.where(t_c < np.minimum(t_a, t_b), None, first)
>>> model = CompetingRisksProportionalHazards.fit(x, Z, e)
>>> model.betas.round(3)
array([[0.985],
       [0.005]])
>>> model.cif([5, 10], [[1]], "a").round(4)
array([0.59  , 0.7369])
classmethod fit_from_df(df: Any, x_col: str, e_col: str, Z_cols: str | list[str] | None = None, c_col: str | None = None, n_col: str | None = None, formula: str | None = None, model: str = 'Cox', tie_method: str = 'efron', center: bool = False) → CompetingRisksProportionalHazards

Fit a competing-risks proportional-hazards model from a pandas DataFrame.

Parameters:

  • df (pandas.DataFrame) – The data.

  • x_col (str) – Column of observed times.

  • e_col (str) – Column of event-type (cause) labels. Use None (or a blank/NaN cell) for a censored observation.

  • Z_cols (str or list of str, optional) – Covariate columns. Either Z_cols or formula must be given.

  • c_col (str, optional) – Column of censoring flags (0 observed, 1 right-censored).

  • n_col (str, optional) – Column of counts per row.

  • formula (str, optional) – A patsy/formulaic formula for the covariates, as an alternative to Z_cols.

  • model ({'Cox', 'Fine-Gray'}, optional) – Cause-specific proportional hazards or Fine-Gray subdistribution hazards. Default ‘Cox’.

  • tie_method (str, optional) – Tie handling for the model='Cox' path, passed to CoxPH.fit(): 'efron' (default), 'breslow', 'exact' or 'kalbfleisch-prentice' (alias 'kp').

  • center (bool, optional) – Report the baselines at the covariate means (model.center) instead of at Z = 0; see fit().

Returns:

The fitted model. Its prediction methods take a DataFrame of the raw covariate columns (the formula is applied to it) or a covariate array in the fitted column order.

Return type:

CompetingRisksProportionalHazards

Examples

A categorical covariate through a formula; the model predicts from a DataFrame of raw covariates, before and after saving:

>>> import numpy as np
>>> import pandas as pd
>>> import surpyval
>>> from surpyval.univariate.competing_risks import (
...     CompetingRisksProportionalHazards,
... )
>>> rng = np.random.default_rng(1)
>>> g = rng.choice(["a", "b", "c"], 300)
>>> rate = 0.1 * np.exp(np.select([g == "b", g == "c"], [0.8, -0.5]))
>>> t_a = rng.exponential(1 / rate)
>>> t_b = rng.exponential(1 / 0.05, 300)
>>> df = pd.DataFrame({
...     "time": np.minimum(t_a, t_b).round(3),
...     "cause": np.where(t_a < t_b, "a", "b"),
...     "g": g,
... })
>>> model = CompetingRisksProportionalHazards.fit_from_df(
...     df, "time", "cause", formula="g"
... )
>>> model.feature_names
['g[T.b]', 'g[T.c]']
>>> new = pd.DataFrame({"g": ["a", "b", "c"]})
>>> model.cif(np.full(3, 5.0), new, "a").round(4)
array([0.3752, 0.6449, 0.2551])
>>> restored = surpyval.from_dict(model.to_dict())
>>> bool(np.allclose(restored.cif(np.full(3, 5.0), new, "a"),
...                  model.cif(np.full(3, 5.0), new, "a")))
True
classmethod from_dict(model_dict: dict) → CompetingRisksProportionalHazards

Rebuild a competing-risks proportional-hazards model from a to_dict() dictionary.

classmethod from_json(fp: str | PathLike) → Any

Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

hf(x: ArrayLike, Z: ArrayLike, event: Any = None, interp: str = 'step') → NDArray

Cause-specific hazard increments at x for covariates Z: one cause’s (event) or the sum over causes (event=None). Not available for a Fine-Gray model.

maximum: str = 'unknown'

What the fit reached, one of MAXIMUM_STATES (surpyval.utils.no_maximum): the worst of the causes’ fits, as their warnings say; "unknown" for a model restored from a dict saved without it.

results: list | None

Each cause’s optimiser result, in event_idx_map order (None for a model restored from a dict: the optimiser objects are not serialised).

sf(x: ArrayLike, Z: ArrayLike, event: Any = None, interp: str = 'step') → NDArray

\(e^{-H}\) at x for covariates Z: the all-cause survival (event=None), which is one minus the sum of the causes’ cif(), or one cause’s net survival (the other causes treated as censoring). For a Fine-Gray model, 1 - cif.

to_dict() → dict

Serialise this fitted model to a plain, JSON-serialisable dict.

Stores the causes (event_idx_map), the per-cause coefficients betas and the per-cause baseline step arrays on the shared time grid x (with the covariate center they are at); for model="Fine-Gray" also each cause’s FineGrayModel (its own to_dict), from which the Fine-Gray predictions come. The reloaded model reproduces every prediction (cif, sf, ff, Hf, hf, df) for any Z. The optimiser results (results) are not stored.

Examples

>>> import numpy as np
>>> import surpyval
>>> from surpyval.univariate.competing_risks import (
...     CompetingRisksProportionalHazards,
... )
>>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> Z = [[0], [1], [0], [1], [0], [1], [0], [1], [0], [1]]
>>> e = ["a", "b", "a", None, "b", "a", "a", None, "b", "a"]
>>> model = CompetingRisksProportionalHazards.fit(x, Z, e)
>>> restored = surpyval.from_dict(model.to_dict())
>>> bool(np.allclose(restored.cif([5, 9], [1], "a"),
...                  model.cif([5, 9], [1], "a")))
True
to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

  • fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.

  • with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.