The study protocol was approved and the requirement for informed consent was waived by the institutional review boards (IRBs) of all participating hospitals. In addition, the use and management of a cloud-based data storage platform for the secondary analysis was approved by the IRB of Samsung Medical Center (SMC 2020-09-100-002). All unique identifiers were removed prior to uploading. All data in the study database were assigned a research specific serial number and deidentified to protect the confidentiality of the study patients.

This study was a nationwide, multicenter, retrospective, prognostic study conducted in South Korea. We collected data for adult patients who were confirmed to have COVID-19 via real-time polymerase chain reaction and were hospitalized at 19 main referral hospitals between January 5, 2020, and August 29, 2022 (Methods S1 in Multimedia Appendix 1). Among 9199 hospitalized patients with COVID-19, we excluded 406 patients diagnosed either more than 15 days before or more than 1 day after the hospitalization date, as well as 2848 patients who had missing data in any variable of interest (Figure 1). A total of 5945 patients (5106 nonsevere and 839 severe) remained for the analysis. The 19 hospitals were divided into 4 strata according to hospital type and location: metropolitan area general hospitals, nonmetropolitan area general hospitals, metropolitan area tertiary care hospitals, and nonmetropolitan area tertiary care hospitals. We then used a random cluster sampling method to partition hospitals in each stratum and construct development and validation cohorts. For reporting and bias-risk assessment, we adhered to the following guidelines: Guidelines for Developing Machine Learning Predictive Models in Biomedical Research [14], Transparent Reporting of a multivariable prediction model for Individual Prognosis (TRIPOD; File S1 in Multimedia Appendix 1), and Diagnosis and Prediction Model Risk of Bias Assessment Tool (PROBAST; File S2 in Multimedia Appendix 1) [15].

A set of data collection guidelines were predetermined by our clinical experts. We developed a standard data collection form and prepared cloud database storage. Adhering to the set guidelines, researchers affiliated with each participating hospital gathered patient data with 32 features from demographic, clinical, laboratory, and radiological findings within the first day of hospitalization. We specified these features based on previous prognostic models and a literature review describing common biomarkers associated with severe COVID-19 [16]. The final severity status of each patient was determined on day 15 of hospitalization. All data collected in each hospital were deidentified and uploaded onto the cloud database storage. The entire data set underwent a quality assurance process, including typo rectification, outlier handling, and double-checking with the electronic health records in each participating hospital.

We declared the COVID-19 severity for patients under one or more of the following conditions during their hospitalization: (1) mechanical ventilation required; (2) extracorporeal membrane oxygenation required; (3) admission to intensive care unit; or (4) patient’s death. This criterion aligns closely with severe status (score of 6 or higher) in the WHO Clinical Progression Scale, which is developed by reaching a consensus among a group of international medical experts [17].

Among the 32 collected features, 27 readily accessible features without missing data remained for prediction modeling (Figure 2A). In order to identify subsets of robust features against feature selection methods, we considered 6 feature engineering methods (FEMs) based on 2 ML algorithms with optimal hyperparameter tuning and 4 feature importance measures: random forest (RF)-based mean decreases in Gini impurity feature importance; RF-based permutation feature importance; RF-based Shapley values; extreme gradient boosting (XGB)-based built-in feature importance; XGB-based permutation feature importance; and XGB-based Shapley values.

All features were ranked by importance measure by each FEM. We considered various criteria (ie, top K; K=5, 6, …, 15) for high-ranking features, termed high rankers. Based on each criterion K, we generated candidate feature subsets in two steps: (1) we selected high rankers with K highest importance rankings by each FEM, and (2) we identified features that were stably selected as high rankers by at least 50% of all FEMs. This process resulted in candidate subsets of robust features. The set of all 27 features was used as the reference model to identify the performance improvements in subset models during model evaluation.

