The AAAI 2023 Workshop on Representation Learning for Responsible Human-Centric AI (R2 HCAI)

Dynamic Outcomes-Based Clustering of Disease Trajectory in MechanicallyVentilated Patients
Emma Rocheteau*1 , Ioana Bica2 , Pietro Liò1 , Ari Ercole3
1

Department of Computer Science and Technology, University of Cambridge, UK
2
Department of Engineering Science, University of Oxford, UK
3
Department of Anaesthesia, University of Cambridge, UK*ecr38@cam.ac.uk
Abstract

may be both more homogeneous and represent a principledsubstrate for individualised intervention. Therefore, there isgreat potential benefit to be gained from phenotype discovery
in order to guide future clinical studies.
Relatively simple clustering techniques have revealed actionable sub-phenotypes by secondary analysis of RCT data
where intervention data is also available. For example, latent trajectory modelling of inflammatory biomarkers hasrevealed pro- and anti-inflammatory sub-types of ARDS (Famous et al. 2017) and these have a differential response toliberal or conservative fluid therapy policies. Clustering oftranscriptomic data has revealed patient populations in whichsteroid therapy may be beneficial in sepsis (Antcliffe et al.
2019). Routinely collected data has also been used to findsimple trajectory clusters for sepsis based on physiologicalparameters (Bhavani et al. 2022) which have been shown tohave differential response to fluids.
Investigating the dynamic clustering of trajectories is appealing in the ICU for reasons touched on above. We knowthat temporal neural network architectures can handle the heterogeneous population in the ICU, yielding high predictive
performance (Rocheteau, Liò, and Hyland 2021). We also
have reason to believe that there is merit in an unsupervised
approach to patient representation (Miotto et al. 2016; Wanget al. 2020). This motivates examining a combined (supervised and unsupervised) approach to patient representation
as a substrate for trajectory clustering.
We choose apply our method to mechanical ventilation, acommon and high-risk intervention, although many clinical
trials (Girbes and de Grooth 2020) have been unrewarding.
Patient heterogeneity is a major stumbling block, becausesome patients can be liberated from mechanical ventilation
quickly, whereas others may need a prolonged period of ventilation and potentially tracheostomy to allow a more gradualtransition to spontaneous breathing. Difficult weaning may
be due to either acute or chronic pulmonary injury, frailty
or deconditioning due to other critical illness (critical illness
may increase physiological demands on the cardiorespiratory
system at the same time as being catabolic i.e. associatedwith muscle breakdown and weakness and this may include
respiratory muscles).
To this end, we have developed a dynamic clustering approach for mechanically ventilated patients in the ICU. Notethat temporal clustering approaches have successfully ap-

The advancement of Electronic Health Records (EHRs) andmachine learning have enabled a data-driven and personalised
approach to healthcare. One step in this direction is to uncover
patient sub-types with similar disease trajectories in a heterogeneous population. This is especially important in the context
of mechanical ventilation in intensive care, where mortality
is high and there is no consensus on treatment. In this work,
we present a new approach to clustering mechanical ventilation episodes, using a multi-task combination of supervised,self-supervised and unsupervised learning techniques. Our dynamic clustering assignment is explicitly guided to reflect thephenotype, trajectory and outcomes of the patient. Experimentation on a real-world dataset is encouraging, and we hopethat we could someday translate this into actionable insights
in guiding future clinical research.

Introduction and Related Work
Patients on mechanical ventilation are a highly heterogeneousgroup, with widely differing outcomes. Some have relatively
healthy lungs e.g. if they are recovering from surgery on
another organ; whereas others have varying degrees of pulmonary failure. Pulmonary failure can be acute e.g. AcuteRespiratory Distress Syndrome (ARDS), or chronic. Whenpatients deteriorate, it can happen rapidly over hours, minutes
or even seconds. Understanding the stability, or otherwise, ofa particular physiological state at a particular time is important in terms of flagging the potentially deteriorating patientand enabling early intervention. Unfortunately, patients onventilators have high mortality rates (Máca et al. 2017; Pooleet al. 2017) and there is no established consensus on optimal treatment strategies from randomised controlled trials(RCTs) (Bein et al. 2016). One important reason for this is theextreme heterogeneity in the ICU population. Unlike diseases
such as ischaemic heart disease or COVID-19 pneumonitis
with potentially devastating but relatively simple aetiologies,
the typical intensive care patient may have a wide variety ofreasons for admission and mechanisms leading to physiological disruption, have complex multi-system disease, complex
and varied comorbidities and treatment histories. There isa great deal of interest in finding stereotypical disease behaviours or ‘phenotypes’ or ‘sub-phenotypes’ in the hopethat these might represent pathobiological ‘endotypes’ which

1plied to other domains e.g. in Parkinson’s (Zhang et al. 2019),
diabetes (Rusanov, Prado, and Weng 2016) and cystic fibrosis (Lee and van der Schaar 2020) and increasingly in intensive care as discussed above. The clusters are designed toshare similarities in phenotype, trajectory and outcomes. We
generate a cluster for each hour of a patient’s stay, meaningthat if an event happens which alters the predicted trajectoryand outcomes, there will be a shift in the cluster assignment.
This is interesting, not only because it can reveal which events
are associated with these shifts, but it may give us insight
into what might have happened if the ventilation strategy
had been different. Experimentation on a real-world dataset
is encouraging, and we hope that this could translate into
actionable insights when guiding future clinical research.

Figure 1: Overview of our model. Only one timestep, t, isshown for simplicity. F and S are the number of time series and static variables respectively. At timestep t, the static
variables (yellow) and preceding time series variables (grey)
and their corresponding decay indicator variables (orange,
explained under ‘Time Series’ in Data Preprocessing in the
Appendix) are given to the encoder, which produces an embedding (green) for timestep t. This is then given to thedecoder networks (yellow), forecasting network (purple) andthe predictor network to obtain the four patient outcomes(red). After training is complete, the test embeddings areused for clustering.

Methods
Broadly, our strategy was to train a temporal encoder toembed the patient data at every timestep (this is analogous toreturning all the hidden states for an LSTM model). We useda mixture of supervised, unsupervised and self-supervised
learning to do this (see ‘Prediction Tasks’ below). Once theencoder training was complete, we used an unsupervised
method to cluster the embeddings, so that we get a cluster
for every timestep in the patient’s ventilation episode. The
code can be found at: https://github.com/EmmaRocheteau/
Mechanical-Ventilation-Clustering.
The data consisted of both timeseries and static features.
The supervised tasks included two binary tasks: predicting
hospital mortality and the risk of receiving a tracheostomy1 ,and two duration tasks: the remaining length of stay (LoS)
from timestep t, and their remaining ventilation duration(VD). This ensured that the patient outcomes are stored within
the embedding. In addition, we trained a decoder to reconstruct timestep t and the static data. This is an unsupervised
approach which encourages the embedding to retain the patient phenotype. Finally, we predicted timestep t + 1, a selfsupervised approach designed to embed the patient trajectory.
See Figure 1 for a schematic.

K-Medoids Clustering We used k-medoids clustering tocluster the learned embeddings. K-medoids is similar to kmeans, except that it operates with medoids rather than centroids. This means that the medoids will always be a true
observation in the data, while that is not usually the casefor centroids. The main advantage is that k-medoids are lesssensitive to outliers than k-means, which is more suitable inthis context where the data is noisy and heavily skewed2 .
Both k-means and k-medoids operate on pairwise similarities. We decided to use Euclidean distance rather thancosine similarity. This is because intuitively, it is not only
the direction that the patient is moving in that matters, butalso the distance along that axis. For example, if a particular‘direction’ represents acute decompensated heart failure, we
also care how severe the decompensation is.
We applied batch normalisation (Ioffe and Szegedy 2015)
to the embeddings, to ensure that the embedding distribution
remained within a reasonable range. The value of k (5 for allmodels) was chosen using the elbow method (see ‘Numberof Clusters’ in the Appendix).

Encoder In recent years, LSTMs have been by far the mostpopular model for predicting clinical outcomes and haveachieved state-of-the-art results (Harutyunyan et al. 2019;
Sheikhalishahi, Balaraman, and Osmani 2019; Rajkomar et al.
2018; Tong et al. 2022). They have also been applied to other
patient prediction tasks e.g. forecasting diagnoses and medications (Choi et al. 2015; Lipton et al. 2015), and mortality
prediction (Che et al. 2018; Harutyunyan et al. 2019; Shickel
et al. 2019). More recently, the Transformer model (Vaswani
et al. 2017) has marginally outperformed the LSTM when predicting LoS (Song et al. 2018). Rocheteau, Liò, and Hyland(2021) showed that Temporal Pointwise Convolution (TPC)
outperformed both the LSTM and Transformer models on
mortality and LoS. Therefore, we chose to investigate thesethree encoders. Details of their implementation are given in
the Appendix.

Data
We used the Amsterdam UMC database version 1.0.2 (Thoral
et al. 2021), which contains 23,106 ICU admissions from
20,109 patients admitted between 2003 and 2016. We selected all of the mechanical ventilation episodes with a minimum duration of 4 hours, capping the maximum duration after 21 days to reduce computational costs. This corresponded
to 14,836 episodes which occurred during 13,502 ICU admissions from a cohort of 12,597 unique patients. Of the 14,836
episodes, 13,783 ended in extubation or death of the patient,2Preliminary experiments revealed that k-means were more
likely to produce small clusters which lay far away from the restof the data, because it is more affected by outliers. This made theclustering process less reliable and reproducible.

