Among hospitalized patients with COVID-19, there is significant interindividual variability. A significant number (20%-67%) progress from moderate illness to life-threatening complications, including acute respiratory distress syndrome (ARDS) [1,2] and septic shock [3], generating a surge in patients who require intensive care unit (ICU)–level respiratory and vasopressor support [4]. Among patients with COVID-19 who are critically ill and require invasive mechanical ventilation, a delay in intubation from the first noninvasive respiratory support is associated with an increase in hospital mortality [4]. Similarly, delayed vasopressor initiation in patients with septic shock has been found to be associated with increased mortality [5]. Acute kidney injury (AKI) is also common among hospitalized patients with COVID-19 and is associated with high mortality [6]. In a recent observational study of 3993 hospitalized patients with COVID-19, AKI occurred in 46% of patients, and 19% required dialysis [7]. A recent meta-analysis of 34 observational studies of hospitalized patients found that delayed ICU admission was remarkably associated with mortality, highlighting the importance of providing timely critical care in non-ICU settings [8].

To support risk stratification among heterogeneous hospitalized patients, recent studies have used machine learning–based clustering [9] to retrospectively analyze routinely available patient electronic health record (EHR) data to identify clinically useful phenotypes [10]. In critical care research, unsupervised machine learning clustering has been used to identify homogeneous subgroups within a broad heterogeneous hospitalized population [11], which elucidates pathophysiology, can predict treatment response, and has the potential to augment clinical trial enrollment [10]. The most common clustering techniques used in medicine are latent class analysis (LCA), an algorithm that derives clusters using a probabilistic model that describes the distribution of the data [12,13], and k-means, which identifies clusters in a data set by using a distance metric to find k centroids (a weighted average) within the n-dimensional space of clinical features [11-14]. Both LCA and k-means have been effectively [15] used to detect homogeneous phenotypes with distinct severities and treatment responses in ARDS [16-18], sepsis [19,20], and COVID-19 [21,22].

In support of point-of-care clinical management, modern predictive machine learning classification algorithms (eg, the gradient-boosting machine [GBM] algorithm [23]) trained using features based on observations recorded early in a new encounter have shown promise in rapidly assigning de novo patients to a clustering-identified phenotype [24]. GBM classifiers are increasingly being applied for prediction in the data science industry and are known to outperform simpler models such as logistic regression in many clinical research fields, including critical care [25,26]. GBM has been used to accurately identify LCA-derived ARDS phenotypes [24], including a hyperinflammatory phenotype characterized by elevated inflammatory biomarkers, higher prevalence of vasopressor use, longer use of ventilation, extended length of stay, higher prevalence of sepsis, and higher mortality [27-31]. In addition, a recent ARDS study observed differential responses to positive end-expiratory pressure strategy by phenotype, with higher positive end-expiratory pressure associated with improved outcomes in the hyperinflammatory phenotype [27].

A recently reported EHR data clustering analysis of a relatively small sample of patients with COVID-19 admitted to a US hospital identified 2 phenotypes designated as cluster 1 and cluster 2 [21]. Patients in cluster 1 were older individuals (mean age 79.5 years) with multiple comorbidities and a higher mortality rate (25.4% vs 8.97%; P<.001) than patients in cluster 2. Patients in cluster 2 were younger individuals (mean age 53.7 years) who were more likely to be male and racial and ethnic minority individuals with higher levels of inflammatory markers and alanine aminotransferase (ALT) and a markedly increased BMI.

In this study, we sought to explore the generalizability of this 2-phenotype finding for COVID-19 using a clustering analysis of EHR data associated with a much larger cohort of hospitalized patients. Analogous to the ARDS study cited previously, we also explored the application of GBM-based phenotype classifier algorithms trained using routinely available clinical data for the rapid identification of clustering-derived COVID-19 phenotypes.

