Prediction of Probability of Survival in Critically Ill Patients Optimizing the Area under the ROC Curve

1 Introduction

The available models for predicting outcomes in intensive care units (ICU) are usually scoring systems that estimate the probability of hospital mortality of critically ill adults. This is the case of APACHE (Acute Physiology And Chronic Health Evaluation) [Knaus et al., 1991], SAPS (Simplified Acute Physiology Score) [Le Gall et al., 1984], and MPM (Mortality Probability Models) [Lemeshow et al., 1993]. The score functions of these predictors were induced from data on thousands of patients using logistic regression. The data required by these systems come from monitoring devices, clinical analysis, and demographic and diagnostic features of patients. So, APACHE III includes age, 16 acute physiologic variables that use the worst value from the first 24 hours in the ICU (temperature, heart rate, blood pressure, respiratory rate, oxygenation, acid-base status, serum sodium, serum blood urea nitrogen, serum creatinine, serum albumin, serum bilirrubin, serum glucose, white cell count, hematocrit, itemized Glasgow Coma Scale score, and urine output), preexisting functional limitations, major comorbidities, and treatment location immediately prior to ICU admission.

∗The research reported here is supported in part under grant TIN2005-08288 from the MEC (Ministerio de Educaci ´on y Ciencia of Spain). The authors acknowledge the work of the Grecia Group in the collection of data.

These prognostic models are mainly used to measure the efficiency of ICU treatments. The risk stratification of patients allows comparison of the observed outcomes versus accepted standards provided by score functions. ICU assessment is very important since it is estimated that end-of-life care consumes 10% to 12% of all healthcare costs. Moreover, in 2001 the average daily cost per patient in ICUs was about $3000 in the USA [Provonost and Angus, 2001]. On the other hand, the literature also shows that prognoses have constituted an important dimension of critical care, as patients and their families seek predictions about the duration and outcome of illness [Lemeshow et al., 1993].

In this paper we propose a new method for learning probabilities that will be tested on the probabilities of survival in ICU patients. The method makes intensive use of the socalled Support Vector Machines (SVM), a powerful family of algorithms for learning classification and regression tasks. When used for binary classification, SVM learn hypotheses that return continuous numbers: positive values for cases of one class, and negative for the other class.

On the other hand, to measure the performance of predictions in medicine, and in general when classes are very unbalanced, the misclassification rate (or accuracy) is usually inadequate. Frequently, the Area Under a receiver operating characteristic (ROC) Curve (AUC for short) is used. This amount can be interpreted as the degree of coherence between a continuous output (such as the probability, or the continuous output of an SVM) and a binary classification. It is important to emphasize that that coherence is established in terms of orderings. For this purpose, continuous outputs or scores are used to rank available cases, while classes in the ICU problem are codified by '+1' when the patient has survived, and '-1' otherwise.

In this context, Hanley and McNeil [1982] showed that the AUC is the probability of a correct ranking; in other words, it is the probability that a randomly chosen subject of class '+1' is (correctly) ranked with greater output than a randomly chosen subject of class '-1'. Therefore, AUC coincides with the value of the Wilcoxon-Mann-Whitney statistic.

Additionally, there are other measures of the goodness of probability estimations; for instance, the Brier score is the average of quadratic deviations of true and predicted probabilities. The relation between AUC and Brier scores was studied in [Ikeda et al., 2002]. However, the relationship found is guaranteed only under very restrictive conditions that are difficult to check in real world cases. Moreover, the relationship is not always direct even in the case study reported in the experimental section of this paper.

To learn a probability distribution using SVM, it is crucial to transform their scores or continuous outputs into probabilities. But this is what a method presented by Platt [2000] does. The core idea is to fit a sigmoid, using a maximum likelihood procedure. The novelty of the proposal reported in this paper is that we postulate that to compute Platt's sigmoid it is better to look for an optimum AUC first than to minimize the error rate with a classification SVM. For this reason, in Section 2.3 we shall discuss how to optimize the AUC with a Support Vector method [Herbrich et al., 2000; Joachims, 2005].

The rationale behind our proposal is that the quality of the sigmoid fit depends on the quality of the ranking of the scores. If most of the cases with a higher score than a given one of class y have a class greater than y, then the task of the sigmoid can be easily accomplished, and the performance of the final probability is nearly optimal.

