In 2020, the clinical research ethics committee of HUF approved the “Database of COVID-19 Patients Treated at Hospital Universitario de Fuenlabrada (FUENCOVID)” project (approval 20/39). Given the retrospective nature of the study and the large number of patients included, obtaining individual informed consent was not feasible. The study relied solely on data previously collected during routine clinical care, and all data were fully anonymized prior to analysis. According to the policies of the ethics committee of HUF, and considering that the study involved no direct patient contact or additional risk to participants, the ethics committee determined that informed consent was not required. No patients received any kind of compensation.

This population-based, observational, and retrospective study covers 3 years of evolution since the beginning of the COVID-19 pandemic, with all the patients attending the only reference hospital in the city of Fuenlabrada (Madrid, Spain). All the processing has been performed using Python (version 3.10.0; Python Software Foundation) and scikit-learn (version 1.3.0) [9].

The objective was to create a large-scale database containing anonymized, structured, and nontraceable medical information of patients with SARS-CoV-2 infection who were treated at the hospital, whether admitted or discharged. The database included records of 11,975 patients.

Baseline epidemiological, clinical, laboratory, therapeutic, vaccine status, and radiological data were collected at admission for all patients. Additional information on clinical evolution and COVID-19 treatment during hospitalization was also recorded, amounting to over 700 variables. The Integrating Biology and Bedside data warehouse was used, enabling standardized data acquisition and incorporating features such as postdischarge disease evolution and mortality records. Deaths occurring during hospitalization or within 30 days after discharge were considered COVID-19–related.

In addition to standard preprocessing procedures, several additional steps were taken: (1) categorical features were encoded while preserving their order relation, (2) out-of-range values in features were removed, and (3) each patient’s initial date was determined as their first contact date with the hospital attention, and all subsequent dates, including vaccination dates, were recalculated relative to this initial date.

Missing values in the dataset were handled using the separate class method. The missing values of each feature are randomly distributed across the records and neither introduce any uninformed bias of the measurement or recording procedure nor provide any insight into the target variable in the prediction [10]. The datasets contained a high number of features relative to the number of instances, which could potentially introduce noise and reduce accuracy. To address this, recursive feature elimination with k-fold cross-validation was applied. This method selected the minimum subset of features that achieved the best model accuracy, avoiding overfitting by using cross-validation [11].

Two datasets were created for training and testing ML models using data from patients hospitalized from February 2020 to April 2022. The first dataset (exitus dataset) focused on predicting patient mortality caused by SARS-CoV-2 infection, while the second dataset (ICU dataset) aimed to predict admission to the ICU during the hospital stay. Both datasets experienced class imbalance, which could negatively impact model performance. Two oversampling techniques—synthetic minority oversampling technique and SMOTETomek [12]—were tested to address this situation. Both are based on generating new synthetic samples (data augmentation) of the minority class along the segment connecting actual minority class samples. SMOTETomek adds a cleaning step by removing the samples in the majority class closest to those in the minority class. Both techniques were used to evaluate whether new synthetic samples could improve the model performance. Since there are no standardized tests for assessing the utility of synthetic data, we relied on the improvement of the models’ performance metrics to evaluate their usefulness. Models’ metrics did not improve. Given the good performance of the models without synthetic data and concerns regarding the realism of the synthetic data when applied to complex, real-world medical scenarios [13], we decided not to use any data augmentation technique.

The models were validated using separate datasets: a validation dataset of hospitalized patients who were ultimately admitted from April 2022 to January 2023 and a validation dataset of nonhospitalized (discharged) patients between February 2020 and January 2023. All datasets share the same features and preprocessing procedures, including only the features known at admission.

Random forest (RF) is a supervised ML technique that uses an ensemble of decision trees for classification and regression. It combines the outputs of multiple trees to provide a probability of positive classification for a given input based on the ratio of decision trees that return a positive classification. The higher the probability, the stronger the confidence that the patient’s data effectively align with patterns associated with a positive outcome [14]. A threshold value should be provided so that a binary yes or no outcome is returned. If risk stratification is required rather than a binary outcome, the probability scores generated by the RF model can be used. RF is favored for three reasons: (1) it achieves accuracy comparable to deep learning models in applications without strong temporal or spatial correlation among data [15], (2) it offers better interpretability compared to deep learning models [16], and (3) previous studies have ranked RF among the top ML techniques for biomedical prediction models [17]. The architecture of an RF model is determined by several hyperparameters, which are empirically tuned using a 10-fold cross-validation grid search to optimize performance and minimize model size [18].

In data with normal distributions, the mean and SD are used, while for nonnormal distributions, the median and IQR are used. Due to the imbalance in the dataset, the accuracy metric is not recommended, as it focuses on improving the true positive rate of the majority class. To address this, the area under the receiver operating characteristic curve (AUC) is chosen as the metric to optimize during model training.