A total of 60 candidate feature subsets were used for the model development, including 59 identified subsets of robust features and the reference set of all 27 features (Figure S1 in Multimedia Appendix 1). We first fine-tuned the hyperparameters of 6 ML-based algorithms, namely, deep neural network (DNN), multivariable logistic regression, RF, XGB, gradient boosting machine, and support vector machine, by applying Bayesian optimization on the development cohort. Then, we simultaneously developed all possible 360 combinations of 6 ML-based algorithms and those 60 feature subsets and evaluated the performances in both internal and external validations (Figure 2B). The model predictive performance was evaluated using the area under the receiver operating characteristic curve (AUROC) score and other cutoff-based measures, such as sensitivity, specificity, positive and negative predicted values, positive and negative likelihood ratios, and diagnostic odds ratio.

Model development and internal validation were done with the development cohort in the following steps. First, we used 5 iterations of stratified 5-fold cross-validation to explore the internal validity of each combination of feature subsets and ML algorithms as a model development procedure. The procedures were evaluated by the mean values of performance metrics and their 95% CIs calculated from the repeated cross-validation process (Methods S2 in Multimedia Appendix 1). Second, we used the entire development cohort to construct prediction models based on each of the internally validated development procedures.

External validation was conducted for the prediction models with the validation cohort. Each prediction model was evaluated by the mean values and 95% CIs of performance measures that were calculated from 100 iterations of bootstrapped sampling (Methods S2 in Multimedia Appendix 1). The final prediction model was proposed in three steps: (1) for each ML algorithm, we selected the optimal feature subset that produced the model with best discriminative performance in both internal and external validations; (2) the 3 prediction models with the best predictivity were chosen to compare their usefulness via calibration, reclassification improvement, and decision curve analysis (DCA); (3) the DNN-based final prediction model, RIETS, was proposed by considering its discriminative ability along with clinical applicability [18].

We used graphical representations to visualize the performance of RIETS contrasted with other ML-based models and provided interpretation for the selected features (Figure 2C; Methods S3 in Multimedia Appendix 1). Receiver operating characteristic curves demonstrated the discriminative model performance. Calibration plots implicated the model’s reliability in practical settings by displaying the correlation between predicted and observed risks. DCA plots indicated the net benefit of incorporating the model in clinical decision-making by quantifying the weighted trade-off between true positive and false positive identifications [19]. Reclassification plots displayed the proportion of patients that were reclassified correctly or incorrectly by RIETS compared to other ML-based models. Lastly, the Shapley additive explanations summary plot interpreted the contributions of individual features in RIETS when classifying severe and nonsevere cases [20].

We used discriminative dimensionality reduction via learning a tree (DDRTree) to cluster and characterize patients based on the features in RIETS. DDRTree is a tree-based unsupervised learning technique that reduces multidimensional features into a 2-dimensional space to visualize patients in the form of a tree structure (see Methods S3 in Multimedia Appendix 1 for procedures). This tool is known to capture cluster information with higher accuracy compared to conventional dimensionality reduction methods [21,22].

Each patient in a tree was colored with dark red to indicate high odds for severity and light green to indicate low odds for severity (Figure 3). Then, in Figure 4, dark blue and light green colors were overlaid to represent high and low concentrations of each laboratory marker, respectively. Severity risk of each patient can be identified through Figure 3A, the risk distribution with or without a pre-existing condition can be seen in Figures 3B and 3C, and feature values associated with each patient can be inferred from Figure 4. We integrated these observations to cluster patients into subgroups and characterize each subgroup (subgroup boundaries are shown in Figure 3A).

The predominant circulating variant at the time of hospitalization was identified through viral whole genome sequencing and could differ across nations [23]. According to predominant circulating variants during the pandemic in South Korea, we segmented our study period into 3 variant-dominant periods and constructed the corresponding patient subcohorts: original Alpha-dominant period (January 5, 2020 to May 1, 2021), Delta-dominant period (May 1, 2021 to November 24, 2021), and Omicron-dominant period (November 24, 2021 to August 24, 2022) [23].