Deidentified EHR data were extracted from the JH-CROWN Registry [32] on patients with COVID-19 who were admitted to the Johns Hopkins (JH) Health System from February 25, 2020, to March 3, 2021. The registry, constructed directly from the JH clinical EHR, was designed to serve as a comprehensive projection of structured clinical data for patients with COVID-19. Diagnosis of COVID-19 was defined as a positive molecular test for SARS-CoV-2 and either a COVID-19 International Classification of Diseases, 10th Revision, diagnosis or an associated diagnosis suggesting that COVID-19 was likely present (eg, pneumonia, ARDS, or anosmia). Patients transferred from other health care institutions were excluded. Given that our goal was to identify phenotypes potentially at high risk of deterioration, we also excluded patients who had initiated critical care treatment (eg, invasive mechanical ventilation or dialysis) or died during the 6-hour time window following admission. The extracted registry data included patient demographics, encounter information, problem lists, diagnoses, flow sheets, laboratory test results, medications, procedures, and outcomes associated with the patients (N=4379).

Data used for clustering included age, BMI, and 36 clinical observations (vitals and laboratory tests) selected based on registry data availability (<25% missingness, as shown in Figure 1) associated with the included patients with COVID-19. Figure 2 is a correlation heat map showing that our selected data elements were mostly uncorrelated except for expected strong associations in observations, such as between creatinine and blood urea nitrogen [33]; among white blood cell count, lymphocytes, and neutrophils; and among red blood cell count, hemoglobin, and hematocrit [34]. Clustering features were generated as minimums or maximums of these vitals and laboratory tests [31] within the context of severe COVID-19 illness recorded during the entire hospital stay (Textbox 1).

The practice of splitting a data set into training and validation data sets toward assessments of generalizability of machine learning–based subgroup discovery is well established [35]. Accordingly, in our study, patient encounters were split into 2 cohorts randomly. Cohort 1 (2179/4379, 49.76%) was used as the training cohort, and cohort 2 (2182/4379, 49.83%) served as the internal validation set. Table 1 shows the basic demographics, comorbidities, vitals, and inflammation biomarkers associated with the full cohort as well as the highly similar training and validation cohorts. Following data cleansing to account for data outliers, as described in more detail in the following sections, missing data imputation and clustering analysis for each cohort were independently performed. Beyond comparing overall clustering result indexes such as the number of phenotypes identified and how Textbox 1 features were statistically distributed across the phenotypes identified in the 2 data sets (ie, internal validation [35]), we also explored the similarity of clinical data distributions across phenotypes detected in the training and validation cohorts not used for clustering (ie, external validation [35]). The features used for external validation included inflammatory biomarkers not used for clustering because of excessive missingness (eg, D-dimer and ferritin), treatment response (eg, the need for critical care treatment such as invasive mechanical ventilation, dialysis, and vasopressors), or outcomes (eg, length of stay, sepsis, and survival).

A recognized challenge in the use of observational clinical data in machine learning analytics is the need to account for potential biases resulting from treatment that can influence patient physiological measurements [36]. For example, in our study, a significant number of inpatients received supplemental oxygen for acute COVID-19 respiratory symptoms in an emergency room setting before admission, thus potentially biasing observations such as oxygen saturation (SpO₂) measured following admission. Another potential source of bias can be treatment for hypotension upon presentation, which is known to occur in patients with chronic hypertension [37], triggering the need for fluid boluses or vasopressors before admission. Fortunately, the JH data recorded the start and end times of critical care therapies (eg, high-flow nasal cannula [HFNC], oxygen flow rate [L/min], mechanical ventilation, fraction of inspired oxygen (FiO₂; %), vasopressors, and dialysis), preadmission treatment, and vital sign information, enabling us to identify the minimums for SpO₂ and systolic blood pressure before treatment. The potential for SpO₂ bias owing to supplemental oxygen was also mitigated by the derived ratio of SpO₂ to FiO₂. This ratio was calculated using either an FiO₂ value recorded contemporaneously with SpO₂ (including cases in which FiO₂ was recorded as 21%, suggesting that SpO₂ was a “room air” measurement) or in cases in which an oxygen flow rate was documented (eg, in cases using nasal cannulas) using an estimated FiO₂ calculated from the oxygen flow rate [38].

