Comparison Tests and Validation Metrics

Tools that compare models or groups rather than fit one: hypothesis tests of whether groups share a survival (or cumulative-incidence) curve, the restricted-mean difference between two groups, automatic selection of the best-fitting distribution, and metrics that score how well a model’s predicted survival matches held-out data. All are importable directly from surpyval. The tests are explained with the estimators they build on in Non-Parametric Estimation (log-rank, restricted mean) and Competing Risks Analysis (Gray’s test); model selection in Parametric SurPyval Modelling; and the validation metrics in Regression Modelling with SurPyval.

Group-comparison tests

The (weighted, optionally stratified) log-rank test for comparing survival distributions across groups, and the result it returns:

surpyval.univariate.nonparametric.logrank.logrank(x: ArrayLike, Z: ArrayLike, c: ArrayLike | None = None, n: ArrayLike | None = None, weighting: str = 'log-rank', rho: float = 0, gamma: float = 0, strata: ArrayLike | None = None) → LogRankResult

The k-sample (weighted) log-rank test for the equality of survival distributions of right censored data.

At each distinct event time the observed number of events in each group is compared with the number expected under the null hypothesis that all groups share the same survival distribution. The weighted sums of these differences form a chi-squared statistic with k - 1 degrees of freedom.

Parameters:

x (array like) – Array of observations of the random variables.
Z (array like) – Array of group labels for each observation. Any hashable values other than NaN or None can be used. The test has one degree of freedom fewer than the number of groups with a positive expected number of events (as R’s survdiff): a group that is never at risk at an event time carries no information and is left out. If fewer than two groups have expected events there is nothing to compare, and the statistic is 0 with dof 0 and a p-value of 1.
c (array like, optional) – Array of censoring flags. 0 is observed and 1 is right censored. Left or interval censored data cannot be used with the log-rank test. If not provided assumes all values are observed.
n (array like, optional) – Array of counts for each x. If None assumes each observation is 1.
weighting (str, optional) –
The weighting to use at each event time. One of:
- ”log-rank”: weight 1 (the standard log-rank test; sensitive to proportional hazards alternatives),
- ”gehan”: weight r (a.k.a. Gehan-Breslow-Wilcoxon; emphasises early differences),
- ”tarone-ware”: weight sqrt(r),
- ”fleming-harrington”: weight S(t-)**rho * (1 - S(t-))**gamma where S is the pooled Kaplan-Meier estimate; rho > 0 emphasises early differences and gamma > 0 late differences.
Defaults to “log-rank”.
rho (scalar, optional) – The parameters of the Fleming-Harrington weighting. Only used when weighting is “fleming-harrington”. Defaults to 0, 0 (which is identical to the log-rank weighting).
gamma (scalar, optional) – The parameters of the Fleming-Harrington weighting. Only used when weighting is “fleming-harrington”. Defaults to 0, 0 (which is identical to the log-rank weighting).
strata (array like, optional) – Array of stratum labels, one per observation. When supplied the test is stratified: the observed-minus-expected numerators and their variances are accumulated within each stratum (risk sets never cross a stratum boundary) and summed before forming the statistic. This removes a nuisance factor – one whose baseline hazard differs across strata – from the comparison, so groups are only ever compared against others in the same stratum. The degrees of freedom are counted as for an unstratified test, from the expected events summed over the strata. NaN or None labels are refused.

Returns:

result – Object with the chi-squared statistic, the degrees of freedom dof, the p_value, the weighting and the number of strata.

Return type:

LogRankResult

Raises:

ValueError – If there are fewer than two groups, Z, c, n or strata does not have one entry per value, a group or stratum label is missing (NaN or None), the weighting is unknown, or the data have left or interval censoring.

Examples

>>> from surpyval import logrank
>>> x = [9, 13, 13, 18, 23, 28, 31, 34, 45, 48, 161,
...      5, 5, 8, 8, 12, 16, 23, 27, 30, 33, 43, 45]
>>> c = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1,
...      0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
>>> Z = [1] * 11 + [2] * 12
>>> res = logrank(x, Z, c=c)
>>> print(round(res.statistic, 2), round(res.p_value, 4))
3.4 0.0653

