Temporal Pointwise Convolutional Networks for Length of Stay Prediction in the Intensive Care Unit
Emma Rocheteau, Pietro Liò, Stephanie Hyland
Introduction
In-patient length of stay (LoS) explains approximately 85-90% of inter-patient variation in hospital costs in the United States (Rapoport et al., 2003). Extended length of stay is associated with increased risk of contracting hospital acquired infections (Hassan et al., 2010) and mortality (Laupland et al., 2006). Hospital bed planning can help to mitigate these risks and improve patient experiences (Blom et al., 2015). This is particularly important in the intensive care unit (ICU), which has the highest operational costs in the hospital (Dahl et al., 2012) and a limited supply of specialist staff and resources.
At present, discharge date estimates are done manually by clinicians, but these rapidly become out-of-date and can be unreliable (for example Mak et al. (2012) found that the average error made by clinicians was 3.82 days). Automated systems drawing on the electronic health record (EHR) have the potential to improve forecasting accuracy using state-of-the-art models that can be updated in light of new data. This has efficiency benefits in reducing the administrative burden on clinicians, and the improved accuracy may enable more sophisticated planning strategies e.g. scheduling high-risk elective surgeries on days with more availability (Gentimis et al., 2017).
In our work, we simulate real-time predictions in retrospective data by updating the patients’ remaining ICU length of stay prediction at hourly intervals during their stay using the preceding data from the EHR (similar to Harutyunyan et al. (2019)). When designing both the architecture and pre-processing, we focus on mitigating the effects of non-random missingness due to irregular sampling, sparsity, outliers, skew, and other common biases in EHR data. Our key contributions are:
A new model – Temporal Pointwise Convolution (TPC) – which combines:
Temporal convolutional layers (van den Oord et al., 2016; Kalchbrenner et al., 2016), which capture causal dependencies across the time domain.
Pointwise convolutional layers (Lin et al., 2013), which compute higher level features from interactions in the feature domain.
Our model significantly outperforms the commonly used Long-Short Term Memory (LSTM) network (Hochreiter and Schmidhuber, 1997) and the Transformer (Vaswani et al., 2017) by margins of 18-68%.
We make a case for using the mean-squared logarithmic error (MSLE) loss function to train LoS models, as it deals more naturally with positively-skewed labels.
By adding in-hospital mortality as a side-task, we demonstrate further performance gains in the multitask setting.
We perform several investigations to improve our understanding of the model, including: an extensive ablation study of the model architecture, a post-hoc analysis of feature importances with integrated gradients (Sundararajan et al., 2017), and a visualisation to show the model reliability as a function of the time since admission and the predicted remaining LoS.
Additionally, we develop a data processing pipeline for the eICU (Pollard et al., 2018) and MIMIC-IV (Johnson et al., 2020) databases that is designed to i) mitigate some of the impact of sparsity (for the diagnoses) and missing data (for time series) in the EHR and ii) extract a wide variety of features semi-automatically such that the approach is generalisable to other EHR databases. Our code is available at: https://github.com/EmmaRocheteau/TPC-LoS-prediction.
Related Work
Despite its importance, LoS prediction has received less attention than mortality prediction. This could be due to its difficulty; LoS depends heavily on operational factors and there is considerable positive skew in its distribution (see Figure 1). While it has been addressed as a regression problem (optimised using the mean-squared error (MSE) (Purushotham et al., 2018; Sheikhalishahi et al., 2019)), it is often simplified into binary classification (short vs. long stay) (Gong et al., 2017; Nestor et al., 2018; Rajkomar et al., 2018), or as a multi-class task (Harutyunyan et al., 2019). This simplification comes at a cost of utility, so we choose to focus on the more challenging regression variant.
Owing to the centrality of time series in the EHR, LSTMs have been by far the most popular model for predicting LoS (Harutyunyan et al., 2019; Sheikhalishahi et al., 2019; Rajkomar et al., 2018). This reflects the prominence of LSTMs in other clinical prediction tasks such as predicting in-hospital adverse events including cardiac arrest (Tonekaboni et al., 2018) and acute kidney injury (Tomašev et al., 2019), forecasting diagnoses, medications and interventions (Choi et al., 2015; Lipton et al., 2015; Suresh et al., 2017), missing-data imputation (Cao et al., 2018), and mortality prediction (Che et al., 2018; Harutyunyan et al., 2019; Shickel et al., 2019). More recently, the Transformer model (Vaswani et al., 2017) been shown to marginally outperform the LSTM on LoS (Song et al., 2018) (and it continues to dominate in many other domains (Mousavi et al., 2020)). Therefore, the LSTM and the Transformer were chosen as key baselines.
Temporal convolution models have previously been applied to the task of early disease detection using longitudinal lab tests (Oh et al., 2018; Razavian and Sontag, 2015; Razavian et al., 2016), yielding similar results to the LSTM. We highlight two main differences in our work: we introduce a set of pointwise convolutions in parallel, and the temporal convolution filters do not share their parameters between features, allowing the model to optimise processing in spite of heterogeneity in the temporal characteristics. We demonstrate via ablation studies how these design choices contribute substantial improvements to the patient state representation, yielding state-of-the-art results on LoS prediction.
Methods
We want our model to extract both temporal trends and inter-feature relationships in order to capture the patient’s clinical state. Consider a patient who is experiencing slowly worsening respiratory symptoms but is otherwise stable. As this patient is unlikely to be weaned from their ventilator in the near future, a clinician might anticipate a long remaining LoS, but how do they come to this conclusion? Intuitively, one of the factors they are evaluating is the trajectory of the patient e.g. they may ask themselves “Is the respiratory rate getting better or deteriorating?”. However, they can obtain a better indication of lung function by combining certain features e.g. the PaO2/FiO2 ratio, and then looking at how these vary over time. A model should therefore be adept at extracting and combining both intra-feature temporal statistics and inter-feature relationships.
2. Temporal Convolution
Temporal Convolution Networks (TCNs) (van den Oord et al., 2016; Kalchbrenner et al., 2016) are a subclass of convolutional neural networks (Fukushima, 1980) that convolve over the time dimension. They operate on two key principles: the output is the same length as the input, and there can be no leakage of data from the future. We use stacked TCNs to extract temporal trends in our data. Unlike most implementations including (Razavian et al., 2016), we do not share weights across features i.e. weight sharing is only across time (like in Xception (Chollet, 2017)). This is because our features differ sufficiently in their temporal characteristics to warrant specialised processing.
We define the temporal convolution operation for the feature in the layer as
We concatenate the temporal convolution outputs for each feature, as follows
We use to denote concatenation i.e. . In our case, the output dimensions are , where is the number of temporal input features. Throughout this section we label terms with numbers , etc. corresponding to objects in Figure 3. We recommend following this alongside the equations.
3. Pointwise Convolution
4. Skip Connections
We propagate skip connections (He et al., 2015) to allow each layer to see the original data and the pointwise outputs from previous layers. This helps the network to cope with sparsely sampled data. For example, suppose a particular blood test is taken once per day. In order not to lose temporal resolution, we forward-fill these data (Section 4) and convolve with increasingly dilated temporal filters until we find the appropriate width to capture a useful trend. However, if the smaller filters in previous layers (which did not see any useful trend) have polluted the original data by re-weighting, learning will be harder. Therefore, skip connections provide a consistent anchor to the input. They are concatenated (like in DenseNet (Huang et al., 2017), and are arranged in the shared-source connection formation (Wang et al., 2018)) as illustrated in Figure 2. The skip connections expand the feature dimension, , to accommodate the pointwise outputs, and also the channel dimension to fit the original data, . This is best visualised in Figure 3.
5. Temporal Pointwise Convolution
Our model – which we refer to as Temporal Pointwise Convolution (TPC) – combines temporal and pointwise convolution in parallel. Firstly, the temporal output is combined with the skip connections to form (Step 3 in Figure 3).
is then concatenated with the pointwise output after it has been broadcast times. We can therefore define the TPC layer as
6. Loss Function
The remaining LoS has a positive skew (shown in Figure 1) which makes the prediction task more challenging. We address this by replacing the commonly-used mean squared error (MSE) loss with mean squared log error (MSLE).
MSLE penalises proportional errors, which is more reasonable when considering an error of e.g. 5 days in the context of a 2-day stay vs. a 30-day stay. The difference can be seen in Figure 4. For bed management purposes it is particularly important not to harshly penalise over-predictions – the model will become overly cautious and regress its predictions towards the mean. This is counter-productive because long stay patients have a disproportionate effect on bed occupancy.
Data
We use the eICU Collaborative Research Database (Pollard et al., 2018), a multi-centre dataset collated from 208 care centres in the United States, available through PhysioNet (Goldberger et al., 2000). It comprises 200,859 patient unit encounters for 139,367 unique patients admitted to ICUs between 2014 and 2015.
We selected all adult patients (>18 years) with an ICU LoS of at least 5 hours and at least one recorded observation, resulting in 118,535 unique patients and 146,671 ICU stays. We selected 87 time series from the following tables: lab, nursecharting, respiratorycharting, vitalperiodic and vitalaperiodic. To be included, variables had to be present in at least 12.5% of patient stays, or 25% for lab variables. As shown in Figure 5, the lab variables tend to be sparsely sampled. To help the model cope with this missing data, we forward-filled over the gaps. This is more realistic than interpolation as the clinician would only have the most recent value. We then added ‘decay indicators’ to specify where the data is stale. The decay was calculated as , where is the time since the last recording. This is similar in spirit to the masking used by Che et al. (2018).
We extracted diagnoses from the pasthistory, admissiondx and diagnoses tables, and 17 static features from the patient, apachepatientresult and hospital tables (see Tables 5 and 16, and Appendix B for the full list of features and further details).
2. MIMIC-IV Database
We verify our results on a second dataset, the Medical Information Mart for Intensive Care (MIMIC-IV v0.4) database (Johnson et al., 2020), a de-identified and publicly available EHR dataset from the Beth Israel Deaconess Medical Center containing 69,619 ICU stays from 50,048 patients admitted between 2008 and 2019.
We use the same cohort selection criteria as in eICU to select 69,609 ICU stays from 50,042 patients. We followed the same feature selection process to obtain a short list of 172 time series from the chartevents and labevents. We manually removed 71 of these from chartevents because the variable did not vary over time, or because the distribution was not found to provide useful discrimination between patients (see Table 17 for the final list of features). We filled the missing data in the same way as in eICU. We extracted 12 flat features from the icustays, admissions, patients and chartevents tables (Table 6). We did not extract diagnoses from MIMIC-IV because they are not associated with reliable timestamps.
Experiments
In this section, we describe the prediction tasks, baseline models and evaluation metrics. As in Harutyunyan et al. (2019) the training and test data was fixed upfront – the patients were divided such that 70% were used for training, 15% for validation, and 15% for testing.
We assign a remaining LoS target to each hour of the stay, beginning at 5 hours and ending when the patient dies or is discharged. We train the models to make a prediction every hour of the stay. We only include the first 14 days of any patient’s stay to protect against very long batches which would slow down training. This cut-off applies to <5% of patient stays, but it does not affect their maximum remaining LoS values.
1.2. In-Hospital Mortality
We also tested the performance of the models on mortality prediction. Unlike LoS, these labels remain static throughout the patient stay. We used the same training procedure as the LoS task i.e. one prediction each hour. However, to reflect the approach taken by Purushotham et al. (2017) and Harutyunyan et al. (2019), we only report the mortality performance once per patient (at 24 hours into the stay). This means that the cohort represented in the mortality metrics in Table 4 is smaller (16,239 of 21,889 test stays in eICU and 8,320 of 10,264 test stays in MIMIC-IV).
1.3. Multitask
Previous work has found merit in a multitask approach to patient outcome prediction (Harutyunyan et al., 2019; Sheikhalishahi et al., 2019). We investigated whether we would see a similar benefit in the TPC model. When combining the LoS and mortality losses, we applied a relative weighting to the mortality loss – dictated by a parameter (which was treated as a hyperparameter). Further information on the hyperparameter search and implementation details is in Appendix C.
2. Baselines
We include the following baselines in our experiments:
These always predict 3.47 and 1.67 days respectively for eICU and 5.70 and 2.70 days for MIMIC-IV (these correspond to the mean and median of the training data). This is to benchmark the level of performance which is achievable ‘for free’ just by predicting in a reasonable range, and to provide points of reference when setting performance expectations for each dataset.
These are generated by a risk assessment scoring model which is evaluated only once per patient at 24 hours. Therefore it cannot be compared directly, but we include it only as a point of reference for a widely used clinical model. APACHE-IV is only present in the eICU dataset.
Our standard LSTM is similar to Harutyunyan et al. (2019).
Again similar to Harutyunyan et al. (2019), this consists of a set of independent LSTMs that process each feature separately before concatenation (note the similarity with the independent temporal convolutions in the TPC model).
This model takes advantage of multi-head self-attention. Like the TPC model, it is not constrained to progress one timestep at a time; however, unlike TPC, it is not able to scale its receptive fields or process features independently.
3. Evaluation Metrics
We report on 6 LoS metrics: mean absolute deviation (MAD), mean absolute percentage error (MAPE), mean squared error (MSE), mean squared log error (MSLE), coefficient of determination () and Cohen Kappa Score. This is important because bad models can ‘cheat’ particular metrics just by being close to the mean or median value (see Appendix D for additional discussion on this).
3.2. In-Hospital Mortality
In the mortality and multitask experiments we report the area under the receiver operating characteristic curve (AUROC) and the area under the precision recall curve (AUPRC).
Results
In this section, we analyse the model in several ways. Firstly, we report overall performance and compare against a set of baselines. Next, we examine the role of the loss function. Finally, we perform a set of ablation studies to find out which components of the model architecture contribute the most to its success.
The TPC model outperforms all of the baseline models on every metric on both datasets (Table 2) – particularly those that are more robust to skewness: MAPE, MSLE and Kappa. Discounting APACHE-IV, the best performing baseline across both datasets is the Transformer (although the channel-wise LSTM (CW LSTM) is similar on eICU). This is consistent with Harutyunyan et al. (2019) (for CW LSTM) and Song et al. (2018) (for Transformers), who found small improvements over standard LSTMs.
Although the pattern of results is remarkably similar between eICU and MIMIC-IV, there are notable differences in the magnitudes of the metrics. These differences can be attributed to their LoS distributions – the positive skew is more severe in MIMIC-IV (Table 1). This skew has a disproportionate impact on the absolute error, which is captured in the MSE and MAD metrics. Interestingly, the Kappa score is higher in MIMIC-IV because the model can assign the longest stay patients to the >8 day bin, whereas eICU has more medium stay patients in the 3-8 day range which need to be precisely placed. The most comparable results are the MSLE and MAPE metrics, both of which penalise the proportional error, making them more robust to shifts in the LoS distribution.
2. Ablation Studies
To understand the impact of each design choice for the TPC model, we study performance under different ablations on the eICU dataset. The results of these ablations are reported in Table 3.
The first two rows of Table 3 show that using the MSLE (rather than MSE) loss function leads to significant improvements in the TPC model, with large performance gains in MAD, MAPE, MSLE and Kappa, while conceding little in terms of MSE and . The MSE results for the other models are in Appendix Table 13; they show a similar pattern to the TPC model.
2.2. Model Architecture
The second subtable shows that the temporal-only model is superior to the pointwise-only model, but neither reaches the performance of the TPC model. The temporal-only model performs much better than its weight-sharing variant, which demonstrates the importance of having independent parameters per feature. Note that the temporal-only model with weight sharing is the most similar to the approach taken by Razavian et al. (2016), and the results are comparable to the LSTM which is consistent with the results presented in the paper. Removing the skip connections reduces performance by 5-25%. Together the ablation studies demonstrate that the superior performance of the TPC model is the culmination of multiple design decisions.
2.3. Data
We also tested the models without the diagnoses or decay indicators. Perhaps surprisingly, we found that the exclusion of diagnoses does not seem to harm the model. This could be because the relevant diagnoses for predicting LoS e.g. Acute Respiratory Distress Syndrome (ARDS), are discernible from the time series alone e.g. PaO2, FiO2, PEEP etc. The decay indicators contribute a small (but statistically significant) benefit. Their contribution is more obvious in the pointwise-only model where all of the metrics see improvements of 5-23%. This difference is expected since they might reveal some of the temporal structure to the pointwise model e.g. reveal links between up-to-date observations and patient deterioration.
In Appendix E we tested the models without the laboratory tests (which are infrequently sampled) and without the other time series (which tend to be regularly monitored). They indicate that the TPC model is able to exploit disparate EHR time series more successfully than the baselines. They also show that the advantage of the CW LSTM over the standard LSTM is only apparent when the model has to process different types of time series simultaneously.
3. Mortality and Multitask Performance
We investigated adding in-patient mortality as a side-task to improve LoS prediction. Table 4 shows the TPC performance both on single-task mortality prediction, as well as the multi-task setting. We observe first that TPC achieves good performance on mortality alone. Comparing the impact on LoS forecasting in the multi-task setting, we see significant improvements on every metric. Multi-task performance for all baselines is reported in Tables 14 and 15 in the Appendix, where the multitask training confers a more modest benefit.
Further Analyses
In this section, we further explore the performance and behaviour of the TPC model for LoS prediction on the eICU dataset. We test its capacity to exploit smaller datasets, explore which features it uses, and provide a visualisation of the reliability of the model. Finally, we simulate the potential use of the model for bed planning.
The TPC model consistently outperforms the baselines when the training data is small, but we noticed even greater potential for big data. We tested the TPC, LSTM, CW LSTM, and Transformer models with 6.25%, 12.5%, 25%, 50%, and 100% of the eICU training data. TPC maintains the best test performance on all data sizes, with an increasing benefit for larger data. Figure 6 shows the effect on MSLE (the full results for all metrics are included in Table 12).
2. Feature Importance
We used the integrated gradients method (Sundararajan et al., 2017) to calculate feature attributions for the LoS estimates in the eICU dataset. This method computes the importance scores by accumulating gradients interpolated between a baseline input b (intended to represent the absence of data) and the current input x:
where the TPC model is represented as . We use the mean feature values as our baseline input vector. We take the absolute attribution values when a single LoS prediction is made for each patient at 24 hours. We aggregate by taking the mean along the time dimension and then the patient dimension to obtain Figure 7. The background and intuition behind the method is explained clearly by Sturmfels et al. (2020).
Analysing Figure 7, we note that the top features are all strong indicators of organ failure: troponin I is a specific biomarker of myocardial infarction; peak inspiratory pressure, O2 L/%, TV/kg IBW, plateau pressure, PEEP and tidal volume indicate mechanical ventilation (on account of respiratory failure); PTT, ALT (SGPT), AST (SGOT) and alkaline phosphatase suggest liver disease; and high BUN and bilirubin levels point towards kidney failure. Additionally we see infection markers such as lactate, basophils and eosinophils which could indicate sepsis. Both multi-organ failure and sepsis are known causes of extended LoS in the ICU (Böhmer et al., 2014).
3. Evaluation by Use-Case
We have reported aggregate performance metrics indicating strong performance of the TPC model for overall LoS forecasting. In this section, we provide further evaluations tailored to two potential users – an individual ICU clinician, and a bed manager for the unit.
Although aggregate measures of performance are typically reported, these can mask underlying variability in model performance. Such variability can undermine trust or result in unsafe application of systems (Sendak et al., 2020b). In this section, we think of a clinician who wishes to interpret the prediction of the system for an individual patient. We break down the aggregate performance metrics based on factors which will be readily-available at the time of the prediction. Specifically, we visualise the MAPE (chosen for its interpretability) as a function of the time since admission and the predicted remaining LoS.
Figure 8 shows an example for the TPC model on eICU. We can see that high predicted remaining LoS on the first day of a patient’s stay can be quite unreliable, with performance rapidly improving over time. Additional investigation revealed these initial predictions to be under-predictions, indicating that it is challenging to accurately forecast very long LoS for patients on their first day. The long tail of LoS in the dataset reflects the abundance of short-stay patients. The model therefore seems to wait for 1-2 days of data to justify a long LoS prediction. The system can therefore be equipped with instructions indicating that a high predicted remaining LoS on the first day should not be acted upon. This could complement information provided on a model card (Mitchell et al., 2019; Sendak et al., 2020b).
3.2. ICU-level Bed Management
From the perspective of a bed manager, aggregate performance of the model is important: an over-prediction for one patient could be offset by an under-prediction for another, resulting in the same net bed availability. To investigate this, we performed a simulation study. We ran 500 ICU simulations by randomly selecting 16 examples from the eICU test set to form a ‘virtual cohort’. The number 16 was chosen because US hospitals have, on average, 24 ICU beds (Wallace et al., 2015) with an occupancy rate of 68% (Halpern and Pastores, 2015). Figure 9 shows the number of patients remaining in the ICU (of the selected cohort; we do not visualise incoming ICU admissions) using their true remaining LoS (blue). We compute the error (red) between the predictions (green) and true values. The model is well calibrated when predicting patients who are going to stay for at least 1 day. After this, the model tends to under-predict the occupancy by approximately 0.8 patients, corresponding to a small bias towards under-estimating the remaining LoS.
Discussion
We have shown that the TPC model outperforms all baseline models in all task settings (LoS, mortality or multitask) on both the eICU and MIMIC-IV datasets. To explain the success of TPC, we start by examining the parallel architectures in the TPC model. Each component has been designed to extract different information: trends from the temporal convolutions and inter-feature relationships from the pointwise convolutions. The eICU ablation studies reveal that the temporal element is more important, but we stress that their contributions are complementary since the best performance is achieved when they are used together.
Next, we highlight that the temporal-only model far outperforms its most direct comparison, the CW LSTM, on all metrics. Theoretically, they are well matched because they both have feature-specific parameters but are restricted from learning cross-feature interactions. To begin to explain this, we consider how the information flows through the model. The temporal-only model can directly step across large time gaps, whereas the CW LSTM is forced to progress one timestep at a time. This gives the CW LSTM the harder task of remembering information across a noisy EHR with distracting signals of varying frequency. In addition, the temporal-only model can tune its receptive fields for improved processing of each feature thanks to the skip connections (which are not present in the CW LSTM).
The difference in performance between the temporal-only model with and without weight sharing provides strong evidence that assigning independent parameters to each feature is important. Some EHR time series are irregularly and sparsely sampled, and can exhibit considerable variability in the temporal frequencies within the underlying data (evident in Figure 5). This presents a challenge for any model, especially if it is constrained to learn one set of parameters to suit all features. The relative success of the CW LSTM over the standard LSTM when processing disparate time series – but not similar – also lends weight to this theory.
However, the assignment of independent parameters to each feature does not explain all the successes of TPC e.g. the TPC model can process disparate time series and gain more marginal performance than the CW LSTM (Table 11). We need to consider that periodicity is a key property of EHR data – this is true in both the sampling patterns and in the underlying biology e.g. medication schedules, sleep cycles, meals etc. The temporal component of the TPC model is the only architecture with an inherent periodic structure (from the stacked temporal filters) which makes it much easier to learn EHR trends. By comparison, a single attention head in the Transformer model does not look at timepoints a fixed distance apart, but can take an arbitrary form. This is a strength for natural language processing, given the variety of sentence structures possible, but it does not help the Transformer to process EHRs.
Additionally, we have shown that the TPC model outperforms baselines on in-hospital mortality both as a standalone task and in combination with LoS. The performance on both mortality and LoS is significantly better in the multitask setting (this is consistent with past works (Harutyunyan et al., 2019; Sheikhalishahi et al., 2019)) because multitask learning helps to regularise the model and reduce the chance of overfitting (Ruder, 2017). Adding further tasks may be a valid strategy to improve LoS performance.
Finally, we reiterate that using MSLE loss instead of MSE greatly mitigates for positive skew in the LoS task, and this benefit is not model-specific (all of the baselines perform better with MSLE – see Table 13). This demonstrates that careful consideration of the task – as well as the data and model – is an important step towards building useful tools in healthcare.
Our work has several limitations. We know that LoS is heavily influenced by operational factors, and clinical practices can change over time (Kalra et al., 2010). Capacity to maintain performance over time is an important consideration before a system could be used in practice. In future work, it would be instructive to test how quickly the models become out-of-date by reserving more recent data as a test set (Nestor et al., 2018). Although we have included a large set of baselines, we acknowledge that a more exhaustive comparison could be performed, for example comparing by Gaussian Processes (Prasad et al., 2017) or ODE-RNNs (Rubanova et al., 2019; De Brouwer et al., 2019) for handling irregularly sampled time-series. Finally, although we have motivated our study by bed management, this work describes a methodological proof of concept and does not constitute a real clinical system. Prospective study and integration into a real-world EHR is necessary to demonstrate real-world benefit, both of which pose their own challenges (Rajkomar et al., 2018; Sendak et al., 2020a).
In future work, we would like to investigate why the TPC model gains more from the multitask setting than the other models. It seems likely that it is related to additional regularisation provided by the mortality task, but further investigation is needed to confirm our speculations.
Conclusion
We have proposed and evaluated a new deep learning architecture, which we call ‘Temporal Pointwise Convolution’ (TPC). TPC combines temporal convolutional layers with pointwise convolutions to extract temporal and inter-feature information. We have shown that the TPC model is well-equipped to analyse EHR time series containing missingness, differing frequencies and sparse sampling. We believe that the following four aspects contribute the most to its success:
The combination of two complementary architectures that are able to extract different features, both of which are important.
The ability to step over large time gaps.
The capacity to specialise processing to each feature (including the freedom to select the receptive field size for each).
The rigid spacing of the temporal filters, making it easy to derive trends.
From a clinical perspective, we have contributed to the advancement of LoS prediction models, a prerequisite for automated bed management tools. Improving the practice of bed management promises cost reduction (Halpern and Pastores, 2015) and better resource allocation (Mathews and Long, 2015) worldwide. From a computational perspective, we have provided key insights for retrospective EHR studies, particularly where LSTMs are the currently model of choice. In the broader context of machine learning for healthcare we have demonstrated that careful consideration of the complexities of health data is necessary to gain state-of-the-art performance in these tasks.
Acknowledgements
The authors would like to thank Alex Campbell, Petar Veličković, and Ari Ercole for helpful discussions and advice. We would also like to thank Louis-Pascal Xhonneux, Seyon Sivarajah, Rudolf Cardinal, Jacob Deasy, Paul Scherer, and Katharina Kohler for their help in reviewing the manuscript. Finally we thank the Armstrong Fund, the Frank Edward Elmore Fund, and the School of Clinical Medicine at the University of Cambridge for their generous funding.
References
Appendix A Model Architecture: Further Details
After TPC layers, we apply two further pointwise convolutions to obtain the final predictions. Formally, these final steps (shown in Figure 10) can be written as
Appendix B Feature Pre-processing
We selected 17 static features from eICU (shown in Table 5) and 12 from MIMIC-IV (Table 6). Discrete and continuous variables were scaled to the interval , using the 5th and 95th percentiles as the boundaries, and absolute cut offs were placed at . This was to protect against large or erroneous inputs, while avoiding assumptions about the variable distributions. Binary variables were coded as 1 and 0. Categorical variables were converted to one-hot encodings.
B.2. Diagnoses
Here we only describe pre-processing for eICU since MIMIC-IV did not contain coded diagnoses with appropriate timestamps.
Like many EHRs, diagnosis coding in eICU is hierarchical. At the lowest level they can be quite specific e.g. “neurologic disorders of vasculature stroke hemorrhagic stroke subarachnoid hemorrhage with vasospasm”. To maintain the hierarchical structure within a flat vector, we assigned separate features to each hierarchical level and use binary encoding. This produces a vector of size 4,436 with an average sparsity of 99.5% (only 0.5% of the data is positive). We apply a 1% prevalence cut-off on all these features to reduce the size of the vector to 293 and the average sparsity to 93.3%. If a disease does not make the cut-off for inclusion, it is still included via any parent classes that do make the cut-off (in the above example we record everything up to “subarachnoid hemorrhage”). We only included diagnoses that were recorded before the 5th hour in the ICU, to avoid leakage from the future.
Many diagnostic and interventional coding systems are hierarchical in nature: ICD-10 classification (World Health Organisation, 2011), Clinical Classifications Software (Elixhauser et al., 2015), SNOMED CT (De Silva et al., 2011) and OPCS Classification of Interventions and Procedures (NHS Digital, 2019), so this technique is generalisable to other coding systems present in EHRs.
B.3. Time Series
For each admission, we extracted 87 time-varying features from eICU (Table 16) and 101 from MIMIC-IV (Table 17) for each hour of the ICU visit, and up to 24 hours before the ICU visit. The variables were processed in the same manner as the static features. In general, the sampling is very irregular, so the data was re-sampled according to one hour intervals and forward-filled. After forward-filling is complete, any data recorded before the ICU admission is removed. Decay indicators are added as described in Section 4.
Appendix C Hyperparameter Search Methodology and Implementation Details
The TPC model and baselines have hyperparameters that can broadly be split into three categories: time series specific, non-time series specific and global parameters (shown in more detail in Tables 7, 8 and 9). The hyperparameter search ranges have been included in Table 10.
First, we ran 25 randomly sampled hyperparameter trials on the TPC model to decide the non-time series specific parameters (diagnosis embedding size, final fully connected layer size, batch normalisation strategy, dropout rate and the parameter ) keeping all other parameters fixed. These parameters (indicated by stars) remained fixed for all the models which share their non-time series specific architecture (NB. the best value for was 100 – not shown in the Tables).
We then ran 50 hyperparameter trials to optimise the remaining parameters for the TPC, standard LSTM, and Transformer models. To train the channel-wise LSTM and the temporal model with weight sharing, we ran a further 10 trials to re-optimise the hidden size (8 per feature) and number of temporal channels (32 channels shared across all features) respectively. For all other ablation studies and variations of each model, we kept the same hyperparameters where applicable (see Table 2 for a full list of all of the models). The number of epochs was determined by selecting the best validation performance from a model trained over 50 epochs. This was different for each model. For eICU this was 8 (LSTM), 30 (CW LSTM), 15 (Transformer) and 15 (TPC). For MIMIC-IV this was 8 (LSTM), 20 (CW LSTM), 15 (Transformer) and 10 (TPC). We noted that the best LSTM hyperparameters (Table 8) were similar to that found in Sheikhalishahi et al. (2019).
All deep learning methods were implemented in PyTorch (Paszke et al., 2019) and were optimised using Adam (Kingma and Ba, 2014). The data (including decay indicators) and the non-time series components of the models were the same as in TPC (Figure 10). We used trixi to structure our experiments and compare different hyperparameter choices (Zimmerer et al., 2017).
The experiments were performed using resources provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service (www.hpc.cam.ac.uk) funded by EPSRC Tier-2 capital grant EP/P020259/1.
The Transformer is a multi-head self-attention model, originally designed for sequence-to-sequence tasks in natural language processing. It consists of both an encoder and decoder, however we only use the former. Our implementation is the same as the original encoder in Vaswani et al. (2017), except that we add temporal masking to impose causality i.e. the current representation can only depend on current or earlier timepoints, and we omit the positional encodings because they were not found to be helpful. This is probably because we already have a feature to indicate the position in the time series (Section B.3).
Appendix D Evaluation Metrics
The metrics we use are: mean absolute deviation (MAD), mean absolute percentage error (MAPE), mean squared error (MSE), mean squared loss error (MSLE), coefficient of determination () and Cohen Kappa Score. We modify the MAPE metric slightly so that very small true LoS values do not produce unbounded MAPE values. We place a 4 hour lower bound on the divisor i.e.
MAD and MAPE are improved by centering predictions on the median. Likewise, MSE and are bettered by centering predictions around the mean. They are more affected by the skew. MSLE is a good metric for this task, indeed, it is the loss function in most experiments, but is less readily-interpretable than some of the other measures. Cohen’s linear weighted Kappa Score (Cohen, 1960) is intended for ordered classification tasks rather than regression, but it can effectively mitigate for skew if the bins are chosen well. It has previously provided useful insights in Harutyunyan et al. (2019), so we use the same LoS bins: 0-1, 1-2, 2-3, 3-4, 4-5, 5-6, 6-7, 7-8, 8-14, and 14+ days. As a classification measure, it will treat everything falling within the same classification bin as equal, so it is fundamentally a coarser measure than the other metrics.
To illustrate the importance of using multiple metrics, consider that the mean and median models are in some sense equally poor (neither has learned anything meaningful for our purposes). Nevertheless, the median model is able to better exploit the MAD, MAPE and MSLE metrics, and the mean model fares better with MSE, but the Kappa score betrays them both. A good model will perform well across all of the metrics.
Appendix E Time Series Ablation
We performed ablations on the type of time series variable that we include: laboratory tests only (labs), which are infrequently sampled, and all other variables (other) which include vital signs, nurse observations, and automatically recorded variables (e.g. from ventilator machines). This shows how well each model can cope with time series exhibiting different periodicity and sampling frequencies. The results are shown in Table 11. The TPC model has the largest percentage gain when the labs and other variables are combined (this is synonymous with the greatest percentage impairment in the ablations). Next are the CW LSTM and Transformer, followed by the LSTM. This suggests that the TPC model is best able to exploit EHR time series with different temporal properties.
When examining the results for LSTM and CW LSTM in more detail, we can see that the CW LSTM only has an advantage when the model has to combine the data types. This supports the hypothesis that the CW LSTM is better able to cope when there are varying frequencies in the data, as it can tailor the processing to each. When the inter-feature variability is small (the same type of time series) they perform similarly.
It is unsurprising that the Transformer does better than the LSTM when combining data types, as it can directly skip over large gaps in time to extract a trend in lab values, while simultaneously attending to recent timepoints for the processing of other variables.
The TPC is the most successful model; its inherent periodic structure helps it to extract useful information from all of the variables. The CW LSTM and Transformer do not have this in their architectures, making the derivation more obscure. The importance of periodicity is discussed in more detail in Section 8.