1
A tracheostomy is a procedure designed for long term mechanical ventilation of a patient.

2648 ended with a tracheotomy procedure occurring within 21days, 399 patients were still on ventilation at 21 days, and 6
patients were converted to a non-invasive ventilation setting.
We selected 31 time series features and 14 static features.
The data were split such that 70%, 15% and 15% were used
for training, validation and testing respectively. These weresplit by patient, not ventilation episode, to avoid data leakage
from the train set. The cohort summary is shown in Table 1.
Further details of the data are provided in the Appendix.

tracheostomy respectively). In order to encourage the model
to prioritise learning these important outcomes, we appliedclass weighting to the task (where the weight is proportional
to the inverse of the frequency of each outcome). We usedbinary crossentropy as the loss function for both tasks. We
report the area under the receiver operating characteristic
curve (AUROC) and the area under the precision recall curve(AUPRC) as metrics.
Reconstruction and Forecasting As shown in Figure 1,
we use the embedding to reconstruct the timestep t, and forecast one timestep (t + 1) ahead. For the reconstruction oft and forecast of t + 1, we apply the mean squared error
since the data can (very approximately) be assumed to be
Gaussian centred at 0 following normalisation. We also reconstruct the following static features: sex, urgency of admission,
agegroup, weightgroup, and heightgroup. The first two arebinary, and so we apply the binary crossentropy loss function. The other three are ordered categorical (as explained in‘Static Features’ in the Appendix), therefore we use the meansquared error loss function. Since these tasks are auxiliary(we are not interested in the performance as an outcome of
the model), we reported their loss function values as ‘metrics’since they do not need to be interpretable.
The relative weightings of all of these tasks are given under‘Hyperparameter Search Methodology’ in the Appendix.

Table 1: Cohort summary for the Amsterdam UMC database.
‘Remaining LoS’ refers to the remaining duration in the hospital after the start of the ventilation episode.
Number of ventilation episodes
Train
Validation
Test

14,836
10,3952,230
2,208

Sex (% male)
Total LoS in days (mean)
Total LoS in days (median)Remaining LoS in days (mean)
Remaining LoS in days (median)
VD in days (mean)
VD in days (median)In-hospital mortality
Tracheostomy patients‘Urgent’ patients

66.6%
8.26
2.13
7.26
2.01
3.95
0.83
14.6%
7.4%
28.1%

Number of input features
Time series
Static45
31
14

Results
In this section, we highlight important performance differences between the three encoders, analyse an ablation study
on the tasks, and provide a detailed analysis of the clusters
produced by the TPC model. A deeper evaluation of theresults can be found in the discussion.

Prediction Tasks
Remaining Length of Stay and Ventilation Duration We
assigned a remaining length of stay (LoS) and remaining ventilation duration (VD) target to each hour of the ventilation
episode, ending when the patient dies or is extubated. The remaining LoS is calculated by subtracting the time elapsed in
the ICU from the total LoS. The remaining VD is calculated
by subtracting the time elapsed in the ventilation episodefrom the total VD. We only trained on data from the first
21 days of the ventilation episode to protect against batches
becoming overly long and slowing down training.
The remaining LoS and VD each have a significant positive skew which makes the duration tasks more challenging.
We partly circumvent this by replacing the commonly usedmean squared error (MSE) loss with mean squared log error(MSLE), as in Rocheteau, Liò, and Hyland (2021). We reported on 2 LoS and VD metrics: mean absolute deviation(MAD) and mean squared log error (MSLE). The MAD wasused as the primary metric in Harutyunyan et al. (2019) but
MSLE is arguably the more holistic metric (Rocheteau, Liò,and Hyland 2021).

Task Performance(a) – Full Task Setting The TPC model performs significantly better than the LSTM and Transformer on the outcome
tasks (Table 2a), which is in line with previous findings in
MIMIC-IV and eICU (Rocheteau, Liò, and Hyland 2021).
The superiority of the TPC model is also evident in the variational and ablation experiments. Interestingly, the Transformer performs poorly on the binary tasks but better on the
duration tasks with respect to the LSTM. Additionally, theLSTM performs the best on the reconstruction and forecasting tasks (Table 14a). Possible reasons for these findings are
explored in the discussion.
(b) – Variational Embedding Spaces We experimented
with making the embeddings ‘variational’, by representing
the embedding as a set of means and standard deviations toallow sampling of embedding coordinates. The rationale wasthat by forcing the embedding space to be smoother, we might
improve the quality of the clustering as the distances between
patients in the embedding space become more reliable. However, this was found to universally hurt performance (Table 2b
and Table 14b) and it produced clusters which were morehomogeneous in terms of outcomes and features, which wascounter to the aim of producing clinically distinct clusters.