Our study explored the detection of phenotypes by clustering routinely available clinical data. However, raw clinical data typically extracted automatically from EHRs can often contain outliers, particularly those associated with observations that may have been manually entered erroneously [20]. Recent studies have confirmed that outliers will negatively affect the quality of derived clusters [39]. Although the JH-CROWN Registry contained syntax error–free structured tables for vitals and laboratory test measurements, unlikely outliers, as recorded in the EHR, were replicated in the registry tables. To cleanse vitals, we adopted reported rules reflecting commonly accepted ranges [40] for human physiology in which outliers were replaced with “null” (ie, treated as “missing”). Table 2 shows the raw total counts of key vitals—with SpO₂, pulse, and respiration having the highest raw counts followed by temperature and blood pressure—and the statistics of the validated (cleansed) vital signs. For laboratory tests, values were replaced with nulls using the statistical “3(IQR)” criteria designed to detect extreme outliers in observational data [41].

Multiple imputation (MI) and weighted clustering analysis were applied to the training and validation cohorts independently. No validation data were used to influence the imputation of the training data. MI, known to reduce bias even when the proportion of missingness is large [42], is an approach to missing data whereby multiple copies of the feature data set are generated with missing values replaced by inferences drawn from the data set. Our approach was based on Bayesian joint models congenial or compatible with k-means clustering [43,44]. The joint-modeling MI was based on the Dirichlet process mixture of multivariate normal distributions to reflect complex distributional features [45]. For each cohort, a total of 100 imputed data sets were created. K-means clustering was then applied to each imputed data set, generating base clusterings for final weighted consensus clustering.

Although an old rule of thumb is that 3 to 10 imputations would typically suffice to ensure precision and replicability [46], recent studies [47] have developed a new formula based on the fraction of missing information (FMI) that estimates how many imputations would be needed for precise and replicable SE, CIs, t statistics, and P values [47]. Although, to achieve a variation of <5% in SE at an FMI of 25%, the estimated number of imputations needed is approximately 20 [48], as the FMI increases, the required number of imputations increases quadratically, and at an FMI of 70%, the estimated number of required imputations is approximately 100 [47]. As, in general, adding more imputations increases precision and replicability, we chose 100 imputations to ensure robustness and accuracy, especially in the context of complex data such as the medical records of patients with COVID-19 with a substantial amount of missing data [49].

Although there are numerous methods that have been proposed to determine the optimal number of partitions or clusters in k-means analysis, clustering stability has emerged as a general model-agnostic evaluation method. In statistical learning terms, if data sets are repeatedly and randomly sampled from the same underlying distribution, a stable clustering algorithm should find similar partitions [50]. The approach used in our study defines a “good clustering” in terms of its instability in response to imputation-related perturbations in the data. Instability was assessed using the bootstrapping method [51]. Accordingly, we selected k as the value that minimizes the instability of the clustering [52]. Instability-based methods are attractive as they are not based on a specific metric for the distance between objects and have been shown to perform at least as well as state-of-the-art distance-based methods [53].

Consensus clustering has the theoretical advantage of minimizing overfitting and optimizing the stability of cluster assignments, as has been shown for identifying subgroups of heterogeneous patients in the ICU [10]. A weighted consensus clustering based on the nonnegative matrix factorization (NMF) framework [54] was obtained for the clustering results from all imputed data sets [55]. Unlike a consensus approach based on an averaging process (wherein all base clusterings are treated with equal weight), the objective of NMF weighted consensus clustering is to aggregate the base clusterings into a final clustering using weights optimized for each base clustering in a manner analogous to least absolute shrinkage and selection operator regression [56]. Under this approach, the solution for the weights is sparse (ie, only a small subset of base clusterings contributes to the final clustering). A wide range of comparative experiments has demonstrated the effectiveness of the NMF-based consensus clustering approach [54,55]. The R package clusterMI (version 0.0.41; R Foundation for Statistical Computing) [43] was used to perform clustering with MI. In addition, it allows for the consensus pooling of results in terms of both partitions and instability [57].

To support the analysis of reproducibility, a statistical analysis of how features and outcomes were distributed across phenotypes in both cohorts, with P values established using the Kruskal-Wallis rank sum test for continuous variables and the chi-square test for categorical values [58], was performed. An overview of the clustering process flow is shown in Figure 3.

Although the identification of phenotypes via clustering has the potential to inform personalized care [59], the lack of point-of-care testing of key defining inflammatory biomarkers, especially during the early stages of an encounter, limits the clinical utility of phenotypes. A recent related ARDS study explored the application of supervised GBM phenotype classifiers trained using routinely available observational data and clustering-identified labels and achieved a phenotype classifier with an area under the curve (AUC) of 0.95 [24]. Our study extends this modeling effort by deriving a predictive GBM phenotype classifier trained using data observed within the first 6 hours of admission in which, as in the cited ARDS study, model performance was evaluated against the clustering-derived phenotype.