These are the AML maintenance data (aml in R’s survival package), for which R’s survdiff also gives 3.4 on 1 degree of freedom. The Gehan-Breslow weights put more weight on the early failures:

>>> res = logrank(x, Z, c=c, weighting="gehan")
>>> print(round(res.statistic, 3), round(res.p_value, 4))
2.723 0.0989

References

Klein, J. P. and Moeschberger, M. L. (2003), “Survival Analysis: Techniques for Censored and Truncated Data”, 2nd ed., Chapter 7.

class surpyval.univariate.nonparametric.logrank.LogRankResult(statistic: float, dof: int, p_value: float, weighting: str, strata: int | None = None)

Bases: object

Result of a (weighted) log-rank test.

statistic

The chi-squared test statistic.

Type:: float

dof

The degrees of freedom: the number of groups with a positive expected number of events, minus one (R’s survdiff convention).

Type:: int

p_value

The p-value of the test.

Type:: float

weighting

The weighting used for the test; for the Fleming-Harrington weights it includes rho and gamma, e.g. 'fleming-harrington(rho=1, gamma=0)'.

Type:: str

strata

The number of strata of a stratified test; None for an unstratified one.

Type:: int or None

Examples

logrank returns one; here comparing two groups of five:

>>> import numpy as np
>>> from surpyval import logrank
>>> x = np.array([3, 5, 7, 8, 10, 2, 3, 4, 4, 6])
>>> Z = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
>>> result = logrank(x, Z)
>>> result
Log-Rank Test
=============
Weighting        : log-rank
Statistic        : 3.69204
DoF              : 1
p-value          : 0.0546728
>>> round(result.p_value, 4)
0.0547

Gray’s test for comparing cumulative incidence functions across groups under competing risks:

surpyval.univariate.competing_risks.nonparametric.gray_test.gray_test(x: ArrayLike, e: ArrayLike, group: ArrayLike, event: Any, c: ArrayLike | None = None, n: ArrayLike | None = None, rho: float = 0.0) → GrayTestResult

Gray’s k-sample test comparing the cumulative incidence of one cause across groups.

Parameters:

x (array_like) – Event/censoring times (finite).
e (array_like) – Cause label per observation; a missing value (None, NaN or pandas NA) marks a right-censored observation, as for the competing-risks model classes.
group (array_like) – Group label per observation (two or more groups, no missing labels). Labels of different types (0 and "a") may be mixed.
event (scalar) – The cause (a label in e) whose cumulative incidence is compared across groups. The result keeps it as cause.
c (array_like, optional) – Censoring flag (0 event, 1 right-censored; left and interval censoring are not supported). If omitted it is derived from e (a missing cause is censored). If given, every row with c == 1 must have a missing cause and every other row a cause, otherwise a ValueError is raised.
n (array_like, optional) – Count weight per observation (default 1), each positive.
rho (float, optional) – Weight-family parameter: the per-time weight is (1 - F(t-))**rho with F the pooled cumulative incidence of event. 0 (the default) is the standard Gray test.

Returns:

(statistic, df, p_value, cause, groups); df is n_groups - 1 and a small p_value is evidence the groups’ cumulative incidence functions differ.

Return type:

GrayTestResult

Notes

This is Gray’s (1988) construction. Group g’s subdistribution risk set is R_g(t) = Y_g(t) (1 - F_g(t-)) / S_g(t-), from the group’s own at-risk count Y_g, Aalen-Johansen incidence F_g of the cause and all-cause Kaplan-Meier survival S_g; since Y_g / S_g estimates the group size times the group’s censoring survival, each group’s censoring is estimated separately, and the groups may be censored differently. The score is the weighted observed-minus-expected count sum_t w(t) (d_g(t) - R_g(t) d(t) / R(t)), with d the failures from the cause and R = sum_g R_g.

