{'fraction_non_alphanumeric': 0.023877518641916548, 'fraction_numerical': 0.01328732349674758, 'mean_word_length': 5.415521628498728, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A novel interpretable machine learning system to generate clinical risk scores: An application for predicting early mortality or unplanned readmission in a retrospective cohort study

Risk scores are widely used for clinical decision making and commonly generated from logistic regression models. Machine-learning-based methods may work well for identifying important predictors, but such 'black box' variable selection limits interpretability, and variable importance evaluated from a single model can be biased. We propose a robust and interpretable variable selection approach using the recently developed Shapley variable importance cloud (ShapleyVIC) that accounts for variability across models. Our approach evaluates and visualizes overall variable contributions for in-depth inference and transparent variable selection, and filters out non-significant contributors to simplify model building steps. We derive an ensemble variable ranking from variable contributions, which is easily integrated with an automated and modularized risk score generator, AutoScore, for convenient implementation. In a study of early death or unplanned readmission, ShapleyVIC selected 6 of 41 candidate variables to create a well-performing model, which had similar performance to a 16-variable model from machine-learning-based ranking.


INTRODUCTION
Scoring models are widely used in clinical settings for assessment of risks and guiding decision making. 1 For example, the LACE index 2 for predicting unplanned readmission and early death after hospital discharge has been applied and validated in various clinical and public health settings since it was developed in 2010. [3][4][5] The LACE index is appreciated for its simplicity, achieving moderate discrimination power 6 by using only four basic components: inpatient length of stay (LOS), acute admission, comorbidity, and the number of emergency department (ED) visits in past 6 months. There have been ongoing efforts to derive new readmission risk prediction models for improved performance and for specific subcohorts, which, despite the increasing availability of machine learning methods in recent years, are dominated by the traditional logistic regression approach. 7,8 Logistic regression models generate transparent prediction models that are easily converted to risk scores, e.g., by rounding the regression coefficients, 9 where model interpretability is mainly reflected by variable selection. Traditional variable selection approaches (e.g., stepwise selection and penalized likelihood) are straightforward, but may not be sufficient in identifying important predictors and controlling model complexity. 10,11 Machine learning methods are well-performing alternatives, e.g., a study used decision trees and neural networks to heuristically develop a 11-variable logistic regression model from 76 candidate variables, 12 and a recently developed automated risk score generator, the AutoScore, 13 uses the random forest (RF) to rank variables and built well-performing scoring models using less than 10 variables in clinical applications. 13,14 Such machine learning methods rank variables by their importance to best-performing model(s) trained from data, but the 'black box' nature of most machine learning models limits the interpretability of the variable selection steps. The interpretability issue may be alleviated via post-hoc explanation, e.g., the widely used Shapley additive explanations (SHAP). 15 However, another general issue in current variable importance analyses has been largely ignored: restricting the assessment to best-performing models can lead to biased perception on variable importance. 16 Since predictive performance is often not the only concern in practice, a recent study extended the investigation to a group of models that are 'good enough' to generate a variable importance cloud (VIC), which provides a comprehensive overview of how important (or unimportant) each variable could be to accurate prediction. 16 Our recent work combined VIC with the well-received Shapley values to devise a Shapley variable importance cloud (ShapleyVIC), 17 which easily integrates with current practice to provide less biased assessments of variable importance. In a recidivism risk prediction study, both ShapleyVIC and VIC found an overclaimed importance of race from single models, 16,17 but ShapleyVIC could easily compute uncertainty intervals for overall variable importance for statistical assessments. In a mortality risk prediction study using logistic regression, 17 ShapleyVIC handled strong collinearity among variables and effectively communicated findings through statistics and visualizations. These applications suggest ShapleyVIC as a suitable method for guiding variable selection steps when developing clinical scores from regression models.
In this paper, we explore the use of ShapleyVIC as an interpretable and robust variable selection approach for developing clinical scores. Specifically, we use an electronic health record dataset from ED to predict death or unplanned readmission within 30 days after hospital discharge. We interpret the importance of candidate variables with respect to a group of 'good' logistic regression models using visualizations of ShapleyVIC values, combine information across models to generate an ensemble variable ranking, and use the estimated overall importance to filter out candidate variables that will likely add noise to the scoring model. To develop scoring models from ShapleyVIC-ranked candidate variables, we take advantage of the disciplined and modularized AutoScore framework, which grows a relatively simple, interpretable model based on a ranked list of variables until the inclusion of any additional variable makes little improvement to predictive performance. 13 The current AutoScore framework ranks variables using the RF, which is a commonly used machine learning method in clinical applications that performs well and is easy to train, and its preferrable performance in identifying important variables over traditional statistical methods has been previously demonstrated. 13 Hence, in this work we include RF-based scoring models as a baseline to evaluate the ShapleyVIC-based model and to demonstrate the benefits of ShapleyVIC that are not available from existing variable importance methods. We also compare these scoring models to the LACE index to show improvements over existing models.

Study cohort
This study aimed to develop a scoring model to predict unplanned readmission or death within 30 days after hospital discharge. The data was derived from a retrospective cohort study of cases who visited the ED of Singapore General Hospital and were subsequently admitted to the hospital. The full cohort consists of 411,137 cases, where 388,576 cases were eligible. 63,938 (16.5%) eligible cases developed the outcome of interest.
As summarized in Table 1, cases with and without the outcome were significantly different (i.e., had p-value<0.05) in all 41 candidate variables except the use of ventilation.
Specifically, compared to those without the event, cases with events tended to be older and have shorter ED LOS, higher ED triage, longer inpatient LOS, more ED visits in the past 6 months and more comorbid. The training, validation and test sets consisted of 272,004 (70%), 38,857 (10%) and 77,715 (20%) cases randomly drawn from the full cohort, respectively.

RF-generated models
We first describe scoring models generated using RF-based variable ranking, which ranks the 41 variables based on their importance to the RF trained on the training set and uses it in the variable ranking module of the AutoScore framework. Figure 1

ShapleyVIC-generated model
Next, we describe the assessment of overall variable importance using ShapleyVIC values and the implications, the resulting ensemble ranking that accounts for variabilities in variable importance across models, the simplification of model building steps based on the significance of overall variable importance, and the final scoring model developed. Using the training set, we trained a logistic regression model on all 41 candidate predictors by minimizing model loss, and generated 350 nearly optimal models whose loss was no more than 5% higher than the minimum loss. We used the first 3500 cases in the validation set to evaluate the ShapleyVIC values of these 350 models, which we found enough for the algorithm to converge and generate stable values in preliminary experiments.
Unlike machine-learning-based methods (e.g., RF as described in the previous subsection) that focus on relative importance of candidate variables for ranking purpose, ShapleyVIC quantifies the extent of variable importance for more in-depth inference, and communicates the findings through different forms of visualizations to facilitate interpretation. Figure 2 visualizes the average ShapleyVIC values across 350 models and the 95% prediction intervals as estimates for overall variable importance and its variability, where non-positive values indicate unimportance and are visualized using grey bars. Figure 3 presents the distribution of ShapleyVIC values from individual models to visualize the relationship between variable importance and model performance among the nearly optimal models. While both the RF (see Figure 1) and the average ShapleyVIC values (see Figure 2) suggested the number of previous ED visits as the most important variable among the 41 candidates, the 95% PI of the average ShapleyVIC value (see Figure 2) further concluded the significance of its overall importance, and the violin plot in Figure 3 revealed its much higher importance than other variables in all nearly optimal models studied. Metastatic cancer had the second highest average ShapleyVIC value among all variables, consistent with its considerable contribution to model performance observed from the RF-based variable ranking.
Both metastatic cancer and age were important to all nearly optimal models studied (indicated by positive ShapleyVIC values), but the wider spread of ShapleyVIC values for age resulted in a wider 95% PI and hence a higher uncertainty for the average value.
Similarly, although admission type had similar overall importance (indicated by the average ShapleyVIC value) as neighboring variables (cancer and systolic blood pressure), the wider spread and long left tail for admission type resulted in a 95% PI containing zero, indicating a non-significant overall importance for this variable. Among the 41 candidate variables, 20 variables had non-significant overall importance and were excluded from subsequent model building steps.
To rank variables using ShapleyVIC values, we ranked the 41 variables within each of the 350 nearly optimal model, averaged across models to generate an ensemble ranking for all variables and subsequently focused on the 21 variables with significant overall importance. Unlike the parsimony plot from RF-based variable ranking (see Figure 1), the parsimony plot based on the ensemble ranking (see Figure 4) increased smoothly until the 6 th variable entered the model, and the increments in model performance became small afterwards. The parsimony plot suggested a feasible model, Model 2, also with 6 variables: number of ED visits in past 6 months, metastatic cancer, age, sodium, renal disease, and ED triage. Four of these six variables were in common with the final models built using RF, but now metastatic cancer was selected automatically without the need for manual inspection. As a sensitivity analysis, we also generated a parsimony plot using the ensemble ranking of all 41 variables, visualized in Supplementary Figure S1. Although some variables with nonsignificant overall importance ranked higher than those with significant overall importance (e.g., admission type ranked 14 th among all 41 variables, higher than ED LOS and the next seven variables in Figure 4), these non-significant variables had small incremental change to model performance and did not affect the development of Model 2. Hence, the exclusion of variables with non-significant overall importance simplified model building steps without negative impact on the performance of the final model. The scoring table of Model 2 (after fine-tuning) is shown in Table 4, and the AUC evaluated on the test set (0.756, 95% CI: 0.751-0.760) was comparable to that of the 16-variable Model 1B and significantly higher than the LACE index (see Table 2).

DISCUSSION
Identifying important variables to the outcome is crucial for developing wellperforming and interpretable prediction models, 18 especially when developing clinical risk scores based on the simple structure of logistic regression models. The recently developed ShapleyVIC method 17 comprehensively assesses variable importance to accurate predictions, and quantifies the variability in importance to facilitate rigorous statistical inference. In this study, we propose ShapleyVIC as a variable selection tool for developing clinical risk scores, and illustrated its application in conjunction with the modularized AutoScore framework 13 by developing a risk score for unplanned readmission or death within 30 days after hospital discharge. Based completely on the scoring-generating logistic regression models, ShapleyVIC avoids bias in importance assessments due to the choice of variable ranking methods, and presents rich statistics on variable contributions through effective visualizations to support in-depth interpretation and guide variable selection. The ShapleyVIC-generated model used six of the 41 candidate variables, which outperformed the widely used LACE index and had comparable performance as a 16-variable model developed from RF-based variable ranking. Our work makes a novel contribution to the development of clinical risk scoring systems by providing an effective and robust variable selection method that is supported by statistical assessment and of variable importance and is tailored to scoregenerating regression models.
Models that are 'good enough' can be as relevant in real-life clinical applications as the 'best-performing' model, especially with respect to practical considerations such as clinical meaningfulness and costs. 16,17 The variable importance cloud proposed by Dong and Rudin 16 was the first to formally extend variable importance to such a group of models with nearly optimal performance. ShapleyVIC creates an ensemble of variable importance measures from a group of regression models using the state-of-the-art Shapley-based interpretable machine learning method, explicitly evaluates the variability of variable importance across models to support formal assessment of overall importance, and conveys such information through effective visualizations. In this work, we further propose a disciplined approach that takes into account such rich information on variable importance when ranking variables, and demonstrated its easy integration with the existing AutoScore framework 13 for automated development of scoring systems. In addition to using the average ShapleyVIC values across models as a measure of overall variable importance, we propose to use them as screening tools to exclude candidate variables that are less useful.
Although machine learning methods (e.g., RF and neural network) can be more efficient than traditional variable selection methods, 12 Four of the six variables used in Model 2 (i.e., number of ED visits in past 6 months, metastatic cancer, age, and sodium) were also used in the 16-variable Model 1B. Both models outperformed the LACE index by using fewer variables, as the CCI used in LACE requires information on 17 comorbidities. A notable difference is the role of inpatient LOS, which is included in Model 1B and LACE index but not in Model 2. The two variables used in Model 2 but not the other models, i.e., renal disease and ED triage, were ranked 22 nd and 23 rd by RF, respectively. Creatinine, which was used in Model 1B but not Model 2, is reflective of renal function to some extent. 20 The use of ED triage in Model 2 is reasonable in our setting, where the full cohort consists of emergency admissions (i.e., inpatient admissions after ED visits), and ED triage is an important predictor for short-term clinical outcomes such as 30-day mortality or readmission. 21 However, scoring models discussed above were developed for illustrative purpose and should not be use in clinical applications without further investigations and validation.
Interestingly, when using the same variables, RF had a slightly lower AUC than Model 1B (AUC=0.754, 95% CI: 0.749-0.758) and a much lower AUC than Model 2 (AUC=0.663, 95% CI: 0.659-0.668). This again highlights the drastic difference between RF and logistic regression models and their perception on variable importance, and that complex and robust machine learning models do not necessarily outperform simple scoring models.
However, compared to variable ranking based on a single RF, ShapleyVIC requires much longer run-time and larger memory space, which would benefit from the use of highperformance computers and parallel computing (option available in the ShapleyVIC R package 22 ). Another limitation of ShapleyVIC is that the sampled set of models generated using our pragmatic sampling approach may not fully represent the entire set of near-optimal models (which is formally referred to as the Rashomon set), 17 which remains a challenge in this area of research. 1,16 In addition, although collinearity was not present in our example (where generalized VIF for all 41 candidate variables were below 2), it is generally an unresolved difficulty in variable importance assessments. 16,17,[23][24][25] Our pragmatic solution of using absolute SAGE values as model reliance measures for variables with VIF>2 worked well in a previous empirical experiment, 17 and future work aims to devise more formal solutions.
Our work contributes to the recent emphasis on interpretability and transparency of prediction models for high-stakes decision making, 26 by devising a robust and interpretable approach for developing clinical scores. Unlike existing statistical or machine learning methods for variable selection, which focus on ranking variables by their importance to the prediction, ShapleyVIC rigorously evaluates variable contributions to the scoring model and visualizes the findings to enable a robust and fully transparent variable ranking procedure.
The ShapleyVIC-based variable ranking, which is tailored to the score-generating logistic regression, is easily applied to the AutoScore framework for generating sparse scoring models with the flexibility to adjust for expert opinions and practical conventions.
ShapleyVIC and AutoScore combines into an integrated approach for future interpretable machine learning research in clinical applications, which provides a disciplined solution to detailed assessment of variable importance and transparent development of clinical risk scores, and is easily implemented by executing a few commands from the respective R packages. 22,27 Although we have illustrated our approach for predicting early death or unplanned readmission, it is not limited to any specific clinical application. Moreover, the model-agnostic property of ShapleyVIC 17 and the recent extension of the AutoScore framework to survival outcomes 28 enables the use of our approach for a wider range of applications.

AutoScore framework for developing risk scoring models
The AutoScore framework provides an automated procedure for developing risk scoring models for binary outcomes. Interested readers may refer to the original paper 13 and the R package documentation 27 for detailed description of AutoScore methodology, and to a clinical application for an example use case. 14 In brief, AutoScore divides data into training, validation and test sets (typically consisting of 70%, 10% and 20% of total sample size, respectively), and ranks all candidate variables based on their importance to a RF trained on the training set. Continuous variables are then automatically categorized (using percentiles or k-means clustering) 13 for preferrable clinical interpretation and to account for potential nonlinear relationships between predictors and outcomes. AutoScore creates scoring models

Shapley variable importance cloud (ShapleyVIC)
In practical prediction tasks, researchers are often willing to select a simpler and interpretable model even when its performance is marginally lower than a more complex model. Hence, conventional variable importance assessments based on a single best model can be limiting. ShapleyVIC 17 is suitable for such purposes and assesses the overall contribution of variables by studying a group of models with near-optimal performance. In the case of developing scoring models for binary outcomes, near-optimal models can be characterized by vectors of regression coefficients corresponding to logistic loss no more than 5% higher than the minimum logistic loss.
ShapleyVIC consists of three general steps. First, a suitable number (e.g., 350) of near-optimal models are sampled from a multivariable normal centered by the regression coefficients of the optimal logistic regression model. The sampling steps are implemented in the ShapleyVIC package 22 and described in detail in the original paper. 17 Next, ShapleyVIC measures reliance of each nearly optimal model on variables using the Shapley additive global importance (SAGE) method, which measures variable impact on a model using Shapley values and is closely related to the commonly used SHAP method. 24 As explained in the original paper on ShapleyVIC, we use the absolute SAGE values instead of the unadjusted values to measure model reliance for variables involved in collinearity, identified by variable inflation factor above two from the optimal logistic regression model. Finally, the ShapleyVIC values are pooled across all nearly optimal models using a random-effect metaanalysis approach, which quantifies the overall importance of variables and explicitly evaluates the variability of these measures across well-performing models (in terms of a 95% PI for a new model) to support formal statistical inference. The average ShapleyVIC values and 95% PIs are visualized using a bar plot, and the variability of ShapleyVIC values across models is highlighted using a colored violin plot for additional insights.

Ensemble variable ranking from ShapleyVIC
To apply ShapleyVIC in practical development of scoring models based on logistic regression models, we use an ensemble variable ranking approach from ShapleyVIC that is tailored to score-generating models and filter out variables that are not likely to make significant contributions. First, we describe our proposed ensemble variable ranking for variables for a single nearly optimal model . Let ̂( ) denote the ShapleyVIC value of the -th variable ( = 1, … , ) for model , which, as described above, is either the SAGE value or the absolute SAGE value. The variability of ̂( ) can therefore be measured using the standard error of the SAGE value, denoted by ( ), that is readily available from the SAGE algorithm. Assuming independent normal distributions for SAGE values, 17 Instead of considering all candidate variables in subsequent model building steps, we propose to filter out variables that are not likely important contributors to accurate predictions based on their overall importance, corresponding to variables where the 95% PI for the average ShapleyVIC value contains or is entirely below zero. In this work, we applied the proposed ensemble variable ranking to the AutoScore framework, using it to replace RF when ranking variables.

Study design and data preparation
Our empirical experiment was based on a retrospective cohort study of patients who

DATA AVAILABILITY
Data are not available because they contain sensitive patient information from electronic health records. De-identified data is available from the corresponding author upon reasonable request.

AUTHOR CONTRIBUTIONS
NL and YN conceptualized the study and designed the analytical method. YN developed the software package for the method, conducted data analyses and drafted the manuscript. SL participated in the development of software package and data analyses. All authors interpreted the data, discussed the results, and critically revised the manuscript for intellectual content. NL supervised the study.
indicates an interval inclusive of the lower limit and exclusive of the upper limit. "--" indicates variables not included in a model. h: hours; min: minutes; ED: Emergency department; LOS: length of hospital stay.   Number of ED visit is within 6 months before current inpatient stay. Number of surgery, ICU stay and HD stay are within 1 year before current inpatient stay. CHF: Congestive heart failure; DBP: diastolic blood pressure; ED: emergency department; HD: high dependency ward; ICU: intensive care unit; LOS: length of stay; MI: Myocardial infarction; PVD: Peripheral vascular disease; PUD: Peptic ulcer disease; SBP: systolic blood pressure; SpO2: blood oxygen saturation.