A function for calculating the Brier score for a survival model.
brier_score(y_true = NULL, risk = NULL, surv = NULL, times = NULL)
loss_brier_score(y_true = NULL, risk = NULL, surv = NULL, times = NULL)a survival::Surv object containing the times and statuses of observations for which the metric will be evaluated
ignored, left for compatibility with other metrics
a matrix containing the predicted survival functions for the considered observations, each row represents a single observation, whereas each column one time point
a vector of time points at which the survival function was evaluated
numeric from 0 to 1, lower scores are better (Brier score of 0.25 represents a model which returns always returns 0.5 as the predicted survival function)
Brier score is used to evaluate the performance of a survival model, based on the squared distance between the predicted survival function and the actual event time, weighted to account for censored observations.
[1] Brier, Glenn W. "Verification of forecasts expressed in terms of probability." Monthly Weather Review 78.1 (1950): 1-3.
[2] Graf, Erika, et al. "Assessment and comparison of prognostic classification schemes for survival data." Statistics in Medicine 18.17‐18 (1999): 2529-2545.
library(survival)
library(survex)
cph <- coxph(Surv(time, status) ~ ., data = veteran, model = TRUE, x = TRUE, y = TRUE)
cph_exp <- explain(cph)
#> Preparation of a new explainer is initiated
#> -> model label : coxph ( default )
#> -> data : 137 rows 6 cols ( extracted from the model )
#> -> target variable : 137 values ( 128 events and 9 censored , censoring rate = 0.066 ) ( extracted from the model )
#> -> times : 50 unique time points , min = 1.5 , median survival time = 80 , max = 999
#> -> times : ( generated from y as uniformly distributed survival quantiles based on Kaplan-Meier estimator )
#> -> predict function : predict.coxph with type = 'risk' will be used ( default )
#> -> predict survival function : predictSurvProb.coxph will be used ( default )
#> -> predict cumulative hazard function : -log(predict_survival_function) will be used ( default )
#> -> model_info : package survival , ver. 3.5.8 , task survival ( default )
#> A new explainer has been created!
y <- cph_exp$y
times <- cph_exp$times
surv <- cph_exp$predict_survival_function(cph, cph_exp$data, times)
brier_score(y, surv = surv, times = times)
#> [1] 0.014504694 0.033959034 0.050857211 0.072372854 0.082965051 0.100349471
#> [7] 0.108958586 0.121159215 0.125346710 0.128121958 0.128948438 0.134527428
#> [13] 0.138530900 0.142016610 0.150285473 0.158946515 0.158878661 0.171125620
#> [19] 0.157656684 0.153866264 0.155959101 0.164717177 0.163641669 0.167682707
#> [25] 0.156722788 0.158009982 0.151887243 0.151107331 0.152375352 0.157368703
#> [31] 0.164691632 0.167253185 0.170942290 0.165724502 0.160732066 0.159359669
#> [37] 0.148453543 0.146427436 0.136663232 0.135249898 0.125593998 0.107657944
#> [43] 0.099012796 0.090316700 0.077607024 0.057225102 0.045652900 0.036427095
#> [49] 0.017892650 0.000299962
loss_brier_score(y, surv = surv, times = times)
#> [1] 0.014504694 0.033959034 0.050857211 0.072372854 0.082965051 0.100349471
#> [7] 0.108958586 0.121159215 0.125346710 0.128121958 0.128948438 0.134527428
#> [13] 0.138530900 0.142016610 0.150285473 0.158946515 0.158878661 0.171125620
#> [19] 0.157656684 0.153866264 0.155959101 0.164717177 0.163641669 0.167682707
#> [25] 0.156722788 0.158009982 0.151887243 0.151107331 0.152375352 0.157368703
#> [31] 0.164691632 0.167253185 0.170942290 0.165724502 0.160732066 0.159359669
#> [37] 0.148453543 0.146427436 0.136663232 0.135249898 0.125593998 0.107657944
#> [43] 0.099012796 0.090316700 0.077607024 0.057225102 0.045652900 0.036427095
#> [49] 0.017892650 0.000299962