Mortality and Tracheostomy Unlike the duration tasks,
these tasks are static, i.e. the labels do not change during theventilation episode. Both tasks have significant class imbalance (only 14.6% and 7.4% of patients died or received a3

Table 2: Encoder performance on the prediction tasks averaged over 5 independent training runs. The error margins are 95%
confidence intervals. For mortality and tracheostomy, higher AUROC and AUPRC is better; for LoS and VD, lower MAD andMSLE is better. (a) shows the full multi-task setting as shown in Figure 1, (b) is a variational alternative to the full task setting.
Statistically significant differences are indicated by daggers († = p < 0.05, ‡ = p < 0.001). If the result is significantly better thanthe comparison models*, it is highlighted in blue, if it is significantly worse it is highlighted in pink. *In (a) the statistical testing
compares the three model types, in (b) each model type is compared to its corresponding ‘non-variational’ model in table (a).

Model

In-Hospital Mortality
AUROC
AUPRC

Tracheostomy
AUROC
AUPRC

Length of Stay
MAD
MSLE

Vent. Duration
MAD
MSLE

TPC
(a) Transformer
LSTM

0.833±0.010† 0.644±0.013‡ 0.804±0.007‡ 0.507±0.020† 7.20±0.13‡ 0.359±0.010‡ 3.24±0.07‡ 0.210±0.008‡
0.697±0.012 0.434±0.019 0.760±0.012 0.419±0.033 8.46±0.07 0.495±0.007 3.95±0.20 0.256±0.016
0.823±0.002 0.608±0.008 0.774±0.002 0.473±0.015 9.16±0.06 0.663±0.008 5.57±0.04 0.681±0.011

TPC(b) Transformer
LSTM

0.807±0.006‡ 0.584±0.014‡ 0.775±0.008‡ 0.437±0.012‡ 9.06±0.10‡ 0.555±0.018‡ 4.42±0.03‡ 0.347±0.006‡
0.660±0.023† 0.373±0.039† 0.714±0.020‡ 0.353±0.018† 9.42±0.27‡ 0.623±0.020‡ 4.63±0.27‡ 0.359±0.030‡
0.803±0.004‡ 0.555±0.006‡ 0.748±0.005‡ 0.411±0.010‡ 10.2±0.1‡ 0.813±0.016‡ 5.95±0.04‡ 0.775±0.007‡

Ablation Studylung. Furthermore, they have the least compliant lungs
with higher peak inspiratory pressure and higher PEEP.
This means that the physical properties of the lung havebeen damaged and it is harder to inflate the lung. The
P/F ratio is also very low, which indicates that the lungfunction is also impaired i.e. the patients cannot absorb asmuch oxygen. This is in keeping with severe respiratory
distress. We could describe this phenotype as a ‘early,
life-threatening pulmonary injury’ patient group.
• Cluster 2 also display substantial mortality and, again,from Table 4 it becomes clear that these patients also represent a group with severe pulmonary dysfunction likecluster 1. However this phenotype is characterised byvery long LoS and VD, with consequent high rates of tracheostomy: this phenotype represents patients who are difficult to wean from mechanical ventilation with resultant
long ICU stay. This might be described as a ‘pulmonary
critical illness’ phenotype.
• Cluster 3 have the best outcomes, with short LoS andlow mortality. They are extubated without tracheostomy,
because they have less pulmonary injury. This appears to
be a ‘short stay’ phenotype who require a brief period oforgan support, perhaps after significant surgery.
• Cluster 4 have relatively low mortality but high rates oftracheostomy. Table 4 shows modest levels of respiratory
failure and good lung compliance. Thus, whilst these patients are difficult to wean from mechanical ventilation(like cluster 2), this is due to factors that are not primarily related to pulmonary pathology: once they receive a
tracheostomy they are rapidly liberated from mechanical ventilation. We could therefore describe them as a‘general critical illness’ phenotype.
• Cluster 5 shows a moderate to severe group, who are notas acutely unwell as cluster 1, but are still high-risk. From
Table 4 we see that pulmonary injury is not a prominent
feature so we could characterise these patients as ‘early,
life-threatening non-pulmonary injury’ patients.

We performed an ablation study on the tasks used to train the
representation space. The full set of results and analysis areincluded in the Appendix. However, the trend is such that thebest results for all tasks (except for the duration tasks) areachieved when all tasks are included (Table 6). The reason
for the exception is in the duration only task setting (g), is
explored further in the discussion. Overall, our ablation study
indicates that having multiple competing learning objectives
has a stabilising effect on learning the representation.

Cluster Analysis
As the best performing encoder, we have focused on
analysing the clusters produced by the TPC model. In order to analyse the average differences between the patientsin each cluster, it was necessary to flatten the clustering intoone ‘primary’ cluster per patient. This was to prevent confusion, since patients can enter multiple clusters during their
ICU stay (sometimes only for one or two timepoints), andthis is disproportionately true of the long stay patients. The
cluster in which each patient spent the majority of their timein was assigned its primary cluster. If there were multiple
modes, then the mode experienced later in the sequence waschosen. The next two sections characterise the behaviour ofthe primary clusters. Subsequently, we analyse the dynamic
aspects of the clustering from multiple different perspectives.
Differences in Phenotype and Outcomes Table 3 shows
the mean outcomes for each cluster. We also analysed somekey features in the original data, to visualise differences inpatient phenotype that the model identified. The average values of key features in patients divided by primary cluster are
shown in Table 4. Broadly we can say that:• Cluster 1 contains the sickest patients, with an average
mortality of 72.0%. They are short stay patients and unsurprisingly they have the lowest rate of tracheostomy
as most do not survive or stay long enough to requirecomplex respiratory weaning. Table 4 shows they are primarily ventilated with ‘mandatory’ ventilation settings,
meaning the machine is deciding the respiratory rate, andeither the tidal volume or pressure differential across theOverall, the findings from Tables 3 and 4 show that there
are statistically significant and clinically meaningful differences between the clusters. These can be visualised in Fig4

Table 3: Average outcomes by cluster ± 95% confidence intervals for the TPC model. Each patient has been classified into a
primary cluster, which is the cluster that they spent the majority of their time in. LoS and VD are shown in days.
Cluster

12
3
45

Patients

Mortality (%)

Tracheostomy (%)

Length of Stay

Vent. Duration

232
133
1,292
347
227

72.0±5.8
34.6±8.2
1.9±0.7
4.0±2.1
26.0±5.7

1.3±1.538.3±8.4
1.5±0.731.1±4.9
8.4±3.6

3.8±0.830.0±3.6
2.8±0.3
22.0±1.8
13.0±1.6

2.4±0.3
21.4±2.2
0.7±0.0
7.4±0.9
7.2±0.9

Table 4: Key features averaged by cluster ± 95% confidence intervals. ‘Urgency’ is a flag given to the patient at admission.
Mandatory Ventilation (MV) settings are provided in Table 13. The peak inspiratory pressure, P/F Ratio and PEEP are expressedin mmHg. A normal P/F ratio at sea level is ≈400-500mmHg; whereas 200-300mmHg is consistent with mild ARDS (Forceet al. 2012). Lung compliance is expressed in ml/cmH2 O (normal for a mechanically ventilated patient is 50-100ml/cmH2 O).
Cluster12345

Age 70+ (%)

Sex (% male)Urgency (%)

MV (%)

Peak Insp. Pressure

Lung Comp.

P/F Ratio

PEEP

52.2±6.5
54.1±8.6
39.8±2.725.9±4.640.1±6.4

59.7±6.3
65.8±8.1
69.7±2.568.4±4.969.6±5.963.4±6.3
39.1±8.4
14.9±1.941.5±5.343.2±6.5

68.3±0.843.2±0.438.6±0.6
22.1±0.441.8±0.525.3±0.2
23.2±0.1
16.1±0.1
17.8±0.1
20.3±0.1

32.7±0.536.8±0.358.8±0.757.5±0.447.1±0.4

217±2
220±1
260±1
237±1
243±1

10.09±0.07
9.97±0.03
6.78±0.03
8.19±0.28
8.83±0.38

ure 2, which provides maps of the representation space, showing how the clusters relate to one another as well as the
distribution of outcomes and features.
Medoid Analysis The medoids produced by the clustering
algorithm are shown in Figure 3 and give a description of arepresentative patient in each cluster. Note that each medoid
corresponds to both a patient and a specific time-point in their
ventilation episode3 .
• The medoid patient for cluster 1 (female, age 60-69) died
4 hours after the episode shown without a tracheostomy.
Her deterioration is predictable from multiple parameters in the data but infection (high WBC) and pulmonary
dysfunction are particularly noteworthy.
• The typical medoid patient representing cluster 2 (male,80+ years old) received a tracheostomy 19 days after
the episode shown, and was discharged at 23 days. Thispatient required late as well as early mandatory ventilationsuggesting that he ran into later pulmonary complications– mostly likely infectious since his CRP remains highthroughout the stay.
• The medoid patient in cluster 3 (female, age 60-69) was
discharged from hospital the day after her brief window ofventilation. She does not display substantial physiologicalderangement.
• The patient in cluster 4 (female, age 60-69) received a
tracheostomy 3 days after the sequence shown. Her lungcompliance and P/F ratio are both high compared to clusters 1 and 2, indicating better lung function. Therefore, we
can conclude that she needed a tracheostomy for reasons
other than lung injury. This is also supported her ventilator
settings: low peak inspiratory pressure and PEEP and no

Figure 2: t-SNE plots for the embeddings produced by the
TPC model. For these figures, 1500 random samples were selected from the test set and projected. In each plot, a different
attribute has been highlighted.

3
Not all of the ventilation episode for each patient is shown in
Figure 3 – only from the start until the medoid time.

5

mandatory ventilation modes. Her end tidal CO2 levels
suggest a degree of hypoventilation which improves as
she is slowly weaned from mechanical ventilation.
• Lastly, the patient in cluster 5 (female, 80+ years old)stayed for 9 further days in hospital before being discharged. The short duration of ventilation and relatively
normal pulmonary physiology is again consistent with a
non-pulmonary phenotype.
Cluster 1

Cluster 2

Cluster 3

Cluster 4

Cluster 5

Figure 5: A sankey plot showing the evolution of the clustering across time. We begin at 4 hours to allow the clustering
to stabilise at the start of the time series. At 21 days there
are still some patients without a final outcome (mostly fromcluster 2) but this is because they are ventilated for longer
than 21 days and have been right censored.

Time since ICU admission (hours)

Figure 3: Raw data from each of the medoids. The data havebeen standardised around the mean value for each feature.
Red means the value is high and blue means low. We can seethat each medoid largely follows the average pattern for thecluster shown in Table 4. WBC is white blood count, CRP is
C Reactive Protein, ABP is arterial blood pressure.

contains patients with the longest ventilation episodes, which
can be seen by its low rate of attrition over time.
Number of Clusters per Patient Figure 6 shows that mostpatients remain in only one cluster during their ventilation
episode. However, when the distribution is broken down byprimary cluster, we can see that this is heavily driven by the
behaviour of cluster 3 patients, which tend to remain in cluster 3 for their entire ventilation duration (note that they tend
to have short VDs so this is not so surprising). In contrast,
clusters 2 and 5 most commonly appear alongside other clusters during a single ventilation episode. This means that formost episodes attributed to cluster 2 or 5, there are transitions
either into or out of these clusters. These are explored next.

Temporal Analysis Broadly, there are two perspectives
when evaluating the dynamic aspects of this clustering.
One is the ‘Markovian’ perspective, where we can examine
the transition function between clusters. This is shown in
Figure 4. Unsurprisingly, this reveals that the patient is always
most likely to remain in the same cluster. However the mostcommon inter-cluster transitions are from cluster 5 to cluster
4, and cluster 1 to cluster 5. Note that these clusters are next
to one another and share lengthy borders in Figure 2. Mostof the patients who transition to ‘Died’ come from cluster 1,
and most of the ‘Discharged’ patients come from cluster 3.

Figure 4: A transition matrix for the TPC model, showing
the probability of entering each cluster at time t + 1, plus thecategories ‘discharged’ or ‘died’, given their cluster at time t.

Figure 6: Distribution of the number of clusters that the patient enters during their ventilation episode, separated byprimary cluster (shown by the colour key). For example, cluster 3 (purple) mainly appears on its own i.e. the patient starts
the episode in cluster 3 and remains in cluster 3 for the whole
duration, whereas cluster 5 (red) rarely appears on its own.

The other perspective is to look at the number of patientsin each cluster at different time points after admission, andobserve the transitions between them (Figure 5). Transitions
from cluster 3 to ‘extubated’ are very common within the firstday, but then they almost disappear by 3 days. This cannot beseen with the Markovian perspective in Figure 4. Cluster 2

Cluster Transitions The clusters produced by the TPC
model are remarkably stable over time, given that there is no
explicit loss incentive to constrain the representation to behave in this way. Figure 7 shows the distribution of timepoints
that the patients first enter their primary cluster. Clusters 2
6

Table 5: Stable cluster transitions (origin cluster → destination cluster) with a count of ≥15, sorted by destination cluster. The
median rather than the mean time is displayed to show a more representative time of transition (as there is positive skew).
Transition

Count

Median Time

Mortality (%)

Tracheostomy (%)

Urgency (%)

VD

LoS

3→1
5→1

17
29

3
16

76.5
51.7

0.0
10.3

47.1
55.2

0.5
4.3

0.7
5.3

1→3
5→3

28
46

11
9

10.7
15.2

0.0
4.3

67.9
41.3

1.0
1.2

2.6
6.5

2→4
5→4

28
27

17
10

10.7
11.1

21.4
7.4

42.9
48.1

6.2
3.4

12.8
9.1

1→5
3→5
4→5

25
15
15

3
4
56

44.0
13.3
26.7

4.0
13.3
26.7

68.0
53.3
46.7

3.9
1.9
6.6

6.5
4.6
11.5

in Figure 5, is that the average predicted risk of death drops
from 56.4% 5 hours prior to the transition, to 41.7% at the
point of transition. There is also a corresponding reduction
in tracheostomy risk (-13%), LoS (-17.1% after adjustment4 )and VD (-26.4% after adjustment) as predicted by the model,
and dramatic improvements in physiological parameters suchas lung compliance (+35%) and P/F ratio (+15%).
Another interesting transition is 3→1, which happens
when the model initially believes the patient to be relatively
healthy, but then quickly re-adjusts to predict poor outcomes.
Looking in more detail at the raw data, we discovered thatthese patients are younger (only 23.5% are 70+), which could
explain why the model was initially optimistic and why thedeterioration is so rapid5 . We also observed a deterioration in
the lung compliance (-26.3%) and P/F ratios (-12.8%), and achange in the ventilator settings – namely higher PEEP andpeak inspiratory pressure and lower tidal volumes – reflecting a drop in lung compliance of the patients. Most of thesepatients died within 12 hours of the transition to cluster 1.

Figure 7: Percentage of patients who enter their primarycluster, by time since the start of the ventilation episode.

and 3 are particularly likely to accurately assigned during thefirst hour of ventilation (87% and 89% respectively), whilecluster 4 is the least likely to be identified early (64%).
Next, we investigated what we will refer to as ‘stable’ transitions between clusters. In order to be characterised as stable,
the origin cluster needed to remain stable in the 5 hours preceding the transition, and the patient was not permitted to
re-enter the origin cluster for 5 hours following the transition.
This was primarily to screen out patients who were at theboundary between two clusters, continually crossing backand forth but not representing a true transition from one tothe other. Before screening, there were 22,036 cluster transitions, corresponding to 870 separate ventilation episodes(39% of the total in the test set). Of these transitions, only 291
represented stable movement between clusters. We further
removed any transitions between two clusters that had fewer
than 15 transition examples, as this would be insufficient toanalyse. The remaining 230 transitions are shown in Table 5.
Firstly, it is noteworthy that the outcomes reflect the destination cluster, not the origin cluster. The exception to this isthe ‘urgency’ column, which is not an outcome, but a label
assigned at admission and hence is more likely to reflect theorigin cluster (although it is worth mentioning that the origin
cluster is not necessarily the cluster at admission).
Cluster 5 stands out as being disproportionately involved
in inter-cluster transitions. Of these, the most common is
5→3, which occurs when the model overestimates the risk
to the patient early on in the ventilation episode. Not shown

Reliability We investigated the reproducibility of these phenotypes. We chose to analyse the clusters in the following
settings: i) alternative encoder models, ii) retraining the TPC
model with different random seeds and iii) varying the valueof k. The clusters were found to be surprisingly stable, withkey features of the extracted phenotypes remaining similar
between models. With increasing value of k, we noticed thatrather than completely rearranging the position of the clusters, increasing k progressively subdivides existing clusters,
hinting that the clusters are hierarchically organised (more in
the discussion). The full analysis is included in the Appendix.

