Interpretable machine learning prediction of all-cause mortality

Tree-based models make minimal assumptions about the data distribution

Several assumptions associated with linear models (e.g., linearity, independence, normality, etc.) restrict the features linear models can use. To satisfy these assumptions, scientists often manually transform non-linear variables before fitting a model (e.g., log-transformation, discretization of continuous variables, etc.). For instance, to explore the effect of blood lead on mortality, researchers first discretized blood lead using different thresholds. They observed that individuals with blood lead levels higher than the threshold had increased mortality risk compared to those with lower blood lead levels [30, 34, 47]. In comparison, tree-based models make minimal assumptions about the data distribution and need no data transformations. Figure 3D shows a positive relationship between blood lead and 5-year mortality risk. Tree-based models can capture complex relationships directly without the need of manually transforming the variables.

Tree-based models capture non-linear relationships and important thresholds

Discovering non-linear relationships is important but challenging for epidemiological research using traditional linear models. J-shaped and U-shaped associations are two common and meaningful non-linear relationships [32]. However, linear models must use manually transformed features to capture non-linear relationships. As an example, Suliman et al. showed a J-shaped relationship between uric acid levels and mortality in patients with stage 5 chronic kidney disease (CKD) using a linear model by dividing uric acid level into three categories and calculating the hazard ratio for each. Unlike linear models, tree-based approaches can directly capture non-linear relationships. We observe a U-shaped relationship between uric acid level and all-cause 5-year mortality predictions in Figure 3B. This relationship differs from the J-shaped one in previous work, possibly because of categorization, which loses essential information about values within the categories.

Additionally, discovering thresholds (i.e., inflection points beyond which changing a feature’s value has diminishing returns) is significant in epidemiological analysis. Figure 3C shows that 250 μg/mL is an important threshold: according to our model, increasing urine albumin generally increases 5-year mortality risk; however, urine albumin higher than this threshold has almost the same impact on mortality risk.

Tree-based models measure feature interaction effects

Feature interaction examines how the effect of one feature on the outcome differs across strata of another feature and shows the complex relationship of two features on the outcome [7]. Tree-based models can naturally capture interaction effects by splitting on different features in the same tree. As shown in Figure 3D-F, SHAP dependence plots can be decomposed into main effects and interaction effects for each sample. Figure 3F highlights a specific interaction: the relationship of blood lead level to mortality presents differently for young and old individuals. Specifically, for those with blood lead higher than 0.1 μmol/L, younger individuals have a higher 5-year mortality risk than older individuals. Figure 3G shows the SHAP interaction effects of gender with blood lead level: females have a higher 5-year mortality risk than males with blood lead levels higher than 0.24 μmol/L. The interaction effects of age and gender with blood lead level cannot be clearly identified without SHAP interaction values because being male or older generally increases mortality risk. These findings highlight how the interaction effects detected by our model open opportunities for further research.

IMPACT verifies well-studied features associated with mortality

Some of the top 20 most important features for our 5-year mortality prediction models have been previously identified. For example, red cell distribution width (RDW), the second most important feature of the 5-year IMPACT model, has been shown to have a strong positive relationship with mortality by many studies under several conditions [10, 39, 40, 41]. We also observe the positive relationship between RDW and risk of mortality (Figure 4B); moreover, 12.7% is an important threshold over which RDW manifests a positive effect on mortality. Serum albumin level’s relation to mortality is also well-studied. Previous studies show that serum albumin is negatively associated with mortality risk [5, 13, 42]. The relationship shown in Figure 4C matches this trend. Furthermore, Corti et al. showed that serum albumin level<35 g/L was associated with a significantly increased risk of mortality compared to serum albumin levels greater than 43 g/L [5]. We observe that 35 g/L and 43 g/L are indeed key inflection points (Figure 4C): serum albumin levels lower than 43 g/L have a positive relationship with mortality prediction, while those around 35 g/L are associated with a dramatically increased mortality risk.

IMPACT identifies less well-studied features associated with mortality

Some of the top 20 most important features identified by IMPACT are less appreciated as mortality risk factors in the existing epidemiological literature. Three of these are arm circumference, platelet count, and serum chloride level. Figure 4D shows a negative relationship between arm circumference and 5-year mortality, especially for older people. This negative relationship is consistent with previous work [2, 64]. IMPACT ranks arm circumference as the fourth most important feature for 5-year mortality prediction, with an importance ranking that significantly exceeds that of BMI (the 56th). This suggests that smaller arm circumference is more predictive than BMI for modeling mortality, as in [52].