Its variance is Gray’s asymptotic estimate, computed as R’s cmprsk::cuminc computes it, and the statistic agrees with cuminc’s to rounding, ties included (#380). The score is linearised in each group’s counting-process martingales, of the cause and of the competing causes (the latter enter through F_g and S_g in the risk set), under the null hypothesis: the martingale variances are the failures each group is expected to have from the cause, and those it had from the competing causes, each with a correction for tied failures. It is not the hypergeometric (log-rank) variance, which ignores the variability of the estimated risk sets. The pooled incidence F^0 in the rho weight and the variance steps by d(t) / sum_g Y_g(t) / S_g(t-), the failures over the whole sample’s censoring-weighted size. cmprsk forms R_g the same way: Y_g(t) counts the rows censored at t and S_g is taken just before t, so Y_g(t) / S_g(t-) is the group size times its censoring survival just before t – a censoring tied with a failure counts after it, as in FineGray’s weights. Counts n enter as frequency weights (cmprsk has none).

Examples

Group 1 has twice group 0’s hazard of cause a:

>>> import numpy as np
>>> from surpyval import gray_test
>>> rng = np.random.default_rng(0)
>>> group = rng.binomial(1, 0.5, 200)
>>> t_a = rng.exponential(1 / (0.1 * np.exp(0.7 * group)))
>>> 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)
>>> res = gray_test(x, e, group, event="a")
>>> round(res.statistic, 3), res.df
(20.113, 1)
>>> bool(res.p_value < 0.001)
True

class surpyval.univariate.competing_risks.nonparametric.gray_test.GrayTestResult(statistic, df, p_value, cause, groups)

Bases: NamedTuple

Create new instance of GrayTestResult(statistic, df, p_value, cause, groups)

cause: Any: Alias for field number 3

df: int: Alias for field number 1

groups: list: Alias for field number 4

p_value: float: Alias for field number 2

statistic: float: Alias for field number 0

Restricted mean survival time

The two-group restricted-mean-survival-time difference (the per-model rmst method lives on the non-parametric model class):

surpyval.univariate.nonparametric.nonparametric.rmst_diff(model_a: NonParametric, model_b: NonParametric, tau: float | None = None, alpha_ci: float = 0.05) → dict

Compare the restricted mean survival time (RMST) of two groups.

The RMST-difference is the standard, assumption-light alternative to the hazard ratio when proportional hazards fails: it needs no PH assumption and reads directly as a difference in expected event-free time within the horizon tau.

Parameters:

model_a (NonParametric) – Two fitted non-parametric estimators (e.g. KaplanMeier.fit per group). They must carry a variance estimate (Greenwood).
model_b (NonParametric) – Two fitted non-parametric estimators (e.g. KaplanMeier.fit per group). They must carry a variance estimate (Greenwood).
tau (scalar, optional) – Common horizon. Defaults to the smaller of the two groups’ largest observed times, so both survival curves are supported by data up to tau (the standard choice). A larger tau is accepted without a warning: a curve is then held at its final value beyond its last observation, an extrapolation that is left to the caller to judge.
alpha_ci (scalar, optional) – Significance level for the interval and test (default 0.05).

Returns:

{"difference", "se", "lower", "upper", "p_value", "ratio", "rmst_a", "rmst_b", "tau"}. difference is RMST_a - RMST_b and ratio is RMST_a / RMST_b; se is the square root of the sum of the two groups’ variances (the groups are independent); lower and upper are difference +- z se; p_value is the two-sided z-test of difference == 0.

Return type:

dict

Raises:

ValueError – If either model has no variance estimate (fit_from_ecdf).

Examples

>>> import surpyval as sp
>>> a = sp.KaplanMeier.fit([2, 3, 4, 5, 6, 7])
>>> b = sp.KaplanMeier.fit([1, 2, 2, 3, 4, 5])
>>> res = sp.rmst_diff(a, b)
>>> print(res["tau"], round(res["difference"], 4), round(res["se"], 4))
5.0 1.1667 0.7233
>>> print(round(res["p_value"], 4))
0.1067

Automatic distribution selection

Fit every candidate continuous distribution and keep the one with the best information criterion:

Fit every candidate continuous distribution to the data and return the fitted model with the best value of metric.

The candidates are the fittable continuous univariate distributions with a regular likelihood (Beta, Exponential, ExpoWeibull, Gamma, Gumbel, Logistic, LogLogistic, LogNormal, Normal, Rayleigh and Weibull). A candidate the data lie outside the support of (a Beta for data outside (0, 1)) is passed over quietly; candidates whose fit fails are skipped and named in one warning. If every candidate fails, None is returned.