Discussion
We evaluated the use of TPC model, trained using supervised,unsupervised and self-supervised learning techniques, for the
purposes of phenotype discovery in mechanically ventilated
patients. We will discuss the most important findings in turn.
4
There is a 5 hour gap between these predictions, therefore thistime difference needs to be removed from the first prediction.
5
This is because younger patients can mask a problem by compensating deceptively well, until they reach a point where the homeostatic mechanisms can no longer cope.

7

Firstly, we reaffirmed that the TPC model performs better
than alternative encoders on EHR data for patient outcome
prediction. This time on the Amsterdam UMC database (Thoral et al. 2021), and with added tasks.
Secondly, we found that the Transformer outperformed
LSTM on LoS and VD, but performed much worse on themortality task, and slightly worse on the tracheostomy task.
This may be because the task weighting was more favourouble to the LSTM and TPC models, whereas the Transformer
would have benefited from greater weighting towards thebinary tasks. Another possibility is that the binary tasks particularly benefit from biases in the LSTM and TPC encoders,
because these models naturally emphasise recent timepoints(and these are especially important for solving the mortalitytask). As for the reason that the Transformer performs relatively better on tracheostomy compared to mortality, it could
be because there is positive correlation between long LoS,
VD and the chance of acquiring a tracheostomy. Solving theduration tasks makes the tracheostomy task easier, whereas
the relation to mortality is more complex (Figure 2).
To briefly comment on the reconstruction results in the
Appendix; it may seem surprising that the LSTM model performs best on the reconstruction and forecasting tasks. However, this could be explained if the LSTM is creating simpler‘lower level’ representations that are easier to translate backto the original data using two-layer decoder networks.
Thirdly, Table 6 reveals a general trend that the more tasks
that are added, the better results across all the tasks, with particular benefits to the tracheostomy task. The exception to thiswas the duration only setting, which achieved stronger results.
There are two possible explanations for the discrepancy:

LoS. This theory would be formally tested by accompanying
the duration tasks with the mortality task only.
Finally, regarding the repeatability of the clustering, we
demonstrated that there are key aspects of the learned representations (both of different encoders and TPC instances)that translate to core phenotypic traits that are consistently
recognised. The separation on other traits, especially when
distinguishing the sickest patients from the moderately ill,
was more malleable. This suggests that perhaps there is not awell defined distinction between these, but rather a scale ofdeterioration, through which an arbitrary line can be drawn.

Limitations and Future Work
Hierarchical Clustering It is evident from the results thatcertain clusters are more related than others. A tree basedhierarchy of clusters seems more natural than a flat structure.
It would be interesting to visualise the shifts in cluster assignment within a tree rather than as discrete jumps. We areparticularly interested in modifying an approach being used
for genetics data (Patel et al. 2022; Chami et al. 2020; Corso
et al. 2021), and are hoping to apply it to the ICU.
Contrastive Learning We are investigating the use of contrastive learning to regularise the embedding space e.g. in(Yèche et al. 2021). Currently, there is no explicit loss toenforce relative positioning of the embeddings. Despite this,we have empirically found the clusters presented in this paper
to be very stable, both temporally and to encoder type. Thisis likely to be because our predictor and decoder networks
are very simple, meaning that the embedding space does nothave the freedom to model similar patient trajectories (whichshould be in the same cluster) in disparate parts of the embedding space. Nevertheless, contrastive learning could provide
further regularisation to reduce the chance of any unexpectedbehaviour.

1. The weighting of the duration task was not sufficient.
2. The tracheostomy task (but not the mortality task) reduces
the performance on the duration tasks.
The former does not seem likely, because as we alreadydiscussed, the Transformer is likely to be over-weighting theduration tasks, and yet, it follows the same trend as the LSTM
and TPC. The latter may appear to be counter-intuitive atfirst, especially in light of the correlation that exists between
the duration tasks and tracheostomy procedures. However,strong associations between tasks may not always lead tobetter performance. For example, one task can interfere withcrucial aspects of the representation required for solving theother. Usually this is an advantage of multitask learning,
because it enhances the signal:noise ratio when certain typesof noise only apply to one task. However, when there is
non-uniform correlation between tasks in the representation
space, it could harm the performance. Looking closely at
Figure 2, we can see that there is an area of patients near thebottom of the figure, in cluster 5. These patients have long
VD and LoS but have been separated from the other long staypatients in clusters 2 and 4. The separation can be attributed
to these patients never receiving tracheostomies, therefore thetracheostomy task forces the representation space to separate
these groups when they would be otherwise be aligned. Given
the simple nature of the predictor networks, this may harm
the performance on the duration tasks because the predictor
cannot effectively map these patients to appropriately long

Generalisability While we have shown that the clusters
are surprisingly stable, repeating the work on another dataset(e.g. MIMIC-IV or eICU) would strengthen this assessment.

Summary
While we acknowledge important limitations in our work, wehave shown that:
1. The TPC model outperforms alternative encoders on patient outcome prediction tasks.
2. It can also be used to generate clinically meaningful andinterpretable clusters using this technique.
3. The phenotypes are similar across choices of encoder andnumber of clusters.
4. The cluster assignment is remarkably stable over time, andmembership is determined early on. This is particularly
encouraging as a substrate for future intervention studies,because they rely on phenotyping before any intervention.
5. Stable cluster transitions do occur but they are infrequent.
Studying these transitions with a view towards understanding the causes is an important avenue for future work.
8

Acknowledgements

He, K.; Zhang, X.; Ren, S.; and Sun, J. 2015. Deep Residual
Learning for Image Recognition. CoRR, abs/1512.03385.

The authors would like to thank Petar Veličković, Sophie
Xhonneux, Stephanie Hyland and Mihaela van der Schaar for
helpful discussions and advice. We also thank the Armstrong
Fund, the Frank Edward Elmore Fund, and the School ofClinical Medicine at the University of Cambridge for theirgenerous funding.