The SHAP technique is used to interpret the importance of features in predictions. SHAP is an agnostic technique that uses Shapley values from coalitional game theory to explain the significance of features in predicting an instance’s outcome [19]. SHAP calculates the difference in model output by including or excluding individual features. This process yields a Shapley value for each feature, providing insights into both its contribution and its impact on the model’s predictions. The SHAP method assesses the contribution of each feature value to the classification, taking into account the removal of the prevalence. This knowledge becomes crucial for clinicians to identify the most influential variables that significantly affect the prediction. It is essential to highlight that the significance of SHAP values does not rely solely on their magnitude, as lower SHAP values do not necessarily imply less importance. If the difference between the SHAP base value (prevalence) and the classification threshold is small, even a low SHAP value for a feature can contribute to a positive classification.

Furthermore, the SHAP value for a specific instance and feature can be decomposed into the interactions of the feature with all other feature values for that instance. The SHAP values of these interactions quantify the importance of combinations of feature values for classifying an instance.

The analysis of hospitalized patient records used for model training and testing finds that 2729 (54.6%) were male with an average age of 61.0 (SD 16.7) years, average BMI of 29.3 (SD 5.2), and average disease evolution time of 7.2 (SD 4.2) days. The remaining 2269 (45.4%) were female with an average age of 62.1 (SD 19.2) years, average BMI of 30.2 (SD 6.3), and average disease evolution time of 7.4 (SD 4.5) days. The median time to death since hospital admission was 11.9 (IQR 5.6-18.2) days, and the median time to ICU admission was 2.8 (IQR 0.8-4.8) days.

Figure 1 [20,21] illustrates the trends in admissions, mortality rate, and age over the study period. Data have been smoothed using a moving 14-day sliding window average to filter high-frequency stochastic behavior. Spain has experienced 7 waves of SARS-CoV-2 infections, leading to 7 hospital admissions waves in HUF [20].

The exitus dataset contains 4710 patients with 126 features known at admission per patient, with 472 (10%) patient deaths during admission or in the first 30 days after discharge. The ICU dataset has 4998 patients with 124 features per patient and 409 (8%) ICU admissions. The prevalence of COVID-19 infections for each SARS-CoV-2 main variant is presented in Table 1.

Feature selection using recursive feature elimination with k-fold cross-validation and AUC resulted in 60 selected features for both mortality and ICU admission prediction. Figure 2 illustrates how models’ AUCs remain stable, with values equal to or greater than 0.90, up to 60 features. Note that while the number of features is the same, the selected features differ depending on whether the prediction is for mortality or ICU admission. Finally, the exitus dataset was reduced to 35 features, while the ICU dataset contained 61 features.

The validation dataset (April 2022 to January 2023) contains 240 hospitalized patients with 43 deaths and 5 ICU admissions and 6737 nonhospitalized patients with 34 deaths within 30 days after discharge. Patients who died more than 30 days after discharge were eliminated from the validated dataset to avoid the class noise [22]. The datasets’ missingness is random and unrelated to the outcomes. Of all the patients predicted by our model to die, 21%-25% actually did. Conversely, among the patients predicted by our model to survive, 98%-99% did not die.

The best RF parameters are (1) 600 estimators, (2) Gini metric, (3) maximum depth of 20, and (4) 8 features per tree. These values have been obtained using Optuna, a framework for searching optimal hyperparameters using Bayesian optimization.

The evolution of SARS-CoV-2 variants with different clinical metrics raises the question of whether predictive models per variant can improve prediction performance. To investigate this, we created 5 models per dataset: 1 global model trained with all variants (global) and 4 per-variant models (Wuhan, Alpha, Delta, and Omicron). Figure 3 compares the AUC, true positive rate, and true negative rate—metrics not affected by the prevalence of the positive class—of the 5 models per variant.

No clear winner emerged among the models. The global model performs better in predicting mortality for the latest variants but performs worse in predicting ICU admission. This behavior is also observed in the Interpretation section.

Since models perform similarly, we present the detailed test metrics for the global model across all variants after evaluating thirty 10-fold cross-validations. The exitus predictor has a mean AUC of 0.92 (SD 0.01), mean best threshold of 0.14 (SD 0.03), mean sensitivity of 0.89 (SD 0.05), mean specificity of 0.82 (SD 0.04), mean positive predictive value (PPV) of 0.35 (SD 0.04; 134/439), and mean negative predictive value (NPV) of 0.98 (SD 0.06; 956/974). ICU predictor has a mean AUC of 0.89 (SD 0.01), mean best threshold of 0.14 (SD 0.05), mean sensitivity of 0.75 (SD 0.05), mean specificity of 0.88 (SD 0.05), mean PPV of 0.37 (SD 0.10; 92/223), and mean NPV of 0.98 (SD 0.01; 1242/1275). Figure 4 presents the receiver operating characteristics (ROCs) of both predictors.