At the end of the paper we provide experimental evidence in favor of our proposal, comparing it with other alternative approaches. For this purpose, we used data collected in general ICUs at 10 hospitals in Spain, 6 of which include coronary patients, while the other 4 do not treat coronary diseases. The total number of patients considered in our study was 2501, 19.83% of whom did not survive.

2 Predicting probabilities

In this section we shall start off by reviewing a standard method for learning probabilities based on Support Vector Machines to then present our proposal. But first of all we must realize that the performance of classification learners is not satisfactory in the ICU problem; otherwise, nobody would turn to probabilities. This is a general situation in medicine, as well as in other fields; accurate crisp predictions are difficult to make, but some useful knowledge can be drawn from data.

The section will end with the description of a straightforward approach for learning probabilities using regression. This method will be used as a baseline for measuring the merits of the other options.

2.1 The goodness of probability predictions

Let S = {(x1, y1), . . . , (xn, yn)} be a training set for a learning task in which a function (or hypothesis) is sought that is able to return outputs yi from entries xi of an input space X. An important issue when we are learning is to fix the way in which we are going to measure the quality of the result. In fact, given S formally, the aim of learning is to find a hypothesis h (from a given hypothesis space) that minimizes the average loss extended over the set of independently identically distributed (i.i.d.) test sets S′ , usually represented by Δ(h, S′ ).

In the ICU problem, training and test examples have no probability attached, they are labeled with +1 or -1.

Therefore, we shall assume that the true class probability, P rtrue (y = +1|x), is 1 when the class of x, y, is +1 and 0 otherwise.

In general, when predictions are discrete probability distributions, there is basically one standard loss function: the average quadratic deviation. If there are two possible outputs, the probability loss is given by

where the hypothesis h returns the estimation of the probability h(x) = Pr(y = +1|x), and pi stands for the observed probability of the i-th case, pi = Prtrue (y = +1|xi).

The measurement in Equation (1) is frequently used in medicine and meteorology, and is known as the Brier [1950] index or score. If the number of possible outputs is greater than two, the estimated probabilities can be seen as a vector, and the Mean Square of the Euclidean (MSE) distance from predicted and observed probabilities is then used; see, for instance [Melville et al., 2005]. It can be seen that, in the ICU problem, MSE is 2 times the Brier score.

2.2 Optimizing accuracy plus a sigmoidal transformation

The straightforward approach to the ICU problem is a binary classification SVM followed by a sigmoid estimated using Platt's method [2000]. Thus, given the training set S, we can use a transformation φ defined from input entries in X into a feature space H, where classes should be mostly separable by means of a linear function. As is well known, H must have an inner product 〈, 〉, and

is called the kernel function of the transformation. We shall use the rbf kernel that is defined by

The work of the SVM consists in solving the following convex optimization problem:

Then, the classification is accomplished by the hypothesis

It can be seen that the kernel and the vector α = (αi : i = 1, . . . , n) of Lagrange multipliers define the implementation of Function (5) computed from input space entries x as follows:

According to (4), the aim of this function is to maximize the margin (between classes) and to minimize the training loss. In fact, the sum of the so-called slack variables, ∑ni=1 ξi, is an upper bound of misclassifications of (6) on the training set. It is acknowledged that the Function (6) so achieved has good classification accuracy on unseen cases.

In order to compute the probabilistic outputs, we get rid of the sign function, and we only consider the continuous outputs

Platt's method then fits a sigmoid to estimate probabilities:

Figure 1 depicts the fit of this sigmoid to the dataset of all patients (2501) at all the available hospitals. Notice that the fac values follow a bell-shape distribution with most individuals having positive values, which means that they have a survival prediction.

2.3 Optimizing the AUC first

When classification predictions are made comparing the values returned from patients' descriptions x by a rating function with a threshold, as in classification SVM (see Equation (5)), then the performance of these predictions can be assessed using the AUC. According to its probabilistic interpretation, the complementary of this amount (1-AUC) can be used as a loss function. Thus, if g is a hypothesis, its loss evaluated on a test set S′ is

Let us stress that the explicit objective of SVM presented in the preceding section is not to minimize Equation (9). [Cortes and Mohri, 2004] provide a detailed statistical analysis of the difference between maximizing the AUC and minimizing the error rate in binary classification tasks.