AIC, AIC_c and BIC compare maximised log-likelihoods, and their penalties assume a regular maximum: an interior point of the parameter space with a zero gradient near which the log-likelihood is quadratic, and a support that does not depend on the parameters. Two kinds of candidate break that, and are set aside – ranked only when no regular candidate fitted:

a family whose support ends are parameters (Uniform, Beta4). Its likelihood is highest with an end on the extreme observations, where the 2k penalty undercounts: on 50 Weibull(100, 2) draws the Uniform’s AIC beat the Weibull’s by 14. These are left out of the default candidates and tried only when named in include;
a fit that is not a verified maximum, as its maximum attribute records (see Parametric): "no finite maximum" (a Beta4 whose shape falls below 1, say), whose likelihood has no maximum, so its value, and every criterion made from it, means nothing – on [1, ..., 7] the Beta4 “won” with a log-likelihood of +24.8 against the Weibull’s -14.6 – or "unverified", a search that did not reach a verified maximum (an ExpoWeibull running towards a limit of its shapes, say), whose value is only where the search stopped.

When a candidate is set aside, one warning names it and says why; the warning its own fit would give that it is not a verified maximum is held back, replaced by that one.

Parameters:

x (array_like) – The observed event times (or intervals), in any of the formats fit accepts.
c (array_like, optional) – The censoring indicators.
n (array_like, optional) – The counts for each observation.
t (array_like, optional) – The truncation intervals.
metric (str, optional) – The model-selection criterion to minimise: "aic" (default), "aic_c", "bic" or "neg_ll".
include (iterable of str, optional) – Only try distributions with these names (matched without regard to case; a name that is not a candidate raises a ValueError). The only way to try the Uniform or the Beta4, which are then set aside (see above). Mutually exclusive with exclude.
exclude (iterable of str, optional) – Try every candidate except distributions with these names, checked in the same way. Mutually exclusive with include.

Returns:

The fitted model that minimises metric, or None when no candidate converged.

Return type:

Parametric or None

Raises:

ValueError – If candidates fitted but none has a finite metric – for "aic_c", when every candidate has at least as many parameters as observed failures minus one.

Examples

>>> from surpyval import fit_best
>>> import numpy as np
>>> np.random.seed(1)
>>> from surpyval import Weibull
>>> x = Weibull.random(50, 10, 2)
>>> model = fit_best(x, metric="bic")

Prediction-validation metrics

Right-censored-standard metrics for scoring a predicted survival function (Brier / integrated Brier score and Uno’s time-dependent AUC), plus the helper that builds a predicted-survival matrix from a fitted regression model or forest.

Prediction-validation metrics for right-censored survival predictors.

These score a predicted survival function against right-censored outcomes, handling censoring by inverse-probability-of-censoring weighting (IPCW):

brier_score() / integrated_brier_score() – the time-dependent Brier score of Graf et al. (1999): the IPCW-weighted squared error between the predicted S(t | Z) and the survival indicator, and its integral over a time grid. The standard scalar summarising calibration and discrimination together (lower is better).
auc_td() – Uno’s (2007) cumulative/dynamic time-dependent AUC: discrimination between subjects who have had the event by t and those still event-free, as a function of the horizon t (0.5 is chance, 1 is perfect).

All three are model-agnostic: they take a matrix of predicted survival probabilities. survival_probability() builds that matrix from a fitted model whose sf(x, Z) pairs x with the rows of Z (the parametric regression families, CoxPH, AdditiveHazards and the beta.ml forest).

The censoring survival G behind the weights is the reverse Kaplan-Meier with the events-first convention at ties (surpyval.utils.ipcw.censoring_survival() with ties="events_first"): an event and a censoring at the same time are ordered event first, as the data record them, so the event is not at risk of censoring there. An event at x_i is weighted by 1 / G(x_i-), the probability of its being observed, and a subject still event-free at the horizon t by 1 / G(t), as in Gerds and Schumacher (2006) and R’s pec. Without ties between event and censoring times the results agree exactly with scikit-survival’s brier_score, integrated_brier_score and cumulative_dynamic_auc. With such ties scikit-survival weights an event by 1 / G(x_i), which also discounts the censorings at x_i and over-weights the event.