Our approach to phenotype prediction model development adheres to the TRIPOD (Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis) guidelines [60]. For prediction modeling, routinely available observational data, as shown in Textbox 1 (except for the C-reactive protein [CRP], which was excluded because of >40% missingness), for included patients (N=4379) recorded within 6 hours following admission were used for training. MI based on Bayesian joint models was applied to create 100 complete feature data sets of the remaining 37 features used for phenotype identification. On each imputed data set, 69.99% (3065/4379) of the included patients were randomly selected as the training set, whereas the remaining 30.01% (1314/4379) were reserved as the test data set. The imputeData() function from the aforementioned clusterMI R package [57] was used to perform MIs of the early feature data before random 70/30 splitting.

The GBM was trained with 10-fold cross-validation for hyperparameter tuning using a grid search to optimize the models using the 100 training sets. Prediction performance (eg, AUC, sensitivity, and specificity) was assessed on the held-out test sets. The final performance metrics were estimated by averaging the performance estimates obtained from each imputed data set. An overview of the phenotype prediction process flow is shown in Figure 4. The classification and training R package caret (version 6.0-93) [61] was used for prediction phenotype classifier development.

This study was approved by the JH institutional review board (IRB00250903).

By examining the total instability over different numbers of clusters [52], 2 clusters were found to be optimal in both the training and validation cohorts as, in both cases, k=2 exhibited the least instability (Figure 5). The final assignment of each patient to 1 of the 2 phenotypes in each cohort (phenotype 1: 1284/4379, 29.32% and 1258/4379, 28.73%; phenotype 2: 898/4379, 20.51% and 939/4379, 21.44% in the training and validation cohorts, respectively) was determined by NMF consensus clustering using 2 clusters.

Figure 6 depicts rank plots in which the 38 features used for training and validation cohort clustering are normalized with respect to the mean and SD of the population of the underlying paired phenotypes. Between-phenotype comparisons through nonparametric statistical methods indicate that, among the considered features in both cohorts, the most significant phenotype-defining features include age, blood urea nitrogen, creatinine, and elevated inflammatory laboratory values (neutrophils, neutrophil-to-lymphocyte ratio, red blood cell count, and albumin).

Figure 7 depicts violin plots of the clustered training and validation data features. In this display of the summary statistics, distribution, and density of each variable, it appears that features across phenotypes share similar distributions and densities.

Figure 8 shows the differences in inflammatory biomarkers (CRP, interleukin 6, D-dimer, ferritin, lactate dehydrogenase, procalcitonin, and fibrinogen) in both cohorts associated with poor COVID-19 outcomes, as reported in previous studies [62,63]. In both the training and validation cohorts, phenotype 2 was associated with elevated inflammatory markers.

Table 3 and Figure 9 show the odds ratios of phenotype 2 versus phenotype 1 associated with comorbidities adjusted for age, race, gender, and ethnicity in both cohorts. These results suggest that patients in phenotype 2 have a higher likelihood of anemias, lymphoma, coagulopathy, congestive heart failure, preexisting renal failure, peripheral vascular disease, AIDS, complicated hypertension, bleeding peptic ulcers, cancer, electrolyte disorders, and diabetes with chronic complications. Figures 10 and 11 are principal-component analysis biplots including a scatterplot that shows the similarity of 2D projections of clustered observations or patients. These figures have a superimposed loading plot that shows how strongly features influence a phenotype (eg, strong associations between phenotype 1 and lymphocytes, SpO₂/FiO₂, and albumin and between phenotype 2 and systolic blood pressure, age, creatinine, and pulse pressure).