Notably, both predictors exhibit very high NPVs, indicating effective classification for most (>90%) hospitalized patients. In addition, PPVs and specificities allow for the detection of patients who will not die or require ICU admission, advising clinicians to pay more attention to those patients who could do it. These results are consistent regardless of the predominant SARS-CoV-2 variant.

On average, across all model evaluations, the mortality prediction requires at least 70 trees returning a positive outcome to classify a case as positive, whereas the ICU prediction requires 94 trees to do the same. The lower number of features needed for mortality prediction (35 features) compared to ICU prediction (61 features) aligns with this finding (Figure 2).

Model validation using new records confirms the performance of the global models. With validation hospitalized dataset, the exitus global model has a sensitivity of 0.86, specificity of 0.74, PPV of 0.24 (18/79), NPV of 0.98 (160/163), and AUC of 0.85. The ICU global predictor has a sensitivity of 1.00, specificity of 0.59, PPV of 0.05 (5/106), NPV of 1.00 (0/151), and AUC of 0.82.

With validation nonhospitalized dataset, the exitus global model has sensitivity of 0.88, specificity of 0.98, PPV of 0.21 (30/113), NPV of 0.99 (4/6590), and AUC of 0.95. Especially relevant is that the model was able to predict the mortality of 30 of 34 patients who were not hospitalized.

In this work, we have developed and validated predictive models of mortality and admission to the ICU to be used by clinicians in the first patient assessment using records of nearly 11,000 patients of a population of 220,000 inhabitants in the Fuenlabrada area. We have trained the AI models with nearly 5000 patients who needed hospital admissions throughout the pandemic. Subsequently, we have validated the models using data from over 6000 patients from the city population collected in a different period. Our AI-based models enable us to predict ICU admission and mortality on a case-by-case basis, taking into account different factors such as age, medical history, analytical data, vaccination status, and the date of infection (that indirectly estimates the circulating strain of SARS-CoV-2) at the moment of diagnosis of SARS-CoV-2 infection. The model’s originality is that it can help clinicians to predict “yes or no,” whether a concrete patient will die or need ICU admission at its first evaluation at the hospital.

RF provides a probability score based on the ratio of decision trees that return a positive classification. A threshold value should be provided so that a binary yes or no outcome is returned. The ROC curves in Figure 4 illustrate how the sensitivity and specificity of the model change with the chosen threshold. Therefore, the model does not inherently produce a strict yes or no classification but rather assigns a probability that can be converted into a binary decision based on a predefined threshold. The models in this work use a threshold of 0.14 (Figure 4), determined by the elbow of the ROC curve. This point indicates the optimal trade-off between sensitivity and specificity, making it a reference for selecting the best classification threshold.

There are many pre-existing prediction models for mortality or ICU admission among patients with COVID-19, but most give probabilities and highlight the population characteristics that predict the patient’s outcome. Few of them apply to a specific concrete patient, as our model does. Likewise, our model is also a pandemic and a population model.

Some of the analytical variables we found to be associated with mortality have been described in other studies. Hence, several studies have identified age as the strongest predictor of mortality [23].

An analytical value that deserves to be mentioned is EDW, which is associated with higher mortality in our work. Our findings are consistent with those in the review of Lee et al [24], which included a meta-analysis of 14,866 patients from 10 studies and found that higher levels of EDW were associated with adverse outcomes in patients with COVID-19.

In this sense, there are works, in which red cell distribution width is associated with mortality in patients with COVID-19 who are nonanemic. These associations were independent of age, sex, race, and the Charlson index [25].

Furthermore, we observed a shift toward older patients requiring hospital admission since January 2022, likely due to more than 5,390,000 (80%) of the Madrid population being fully vaccinated, and admissions occurred mainly in older patients with more comorbidities and patients with immunosuppression. While overall admissions decreased, a higher proportion of hospitalizations occurred among individuals aged 65 years and older [26].

Several predictive models related to the severity and prognosis in SARS-CoV-2 infections have been published. Hudda et al [27] reviewed 19 studies, including 7 (37%) focused on model development, 2 (11%) involved external validation of existing models, and 10 (53%) addressed both the development and external validation of the same model.

Asteris et al [28] developed an artificial neural network model based on laboratory indices at admission, achieving high accuracy in predicting ICU hospitalization. They used data from 248 adult patients with COVID-19 for database creation, training, and model development. The 5 most important laboratory indices were the neutrophil to lymphocyte ratio, LDH, fibrinogen, albumin, and D-dimer. The best artificial neural network model achieved an accuracy of 95.97%, a precision of 90.63%, a sensitivity of 93.55%, and an F₁-score of 92.06%, as verified in the validation cohort [28].