Figure 4E shows a negative relationship between platelet count, the 13th most important feature, and 5-year mortality. 175 × 1, 000 cells/μL is an important threshold; platelet count lower than that level is associated with dramatically increased mortality risk. Serum chloride is also inversely related to 5-year mortality (Figure 4F). The normal adult value for chloride is 98-106 mmol/L. We observe that serum chloride lower than 98 mmol/L is associated with sharply increased mortality risk. In Figure 4G–H, we plot the interaction effect of age and sex with serum chloride level. This analysis reveals that younger people and females with low serum chloride have a higher mortality risk than older people and males. The interaction effect of age and serum chloride shows that early rather than late-onset low chloride level has a greater effect on the model.

IMPACT can purovide additional perspective to laboratory reference intervals

Reference interval (RI) is the range of values that is deemed normal for a physiologic measurement in healthy persons [22] It is the most common decision support tool to interpret patient laboratory test results. RIs enable differentiation of healthy and unhealthy individuals [38, 21]. Hence, the quality of the RIs is as crucial as the quality of the result itself. RIs in use today are most commonly defined as the central 95% of laboratory test results in a reference population. Unfortunately, this definition does not consider mortality or disease risk, which may lead to misdiagnosis since RIs are often used to identify unhealthy individuals. The partial dependence plots of IMPACT models directly reflect the effects of the features on mortality risk, which provides an alternative perspective for identifying inappropriate reference intervals with mortality/disease relevance.

We define the relative risk percentage (RRP) that measures the relative risk (Method 1 3.3) of the feature values within the reference interval compared to the relative risk of all values (Table 2). A higher RRP indicates that the feature values within the reference interval may lead to high mortality risk, which we need to pay special attention to. The first four features in Table 2 have relatively low 5-year mortality RRP. From Figure 5A-D, we observe that the values of these features within the reference interval have low 5-year relative mortality risk; the values outside reference interval may lead to increased 5-year mortality risk. Therefore, IMPACT confirms the reference intervals of these four features as optimal for mortality risk. In contrast, the RRP of the last four features in Table 2 are high. Figure 5E-H also shows that the relative 5-year mortality risk of the values within the reference interval is high compared to the maximum relative risk of all values. Hence, IMPACT identified the divergence where reference intervals appear to be poorly tuned to mortality risk, suggesting that these reference intervals may in fact be sub-optimal for health.

• Download figure
• Open in new tab

Figure 5: Effect of varying laboratory feature values on 5-year mortality risk.

These partial dependence plots show the change in relative 5-year mortality risk (Method 1 3.3) for all values of a given laboratory feature. The grey histograms on each plot show the distribution of values for that feature in the test set. The green shaded region shows the reference interval of each feature.

IMPACT identifies trends for 1-year, 3-year, 5-year and 10-year mortality prediction models

Figure 6A shows the top 20 most important features and relative importance of input features in IMPACT’s 1-year, 3-year, 5-year, and 10-year mortality prediction models. Feature importance rankings change significantly between these four models. Some features are important for all four (e.g., age, RDW, and urine albumin level). Some features become more important over time (e.g., platelet count, whose importance ranking is 75 for the 1-year model and 12 for the 10-year model). Other features become less important over time (e.g., serum potassium, whose importance ranking is 17 for the 1-year model and 42 for the 10-year model). These results provide a more comprehensive understanding of shorter- and longer-term mortality risk, which can facilitate the investigation of mechanisms underlying risk predictors and potentially help validate interventions.

• Download figure
• Open in new tab

Figure 6: Identifying important discoveries for mortality prediction from tree-based models using different follow-up times.

(A) Relative importance of input features in 1-, 3-, 5- and 10-year mortality models. For each model, the figure shows the 20 most important features of prediction (ordered by the importance). The purple line indicates that the feature is in the top 20 features of two models. Blue and red lines indicate the feature is in the top 20 features of one model, but not in the top 20 features of the other. (B) The SHAP value of serum potassium in the 1-year mortality model. (C) The SHAP value of serum potassium in the 5-year mortality model. (D) The SHAP value of serum sodium in the 1-year mortality model. (D) The SHAP value of serum sodium in the 5-year mortality model.