Herbrich et al. [2000] presented a direct implementation that solves a general ranking problem that is applicable to maximizing the AUC. The core idea is that if a hypothesis f : φ(X ) → R is linear and has to fulfill that f(φ(xi)) > f(φ(xj)), since yi > yj, then

Notice that this statement converts ordering constraints into classification constraints (with one class), but now the input space is X ×X and each pair (xi, xj) is represented by the difference φ(xi) - φ(xj). According to this approach, the aim is to find a hypothesis f(x) = 〈w, φ(x)〉 such that w solves the following convex optimization problem:

For each x of the input space, the hypothesis so found returns

where αi,j are again the Lagrange multipliers computed by the optimizer.

Unfortunately, this approach leads to dealing with one constraint for each element of the dataset

whose size is the number of positive (class +1) examples times the number of negatives, #pos×#neg, i.e. O(n2 ) when the size of S is only n. This means that some applications become intractable, although the approach (or a simplified version of it) has been successfully used on other occasions [Joachims, 2002; Bahamonde et al., 2004].

To alleviate the difficulties caused by the size of data sets, it is not straightforward to reformulate Herbrich's approach as an optimization problem with a small number of constraints. The main problem is that the loss function (1-AUC) (see Equation (9)) cannot be expressed as a sum of disagreements or errors produced by each input xi.

Following a different procedure, Joachims [2005] recently proposed a multivariate approach to solve this problem with a convex optimization problem that converges using only a few constraints.

The optimization problem is:

Despite the enormous potential number of constraints, the algorithm proposed in [Joachims, 2005] converges in polynomial time. Moreover, it only requires a small set of constraints. However, the most interesting result is that the solution w of problem (14) is also the same as that of the optimization problem (11). Additionally, the slack variables in both cases are related by

Finally, the multivariate SVM returns a function fAUC of the form

Then Platt's method can fit a sigmoid to transform the output of fAUC into a probability.

2.4 Regression is a baseline approach

Considering that probabilities are real numbers, regression algorithms must be a first attempt to learn them. For this purpose, all training examples of class -1 are labeled as 0.

In order to maintain the uniformity of approach with preceding subsections, we considered the regression based on support vectors, therefore we used the so-called Support Vector Regression (SVR). Although there are least squares SVR, we used the standard version; i.e. a learner of a function

where K is once again the rbf (3) kernel, and αi are the Lagrange multipliers of the solution to the convex optimization problem:

However, given that nothing forces fRe (18) outputs to be in [0, 1], we set the hypothesis output to 1 whenever fRe returns values above 1, and 0 for fRe values below 0. In symbols, finally we have the hypothesis

3 Experimental results

Using a collection of data sets of survival probabilities in critically ill patients, we carried out an experimental comparison of four different learning approaches. SVM followed by Platt's fit of a sigmoid: the accuracy optimizer described in subsection 2.2, which will be represented by SVM(Accu); the multivariate version, aimed at optimizing the AUC (subsection 2.3), for short SVM(AUC); and finally the regression approach, SVR (subsection 2.4). The fourth predictor used was the commercial system APACHE III; we used the customization described in [Rivera-Fern´andez et al., 1998] that was developed to improve its performance in Spain.

First of all, we have to point out that this is an unfair comparison, since APACHE III was trained with a cohort of 17440 patients from 40 different hospitals in the USA [Knaus et al., 1991]; the Spanish version used records of 10929 patients from 86 ICUs; while the available data sets in our experiments only included 2501 patients. Nevertheless, this comparison is useful to test whether or not the scores achieved by SVM methods are good enough to be considered for future learning tasks.

To estimate the performance of the algorithms described in the preceding section, we used data collected from ICUs at 10 different Spanish hospitals, 6 of which include coronary patients. It is acknowledged among the medical community that coronary diseases generally have a lower mortality risk than other critical illnesses. So from a learning perspective, it makes sense to differentiate between ICUs with and without coronary patients.