References

Graf, E., Schmoor, C., Sauerbrei, W. and Schumacher, M. (1999), “Assessment and comparison of prognostic classification schemes for survival data”, Statistics in Medicine 18, 2529-2545.

Uno, H., Cai, T., Tian, L. and Wei, L. J. (2007), “Evaluating prediction rules for t-year survivors with censored regression models”, JASA 102, 527-537.

Gerds, T. A. and Schumacher, M. (2006), “Consistent estimation of the expected Brier score in general survival models with right-censored event times”, Biometrical Journal 48, 1029-1040.

surpyval.metrics.validation.auc_td(x: ArrayLike, c: ArrayLike, risk: ArrayLike, times: ArrayLike, x_train: ArrayLike | None = None, c_train: ArrayLike | None = None) → tuple[NDArray, NDArray]

Uno’s cumulative/dynamic time-dependent AUC (2007).

At each horizon t a case is a subject with the event by t (x_i \le t, \delta_i = 1) and a control is a subject still event-free (x_j > t). The AUC estimates the probability that a case is assigned a higher risk than a control, with cases IPCW-weighted by 1 / \hat G(x_i) to correct for censoring:

\[\widehat{AUC}(t) = \frac{\sum_{i,j} w_i\, \big(I(r_i > r_j) + \tfrac12 I(r_i = r_j)\big)\, I(\text{case}_i)\, I(\text{control}_j)} {\big(\sum_i w_i I(\text{case}_i)\big)\, \big(\sum_j I(\text{control}_j)\big)}.\]

Parameters:

x (array_like) – Observed times and censoring flags (0 event, 1 right censored).
c (array_like) – Observed times and censoring flags (0 event, 1 right censored).
risk (array_like, shape (n_samples, n_times)) – Risk scores where higher means earlier event. For a survival predictor use 1 - survival (see survival_probability()). A single column (or 1-D array) is broadcast across all times.
times (array_like) – Horizons at which to evaluate the AUC.
x_train (array_like, optional) – Data used to estimate the censoring distribution G. Defaults to the evaluation x / c.
c_train (array_like, optional) – Data used to estimate the censoring distribution G. Defaults to the evaluation x / c.

Returns:

times, auc – The horizons and the AUC at each. A horizon with no cases or no controls, or with a case whose weight is not identified (G = 0, see brier_score()), yields nan.

Return type:

ndarray

Examples

At t = 1.5 the one case (risk 0.9) outranks all four controls. At t = 2.5 the second case (risk 0.4) outranks only one of the three controls (risks 0.7, 0.5, 0.2), so the AUC is (3 + 1) / 6:

>>> from surpyval.metrics import auc_td
>>> x = [1.0, 2.0, 3.0, 3.0, 4.0]
>>> c = [0, 0, 0, 1, 0]
>>> risk = [0.9, 0.4, 0.7, 0.5, 0.2]
>>> auc_td(x, c, risk, [1.5, 2.5])[1].round(4)
array([1.    , 0.6667])

surpyval.metrics.validation.brier_score(x: ArrayLike, c: ArrayLike, survival: ArrayLike, times: ArrayLike, x_train: ArrayLike | None = None, c_train: ArrayLike | None = None) → tuple[NDArray, NDArray]

Time-dependent Brier score (Graf et al. 1999).

At each horizon t the Brier score is the IPCW-weighted mean squared error between the survival indicator I(T_i > t) and the predicted survival S(t | Z_i):

\[BS(t) = \frac1n \sum_i \Big[ \frac{S(t\mid Z_i)^2\, I(x_i \le t,\ \delta_i=1)}{\hat G(x_i-)} + \frac{(1-S(t\mid Z_i))^2\, I(x_i > t)}{\hat G(t)} \Big],\]

where \(\hat G\) is the Kaplan-Meier estimate of the censoring survival with the events-first convention at ties (see the module notes). Subjects censored before t contribute nothing (their status at t is unknown), nor do those censored at t; the IPCW weights correct for that loss. Lower is better.

Parameters:

x (array_like) – Observed times and censoring flags (0 event, 1 right censored) of the evaluation set.
c (array_like) – Observed times and censoring flags (0 event, 1 right censored) of the evaluation set.
survival (array_like, shape (n_samples, n_times)) – Predicted survival S(times[k] | Z_i); see survival_probability().
times (array_like) – Horizons at which to score, matching the columns of survival.
x_train (array_like, optional) – Data used to estimate the censoring distribution G. Defaults to the evaluation x / c. G is held at its last value beyond the largest training time. If it has reached 0 there, a horizon at which an evaluation row needs that zero (an event past the training censoring support, or a survivor at such a horizon) scores nan.
c_train (array_like, optional) – Data used to estimate the censoring distribution G. Defaults to the evaluation x / c. G is held at its last value beyond the largest training time. If it has reached 0 there, a horizon at which an evaluation row needs that zero (an event past the training censoring support, or a survivor at such a horizon) scores nan.

Returns:

times, bs – The horizons and the Brier score at each.

Return type:

ndarray

Examples

An event and a censoring tie at t = 3. The censoring survival is G = 1 before 3 and 1 - 1/2 from 3 (at risk of censoring at 3: the one censored there and the one still under observation). At the horizon 3.5 the three events are weighted by 1/G(x_i-) = 1 and the survivor by 1/G(3.5) = 2:

>>> from surpyval.metrics import brier_score
>>> x = [1.0, 2.0, 3.0, 3.0, 4.0]
>>> c = [0, 0, 0, 1, 0]
>>> S = [[0.8], [0.6], [0.5], [0.5], [0.3]]
>>> brier_score(x, c, S, [3.5])[1]  # (.64 + .36 + .25 + 2 * .49) / 5
array([0.446])

surpyval.metrics.validation.integrated_brier_score(x: ArrayLike, c: ArrayLike, survival: ArrayLike, times: ArrayLike, x_train: ArrayLike | None = None, c_train: ArrayLike | None = None) → float

Integrated Brier score: the Brier score averaged over times.

The trapezoidal integral of brier_score() over the time grid divided by its span, as scikit-survival’s integrated_brier_score. A single number summarising a survival predictor’s accuracy (lower is better); a model that predicts the true S(t | Z) scores below the marginal Kaplan-Meier reference. The grid need not be sorted (the columns of survival follow times); a single time, or a grid with no span, returns the mean Brier score. Parameters as for brier_score().

Examples

A Cox model of the Rossi recidivism data scored over the first 39 weeks, against the Kaplan-Meier curve that ignores the covariates:

>>> import numpy as np
>>> from surpyval import CoxPH, KaplanMeier
>>> from surpyval.datasets import load_rossi_static
>>> from surpyval.metrics import (
...     integrated_brier_score,
...     survival_probability,
... )
>>> df = load_rossi_static()
>>> x, c = df["week"].values, 1 - df["arrest"].values
>>> Z = df[["fin", "age", "prio"]].values
>>> times = [13, 26, 39]
>>> cox = CoxPH.fit(x, Z, c=c)
>>> S_cox = survival_probability(cox, Z, times)
>>> round(integrated_brier_score(x, c, S_cox, times), 4)
0.0998
>>> S_km = np.tile(KaplanMeier.fit(x, c).sf(times), (len(x), 1))
>>> round(integrated_brier_score(x, c, S_km, times), 4)
0.1038

surpyval.metrics.validation.survival_probability(model: Any, Z: ArrayLike, times: ArrayLike) → NDArray

Predicted survival matrix S(t | Z_i) from a fitted model.

Parameters:

model (object) – Any fitted model exposing sf(x, Z) where x is paired element-wise with the rows of Z (the parametric regression families, CoxPH, AdditiveHazards), or returning an (n_samples, n_times) grid (the beta.ml SurvivalTree and RandomSurvivalForest). Models whose sf takes a single covariate vector (BuckleyJames) are not supported: build their matrix row by row with model.sf(times, Z[i]).
Z (array_like or pandas.DataFrame) – Covariate matrix, one row per subject. A DataFrame is passed to model.sf as it is, so a model fitted with fit_from_df (Z_cols or a formula, string levels included) reads it by column name; an array is taken as numbers, one column per covariate.
times (array_like) – Evaluation times.

Returns:

survival – survival[i, k] is the predicted survival of subject i at times[k]. A subject with a missing covariate, or a missing time, gets nan where the model’s sf gives it.