The relationship between each feature and mortality may change for models that predict different mortality outcomes or utilize different subsamples of the general population. For instance, Figure 6B-C shows the SHAP value for serum potassium in IMPACT’s 1-year and 5-year mortality prediction models. The finding that serum potassium lower than 3.5 mmol/L and higher than 4.0 mmol/L are associated with increased mortality risk has been previously observed [1, 14, 36]. Interestingly, for the 1-year model, hyperkalemia (high potassium) has a higher mortality risk than hypokalemia (low potassium). For the 5-year model, hypokalemia has the same or higher mortality risk than hyperkalemia. Figure 6D shows that serum sodium higher than 139 mmol/L increases 1-year mortality risk, and low serum sodium with negative SHAP values decreased mortality risk. However, the relationship differs completely in the 5-year mortality prediction model (Figure 6E): hyponatremia (serum sodium <135 mmol/L) is associated with a higher 5-year mortality risk. This type of insight, especially regarding the differences of non-linear trends, is not apparent using linear models.

Likewise, we can compare models trained on distinct subpopulations (e.g., samples in different age groups). The differences between these models can help researchers identify risk predictors relevant to each subpopulation. Comparing models in this way can provide epidemiological insights that may guide policy for specific at-risk populations. The discoveries for mortality prediction using different age groups are discussed in Supplementary appendix 1.

Supervised distances measures feature redundancy and identifies redundant groups of features

Researchers often use unsupervised methods such as some form of correlation-based clustering to find dependent features [61, 54]. However, when we have a specific prediction task in mind, we would like to measure the feature redundancy with respect to the task. The supervised distance can be an accurate measure of this feature redundancy. (Method 1 5.1). Specifically, supervised distance measures the similarity of the two features’ information about the prediction task. It is scaled roughly between 0 and 1, where 0 distance means the features are perfectly redundant regarding the prediction task and 1 means they are not redundant at all.

To identify groups of redundant features, we hierarchically cluster all features according to supervised distance (Supplementary Figure 3; Method 1 5.1). Redundant features that have the same information about the output group together. For example, arm circumference, the fourth most important feature of 5-year IMPACT model, is grouped with weight-related features: BMI, waist circumference, weight, etc. These weight-related features all contain similar information about 5-year mortality. To further explore the predictive ability of the features, we train models using one weight-related feature and all non-weight-related features (reducing redundancy models) and models using one weight-related feature in addition to age and gender (single feature models) (Method 1 5.2). Arm circumference is the most predictive weight-related feature across all settings (Figure 7A). These results indicate that arm circumference may be more informative than other weight-related features with respect to all-cause mortality. Another example would be the cluster that includes many blood test features (Figure 7B). Similar to arm circumference, serum albumin is the most predictive feature among these blood test features. In summary, using supervised distance, we can easily identify redundant feature groups and select the most representative feature based on predictive power. These selected features can be the strongest risk predictors because they have strong predictive power and can represent a number of features.

• Download figure
• Open in new tab

Figure 7: Exploring features redundancy using supervised distance.

(A) The feature importance ranking of the BMI-related features in original models and reducing redundancy models, and the AUC of the single feature models confounded by age and gender. (Method 1 5.2) (B) The feature importance ranking of the selected laboratory features in original models and reducing redundancy models and the AUC of the single feature models confounded by age and gender. (C) The AUC of the models using the selected feature sets and the minimum feature redundancy within the selected feature sets when running supervised distance-based feature selection. The purple dashed line shows the AUC of the model trained on age and gender. The pink dashed line indicates the feature set we select for further analysis. (D) SHAP summary plot for the gradient boosted trees trained on the selected 90 features for the 5-year mortality prediction.

Supervised distance-based less-redundant feature selection

To better address the redundancy in the dataset, we propose a recursive feature selection method based on the supervised distance to select the predictive and less-redundant feature sets (Method 1 5.3; Supplementary appendix 2). Figure 7C shows the predictive power and minimum supervised distance of subsets of features refined by the feature selection approach. We can see that as the number of features reduces, the predictive performance drops, and the feature redundancy reduces (as indicated by an increasing minimum supervised distance). The figure shows that when using 90 features, the model can achieve good predictive performance (AUROC = 0.8845) and the minimum supervised distance within the features is high (0.9301). Figure 7D shows the summary plot of the top 10 features in the 5-year mortality prediction model using the selected 90 features. Since there is less redundancy in the selected features, we mitigate the issue of redundant features splitting credit. It allows us better to explore the effect of important risk predictors on mortality. In our low redundancy model, arm circumference is selected to represent the weight-related features and still receives high importance. Furthermore, we find that “requiring special healthcare equipment”, one of the top 10 features in the model trained on all features, is removed from the feature list because it is redundant with “general health condition”. In summary, our supervised distance-based feature selection method helps remove the redundant features and select the predictive and less-redundant feature set.