The detailed demographics, clinical characteristics, and statistical significance of the feature and outcome distribution across phenotypes in both the training and validation cohorts are shown in Tables 4 and 5. In these tables, P values suggest statistically significant associations between the need for ICU-level care and poor outcomes associated with phenotype 2 in both cohorts. Specifically, as shown in Table 4, phenotype-2 patients were associated with advanced age (mean 76, SD 14.2 years in phenotype 2 vs 52 years in phenotype 1) and with statistically significant (P<.001) increased risk of developing sepsis (34% in phenotype 2 vs 21% in phenotype 1), requiring mechanical ventilation (11% in phenotype 2 vs 4.5% in phenotype 1), using vasopressors (10.5% in phenotype 2 vs 3.5% in phenotype 1), requiring HFNC (16% in phenotype 2 vs 8% in phenotype 1), requiring continuous renal replacement therapy (CRRT; 2.4% in phenotype 2 vs 0.5% in phenotype 1), requiring dialysis (8% in phenotype 2 vs 0.6% in phenotype 1), and mortality (17% in phenotype 2 vs 2.5% in phenotype 1).

Figure 12 shows survival curves. In both cohorts, survival between phenotypes diverged on day 1 from admission, and the divergence was sustained over 60 days, with significantly lower survival in the phenotype 2 group than in the phenotype 1 group.

The 4379 patients used for phenotype prediction analysis did not include those who initiated some form of critical care therapy (eg, HFNC or mechanical ventilation) or died within the 6 hours following admission. Clustering-identified phenotypes and clinical features (Textbox 1) based on observations recorded within the first 6 hours associated with the included patients were used for GBM predictive classifier derivation except for CRP, which was excluded because of excessive missingness (>40%). The classifier performance was based on comparing clustering-identified phenotype labels with predicted labels in the held-out test sets. Table 6 shows predictive metrics with 95% CIs of the classifier performance over the 100 imputed test sets (1314/4379, 30.01%); the mean AUC to accurately predict a test patient’s clustering-derived phenotype was at 0.89 (95% CI 0.887-0.893).

The key aim of our study was to develop the foundations of an EHR data-screening tool that may assist clinicians in the early identification of patients among a population of highly heterogeneous hospitalized patients with COVID-19 likely to deteriorate to hyperinflammation and require ICU-level care that may include respiratory support across a broad spectrum of modalities, such as HFNC, Nasal intermittent positive pressure ventilation, intubation or mechanical ventilation, and extracorporeal membrane oxygenation. Patients with hyperinflammation may also develop life-threatening comorbidities such as septic shock [64] and AKI [65], driving the need for specialized care, such as vasopressors [66], intermittent dialysis, and CRRT [67].

The heterogeneity of hospitalized patients with COVID-19 [68] suggests the potential benefits of clustering encounter data to identify phenotypes with distinct host response patterns to treatment that may help guide personalized therapeutics. Recent studies in Europe and the United States have identified 2 homogeneous clinical phenotypes in hospitalized patients with COVID-19 using machine learning or clustering algorithms with the potential utility to identify targeted treatment protocols [21,22].

Given that, in its most severe form, SARS-CoV-2 infections lead to life-threatening pneumonia and ARDS, clustering studies identifying phenotypes of patients with ARDS or who are mechanically ventilated and at risk of ARDS [69] are highly relevant. In total, 2 ARDS phenotypes (hyper- and hypoinflammatory) have been consistently identified in previous clustering studies [16,24,31,59,69-72] that have statistically significant similar clinical, physiological, or biomarker traits, including differential responses to treatments, interventions, and mortality rates, supporting the potential utility of machine learning–based ARDS phenotyping [27,59].

Using k-means clustering analysis of training cohort clinical data, we identified 2 distinct phenotypes that differed significantly in demographics, sepsis incidence, inflammatory biomarkers, the need for ICU-level care, and clinical outcomes including mortality. This result was reproduced in an independent clustering analysis of an internal validation cohort. These findings suggest that the early association of a new patient with a clustering-identified phenotype may provide useful prognostic information. For example, a hospitalized patient predicted to be phenotype 2 and not in current need of supplemental oxygen may be viewed as at high risk of progression to requiring ICU-level care. This patient could be flagged by hospital staff so that they might initiate close monitoring, offer empiric use of therapies such as remdesivir [73], and prepare critical care resources. Alternatively, a patient predicted to be phenotype 1 and not in current need of supplemental oxygen may be viewed as low risk, warranting supportive care only. To facilitate the early identification of patient phenotypes, this study developed a predictive GBM classifier with a mean AUC of 0.89, which would be considered an excellent statistical performance [74]. In addition, the GBM classifier used only routinely available vital signs and laboratory results observed within the first 6 hours of admission, enhancing the value of this tool as an early warning system.

