The global spread of COVID-19, the disease caused by SARS-CoV-2, has resulted in the loss of around 700,000 lives. Repercussions of the pandemic have wreaked havoc on the economy, sending billions of people into lockdown to flatten the curve for the healthcare systems and resulting in record-high unemployment around the world. The American Hospital Association estimates a total four-month financial impact of over $200 billion in losses for the U.S. healthcare systems as a result of cancelled hospital services (e.g., cancelled non-elective surgeries and outpatient treatment) due to the COVID-19 pandemic and COVID-19 hospitalizations. 1 Our inability to provide precise estimates of the COVID-19 outcomes such as death, hospitalization, and need for ICU and ventilation has contributed to lost opportunities for saving lives with personalized preventive measures and making smart resource allocation plans. Although our inability to predict was partly due to the novelty of the disease, the critical question is: do we have the data and technology to predict?
Over the past decade, the U.S. federal government has made extensive investments to institute meaningful use of electronic health record (EHR) systems. Clinical data in EHRs, however, are still complex and have important quality issues, impeding their leverage to address pressing health issues that require rapid response. Nevertheless, biomedical researchers are increasingly applying data mining and Machine Learning techniques to clinical data for predicting health outcomes.
Recent studies have shown that COVID-19 disease severity and mortality are associated with a number of comorbidities including cardiovascular disease, diabetes mellitus, hypertension, chronic lung disease, cancer, chronic kidney disease, and obesity, and demographics including age, sex, and race/ethnicity. [2] [3] [4] [5] [6] [7] [8] A number of Machine Learning (ML) models have been developed to predict susceptibility in the general population, the likelihood of a positive diagnosis in a patient with symptoms, and prognosis in those with the disease. 9 Many of these models are based on a combination of demographics, comorbidities, symptoms, and biomarkers, [10] [11] [12] [13] [14] [15] and use data from relatively small cohorts of COVID-19 patients. Jiangfent et al, 10 for example, performed a logistic regression analysis on 299 patients that identified age, lymphocyte count, lactate dehydrogenase, and oxygen saturation as independent predictors for mortality. Jiang et al, 15 used data from 53 patients to develop ML models to identify elevated alanine aminotransferase (ALT), the presence of myalgias, and an elevated hemoglobin as the most predictive features for disease severity, achieving approximately 70-80% accuracy in modeling COVID-19 morality. Huang et al 11 examined clinical data from 125 COVID-19 patients and identified the presence of comorbidities, increased respiratory rate, elevated C-reactive protein, and elevated lactate dehydrogenase as independently associated with a worse prognosis. While the inclusion of vital signs and biomarkers in these models may be highly predictive of some adverse outcomes of COVID-19 infection, they are typically measured after the patient has already started to show signs of the disease and may be at a point that is too late for a useful intervention.
Estiri et al. recently introduced the transitive Sequential Pattern Mining (tSPM) along with a dimensionality reduction algorithm, MSMR, and showed that together, the algorithms can successfully predict different health outcomes. 16, 17 The goal in tSPM is to mine temporal data representations from clinical data for application in downstream ML. The number of transitive sequential representations is quadratic in the number of records, and thus create dimensionality issues. To address this, the MSMR algorithm --short form stands for Minimize Sparsity, Maximize Relevance, applies high performance dimensionality reduction ( Figure 1 ). It takes the initial set of N features mined by the tSPM algorithm and provides a list of <400 features, through a 3-step process including frequency-based and information-based (mutual information, and joint mutual information) filtering of a large number of features. 16, 17 We adapt the tSPM temporal representation mining and MSMR dimensionality reduction algorithms, and make adjustments to architect an end-to-end ML pipeline that enables iterative feature and algorithm selection to predict health outcomes (MLHO, pronounced as melo). MLHO offers an architecture that enables a parallel outcome-targeted calibration of the features and algorithms. The goal in MLHO is to mine relevant data representations from past clinical records and select the most efficient algorithmic solution for accurately predicting future health outcomes. As COVID-19 cases are rising again in many areas of the world, being able to provide personalized predictions of its adverse outcomes that are directly connected to the healthcare systems' responsiveness and mortality can be game-changing. In this study, we focus on predicting four outcomes (hospitalization, ICU, ventilation, and death) in patients with a verified COVID-19 infection, using their past medical records. Using about 600 features mined from patients' past medical records (before contracting COVID-19), we trained and tested predictive models that estimate risks of hospitalization, ICU admission, need for mechanical ventilation, and death.
Implementation of the MLHO pipeline is depicted in Figure 1 . MLHO mines both sequential (temporal) and raw data representations from clinical data and performs iterative feature and algorithm selection in a 2-step evaluation process. MLHO's architecture enables a parallel outcome-targeted calibration of the features and algorithms, in which different statistical learning algorithms and vectors of features are simultaneously tested and leveraged to improve prediction of health outcomes. The four adverse outcomes of interest in this study reflect a hypothetical sequential spectrum of outcome severity in patients with verified COVID-19 infection, ranging from hospitalization to ICU admission, to need for mechanical ventilation, and ultimately to death (Hospitalization → ICU → Ventilation → Death). 
We used electronic health records data from over 14,500 patients with a confirmed case for COVID-19 between March and May 2020 and who had at least 1 year of medical history with Mass General Brigham (MGB), since 2016. Table 1S presents general demographic information about the study cohort. For each patient, we only included clinical records from 14 days prior to the positive COVID-19 test date. This temporal buffer ensures that no COVID-19-related medical conditions are included in the model as a risk factor. The use of data for this study was approved by the Mass General Brigham Institutional Review Board (2020P001063).
We randomly split the data with an 80-20 training-to-testing ratio. We iterated the train-test sampling 10 times to account for possible patient population differences caused by the sampling. We used the tSPM algorithm to mine temporal sequential representations. Given a list of clinical records for patient at times , the tSPM algorithm
mines all avector of transitive sequential patterns from possible pairs of distinct records X ij where by setting (as samples of random variable for patient ) to be 1 R , ),
if and and 0 otherwise. We also mined raw representations from the clinical 1, k k i ≥ j ≥ 1 , t i1 ≤ t ji data, which are composed of all clinical records. To obtain a count of the raw records, for each patient , we sum the frequency of , as samples of a random variable . p , k , ...,
On the training sets, we performed feature and algorithm selection and the final predictive modeling. We feed the combined representations to the feature selection step.
Unlike other studies that begin with a limited set of hypothetical risk factors, we took a primarily inductive approach to selecting clinical representations for predicting the risk of adverse outcomes in COVID-19 patients. First, we apply a filter method for feature selection, using a computational algorithm that minimizes sparsity and maximizes relevance (MSMR 16, 17 ) .
Step 1 in MSMR is to cut the initial combined vector of representations to those that were observed in X′ fewer than 0.2 percent of the patients. On the remaining representations, step 2 in MSMR computes the mutual information with the outcome variable which in this study (y , ..., y ), Y = 1 i includes labels for hospitalization, ICU, ventilation, and death. Mutual information, 18, 19 in this case, measures the amount of information that each remaining representation contains about the outcome. Given the joint probability distribution , the mutual information between
We cut the remaining representations from by ranking based on the mutual information X′ coefficient, and update to a list of around 30,000 representations with the highest mutual X′ information. In the third and final step, MSMR computes the joint mutual information (JMI) 20 score for the updated vector of remaining representations, . The algorithm starts with a set X′ S containing the top feature according to mutual information, then iteratively adds to the features S maximizing the joint mutual information score
Where the random variable corresponds to the joint distribution of and . As a result, X X′ ′ * X′ X′ * JMI also takes into account the redundancy between the features -i.e., reducing multicollinearity among covariates. Second, we combine the feature selection with a preliminary evaluation of algorithms for the prediction task.
Using the <400 features identified by the MSMR algorithm, we train a set of preliminary classification algorithms that perform embedded feature selection. The preliminary classification serves two goals. First, we screened features used in those algorithms to compile a list of the common features used for modeling the 4 outcomes. Second, during each sampling iteration, we computed the Area Under the Receiver Operating Characteristic Curve (AUC ROC) on the held-out test sets to evaluate the algorithms' performance for predicting the outcome labels. We used 10-fold cross-validation to train the prediction algorithm --therefore, a 72-8-20, train-evaluation-test split. We tested 10 classification algorithms --bartMachine: Bayesian Additive Regression Trees, 21, 22 Stacked AutoEncoder Deep Neural Network (dnn), Stochastic Gradient Boosting (gbm), 23, 24 glmboost: Boosted Generalized Linear Model, eXtreme Gradient Boosting (xgb) 25 with DART booster (xgbDART), 26 linear model solver (xgbLinear), tree learning (xgbTree), model-averaged Neural Network (avNNet), Elastic-Net regularized generalized linear model (glmnet), 27, 28 and Oblique random forest (ORFlog). 29 
The iterative feature and algorithm (preliminary classification) selection results in a set of common features and a ranking of classification algorithms. Using the top algorithms (by AUC ROC) and the common features, we perform a final round of classification training with 10 train-test sampling iteration and 10-fold cross-validation. We also measured a model-specific feature importance/influence metric in the final modeling rounds.
As we shall reveal in the results, the gbm algorithm was one of the top 2 performing algorithms. As a result, we computed the relative feature influence metrics from the final gbm models to measure features' importance. In the gbm model, boosting estimates as an 'additive'
Where the expansion coefficients and the base learner's --
are jointly fit to the training data in a feedforward process. A regression tree model specializes tue base learner partitions the feature space into disjoint regions to (X;
)
predict a different constant value for each. 30 The relative influence (or contribution), averages I j 2 the improvement made by each variable when it is permuted from all trees in which the given variable was incorporated. 30, 31 In an additive tree model, the relative influence measure is provided by Friedman and Meulman (2003) (T )
where is the measure of relevance for a single tree . 32 (T ) I j 2 T We summarize and report the AUC ROCs from the final models as well as the features' relative influence values. Table 1S provides a summary of demographic characteristics of the patients. We found that among racial/ethnic groups, the rate of hospitalization, ICU admission, and ventilation was highest among African American/black individuals. The overall mortality rate was three percent, and white individuals had the highest rate (3.9 percent) compared to other racial/ethnic groups. Compared to females, male COVID-19 patients had a significantly higher chance of hospitalization, ICU admission, ventilation, and mortality. While the average age in the cohort of COVID-19 patients was 49.5, the average age of patients who were hospitalized, admitted to ICU, or needed mechanical ventilation was between 61 and 62. The average age of mortality in COVID-19 patients was much higher at 77 years (Table 1S ). Figure 2 summarizes the sequential scenarios in which the four adverse outcomes are observed in COVID-19 patients. Overall, in more than 75 percent of the patients we did not observe any adverse outcomes. Approximately 13 percent of patients were discharged after hospitalization. The cumulative probability of patients needing to be admitted to the ICU was less than eight percent, which sorted in a declining order of sequential events from four percent to 1.4 percent for Hospitalized → ICU → Discharged, Hospitalized → ICU → Ventilation → Discharged, Hospitalized → ICU → Ventilation → Died. This order supports our general hypothesis about the severity spectrum. Of the three percent overall chance of mortality, 2.5 percent fall into the sequential scenarios of Hospitalized → ICU → Ventilation → Died and Hospitalized → Died (without ICU/Ventilation). Iterative feature and algorithm selection MLHO mined over 60 thousand raw (e.g., diagnosis/medication/procedure codes) and 160 million transitive sequential (e.g., medication → diagnosis) representations. Through the iterative feature selection, these features were shrunk to over 2,300 representations, about 900 of which were raw features and about 1,400 were transitive sequential features. As described in the methods, both MSMR filter method and embedded methods (while training preliminary classification algorithms) were utilized in iterative sampling of train-test data. The other outcome of this step was identification of the top predictive algorithms. Figure 3 illustrates the AUC ROC results obtained from the 10 algorithms for estimating the risk of hospitalization. We found that two Boosting algorithms obtained the best overall results: the Stochastic Gradient Boosting 
We fed the shrunken feature set (with ~2,300 features) to the two top algorithms for the final model training --with 10 iterative train-test sampling and 10-fold cross-validation, which means for each outcome, we trained 20 final models. We picked the top 10 models, by AUC ROC. Table 1 presents the median, standard deviation, and the best AUC ROC values obtained from the final models for each outcome. Overall, the median AUC ROCs was the best for mortality, while it ranged near the 0.8 threshold for hospitalization, ICU, and ventilation. Among the top 10 models, the gbm algorithm was the best algorithm for hospitalization and ICU, and the xgbDART was the best algorithm for predicting the need for mechanical ventilation. 
Because the relative influence measures from the gradient boosting models are relative, by default the largest value is assigned to 100, and the remaining values are scaled accordingly. 30 Using only the non-zero relative influence values resulted in a final set of 603 features that were used for predicting at least one of the outcomes. Table 2S provides the list of features with their median and interquartile range values of relative influence by outcome.
The features were then grouped under nine different clusters that correspond to comorbidities and demographics that have been associated with poor outcomes in COVID-19 by the Center for Disease and Control 33 . The categories are cardiovascular disease, chronic renal disease, chronic lung disease, neurological disorder, diabetes mellitus, other chronic illness, age, sex, and race/ethnicity. The "Other chronic illness" cluster includes comorbidities such as hypertension, hyperlipidemia, hypothyroidism, obstructive sleep apnea, and obesity. Table 2 is a sample of features that were assigned to two of these disease clusters. The features used include diagnoses, medications, laboratory tests, imaging studies, and procedures. Diagnoses can easily be identified as a subtype of a particular category. For example, "Type 2 diabetes mellitus with diabetic nephropathy" is grouped with diabetes mellitus and "Unspecified dementia without behavioral disturbance" is grouped with neurological disorder. Medication labels are assigned to a specific cluster if the medication is exclusively used for that disease. For example, "metformin" is assigned to the diabetes mellitus cluster because it is primarily used for treating diabetes mellitus. However, in many cases a medication can be used for multiple diseases and cannot be easily grouped with a specific cluster. If a medication cannot be attributed to a specific cluster, it was left unassigned.
For the sequential features, if one of the two components in the sequence was related to a disease cluster; the whole sequence was assigned to the corresponding disease. Laboratory tests, imaging studies and procedures on their own are not associated with a particular disease. But they were assigned to particular disease clusters if they were in sequence with another label that was specific to a disease cluster. For example, the sequence "CT Head → Altered Mental Status" was associated with a neurological disorder because it includes the diagnosis "Altered Mental Status". Being able to predict the four outcomes can also create possibilities for better preventive measures. For instance, healthcare systems and regional health authorities can use predicted outcomes to identify patients who might be at higher risks and plan for preventive measures such as alerting the patients' primary care providers to prioritize healthcare maintenance visits in order to make sure the patients are up to date on immunizations and health care screening. There could also be an in-person discussion of the importance of taking precaution against contracting COVID-19. The estimated risks provided by MLHO can also help stratify the most at-risk from already known vulnerable populations such as those living in nursing homes, subacute rehabilitation centers, and homeless shelters for potential high risk.
We found that the average age for patients who died was about 77, whereas the average age for hospitalization, ventilation, and ICU admission was around 61. This shows that the patients hospitalized and having critical illness interventions (such as an ICU admission or intubation) in general have a younger age. One possible explanation for this discrepancy is that many of the younger patients with these interventions may survive. This suggests that the younger patients who get intensive therapies such as ventilation or admission into the ICU, see an improvement compared to the older patients who may have less improvement with such therapies.
The relative influence of specific disease clusters for different outcomes suggests certain diseases may have increased predictive value for determining adverse outcomes. Age is by far the best predictor for an adverse outcome. It is far more influential of the prediction than any other comorbidity or demographic including severe diseases such as congestive heart failure, end stage renal failure, or dementia. Another interesting trend between outcomes is that chronic illness (hypertension, hypothyroidism, and obesity) and diabetes mellitus are relatively more influential predictors for hospitalization, ICU admission, and need for ventilation than they are for death. However, the comorbidity with the best prediction for death is chronic lung disease. This could suggest that intensive therapies such as ventilation and admission into the ICU, may improve the outcome of patients with common comorbidities such as hypertension and hyperlipidemia but they do not have as beneficial of an impact on patients with chronic lung disease.
Findings of this study can exhibit limitations due to the potential availability of noise in outcome labels. Outcomes such as hospitalization, ICU admission, and ventilation can be difficult to uniformly measure since they can be correlated with healthcare utilization. For example, hospitalization could be due to a heart failure exacerbation and the COVID-19 diagnosis is only incidentally discovered, but unrelated to the hospitalization.
The COVID-19 pandemic has had devastating health and economic impacts across the globe. Being able to predict adverse outcomes that directly impact the capacity of healthcare systems in a timely manner can lead to better allocation of healthcare resources and more efficient targeted economic and preventive measures. MLHO's architecture enables a parallel outcome-targeted calibration of the features and algorithms, in which different statistical learning algorithms and vectors of features are simultaneously tested and leveraged to improve prediction of health outcomes. As the COVID-19 cases are re-surging around the world, a pipeline like MLHO is crucial to improve our readiness for confronting not only the potential future waves of COVID-19, but also other novel infectious diseases that may emerge.
This work was funded through the National Human Genome Research Institute grants R01-HG009174.