Hochreiter, S.; and Schmidhuber, J. 1997. Long short-term
memory. Neural computation, 9(8): 1735–1780.
Ioffe, S.; and Szegedy, C. 2015. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. In Proceedings of the 32nd International Conference on International Conference on Machine Learning Volume 37, ICML’15, 448–456. JMLR.

References
Antcliffe, D. B.; Burnham, K. L.; Al-Beidh, F.; Santhakumaran, S.; Brett, S. J.; Hinds, C. J.; Ashby, D.; Knight, J. C.;
and Gordon, A. C. 2019. Transcriptomic Signatures in Sepsis
and a Differential Response to Steroids. From the VANISH
Randomized Trial. Am J Respir Crit Care Med, 199(8):980–986.
Bein, T.; Grasso, S.; Moerer, O.; Quintel, M.; Guerin, C.;
Deja, M.; Brondani, A.; and Mehta, S. 2016. The standard ofcare of patients with ARDS: ventilatory settings and rescue
therapies for refractory hypoxemia. Intensive Care Medicine,42(5): 699–711.
Bhavani, S. V.; Semler, M.; Qian, E. T.; Verhoef, P. A.; Robichaux, C.; Churpek, M. M.; ; and Coopersmith, C. M. 2022.
Development and validation of novel sepsis subphenotypes
using trajectories of vital signs. Intensive Care Med, 48(11):1582–1592.
Chami, I.; Gu, A.; Chatziafratis, V.; and Ré, C. 2020. From
Trees to Continuous Embeddings and Back: Hyperbolic Hierarchical Clustering. CoRR, abs/2010.00402.
Che, Z.; Purushotham, S.; Cho, K.; Sontag, D.; and Liu, Y.
2018. Recurrent Neural Networks for Multivariate Time
Series with Missing Values. Scientific Reports, 8(1): 6085.
Choi, E.; Bahadori, M. T.; Schuetz, A.; Stewart, W. F.; andSun, J. 2015. Doctor AI: Predicting Clinical Events via
Recurrent Neural Networks. JMLR workshop and conference
proceedings, 56: 301–318.
Corso, G.; Ying, R.; Pándy, M.; Veličković, P.; Leskovec, J.;
and Liò, P. 2021. Neural Distance Embeddings for BiologicalSequences.
Falcon, W. A. 2019. PyTorch Lightning. GitHub.
Famous, K. R.; Delucchi, K.; Ware, L. B.; Kangelaris, K. N.;
Liu, K. D.; Thompson, B. T.; Calfee, C. S.; and Network,
A. 2017. Acute Respiratory Distress Syndrome Subphenotypes Respond Differently to Randomized Fluid Management
Strategy. Am J Respir Crit Care Med, 195(3): 331–338.
Force, A. D. T.; Ranieri, V. M.; Rubenfeld, G. D.; Thompson,
B. T.; Ferguson, N. D.; Caldwell, E.; Fan, E.; Camporota,
L.; and Slutsky, A. S. 2012. Acute Respiratory DistressSyndrome: The Berlin Definition. JAMA, 307(23): 2526–
2533.
Girbes, A. R. J.; and de Grooth, H.-J. 2020. Time to stop
randomized and large pragmatic trials for intensive care
medicine syndromes: the case of sepsis and acute respiratory
distress syndrome. J Thorac Dis, 12(Suppl 1): S101–S109.
Harutyunyan, H.; Khachatrian, H.; Kale, D. C.; Ver Steeg,
G.; and Galstyan, A. 2019. Multitask Learning and Benchmarking with Clinical Time Series Data. Scientific Data,
6(96).

Kalchbrenner, N.; Espeholt, L.; Simonyan, K.; van den Oord,
A.; Graves, A.; and Kavukcuoglu, K. 2016. Neural MachineTranslation in Linear Time. CoRR, abs/1610.10099.
Kingma, D. P.; and Ba, J. 2014. Adam: A Method for Stochastic Optimization. CoRR, abs/1412.6980.
Lee, C.; and van der Schaar, M. 2020. Temporal Phenotyping
using Deep Predictive Clustering of Disease Progression.
Lin, M.; Chen, Q.; and Yan, S. 2013. Network In Network.
arXiv:1312.4400.
Lipton, Z. C.; Kale, D. C.; Elkan, C.; and Wetzel, R. C.
2015. Learning to Diagnose with LSTM Recurrent Neural
Networks. CoRR, abs/1511.03677.
Miotto, R.; Li, L.; Kidd, B. A.; and Dudley, J. T. 2016. DeepPatient: An Unsupervised Representation to Predict the Future of Patients from the Electronic Health Records. Scientific
Reports, 6(1): 26094.
Máca, J.; Jor, O.; Holub, M.; Sklienka, P.; Burša, F.; Burda,
M.; Janout, V.; and Ševčík, P. 2017. Past and Present ARDS
Mortality Rates: A Systematic Review. Respiratory Care,
62(1): 113–122.
Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.;
Chanan, G.; Killeen, T.; Lin, Z.; Gimelshein, N.; Antiga, L.;
Desmaison, A.; Köpf, A.; Yang, E. Z.; DeVito, Z.; Raison,
M.; Tejani, A.; Chilamkurthy, S.; Steiner, B.; Fang, L.; Bai,J.; and Chintala, S. 2019. PyTorch: An Imperative Style,High-Performance Deep Learning Library. In Advancesin Neural Information Processing Systems 32, 8024–8035.
Curran Associates, Inc.
Patel, A.; Montserrat, D. M.; Bustamante, C.; and Ioannidis,
A. 2022. Hyperbolic geometry-based deep learning methods
to produce population trees from genotype data. bioRxiv.
Poole, J.; McDowell, C.; Lall, R.; Perkins, G.; McAuley,
D. F.; Gao, F.; and Young, D. 2017. Individual patient data
analysis of tidal volumes used in three large randomized
control trials involving patients with acute respiratory distresssyndrome. BJA: British Journal of Anaesthesia, 118(4): 570–
575.
Rajkomar, A.; Oren, E.; Chen, K.; et al. 2018. Scalable andaccurate deep learning with electronic health records. Nature,1(1): 18.
Rocheteau, E.; Liò, P.; and Hyland, S. 2021. Temporal Pointwise Convolutional Networks for Length of Stay Predictionin the Intensive Care Unit. In Proceedings of the Conference
on Health, Inference, and Learning, CHIL ’21, 58–68. New
York, NY, USA: Association for Computing Machinery.
9Rocheteau, E.; Liò, P.; and Hyland, S. 2020. Predicting
Length of Stay in the Intensive Care Unit with Temporal
Pointwise Convolutional Networks.
Rusanov, A.; Prado, P. V.; and Weng, C. 2016. Unsupervised
Time-Series Clustering Over Lab Data for Automatic Identification of Uncontrolled Diabetes. In 2016 IEEE International
Conference on Healthcare Informatics (ICHI), 72–80.
Sheikhalishahi, S.; Balaraman, V.; and Osmani, V. 2019.
Benchmarking Machine Learning Models on eICU Critical
Care Dataset. arXiv:1910.00964.
Shickel, B.; Loftus, T. J.; Adhikari, L.; Ozrazgat-Baslanti, T.;
Bihorac, A.; and Rashidi, P. 2019. DeepSOFA: A Continuous Acuity Score for Critically Ill Patients using ClinicallyInterpretable Deep Learning. In Scientific Reports.
Song, H.; Rajan, D.; Thiagarajan, J.; and Spanias, A. 2018.
Attend and Diagnose: Clinical Time Series Analysis using
Attention Models. In 32nd AAAI Conference on Artificial
Intelligence, AAAI 2018, 4091–4098.
Srivastava, N.; Hinton, G.; Krizhevsky, A.; Sutskever, I.; andSalakhutdinov, R. 2014. Dropout: A Simple Way to Prevent Neural Networks from Overfitting. Journal of Machine
Learning Research, 15: 1929–1958.
Szegedy, C.; Liu, W.; Jia, Y.; Sermanet, P.; Reed, S. E.;
Anguelov, D.; Erhan, D.; Vanhoucke, V.; and Rabinovich,
A. 2014. Going Deeper with Convolutions. CoRR,
abs/1409.4842.
Thoral, P.; Driessen, R.; and Dam, T. 2020. AmsterdamUMCdb. GitHub.
Thoral, P. J.; Peppink, J. M.; Driessen, R. H.; Sijbrands, E.
J. G.; Kompanje, E. J. O.; Kaplan, L.; Bailey, H.; Kesecioglu,
J.; Cecconi, M.; Churpek, M.; Clermont, G.; van der Schaar,M.; Ercole, A.; Girbes, A. R. J.; and Elbers, P. W. G. 2021.
Sharing ICU Patient Data Responsibly Under the Society of
Critical Care Medicine/European Society of Intensive Care
Medicine Joint Data Science Collaboration: The Amsterdam
University Medical Centers Database (AmsterdamUMCdb)
Example. Critical Care Medicine, 49(6).
Tong, C.; Rocheteau, E.; Veličković, P.; Lane, N.; and Liò,P. 2022. Predicting Patient Outcomes with Graph Representation Learning, 281–293. Cham: Springer International
Publishing.
van den Oord, A.; Dieleman, S.; Zen, H.; Simonyan, K.;
Vinyals, O.; Graves, A.; Kalchbrenner, N.; Senior, A. W.; andKavukcuoglu, K. 2016. WaveNet: A Generative Model for
Raw Audio. CoRR, abs/1609.03499.
Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones,L.; Gomez, A. N.; Kaiser, u.; and Polosukhin, I. 2017. Attention is All You Need. In Proceedings of the 31st International Conference on Neural Information Processing Systems,
NIPS’17, 6000–6010. Curran Associates Inc.
Wang, Y.; Zhao, Y.; Therneau, T. M.; Atkinson, E. J.; Tafti,
A. P.; Zhang, N.; Amin, S.; Limper, A. H.; Khosla, S.; andLiu, H. 2020. Unsupervised machine learning for the discovery of latent disease clusters and patient subgroups using
electronic health records. Journal of Biomedical Informatics,
102: 103364.

