Skip to content

shi-ang/SurvivalEVAL

Repository files navigation

SurvivalEVAL

PyPI


PyPI PyPI - Python Version License Maintenance

This python-based package contains the most completeness evaluation methods for Survival Algorithms (see paper). These evaluation metrics can be divided into 3 categories:

Visualization of the evaluation metrics

Installation

You can install the package via pip.

pip install SurvivalEVAL

Or if you want to do some modification by yourself. Clone the repo, cd into it and install it in editable mode (-e option). That way, these are no more need to re-install the package after modification.

git clone https://github.com/shi-ang/SurvivalEVAL.git
cd SurvivalEVAL
pip install -r requirements.txt
pip install -e . 

Quickstart Example

Install a survival analysis package, such as lifelines, and load the data. Then, you can use the following code to evaluate the model.

from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi

from SurvivalEVAL.Evaluator import LifelinesEvaluator

# Load the data
rossi = load_rossi()
rossi = rossi.sample(frac=1.0)

# Split train/test set
train = rossi.iloc[:300, :]
test = rossi.iloc[300:, :]
train_event_times = train.week.values
train_event_indicators = train.arrest.values
test_event_times = test.week.values
test_event_indicators = test.arrest.values

# Fit the model
cph = CoxPHFitter()
cph.fit(train, duration_col='week', event_col='arrest')

survival_curves = cph.predict_survival_function(test)

# Make the evaluation
evl = LifelinesEvaluator(survival_curves, test_event_times, test_event_indicators,
                          train_event_times, train_event_indicators)

cindex, _, _ = evl.concordance()

mae_score = evl.mae(method="Pseudo_obs")

mse_score = evl.mse(method="Hinge")

# The largest event time is 52. So we use 53 time points (0, 1, ..., 52) to calculate the IBS
ibs, (fig, ax) = evl.integrated_brier_score(num_points=53, draw_figure=True)

d_cal, details = evl.d_calibration(return_details=True)

# The target time for the single time probability prediction is set to 25
auc_score = evl.auc(target_time=25)
bs_score = evl.brier_score(target_time=25)
one_cal = evl.one_calibration(target_time=25)

See the Examples for more usage examples.

Point Prediction

Discrimination for point prediction

Concordance index (CI) identifies the “comparable” pairs of patients and calculates the percentage of correctly ranked pairs to assess a survival model’s performance. Given two predicted survival curves of a paired patients, it compares the predicted median/mean survival times and marks that pair as correct if the model's prediction about who died first matches with the reality.

Metric Name Description Code Paper Link
Harrell's CI Original Harrell's CI, using comparable pair for concordant checking. evl.concordance(method="Harrell") Harrell et al.
Uno's CI Uno's CI adds a (squared) IPCW weighting to each comparable pair to the overall calculation. To be implemented Uno et al.
Margin CI Calculate the marginal survival time (and proxy time) to substitute the censoring time for censored patients, then calculate CI as everyone is uncensored. eval.concordance(method="Margin") Kumar et al.

Mean Absolute Error, Mean Squared Error and Root Mean Squared Error

One straightforward metric would be “MAE” – the absolute difference between the actual and predicted survival times (e.g. median of a curve). This requires using the “actual survival time”, which is trivial for uncensored instances, but problematic for censored individuals. This python package implemented MAE loss metrics using different ways of handling censored instances. Here we list three of them:

  1. Uncensored simply discards all the censored individuals and compute the MAE for all the uncensored instances.
  2. Hinge calculates the early prediction error. For a censored instance, if the predicted survival time is smaller than the censored time, then MAE = censor_time - predict_time. If the predicted survival time is equal or larger than the censored time, then MAE = 0. Note that the standard Hinge method requires the Weighted parameter to be set to False.
  3. Pseudo_obs “de-censors” the censored patients, using pseudo-observation method (by estimating the contribution of a censored subject to the whole Kaplan-Meier distribution). Then it calculates the MAE between de-censoring time and the predicted survival time, just like the normal way. Note that the standard Pseudo_obs method requires the Weighted parameter to be set to True.

Mean squared error (MSE) is another metric to measure the difference between the actual and predicted survival times. Similar to MAE, mean squared error (MSE) also has multiple ways to handle censored instances, similar to MAE. We also have root mean squared error (RMSE) for each of the different ways.