To our knowledge, this is the largest clustering study identifying homogeneous phenotypes in hospitalized patients with COVID-19 using routinely available early clinical data. Independent clustering analysis of randomly selected patients in the training and validation cohorts identified a hyperinflammatory phenotype (phenotype 2) characterized by higher plasma levels of inflammatory biomarkers that were associated with a higher prevalence of HFNC, invasive mechanical ventilation, extracorporeal membrane oxygenation, CRRT, dialysis, vasopressor use, diagnosis of sepsis or ARDS, and increased mortality compared with a hypoinflammatory phenotype (phenotype 1).

Recent reports have concluded that the imbalance between hyperinflammation and immune paralysis is a hallmark of sepsis [75] and that the levels of inflammatory biomarkers such as interleukin 6 in patients with COVID-19 are associated with mortality [76]. As inherent characteristics and genetic predisposition are likely key to the heterogeneity of individual immune responses, the ability to categorize patients based on the risk of hyperinflammation allows for risk stratification and personalized treatment using targeted therapeutic regimens. The ability to identify 2 separate phenotypes based on immune condition also allows for specific treatment approaches using immunomodulators.

Our findings are in concert with those of other studies that have associated worse COVID-19 outcomes with comorbid conditions that include depression [77], anemia [78], hypertension [79], congestive heart failure [80], preexisting renal failure [81], peripheral vascular disease [82], cancer [83], paralysis or spinal cord injuries [84], chronic obstructive pulmonary disease [85], obesity [86], electrolyte imbalance [87], and diabetes [88].

Although our 2-phenotype findings are similar to those of the 2021 study of 483 patients with COVID-19 at Yale New Haven Health [21], there are several differences. In the Yale study, among the 2 identified phenotypes, the phenotype with the higher risk of mortality comprised older individuals with more comorbidities, whereas patients in the group with a lower risk of mortality comprised younger individuals who were more likely to be obese, male, and racial and ethnic minority individuals with higher levels of the CRP and ALT inflammatory markers. In contrast, our analysis identified a hyperinflammatory phenotype associated with age and comorbidities highly relevant to the development of ARDS and sepsis, both leading to increased mortality rates. However, in both the training and validation cohorts, we found that elevated ALT and BMI, male gender, and racial and ethnic minority individuals were associated with the hypoinflammatory phenotype 1 (Table 3), which is consistent with the results of the Yale study. Although the distribution of mortality between the 2 phenotypes was similar (Yale study: 25% vs 9%; this study: 23% vs 3%), the Yale study did not find statistically significant differences in the use of critical care treatments (eg, dialysis or mechanical ventilation) between the 2 phenotypes. Overall, the Yale study showed that patients who were admitted for COVID-19 were found to be classified into 2 cohorts mostly based on age-related comorbidities and specific demographics. It should be noted that both the Yale study and our study support the recent finding that, although there may be an increased incidence of severe COVID-19 among Black and Hispanic patients, this is not due to an inherent susceptibility to progression [89].

A recent systematic review of prediction models for COVID-19 [60] enumerated common weaknesses. These include a high risk of bias from inadequate sample sizes and inappropriate or incomplete evaluation of model performance with insufficient internal or external validation. In addition, calibration was often incomplete or performed using inappropriate statistics. Finally, inappropriate handling of missing data was common, including the omission of how missing data were handled. The authors summarily recommended that prediction modelers “should adhere to the TRIPOD (Transparent Reporting of a multivariate prediction model for Individual Prognosis Or Diagnosis) reporting guideline” [60]. We believe that a key strength of this study is the proactive adherence to the TRIPOD guidelines [90], thereby avoiding the weaknesses described in previous publications [91].

Our study serves as a proof of concept that combines unsupervised clustering for COVID-19 phenotype identification in historical data and supervised machine learning for phenotype prediction model derivation using routine clinical data, which is feasible as a basis for an “early warning” bedside COVID-19 screening tool [92]. If validated prospectively, such EHR data–derived and embedded models could automatically incorporate and analyze clinical data to provide real-time COVID-19 critical care decision support while minimizing disruptions to the workflow.