Yèche, H.; Dresdner, G.; Locatello, F.; Hüser, M.; and Rätsch,
G. 2021. Neighborhood Contrastive Learning Applied toOnline Patient Monitoring. In Meila, M.; and Zhang, T.,
eds., Proceedings of the 38th International Conference on
Machine Learning, volume 139 of Proceedings of Machine
Learning Research, 11964–11974. PMLR.
Zhang, X.; Chou, J.; Liang, J.; Xiao, C.; Zhao, Y.; Sarva, H.;
Henchcliffe, C.; and Wang, F. 2019. Data-Driven Subtypingof Parkinson’s Disease Using Longitudinal Clinical Records:
A Cohort Study. Scientific Reports, 9(1): 797.

10

Further Analyses
Ablation Study
We performed an ablation study on the tasks used to train
the representation space. The results are shown in Table 6.
Firstly, we see that the best results for all tasks (except for theduration tasks) are achieved in the full multi-task setting. Nota single metric improves in the other ablation settings, and yet
at least one metric showed a deterioration in performance (theexception in task setting (g) is discussed below). Overall thisindicates that having multiple competing learning objectives
has a stabilising effect on learning the representation.

Figure 8: A comparison of the cluster assignments produced
by different encoders. We can see that there is strong cohesion
between some of the clusters. For example, cluster 3 appears
to be the same in all three models.

(c) – No Forecasting Experiment (c) included all the tasks
except forecasting one timestep ahead. When we compare experiment (c) to (a), we see that the results are mostly similar,
but there is a consistent decrease in performance, which is statistically significant at the p<0.05 level on the tracheostomy
task (AUPRC in the TPC model and AUROC in the Transformer model). On the reconstruction task, again the performance is similar but statistically worse in the last timestep
reconstruction in the LSTM model. This means that the forecasting task is contributing slightly to the performance in (a),but the benefit is small.

Figure 9: A comparison of the cluster assignments producedby TPC models which have been trained with different random seeds.
If we examine the TPC/LSTM comparison (far left in
Figure 8) we can see that clusters 2, 3 and 4 in the TPC correspond to 5, 3 and 1 respectively in the LSTM. In addition,clusters 1 and 5 (TPC) imperfectly map to 2 and 4 (LSTM)– the main difference being that some additional patients incluster 5 in the TPC map to cluster 2 in the LSTM. Thismeans that the k-medoids algorithm has placed a different
boundary between these groups in the clustering process.
Looking back to Figure 2, clusters 1 and 5 are revealed to beneighbours – in fact, cluster 5 appears to envelope cluster 1,
suggesting that 1 is a sub-cluster of 5. This is also consistent
with the clinical picture shown in Table 3, where the maindifference is that cluster 1 appear to have acute pulmonary
dysfunction, most likely in addition to other organ failures.
Therefore, the explanation for this discrepancy is that the
LSTM has a more generous threshold than the TPC for inclusion in its highest risk ‘early, life-threatening pulmonary
injury’ category (cluster 2).
In the TPC/Transformer comparison, the clusters largely
correlate, except that cluster 5 patients in the TPC have beenplaced into cluster 4 in the Transformer. Further investigation
revealed these to be patients with increased risk of receivinga tracheostomy (i.e. the patients which lay closest to thedecision boundary between the clusters).
The LSTM/Transformer comparison mirrors some of thecorrelations in TPC/LSTM but the mapping is less precise.
This could be because the models do not perform as well,making the representations less reliable. It could also bebecause the models have different affinities for the various
tasks, creating divergent biases in the representation space.
It is worth noting that in all three models, cluster 3 is themost distinct. This is unsurprising because it corresponds
to patients whose physiology is closest to ‘normal’. Therefore this group is the most homogeneous and can always be
identified easily.

(d) – No Reconstruction Experiment (d) removes both thetimestep t reconstruction and the static data reconstruction
tasks, but keeps the forecasting task. The effect size is larger
than in (c), but again is only statistically significant on the tracheostomy task. The forecasting task performs significantly
worse in the Transformer and LSTM models without thereconstruction.
(e) – Prediction Tasks Only Experiment (e) includes thebinary and duration prediction tasks, but no reconstruction orforecasting. The performance again deteriorates, particularlyon the tracheostomy task, we also start to see a more noticeable deterioration in the duration tasks, although this is not
yet statistically significant.
(f) – Binary Tasks Only Experiment (f) follows the trendof worsening performance as tasks are removed. This meansthat the mortality and tracheostomy tasks consistently benefit
from supplementary tasks which help to distinguish signal
from noise.
(g) – Duration Tasks Only Experiment (g) shows unexpected results; all of the models return better results when
only predicting LoS and VD. This is not what has been observed previously in multitask settings ((Rocheteau, Liò, andHyland 2021; Harutyunyan et al. 2019)). This is exploredfurther in the discussion.
However, overall the trend is such that the more tasks thatare included, the better the average results across tasks.

Cluster Reliability
Choice of Encoder Figure 8 compares the cluster assignments using different encoder models. It is encouraging that
there is a strong cohesion between some of the clusters, meaning that the models are picking out genuine and consistentpatterns in the data.

Retraining TPC We retrained the TPC model 5 times withdifferent initialisation generated by different random seeds
11Table 6: Prediction task results for the task ablation study. The full task setting from Table 2a has been repeated for ease ofcomparison. Various task ablations are compared to (a): (c) includes all tasks except for the forecasting task, (d) includes alltasks except for the reconstruction tasks, (e) includes only the prediction tasks, (f) is only the binary tasks, and (g) is only the
duration tasks. The colour scheme, metrics and statistical test comparisons are explained in the legend to Table 2.
In-Hospital Mortality
AUROC
AUPRC

Tracheostomy
AUROC
AUPRC

Length of Stay
MAD
MSLE

Vent. Duration
MAD
MSLE

TPC
(a) Transformer
LSTM

0.833±0.010
0.697±0.012
0.823±0.002

0.644±0.013
0.434±0.019
0.608±0.008

0.804±0.007
0.760±0.012
0.774±0.002

7.20±0.13
8.46±0.07
9.16±0.06

0.359±0.010
0.495±0.007
0.663±0.008

3.24±0.07
3.95±0.20
5.57±0.04

0.210±0.008
0.256±0.016
0.681±0.011

TPC(c) Transformer
LSTM

0.831±0.006
0.675±0.052
0.820±0.003

0.645±0.009
0.399±0.079
0.608±0.003

0.796±0.006 0.499±0.016† 7.24±0.120.743±0.011† 0.406±0.022 8.44±0.29
0.773±0.005 0.473±0.014 9.16±0.04

0.360±0.005
0.492±0.024
0.663±0.005

3.26±0.07
3.95±0.25
5.60±0.05

0.210±0.004
0.251±0.026
0.685±0.010

TPC(d) Transformer
LSTM

0.832±0.005
0.698±0.017
0.820±0.003

0.645±0.016
0.431±0.041
0.608±0.007

0.796±0.007 0.483±0.020† 7.28±0.09
0.743±0.008† 0.391±0.008 8.44±0.23
0.773±0.002 0.464±0.011 9.19±0.04

0.362±0.007
0.492±0.019
0.669±0.006

3.29±0.06
3.91±0.39
5.59±0.03

0.213±0.002
0.253±0.033
0.688±0.010

TPC(e) Transformer
LSTM

0.828±0.004 0.643±0.010
0.676±0.019† 0.410±0.034
0.819±0.005 0.604±0.013

0.798±0.005 0.480±0.020† 7.38±0.20
0.736±0.021† 0.383±0.026 8.67±0.27
0.773±0.002 0.475±0.008 9.20±0.04

0.367±0.020
0.509±0.024
0.669±0.008

3.24±0.07
4.12±0.22
5.61±0.04

0.212±0.012
0.268±0.017
0.691±0.012

TPC(f) Transformer
LSTM

0.823±0.006† 0.626±0.014† 0.793±0.002† 0.477±0.017† 0.669±0.036 0.373±0.048† 0.737±0.021† 0.400±0.038 0.817±0.003† 0.597±0.007† 0.767±0.003‡ 0.458±0.016 -

-

-

-

TPC
(g) Transformer
LSTM

-

Model

-

-

0.507±0.020
0.419±0.033
0.473±0.015-

6.99±0.10† 0.341±0.007† 3.08±0.09† 0.180±0.004‡
8.18±0.12‡ 0.472±0.012† 3.68±0.18† 0.224±0.009†
9.05±0.05† 0.644±0.006‡ 5.55±0.01 0.668±0.003†