Metric Name Description Code Paper Link
MAE-Uncensored Calculate MAE for uncensored instances only. evl.mae(method="Uncensored") N/A
MAE-Hinge For censored instances, Hinge calculates the early prediction error. For a censored instance, if the predicted survival time is smaller than the censored time, then MAE = censor_time - predict_time. If the predicted survival time is equal or larger than the censored time, then MAE = 0. Note that the standard Hinge method requires the Weighted parameter to be set to False. evl.mae(method="Hinge") Shivaswamy et al.
MAE-Margin It "de-censors” the censored patients, using the margin time (mean of the conditional KM curve given the censoring time). evl.mae(method="Margin") Haider et al.
MAE-PO It "de-censors” the censored patients, using pseudo-observation method. evl.mae(method="Pseudo_obs") Qi et al.
MSE MSE has the same user-interface with MAE. evl.mse(method=method) N/A
RMSE RMSE has the same user-interface with MAE. evl.rmse(method=method) N/A

Reliability

Ideally, a reliable model should produce similar predictions (in this context, similar point prediction) to the true labels. Here we use the log-rank test, to test whether the predicted survival times are statistically different than the observed labels (subject to censoring).

Metric Name Description Code Paper Link
Log-rank Test Use log-rank test to check whether the predicted survival times are statistically different than the observed labels (subject to censoring). A model with p-value higher than 0.05 can be considered as reliable. evl.log_rank() Mantel
Weighted Log-rank Test Weighted log-rank test is a generalization of the log-rank test that accounts for different weights for different time regions. This can be useful when the proportional hazard assumption is violated, or you care about specific region of the data. evl.log_rank(weighting=["wilcoxon"/"tarone-ware"/"peto"/"fleming-harrington"]) Fleming&Harrington

Single Time Probability Prediction

Discrimination for single time probability prediction

Area Under the Receiver Operating Characteristic (AUROC) is a metric to measure the performance of a single time probability prediction. It is the area under the receiver operating characteristic curve, which is the plot of the true positive rate against the false positive rate at various threshold settings. In the survival analysis, the single time probability prediction is the prediction of the survival probability at a specific time point. And the true label is whether the patient died at that time point. AUROC excludes the censored instances whose censoring time is earlier than the target time point.

Metric Name Description Code Paper Link
AUC Calculate AUC/AUROC for uncensored instances only. evl.auc() N/A
Uno AUC Uno's CI adds a IPCW weighting to each pairs to the overall calculation. To be implemented Uno et al.

Errors between the predicted probability and survival status

The Brier score (BS), at a specific time-point, is computed as the mean squared error between the observed event (binary indicator variable) and the predicted event probability at that time-point. It is meaningful in the sense that the square root of the Brier score is the distance between the observed and predicted event on the probability scale.

Metric Name Description Code Paper Link
Plain BS Calculate BS for uncensored instances only. evl.brier_score(IPCW_weighted=False) N/A
IPCW BS Adding a IPCW weighting to uncensored instances after the target time point. evl.brier_score(IPCW_weighted=True) Graf et al.

Calibration for single time probability prediction

Calibration measures the alignment between the predicted survival probability and the observed survival status at a specific time point. The most common Hosmer-Lemeshow (HL) goodness-of-fit test is used to check the calibration of a model at a specific time point. Models with p-value higher than 0.05 can be considered as well-calibrated model at that time.

Metric Name Description Code Paper Link
Uncensored HL test Standard Hosmer-Lemeshow test for uncensored part of the dataset. evl.one_calibration(method="Uncensored") Hosmer&Lemeshow
DN HL test DN's extension of HL test, using Kaplan-Meier estimation to approximate the observed probability. evl.one_calibration(method="DN") D'Agostino&Nam
Integrated Calibration Index Use non-linear model to fit a smooth calibration curve and calculate the integrated error. eval.integrated_calibration_index() Austin et al.

Survival Distribution Prediction

Discrimination for survival distribution prediction

Time-dependent concordance index (td-CI) is a generalization of the concordance index to the survival distribution prediction. It is calculated by integrating the concordance index over times.

Metric Name Description Code Paper Link
td-CI Calculate td-CI for uncensored instances only. To be implemented Antolini et al.

Errors for individual survival function and the heaviside step function