We developed modifications of RIETS to explore its prediction transportability across different variant-dominant periods. Each modified model was constructed using the variant dominant subcohorts in the development cohort. For instance, the “RIETS-All” model was based on the entire development cohort and the “RIETS-Omicron” model was based on the Omicron-dominant development cohort. We evaluated all possible combinations of modified RIETS and compared their discriminative performances among patients in the external validation cohort. Consequently, we identified the best-performing model, named “RIETS-Ensemble,” that integrates the 3 models based on original Alpha-Omicron, Delta-Omicron, and Omicron. Then, the “RIETS-Ensemble” model was contrasted to the “RIETS-All” and “RIETS-Omicron” models to visualize marginal improvements. All developed models were compared using the AUROC as a measure for discriminative performance.

Patient characteristics were summarized as median (IQR) and number (%) for continuous and categorical variables, respectively, and compared between the development and validation cohorts via absolute standardized mean difference (ASMD). The ASMD was calculated using Cohen D and H formulas for continuous and categorical variables, respectively. No considerable difference was identified with an ASMD below 0.2. For cutoff-based performance measures, Youden index was used to find an optimal threshold at which the average of sensitivity and specificity was maximized. The integrated calibration index (ICI), derived from the weighted mean difference between observed and predicted probabilities for the outcome, was used to quantify and assess calibration. ICI was preferred over other calibration metrics (eg, calibration-in-the-large and slope) due to its high stability from capturing the entire range of predicted probabilities during its computation [24]. A 2-sided P value below 0.05 was set to declare statistical significance. All statistical analyses were performed using Python (Python Software Foundation, version 3.9).

Among the 5945 hospitalized patients with COVID-19 used for the development and validation of RIETS, 4019 (67.6%) and 1926 (32.4%) were allocated into the development and validation cohorts, respectively (Table 1). The median age was higher in the development cohort than in the validation cohort (mean 60, SD 45-70 years vs mean 55, SD 35-65 years, respectively; ASMD=0.333). The proportion of male patients was similar in both the development (n=2130, 48.8%) and validation (n=757, 47.9%) cohorts. Hypertension was the most prevalent comorbidity for both the development (n=1622, 37.2%) and validation (n=492, 31.1%) cohorts. While the most frequent symptoms across both cohorts were cough (n=2623, 44.1%) and fever (n=2366, 39.8%), the rankings of observed symptoms were similar in both cohorts. All variables pertaining to vital signs and blood biochemistry showed no considerable difference between cohorts (ASMD<0.2), except for the absolute neutrophil count (ANC; ASMD=0.319).

The baseline characteristics were also compared between 839 (14.1%) patients with nonsevere COVID-19 and 5106 (85.9%) patients with severe COVID-19 (Table S1 in Multimedia Appendix 1). Patients with severe COVID-19 were older (mean 70, SD 60-75 years vs mean 55, SD 40-70 years; ASMD=0.733). A larger proportion of patients with severe COVID-19 had dyspnea (ASMD=0.902) and diabetes mellitus (DM; ASMD=0.524). Patients with severe COVID-19 were more likely to have an increased respiratory rate (RR; ASMD=0.911) and decreased saturation of peripheral oxygen (SPO2; ASMD=0.705) upon hospital admission. In addition, patients with severe COVID-19 presented with higher ANC (ASMD=0.971), higher lactate dehydrogenase (LDH; ASMD=0.726), and higher white blood cell (WBC) count (ASMD=0.693).

RIETS is a DNN-based final model with the subset of 11 features that demonstrated the highest discriminative power (AUROC=0.937, 95% CI 0.935-0.938; diagnostic odds ratio=46.14, 95% CI 43.40-48.87; specificity=0.867, 95% CI 0.865-0.869; sensitivity=0.869, 95% CI 0.864-0.875) amongst the 6 ML-based models (AUROC=0.862-0.929) (Figure 5A and Table 2). RIETS also exhibited a superior discriminative ability compared to the existing low risk of bias models (AUROC=0.60-0.80; Table S2 in Multimedia Appendix 1).

