Early prediction of in-hospital death of COVID-19 patients: a machine-learning model based on age, blood analyses, and chest x-ray score

An early-warning model to predict in-hospital mortality on admission of COVID-19 patients at an emergency department (ED) was developed and validated using a machine-learning model. In total, 2782 patients were enrolled between March 2020 and December 2020, including 2106 patients (first wave) and 676 patients (second wave) in the COVID-19 outbreak in Italy. The first-wave patients were divided into two groups with 1474 patients used to train the model, and 632 to validate it. The 676 patients in the second wave were used to test the model. Age, 17 blood analytes, and Brescia chest X-ray score were the variables processed using a random forests classification algorithm to build and validate the model. Receiver operating characteristic (ROC) analysis was used to assess the model performances. A web-based death-risk calculator was implemented and integrated within the Laboratory Information System of the hospital. The final score was constructed by age (the most powerful predictor), blood analytes (the strongest predictors were lactate dehydrogenase, D-dimer, neutrophil/lymphocyte ratio, C-reactive protein, lymphocyte %, ferritin std, and monocyte %), and Brescia chest X-ray score (https://bdbiomed.shinyapps.io/covid19score/). The areas under the ROC curve obtained for the three groups (training, validating, and testing) were 0.98, 0.83, and 0.78, respectively. The model predicts in-hospital mortality on the basis of data that can be obtained in a short time, directly at the ED on admission. It functions as a web-based calculator, providing a risk score which is easy to interpret. It can be used in the triage process to support the decision on patient allocation.


Introduction
Starting from late February 2020, the COVID-19 outbreak struck the north of Italy causing more than 30,000 deaths in Lombardy alone, up to the end of March 2021. At the beginning of the outbreak, the Spedali Civili di Brescia (SCBH), the university hospital of one of the hardest hit cities in Europe, was faced with a 'flash flood' of severely ill patients seeking admission to the emergency department (ED). For several weeks, their number exceeded the available resources, obliging a continuous organizational restructuring of the hospital wards (Garrafa et al., 2020b).
In those weeks, given the limited evidence of clinically proven predictors (Marengoni et al., 2021;Wynants et al., 2020;Sperrin et al., 2020), prioritizing hospital admission of non-critical patients was an arduous task. Essentially, the criteria were based on the presence of fever, respiratory symptoms, and the level of blood oxygenation. A significant drawback of this approach was that patients referring to the ED with very similar clinical findings underwent inconsistent assessments. In this scenario, the availability of predictors would have been extremely beneficial, not only to triage patients, but also to monitor hospitalized patients and warn of exacerbation of the outbreaks.
Starting from March 2020, all patients referred to EDs underwent a chest X-ray at admission or within a few hours. With the purpose of grading pulmonary involvement and tracking changes objectively over time, a chest X-ray severity score was developed (Brescia X-ray score) (Borghesi and Maroldi, 2020;Maroldi et al., 2021;Borghesi et al., 2020a;Borghesi et al., 2020b). The score was able to predict in-hospital mortality in 302 patients. In addition to the chest X-ray severity score, a dedicated blood sampling profile was included in the COVID-19 ED work-up (Garrafa et al., 2020a). Among its 17 blood analytes, the sampling profile encompassed hemachrome, inflammation biomarkers such as C-reactive protein (CRP), lactate dehydrogenase (LDH), and ferritin, and coagulation markers (fibrinogen and D-dimer). Since that time, the medical literature began to encompass an increasing number of studies advocating the prognostic value of single or grouped blood parameters (Bonetti et al., 2020;Borghi et al., 2020;Avouac et al., 2021;Knight et al., 2020). All these parameters were present in our COVID-19 sampling profile.
This study aims to develop and validate an early-warning model (BS-EWM), predictive of in-hospital death, based on data that could easily be acquired on admission to the ED: age, simple blood biomarkers, and chest X-ray. The model was constructed based on the analysis of a cohort of 2872 COVID-19 patients treated in a single reference center over a 10 -month period.
This paper adheres to the TRIPOD checklist for predictive model development and validation (Collins et al., 2015).
The study was approved by the local ethics committee (COVID-SURG-BS; NP 4059).

Results
Description of the sample The descriptive statistics for all variables in the dataset are presented in Supplementary file 1b and were computed and stratified by the two waves (MA vs. MD) and by outcome (alive vs. dead). The two subsets were similar for most variables.
The correlations between the 17 analytes and the Brescia X-ray score were investigated using Spearman correlation coefficients and visualized using a correlation plot (Figure 2). The Brescia X-ray score was positively correlated with neutrophil to lymphocyte ratio (NLR), CRP, LDH, standardized ferritin, and D-dimer, and was negatively correlated with lymphocyte %, monocyte %, and basophil %.

BS-EWM
A machine-learning model (BS-EWM) was developed by inputting a dataset of 2782 COVID-19 patients admitted to the ED and hospitalized at SCBH from March to December 2020. The majority of patients (2106/2782, 75.70%) belonged to the first wave (MA), the remaining fraction (676/2782, 24.30%) to the second wave (MD). As outcome, the machine-learning model had the condition dead/ alive, and, as covariates: age, Brescia X-ray score, and 17 blood sample analytes. Figure 2. Correlation plot on biomarkers and Brescia chest X-ray score. The relationships between 17 analytes and Brescia chest X-ray score are inspected with the Spearman correlation coefficients, ρs, which are represented in this correlation plot by means of blue and red circles (positive and negative correlation, respectively). The diameter of the circle is proportional to the magnitude of ρs and black crosses on them identify correlation not significantly different from zero (p-values > 0.05). The correlation matrix is reordered according to the hierarchical cluster analysis on the quantitative variables. Figure 1 reports the flowchart that describes how data were divided for training, validation, and testing the BS-EWM.
The synthetic minority oversampling technique (SMOTE) procedure, rebalancing the dead/alive ratio (50% vs. 50%) from the original 20.09%, improved accuracy, specificity, and sensitivity of the random forest applied on it (see Supplementary file 1c which compares performance metrics with/ without the SMOTE method).
The relative variable importance measure (rel VIM) and partial dependence plot (PDP) were extracted from the random forests (Figures 3 and 4, respectively). In Figure 3, the rel VIM of BS-EWM based on age, Brescia X-ray score, and 17 blood analytes are reported on a bar plot. Since age was strongly associated with the risk of death, it masked the role of the other covariates. For completeness, the relevance of the 17 analytes and Brescia X-ray score was estimated in an additional EWM, in which the covariate 'age' was excluded. In Figure 4, 9/17 analytes and the Brescia X-ray score were noted as being important in predicting the risk of death (rel VIM >60). The effects of changes in covariate values on the risk of death threshold of the EWM were reported by means of a PDP (a 2D plot in the x-y plane) (Figure 4). Only fibrinogen was excluded from this graphical representation since in Table 1, it was not significantly different in the two subpopulations deceased/alive. Most PDPs showed nonmonotonic increasing relationships between the x-variable and the EWM, resulting in a plateau corresponding to high values of x.
When compared to other models such as gradient boosting machine (GBM) and logistic regression, the random forest showed better performance in terms of area under the curve (AUC), sensitivity, and specificity. The in-sample sensitivity (0.93) yielded by the model was the highest, and it maintained an important 0.82 in validating the out-of-sample sensitivity, and this decreased to 0.73 when testing the MD subgroup (see Table 2 which contains details on all the metrics extracted from the ROC analysis). ROC curves are visualized in Figure 5 where, for each model (random forest, GBM, and logistic regression), the performances in training, validating, and testing are compared in a unique graph.
In order to compare the BS-EWM score with univariate models based on single biomarkers, three random forest (on training, validation, and testing) are estimated on the most important biomarkers (LDH, D-dimer, neutr/lymph, neutrophils %, fibrinogen, CRP, Brescia chest X-ray, lymphocytes %, ferritin std, monocytes %). Results of these 30 models are reported in Supplementary file 1d. It is Figure 3. Relative variable importance measure (rel VIM). In , Figure 3, A1, there is the rel VIM based on Gini index. It was extracted from a random forest where the outcome is dead/alive and covariates are: the 17 biomarkers, Brescia X-ray score, and age. The algorithm grows 10,000 trees where the number of splitting variables at each tree node is √(# covariates in the model). Missing values are imputed with the 'on-the-fly-imputation' algorithm. A model with the same features was run (Figure 3, A2) excluding the covariate 'age' since it was strongly associated with the risk of death, masking the role of remaining covariates.  . Partial dependence plot (PDP) of random forest grown on the 17 biomarkers and Brescia X-ray score. Considering the random forest that excludes the 'age' variable, the PDPs were computed in correspondence of covariates with relative variable importance measure (rel VIM) of Appendix 1-figure 2 > 60 (cut-off identified by the red dashed line) and p-value in  Comparison between the performances of three methods: random forest, GBM, and logistic regression model applied on the rebalanced dataset obtained with SMOTE methodology.
Logistic regression predictions are computed using the 10-fold cross-validation in order to be comparable with random forest and GBM predictions (which use out-of-bag and 10-fold cross-validation, respectively). evident that considering one biomarker per time, the model provides good predictions in training, but bad performances out of sample (contrary to the BS-EWM score it loses its predictive power on fresh data). Hence, it is evident that a score based on a multivariate model provide better results since it also considers interactions among variables.

Discussion
The dataset for the development, validation, and testing of the BS-EWM originated entirely from an Italian region, potentially limiting the generalizability of the risk score in other areas of the world. Additional validation studies from different geographic areas are welcomed. Furthermore, though the BS-EWM has been validated using blood sample values obtained by instruments that satisfy internal and external quality control, different equipment could lead to divergent results (Martens et al., 2021;Lippi et al., 2020). Therefore, it would be appropriate to harmonize the results. Another limit could have been the presence of missing values, though the BS-EWM has also performed adequately in this condition since it used a multiple imputation technique to overcome the problem. Finally, it is important to point out that the BS-EWM risk score should not be used for asymptomatic COVID-19 patients or for the pediatric population. It will be interesting, in the future, to verify if the BS-EWM could be applied by general practitioners to the unhospitalized population. This could allow the generalizability of the model outside the hospital context to be tested. Though the BS-EWM has been developed on a cohort of 2106 patients belonging to the COVID-19 first wave, the model also demonstrated a sensitivity greater than 70 % in the early prediction of high risk in patients in the second wave, when in-hospital mortality was 40 % lower.
While these models are mostly based on age and a set of vital (clinical) parameters, in addition to age, the BS-EWM depends on blood parameters. It is conceivable that blood analytes capture a snapshot at hospital admission signaling a specific bodily reaction to viral infection in terms of hyperinflammation, immune response, and thrombophilia. On the other hand, the other models are more influenced by the general status of the patient, which may be determined by concomitant and preexisting diseases. According to the International Federation of Clinical Chemistry (Bohn et al., 2020), no single biochemical or hematological marker is sufficiently sensitive or specific to predict the outcome of SARS-CoV-2 infection. Notably, the IFCC recommends that the interpretation of laboratory abnormalities should be based on groups of analytes (Bohn et al., 2020). In the BS-EWM, three analytes reached a significant value in predicting death: LDH, D-dimer, and NLR. LDH is a non-specific indicator of tissue damage (Bohn et al., 2020;Liang et al., 2020) that emerges as one of the most consistently elevated markers in patients at higher risk of developing adverse outcomes, probably because COVID-19 infection is characterized by systemic tissue damage. Another key feature of SARS-CoV-2 is the coagulopathy: high levels of D-dimers have been reported to correlate with unfavorable disease progression in several cohorts of patients. The coagulopathy linked to COVID-19 infection is likely to involve a complex interplay between pro-thrombotic and inflammatory factors, thus the combined analysis of both inflammatory and thrombophilic markers could play an important role in the early identification of patients at higher risk of unfavorable progression (Bohn et al., 2020;Lazzaroni et al., 2021). Finally, lymphopenia has become a hallmark of SARS-CoV-2. It has been demonstrated in almost all symptomatic patients, though in varying degrees. Disease severity has been correlated with the level of lymphocyte count reduction. A direct infection of lymphocytes, which express the coronavirus receptor ACE-2, is among the mechanisms proposed. A poor prognosis is also associated with an elevated neutrophil count combined with lymphopenia, resulting in a high NLR. The increase in granulocytes is the result of the cytokine storm induced by the virus and is responsible for tissue damage (Bonetti et al., 2020;Bohn et al., 2020).
A further remark concerning the blood analytes is that, in the BS-EWM, the thresholds of the single analytes (namely the point where the functions in Figure 4 become constant and the probability of death no longer increases/decreases) closely overlap with the values recently proposed by other authors (Webb et al., 2020;Caricchio et al., 2021). For completeness, the optimal threshold (computed through the Youden index) for each biomarker to predict the outcome (dead or alive) are reported in Supplementary file 1e.
The present study is not unique in encompassing radiological findings combined with blood analysis. The study by Schalekamp et al., 2021, integrated blood analysis parameters and radiological information derived by grading chest X-rays (0-8 scale points). Unlike the cited study, with the BS-EWM in this study, the radiological score did not reach a high relevance (rel VIM) in predicting high risk. This difference can be explained by the different approaches used to build the model (logistic regression vs. random forests) and by the high degree of correlation of the X-ray score with multiple blood analytes: 'collinearity' thus could have 'stolen' importance from the information provided by imaging. Nevertheless, at admission, the chest X-ray score of patients who subsequently died was significantly higher than for patients who survived. Furthermore, the chest X-ray score may provide additional stability to the model, playing an important role in the case of missing data in the blood sample counterpart.
Further, the BS-EWM delivers high prediction performance and only requires a limited number of readily available variables with easy operability, no time consuming, no extra money since these analytes are required for COVID-19 diagnosis and monitoring. An important and pragmatic aspect offered by the BS-EWM is that the biomarkers employed may be obtained by the emergency laboratory in less than an hour (Garrafa et al., 2020a) and, differently from other biomarkers (Kyriazopoulou et al., 2021), they are non-expensive and frequently used also in developing countries. It is important to note that the same methodology could be applied to other infections and be practical to triage people.
Most laboratories, including the small or peripheral ones, may provide results in a short time. At the Spedali Civili of Brescia, the BS-EWM is integrated within the Laboratory Information System (LIS). It works as a web-based calculator and is easy to interpret. The online calculator allows an easy assessment of the EWM, requiring only the entering of analyte values and X-ray score. The score is calculated even if some of the values are missing. Furthermore, in our center the system may be integrated with the electronic health record and the radiology information system, allowing a completely automatic data retrieval and entering, without any operator interaction. It provides a risk threshold of 0.5, above which patients are graded as having a potentially high death-risk, thus supporting closer clinical observation or admission to a high-intensive care ward. In patients yielding a low risk (score 0-0.49), the decision by clinicians to allocate them to a low-intensive care ward or to monitoring is further sustained. The online calculator allows an easy assessment of the EWM, requiring only the entering of analyte values and X-ray score. The score is calculated even if some of the values are missing. Furthermore, in our center the system may be integrated with the electronic health record and the radiology information system, allowing a completely automatic data retrieval and entering, without any operator interaction.
Finally, the need to regularly update models and closely monitor their performances over time and geographically should be underlined, given the rapidly changing nature of the disease and its management.
According to the two temporal peaks of incidence of the COVID-19 outbreak in Lombardy, the 2782 patients were divided into two groups: (i) MA including 2106 patients admitted during the first wave; (ii) MD including 676 patients in the second wave. Quantitative variables were described using mean (SD), median (IQR), and range (min-max), while categorical variables were reported as counts and percentages. The comparisons between groups were performed using the Wilcoxon rank-sum test for quantitative variables and Fisher's exact test for qualitative variables.
The relationships between the 17 analytes and the Brescia X-ray score were inspected using the Spearman correlation coefficient, ρs, and visualizing results using a correlation plot (Dancelli et al., 2013;Marziano et al., 2019;Figure 2).
To estimate the BS-EWM, the outcome (alive/dead) was modeled using as covariates: (i) Brescia X-ray score, (ii) 17 analytes, (iii) age. Since most of the covariates analyzed were strongly correlated (multi-collinearity) ( Figure 2) and their relationships with the outcome were non-linear, the BS-EWM was estimated using random forests (Breiman, 2001;Carpita and Vezzoli, 2012), a non-parametric machine-learning method (Vezzoli, 2011;Vezzoli et al., 2017). Moreover, the algorithm is able to manage missing values which are common in clinical studies. The 'on-the-fly-imputation' algorithm (Hong and Lynn, 2020) imputes data when it grows the forest handling interactions and non-linearity in the dataset.
Since the prevalence rate of death in the two waves was different (20 % in MA vs. 12 % in MD), a strategy to generalize results in unbalanced datasets was applied, adopting a rebalancing method able to improve the detection of patients with a high death-risk.
The EWM was developed using the 2106 patients in the first COVID-19 wave (MA 2020) when in-hospital death prevalence was 20 %. Seventy percent of them (1474 patients) were used for training the model and the remainder (632 patients) for testing it. Patients were randomly assigned to the two subgroups, and further stratified according to the outcome (alive/dead). Consequently, both the training and testing subgroups included the same rate of deaths (20.09%) as the full sample (2106 patients). With such a 'moderate' incidence of death, the dataset was statistically unbalanced. This limitation could have implied the development of a model yielding unsatisfactory results in predicting new observations for the minority class, that is, patients with death as outcome. An approach to address this limitation is to oversample the minority class (deceased patients) and, subsequently, create the predictive model (BS-EWM). The SMOTE (Chawla et al., 2002) was chosen. The SMOTE function oversamples the minority class by using bootstrapping and k-nearest neighbor to synthetically create additional observations belonging to that class (dead). The procedure is combined with under-sampling of the majority class (alive). To determine the optimum number of k-groups into which to assign the dataset, a matrix containing the 17 analytes and the Brescia X-ray score was used to compute the hierarchical cluster analysis (Salvi et al., 2019;Codenotti et al., 2016). By means of silhouette analysis, k = 2 was determined as the optimal number of clusters into which to assign the dataset. Hence, a synthetic rebalanced dataset was obtained with an equal number of living and deceased patients (888 + 888). The rebalancing procedure enabled a risk score to be devised ranging from 0 to 1 with a threshold of 0.5 to separate non-severely affected from severely affected patients. Subsequently, the model was tested on the subgroup of 632 patients in the first wave excluding the training set. A further validation of the EWM was conducted on the 676 COVID-19 patients in the second wave (Wynants et al., 2020).
The rel VIM (Carpita andVezzoli, 2012, Doglietto et al., 2020b) and the PDP (Friedman, 2001;Doglietto et al., 2020a) were extracted from the model for a better understanding of the relationship between outcome and covariates.
The predictions extracted from the random forests classification were interpreted as in-hospital death probability conditional on the combination of the values of analytes, Brescia X-ray score, and age in COVID-19 patients at admission to the ED.
The BS-EWM performance was evaluated by AUC of an ROC curve. The robustness of the model was compared to other models by running GBM (a machine-learning approach and competitor to random forests), and logistic regression, and computing the same metrics.
The BS-EWM score is available for use online (https://bdbiomed.shinyapps.io/covid19score). In the SCBH it is integrated within LIS returning the death-risk score directly from the medical report.
All the analyses were performed by R, version 4.0.0 (R Development Core Team, 2020). The code is available at 288 https://github.com/biostatUniBS/BS_EWS copy archived at swh:1:rev:7416ba71075402e6a0ed997e7aa6a527e93247b2 . The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Data availability
We are unable to share the dataset as it contains sensitive personal data collected during the pandemic at Spedali Civili di Brescia. We cannot share the full data since are data from patients. Interested researchers should contact the authors for any query related to data sharing and submit a project proposal Once defined the goal of the study, and the data needed authors will submit the potential project of collaboration to the IRB of Spedali Civili di Brescia to receive approval to access a deidentified dataset. Please note that other informations related to patients can be acquired, always following approval of IRB of Spedali Civili di Brescia, not only the ones studied in the paper. Anyway, following request to the authors, it will be possible to share processed version of the dataset ( e.g. an Excel sheet with numbers used to plot the graphs and charts of the manuscript). All code used to analyse the data can be found on GitHub at https://github.com/biostatUniBS/BS_EWS, (copy archived at swh:1:rev:7416ba71075402e6a0ed997e7aa6a527e93247b2).