Moslehi et al [29] implemented a generalized neural additive model to predict mortality in more than 2000 hospitalized patients. A total of 22 baseline features, including demographic information and clinical biomarkers, were collected. The model achieved an accuracy of 0.84 [29].

Yu et al [30] developed various models to predict in-hospital mortality in patients with COVID-19 during their first day of hospitalization using 138 patient records. Advanced age and the abnormal expression of certain cytokines were found to be closely associated with a worse prognosis. The accuracy of the logistic regression, RF, and decision tree models was 87%, 80%, and 86%, respectively [30].

Our mortality prediction model achieved an AUC of 0.90, with a sensitivity of 0.89, specificity of 0.76, PPV of 0.31, and NPV of 0.98. The model for predicting ICU admission yielded an AUC of 0.89, with a sensitivity of 0.74, specificity of 0.90, PPV of 0.41, and NPV of 0.97. Notably, both models exhibited high NPV, suggesting their effectiveness in identifying patients at low risk of mortality or ICU admission. This is particularly relevant, given that 4238 (90%) of hospitalized patients survived, and 4589 (92%) did not require ICU admission.

Furthermore, the model’s performance is independent of the circulating SARS-CoV-2 variant, demonstrating stability in both testing and external validation. A key advantage is its ability to handle missing values natively, allowing patient classification even when some data are unavailable.

To implement the model into a hospital triage system, an application should be developed that is accessible from the hospital’s digital interface and automatically collects the clinical, analytical, and epidemiological data defined by our model in the evaluation of the patient in the emergency department. This application would perform the calculation used by the model for its prediction, providing an automated response with PPV and NPV, for example, for patient mortality—risk of death: yes or no.

Patients with a no risk of death result could be discharged thinking about the high NPV of the model prediction. Patients with a “yes” risk of death result should be evaluated for hospital admission, or they could be discharged only if they are clinically stable with no supplemental oxygen requirements and with daily clinical follow-up.

This study stands out for its robust validation, which is uncommon among other AI-based predictive models for COVID-19. The validation cohort includes hospitalized patients, potentially infected with Omicron variants, as well as a significant number of patients who were discharged from the emergency department without hospitalization. The model successfully identified most patients who died within 30 days after hospital discharge in the validation cohort (30/34 deaths). We do not know why these patients were not admitted, as hospitalization decisions were made by their personal physicians. Most of these cases occurred during the first wave when the overwhelming health care crisis and epidemic conditions may have influenced such decisions. Nevertheless, we believe this does not impact our analysis, as the model is trained on patients’ initial clinical assessment. As hospitalization decisions are independent of the model, they do not influence its training process or predictions.

This study is a single-center retrospective study. The results obtained may be influenced by the number of patients hospitalized during each wave and could be specific to periods dominated by a particular SARS-CoV-2 variant. However, our global model has demonstrated stability in these contexts and has also been validated using data from a different period and with patients distinct from those used to create the predictive model.

Another limitation is that our model was developed using data from 5000 patients and validated with another 6000 patients, all belonging to the same region of Spain (the city of Fuenlabrada). It would be of interest to apply this model to other regions in Spain or other countries worldwide to assess its generalizability.

The scientific community has increasingly focused on integrating AI into COVID-19 predictive modeling. RF techniques were predominantly used in papers published from 2019 to 2023. However, most studies lacked external validation, which hindered reproducibility and applicability to diverse datasets. None of these studies covered the entire pandemic period, and nearly all involved fewer than 5000 patients. While the models demonstrated moderately good performance, with AUC values above 0.7, data collection often included demographics, clinical records, laboratory results, and pharmacological studies, but the applicability settings of these models were often unspecified.

Our study is unique, as it represents a population-based study in a city reflective of the Spanish population. Unlike other studies, it covers a 3-year pandemic period. We developed an AI predictive model for COVID-19–related mortality and ICU admission, which was validated with patients different from those used to create the model. This model covers a 3-year pandemic period, demonstrating its usefulness during both pre- and postvaccination periods. Ethnicity could likely enhance the model’s utility, especially in Southern Europe. We discuss ethnicity because this study has been conducted in a Mediterranean country, specifically in a city with 217,184 inhabitants, where 195,465 (90%) of the population is Spanish and 23,337 (10.7%) are immigrants (n=5991, 2.8% from Central and South America). The model, based on ML, can be automatically applied when a patient seeks medical attention. It is valuable for identifying patients requiring hospitalization and those eligible for home care. This multidisciplinary study, conducted by clinicians and engineers, used data from over 12,000 patients and extensive covariables. Our model offers practical applicability in daily clinical practice, providing a clear risk assessment for mortality and ICU admission.