IMPACT develops highly accurate and efficient 5-year mortality risk scores

The predicted probability of IMPACT models can be directly used as mortality risk scores (IMPACT risk scores). We did a temporal validation of the risk scores by assessing their performances in the samples collected in NHANES 2009-2014. For comparison, we train linear and tree-based Cox proportional hazard models widely used in previous work. (Method 1 6.1) To find less costly but nearly as accurate models, we select the features using recursive feature elimination (RFE; Method 1 6.2). Moreover, we compare IMPACT risk scores with Intermountain sex-specific risk scores [18] ¹⁰ (Method 1 6.3). The models are evaluated on different gender groups.

In Figure 8A and B, we show the AUROC of the 5-year mortality risk scores of female samples (See Supplementary Figure 4 for male results) in the test set and the temporal validation set. We can see that the IMPACT model with only 20 features obtains an AUROC of 0.8971, which is almost as same as the performance of the model using all features (AUROC = 0.9030), and using fewer than 20 features leads to a dramatic accuracy drop. Figure 8A and B also show that IMPACT models achieve better performance than linear and tree-based Cox proportional hazard models. Furthermore, we can see that IMPACT risk score using the laboratory features (AUROC = 0.8881) and the risk score using the questionnaire and examination features (AUROC = 0.8835) both get acceptable predictive performance. The IMPACT risk score using the features from common test panels can achieve higher AUROCs than the intermountain risk score, which uses CBC and BMP panels features. With the models trained using different numbers of features, users can measure more risk predictors to use more accurate mortality prediction models. Figure 8B shows that the performance of our models only drops a little on the temporal validation set, which can indicate that our risk scores have generalization ability to some extent. The selected top 20 features and features from CBC and BMP panels are listed in Supplementary Table 2. In summary, we build IMPACT risk scores that are applicable to professional and non-professional individuals with different cost-vs-accuracy tradeoffs.

• Download figure
• Open in new tab

Figure 8: Developing highly accurate and efficient interpretable 5-year mortality risk scores.

(A)–(B) The AUC of the models using different feature sets after recursive feature elimination. Lines are mean performance over 1000 random train/test splits, and shaded bands are 95 percent normal confidence intervals. (A) The AUC of the models tested on the female group in the test set of NHANES 1999-2008. (B) The AUC of the models testing on the female group in the temporal validation set (NHANES 2009-2014). (C)–(D) IMPACT can analyze individualized mortality risk scores. (C) The individualized explanation for an individual who is alive after 5 years. The output value is the risk score for that individual. The base value is the mean risk score, i.e., the score that would be predicted if we did not know any features for the current output. The features in red increase mortality risk, and those in blue decrease it. (D) The individualized explanation for a sample who is deceased after 5 years.

IMPACT exposes individualized mortality risk score explanations

TreeExplainer can help researchers analyze the prediction for each individual and illustrate each features’ contribution to the mortality risk score. We explain the mortality prediction model in terms of its probability predictions (risk scores) Method 1 6.1. Figure 8C,D shows individualized explanations for two individuals from the model using the top 20 features. The first individual (Figure 8C) was alive after 5 years. From the figure, we observe that IMPACT predicted that the individual’s 5-year mortality risk score was 0.02, lower than the average predicted risk (i.e., base value). There are features that increase mortality risk, such as red cell distribution width, and features that decrease mortality risk, such as urine albumin level. For this individual, the features that drive down mortality risk outweigh those that increase it. The second individual (Figure 8D) was deceased after 5 years, and the model’s predicted mortality probability is 0.61, much higher than the average predicted risk. The top three features that increase this individual’s risk are high age, high red cell distribution width, and high urine albumin concentration. These individualized explanations can help individuals understand their health status, adjust their lifestyle, and help doctors give personalized treatment and implement precision medicines.

Method 1 3.1 SHAP (SHapley Additive exPlanation) values

