Dataset collection and preparation.
We extracted a dataset of 19,957 chest X-ray exams collected from 4,772 patients who tested positive for COVID-19 between March 2, 2020, and May 13, 2020. We applied inclusion and exclusion criteria that were defined in collaboration with clinical experts, as shown in Figure 2.b. Specifically, we excluded 783 exams that were not linked to any radiology report, nine exams that were not linked to any encounter information, and 5,213 exams from patients who were still hospitalised by May 13, 2020. To ensure that our system predicts deterioration prior to its occurrence, we excluded 6,260 exams that were collected after an adverse event and 187 exams of already intubated patients. The final dataset consists of 7,502 chest X-ray exams corresponding to 4,204 unique patients. We split the dataset at the patient level such that exams from the same patient exclusively appear either in the training or test set. In the training set, we included exams that were collected both in the ED and during inpatient encounters. Since the intended clinical use of our model is in the ED, the test set only includes exams collected in the ED. This resulted in 5,224 exams (5,617 images) in the training set and 770 exams (832 images) in the test set. We included both frontal and lateral images, however there were less than 50 lateral images in the entire dataset.
The data used to evaluate the models during deployment consist of 375 exams from 217 patients collected between May 22, 2020 and June 24, 2020. The exams were filtered based on the same criteria described above. Among the 375 exams, 25 chest X-ray exams were collected from patients who were admitted to the ICU within 96 hours, and six exams were collected from patients who were intubated within 96 hours.
After extracting the images from DICOM files, we applied the following preprocessing procedure. We first thresholded and normalized pixel values, and then cropped the images to remove any zero-valued pixels surrounding the image. Then, we unified the dimensions of all images by cropping the images outside the center and rescaling. We performed data augmentation by applying random horizontal flipping (p=0.5), random rotation (-45 to 45 degrees), and random translation. Supplementary Figure 1 shows the distribution of the size of the images prior to data augmentation, as well as examples of images before and after preprocessing.
In addition to the chest X-ray images, we extracted clinical variables for patients including patient demographics (age, weight, and body mass index), vital signs (heart rate, systolic blood pressure, diastolic blood pressure, temperature, and respiratory rate), and 25 lab test variables listed in Supplementary Table 1. All vital signs were collected prior to the chest X-ray exam.
Adverse event prediction.
Our main goal is to predict clinical deterioration within four time windows of 24, 48, 72, and 96 hours.
We frame this as a classification task with binary labels y=[y24,y48,y72,y96] indicating clinical deterioration of a patient within the four time windows. The probability of deterioration is estimated using two types of data associated with the patient: a chest X-ray image, and routine clinical variables. We use two different machine learning models for this task: COVID-GMIC to process chest X-ray images, and COVID-GBM to process clinical variables. For each time window t∈Ta={24,48,72,96}, both models produce probability estimates of clinical deterioration, ^ytCOVID-GMIC,^ytCOVID-GBM∈[0,1].
In order to combine the predictions from COVID-GMIC and COVID-GBM, we employ the technique of model ensembling dietterich2000ensemble . Specifically, for each example, we compute a multi-modal prediction ^yENSEMBLE as a linear combination of ^yCOVID-GMIC and ^yCOVID-GBM:
|
^yENSEMBLE=λ^yCOVID-GMIC+(1−λ)^yCOVID-GBM, |
|
(1) |
where λ∈[0,1] is a hyperparameter. We selected the best λ by optimizing the average of the AUC and PR AUC on the validation set. In Supplementary Figure 2.b, we show the validation performance of ^yENSEMBLE for varying λ.
Chest X-ray image model.
We process chest X-ray images using a deep convolutional neural network model, which we call COVID-GMIC, based on the GMIC model shen2019globally ; shen2020interpretable . COVID-GMIC has two desirable properties. First, COVID-GMIC generates interpretable saliency maps that highlight regions in the X-ray images that correlate with clinical deterioration. Second, it possesses a local module that is able to utilize high-resolution information in a memory-efficient manner. This avoids aggressive downsampling of the input image, a technique that is commonly used on natural images he2016deep ; huang2017densely , which may distort and blur informative visual patterns in chest X-ray images such as basilar opacities and pulmonary consolidation. In Supplementary Table 2, we demonstrate that COVID-GMIC achieves comparable results to DenseNet-121, a neural network model that is not interpretable by design, but is commonly used for chest X-ray analysis rajpurkar2017chexnet ; allaouzi2019novel ; liu2019sdfn ; guan2020multi .
The architecture of COVID-GMIC is schematically depicted in Figure 1.b. COVID-GMIC processes an X-ray image x∈RH,W (H and W denote the height and width) in three steps. First, the global module helps COVID-GMIC learn an overall view of the X-ray image. Within this module, COVID-GMIC utilizes a global network fg to extract feature maps hg∈Rh,w,n, where h, w, and n denote the height, width, and number of channels of the feature maps. The resolution of the feature maps is chosen to be coarser than the resolution of the input image. For each time window t∈Ta, we apply a 1×1 convolution layer with sigmoid activation to transform hg into a saliency map At∈Rh,w that highlights regions on the X-ray image which correlate with clinical deterioration. Each element Ati,j∈[0,1] represents the contribution of the spatial location (i,j) in predicting the onset of adverse events within time window t. In order to train fg, we use an aggregation function fagg:Rh,w↦[0,1] to transform all saliency maps At for all time windows t into classification predictions ^yglobal:
|
fagg(At)=1|H+|∑(i,j)∈H+Ati,j, |
|
(2) |
where H+ denotes the set containing the locations of the r% largest values in At, and r is a hyperparameter.
The local module enables COVID-GMIC to selectively focus on a small set of informative regions. As shown in Figure 1, COVID-GMIC utilizes the saliency maps, which contain the approximate locations of informative regions, to retrieve six image patches from the input X-ray image, which we call region-of-interest (ROI) patches. Figure 3 shows some examples of ROI patches. To utilize high-resolution information within each ROI patch ~x∈R224,224, COVID-GMIC applies a local network fl, parameterized as a ResNet-18 he2016deep
, which produces a feature vector
~hk∈R512 from each ROI patch. The predictive value of each ROI patch might vary significantly. Therefore, we utilize the gated attention mechanism
ilse2018attention to compute an attention score
αk∈[0,1] that indicates the relevance of each ROI patch
~x for the prediction task. To aggregate information from all ROI patches, we compute an attention-weighted representation:
The representation z is then passed into a fully connected layer with sigmoid activation to generate a prediction ^ylocal. We refer the readers to Shen et al. shen2020interpretable for further details.
The fusion module combines both global and local information to compute a final prediction. We apply global max pooling to
hg, and concatenate it with
z to combine information from both saliency maps and ROI patches. The concatenated representation is then fed into a fully connected layer with sigmoid activation to produce the final prediction
^yfusion.
In our experiments, we chose H=W=1024. Supplementary Table 2 shows that COVID-GMIC achieves the best validation performance for this resolution.
We parameterize fg as a ResNet-18 he2016deep which yields feature maps hg with resolution h=w=32, and number of channels n=512
. During training, we optimize the loss function:
|
l(y,^yglobal,^ylocal,^yfusion)=1|Ta|∑t∈TaBCE(yt,^ytglobal)+%
BCE(yt,^ytlocal)+BCE(yt,^ytfusion)+β|At|, |
|
(4) |
where BCE denotes binary cross-entropy and β is a hyperparameter representing the relative weights on an ℓ1-norm regularization term that promotes sparsity of the saliency maps. During inference, we use ^yfusion as the final prediction generated by the model.
Estimation of deterioration risk curves.
The deterioration risk curve (DRC) represents the evolution of the deterioration risk over time for each patient. Let T denote the time of the first adverse event. The DRC is defined as a discretized curve that equals the probability P(T≤ti) of the first adverse event occurring before time ti∈{ti|1≤i≤8}, where t1=3, t2=12, t3=24, t4=48, t5=72, t6=96, t7=144, t8=192 (all times are in hours).
Following recent work on survival analysis via deep learning gensheimer2018scalable , we parameterize the DRC using a vector of conditional probabilities ^p∈R8. The ith entry of this vector, ^pi, is equal to the conditional probability of the adverse event happening before time ti given that no adverse event occurred before time ti−1, that is:
|
^pi={P(T≤t1),i=1,P(T≤ti|T>ti−1),2≤i≤8. |
|
(5) |
Given an estimate of ^p
, the DRC can be computed applying the chain rule:
|
DRC(ti)=P(T≤ti)=1−P(T>ti)=1−i∏j=1P(T>tj|T>tj−1)=1−i∏j=1(1−^pj). |
|
(6) |
We use the GMIC model to estimate the conditional probabilities ^p from chest X-ray images. We refer to this model as COVID-GMIC-DRC. As explained in the previous section, the GMIC model has three different outputs corresponding to the global module, local module and fusion module. When estimating conditional probabilities for the eight time intervals, we denote these outputs by ^pglobal, ^plocal, and ^pfusion. During inference, we use the output of the fusion module, ^pfusion, as the final prediction of the conditional-probability vector ^p. We use an input resolution of H=W=512 and parameterize fg as ResNet-34 he2016deep . The resulting feature maps hg have resolution h=w=16 and number of channels n=512. The results of an ablation study that evaluates the impact of input resolution and compares COVID-GMIC-DRC to a model based on the Densenet-121 architecture, are shown in the Supplementary Tables 2 and 5.
During training, we minimize the following loss function defined on a single example:
|
l(T,^pglobal,^plocal,^p%
fusion)=ls(T,^pglobal)+ls(T,^plocal)+ls(T,^pfusion)+8∑m=0β|Am|, |
|
(7) |
where ls is the negative log-likelihood of the conditional probabilities. For a patient who had an adverse event between ti−1 and ti (where t0=0), this negative log-likelihood is given by
|
ls(T,^p)=−lnP(ti−1≤T≤ti)=−lni−1∏j=1P(T>tj|T>tj−1)P(T≤tj|T>ti−1)=−i−1∑j=1ln(1−^pj)−ln^pi. |
|
(8) |
The framework can easily incorporate censored data corresponding to patients whose information is not available after a certain point. The negative log-likelihood corresponding to a patient, who has no information after ti and no adverse events before ti, equals
|
ls(T,^p)=−lnP(T>ti)=−lni∏j=1P(T>tj|T>tj−1)=−i∑j=1ln(1−^pj). |
|
(9) |
Note that each ^pi is estimated only using patients that have data available up to ti
. The total negative log-likelihood of the training set is equal to the sum of the individual negative log-likelihoods corresponding to each patient, which makes it possible to perform minimization efficiently via stochastic gradient descent. In contrast, deep learning models for survival analysis based on Cox proportional hazards regression
cox1984analysis require using the whole dataset to perform model updates
ching2018cox ; katzman2018deepsurv ; liang2020early , which is computationally infeasible when processing large image datasets.
Model training and selection.
In this section, we discuss the experimental setup used for COVID-GMIC, COVID-GMIC-DRC, and COVID-GBM. The chest X-ray image models were implemented in PyTorch
NEURIPS2019_9015 and trained using NVIDIA Tesla V100 GPUs. The clinical variables models were implemented using the Python library LightGBM
ke2017lightgbm .
We initialized the weights of COVID-GMIC and COVID-GMIC-DRC by pretraining them on the ChestX-ray14 dataset wang2017chestx (Supplementary Table 3 compares the performance of different initialization strategies). We used Adam kingma2014adam with a minibatch size of eight to train the models on our data. We applied data augmentation during training and testing, but not during validation. During testing, we augmented each image ten times and averaged the corresponding outputs to produce the final prediction.
We optimized the hyperparameters using random search bergstra2012random . For COVID-GMIC, we searched for the learning rate η∈10[−6,−4] on a logarithmic scale, the regularization hyperparameter β∈4×10[−6,−3] on a logarithmic scale, and the pooling threshold r∈[0.2,0.8] on a linear scale. For COVID-GMIC-DRC, based on the preliminary experiments, we fixed the learning rate to 1.25 × 10−4. We searched for the regularization hyperparameter, β∈10[−6,−4] on a logarithmic scale, and the pooling threshold r∈{0.2,0.5,0.8}. For COVID-GBM, we searched for the learning rate η∈10[−2,−1] on a logarithmic scale, the number of estimators e∈10[2,3] on a logarithmic scale, and the number of leaves l∈[5,15] on a linear scale. For each hyperparameter configuration, we performed Monte Carlo cross-validation xu2001monte (we sampled 80% of the data for training and 20% of the data was used for validation). We performed cross-validation using three different random splits for each hyperparameter configuration. We then selected the top three hyperparameter configurations based on the average validation performance across the three splits. Finally, we combined the nine models from the top three hyperparameter configurations by averaging their predictions on the held-out test set to evaluate the performance. This procedure is formally described in Supplementary Algorithm 1.
share
