The task of predicting hospitalisation outcomes from multivariate time-series is well-studied in the literature. Existing models range from ad-hoc early warning scoring systems based on aggregates of manually-selected physiological measurements [34] and machine learning models [8] . A prominent example of early warning scores is the National Early Warning Score (NEWS2) [44] , which has received formal endorsement by the United Kingdom's National Health Services to become the early warning system for identifying acutely ill patient. NEWS2 is representative of early warning scores in assuming independence among the predictive variables measured [22] , subsequently failing to capture the dependencies among the temporal signatures of the patient's physiology. In contrast, machine learning models overcome this limiting factor using sophisticated architectures to capture the non-linearities within the multivariate temporal data [8] . However, we find that existing machine learning approaches suffer from several problems. First, most of the current models have either been tested on a single condition (e.g. sepsis [23] , cardiac patients [26] , COVID-19 [57] , brain injury [43] ), or solely target Intensive Care Unit (ICU) settings [3] , [53] , where the magnitude of measurements is high, and the population is more uniform in the level of acuity (the interested reader can refer to [48] for a comprehensive review). Second, none of the existing models has been evaluated using metrics that accurately describe the model's ability to predict (the less frequent) outcomes (e.g. mortality). To clarify, consider the United Kingdom's inhospital mortality rates, which are around 23% in ICU settings [2] and 4% in secondary hospital wards [17] . If one thinks of outcome prediction as a binary classification problem, then it is undoubtedly one with a highly skewed distribution of the target variable (mortality in this case); with samples corresponding to survival being orders of magnitude larger than those corresponding to a positive outcome. In such a problem, the classifier's predictions should align with the user preference bias towards performance on poorly represented cases [5] . Despite this, the majority of the general outcome prediction models, i.e. those that are not constrained by a condition or patient population type [1] , [7] , [33] , [49] , generally rely on achieving a high Receiver-Operator Curve Area Under the Curve (ROC-AUC) without considering the metric's impartial assessment towards the imbalanced outcomes [13] . Finally, in contrast to medical practice, where the combined view of physiological sign changes and patient characteristics (e.g. demographics, pre-existing conditions) used to make a prognosis, the combination of multivariate temporal signatures and variables underlying patient characteristics has either been examined descriptively and not in the context of a predictive model [4] , or distinctly without consideration of the interplay between the two [47] . There is, therefore, a missed opportunity in developing forecasting models that combine dynamic timeseries with static features as done in other domains, e.g. [32] .
To overcome the above difficulties, we propose the reformulation of the task of outcome prediction from multivariate time series hospital data into one of outlier detection, whereby positive outcomes (e.g. mortality =1) are regarded as outliers. We subsequently develop a 2-level stacked machine learning architecture centred around an unsupervised LSTM-AutoEncoder module trained on the majority (negative outcome) samples, using it to learn a compressed representation of the temporal context that distinguishes those from 'outliers'. The learned context captures the multivariate interactions in the time-series. It is also used to guide a gradient boost module to estimate the risk of adversity using static data comprising summary statistics of vital signs and laboratory tests, demographics and admission details. The two modules, therefore, jointly make predictions using a combined view of a patient's state, based on the temporal signatures and static data. The stacked architecture is equipped with a justification module to generate visual justifications of the predicted outcomes.
LSTM AutoEncoders have shown competitive accuracy in detecting outliers from multivariate time-series data [28] , [29] . They do so by encoding the time-series in a low dimension to capture its most representative features. Since the encoded representation is compact, it is only possible to reconstruct representative features from the input, not the specifics of the input data, including any outliers [28] .These models have been effectively used in fall detection [40] , sensor failure prediction [36] , fraud detection [15] and video surveillance [55] . Despite their potential as solutions for healthcare problems, their use has been limited to retinal eye research [45] and fraud detection in healthcare settings [51] .
The problem of outcome prediction from hospitalisation data has recently witnessed the design of numerous ensemble architectures to enhance model performance and to support interpretability. Ensemble models have shown superior performance as compared to single-classifier architectures, either by enhancing predictive capability by different data modalities (i.e. a static view and a dynamic view) [47] , or by consolidating predictions of several 'weak' predictors in an additive way (e.g. voting) [54] . However, model combination in an ensemble is often performed during the prediction stage, which is problematic regardless of the type of ensemble used. In ensembles operating on different data modalities, the final prediction is not representative of the possible dependencies between the static view (e.g. demography, preexisting conditions, etc.) and the dynamic view (the temporal interactions embedded within the time-series) [32] . Similarly, ensembles of several weak classifiers have been shown to fail in alleviating the bias of the individual learners, and are generally outperformed by alternative models that stack strong classifiers into an ensemble [10] . These two observations are related to how existing models are evaluated in the literature. Specifically, all existing ensemble models unanimously use the Receiver-Operator Curve (ROC) space for evaluation [7] , [30] , [35] , [42] , [49] . However, the ROC space is known to be highly optimistic of classification performance when using highly-skewed datasets and where the interest lies in the correct prediction of the minority class [9] . Such conditions are representative of the problem of predicting (the less likely) adverse outcomes from Electronic Health Records (EHR) data.
The contribution of this work is the design and validation of an end-to-end machine learning framework to predict hospitalisation outcomes over increasing intervals. The framework possesses several attractive features that complement the literature in the following ways:
1) The framework is robust to the marginal representation of positive outcomes in a given population; this is achieved by regarding hospitalisation outcomes as outliers and making use of a successful outlier detection architecture as the base for prediction. 2) The framework captures the true interplay between the temporal signatures of a patient's physiology and static clinical summaries commonly used by practitioners in clinical settings. This is achieved via a stacked ensemble architecture whereby a compressed representation of the temporal signatures drives a classification model based on static features.
3) The framework is capable of justifying its predictions through visual illustrations. The justifications highlight the most contributing temporal signatures of the patient's physiology, as well as the static features, to the predicted outcome. 4) Evaluating the framework using PR-AUC and macroaveraged precision, recall and F1-score on real patient time-series data shows robustness over the diversity of patient populations (ICU data and secondary hospital wards; pneumonia, chronic kidney disease (CKD) and COVID-19 outcomes), outcome distributions, sample size, and short and long-term outcomes.
In this work, vectors are represented by boldface lower-case letters and matrices by boldface upper case letters. Unless specified, all vectors are assumed to be row vectors. The multivariate time-series input data
for n patients is observed as:
Where X d p contains the totality of the dynamic observations for patient p in the sample. Furthermore:
Where T is the number of observations per patient and each vector x p ∈ R v denotes the feature space consisting of a sequence of v vital signs or laboratory results dynamically accrued during the first 24 hours of admission. Furthermore, the static input data is D s = {x s p } n p=1 where x s p ∈ R u denotes the feature space consisting of a sequence of u static variables which may correspond to demography, pre-existing conditions or summary statistics of clinical variables. Finally,
is a column vector corresponding to the true incidence of outcome o during interval i for each patient p.
The goal of the framework is to predict:
is a binary variable representing the predicted onset of adverse outcome o during interval i for a single patient;p i o is the predicted probability of the outcomeŷ i o , which we retain for use during the interpretation stage. In this paper o ∈ {M ortality, ICU Admission, ICU ReAdmission} and i ∈ {5Day, 7Day, 14Day, 30Day}. Naturally, the class distribution ofŷ i o is highly imbalanced in favour of the negative outcome, as will be demonstrated in the evaluation sections.
In the framework presented here, the probability of outcome o during interval i,p i o , is learned using the two views of the clinical data D d and D s and is used to estimateŷ i o .
As the resulting time-series classification problem is one with imbalanced samples, the framework is designed to capitalise on the lower frequency of positive outcomes (e.g. ITU admission = 1) with respect to negative outcomes (e.g. no ITU admission) in the overall population. This is done by first 'self-training' using a subset of the negative samples of D d in order to learn a compressed representation describing the multivariate temporal dependencies within the majority (no adversity) population, enabling the evaluation of any incoming time-series for deviation from the majority (i.e. positive outcomes). The learned representation is then used to complement outcome prediction from the static view of the data D s .
The overall framework is shown in Figure 1 and comprises two modules stacked to predictp o i andŷ o i . Each of the two modules is designed to learn a view of the data: Dynamic-KD (KD: Knowledge Distiller) learns a uni-dimensional representationp i oDynamic ∈ [0 − 1] from the time-series such that for two individuals a and b whose time-series observations X d a and X d a ∈ D d , ifŷ i oa >ŷ i o b then :
That is, the vector representationp i oDynamic creates a separation between negative and positive-outcome time-series. The joint use of log transformation and flooring captures the difference in the order of magnitude of the resulting representations, rather than possibly non-significant fluctuations within the actual values [11] .
The second module, Static-OP (OP: Outcome Predictor) combinesp i oDynamic with D s to estimate the final prediction probabilityp i o . Static-OP is a classification ensemble based on gradient boost trees. The two modules form a bi-level stacked classification system we call KD-OP (Knowledge-Distillation Outcome Predictor). The output for the ensemble is:
Where γ ∈ [0 − 1] is a learned parameter corresponding to the optimal threshold for classification selected by optimising the mean Area Under the Precision-Recall Curve (PR-AUC) in the validation set.
KD-OP provides a visual justification of its predictedŷ i o for a given outcome o during interval i. The justification component combines the attention vector generated by Dynamic-KD and the relative importance of the static features generated by Static-OP along with the relative contribution of the two modules to the final prediction. The remainder of this section details the design of KD-OP's individual modules.
The multivariate time-series are first processed by Dynamic-KD, which consists of an unsupervised LSTM-Autoencoder architecture trained to reconstruct the input sequence by minimizing a distance-based objective function J . J measures the difference between the input vectors in the original series X d p and the vectors of the reconstructed seriesX d p for each batch (patient) p in D d . J is defined as below:
Where m is the number of multivariate observations for each patient p and . 2 is the L2-norm of a given vector.
Dynamic-KD adopts an attention mechanism over the time steps to capture the most important features in each sequence as proposed by [19] and successfully implemented in [46] , [49] , [50] . Figure 2 shows a feature-level representation of the attention mechanism in the encoder-decoder architecture of Dynamic-KD, reconstructing a multi-variable sequences over T time step batches (i.e. T ordered sequences per patient). For each feature j, a soft attention mechanism is implemented over the encoder's hidden states to obtain a distinct context vector c j . c j attenuates the most informative hidden states in s j,1 , ....s j,T of the decoder and is computed as follows:
For each feature j , the attention probabilities based on the encoded sequence α = (α 1 , ...., α T ) are calculated using the encoded sequence and the encoder's internal hidden states. First, the importance of information at each time step for feature j is calculated: Where U j and W | are trainable parameters capturing the input-to-hidden and hidden-to-hidden transitions for a given layer j respectively. Terms W j h t−1 and U j x t respectively capture the update from the hidden states at the previous step and the new input. a is the activation function.
In the decoder layers, we can measure the importance of the information at each time step for each feature j denoted by e j,t using proximity to U j . Then α j,t is obtained by normalising e j,t using the softmax operation:
T t=1 exp(e t ) Finally, the context vector for each feature c j is calculated using the weighted sum of the encoded sequence with the attention probabilities. Intuitively, this vector summarizes the importance of the encoded features in predicting t th sequence:
The Dynamic-KD LSTM-Autoencoder is trained in batches as shown in Algorithm 1. Training is done strictly on the negative-outcome data (the majority class), which we term D d T rain,0 (line 1), to minimise the reconstruction loss J (line 2). The resulting loss is therefore representative of the training errors associated with the negative (majority) class.
Receives: Training and validation subsets of the multivariate, regularly-sampled and batched time-series,
p=1 for z 2 patients in the validation set. Returns : Validation reconstruction lossp Dynamic , and attention matrix C V alid
The LSTM-Autoencoder is validated using a nonoverlapping subset of the time-series D d V alid (line 3), which contains mixed data (positive and negative outcomes), using the optimal loss obtained during training. Validation yields a reconstructionD d V alid of D d V alid and an attention matrix C V alid . At the end of the procedure, the difference between the original and reconstructed validation sets augments the original highly-dimensional feature space into a linear representation which is descriptive of the deviation from normality (no adversity) with respect to the temporal interactions embedded within the data. The validation reconstruction losŝ p Dynamic (line 4) is therefore discriminatory between the two classes and corresponds to the likelihoods of each batch (patient) p, where 1 ≤ p ≤ z 2 , belonging to the positive class. p Dynamic is used to complement the learning from static features performed by Static-OP, while C V alid is fed into the explainer component of the framework.
The goal of this module is to complement the predictions made through the temporal representation learned by Dynamic-KD using static information routinely employed by healthcare practitioners to assess a patient's risk factors (e.g. demographics, symptoms, summary statistics of physiology). In other words, instead of usingp Dynamic as a predictor of the outcome, it is instead used to drive further classification using static features via a gradient boosting model [27] . This way, the overall pipeline has the advantage of capturing the interplay between dynamic physiological measurements and static features in making the final predictions. The overall structure of Static-OP is given in Algorithm 2. First, Static-OP is trained using D s V alid and y V alid , with the reconstruction errorsp Dynamic serving as sample weights ω (line 1). Because Dynamic-KD ensures thatp Dynamic creates a separation between positive and negative classes, the minority samples of D s V alid will be the determinant of the decision threshold the model is trained to discover. Using the model to predict the outcome probabilities from y T est produces the predicted probabilitiesp, and will also produce the variable importance vector (line 2). The class labelsŷ, are obtained as in Equation 1 (line 4), using a prediction threshold learned by maximising the precision-recall area under the curve from the predicted probabilities (line 3). for z 3 patients in the test set. 2) True validation and testing class labels y V alid and y T est 3) Validation reconstruction loss,p Dynamic for z 2 patients in the validation set, obtained from Dynamic-KD.
We critically assess the model's performance in predicting mortality, as well as unplanned ICU admission and readmission via three different case studies: COVID-19 outcomes using general ward hospital data, and pneumonia and chronic kidney disease (CKD) outcomes using ICU data. The three use cases were chosen to represent different demographical and outcome distributions, as will be detailed in Section VI-A. The model is used to predict mortality in all studies but is only used to predict ICU admission in the COVID-19 dataset. Instead, since pneumonia and CKD use-cases are based on ICU time-series, the model is used to predict ICU readmission. The risks of all adverse outcomes are predicted at intervals of 5, 7, 14 and 30 days within hospital admission.
A. Datasets 1) COVID-19 Case Study : Data was collected from 1,276 adults (≥ 18 years old) inpatients of two acute hospitals: King's College Hospital and Princess Royal University Hospital in South East London, UK. All patients in the dataset tested positive for SARSCov2 between the 1 st of March and 31 st April 2020. Static data collected include age, gender, ethnic background, the length of the period from symptoms onset to hospital admission, and pre-existing conditions (specifically, chronic obstructive pulmonary disease (COPD), Asthma, heart failure, diabetes, ischemic heart disease (IHD), hypertension and chronic kidney disease). For training and risk prediction, pre-existing conditions were aggregated into one ordinal feature describing the number of comorbidities at the time of admission. The dynamic features included 14 routinely collected vital signs and laboratory tests and are available on our online repository 1 .
2) Pneumonia and CKD Case Studies: We used the data of ICU stays between 2001 and 2012 obtained from the anonymized Medical Information Mart for Intensive Care III (MIMIC-III) database, which is a freely-available anonymised ICU database and is the largest resource of time-series hospital data available worldwide [25] . We extracted admission details, demographics, time-stamped vital signs and laboratory test results obtained over the first 24 hours of admission of ICU stays of adults having ICD-9 code = 482.9 (pneumonia, cause not otherwise specified) and 585.9 (CKD, cause not otherwise specified) recorded as the primary diagnoses in the ICU admission notes. Since the MIMIC-III database is structured such as each hospital admission may correspond to multiple ICU stays, we extract the time-series pertaining to the first ICU stay of each admission, and used subsequent ICU admission to describe readmission outcomes. The resulting datasets comprise 509,323 records corresponding to 2,798 pneumonia ICU stays and 702,813 records corresponding to 2,822 CKD ICU stays (SQL and python scripts for recreating the dataset using the MIMICIII database are available on our online repository 2 ).
3) Data Description and Characteristics: Table I provides statistical summaries of the three datasets. The datasets vary in size, where pneumonia and CKD are much larger than COVID-19. The difference in size is a direct consequence of the mode of collection. The pneumonia and CKD datasets were extracted from the largest publicly-available ICU timeseries database [25] , while the COVID-19 data was locally collected over a short time span. Females were only the majority of cases in the COVID-19 dataset (females = 57.6 %), but gender distribution only significantly differed from the CKD dataset (females = 37.27%) and not from the pneumonia dataset (females = 45.58%). The pneumonia cohort was significantly younger and less co-morbid than the other two. The pneumonia cohort also showed a wider distribution of age compared to COVID-19 and CKD . In addition, the table shows that the number of pre-existing conditions varied greatly in the pneumonia and CKD cohorts, while the distribution of pre-existing conditions was more uniform in the COVID-19 dataset. The different distributions in age and pre-existing conditions is quite reasonable and align with the nature of the use cases: CKD is an age-related chronic illness [37] , with previous studies showing that the rate of comorbidities is around 41% [31] . COVID-19 hospital admissions are more likely in the elderly with pre-existing conditions such as hypertension and diabetes, where symptoms are likely to be more severe as opposed to the young healthy individuals [41] . In contrast, although both older age and pre-existing conditions increase the risk of acquiring pneumonia, they have not been found to be associated with the severity of the condition and subsequent intensive care needs [39] , [52] .
Across all prediction intervals, the COVID-19 dataset had higher rates of mortality, while mortality rates of CKD were significantly lower than the other two cohorts. In addition, the CKD cohort had significantly lower rates of ICU admissions across all prediction intervals compared to the other two cohorts. However, the time to ICU admission was much lower in COVID-19, where the average duration from admission to ICU admission was 4.35 days, compared to 12.23 days in pneumonia and 11.18 days in CKD. 
The three data extracts comprised irregularly-sampled time sequences of the vital signs. Using the extracted data to train KD-OP required transforming the irregular time series into D d = {X d 1 , ...X d n } for the n patients, where each X d p is a T × v matrix, T is a fixed number of observations and v is the number of vital signs measured at each time window t ∈ T . To alleviate missingness and to overcome non-uniform sampling, we deployed a patient-specific intervalbased aggregation of the multivariate time-series whereby the length of the aggregation window was chosen via an iterative process aiming to maximise completeness while minimising the length of the aggregation window. For each aggregation window, if a vital sign has multiple observation, then those are aggregated using a knowledge-based approach that mimics the summaries used by clinicians for each variable in practical settings and are available on our online repository 3 . The procedure produced T = 48 for pneumonia and CKD (aggregation into half-hourly intervals) and T = 12 for the COVID-19 use-case (aggregation into 2-hourly intervals). The resulting datasets comprised 30,624 samples with 12 variables for the COVID-19 dataset, 134,304 and 135,456 samples with 30 variables for the pneumonia and CKD cases respectively. The datasets were further imputed using a batch-aware Gaussian Process Regression model to interpolate the posterior mean and variance for the missing observations for each feature. The datasets have two distinguishing properties: 1) the outcomes are skewed, with positive outcomes being highly underrepresented in the time-series, and 2) the temporal ordering is defined over T batches, where each batch corresponds to a single patient. In order to retain the natural distribution of outcomes and temporal ordering during training and validation and to prevent information leakage, we used stratified grouped k-fold cross validation 4 , with k=3, to split the data for training, validation and testing sets as shown in Figure 3 . At each iteration, the data used to train Dynamic-KD was obtained using one fold, discarding the samples corresponding to patients with positive outcomes to yield D d T rain,0 . The second fold was used as a validation set of Dynamic-KD, in order to obtainŷ DynamicV al for the validation set and further to train Static-OP, usingŷ DynamicV al as sample weights. The third fold was first used to predict the testingŷ DynamicT est using Dynamic-KD and to test the performance of Static-OP, using the testingŷ DynamicT est as sample weights.
Furthermore, the data was normalised prior to training using a number of scaling techniques including absolutevalue scaling, standard scaling and min-max scaling. Using a min-max scaler with a range of [0 − 1] achieved the best classification performance and was therefore adopted for all the models. It should be noted that the data was scaled at each split as opposed to prior to splitting the data in order to ensure that the testing data is completely novel to the models. At each split, the training data was used to fit the scaler to obtain the shift and scale statistics used to further normalise the validation and testing data. 
We used the Python language and the Keras library with Tensorflow backend 5 . For Dynamic-KD, the hyperparameters used to train the model were optimised through empirical evaluation, by carefully observing the prediction performance using a set of candidate values of the hyperparameters; those included the number of neurons per layer, the number of hidden layers, dropout rates and the activation function used. The final design included bi-layered encoder and decoder, with the outmost layers having neurons in the order of 2 × n f eatures, where n f eatures is the number of dynamic predictor variables used (14 for COVID-19 and 30 in the pneumonia and CKD studies). A dropout rate of 0.5 was used between each two layers to prevent the autonecoder from overfitting the training data and an adaptive learning rate was used using the Adam optimizer and starting at 0.001. The number of epochs was 1,000, which was selected via cycles of experiments and careful monitoring of the validation loss. An early stopping criteria was used to retain the best model by minimising the validation loss with a patience parameter of 50 epochs. All layers of the autoencoder used ReLU as their activation function, which performed best during our evaluation. The Static-OP module was implemented using the XGBoost algorithm. The parameters were chosen through a grid-search over the hyperparameter space. Static-OP's sample weights were set to Dynamic-KD's prediction errors.
We evaluate KD-OP's performance across four dimensions. First, we evaluate the performance under different settings presented by the three datasets; these include cohort heterogeneity with respect to individual characteristics and outcome distribution with respect to the minority (positive) cases. Here, we initially report metrics averaged across the different prediction intervals for each setting to obtain an overall view, and subsequently evaluate the model's robustness across different prediction intervals. We then evaluate the contributions of the two modules Dynamic-KD and Static-OP to the overall performance, validating those empirically and against clinical knowledge. Finally, we compare the predictive power of KD-OP with existing outcome prediction models as reported in the literature. After evaluating KD-OP's performance, we demonstrate its visualisation capability in section VI-D. 5 .
Throughout the experiments, we report the Precision-Recall Area Under the Curve (PR-AUC) to capture the model's performance with respect to the minority cases, as well as the widely-used Receiver-Operator Area Under the Curve (ROC-AUC). Despite our knowledge of ROC-AUC's impartial assessment of the model's performance under positive and negative outcomes [20] , we choose to show it here due to its wide usage in the literature. Specifically, we use ROC-AUC to compare our model's performance with state-of-the-art models in section VI-D.4. We also report the macro-averaged precision, recall and F1-score. We used macro averages to understand the modules' true performance with respect to the under-represented outcomes [12] .
1) Overall Performance and Sample Diversity: We first evaluate the overall performance across the three case studies. For each dataset, Table II shows the model's performance averaged across the prediction intervals of 5, 7, 14 and 30 days for each outcome. As the table shows, the performance is high overall. However, better performance across prediction intervals was obtained using the COVID-19 dataset compared pneumonia and CKD, despite the latter two being larger datasets with a higher resolution of observations (half-hourly intervals as opposed to two-hourly intervals used to construct the COVID-19 time-series). A close evaluation is shown in Figure 4 , where higher performance ranges appear to be closely correlated with sample homogeneity (lower standard deviation) in age (shown in the legends) and the number of pre-existing conditions (reflected by the lightness of the plot circles).
COVID-19 admissions show a higher uniformity in age and pre-existing conditions, which consequently influences the patterns of changes in their physiological states. In both pneumonia and CKD, the diversity in the number of pre-existing conditions is starker (darker circles) in younger patients where lower-performance is achieved by the model.
Distribution Settings: Figure 5 shows KD-OP's performance in predicting mortality and ICU admission/re-admission over 5, 7, 14 and 30-day intervals on COVID-19, pneumonia and CKD. Two observations can be made when examining this figure in conjunction with the distribution of the outcomes of Table I . Apart from 5-day ICU re-admission in pneumonia and CKD, (a) KD-OP shows high performance across short and long-term intervals, and (b) KD-OP's performance remains high when the minority (positive) samples constitute < 10% of the overall population, which confirms the merit of relying on outlier detection to construct the temporal representation used in the pipeline. These findings are in-line with the demographic diversity results of Section VI-D.1, as the mean and standard deviations of age and number of pre-existing conditions of those re-admitted to the ICU within 5 days were and CKD populations as shown in table I. In contrast, the 5-day ICU admission sample in the COVID-19 study had a mean and standard deviation of age and number of pre-existing conditions being 63.3 (9.79) and 0.32 (0.05) respectively, showing a narrow range of demographical variation compared to pneumonia and CKD.
3) The Contribution of Static-OP vs Dynamic-KD: We now turn to compare the relative contribution of the two modules to KD-OP's overall prediction across the three use cases, outcomes and four intervals. The detailed comparison is provided in table III. In the table, we list the contribution of each module per prediction interval for each outcome using macroaveraged precision, recall and F1 score, as well as PR-AUC and ROC-AUC. We also show avg ∆, the average change in each metric's value between Dynamic-KD and the final prediction made by KD-OP. It is clear that the two modules complement each other to reach a high performance that is not otherwise achievable by the time-series predictor alone. This effect is especially noticeable in recall, where Static-OP significantly increases avg ∆. In mortality outcomes, the stacked model slightly decreases the precision of Dyanmic-KD, but the magnitude of the decrease (avg ∆) is insignificant compared to the increase in recall.
Examining the performance from a domain angle, Static-OP's contribution to the overall performance appears to be more pronounced in short-term outcomes. A highly noticeable difference is in the case of COVID-19 5-day mortality, where the average macro F1 score increases by 0.127 (from 0.786 using Dynamic-KD alone to 0.913 using the full pipeline). In contrast, the increase in F1 goes down to 0.021 (from 0.922 to 0.943) when examining 30-day mortality. This observation is consistent with current knowledge and recent findings that demographic information (e.g. age, pre-existing conditions) are highly predictive of short-term mortality in COVID-19 patients [56] . Similarly, for ICU readmission, replicated studies have found co-morbidities to be highly predictive of intensive care readmission during the same hospitalization [21] .
Here, we compare KD-OP's performance with the reported performance of relevant models and studies for outcome detection, showing the results in Table IV . Having gone through the literature, the only machine learning frameworks found to have been validated in nonICU settings are DEWS [49] , eCART [7] , and LightGBM [30] so we list those first. As the NEWS2 score is widely used to predict deterioration, we include the latest evaluation of its performance (from [22] ) in the table for comparison. For these models, we compare their performance against the average performance of KD-OP when applied to the COVID-19 use-case, as it is a general ward population. For each model, we highlight the class distribution of the target outcome as reported by each model's experimental evaluation settings. For KD-OP, the class distribution was taken as the average distribution of the outcomes over the intervals evaluated (5, 7, 14 and 30 days), as shown in Table I .
As none of the models reports PR-AUC, we compiled a list of all reported metrics, using them to compare against KD-OP's performance in conjunction with the widely-used ROC-AUC. For mortality, LightGBM achieved a high ROC-AUC of 0.961 with a class distribution of 5.1% of the outcome. However, the only other available metric for LightGBM is specificity, which is at a low 0.641 and entails a high rate Deterioration (14 days) -0.78 (0.74 -0.82) Guo [18] Deterioration (14 days) -0.67 (0.61 -0.73) Liu [24] Mortality & Deterioration -0.74 (0.69 -0.79) Galloway [14] Deterioration -0.72 (0.68 -0.77) Gong [16] Deterioration -0.853 (0.790-0.916) Overall, KD-OP shows the highest performance stability across the two outcomes, rendering it a better candidate for general hospitalisation outcome prediction; especially given the lack of thorough assessment of competitive models using metrics suitable for the problem under study. We also list high-performing machine learning models that have only been strictly validated in ICU settings; those include SANMF [35] , SICULA (a.k.a. the super learner) [42] and [38] . It is worth noting that none of the models predicts ICU readmission. We, therefore, resort to comparing with KD-OP's average performance in predicting mortality when applied to pneumonia and CKD using the MIMIC-III ICU dataset. Also, apart from [38] , which reports sensitivity, the models strictly rely on ROC-AUC in reporting their performance. We will, therefore resort to comparing with KD-OP's performance using ROC-AUC. As the table shows, KD-OP is the best predictor of mortality in an ICU setting, marginally exceeding SICULA's performance (ROC-AUC of 0.881 vs 0.880). Given that the SICULA's performance is the current benchmark for mortality prediction in the ICU, KD-OP's performance is wellaligned with existing prediction potential.
Finally, as the literature now contains several statistical models aiming to make prognostic predictions of COVID-19 hospital admissions, we compare those with KD-OP applied to the COVID-19 case. It is worth noting that all of the listed under the COVID-19 section of table IV are scoring systems aiming to mimic or exceed the performance of NEWS2 in predicting COVID-19 deterioration. Hence, KD-OP presents a novel contribution to the COVID-19 use case in being a scalable end-to-end machine learning architecture for predicting hospitalisation outcomes for COVID-19 admissionS. 5) Visual Justification of Predicted Outcomes: The stacked nature of KD-OP naturally enables visualising its predictions using the built-in visualisation properties of each module and obtaining the relative contributions of each module's predic-tion to the outcome. For Dynamic-KD, the feature attention weights at each time interval make up the relative importance of the temporal signatures of each feature. On the other hand, the gradient boost implementation of Static-OP provides a feature importance capability, which we use to understand the relative contribution of each static feature. Since the relative contribution of each module to the final prediction is outcome and interval dependent, including it in the visualisation of the output is highly essential for clinical utility as it directs the attention to the most contributing view (static or dynamic) of the patient. We define the contribution of each module using the ratios of the respective modules' PR-AUC.
An example of the generated visualisation of a positive 30day mortality outcome of a COVID-19 patient is shown in figure 6 . In the figure, the left-most bar shows the relative contribution of the individual modules. In this scenario, Dynamic-KD is a significant contributor (≈ 93% of the overall contribution). Examining the attention weights generated by the framework shows that the highest weights are of lymphocytes and neutrophils (Lymph and Neut in the figure) at hours 36-48 as well as C-reactive protein (CRP in the figure) at time-step 48 (24th hour). We use this information as a justification for the prediction made by Dynamic-KD. On the other hand, the patient's age and the mean NEWS2 score show the highest importance among the static features, followed by the length of the period from symptoms to admission (SxToAdmit) and the maximum C-reactive protein level over the 24 hours.
It is essential to view this justification in relation to current findings. C-reactive protein, lymphocytes and lactic dehydrogenase have been recently found to be highly correlated with adverse outcomes in COVID-19 patients [57] . Although lactic dehydrogenase was not part of our COVID-19 dataset, the temporal signatures of both C-reactive protein and lymphocytes have been accurately identified by Dynamic-KD as predictors of the particular patient's outcome. This, in addition to age being marked as an important static feature, agrees with recent findings [56] , showing that the signals jointly picked up by the framework's modules are coherent and well-aligned with clinical findings. 
We developed and validated KD-OP, an end-to-end pipeline for predicting adversity during hospitalisation. The pipeline comprises two stacked modules, each making predictions from a view of the patient's data: dynamic time-series and static features. The stacking of the pipeline's modules enables mimicking a clinician's approach to making prognostic decisions, by taking into account the interplay between the temporal signatures of a patient's physiology as well as time-invariant characteristics. By design, the pipeline is cognizant of the class imbalance natural to hospitalisation outcome data. It is trained and validated using stratified data that retains the original distribution of the real outcomes within the population. The pipeline's visualisation component complements its prediction by providing visual interpretations to aid clinical decision making. The visual interpretation provided by KD-OP is unique in that it accounts for the interplay between dynamic and static features in justifying the predictions made by the pipeline; a feature that derives directly from the stacked architecture. To our knowledge, this feature is not available in any existing hospitalisation outcome predictor.
We evaluated KD-OP's performance using real hospital data on three use cases representative of the diversity of electronic health records data. Using the pipeline to predict mortality and ICU admission/re-admission over 5-day, 7-day, 14-day, and 30-day intervals resulted in prediction accuracies exceeding 90% in all mortality outcomes and most of the ICU admission/re-admission outcomes.
A disease-agnostic model such as KD-OP could be built into the visual display of an EHR for all clinicians to use. The challenge at the moment is that each hospital department has its outcome prediction scoring system, subsequently making it unrealistic to build over 30 distinct models into an EHR system. The generic nature of KD-OP, coupled with high performance and visualisation capability, gives it a broader potential for integration in ICU and nonICU settings.
There are several possible avenues to build on the existing framework. First, it would be interesting to design a modified platform which projects the progression of the risk of adversity; we are currently developing a temporal risk score model to predict and visualise the risk of a given outcome on an individual level over time, using KD-OP as the base model. Second, the pipeline currently only supports classification outcomes, which limits its utility. Existing targets include the prognosis of continuous outcomes such as worsening oxygenation and worsening of cardiac function. Also, the current framework strictly uses routinely collected clinical variables as predictors. Other types of data can be of high relevance to a given use case. For example, ECG signals are the predictors of choice for cardiology-related outcomes; X-ray images can positively improve predictive power in the case of COVID-19, etc. Although the stacked architecture has proven to be highly robust compared to parallel ensembles, it is intrinsically less flexible towards extensions to incorporate additional models. It is, therefore, an interesting research problem to address, exploring the avenues of further extending the stacked model.
Our final word pertains to all outcome prediction models. The development of a useful clinical tool requires full engagement with stakeholders and ongoing clinical assessment of its validity. Numerous scoring systems have been developed over the decades, but few are routinely used because many are developed in isolation of the clinical teams. Therefore, strengthening academic and clinical collaboration is key to the success of any model.