SHAP (SHapley Additive exPlanation) values attribute to each feature the change in the expected model prediction when conditioning on that feature. The change of the model’s prediction when the feature is masked is recorded across all possible subsets of features, yielding an average change in prediction resulting from the inclusion of a feature in the model: where ϕ_i is the feature attribution (SHAP value) of feature i in model f for data point x, ℛ is the set of all feature permutations, is the set of all features before i in the ordering R, M is the number of input features, and f_x is an estimate of the conditional expectation of the model’s prediction: f_x(S) ≈ E[f(x) | x_S] where x_S is the set of observed features.

SHAP values which guarantee a set of desirable theoretical properties, including additivity and consistency. Additivity states that when approximating the original model f for a specific input x, the SHAP values sum up to the output f(x):

The sum of feature attributions (SHAP values) matches the original model output f(x), where ϕ₀(f) = E[f(z)] = f_x(Ø). Consistency states that if a model changes so that some feature’s contribution increases or stays the same regardless of the other inputs, that input’s attribution should not decrease. Therefore, SHAP values are consistent and accurate calculations of each feature’s contribution to the model’s prediction.

Method 1 3.2 SHAP interaction values and main effects

The SHAP interaction effects is based on the Shapley interaction index from game theory. While standard feature attribution results in a vector of values, one for each feature, attributions based on the Shapley interaction index result in a matrix of feature attributions. The main effects are on the diagonal and the interaction effects on the off-diagonal. The SHAP interaction values are defined as: when i ≠ j, and where ℳ is the set of all M input features. In Equation 3 the SHAP interaction value between feature i and feature j is split equally between each feature so Φ_i,f (f, x) = Φ_j,i(f, x) and the total interaction effect is Φ_i,f (f, x) + Φ_j,i(f, x).

The main effects for a prediction can then be defined as the difference between the SHAP values and the off-diagonal SHAP interaction values for a feature:

Method 1 3.3 Partial dependence plots and additional perspective to reference interval

We use the partial dependence plots to show the change in mortality risk for all values of a laboratory feature. The partial dependence plot shows the marginal effect one features have on the predicted outcome of a machine learning model. The relative mortality risk is defined as the average value of the model predicted probability when we fix a specific feature to a given value divided by the average value of the model predicted probability. The relative risk percentage is the maximum relative risk for the values within the reference interval divided by the maximum relative risk for all values of a laboratory feature. High relative risk percentage indicates that the values within the reference interval have a relatively high mortality risk. The partial dependence plots of selected laboratory feature values on 1-, 3-, and 10-year mortality risk are shown in Supplementary Figure 5.

SHAP value, SHAP main effect value and SHAP interaction value plots

In SHAP value/SHAP main effect value/SHAP interaction value plots, every point corresponds to a single sample where the x-axis is the value of the feature and the y-axis is the SHAP value/SHAP main effect value/SHAP interaction value. The coloring of the points often denotes the value of a separate feature.

Summary plot

Summary plots show the feature attributions (SHAP values) for many samples and multiple features in order of global feature importance (the mean absolute SHAP values). Summary plots stack multiple subplots plots for each feature. For the feature plots, every point corresponds to a single sample where the x-axis is the feature attribution value and the y-axis is vertical dispersion representing the frequency of samples with a particular feature attribution value. Finally, the color of each point represents the normalized feature value, with red representing a high value and blue representing a low one. Intermediary feature values are interpolations between red and blue.

Individualized explanation plot

Individualized explanation plot shows the feature attributions (SHAP values) for an individual in terms of how they drive the model’s prediction for the individual away from the average model prediction across the baseline distribution. The width of the bars indicate the SHAP value with red indicating a positive affect and blue indicating a negative one. The features corresponding to the largest bars are below with their actual values for the individual.

Method 1 5.1 Supervised distance and hierarchical clustering

Supervised distance can accurately measure feature redundancy based on a specific prediction task. As Supplementary Figure 2 shows, to calculate the supervised distance between feature i and feature j, we firstly train a uni-variate GBTs model to predict the label (e.g. 5-year mortality in our study) using feature i. Then, we can obtain the Prediction_i which is the output of the fitted uni-variate GBTs. Next, we fit another uni-variate GBTs to predict Prediction_i using feature j. We define the output of the new GBTs as . All hyperparameter values of the uni-variate GBTs are set to their default values. Following the same above steps, we can obtain . The supervised distance between feature i and feature j (supervised distance(i,j)) is defined as: where var(x) is the variance of the vector x, mean(x) is the average of the vector x. Supervised distance is scaled roughly between 0 and 1, where 0 distance means the features perfectly redundant and 1 means they are completely independent.

