Abstract-X-ray images may present non-trivial features with predictive information of patients that develop severe symptoms of COVID-19. If true, this hypothesis may have practical value in allocating resources to particular patients while using a relatively inexpensive imaging technique. The difficulty of testing such a hypothesis comes from the need for large sets of labelled data, which not only need to be well-annotated but also should contemplate the post-imaging severity outcome. On this account, this paper presents a methodology for extracting features from a limited data set with outcome label (patient required ICU admission or not) and correlating its significance to an additional, larger data set with hundreds of images. The methodology employs a neural network trained to recognise lung pathologies to extract the semantic features, which are then analysed with a shallow decision tree to limit overfitting while increasing interpretability. This analysis points out that only a few features explain most of the variance between patients that developed severe symptoms. When applied to an unrelated, larger data set with labels extracted from clinical notes, the method classified distinct sets of samples where there was a much higher frequency of labels such as Consolidation, Effusion, and alveolar. A further brief analysis on the locations of such labels also showed a significant increase in the frequency of words like bilateral, middle, and lower in patients classified as with higher chances of going severe. The methodology for dealing with the lack of specific ICU label data while attesting correlations with a data set containing text notes is novel; its results suggest that some pathologies should receive higher weights when assessing disease severity.
Index Terms-Covid-19, deep learning, ICU, severity, X-ray.
C OVID-19 remains a serious worldwide health pandemic with infections numbering 21,294,845 as of 16 August 2020 with 761,779 deaths [1] . One emerging characteristic of COVID-19 disease is the wide range of symptoms experienced by infected persons ranging from entirely asymptomatic through admission to the general ward for a range of symptoms including fever, cough, fatigue, headache and diarrhoea [2] to severe pneumonia requiring admission to ICU with mechanical ventilation [3] . Reported case-fatality rates vary from 1% to greater than 7%, usually due to respiratory failure [4] .
Whilst care in modern ICUs will result in death rates towards the lower end of the range, life-sustaining therapies will, in practice, be limited by lack of personnel, infrastructure and materials and equipment. This limitation of resources leads clinicians to make prognostic decisions based on criterion such as old age, fragility and comorbidity that can lead to the death of patients with poor prognosis in favour of patients with better prognosis [4] .
The most widely used testing methodology for COVID-19 is the real-time reverse-transcriptase-polymerase-chain-reaction (RT-PCR) test [2] , [5] . This test provides a binary indication of whether a person is infected with COVID-19 or not. Similarly, serology point of care tests that detect COVID-19 antibodies provide a retrospective measure of infection [5] but do not provide information relating to the severity or progression of the disease and are not useful in helping clinicians to devise an individual plan of treatment for a patient.
Chest medical imaging has proven to be useful in managing more serious COVID-19 infections since respiratory dysfunction is the primary source of COVID-19 morbidity and mortality. Several researchers have shown Deep Learning to be useful in classifying COVID-19 cases from medical images including CXR [6] - [9] , CT [10]- [12] with some research also achieving promising results with the Ultrasound imaging mode [13] , [14] . However, using Deep Learning as a diagnostic tool can be problematic as it is hard to assess biases, risk, potential overfitting, and ability to generalize in clinical settings [15] , [16] . Its value resides more on the prognosis and treatment side than in actual diagnostic use [17] . Chest X-rays (CXR) and Computed Tomography (CT) imaging are useful tools in the management of moderate to severe COVID-19 cases since these methods help clinicians to establish a baseline pulmonary status and identify underlying pulmonary conditions that may contribute to the patients risk. Chest X-rays, in particular, are also non-invasive and can be potentially used as bed-side patient monitoring tool [18] . Compared to CT, CXR is less expensive, more available, and require less technical expertise to perform and interpret than Ultrasound. In the case of clinical worsening, imaging can be useful for assessing COVID-19 progression and secondary abnormalities [19] .
This paper therefore focuses on the prognosis and provides a machine learning-based method for assessing features with predictive potential to classify patients who are more likely to develop severe COVID-19 symptoms and be eventually admitted in ICU. The method proposed here leverages the existence of a data set with patient outcome labels (only publicly accessible for CXR images) and the fact that such images were taken before the patients were admitted to ICU. As such a data set has a small sample size, significant measures were taken to aggressively limit potential overfitting; from image pre-processing techniques and lung segmentation to model and feature selection. Measures like lung segmentation and cropping are not only needed given the data set size but also to prevent the method to rely on confounding image features, as already pointed as a known problem in a recent work [20] . For these tasks, lung segmentation and cropping, a model achieving state-of-the-art results is also proposed and part of the contribution herein. The validation is not only attested from separability metrics using the selected semantic features in the ICU data set but also by correlating the results to a more extensive external data set. Despite not having ICU-related labels, such a data set contains pathologies and localizations labels in its metadata that allows for assessing the significance of the selected semantic features. The authors present such a machine learning system in the hopes that it will inform practitioners of the weight of particular features in patient prognosis, thereby helping early intervention to improve patient outcomes and to aid the allocation of critical resources. A second significant benefit of this novel framework is that it can be used as a method to assist the medical community in quickly identifying features in medical images that are associated with severe progression of COVID-19 and other novel infections.
Given the limited amount of data to learn the severitycorrelating features, a focus on limiting potential overfitting was central to most methodological decisions. Instead of disregarding such investigation simply because the data is limited, the idea proposed here is to use low-complexity, interpretable methods to test features for potential predictive value. Therefore, in an attempt to increase the features generalization potential, a conservative bias for complexity was adopted when choosing models and hyper-parameters. Some of this concern for overfitting manifested as procedures adopted in the pre-processing of images: histogram equalizations (adaptive and standard) and lung segmentation and cropping. In the learning stage, more importantly, it inspired the choice to use a low complexity (low VapnikChervonenkis dimension), shallow decision tree, as to limit the ability of the method to simply shatter the data set. Such an approach also presents the considerable advantage of resulting in a highly interpretable model, which is much desired in systems designed to support human decision-making. To improve understandability, Fig.  1 is set to depict a high-level graphical summary of the methodology.
The authors nevertheless understand that, despite the care taken for overfitting and model selection, it is highly desirable to have additional validation methods that can attest for the model's generalization ability. To this end, this study also contemplates the application of such a model to a larger, different dataset (1000+ images). Although lacking the metadata regarding ICU admission, such X-ray images have clinical annotations that can be tokenized and correlated to the semantic features learned by the low-complexity model. Given the features semantic nature, one can also correlate them to other works in the related literature that report on the lung pathologies present in patients that developed severe symptoms.
The data set used to learn the features present in patients with higher chances of being admitted to ICU was a subset of the data presented by Cohen et al. [21] . It is one of the most popular data sets on the literature, favourited more than 2000 times on Github. One of its most positive characteristics is its rich metadata containing categories such as sex, age, location, patient condition, and outcome (ICU admission), allowing for the investigation presented here. However, it should be noted that many of such interesting labels are sparsely distributed in the data set and not present for every every image. For the images that had rich descriptors, two prominent labels were of interest when setting up this ICU data set: went-icu and in-icu. An image taken from a patient marked with Y on the former label and N on the latter is a sample from a patient that eventually developed severe symptoms before they were admitted in ICU. Therefore, one can reasonably form the hypothesis that there might be features in these images that are associated with patients that were eventually admitted in ICU. In total, 100 images contained these two labels were used in the analysis, 65 of them labelled as not-icu (0) and 35 as icu (1).
The selected images are assigned to a class (ICU or not) and fed through a pre-processing pipeline established to normalize and remove potential bias-inducing artifacts. The sequential blocks, illustrated in Fig. 1 , are responsible for histogram equalizations, lung segmentation and lung-area cropping. Histogram equalization is set to impose a normalizing effect on the contrast of images by equalizing the distribution of pixel intensities that might be concentrated in a narrow range. Such a procedure usually improves the visual contrast in the X-ray images, which can come from particular distribution depending on the equipment manufacturer.
Subsequent to the normalization step, a lung-segmentation model is applied to each X-ray image to improve the signalto-noise ratio for machine learning classification. Such a step is potentially the most important regarding the mitigation of potential bias as all areas that are not lung-related are excluded, allowing the model the focus on lung-specific features. Given that COVID-19 is a respiratory disease primarily affecting the lung, non-lung field pixels are considered here to be noise and therefore segmented from the lung field. The utility of this approach in medical image classification is supported in past studies [22] , and recent works have pointed the fact that that lack of such practice can make the model focus on confounding, not disease-related features [20] . A comprehensive review of lung area segmentation techniques may be found in [23] , which notes that Deep Learning methods provide similar accuracy as in inter-observer performance at the expense of a long training process.
With the assumption that one could improve the U-Net architecture to yield better dice similarity co-efficient results, as well as using training data from two additional data sources, a lung-segmentation model was proposed and trained. The additional images are from Montgomery and Shenzen datasets [24] and resulted in a combination corpus of 1185 CXR image/mask pairs. As typical U-net designs are based on simple convolution stacks similar to the VGG network architecture [25] , a U-Net based on a skip-connection architecture following [26] , [27] is also proposed. Such an architecture allows the network to re-use features learned in earlier layers, thereby improving its ability to learn the identity function and increase generalization. As an additional benefit, this architecture also allows the U-Net to train on a smaller number of epochs, thereby reducing the training time. Given the state-of-the-art dice similarity coefficient achieved by the lung-segmentation model proposed here, validated in three independent data corpora, it is expected that the model will be useful in unseen images. Artefacts in the generated lung field masks were removed by a combination of morphological closing, contour filling and flood-filling, resulting in a set of automatically segmented lung-field images. As a final step, the images were also automatically cropped to fully contain the segmented lung field, resulting in a uniform image set to the downstream classifier and normalizing for different lung sizes. In examples of the result of applying the pre-processing pipeline, Fig. 2 illustrates three images in the input and output stage, with and without normalization.
C. Feature extraction 1) TorchXRayVision: The feature extraction model is a DenseNet [28] pre-trained on 80,000+ lung X-ray images. The model, part of the TorchXRayVision library, was firstly presented in [29] and also later used in a work proposing a COVID-19 pneumonia severity score [18] . A particular aspect of this model is that the author trained it to detect specific lung pathologies from large data sets such as the CheXpert (Stanford) [30] , ChestX-ray8 (NIH) [31] and MIMIC-CXR (MIT) [32] . The semantic labelling was performed by adding an 18-node layer at the end of the network and training it to classify different pathology labels through a sigmoid activation layer. The labels that each of the nodes was trained to classify were the ones found in the large data sets used: Feature extraction from such a model is feasible from many parts of its architecture. In a way, any activation value from the forward pass through the network can be used as features. In this work, specifically, three modes were useful to attest separability between classes, assess the semantic meaning, and aid the pathology location. These methods were labelled mid-layer, last-layer, and class-activation map features. The first method, extracting features from an intermediate layer, extracts the activation values from a layer previous to last, which contains 1024 nodes. These features were helpful to attest the separability between classes but, as they are not set to translate any semantic meaning, they were not used to infer any correlation with particular pathologies. The last-layer features were extracted from the activations of the 18-node layer, pre-sigmoid, and were useful at correlating some of the mentioned pathologies to the defined classes (ICU admission or not). The class-activation map features, or gradient features for short, were hand-engineered features soon to be described, set with the intention of aid the localization of the pathologies in the lung. Fig. 3 is set to illustrate potential feature-extraction branches in the network architecture. It is worth mentioning that none of these features was combined or even used in their totality; they are presented as different classification scenarios where just a few are used, complying with the discussed initial intention to limit potential overfitting.
2) Gradient features: The mid-layer and last-layer features, despite sufficient when attesting separability, do not translate any information regarding the location of the pathologies in the lung. Therefore, gradient features were developed as a way to aid the correlation of the location of the lung pathologies with the patients with higher chances of developing severe symptoms. The calculation of the gradients in the network graph leaves in respect to an input image was performed by the autograd torch class, which is often used to create saliency maps and associate explainability to models [18] , [33] , [34] . In this study, however, a procedure to create features from the gradient maps is proposed. First, the maps are construed with the energy of the accumulated gradients and then sectioned in two cuts: longitudinal and transversal. For each of these cuts, an entropy measure is calculated. The resulting feature is proposed here as a measure of the spread of the activations in their respective cuts. Such a single value measure is inspired by the Shannon entropy [35] can be trivially calculated with Eq. 1. Examples of such features and how the images were sectioned for spatial correlation are illustrated in Fig. 4 .
Given the advantages of interpretability and the commitment to limit the model complexity as primary factors in model selection, it is hard to argue for the application of another method than decision trees [36] . Besides being human interpretable, one can easily limit the effects of overfitting by mechanisms of pruning, such as setting the minimum amount of samples in the tree leaves. For the experiments performed here, for example, trees were pruned to contain at least 10 samples in each leaf, which represent one-tenth of the data. This constraint resulted in very shallow trees that only relied on a few features to make the classification (usually three to five). Such a pruning mechanism was the only hyper-parameter set when fitting the trees; the rest were all default values from the scikit learn package [37] . It is worth noting here that the goal, despite trying to be as accurate as possible, is not claiming the ability to precisely predict the patients that will be admitted to ICU but discover potential semantic predictors that are correlated with such patients.
The results are presented in three distinct scenarios with respective feature sets that the tree-fitting algorithms could greedily choose from mid-layer, last-layer, and gradient features. In each of these scenarios, the results are presented by fitting the tree in the whole data set and by performing cross-validation. To the latter, given the size of the data set, a simple 5-fold would probably not be representative of the predictive potential (especially when limited to leaf size of 10 samples). The number of folds was given by a leave-two-out approach where all the data set is used for training, except for 2 samples (one of each class), which are then used for testing. This procedure is repeated until all the samples in the dataset are used for testing. The shallow trees resulting from fitting the whole data set are also illustrated to depict their simple structure and interpretability.
In an attempt to add to the validation and address potential generalization concerns, a larger, external data set was also used in testing the fitted decision trees. Such a data set could only be used to attest correlations since its metadata does not contain labels regarding patient outcome (went to ICU or not). The absence of these labels means that it could not be used for validation (in the strict sense), but the information it contains could be correlated to results from the learned semantic features. This data set of X-ray images from COVID-19 patients, named BIMCV [38] and publicly available, contains doctors annotations for each subject that were used here as discrete labels and descriptive terms. Besides numerous, such labels refer not only to the pathologies observed by the doctors in the patients lungs but also points to their location. Therefore, because they have a semantic meaning, they can then be correlated to the features from the last layer, which are also descriptive.
The method in which the features are correlated to the external data set is also presented as a contribution and, to the best of the authors' knowledge, not previously presented. As the quality of images in the data set varies a lot, the procedure starts by running the images through the same pre-processing pipeline described before with histogram equalization, lung segmentation and cropping, and feature extraction. Some images that did not have the lungs successfully segmented or were very low-quality are filtered out, but a total of 1312 images were labelled as valid. Subsequent to the processing pipeline, the decision tree performs the classification of features in the binary classes 0 and 1, and the samples pertaining to each class are assigned to two different sets. At this point, the correlation method enters the picture by creating a dictionary with the pathologies and locations labels and their respective frequency for each group. Such a frequency dictionary is then normalized by the set size and compared between classes. The correlation is a measure of how the features relevant to the classification in the first data set are represented (frequency) in the different classes assigned in the second data set. The null hypotheses, in this case, is that if the images were sampled at random, the frequency of labels would be equal in both sets. Such a procedure will become clearer in the presentation of the results in the following section.
III. RESULTS Prior to discussing the classification metrics and feature work, it is worth noting the outcomes of training the lung segmentation network, also part of the contribution herein. Trained on large chest X-rays image data sets, the designed architecture achieved a maximum validation dice similarity coefficient of 0.988 at epoch 93, which is comparable to the best result encountered thus far. The best performance observed in the literature was achieved using a complex CNN-based Deep Learning system, achieving a dice similarity coefficient of 0.980 [39] . A much simpler approach leveraging a U-Net architecture [40] , similar to the the one here, was trained on the JSRT dataset [41] consisting of 385 CXR images with gold-standard masks to achieve a dice similarity co-efficient of 0.974 [42] . Table I is set to contrast the performance of previously methods and the one proposed here.
The separability metrics were calculated for the final semantic features and other two different scenarios. The first refers to the non-task-specific features in the intermediate layer, which although contemplating 1024 nodes, had only 4 features used to attest invariance between classes. The second scenario regards the separability of the actual semantic (pathology labels) features output by the 18-node last layer of the network. The third, gradient features, will have its discussion delayed to a further sub-section regarding the pathology localizations, which is where its value really lies. The metrics presented as whole set refers to using all samples to fit the tree while cross-val indicates the ones resulting from the leave-two-out cross-validation approach. Such an approach means fitting the same number of trees as in the number of samples in the smallest class. For example, if one class has n samples and another has m, n decision trees will be fitted if n ¡ m.
The resulting metrics showed that the mid-layer features are better at separating the classes but not for a large margin when compared to the semantic features from the last layer. While the mid-layer whole set features resulted in 0.85 accuracy, the last-layer features' accuracy was 0.82. Such metrics were obtained when limiting the minimum leaf size to 15 samples in the mid-layer scenario and 10 samples in the other; both values were arbitrarily set by trial and error attempts guided by crossvalidation. The decision trees (whole set) fitted on the midlayer and last-layer scenarios only had 4 and 5 nodes, respectively. The resulting trees, along with the confusion matrices for the cross-val metrics, are illustrated in Fig. 5 . The wholeset tree is the result to be highlighted here as it has allegedly semantic features relating to the lung pathologies. The features corresponding to the 5 splits in the tree with the pathological features are in hierarchical order: 'Effusion', 'Consolidation', 'Pneumonia', 'Lung Lesion' and 'Pleural thickening'. Such labels will be important in the following subsection when correlating with the external dataset and in the discussion section where they will be qualitative compared to the pathologies referred to in the literature regarding patients that develop severe symptoms. As for the cross-val scenario, mid-layers and last-layer features resulted in 0.78 and 0.73 accuracies, respectively. The distribution of predicted versus real labels can be inspected by the decision matrices presented in Fig.  5 . As an additional worth-mentioning note from preliminary experiments, it was observed that such accuracies would drop to values close to 0.5 if the normalization techniques were not applied, which was significant evidence of their importance.
The many labels contained in the BIMCV data set allowed for an interesting comparison between the classes after classification was applied to the images. From the 1312 images selected after filtering, 73 did not have text labels assigned, but most had at least a few pathology and location labels. After the images went through the processing and feature extraction pipeline, the fitted decision tree (whole set, last layer) classified the features and their respective text labels were assigned to sets given by their predicted class. In this manner, two dictionaries could be construed with the keys being the unique words in the set, and the values representing the frequency of such labels. The five most frequent pathological labels in the data set can be listed as examples of the types of terms describes: pattern (467), COVID (443), pneumonia (301), increased density (368) and unchanged (272).
The comparison presented here was set to be the normalized ratio between the frequency of words in the potentially severe class (1) and not (0). By normalized ratio, it is meant that the frequency of words was multiplied by the inverse of the number of samples in the respective set before the ratio calculation. The number of samples in each set after classification was 694 images in class 0 and 618 in class 1. The ratio was calculated only for words that had the minimal arbitrary number of appearances of 20 in each class. It may be worth noting that, if the samples were picked at random, despite the selection of images before classification or the number of images in a set, such a ratio would be equal to 1, given the normalization. To the result of this comparison, presents the pathological labels and respective ratios higher than 1.2 and lower than 0.8 in descending order. Such values relate to the labels that were over-represented in the class of potentially severe patients, which correlates to the features used in the decision tree. The label named Normal appeared much less in the images classified by the tress as class 1. The same analysis of over-and under-represented labels can also be performed for the localization metadata; Table III presents the relevant labels with ratios higher than 1.2 and lower than 
While the pathological features presented in the previous subsections are expressive and correlate with the frequency of labels in the external data set, they cannot be correlated to any information regarding the locations of the pathologies. The features are scalars given by their respective layers activation function, which does not translate information about their spatial distribution. The simple method proposed here uses hand-engineered features on gradient maps to check for areas where the predictive features are prevalent. The decision trees were trained in the same manner as other features but, in this case, every feature relates to an entropy metric in one of the segments extracted from the gradient maps: left and right longitudinal, and superior and inferior transversal cuts (see Fig. 4 ). The entropy features extracted from the gradient maps segments and selected as splits in the decision tree presented some correlation with the external data set and other reported findings. Regarding the separability metrics, fitting a shallow decision tree with only three features in the whole data set resulted in 83% accuracy (separability), while with 71% accuracy score in cross-validation. Shown in Fig. 6 , the three features were Fibrosis and Consolidation, on the longitudinal left, and Infiltration, on the transversal superior.
Due to the fact that the methodology includes the segmentation and cropping of the lung area as an important normalization measure, an interesting way to visualize the spatial distribution of some features becomes possible. Such a way proposed here differs from the usual by not just plotting the energy of the gradient maps of a particular feature but averaging all the maps of a certain class and then platting that as a surface. An example is presented in Fig. 7 where the chosen feature is the one in the top of the tree of semantic features: Effusion. The two images are the average of activations of the two classes, normalized by class size. The plotted surfaces show that such a feature not only has gradients with higher energy in the images but also are more spatially distributed throughout the lung (smaller entropy).
It is worth iterating that the contribution here is not the development of a deployment-ready classifier but the presentation of pathological features with similar characteristics as ones observed in patients with higher chances of developing severe symptoms. As works in the literature have been discussing [15] - [17] , using Deep Learning methods in clinical settings is highly inadvisable, even when interpretable. Practical prediction claims would already be restrained by the classification metrics presented here, but they still cant take away the fact that some features showed interesting correlations and significance to patient severity. Moreover, the contribution here is not constrained to the predictive potential of features; it also includes the methodology developed to evaluate and select them. More than informing practitioners of evidence for setting heavier weights to the importance of some pathologies in disease development, one can hope that the method presented could be useful in for problems in similar domains. Given that computer vision-based COVID-19 diagnostic tools are not recommended for practical use due to their potential bias and risk [15] - [17] , developing systems that can support human decisions in critical resource management becomes one of the most promising uses of Machine Learning [17] . Regarding significance, it is arguable that the most significant aspect of the results is how the decision tree, fitted by a limited dataset, selected images on an external data set that had an over-representation of labels relating to disease severity. The two most over-represented features in the class of higher chances to go severe Consolidation and Bilateral are often cited as pathologies relating to severity in the literature [47] - [49] . For example, the authors in [50] reported that the evolution from ground-glass opacities to consolidation was present in some severe patients. They also pointed out that a systematic review on the subject showed that the prevalent locations were bilateral and peripheral (another over-represented term from results presented here). Another highly cited work [2] reported that most of the patients had bilateral involvements and that ICU-admitted patients showed bilateral multiple lobular and subsegmental areas of consolidation.
Nevertheless, one should also note that these pathologies were not present in all severe patients, neither in the analysis presented here and in the literature. This observation should be evident since lung pathologies are not the only factor in severity and ICU admission. Several patient conditions and comorbidities should play a role in this picture, thus representing the primary limitation of the present work. Fig. 8 summarizes and illustrates some of the worth-noting correlations found here are for improved clarity and comprehension.
The results presented in this paper point to the fact that some features have reasonable predicting potential to patients that eventually develop severe symptoms and are admitted to ICU. Extracted from a limited data set of chest x-ray images, the features only presented any predictive potential when enough care was taken to prevent overfitting. Such measures included training a state-of-the-art lung segmentation network, using a pre-trained, specialized network to extract the features, image normalization techniques, and strictly selecting an interpretable classifier of low complexity. Results from running the fitted classifier on an external data set showed that it indeed selects for the semantic features chosen when learning: consolidation, effusion, and pleural. These features not only appeared more frequently in the set of patients classified as with higher chances of going severe but a frequency analysis on all labels showed that descriptors often used to describe COVID-19 pathologies were over-represented in such patients. Other worth-noting labels included alveolar and atelectasis in the over-represented features and normal and unchanged in the under-represented ones. This observation was also attested in localizations labels, which showed overrepresentation of important terms like bilateral, peripheral, and lower. Although some original visualizations were also shown, as the average energy of class-specific activations, it is clear that this paper only shows the potential of interpretable methods to overcome data size limitations and many paths could be followed for further improvements. The authors hold the opinion that the contribution here not only lies in attesting the potential of the features to aid human decision-making in resource management but also as a template method that could be used in similar problems from other domains.