Year: 2018

Abstract:
Personalized probabilistic forecasts of time to event (such as mortality) can be crucial in decision making, especially in the clinical setting. Inspired by ideas from the meteorology literature, we approach this problem through the paradigm of maximizing sharpness of prediction distributions, subject to calibration. In regression problems, it has been shown that optimizing the continuous ranked probability score (CRPS) instead of maximum likelihood leads to sharper prediction distributions while maintaining calibration. We introduce the Survival-CRPS, a generalization of the CRPS to the time to event setting, and present right-censored and interval-censored variants. To holistically evaluate the quality of predicted distributions over time to event, we present the scale agnostic Survival-AUPRC evaluation metric, an analog to area under the precision-recall curve. We apply these ideas by building a recurrent neural network for mortality prediction, using an Electronic Health Record dataset covering millions of patients. We demonstrate significant benefits in models trained by the Survival-CRPS objective instead of maximum likelihood.

heading: 1 INTRODUCTION
text: Having patient-specific predictions of time to an event such as mortality or bone fracture allows caregivers to make better informed decisions around patient care. Historically, prognosis scores have served as simple tools to stratify patient risk within a predefined time window (Lau et al., 2006; Cardona-Morrell & Hillman, 2015). However, such models tend to be too simplistic to be widely useful. They are often estimated from a large population of patients, and do not take into account patient-specific information to make individualized predictions (Yu et al., 2011). Meanwhile, the adoption of Electronic Health Record (EHR) systems over the past few decades has resulted in the collection of observational data on millions of patients spanning multiple years. This data enables development of patient-specific prediction models using machine learning. Such models are applicable to the larger patient population without being specific to a disease type or demographic, and this makes it possible to develop novel workflows in care delivery. For example, a high predicted probability of 3-12 month mortality could proactively notify palliative care teams of otherwise overlooked patients with end-of-life needs (Avati et al., 2017).
One way to obtain patient-specific survival predictions is to treat the problem as probabilistic classification; that is, training a binary classifier to predict outcomes of event by a particular time of interest (Avati et al., 2017; Rajkomar et al., 2018). However, such an approach has drawbacks. First, the model is specific to the time of interest it was trained upon – it is not straightforward how to take a model that was trained to predict probabilities of 1-year mortality and obtain predictions of 6-month mortality from it. Second, it is not usually possible to use data on all patients – for example, if a patient has only 3 months of history in the EHR system, it is neither possible to include that patient as a positive case nor a negative case in the 1-year mortality prediction task. Third, the process of constructing the data set implicitly conditions on the future outcome to select prediction times – evaluation is performed only at times looking backward from the event of interest. It has been shown that evaluation metrics can be overly optimistic relative to real world performance as a result (Sherman et al., 2017).
An alternative approach to the problem is survival prediction; that is, predicting time to event by estimating a distribution over future time. In this setting, traditional survival analysis methods such
as the Cox proportional hazards model (R., 1972) or accelerated failure time models (J.) are capable of handling data with censored observations (cases in which the event was not observed, but we know that the event did not occur up to a certain time). This addresses concerns raised by the classification approach, but there are a few nuances. First, traditional models typically make strong assumptions, such as proportional hazards or linearity. Second, challenges of low prevalence often arise when these methods are applied to large-scale observational datasets with heavy censoring, which is the case in real EHR data. Third, these survival analysis methods are typically evaluated as point estimates of risk, such as 10-year probabilities of events, rather than holistic measures of quality of the predicted distributions (Goff et al., 2014; Ranganath et al., 2016; Lee et al., 2018). Common metrics of evaluation include the C-statistic (Uno et al., 2007), log-`1 loss (Yu et al., 2011), and mean-squared-error (Katzman et al., 2018). While useful for the purposes of relative risk stratification, model comparisons made using point estimates leaves the quality of uncertainty in predicted distributions left unmeasured. If a point prediction is way off, it is penalized by the same amount whether the model was confident or not (that is, whether the predicted distribution had low or high variance).
In contrast, forecasts in the field of meteorology are typically made as full prediction distributions over all weather conditions given past and current observations (Gneiting et al., 2008). Evaluation of predictive performance is assessed by the paradigm of maximizing the sharpness of the predictive distribution, subject to calibration (Gneiting & Katzfuss, 2014). The intuition behind this paradigm is that probabilities have to be calibrated in order to be correct. However, that does not necessarily make them useful (one could always predict the marginal probability of an outcome without looking at the data, and still be well calibrated). The usefulness of a prediction distribution lies in its sharpness, or how well its mass concentrates. In summary, uncalibrated predictions (sharp or not) are useless, calibrated but non-sharp predictions are correct but less useful, and calibrated and sharp distributions are most useful.
To improve the sharpness of prediction distributions in the survival setting, we propose the use of proper scoring rules beyond maximum likelihood as the training objective. Proper scoring rules are known to measure calibration, and any model trained with a proper scoring rule will tend to maintain calibration (Gneiting & Katzfuss, 2014). For our purposes, we focus on the continuous ranked probablity score (CRPS) which has been used as an objective in the regression setting (Gneiting et al., 2008; Mohammadi et al., 2016; 2015). We generalize the CRPS for the survival setting, called Survival-CRPS, with right-censored and interval-censored extensions. Our work is one among many recent works Chapfuwa et al. (2018) in using non-MLE training objectives for survival prediction.
Summary of contributions. (1) We introduce the proper scoring rule Survival-CRPS, a generalization of CRPS, as an objective in survival prediction. We present its right-censored and interval-censored variants. (2) We propose a new metric, Survival-AUPRC, inspired by the paradigm of maximizing sharpness subject to calibration, to holistically measure the quality of a prediction distribution with respect to a possibly censored outcome. (3) We give practical recommendations for the mortality prediction task, by recommending use of the log-normal parameterization and interval censoring when training. (4) We employ the above techniques and demonstrate their efficacy by training a deep recurrent neural network model for accurate survival prediction of patient mortality using EHR data.

heading: 2 COUNTDOWN REGRESSION
text: Parametric survival prediction methods model the time to an event of interest with a family of probability distributions, indexed by the distribution parameters. The survival function, denoted S(t) : [0,∞)→ [0, 1], is a monotonically decreasing function over the positive reals with S(0) = 1 and limt→∞ S(t) = 0. The survival function represents the probability of an individual not having the event of interest up to a given time. Every survival function has a corresponding cumulative density function (CDF), denoted F (t) = 1− S(t), and probability density function (PDF), denoted f(t) = ddtF (t). The choice of the family of probability distributions implies assumptions made about the nature of the data generating process.
We denote the medical record of a patient i as ( {(x(i)t , a (i) t )}T (i) t=1 , d (i), c(i) ) , where t ∈ {1 . . . T (i)}
denotes the interaction number of this patient with the health record, x(i)t ∈ RD is the set of features corresponding to the t-th interaction, a(i)t ∈ R+ is age at time t, d(i) ∈ R+ is the age of death or age of last known (alive) encounter, and c(i) ∈ {0, 1} is a censoring indicator where c(i) = 0 means the age of death is d(i), and c(i) = 1 means the age of death is at least d(i). For each x(i)t we define the quantity y(i)t = d
(i) − a(i)t which represents the corresponding time to event or time to censoring. Traditional methods in survival analysis are designed to handle right-censored outcomes, but we observe that in many common scenarios outcomes are actually interval-censored. In the context of mortality prediction, for example, we know that humans almost never live past 120 years of age. Therefore, we assume that the true age of death lies between d(i) and A = 120 years, implying that the true time to death lies between 0 and T (i)t = A− a (i) t . We omit patient superscripts i and interaction subscripts t for succinctness where possible. We note that although our notation focuses on the problem of mortality prediction, our techniques generalize to any time to event task of interest.

heading: 2.1 SURVIVAL-CRPS: PROPER SCORING RULE OBJECTIVES
text: A scoring rule is a measure of the quality of a probabilistic forecast. A forecast over a continuous outcome is a probability density function over all possible outcomes, f̂ with corresponding cumulative density function F̂ . In reality, we observe some actual outcome, y. A scoring rule S takes a predicted distribution and an actual outcome, and returns a loss S(F̂ , y). It is considered a proper scoring rule if for all possible distributions G,
Ey∼F̂ [S(F̂ , y)] ≤ Ey∼F̂ [S(G, y)],
and strictly proper when equality holds if and only if F̂ = G (Gneiting et al., 2008). A proper scoring rule is one in which the expected score is minimized by the distribution with respect to which the expectation is taken. Intuitively, it encourages a model for being honest by predicting what it actually believes (Savage, 1971). When a proper scoring rule is employed as a loss function, it naturally forces the model to output calibrated probabilities (Gneiting & Katzfuss, 2014).
There are many commonly used proper scoring rules. Perhaps the most widely used is the logarithmic scoring rule, equivalent to the maximum likelihood objective:
SMLE(F̂ , y) = − log f̂(y).
In the presence of possibly censored data, we maximize the density for observed outcomes, and tail or interval mass for censored outcomes, and this is a proper scoring rule (Dawid & Musio, 2014).
SMLE-RIGHT(F̂ , (y, c)) = − log ( (1− c)f̂(y) + cŜ(y) ) SMLE-INTVL(F̂ , (y, c, T )) = − log ( (1− c)f̂(y) + c(F̂ (T )− F̂ (y))
) However, the logarithmic scoring rule is asymmetric, and harshly penalizes predictions that are wrong yet confident. This results in the training process becoming sensitive to outliers, and in general conservative in prediction-making (that is, hesitant to make sharp predictions) (Gneiting & Raftery, 2007).
(a) Uncensored
(b) Right-censored
(c) Interval-censored
Another proper scoring rule for forecasts over continuous outcomes is the CRPS (Gneiting et al., 2007), defined as
SCRPS(F̂ , y) = ∫ ∞ −∞ ( F̂ (z)− 1{z ≥ y} )2 dz = ∫ y −∞ F̂ (z)2dz + ∫ ∞ y (1− F̂ (z))2dz.
The CRPS has been used in regression as an objective function that yields sharper predicted distributions compared to maximum likelihood, while maintaining calibration (Gneiting et al., 2008). Intuition for the CRPS is better understood by analyzing the latter expression and noting that the two integral terms correspond to the two shaded regions in Figure 2a. The CRPS score is completely reduced to zero when the predicted distribution places all the mass on the point of true outcome, or equivalently, when the shaded region completely vanishes.
In the context of time to event predictions we propose the Survival-CRPS which accounts for the possibility of right-censored or interval-censored data:
SCRPS-RIGHT(F̂ , (y, c)) = ∫ y 0 F̂ (z)2dz + (1− c) ∫ ∞ y (1− F̂ (z))2dz,
SCRPS-INTVL(F̂ , (y, c, T )) = ∫ y 0 F̂ (z)2dz + (1− c) ∫ T y (1− F̂ (z))2dz + ∫ ∞ T (1− F̂ (z))2dz.
Note that when c = 0, both of the above expressions are equivalent to the original CRPS. Again, the intuition behind the Survival-CRPS is better understood by mapping each of the integral terms to the corresponding shaded region in Figure 2b and Figure 2c. The Survival-CRPS behaves like the original CRPS when the time of event is uncensored. For censored outcomes, it penalizes the predicted mass that occurs before the time of censoring and, if interval censored, also the mass after time by which the event must have occurred.
Both variants of the Survival-CRPS are proper scoring rules. They are special cases of the threshold weighted CRPS (Gneiting & Ranjan, 2011), where the weighting function is an indicator over the uncensored regions.

heading: 2.2 EVALUATION BY SHARPNESS SUBJECT TO CALIBRATION
text: Calibration assesses how well forecasted event probabilities match up to observed event probabilities. It is crucial in development of useful predictive models, especially for clinical decision-making. In binary prediction tasks without censoring, the Hosmer-Lemeshow test statistic (Hosmer et al., 2011) is commonly used to assess goodness-of-fit by comparing observed versus predicted event probabilities at quantiles of predicted probabilities. Extensions to account for censoring have been proposed (Grønnesby & Borgan, 1996; D’Agostino & Nam, 2003; Demler et al., 2015), but these methods apply only to predictions of dichotomous outcomes within a particular time frame (for example, 1-year risks of mortality).
There is no widely accepted method for evaluating the calibration of a set of entire prediction distributions, over multiple time frames, in the survival setting. D-calibration has been recently
proposed as a method for holistic evaluation (Andres et al., 2018), but relies on handling censored observations by assuming the true times to death are uniformly distributed past the times of censoring in the predicted distributions. When censored observations far outnumber the uncensored observations, this can lead to overly optimistic assessments of calibration. Another option is to evaluate observed event times on the cumulative density scale of predicted distributions, using a Kaplan-Meier estimate to account for censoring (Harrell, 2006). Again, this method has limitations in the heavily censored setting, as the quantiles in the tail of predicted cumulative densities have few uncensored observations, and will rarely yield well calibrated values.
We instead employ the following method to measure calibration. We compare predicted cumulative densities against observed event frequencies, evaluated at quantiles of predicted cumulative density. Right-censored observations are removed from consideration in quantiles that correspond to times after their points of censoring. Interval-censored observations are similarly removed from consideration in quantiles that correspond to times after censoring, but are additionally re-introduced in quantiles that correspond to times past the time by which the event must have occurred (in the mortality prediction task, this corresponds to 120 years of age).
Subject to calibration, we strive for prediction distributions that are sharp (i.e, concentrated). There are several metrics that could be used for measuring sharpness, such as variance or entropy. In the context of time to event predictions, holding two distributions with vastly different means to the same standard of variance or entropy would be unfair (for example, we would want lower variance for a prediction distribution with a mean of a day, compared to a mean of a year). Instead, we use the coefficient of variation (CoV) as a reasonable measure of sharpness. The CoV is defined as the ratio of one standard deviation to the mean, CoV(F̂ ) = √
Var[F̂ ] E[F̂ ] .

heading: 2.3 SURVIVAL-AUPRC: HOLISTIC EVALUATION OF TIME TO EVENT PREDICTIONS
text: Since sharpness is only a function of the predicted distributions, a measure of sharpness is only meaningful if the model is sufficiently calibrated. We now propose a metric that measures how concentrated the mass of the prediction distribution is around the true outcome, robust to miscalibration. The idea is similar to the area under a precision-recall curve, except here it is with respect to only one predicted distribution and one outcome. We first consider the uncensored case. As an analog to precision, we consider intervals relative to the true time of event, defined by ratios. For example, a region of precision 0.9 around an event that occurs at time y is the interval [0.9y, y/0.9]. Corresponding to this region of precision, the analogy to recall is the mass assigned by the predicted distribution over this interval, F̂ (y/0.9)− F̂ (0.9y). By exploring the full range of precision from 0 to 1, we obtain the Survival Precision Recall Curve. The area under this curve measures how quickly predicted mass concentrates around the true outcome as we expand the precision window.
Survival-AUPRCUNCENS(F̂ , y) = ∫ 1 0 (F̂ (y/t)− F̂ (yt))dt
The highest possible score is 1, when the predicted distribution is a Dirac δ function centered over the time of outcome. The lowest possible score is 0, when the predicted distribution is infinitely dispersed. The mean of all Survival-AUPRC scores across examples provides an overall measure of the quality of the predictions.
The aforementioned metric only applies when the event outcome is uncensored. In the case of censored observations, we use the same analogy but with the right end of precision intervals defined with respect to the time by which the event must have occurred in the interval-censored case, or infinity in the right-censored case.
Survival-AUPRCRIGHT(F̂ , y) = ∫ 1 0 (1− F̂ (yt))dt
Survival-AUPRCINTVL(F̂ , y, T ) = ∫ 1 0 (F̂ (T /t)− F̂ (yt))dt

heading: 2.4 RECURRENT NEURAL NETWORK MODEL
text: We apply our techniques to the mortality prediction task by building a multilayer recurrent neural network (RNN) with parameters θ, denoted RNNθ, that takes as input a sequence of features (in our case, information about a patient recorded in the EHR, for each interaction they had with the hospital) to predict parameters of a parametric probability distribution F̂ over time to death at each timestep (Figure 3). The network depends only on data from the current and previous timesteps, and not the future. The approach here is similar to the recently proposed Weibull time to event RNN (Martinsson, 2016), though we generalize to any choice of noise distribution. The distributions that are output in each timestep are used to construct an overall loss,
LRIGHT = N∑ i=1 T (i)∑ t=1 SRIGHT ( F̂ RNNθ { x (i) 1:t }, (y(i)t , c(i))) LINTVL = N∑ i=1 T (i)∑ t=1 SINTVL ( F̂ RNNθ { x (i) 1:t
}, (y(i)t , c(i), T (i)t )) , where N is the total number of patients in the training set, T (i) is the sequence length for patient i, and F̂RNNθ denotes the distribution parameterized by the output of the RNN. It is the sequential and monotonically decreasing predicted times to event that inspires the name Countdown Regression.

heading: 2.5 CHOICE OF LOG-NORMAL NOISE DISTRIBUTION
text: Common parametric distributions over time to event used in traditional survival analysis models include the Weibull, log-normal, log-logistic, and gamma (in order to be sufficiently expressive in model space, we seek distributions with at least two parameters). We choose the log-normal distribution because other distributions either involve the Gamma function in their density, or involve the pattern (y/p1)p2 , where p1 and p2 are parameters output from the neural network. We found these patterns to be highly sensitive to the inputs and to suffer from numerical instability issues.
For the log-normal distribution, a closed form expression for the CRPS is well known (Baran & Lerch, 2015). However, a closed form expression for the Survival-CRPS does not exist. We perform a change of variable to express the integral terms as finite integrals, and numerically approximate with the trapezoid rule. When training, we then back-propagate through the trapezoidal approximation. Details are given in Appendix B and C. We note that the approximation formulas are themselves proper scoring rules, as they are just weighted sums of Brier scores. Closed form expressions for the log-normal Survival-AUPRC are also given in Appendix D, E, and F.

heading: 3 EXPERIMENTS
text: We run experiments for the mortality prediction task to evaluate four different training objectives: maximum likelihood SMLE-RIGHT and SMLE-INTVL, and our scoring based loss SCRPS-RIGHT and SCRPS-INTVL. For interval censoring we assume a maximum lifespan of A = 120 years.
The neural network architecture is kept identical for all four experiments and implemented in PyTorch (Paszke et al., 2017). The input at each timestep consists of both real valued (for example, age of patient) and discrete valued (for example, ICD codes) data. Discrete data is embedded into a trainable real-valued 126-dimensional vector space, and vectors corresponding to the codes recorded at a given timestep are combined into a weighted mean by a soft self-attention mechanism. All real valued inputs are appended to the averaged embedding vector. We also provide the real valued features to every layer by appending them to the output of previous layer. The input vector feeds into a fully connected layer, followed by multiple recurrent layers. We use the Swish activation function (Ramachandran et al., 2017) and layer normalization (Ba et al., 2016) at every layer. Recurrent layers are defined using GRU units (Chung et al., 2014) with layer normalization inside. After the set of recurrent layers, the network has multiple branches, one per parameter of the survival distribution (for the lognormal, µ and σ2). The final layer in each branch has scalar output, optionally enforced positive with the softplus function, Softplus(z) = log(1 + exp(z)). We use Bernoulli dropout (Srivastava et al., 2014) at all fully connected layers, and Variational RNN dropout (Gal & Ghahramani, 2015) in the recurrent layers, with a dropout probability of 0.5. Optimization is performed using the Adam optimizer (Kingma & Ba, 2014), with a fixed learning rate of 1e-3.

heading: 3.1 DATA
text: We use electronic health records, with IRB approval, from the STARR Data Warehouse (previously known as STRIDE) for training and evaluation (Lowe et al., 2009). The Warehouse contains deidentified data for over 3 million patients (about 2.6% having a recorded date of death), spanning approximately 27 years. Each timestep in the sequence for a patient corresponds to all the data in the EHR for a given day. Only days having any data have a corresponding timestep in the sequence for each patient. We use diagnostic codes, medication order codes, lab test order codes, encounter type codes, and demographics (age and gender). Each code has a randomly initialized embedding vector as a trainable parameter. The set of 3 million patients, correspond to 51 million overall timesteps, and was randomly split in the ratio 8:1:1 into train, validation and test splits.

heading: 3.2 RESULTS
text: We first verify that all models are reasonably well-calibrated (Figure 4). Both the coefficient of variation and the Survival-AUPRC metrics suggest that the Survival-CRPS with interval censoring yields the sharpest prediction distributions (Table 1). Inspecting the mass past 120 years of age shows
that a naively trained prediction model with maximum likelihood can assign more than 75% of the mass to unreasonable regions, which is highly undesirable for the purpose of prediction. We note that this behavior is largely due to low prevalence of uncensored examples, which is typical in real world EHR data sets. As a result, the loss for the censored examples, which can be minimized by pushing mass as far away to the right as possible, dominates the small number of uncensored examples.
By predicting an entire distribution over time to death, the same model can be used to make classification predictions at various time points, highlighting the flexibility of our approach. When evaluated at 6 month, 1 year, and 5 year probabilistic predictions of mortality, our model remains well-calibrated with high discriminative ability (Appendix G, Figure 6).
In the interest of reproducibility, we run similar experiments on the publically available MIMIC-III dataset (Johnson et al., 2016) (Appendix H), and have published our corresponding source code 1.

heading: 4 RELATED WORK
text: Recent works have demonstrated potential to significantly improve patient care by making predictions with deep learning models on EHR data (Avati et al., 2017; Rajkomar et al., 2018), but these have been limited in treating the task as binary classification over a fixed time frame. Predicting survival curves instead of dichotomous outcomes has been explored (Yu et al., 2011; Lee et al., 2018), but only over finite length horizons. Deep survival analysis (Ranganath et al., 2016) has been proposed, but is limited to a fixed shape Weibull (bypassing the concerns we raised about stability, but limited in expressivity). The work by Yang et al. (2017) is similar to ours in terms of using an RNN and log-normal noise distribution, but limited to MLE training. DeepSurv (Katzman et al., 2018) uses a Cox proportional hazards model, which similarly makes a set of inflexible assumptions. The WTTERNN (Martinsson, 2016) model has a similar network architecture to ours, but is also limited to a Weibull distribution. All aforementioned models have only been optimized for maximum likelihood, instead of more robust proper scoring rules. The CRPS scoring rule has been used with Neural Networks in (Rasp & Lerch, 2018). Work in (Alaa & van der Schaar, 2017) also predicts full survival curves specific to a patient, but the use of GPs makes it difficult to scale to millions of patients. The work in Miscouridou et al. (2018) predict survival curves (both non-parametric, and flexible flow based parametric curves) while also handling missing covariates. Another recent work Chapfuwa et al. (2018) uses adversarial training for survival prediction. It has been shown that modern neural networks can be uncalibrated, and the work by Guo et al. (2017) suggests ways to improve calibration (though the work focuses on classification).

heading: 5 CONCLUSION
text: Better survival prediction models can be built by exploring objectives beyond maximum likelihood and evaluation metrics that assess the holistic quality of predicted distributions, instead of point estimates. We introduce the Survival-CRPS objective, motivated by the fact that the CRPS scoring rule is known to yield sharp prediction distributions while maintaining calibration. There are perhaps others scoring rules that work better, leaving avenues for future work. To evaluate, we introduce the Survival-AUPRC metric, which captures the degree to which a prediction distribution concentrates around the observed time of event. We demonstrate success in large-scale survival prediction by using a deep recurrent model employing a log-normal parameterization. By predicting an entire distribution for time-to-event, we circumvent issues associated with binary classification. Meanwhile, our model still yields accurate predictions when evaluated as dichotomous outcomes at particular times. The impact of having meaningfully accurate survival models is tremendous, especially in healthcare. We hope our work will be useful to those looking to build and deploy such models.

heading: ACKNOWLEDGMENTS
text: We thank the PyTorch team, particularly for the erf implementation that allowed use of the lognormal distribution. We thank Baran Sandor, Sebastian Lerch, Alejandro Schuler, Jeremy Irvin, and Russell Greiner for valuable feedback.
1http://github.com/anonymoususer/anonymous.git

heading: B. SURVIVAL-CRPS FOR LOG-NORMAL (RIGHT-CENSORED)
text: For a general continuous prediction distribution F , with actual time to outcome y ∈ R+, and censoring indicator c, we generalize the CRPS to the Right Censored Survival CRPS score as:
SCRPS-RIGHT(F, (y, c)) = ∫ ∞ −∞ (F (z)1{z ≤ log y ∪ c = 0} − 1{z ≥ log y ∩ c = 0})2dz
= ∫ ỹ −∞ F (z)2dz + (1− c) ∫ ∞ ỹ (F (z)− 1)2dz.
In the above expression F would generally be in the family of continuous distributions over the entire real line (eg. Gaussian). Alternately, one could also use a family of distributions over the positive reals (e.g log-normal), in which case the Survival CRPS becomes:
SCRPS-RIGHT(F, (y, c)) = ∫ ∞ 0 (F (z)1{z ≤ y ∪ c = 0} − 1{z ≥ y ∩ c = 0})2dz
= ∫ y 0 F (z)2dz + (1− c) ∫ ∞ y (F (z)− 1)2dz.
For the case of F being log-normal, the expression becomes SCRPS-RIGHT(FLN(µ,σ2), (y, c)) = ∫ y 0 Φµ,σ2(log z) 2dz + (1− c) ∫ ∞ y (1− Φµ,σ2(log z))2dz
= ∫ y 0 Φµ,σ2(log z) 2dz + (1− c) ∫ ∞ y Φ−µ,σ2(− log z)2dz
= ∫ y 0 Φµ,σ2(log z) 2dz + (1− c) ∫ 1/y 0 Φ−µ,σ2(log z) 2(1/z)2dz = I+µ,σ2(y) + (1− c)I − µ,σ2(y).

heading: C. SURVIVAL-CRPS FOR LOG-NORMAL (INTERVAL-CENSORED)
text: We further extend the Right Censored Survival CRPS to the case of interval censoring. This is particularly useful for all-cause mortality prediction where we assume a particular event must occur by time T . Using the same notations as before, the Interval Censored Survival CRPS is:
SCRPS-INTVL(F, (y, c, T )) = ∫ ∞ 0 (F (z)1{{z ≤ y ∪ c = 0} ∪ z ≥ T } − 1{{z ≥ y ∩ c = 0} ∪ z ≥ T })2dz
= ∫ y 0 F (z)2dz + (1− c) ∫ T y (F (z)− 1)2dz + ∫ ∞ T (F (z)− 1)2dz.
For the case of F being log-normal, the expression becomes SCRPS-INTVL(FLN(µ,σ2), (y, c, T )) = ∫ y 0 Φµ,σ2(log z) 2dz + (1− c) ∫ T y (1− Φµ,σ2(log z))2dz
+ ∫ ∞ T (1− Φµ,σ2(log z))2dz
= ∫ y 0 Φµ,σ2(log z) 2dz + (1− c) ∫ 1/y 1/T Φ−µ,σ2(log z) 2(1/z)2dz
+ ∫ 1/T 0 Φ−µ,σ2(log z) 2(1/z)2dz
= I+µ,σ2(y) + I − µ,σ2(T ) + (1− c) [ I−µ,σ2(y)− I − µ,σ2(T ) ] .

heading: D. SURVIVAL-AUPRC FOR LOG-NORMAL (INTERVAL-CENSORED)
text: We start with the most general case (interval censoring). For a general continuous prediction distribution F with an interval outcome [L,U ], we define the Survival-AUPRC as
Survival-AUPRC(F,L, U) = ∫ 1 0 [F (U/t)− F (Lt)] dt.
Specifically for the case of log-normal, where φ and Φ are PDF and CDF of N (0, 1) respectively, and L̃ = logL and Ũ = logU :
Survival-AUPRC(FLN(µ,σ2), L, U) = ∫ 1 0 [ FLN(µ,σ2)(U/t)− FLN(µ,σ2)(Lt) ] dt
= ∫ 1 0 [ FN (µ,σ2)(Ũ − log t)− FN (µ,σ2)(L̃+ log t) ] dt
(substituting s = log t) = ∫ 0 −∞ [ FN (µ,σ2)(Ũ − s)− FN (µ,σ2)(L̃+ s) ] esds
= [ FN (µ,σ2)(Ũ − s)− FN (µ,σ2)(L̃+ s) ] es ∣∣s=0 s=−∞
− ∫ 0 −∞ [ −fN (µ,σ2)(Ũ − s)− fN (µ,σ2)(L̃+ s) ] esds
= ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) +
∫ 0 −∞ [ fN (µ,σ2)(Ũ − s) + fN (µ,σ2)(L̃+ s) ] esds
= ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) +
∫ 0 −∞ fN (µ,σ2)(Ũ − s)esds+ ∫ 0 −∞ fN (µ,σ2)(L̃+ s)e sds
= ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) +
∫ 0 −∞ 1 σ φ ( Ũ − s− µ σ ) esds+ ∫ 0 −∞ 1 σ φ ( L̃+ s− µ σ ) esds(
substituting u = Ũ − s− µ
σ
) = ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) + ∫ Ũ−µ σ
∞
1 σ φ (u) eŨ−σu−µ(−σ)du+ ∫ 0 −∞ 1 σ φ ( L̃+ s− µ σ ) esds
( substituting v =
L̃+ s− µ σ
) = ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) + ∫ Ũ−µ σ
∞
1 σ φ (u) eŨ−σu−µ(−σ)du+
∫ L̃−µ σ
−∞
1 σ φ (v) evσ−L̃+µσdv
= ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) − eŨ−µ ∫ Ũ−µ σ
∞ φ (u) e−σudu+ e−L̃+µ
∫ L̃−µ σ
−∞ φ (v) evσdv( using ∫ ecxφ(x)dx = e c2 2 Φ(x− c) ) = ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) + U
eµ
[ e σ2 2 Φ(u+ σ) ]u= Ũ−µσ u=∞ + eµ L [ e σ2 2 Φ(v − σ) ]v= L̃−µσ v=−∞
= ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃) ) + U
eµ
[ e σ2 2 Φ ( Ũ − µ σ + σ ) − eσ 2 2 ] + eµ L [ e σ2 2 Φ ( L̃− µ σ − σ )] = ( FN (µ,σ2)(Ũ)− FN (µ,σ2)(L̃)
) + e σ2 2 [ eµ
L Φ ( L̃− µ σ − σ ) + U eµ ( 1− Φ ( Ũ − µ σ + σ ))] = Φµ,σ2(logU)− Φµ,σ2(logL)
+ e σ2 2
[ eµ
L Φ
( logL− µ
σ − σ
) + U
eµ Φ
( − logU − µ
σ − σ
)] .

heading: E. SURVIVAL-AUPRC FOR LOG-NORMAL (RIGHT-CENSORED)
text: For a general continuous prediction distribution F with an interval outcome [L,∞), we define Survival-AUPRC as
Survival-AUPRC(F,L) = ∫ 1 0 [1− F (Lt)] dt.
Specifically for the case of log-normal, where Φ is the CDF of N (0, 1), and L̃ = logL (following Appendix-D),
Survival-AUPRC(FLN(µ,σ2), L) = ∫ 1 0 [ 1− FLN(µ,σ2)(Lt) ] dt = 1− Φµ,σ2(L̃) + eµ+ σ2 2 L Φ ( L̃− µ σ − σ ) .

heading: F. SURVIVAL-AUPRC FOR LOG-NORMAL (UNCENSORED)
text: For a general continuous prediction distribution F with a point outcome y, we define Survival-AUPRC Survival-AUPRC(F, y) = ∫ 1 0 [F (y/t)− F (yt)] dt.
Specifically for the case of log-normal, where Φ is the CDF of N (0, 1), and ỹ = log y (following Appendix-D),
Survival-AUPRC(FLN(µ,σ2), y) = ∫ 1 0 [ FLN(µ,σ2)(y/t)− FLN(µ,σ2)(yt) ] dt
= e σ2 2
[ eµ
y Φ ( ỹ − µ σ − σ ) + y eµ Φ ( − ỹ − µ σ − σ )] .

heading: G. EVALUATION AS BINARY OUTCOME

heading: H. EXPERIMENTS ON THE MIMIC-III DATASET
text: On the MIMIC-III dataset (Johnson et al., 2016), we built a feed forward neural network that takes in 51015 hospital admissions in the dataset (70.1% censored) and makes predictions at the time of discharge. There is only one time of prediction per patient, so a recurrent model was not used. We removed admissions where the patient’s age was obfuscated or where the patient’s discharge time occurred after their recorded date of death. As features, we used demographics (age and gender) and embedded diagnostic codes into a 128-dimensional space.

Review summaries:
The paper proposes the use of Survival Continuous Ranked Probability score instead of maximum likelihood estimation for personalised probabilistic forecasts of time-to-event data, thus estimating a distribution over future time. The authors describe the evaluation their method using (1) proper scoring rule objectives; (2) evaluation of calibration using sharpness as a metric; (3) the survival precision recall curve. The authors then apply these techniques to predicting time-to-mortality using an RNN that takes EHR patient records to predict the probability of death at a given time point. It’s not clear how this is related to the Survival CRPS model or how this model is incorporated into the RNN.