The prospective validation must address 2 factors. First, it is recognized that COVID-19 populations may differ significantly across time and geography with changing availability or use of vaccines, circulating SARS-CoV-2 variants, treatments, and the influence of comorbidities such as seasonal influenza and respiratory syncytial virus. Hence, the models derived in our study to identify or predict phenotypes need to be routinely “retrained” to reflect hospitalized populations with varying characteristics. With continuous access to EHR data, this issue could be addressed through machine learning models that are routinely updated with changing inpatient population characteristics.

Second, it must be prospectively demonstrated that the models can classify phenotypes robustly and consistently in real-time clinical scenarios in diverse settings. Before their clinical implementation, the models will need rigorous evaluation of their interaction with missing data frequently encountered in the real-world setting of critical care. Although we used a robust set of 38 features combined with imputation, it may be that a more effective approach would involve fewer readily available features that might predict membership to a phenotype with sufficient accuracy. The development and validation of such parsimonious models [92] require a careful analysis of the most important phenotype-defining features that would also most likely be reliably available during the early stages of an encounter. Moreover, although we used observations recorded within the first 6-hour window following admission to derive our predictive model, multisite studies have shown that the mean length of stay for patients with COVID-19 requiring ICU-level care ranges from 12 to 19 days. This suggests that predictive models trained using data over longer intervals (eg, recorded within 24 hours following admission, decreasing the prediction horizon [93]) or updating the prediction longitudinally [94] may lead to clinically useful models with improved prediction performance [95].

Concerning potential methodological weaknesses of this study, it should be noted that, in a head-to-head comparison of LCA versus k-means clustering in a relatively small sample of pediatric patients with sepsis (n=151), LCA was found to be somewhat more useful in identifying homogeneous phenotypes. However, both approaches identified at least one distinct high-severity phenotype [96]. Given that LCA is computationally challenging whereas k-means is better scaled to large data sets [97], most critical medicine clustering studies involving large cohorts (N>1000) in sepsis [20], ARDS [98], and COVID-19 [99,100] have effectively used k-means to identify well-separated phenotypes, leading to early detection of those who would benefit from certain treatments and close monitoring. Notably, a large cohort study by Seymour et al [20] identified and validated 4 clinical phenotypes of sepsis through k-means clustering analysis that were positively correlated with host response patterns and clinical outcomes. Most recently, Duggal et al [98] reported a k-means analysis of routine clinical data associated with a large cohort (4773 patients) that successfully identified 2 distinct ARDS phenotypes that included a phenotype with increased levels of proinflammatory markers, higher mortality, and longer duration of ventilation compared with patients in the second phenotype [99]. These studies support the validity of k-means as an effective machine learning technology for the identification of clinically useful phenotypes in “big EHR data” studies.

Another methodological weakness concerns the dependence on instability analysis (Figure 5) to identify the optimal number of phenotypes. Studies have shown that stability-based methods can be sensitive to underlying data distributions and may not always provide a valid and meaningful choice of the optimal number of k-means–derived clusters [101]. Although instability-based methods compare favorably with commonly used distance-based methods to identify the optimal k (eg, elbow and silhouette), alternatives such as the Calinski-Harabasz [102] evaluation metric that measures the compactness and separation of clusters, thereby providing a measure of the quality of the clustering results, would be a useful addition to the analysis. As this metric can be sensitive to the density and shape of clusters, in future studies, it may be beneficial to consider both stability and evaluation metrics when selecting an optimal k [103]. However, the validity of our finding that k=2 identifies the true number of COVID-19 phenotypes is bolstered by our use of a recently improved instability metric that corrects for the distribution of cluster sizes [52] and the fact that independent studies in other related populations (ARDS and other populations with COVID-19) have also identified 2 phenotypes using totally different clustering algorithms (eg, LCA).

In summary, k-means clustering was effective in identifying phenotypes with distinct treatments or intervention responses and outcomes in a large cohort of hospitalized patients with COVID-19. In addition, a GBM machine learning classifier model using readily available early encounter data accurately assigned patients to phenotypes, suggesting that the application of these models in a clinical setting may provide valuable prognostic information that could inform personalized COVID-19 management. Although future studies and trials are needed to validate the clinical utility of phenotype assignment, it would seem reasonable to implement successfully validated machine learning algorithms in extant EHR systems as a tool to support those trials.