Return type:

ndarray, shape (n_samples, n_times)

Examples

A formula fit with a string-valued factor is scored from a DataFrame:

>>> import numpy as np
>>> import pandas as pd
>>> from surpyval import WeibullPH
>>> from surpyval.metrics import survival_probability
>>> rng = np.random.default_rng(1)
>>> g = rng.choice(["a", "b"], 60)
>>> x = rng.weibull(1.5, 60) * np.where(g == "b", 5.0, 10.0)
>>> df = pd.DataFrame({"x": x, "g": g})
>>> model = WeibullPH.fit_from_df(df, x_col="x", formula="g")
>>> new = pd.DataFrame({"g": ["a", "b"]})
>>> survival_probability(model, new, [2.0, 5.0]).shape
(2, 2)

Harrell’s concordance index of any risk score, with Therneau’s treatment of tied event times by default, as R and lifelines (every regression model also has a concordance method that picks its family’s score):

surpyval.metrics.concordance.concordance_index(x: ArrayLike, c: ArrayLike, risk: ArrayLike, tie_tol: float = 1e-08, ties: str = 'therneau') → float

Harrell’s concordance index (C) of risk scores against right-censored outcomes.

A pair of subjects is usable when the one with the earlier time had the event; it is concordant when that subject also has the higher risk. C is the proportion of usable pairs that are concordant: 1 is a perfect ranking, 0.5 is chance. It measures discrimination only (a model can rank perfectly and still be miscalibrated; see brier_score()).

Parameters:

x (array_like) – Observed times.
c (array_like) – Censoring flags: 0 event, 1 right censored. Left and interval censoring are not supported.
risk (array_like) – Risk scores, higher meaning an earlier event: a hazard ratio or linear predictor, a cumulative hazard, 1 - sf at a horizon. To score predicted times or survival probabilities (higher meaning a later event), pass their negative.
tie_tol (float, optional) – Two scores within tie_tol of each other (or within a relative 1e-9, as math.isclose()) are tied. Default 1e-8.
ties ({"therneau", "harrell"}, optional) – How a pair of events at the same time counts. "therneau" (the default, as R’s survival::concordance and lifelines): it is not usable, since neither subject outlived the other. "harrell" (Harrell’s original definition): it is usable, and counts 1 if the scores are tied, else 0.5. Every other pair is treated the same way by both (see Notes).

Returns:

The concordance index; nan if a time or a score is missing.

Return type:

float

Raises:

ValueError – If the arrays differ in length, a flag is not 0 or 1, ties is not one of the conventions, or no pair is usable (every event is tied with, or later than, every other subject’s time).

Notes

Each pair is scored by the pairwise definition this function reproduces exactly (#276):

times x_i < x_j with an event at x_i (whatever c_j): 1 if risk_i > risk_j, 0.5 if the scores are tied, else 0;
equal times, both events: not usable under ties="therneau"; under ties="harrell" usable, 1 if the scores are tied, else 0.5;
equal times, one event and one censored: usable (the censored subject outlived the event), 1 if the event has the higher score, 0.5 on a tie, else 0;
equal times, both censored, or an earlier censored time: not usable.

The default is R’s (survival::concordance, Therneau) and lifelines’ (lifelines.utils.concordance_index, which scores predicted times, so its value for -risk is this one). The two conventions agree on data without tied event times.

The pairs are counted in \(O(n \log^2 n)\) array operations, not one by one: 50,000 subjects take a fraction of a second.

Examples

>>> from surpyval.metrics import concordance_index
>>> x = [1.0, 2.0, 3.0, 4.0, 5.0]
>>> c = [0, 0, 1, 0, 0]
>>> risk = [0.9, 0.5, 0.7, 0.6, 0.2]
>>> concordance_index(x, c, risk)
0.75

Two deaths at the same time are a pair only under Harrell’s convention:

>>> x = [1.0, 1.0, 2.0, 3.0]
>>> c = [0, 0, 0, 1]
>>> risk = [0.9, 0.5, 0.7, 0.2]
>>> concordance_index(x, c, risk)
0.8
>>> concordance_index(x, c, risk, ties="harrell")
0.75