Integrated Brier Score (IBS) measures the squared difference between the predicted survival curve with the Heaviside step function of the observed event. IBS can be viewed as the integration of the (single-time) Brier score across all the time points. A smaller IBS value is preferred over the larger value. This python implementation uses IPCW weighting to handle the censored instances. Please refer to Assessment and Comparison of Prognostic Classification Schemes for Survival Data for the detail of IPCW weighting. Please also note that IBS is also similar to the Continuous Ranked Probability Score (CRPS), except (1) the IPCW weighting, and (2) squared error instead of absolute error.

Metric Name Description Code Paper Link
Plain IBS Calculate IBS for up to the censoring time, for uncensored instances. evl.integrated_brier_score(IPCW_weighted=False) N/A
CRPS Calculate integrated absolute error between survival function and heaviside function. evl.crps() Avati et al.
IPCW IBS Adding a IPCW weighting to uncensored instances after the target time point. evl.integrated_brier_score(IPCW_weighted=True) Graf et al.

Calibration for survival distribution prediction

Calibration for the entire survival distribution is a more complex task than the single time calibration. Some residuals are proposed to check the calibration of the entire survival distribution. And Haider et al. proposed distribution calibration (D-calibration) test for determining if a model that produces ISDs is meaningful. D-calibration splits the time-axis into a fixed number of intervals and compares the actual number of events with the predicted number of events within each interval. A well D-calibrated model is the one where the predicted number of events within each time interval is statistically similar to the observed number. Models with p-value higher than 0.05 can be considered as well-calibrated model across the survival distribution.

Metric Name Description Code Paper Link
Cox-Snell Residual Calculate the cumulative hazard function at the observed time, and then fit a Nelson-Aalen estimator to check the goodness-of-fit. evl.residuals(method="CoxSnell") Cox&Snell
Modified Cox-Snell Residual v1 Calculate the cumulative hazard function at the observed time, add 1 for censored subjects. evl.residuals(method="Modified CoxSnell-v1") Collett, Chapter 4
Modified Cox-Snell Residual v2 Calculate the cumulative hazard function at the observed time, add ln(2) for censored subjects. evl.residuals(method="Modified CoxSnell-v2") Crowley&Hu
Martingale Residual Calculate the Martingale residual, which should have a mean of 0 evl.residuals(method="Martingale") Barlow&Prentice
Deviance Residual Calculate the Deviance residual, which should follow a normal distribution N(0, 1). evl.residuals(method="Deviance") Therneau et al.
Distribution Calibration Check whether the survival probability at the event time follows the Uniform distribution betwen 0 and 1. eval.d_calibration() Haider et al.
KM Calibration Check whether the averaged predicted survival distribution is aligned with KM estimator. eval.km_calibration() Chapfuwa et al.

Non-Parametric Methods

This package also provides some non-parametric methods to estimate the marginal survival function, without considering the covariates.

The non-parametric methods are included in the SurvivalEVAL.NonparametricEstimator module.

Method Description Code Paper Link
Kaplan Meier The well known Kaplan Meier (KM) estimator that directly estimate the survival function. SingleEvent.KaplanMeier() Kaplan&Meier
Nelson Aalan Nelson Aalan (NA) estimate the cumulative hazard function directly, then transform to survival function. SingleEvent.NelsonAalen() Aalan
Copula Graphic The Copula Graphic (CG) estimator estimate the survival function under the dependent censoring assumption, with known type of copula and parameter. SingleEvent.CopulaGraphic() Emura&Chen
Turnbull The Turnbull estimater estimate the survival function for interval censoring dataset. SingleEvent.TurnbullEstimatorLifelines() Turnbull
Aalan Johansen The Aalan Johansen (AJ) estimator the cumulative incidence function (CIF) for competing risk dataset. To be implemented Aalan&Johansen

Citing this work

We recommend you use the following to cite SurvivalEVAL in your publications:

@article{qi2024survivaleval,
year = {2024},
month = {01},
pages = {453-457},
title = {{SurvivalEVAL}: A Comprehensive Open-Source Python Package for Evaluating Individual Survival Distributions},
author={Qi, Shi-ang and Sun, Weijie and Greiner, Russell},
volume = {2},
journal = {Proceedings of the AAAI Symposium Series},
doi = {10.1609/aaaiss.v2i1.27713}
}

About

The most comprehensive Python package for evaluating survival analysis models.

Topics

Resources

License

Stars

Watchers

Forks

Packages

No packages published

Languages