Coronavirus disease 2019 (COVID-19), an illness caused by novel severe acute respiratory syndrome coronavirus 2, has spread rapidly across the world with over 9.5 million cases internationally and over 2.4 million cases in the United States as of June 25, 2020 [1] .
While the popularity of chest computed tomography (CT) in the early detection and monitoring of COVID-19 is growing in many countries such as China, hospitals in the United States employ chest radiographs (CXRs) as the primary imaging modality for the monitoring of the disease [2] [3] [4] [5] [6] [7] [8] . CXRs are useful due to the speed, portability, and easy disinfection of radiography units. The American College of Radiology has suggested that CT be reserved for only severe cases of COVID-19 and strongly recommends the usage of portable radiography units to minimize spread of infection [8] . However, CXRs have lower resolution than CT images and provide 2-Dimensional (2D) rather than 3D representations of the lungs. Additionally, portable radiography units may result in non-uniform orientations and partial visual fields on images. These features of CXRs make them more difficult to interpret than CTs. Current reports suggest that radiologist diagnosis of COVID-19 from CXR has a sensitivity of 69% compared to a sensitivity of up to 97% on CT [4, 5, 7] .
There is a growing need for methods to monitor and predict disease progression in COVID-19. In severe cases of the infection, patients may progress to hypoxemic respiratory failure and acute respiratory distress syndrome requiring mechanical ventilation [9] . In the United States, early projections forecasted a deficit in ventilators during the COVID-19 pandemic, indicating the importance of efficient resource management [10] . The ability to identify patients that might progress to critical illness beginning at clinical presentation will be invaluable in potential ventilator shortages, and a few studies have demonstrated that radiologic imaging may be of use in this regard [3, 6, [11] [12] [13] . Recent studies have qualitatively described the association of groundglass opacities and lung consolidations with disease presence and progression on CXR and CT [2-4, 6, 7] . Specifically, the presence of opacities in multiple lobes has been shown to predict severe illness [3] . Studies have also evaluated various clinical biomarkers and comorbidities as predictors of disease progression, and there is some evidence that imaging data might complement these models [6, [13] [14] [15] [16] [17] . However, uses of CXRs have been largely qualitative, and there has been little work linking quantifiable CXR findings with patient outcomes in COVID-19. Nevertheless, portable radiography is an appealing, and often the only modality in high-volume hospital settings. In this study we utilize computational techniques to further evaluate the role of CXR in predicting patient outcomes.
Computational radiology is a rapidly advancing field that employs machine learning to interpret medical images. Two general approaches include deep learning (DL) and radiomic analysis [12, 18] . Recent studies have used these techniques in order to study COVID-19, but few have applied them to multi-institutional CXR cohorts [2, 12, 19, 20] . Furthermore, there has not been extensive work utilizing computational models to predict COVID-19 outcomes using CXR.
In this study, we use computational techniques to identify clinically actionable information from CXR. We first predict both mechanical ventilation requirement and mortality from baseline CXR using a DL model. Our second experiment uses machine learning classifiers to predict these same outcomes based upon computer-extracted, pre-defined radiomic features. Third, we propose a combined DL model using both processed CXRs and corresponding radiomic-embedded feature maps as inputs to predict outcomes. This synergistic approach utilizing radiomic-embedded maps for DL has not been explored in evaluating CXRs and may offer insights into novel interpretations of pre-defined radiomic features. Figure 1 displays a general flowchart of experiments. Figure 1 : Study pipeline. Visualized here is the schema for the experiments performed in this study. Experiment 1 uses a CNN deep learning model to predict COVID-19 patient outcomes using segmented CXRs as inputs. In Experiment 2 we extract pre-defined radiomic features from segmented CXRs and input them into machine learning models such as Random Forests. In Experiment 3 we use extracted radiomic features to generate radiomic-embedded maps which are inputted with segmented CXRs into a CNN deep learning model.
In this two-center, IRB approved study, anonymized coronal CXRs were obtained from patients suspected of COVID-19 on presentation at Stony Brook University Hospital (SBUH) and Newark Beth Israel Medical Center (NBIMC) between March and June 2020 (Figure 2 ). At SBUH, 484 baseline CXRs for 463 patients were analyzed. Among these, 17 CXRs from 16 patients were discarded due to being images of pediatric patients or due to poor image quality. Here, a baseline CXR refers to any CXR taken on the first day for which CXR data exists for a patient. A total of 72 baseline CXRs obtained from 72 patients prior to ventilation were included from NBIMC. Of these, 5 CXRs were discarded due to indistinguishable lung fields.
In total, 534 CXRs taken from 514 patients were analyzed in this study. 305 of these images were from 290 male patients and 229 were from 224 female patients. The mean age of patients studied was 54 years old (median=55 years, standard deviation=18.055 years, Table 1 ). COVID-19 positivity was tested for each patient via reverse transcriptase polymerase chain reaction (RT-PCR). 353 CXRs were taken from 338 patients who tested positive for COVID-19 (Table 2 ) and 181 CXRs were from 176 patients who were found to be negative. CXRs taken from COVID-19 positive patients were used in outcome prediction experiments whereas those from both COVID-19 positive and negative patients were used to build lung and artifact segmentation models. Of the 353 CXRs from positive patients, 113 baseline CXRs were taken for 109 patients that later required mechanical ventilation. 89 CXRs were from 84 patients who later died from the disease. Representative CXR images are displayed in Figure 3 . Training of all machine learning models was performed on CXRs acquired at SBUH and evaluated using a combination of CXRs taken at both SBUH and NBIMC. 
To conduct our analysis, we ensured that CXRs were properly segmented to avoid analysis of features unrelated to lung fields. In order to segment the lungs from CXR images, a Residual U-Net DL model was employed [21, 22] . This architecture was augmented using a multiscale image input pyramid for better intermediate feature representations with deep supervision (Supplementary Figure S1 ) [23] . To train the network, lung fields were first manually segmented for a dataset of 100 CXRs, excluding heart shadows. Additionally, artifacts such as EKG leads, pacemakers, and other non-anatomical objects were manually segmented in this training set. These segmentations were used to train a second multiscale-input Residual U-Net model to generate artifact masks. A focal Tversky loss function was employed (alpha=0.3, gamma=1.0) to ensure a higher penalty to false positive results [24] . This was to avoid misidentification of high-intensity objects as lungs and to mitigate misclassification of lungs as unwanted artifacts. These models were then used to generate lung and artifact masks for the remaining 434 CXRs. Each of these masks was manually reviewed and errors in segmentation, if any, were corrected. Convolutional neural networks (CNNs) were employed to predict future mechanical ventilation requirement and patient mortality from the baseline CXRs. Utilizing the earlier-described image processing pipelines, lung and artifact segmentation was performed. Additional preprocessing steps included cropping of the full CXR to a tight boundary around the lungs, resizing of input images to 224 by 224 pixels, and the application of min-max normalization to rescale image intensity values between 0 and 1. For mortality classification, training was performed on 27 CXRs from COVID-19 patients who did not survive (Group NR) and 27 CXRs from patients who recovered (Group R). Validation was performed on 10 CXRs from (Group NR) and 86 CXRs from (Group R). A ResNet-18 architecture was determined to produce the best results on validation. Testing was performed on 52 CXRs from (Group NR) and 138 CXRs from (Group R). To develop machine learning classifiers to predict future mechanical ventilation requirement and mortality, baseline CXRs were divided into the same training, validation, and testing sets used in Experiment 1. First, the training set was used to train Random Forest (RF), Linear Discriminant Analysis, and Quadratic Discriminant Analysis classifiers [31, 32] . For each of 100 iterations in a 3-fold cross validation setting, feature reduction among radiomic and clinical features was performed on the training set using a Wilcoxon rank sum test, student's t-test, or a maximum relevance minimum redundancy approach [33]. Highly correlated features (identified using a correlation threshold of 0.9) were removed to reduce redundancy. Validation experiments were then performed to identify the best classifier. Finally, the chosen classifiers, trained on the training set were evaluated on the distinct testing set of CXRs.
Based upon the results of validation experiments, RF classifiers were used for both COVID-19 mechanical ventilation requirement and mortality prediction. For ventilation classification, eight radiomic features were used while six were selected for mortality classification (Table S1 ).
The radiomic feature statistics from Experiment 2 were used to create radiomic-embedded feature maps for each CXR. t-Distributed Stochastic Neighbor Embedding (random state=1) was employed to perform feature reduction and to convert radiomic data to a 2D representation [34] . The normalization procedure is detailed in Supplementary Section I.
To assess the predictive capability of a model trained using both these radiomic-embedded feature maps and CXR images as inputs, the same general procedure employed in Experiment 1 was used.
A key difference was a change in the first input convolution filter of each ResNet architecture to receive a 2-channel rather than a 3-channel input. Here our first channel was the CXR and the second was the corresponding radiomic-embedded feature map. This new input layer did not use pre-trained weights. All other network configurations are identical to those described in Experiment 1. A ResNet-50 architecture and ResNet-34 architecture were determined by validation to perform best for mechanical ventilation requirement and mortality prediction, respectively. Dataset splits of each of these classifiers were identical to those detailed in Experiment 1.
Evaluation of all classifiers in Experiments 1 through 3 was performed using a bootstrap resampling approach on the test set over 100 iterations.
Class Activation Maps (CAMs) were also generated for Experiment 1 using network outputs prior to the global average pooling layer in each ResNet architecture. These CAMs enable a degree of visualization of the areas in an image that the model used to make predictions.
Results and 95% confidence intervals for Experiments 1, 2, and 3 are summarized in Table 3 and  Table 4 .
In Experiment 1, a ResNet-34 model trained to predict future mechanical ventilation requirement had an AUC of 0.842, a specificity of 73%, and a sensitivity of 83% on the testing dataset. A ResNet-18 model trained to predict mortality yielded an AUC of 0.639, a specificity of 66%, and a sensitivity of 50% on the testing dataset. Representative CAMs are shown in Figure 4 . An expert reader (J.G, 15 years of experience) noted diffuse, bilateral patchy infiltrates most notable in the middle and lower lung zones on both CXR images. These findings align with our CAM findings of high network activations within lung fields, predominantly in the middle and lower lung zones. For Experiment 2, an RF classifier, trained to predict need for mechanical ventilation, yielded an AUC of 0.905, a specificity of 63%, and a sensitivity of 97% on the testing dataset. An RF classifier was also used to predict mortality in COVID-19 positive patients and had an AUC of 0.926, a specificity of 69% and a sensitivity of 100% on the testing dataset. The top radiomic features used for predicting mechanical ventilation requirement and mortality were the Laws S5L5 filter responses and Haralick correlation, respectively ( Figure 5 ). For Experiment 3, a ResNet-50 model trained to predict future mechanical ventilation requirement had AUC of 0.903, a specificity of 71%, and a sensitivity of 91% on the testing dataset. A ResNet-34 model trained to predict mortality yielded an AUC of 0.756, a specificity of 62%, and a sensitivity of 70% on the testing dataset. For both predictions, these models had higher prediction confidence than those reported for Experiment 1 ( Figure 6 ). 
As the COVID-19 pandemic continues to unfold, there will be a growing need for useful interpretations of CXRs. In this work we have presented models for baseline CXRs that demonstrate high sensitivities in predicting future mechanical ventilation requirement (≤97%) and mortality (≤100%). For these predictions of disease outcomes, radiomic classifiers had higher sensitivities than those yielded by DL. However, we did report higher specificities for prediction of mechanical ventilation requirement using DL classifiers with radiomic embedding (Table S2) . Several non-imaging models have been proposed with high sensitivities for various clinical outcomes using biomarkers such as serum lactate dehydrogenase, C-reactive protein, lymphocyte counts, and coagulation factors [13] [14] [15] [16] [17] . These models might be complemented by novel imagingbased approaches. The quantitative modeling of baseline CXRs is relatively unexplored and will be of particular importance for resource allocation if projected scarcity of ventilators and other supplies occurs [10] .
Previous studies have applied DL to the analysis of COVID-19 CXRs [12, 19, 20, 35] . However, at least one study has reported potential deficiencies in these approaches, including insufficiencies in a commonly used public dataset, a neglect to segment lung fields, and a failure to account for large differences between disparate public datasets [35] . Previous studies have also not explicitly accounted for foreign objects in the lung fields such as buttons, clips, wires, etc., which can obscure pathological findings. Here, we have developed a unique pre-processing pipeline to segment lung fields from CXRs. We have further shown that our classifiers are effective in predicting COVID-19 outcomes in patients treated at multiple institutions, reflecting a degree of robustness in their predictive value.
The high middle and lower lung activations observed in network CAMs are consistent with studies reporting lower lobe lung findings on both CT and CXR (Figure 4 ) [4] [5] [6] . In the CAMs presented in Figure 4 , high network activations are visualized in areas where a radiologist had independently identified lung opacities, providing evidence that our DL models are using pathology-relevant features to predict patient outcomes. Using an integrated model that incorporates both baseline imaging and the associated radiomicembedded feature maps, we demonstrated an improvement in outcome prediction results over use of baseline image-based DL models alone ( Figure 6 ). The generation of a radiomic-embedded feature map for DL classification is a novel approach in CXR analysis and may be used in future studies to augment DL approaches. While most learning approaches assume independence between prediction variables, DL considers the 'relative arrangement' of primitives such as that of image pixels. To increase the versatility of CNNs, we first embedded radiomic statistics in a reduced dimensional space, clustered similar features together, and then used the embedding to obtain a pseudo image. This pseudo image was inputted into our classifiers in addition to the original CXR scan. This fusion framework can help assess feature importance and may also be extended to other clinical parameters. Both radiomic and radiomic embedding approaches can give insight into what features of a patient's CXR are significant in making predictions, and the use of either can be more informative to a physician than exclusively DL approaches.
There are certain limitations in our work. First, we used baseline CXRs that are likely to be nonuniform in the interval between COVID-19 infection and image acquisition. While our data is representative of the clinical reality of patients receiving baseline CXRs at varying timepoints in the course of their disease, future studies might build improved time-to-event prediction models using data with a more uniform temporal distribution. Furthermore, we are limited in the number of clinical features studied and our models might benefit from including co-morbidities such as a history of cancer, chronic obstructive pulmonary disease, hypertension, etc. The inclusion of comorbidities and other laboratory variables might further improve our models [13] [14] [15] [16] [17] . Finally, future work will be helpful in demonstrating the robustness of these classification models in the broader context of COVID-19 treatment in other hospitals and locations.
In summary, we have presented a complete pipeline for computational evaluation of CXR in COVID-19 patients. Both radiomic and DL classification models enable us to predict mechanical ventilation requirement and mortality from baseline CXRs. Furthermore, we demonstrate the improvement that a novel radiomic embedding approach has on DL predictions of COVID-19 outcomes. We posit that the ability to make early predictions of disease outcomes may aid in triage, clinical decision-making, and efficient hospital resource allocation as the COVID-19 pandemic progresses. 