To explore the redundant feature groups, we hierarchically cluster all features according to the supervised distance. Specifically, we use complete linkage hierarchical clustering which merges in each step the two clusters whose merger has the smallest diameter. The hierarchical clustering tree is shown in Supplementary Figure 3.

Method 1 5.2 Redundant feature groups experiments training details

Method 1 6.1 Mortality risk scores training details

IMPACT mortality risk scores is defined as the prediction of the 5-year mortality prediction models. For comparison, we train linear ¹³ and gradient boosted tree-based Cox proportional hazard models ¹⁴. We did a temporal validation of the risk scores by assessing their performances in the samples collected in 2009-2014 (N = 7, 034). Specifically, the samples collected in 1999-2008 (N = 28, 820) are randomly divided into training (80%) and testing (20%) sets. To compare with Intermountain gender-specific risk scores, we evaluate the models on different gender groups. The models are trained on the whole training set and evaluate on different gender groups in the testing set. Furthermore, considering the different feature collection cost for general public and medical professionals, we build the risk scores using different feature sets. For general public, the models are trained on all demographics, questionnaire features and examination features that are accessible at home for general public, For medical professionals, the models are trained on all demographics and laboratory features. All trained models are evaluated on different gender groups of the samples collected in 2009-2014 for temporal validation.

The hyperparameters are chosen by GridSearch and 5-fold cross validation. For XGBoost 5-year mortality prediction models, the max depth is selected from {1, 3, 5, 7, 9} and the subsample ratio is selected from {0.2, 0.5, 0.8, 1.0}. The max number of trees is set to 1000. We use 20% of the training samples as validation set for early stopping. The number of early stopping rounds is set to 100. For linear Cox proportional hazard models, the regularization parameter α is selected from {0.01, 0.1, 1, 10, 100}. For tree-based Cox proportional hazard models, the max depth is selected from {1, 3, 5, 7, 9} and the subsample ratio is selected from {0.2, 0.5, 0.8, 1.0}. Other hyperparameter values are left at their default values.

We explain the mortality prediction model in terms of its probability predictions. Specifically, we rescales the SHAP values (in the log-odds space) to be in the probability space directly. The rescaled SHAP values now sum to the probability output of the model.

Method 1 6.2 Recursive feature elimination

Recursive feature elimination works by searching for a subset of features by starting with all features in the training dataset and successfully removing features until the desired number of feature remain. Firstly, we train a model on the full dataset with all features. Then we rank features by importance (mean absolute SHAP values) and remove the least important features. Another model is trained on the resulting feature set, and the process iterates until only the desired number of features are left. Starting from 151 features, we remove 6 features at the first iteration Then, We remove 5 features in each iteration until only one feature left. We bootstrap the test set and assess the predictive performance. Specifically, we resample with replacement from the test set 1,000 times and report the average and the 95% confidence interval of the AUCs. The selected features in each iteration are listed in Supplementary appendix 2.

Method 1 6.3 Intermountain mortality risk scores and exhaustive feature selection

Intermountain mortality risk scores [18] are built using complete blood count and basic metabolic profile. Specifically, 13 laboratory features are used to predict 30 days, 1-year and 5-year mortality. Logistic regression was used to model the risk prediction equations with adjustment for age and sex. Dummy variables modeled each category, with the referent defined as the lowest risk group (except for age categories: 18-29, 30-39, 40-49 [referent], 50-59, 60-69, 70-79, and ≥80 years). A scalar score value was derived for each variable category by multiplying its β-coefficient by 3 and rounding to the nearest integer (referent value = zero). Each individual’s risk score became the sum of the score values based on his or her individual data.

We implement exhaustive search to select features of Intermountain risk scores. The number of features ranges from 1 to 14 (including age). Given the number of features, we search all possible feature combinations. The risk score becomes the sum of the score values of the selected features. The 5-year mortality risk scores are evaluated on the training set. We select the feature combination that achieve the highest AUC on the training set. Then, the risk scores of the selected feature combinations are evaluated on the testing set with bootstrapping for 1,000 times.