and compared the resulting clusters to the original (Figure 9).
Overall, there is strong cohesion between the models, butsometimes there are shifts in the boundaries between neighbouring clusters in Figure 2. Especially between cluster 5(moderate-severe) and cluster 1 (severe), where some models(TPC 2, TPC 3, TPC 4) allow patients in cluster 5 to enter
their ‘cluster 1’ equivalent. Nevertheless, very distinct phenotypes are almost never mixed e.g. clusters 2 and 3, 1 and 4, or
1 and 3 (as defined by the TPC 1 model). As in the encoder
comparisons, cluster 3 is always well characterised.

Figure 10: Cluster labels with increasing number of clusters(from k=2 to k=7).
Number of Clusters The value of k was determined using
the elbow method. In Figure 10, we show how the clusters
would appear with increasing value of k. What is most striking is that each time a cluster is added, the new cluster either
inserts itself within an existing cluster, or it appears at theintersection between existing clusters. For example, as wemove from 2 to 3 clusters, the new cluster 3 is almost completely contained within the old cluster 1. This pattern of
sub-dividing an existing cluster generally continues until wereach 6 and 7 clusters, when the new cluster inserts itself
at the boundary between two or more old clusters. In otherwords, increasing the value of k does not completely shift
the position of all of the clusters, but rather it carefully subdivides them. The importance of this behaviour with increasing
value of k is discussed in the next section.

Data Preprocessing
Static Features
We extracted 14 static features from the admissions table(Table 7). Discrete and continuous variables were scaled tothe interval [-1, 1], using the 5th and 95th percentiles as
the boundaries, and absolute cut offs were placed at [-4, 4].
Binary variables were coded as 1 and 0. Categorical variables
were converted to one-hot encodings, with the exception of‘agegroup’, ‘heightgroup’ and ‘weightgroup’. These appear
as ordered categories e.g. [18-39, 40-49, 50-59, 60-69, 7079, 80+] for agegroup. We converted these to an ordered
set centred on 0, [-1, -0.6, -0.2, 0.2, 0.6, 1], to preserve the12Decay Indicators With the forward-filling method alone,
the model would not know whether a particular data point
was genuine or whether the data had been imputed. This isimportant because the sampling itself may be informative,
for example a deteriorating patient may have more frequent
investigations. To mitigate for this, we added ‘decay indicators’ to specify where the data had been imputed, and if it
had, how long it had been since the genuine measurement
was taken. The decay was calculated as 0.8j , where j is thetime since the last recording. This is similar in spirit to themasking used by Rocheteau, Liò, and Hyland (2021); Cheet al. (2018).

quantitative significance of each category.
Table 7: Static features used in the model. ‘Null Height’and ‘Null Weight’ were added as indicator variables to indicate when the height or weight were missing and have been
imputed with the mean value. We added the variables ‘Admission Count’ and ‘Ventilation Episode Count’ based on
previous admissions and ventilation episodes.
FeatureType

Source Table

Sex
Age Group
Height Group
Weight Group
Admission Count
Ventilation Episode Count
Urgency
Previous Ward
Specific Location in ICU
Physician Speciality
Weight Source
Height Source
Null Height
Null Weight

BinaryDiscrete
Discrete
Discrete
Discrete
Discrete
BinaryCategorical
Categorical
Categorical
Categorical
Categorical
BinaryBinaryadmissions
admissions
admissions
admissions

Additional Implementation Details
We tested three different encoder models to generate the embeddings. They were all trained as follows. Firstly, the timeseries are given to the encoder network which processes andthen combines them with the static features. These are thenpassed through a small two-layer pointwise convolution to
generate the embeddings (shown in green on Figure 1). These
are given to a predictor network, a reconstruction network
and a forecasting network.
The predictor network is one layer, with four outputs, corresponding to the four outcome tasks – tracheostomy, mortality,
LoS and VD. For the binary predictions, we apply a sigmoid
activation function to generate a prediction between 0 and
1 and for the duration predictions we apply an exponentialfunction. This is intended to help to circumvent a common
issue seen in previous models (e.g. Harutyunyan et al. (2019),
as they struggle to produce predictions over the full range of
durations when the data is very skewed) because it effectively
allows the upstream network to model log(LoS) instead of
LoS. The log(LoS) distribution is much closer to a Gaussian
distribution than the remaining LoS. No activation functionis placed on the outputs of the forecasting or time series reconstruction networks, because the variables are continuous.
Batch normalisation (Ioffe and Szegedy 2015) and dropout(Srivastava et al. 2014) is used to regularise the model.
The LSTM (Hochreiter and Schmidhuber 1997) is very
similar to the one used in a recent eICU benchmark paper
including LoS prediction (Sheikhalishahi, Balaraman, andOsmani 2019). The Transformer (Vaswani et al. 2017) is very
similar to its original implementation except that we added
temporal masking to impose causality6 (see ‘Hyperparameter
Search Methodology’ for their hyperparameters).

admissions
admissions
admissions
admissions
admissions
admissions

Time Series
For each ventilation episode, we selected 31 time series variables, mostly from the numericitems table (these are shownin Table 12). We used a semi-automatic process for feature
selection. To be included, the variable had to be present in
at least 25% of patient stays, and these were further narrowed down with advice from Dr Ari Ercole. We extracted‘diagnosissubgroups’ using a query from the AmsterdamUMCdb github repository (Thoral, Driessen, and Dam 2020) andventilator settings from listitems. The ventilator settings classification is given in Table 13. We engineered the features‘lung compliance’ and ‘P/F ratio’ because they are clinicallyimportant, and we have previously noted that neural networks
are unreliable when performing divisions. We calculated lungcompliance as:
0.73556 × Expiratory Tidal Volume
Peak Inspiratory Pressure − PEEP
(1)
where 0.73556 is a conversion factor to convert lung compliance to its usual unit of ml/cmH2 O. We calculated the P/F
ratio as:PaO2
P/F Ratio =(2)
FiO2
where FiO2 is expressed as a fraction rather than a percentage.
The time series variables were standardised in the samemanner as the static features. To help the model cope withthis missing data, we re-sampled according to one-hour intervals and forward-filled the data over the gaps. Note that
this is more realistic than interpolation as the clinician would
only have the most recent value. After forward-filling wascomplete, any data recorded before the ICU admission wasremoved.
Lung Compliance =TPC Model Architecture
We use a Temporal Pointwise Convolution (TPC) (Rocheteau,Liò, and Hyland 2020) network as one of our encoders. Thisis a model that takes advantage of both temporal convolution(to analyse trends) and pointwise convolution (to look forany important variable interactions). It is inspired by the waythat clinicians would approach an assessment of a patient e.g.
they might check how the respiratory rate is changing overtime, and they may also look at combination features e.g. the6
The processing of each timepoint can only depend on current
or earlier positions in the sequence.

13

PaO2 /FiO2 ratio. The components of the network are briefly
explained below.
Y+1T
F + Zt-1

POINTWISETEMPORAL

add flat features(repeated T times)& decay indicators

add padding
Y+1

Y+1

T

T T+ d(k - 1)
F + Zt-1

F + Zt-1

f

1

1
F

batchnorm,
dropout

TTY*batchnorm,
dropout

T
F + Zt-1
Z*T

repeat Y* times
Y*

Figure 12: Overview of the encoding framework. F, T, Y∗ ,Z∗ , Zt and f are defined in the caption to Figure 11. The
original time series (grey) along with the decay indicators(orange) (explained under ‘Time Series’) are processed by n
TPC layers. If a baseline model were used instead of TPC,
the time series output dimensions would be M x T, where M
is the LSTM hidden size or dmodel in the Transformer (this isin place of the light blue and grey output in the TPC model).
The diagnoses, d, are embedded by a diagnosis encoder –a single fully connected layer of size D. The time series(blue and grey), diagnosis embedding (purple) and staticfeatures (yellow) are concatenated along the feature axis, anda two-layer pointwise convolution is applied to obtain theembeddings (green).

T
F + ZtReLU,
add skip connections
Y* + 1T
F + ZtFigure 11: One layer of the TPC model. F is the numberof time series features. T is the time series length. Y is thenumber of temporal channels per feature in the previous TPC
layer (except for the first layer where Y is 1; decay indicators(explained under ‘Time Series’ in ) make up this channel).
Zt−1 is the cumulative number of pointwise outputs from all
previous TPC layers. Y∗ and Z∗ are the number of temporal channels per feature and pointwise outputs respectively
in the current TPC layer. Zt = Zt−1 + Z∗ . The differently
coloured temporal filters indicate independent parameters. d
is the temporal dilation, k is the kernel size. Decay indicatorfeatures () are shown in orange, f static features are shownin yellow. The skip connections consist of F original features(grey) and Zt+1 pointwise outputs (light blue). We ignore thebatch dimension for clarity.

data. Unlike most implementations, it does not share weights
across features i.e. weight sharing is only across time. Thisis because the features in a typical EHR differ sufficiently intheir temporal characteristics and warrant specialised processing. In TCNs, the receptive field sizes7 are highly adaptable.
They can be increased by using greater dilation, larger kernel
sizes or by stacking more layers. By contrast, recurrent neuralnetworks such as LSTM can only process one time step at
a time, and Transformers have a weaker sense of temporal
structure e.g. periodicity, which is central to understandingtime series in the EHR.
Pointwise Convolution Pointwise (or 1x1) convolution(Lin, Chen, and Yan 2013) is typically used to reduce thedimensions in an input (Szegedy et al. 2014). However in the
TPC model it is used to compute interaction features from