In comparison with other ML-based models, RIETS exhibited net reclassification improvement (0.54%-6.14%), especially on nonsevere cases (2.14%-6.14%) (Table S3 and Figure S2 in Multimedia Appendix 1) and had the most stable prediction during cost sensitivity learning [25] (Figure S3 in Multimedia Appendix 1). RIETS also maintained sustainable prediction transportability (AUROC=0.903, 95% CI 0.897-0.910) on the limited number of cases (n=449, 7.6%) in the Omicron-dominant period when an ensemble learning technique was applied (Figure S4, Table S4, and Methods S4 in Multimedia Appendix 1).

Moreover, a PROBAST evaluation indicated that RIETS has a low risk of bias and minimal concerns regarding applicability (Methods S2 in Multimedia Appendix 1). RIETS also attained the best calibration (ICI=0.041) among comparable ML-based models (ICI=0.052-0.071; Figure 5B). Overall, it showed a higher net clinical benefit than the “intervention for none” and “intervention for all” reference strategies in DCA (Figure 5C).

RIETS comprised 11 clinical and laboratory features: LDH, age, absolute lymphocyte counts (ALC), dyspnea, RR, DM, c-reactive protein (CRP), ANC, platelet counts (PLT), WBC, and SPO2. These features were ordered by their contribution to the severity prediction by using Shapley values (Figure 6). LDH was the highest ranked, followed by age, ALC, and dyspnea. In addition, pre-existing conditions (age, dyspnea, and DM) available at the time of admission were generally ranked higher relative to those of laboratory markers (CRP, ANC, PLT, WBC, and SPO2).

We identified 4 patient subgroups using DDRTree, a tree-based unsupervised learning technique, based on the features in RIETS: the upper-right group (URG), middle-right group (MRG), lower-right group (LRG), and lower-left group (LLG) (Figure 3A). Among the 4 subgroups, the URG comprised the largest proportion of patients at high risk for severity, followed by the MRG, LRG, and LLG. The majority of patients in the URG and MRG had dyspnea, were older than 60 years, and had elevated RR, CRP, and LDH (Figure 3B). High ANC and WBC were additionally observed in the MRG. Those in the LRG and LLG had elevated ANC, WBC, and PLT. There was a negligible variation in SPO2 or DM across the tree. Moreover, we compared the patient distribution per each variant period (original Alpha-dominant, Delta-dominant, and Omicron-dominant) and found no distinguishable pattern (Figure S5 in Multimedia Appendix 1).

We developed and validated RIETS, an ML-based prognostic model for severity among patients hospitalized with COVID-19, based on a temporally and geographically extensive cohort with heterogeneous feature distributions (Figure S6 in Multimedia Appendix 1). RIETS incorporates 11 promptly and routinely available features upon hospitalization and is intended to assist early patient triaging. RIETS provides risk estimates that indicate the odds for severity progression along with feature and patient interpretation. These outputs can support clinicians in making decisions for appropriate medical measures, such as the administration of antiviral medication, transportation to the intensive care unit, and proactive preparation of medical resources. Although several prognostic models with low risk of bias excel in analogous tasks [26-30], RIETS offers substantial improvements in several aspects owing to its discriminative power and novel interpretability (Table S2 and S5 in Multimedia Appendix 1).