The data were organized in 13 different training sets, one for each single hospital, two collecting the data from not coronary/coronary ICUs respectively, and the last one containing all the data. Each patient in these data sets was described by the same set of variables used by APACHE III. However, given that some of these variables have discrete values, we had to transform them to be handled by SVM-based systems. Thus, we codified each discrete variable using as many new binary variables (with values 0 and 1) as the number of possible values of the original variable, setting only the variable corresponding to the discrete value actually taken by the original variable to '1'.

Performance estimations were made using a 10-fold stratified cross-validation on each of the data sets, for all the algorithms except for APACHE III; since it was already trained with a different data set, we used the available data to test its predictions. Additionally, the data was standardized according to the mean and deviation observed on each training fold.

It is important to recall that the AUC achieved by the Spanish version of APACHE III in our experiments, 82.27% (in percentage) is similar to the amount reported by RiveraFern´andez et al. [1998]: 81.82%. This fact supports the representativeness of the sample of critically ill patients considered in the experiments described here.

As usual, when dealing with SVM, the parameter setting stage is very important. To set the regularization parameter C (see optimization problems in Section 2) and the rbf kernel parameter σ (see Eq. (3)) in the three support vector based algorithms, we performed a grid search on a validation set formed by the patients at 3 hospitals: one hospital without coronary patients (8), and 2 with coronary patients (1 and 9); see Table 1. The ranges searched were the following: for C we tested values from 10-4 to 102 varying the exponent in steps of 1; for σ we tested values from 10-2 to 101 varying the exponent in steps of 0.5. We found that the most promising values were C = 101 and σ = 10-2 for SVM(Accu)

and SVR; and C = 10-1 and σ = 10-2 for multivariate SVM(AUC). It is worth noting that for SVM and SVR the parameter search was aimed at minimizing the Brier score, while for multivariate SVM it was aimed at maximizing the AUC.

Table 1 shows the results obtained (Brier score and AUC) in the experimental setting described above. Focusing on the results obtained by the three support vector algorithms, we can observe that, in general, the best performance (lowest Bs and highest AUC) is achieved by multivariate SVM(AUC). The differences are statistically significant according to a one tail t-test with threshold p < 0.01. This should not be surprising for the AUC measure, since this algorithm was specially devised to optimize such a measure. But it also outperforms SVM and SVR in terms of the Brier score, whose parameters were set to optimize this score.

Let us stress that, although the optimization problem posed to SVR is precisely the minimization of the distance between true and predicted probabilities, a large amount of data is required to tie the scores of SVM(AUC) in the Brier score. The underlying reason explaining this behavior may be that the hypothesis space used by SVR is not adequate so as to induce probability distributions from a reduced set of training data, even with an rbf kernel.

As regards the data sets used in the experiments, support vector machines yielded the worst performance on the first three data sets, i.e. the smallest. SVR performance was particularly poor on these data sets. Considering that the rows of Table 1 are in ascending order of size of the data set, the trend indicates that performance could be improved if more training cases were available. In fact, when the data set included all available patients' records, the results obtained were similar to those yielded by APACHE III (recall that it was trained with data sets that were several times bigger). On the other hand, we also observe that survival predictions seem to be slightly harder for ICUs without coronary patients (Hospitals 2, 3, 6 and 8) than for ICUs including coronary patients.

4 Conclusions

We have presented a learning method for estimating probabilities in a real world problem: the prediction of survival in critically ill patients. However, the approach is general enough to be applied to other learning tasks. The method is an alternative to the standard procedure when the learning machine is based on Support Vectors and uses Platt's method [Platt, 2000] to fit a sigmoid. Instead of using an SVM devised to optimize classification accuracy, we propose to use a learner that optimizes the Area Under the ROC Curve (AUC). This can be done using a multivariate SVM described in [Joachims, 2005].

We experimentally compared the results obtained by this method with other approaches, and with a commercial scoring system trained with thousands of cases, APACHE III [Knaus et al., 1991; Rivera-Fern´andez et al., 1998]. In the reported experiments, we used real data from 10 ICUs at hospitals in Spain that contain records from 2501 patients. The medical description of each patient includes monitoring variables, clinical analysis, and demographic and diagnostic features.

The method proposed here outperforms the standard SVM approach, especially when the available data is scarce, which is the usual situation. On the other hand, increasing the number of training examples reduces differences in performance; even between probability predictions of APACHE III and those made by the baseline method, a simple regression with the output trimmed to the interval [0, 1].