Temporal Convolution Temporal Convolution Networks(TCNs) (van den Oord et al. 2016; Kalchbrenner et al. 2016)
are models that convolve over the time dimension. The TPC
model uses stacked TCNs to extract temporal trends in the7‘Receptive field’ refers to the width of the filter. For TCNs thiscorresponds to a timespan.

14

the existing feature set at each timepoint.
Skip Connections Skip connections (He et al. 2015) allow
each layer to see the original data and the pointwise outputs
from previous layers. This helps the network to cope with
infrequently sampled data.
Temporal Pointwise Convolution The full model combines temporal and pointwise convolution in parallel. Figure 11 shows just one layer, however our implementation has
6 layers stacked sequentially (Table 8). With each successive
TPC layer, the temporal dilation is increased by 1.

Hyperparameter Search Methodology
All the encoders 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 8, 9 and 10). The hyperparameter search ranges havebeen included in Table 11. We ran 10 hyperparameter trials
to optimise the remaining parameters for the TPC, LSTM,
and Transformer models. The number of epochs was determined by selecting the best validation performance from amodel trained over 300 epochs (early stopping was then used
for each individual model). All deep learning methods were
implemented in PyTorch (Paszke et al. 2019) using PyTorch
Lightning (Falcon 2019) and were optimised using Adam
(Kingma and Ba 2014).
We also optimised for the weighting between the tasks. We
simply multiplied the loss for each component by a hyperparameter. The best overall learning curves were found when
the task weighting coefficients were: 0.5 for the duration
tasks, 1 for the binary tasks, 0.1 for time series reconstruction
and forecasting, and 0.002 for binary feature reconstruction.
The reason for the small weighting for binary feature reconstruction was that the task appeared very easy for the models,
especially predicting the sex of the patient, and so the representation became dominated with this at the expense of theother tasks.

Number of Clusters
The value of k was determined using an average value fromthe elbow method across various encoders. Specifically we
looked for the point at which the Within Cluster Sum of
Squares (WCSS) started to tail off with increasing values ofk. Figure 13 shows an example elbow plot. We selected thevalue 5 across all the models.

Figure 13: Elbow plot for the TPC model.
15

Table 8: The TPC model has 11 hyperparameters (Main Dropout and Batch Normalisation have been repeated in the table because
they apply to multiple parts of the model). We allowed the model to optimise a custom dropout rate for the temporal convolutions
because they have fewer parameters and might need less regularisation than the rest of the model. The best hyperparametervalues are shown in brackets. Hyperparameters marked with * were fixed across all of the models.

Temporal Specific

TPC Specific
Pointwise SpecificTemp. Channels (6)Point. Channels (14)Temp. Dropout (0.05)
Main Dropout* (0.05)
Kernel Size (3)Batch Normalisation* (True)
No. TPC Layers (6)Non-TPC Specific

Global Parameters

Batch Normalisation* (True)
Main Dropout* (0.05)
Final FC Layer Size* (16)Batch Size (128)
Learning Rate (0.0001)
Embedding Size (128)

Table 9: The LSTM model has 8 hyperparameters. We allowed the model to optimise a custom dropout rate for the LSTM layers.
Note that batch normalisation is not applicable to the LSTM layers. The best hyperparameter values are shown in brackets.
Hyperparameters marked with * were fixed across all of the models.
LSTM Specific

Non-LSTM Specific

Global Parameters

Hidden State (128)
LSTM Dropout (0.05)
No. LSTM Layers (2)Batch Normalisation* (True)
Main Dropout* (0.05)

Batch Size (128)
Learning Rate (0.0001)
Embedding Size (128)

Table 10: The Transformer model has 9 hyperparameters. Note that batch normalisation is not applicable to the Transformer layers(the default implementation uses layer normalisation). The best hyperparameter values are shown in brackets. Hyperparameters
marked with * were fixed across all of the models.
Transformer Specific

Non-Transformer Specific

Global Parameters

No. Attention Heads (2)Feedforward Size (256)dmodel (16)Transformer Dropout (0.05)
No. Transformer Layers (6)Batch Normalisation* (True)
Main Dropout* (0.05)

Batch Size (128)
Learning Rate (0.0001)

16

Table 11: Hyperparameter Search Ranges. We took a random sample from each range and converted to an integer if necessary.
For the kernel sizes (not shown in the table) the range was dependent on the number of TPC layers selected (because large kernel
sizes combined with a large number of layers can have an inappropriately wide range as the dilation factor increases per layer).
In general the range of kernel sizes was around 2-5 (but it could be up to 10 for small numbers of TPC Layers).
Hyperparameter

Lower

Upper

Scale

Batch Size
Dropout Rate (all)Learning Rate
Batch Normalisation
Final FC Layer Size
Point. Channels
Temp. Channels
LSTM Hidden State Size
dmodelFeedforward Size
No. Attention Heads
No. TPC Layers
No. LSTM Layers
No. Transformer Layers

400.0001
True
16
4
4
16
16
16
2
1
1
1

512
0.5
0.01
False
64
16
16
256
256
256
16
12
4
10

log2
Linear
log10
log2
log2
log2
log2
log2
log2
log2
LinearLinearLinearTable 12: Time Series features. The features which do not have a source table were calculated from the other features available inthe data. ‘Mandatory Ventilation’ and ‘Patient Triggered’ were calculated from the ventilator settings as outlined in Table 13.
FeatureType

Source Table

ABP gemiddeld
Ademfreq.
Alb.Chem (bloed)
Bilirubine (bloed)
CRP (bloed)
End tidal CO2 concentratie
Exp. tidal volume
Glucose (bloed)
Hartfrequentie
Ht (bloed)
Kalium (bloed)
Kreatinine (bloed)
Lactaat (bloed)
Leuco’s (bloed)
Natrium (bloed)O2 concentratie
P/F ratio
PC
PEEP (Set)
PO2 (bloed)
Piek druk
Saturatie (Monitor)
Temp.
Thrombo’s (bloed)
TroponineT (bloed)
UrineCAD
lung compliancemandatory ventilation
pCO2 (bloed)
pH (bloed)patient triggered
Time in the ICU

Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
Continuous
BinaryContinuous
Continuous
BinaryDiscretenumericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems

17

numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems
numericitems

Table 13: Ventilator Settings Classification, used to produce the features ‘Patient Triggered’ and ‘Mandatory Ventilation’ in
Table 12.
Patient Triggered Ventilation

Mandatory Ventilation

Bi Vente
NAVA
PRVC
PRVC (trig)
PS/CPAP (trig)
SIMV(PC)+PS
SIMV(VC)+PS
VC (trig)
VS
SIMV_ASB
CPAP
BIPAP-SIMV/ASB
MMV_ASB
MMV/ASBASB
IPPV/ASSISTCPPV/ASSISTCPPV_Assist
IPPV_Assist
SIMV/ASB
CPAP_ASB
PS/CPAPBIPAP/ASB
CPAP/ASB

MMV
VC
PC
Pressure Controled
PC (No trig)
PRVC (No trig)
VC (No trig)
CPPV
IPPV
SIMV
BIPAP

Table 14: Losses for the reconstruction tasks and forecasting task averaged over 5 independent training runs. The error margins
are 95% confidence intervals. See ‘Reconstruction and Forecasting’ for explanations of the losses shown. The meaning of (a),(b), the colour scheme and statistical tests are defined in the legend to Table 2.
Model

Reconstruction Tasks
Forecasting
Last Timestep Static (Binary) Static (Other)

TPC(a) Transformer
LSTM

0.334±0.004 0.013±0.000 0.210±0.038 0.334±0.005
0.351±0.005 0.013±0.000 0.354±0.005 0.347±0.001
0.297±0.006‡ 0.012±0.001† 0.078±0.010‡ 0.299±0.004‡

TPC(b) Transformer
LSTM

0.345±0.002‡ 0.013±0.000 0.332±0.006‡ 0.345±0.003‡
0.355±0.006 0.013±0.000† 0.356±0.001 0.353±0.005†
0.322±0.003‡ 0.012±0.000 0.266±0.004‡ 0.323±0.003‡

Table 15: Reconstruction and forecasting losses for the task ablation study. The full task setting from Table 14(a) has been
repeated for ease of comparison. The following task ablations are compared to (a): (c) includes all tasks except for the forecastingtask, (d) includes all tasks except for the reconstruction tasks. The colour scheme and statistical test comparisons are explained in
the legend to Table 2.
Model

Reconstruction Tasks
Forecasting
Last Timestep Static (Binary) Static (Other)

TPC(a) Transformer
LSTM

0.334±0.004
0.351±0.005
0.297±0.006

0.013±0.000
0.013±0.000
0.012±0.001

0.210±0.038
0.354±0.005
0.078±0.010

0.334±0.005
0.347±0.001
0.299±0.004

TPC(c) Transformer
LSTM

0.334±0.004 0.012±0.000
0.349±0.007 0.013±0.000
0.305±0.003† 0.011±0.000

0.198±0.020
0.358±0.007
0.085±0.010-

TPC(d) Transformer
LSTM

-

-

0.339±0.006
0.355±0.007†
0.309±0.006†-

18