According to PROBAST, RIETS can be regarded as a clinically applicable model with a low risk of bias because of its generalizability and methodologically rigorous procedure. First, RIETS can be generalized across diverse populations because it was developed and validated based on a large data set from a multicenter cohort (19 general and tertiary care hospitals) over the 3-year pandemic period (from January 2020 to August 2022). In contrast, previous prognostic models were either based on a large multicenter cohort during the early pandemic period [28,31,32] or a single center cohort covering a longer pandemic period [29,30,33]. Second, we executed a rigorous modeling procedure to establish RIETS. We exhaustively developed and simultaneously validated all possible combinations of candidate feature subsets and modeling algorithms (Figure S7 in Multimedia Appendix 1). Contrary to our study design, previous prognostic studies relied on a single feature selection approach (clinical consensus, least absolute shrinkage and selection operator regression, recursive feature elimination, and sequential forward selection) [26,28,32,34,35]. Since there is no one-size-fits-all solution in the model fine-tuning process [36], this comprehensive modeling procedure can provide engineering value in attaining optimal prediction with parsimonious feature usage.

RIETS demonstrated superior discriminative performance in contrast to previous prognostic models with a low risk of bias (RIETS: AUROC=0.937, 95% CI 0.935-0.938; previous studies: 95% CI 0.60-0.80) while maintaining comparable calibration (ICI=0.041 vs calibration-in-the-large=0.00; slope=0.96; Table S2 in Multimedia Appendix 1). It has high accuracies both in severe cases (sensitivity=0.869, 95% CI 0.864-0.875) and in nonsevere cases (specificity=0.867, 95% CI 0.865-0.869). This strength can offer considerable benefits in triaging situations because prompt treatment for critically ill patients is facilitated without the resource overutilization on less critical patients [37].

RIETS also can be broadly adaptable across health care systems. Unlike some well-established models based on advanced technology-based, expensive, and time-consuming features [30,32,38], RIETS comprises 11 readily available features obtainable from routine blood tests and patient-reported conditions at admission. Thus, it is interoperable even for health care systems in low- and middle-income countries and may offer significant operational benefits during resource allocation across the global population [12]. In addition, RIETS exhibited sustainable performance on Omicron cases, implicating its potential for transportability across new variant cases with differing virulency [39] and limited case availability (Figure S4 and Table S4 in Multimedia Appendix 1).

Lastly, RIETS offers substantial interpretability that may induce improvements in model reliability and operational workflow. To our knowledge, this is the first attempt in patient clustering and characterization amongst the COVID-19 prognostic models (Figures 3 and 4). Given that bias risks are inevitable in ML systems, the interpretability of RIETS can promote transparent feedback, mitigate those bias risks, and earn trust as a clinical decision support system [40,41]. Moreover, the patient clustering tool in RIETS provides clinicians with useful information for treatment planning and resource preparation. For instance, the graphical representations of patients can enable monitoring of the characteristics of incoming patients and facilitate the identification of representative clusters at the moment. This can be used to plan the customized patient care and to initiate the preemptive preparation of medical resources for those representative patient clusters.

This study has some limitations to be addressed. First, the study participants were patients hospitalized with COVID-19 in South Korea from January 2020 to August 2022. Hence, a further study with other ethnic and variant groups is recommended to validate the generalizability of RIETS. Second, the vaccination records were not accounted for during the analysis due to a high missing rate. Although vaccination often decreases the severity [42,43], a recent study showed that some vaccinated patients with certain chief complaints remained at high risk for severity [44]. This finding implicates that the impact of vaccination on severe case discrimination may not be large as long as the distributions of clinical signs remain similar across different variants. Lastly, the information on SARS-CoV-2 variants confirmed by viral whole genome sequencing were not available for each patient. We used variant dominant periods to define variant subcohorts while anticipating some misclassifications.

We developed and validated RIETS, an ML-based COVID-19 severity prediction system, to promote the early triaging of hospitalized patients with COVID-19. RIETS demonstrated high prediction power and considerable reliability with low bias risk. Model development and validation on a nationwide, multicenter cohort implicated its generalizability. The use of routinely collected features for model construction facilitated its adaptability. Visual interpretations of model parameters and patients improved its usability and applicability. When incorporated into routine clinical practice, we anticipate RIETS to have a direct clinical impact for enabling efficient medical resource allocation as well as proactive patient care.