\section{Abstract}
Coronavirus Disease 2019 (COVID-19) is an emerging respiratory disease caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) with rapid human-to-human transmission and a high case fatality rate particularly in older patients. Due to the exponential growth of infections, many healthcare systems across the world are under pressure to care for increasing amounts of at-risk patients. Given the high number of infected patients, identifying patients with the highest mortality risk early is critical to enable effective intervention and optimal prioritisation of care. Here, we present the COVID-19 Early Warning System (CovEWS{}), a clinical risk scoring system for assessing COVID-19 related mortality risk. CovEWS{} provides continuous real-time risk scores for individual patients with clinically meaningful predictive performance up to 192 hours (8 days) in advance, and is automatically derived from patients' electronic health records (EHRs) using machine learning. We trained and evaluated CovEWS{} using de-identified data from a cohort of \numprint{66430} COVID-19 positive patients seen at over 69 healthcare institutions in the United States (US), Australia, Malaysia and India amounting to an aggregated total of over \numprint{2863} years of patient observation time. On an external test cohort of \numprint{5005} patients, CovEWS{} predicts COVID-19 related mortality from $78.8\%$ ($95\%$ confidence interval [CI]: $76.0$, $84.7\%$) to $69.4\%$ ($95\%$ CI: $57.6, 75.2\%$) specificity at a sensitivity greater than $95\%$ between respectively 1 and 192 hours prior to observed mortality events - significantly outperforming existing generic and COVID-19 specific clinical risk scores. CovEWS{} could enable clinicians to intervene at an earlier stage, and may therefore help in preventing or mitigating COVID-19 related mortality.

\input{1_main.tex}

\section{Results}
\input{2_overview.tex}
\input{3_results_survival.tex}
\begin{rotatepage}
\input{4_cohort_table.tex}
\input{4_subcohort_table.tex}
\end{rotatepage}

\section{Discussion}
\input{5_discussion.tex}

\section{Data availability}
Accredited users may license the TriNetX COVID-19 research and Optum de-identified COVID-19 electronic health record databases used in this study  at TriNetX and Optum, respectively.


\section{Acknowledgements}
SP is supported by the Swiss National Science Foundation under P2BSP2$\_$184359. LAC is funded by the National Institute of Health through NIBIB R01 EB017205. BS is a member of the excellence cluster “Machine Learning in the Sciences” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC number 2064/1 – Project number 390727645. We thank Annika Buchholz for helpful discussions. 

\section{Author Contributions}
AM and PS created the new software used in the work. PS, AM, SP, LAC, JH, MH, SB were involved in the conception and design of the work, and the analysis and interpretation of results. JH, LAC, and MH provided the clinical motivation and interpretation. All authors were involved in reviewing, drafting and/or editing of the manuscript. PS and SB supervised the work.

\section{Competing Interests}
PS is an employee and shareholder of F. Hoffmann-La Roche Ltd.


\section{Introduction}
The coronavirus disease 2019 (COVID-19) pandemic has recently emerged as a major and urgent threat to healthcare systems worldwide. Since early reports of its outbreak in China in December 2019, the number of global cases has risen to over 21 million known infections and resulted in over \numprint{750000} deaths worldwide as of August 16, 2020 \citep{world2020coronavirus}. Despite public health efforts aimed at improving testing \cite{yan2020laboratory}, developing potential vaccines \cite{lurie2020}, and improving prevention strategies \cite{hellewell2020feasibility}, the disease is placing a significant burden on healthcare systems and existing resources in many countries, particularly where its spread has not been mitigated. Efficient early detection of patients likely to develop critical illness is thus crucial to optimise the allocation of limited resources, and monitor overall disease progression \cite{emanuel2020fair,fauci2020}. The use of clinical predictive models from electronic health records (EHRs) can help reduce some of this burden and inform better decisions overall \cite{henry2015targeted, wang2020response, green2018comparison, churpek2017investigating}. For instance, a model able to predict in advance which patients are at higher risk of mortality may help ensure resources are prioritized accordingly for these individuals. In addition, as more observational data are gathered, these models could be used both to discover new risk factors as well as reveal interactions between existing factors, offering better insights and opportunities for appropriate intervention.

Several approaches have been proposed to determine potential risk factors that contribute to COVID-19 mortality. Some of these approaches identify demographics and inflammatory markers associated with increased mortality \cite{Petrillim1966, Ji2020}, but do not account for risk factors potentially changing over time. Moreover, many existing analyses are limited to a single source of data, often from a single hospital, for both learning a model and predicting a patient's risk which may limit the generalisability of these analyses \cite{wynants2020prediction}. Other more traditional measures of patient prognosis such as Sequential Organ Failure Assessment (SOFA) scores \citep{vincent1996sofa} are based on examining a fixed set of risk factors not specifically adapted to COVID-19; such measures fail to account for relevant changes in patient status outside these risk factors, and therefore often do not reach high levels of sensitivity and specificity in identifying high-risk patients. Due to these challenges, to date, there does not yet exist a COVID-19 risk score that (i) makes use of multiple, representative sources of data to account for patient heterogeneity, (ii) includes important short-term and long-term risk factors that have a significant impact on mortality risk, (iii) reacts in real time to potentially rapid changes in patient status, and (iv) is adapted to risk factors relevant to COVID-19. 

To address these issues, we developed the COVID-19 Early Warning System (CovEWS{}), a risk assessment system for real-time prediction of COVID-19-related mortality that we trained on a large and representative sample of EHRs collected from more than 69 healthcare institutions using machine learning. In contrast to existing risk scores, CovEWS{} provides early warnings with clinically meaningful predictive performance up to 192 hours prior to observed mortality events, hence enabling critical time to intervene to potentially prevent such events from occurring. Since CovEWS{} is automatically derived from patient EHRs, it updates in real time without any necessity for manual action to reflect changes in patient status, and accounts for a much larger number of risk factors correlated with COVID-19 mortality than existing risk scores. CovEWS{} is based on a time-varying neural Cox model that accounts for risk factors changing over time and potential non-linear interactions between risk factors and COVID-19-related mortality risk\footnote{While these extensions have been pursued in~\cite{zhang2018time} and~\cite{ching2018cox} separately, they have neither been considered in combination nor in the context of COVID-19 risk scoring using EHRs.}, and was derived from the de-identified EHRs of \numprint{66430} diverse COVID-19 positive patients. We demonstrate experimentally that the predictive performance of CovEWS{} is superior to existing generic risk scores, such as SOFA \citep{vincent1996sofa}, COVID-19 specific risk scores, such as the machine learning models from \citet{yan2020interpretable} and \citet{liang2020early}, and COVER\_F \citep{williams2020seek}, and a time-varying Cox model with linear interactions \citep{cameron2020lifelines}. We additionally show that the gradient information of our differentiable CovEWS{} model can be used to quantify the influence of the input risk factors on the output score in real time. CovEWS{} may enable clinicians to identify high-risk patients at an early stage, and may, therefore, help improve patient outcomes through earlier intervention.






\subsection{Overview}
\label{sec:overview}
The overall pipeline of the method is shown in~\Cref{fig:overall_pipeline}. We refer to~\Cref{sec:nonlinear_time_varying} for a detailed presentation of the predictive model used by CovEWS{} and \Cref{fig:neuralnet_schematic} for a detailed diagram of the model architecture.

\subsection{Data Collection}
\label{sec:data}

We used data collected by two federated networks of healthcare organisations:

\paragraph{Optum.} The Optum de-identified COVID-19 electronic health records database includes de-identified electronic medical records and clinical administrative data including bedside observations and laboratory data from a geographically diverse set of healthcare institutions in the United States (US). The EHR data was sourced from more than 45 provider groups and integrated delivery networks. We used Optum cohort data collected between 21st March and 5th June 2020, and another cohort separated in time from 6th June to 13th July 2020 for our analysis.

\paragraph{TriNetX.} TriNetX is a global health research network providing a de-identified dataset of electronic medical records (diagnoses, procedures, medications, laboratory values, genomic information) including patients diagnosed with COVID-19. The data is de-identified based on standard defined in Section $\S$164.514(a) of the Health Insurance Portability and Accountability Act (HIPAA) Privacy Rule. The process by which Data Sets are de-identified is attested to through a formal determination by a qualified expert as defined in Section $\S$164.514(b)(1) of the HIPAA Privacy Rule. We used TriNetX cohort data collected between 21st March and 25th June 2020 from 24 healthcare organisations in the US, Australia, Malaysia and India for our analysis.

\paragraph{Data Quality.} Both Optum and TriNetX as well as the data providing healthcare institutions applied quality control steps to their data, but these procedures are not standardised neither across the federated networks nor across healthcare institutions. Varying levels of data quality across EHRs collected at different healthcare institutions and networks are therefore expected. However, heterogeneous data quality standards are characteristic for real-world data collected at different healthcare institutions. By evaluating CovEWS{} against an external test cohort from healthcare institutions with data collection policies different from our training cohort, we are able to give a fair assessment as to how robust and transferable CovEWS{} is in presence of realistic variations in data quality.

\begin{table}[t!] 
\setlength{\tabcolsep}{0.91em}
\caption{Comparison of the percentage of COVID-19 diagnosed patients without (w/o) a recorded SARS-CoV-2 test result across the analysed datasets (Optum, Optum Future, TriNetX). The higher fraction of COVID-19 diagnoses without a recorded SARS-CoV-2 test result in the Optum cohort compared to the Optum Future cohort is likely a result of the available testing resources having been more scarce early on in the pandemic. Since the Optum database indicates a generally high testing ratio substantiating COVID-19 diagnoses, the relatively higher percentage of COVID-19 diagnoses without corresponding tests in TriNetX is likely a result of inconsistent coding of SARS-CoV-2 tests early in the pandemic.}
\label{tbl:diagnoses_without_test_results}
\centering
\begin{small}
\begin{tabular}{l@{\hskip 0.85ex}rrr}
\toprule
 & Optum & Optum Future & TriNetX \\
\midrule
COVID-19 diagnoses w/o test results [\%] & 6.31 & 0.00 & 37.92 \\
\bottomrule
\end{tabular}
\end{small}
\end{table}

\paragraph{Inclusion Criteria.} We only included patients that were COVID-19 positive in our analysis. In both datasets, we considered patients COVID-19 positive if they either (i) were diagnosed with any of the International Statistical Classification of Diseases and Related Health Problems 10$^\text{th}$ revision (ICD-10) codes J12.89, J20.8, J40, J22, J98.8, and J80 together with B97.29\footnote{The listed criteria correspond to the Centers for Disease Control and Prevention (CDC) COVID-19 coding guidelines effective February 20, 2020 (see \url{https://www.cdc.gov/nchs/data/icd/ICD-10-CM-Official-Coding-Gudance-Interim-Advice-coronavirus-feb-20-2020.pdf}.}, or (ii) had a positive COVID-19 lab test result (\Cref{tbl:diagnoses_without_test_results}). For patients identified as COVID-19 positive via ICD diagnosis codes, we used the date of diagnosis as the reference diagnosis date for our analyses. For those patients identified as COVID-19 positive via a positive lab test, we used the date of the test sample collection as the diagnosis date. For patients with both a positive COVID-19 lab test and diagnosis, the available diagnosis date took precedence. For the subgroup of patients that were not hospitalised, we included all patients that were neither admitted to a hospital as inpatients nor an intensive care unit (ICU) at any point according to their EHRs. We note that it is possible that hospitals did not in all cases record inpatient hospital admissions and ICU admissions in their respective EHRs - which may  explain the observed non-zero rate of intubations in the non-hospitalised group. Membership in the Asian, Caucasian and Black or African American subgroups was mutually exclusive in the underlying EHR data model, and a patient could therefore only be assigned to one of the subgroups. In contrast, hispanic ethnicity was assigned in conjunction with any of the previous race categorisations.

\paragraph{Feature Selection.} We selected EHR-derived covariates for inclusion as input variables for CovEWS{} based on (i) previously published research on clinical risk factors for COVID-19 \cite{zheng2020risk,zhou2020clinical,richardson2020presenting}, and (ii) expert input from several medical professionals involved in the treatment of COVID-19 patients. In addition, we aimed to include both short-term and long-term risk factors for COVID-19 related mortality due to the continuous real-time evaluation of CovEWS{}. We present the list of all included model input covariates including their p-values in \Cref{tbl:feature_pvalues}, and their distributions across the datasets in \Cref{tbl:cohort_description}.

\begin{longtable}{l@{\hskip 4ex}p{33.5ex}r}
\caption{Descriptions, Logical Observation Identifiers Names and Codes (LOINC) codes (if available), and p-values for each input covariate used by CovEWS{}. The p-values were calculated using Wald's $\chi^2$ tests for effect (H$_0$ = no effect) of the coefficients in CovEWS (linear) corresponding to the respective input covariates on the Optum training cohort.}
\label{tbl:feature_pvalues}
\centering
\endfirsthead
\caption*{\Cref{tbl:feature_pvalues} cont.}
\endhead
\toprule
Covariate & Description & p-value \\
\midrule
Sex & - & $0.41$ \\
Age & - & $<0.005$\\
Weight & - & $0.92$ \\
Height & - & $0.80$\\
Body Mass Index (BMI) & - & $0.85$\\
Intubation & Whether or not the patient is intubated & $0.04$ \\
Temperature & Body temperature (LOINC: 8310-5 on TriNetX) & $0.89$ \\
SpO$_2$ & Oxygen saturation by pulse oximetry (LOINC: 59408-5) & $<0.005$\\
Heart rate & - & $0.75$\\
Respiratory rate & - & $0.32$ \\
Systolic blood pressure & - & $0.01$\\
Diastolic blood pressure & - & $<0.005$\\
Kidney disease & see \Cref{tbl:icd_codes} & $0.05$ \\
Ischemic heart disease & see \Cref{tbl:icd_codes} & $0.04$ \\
Other heart disease &see \Cref{tbl:icd_codes} & $0.01$ \\
Cerebrovascular disease &see \Cref{tbl:icd_codes} & $0.20$\\
Hypertension &see \Cref{tbl:icd_codes} &  $0.03$ \\
Diabetes &see \Cref{tbl:icd_codes} &  $0.15$\\
Hyperlipidemia &see \Cref{tbl:icd_codes} &  $0.17$ \\
Cancer & see \Cref{tbl:icd_codes} &  $0.65$ \\
Dyspnea & see \Cref{tbl:icd_codes} & $0.42$ \\
COPD & see \Cref{tbl:icd_codes} &  $0.18$ \\
Asthma & see \Cref{tbl:icd_codes} &  $0.80$ \\
Pulmonary embolism & see \Cref{tbl:icd_codes} & $0.58$ \\
Connective tissue disease & see \Cref{tbl:icd_codes} &  $0.84$ \\
Inflammatory bowel disease & see \Cref{tbl:icd_codes} & $0.99$ \\
Osteoarthritis & see \Cref{tbl:icd_codes} & $0.37$ \\
Rheumatoid arthritis & see \Cref{tbl:icd_codes} & $0.79$ \\
HIV & see \Cref{tbl:icd_codes} & $0.97$ \\
Smoking (never) & - & $0.81$ \\
Smoking (previous) & - & $0.61$ \\
Smoking (current) & - & $0.84$ \\
Smoking (unknown) & - & $0.85$ \\
White blood cells & White blood cell count (LOINC: 26464-8) & $<0.005$ \\
Neutrophil & Neutrophils per 100 leukocytes in blood (LOINC: 26511-6) & $<0.005$\\
Lymphocytes & Lymphocytes per 100 leukocytes in blood (LOINC: 26478-8) & $<0.005$ \\
Eosinophil & Eosinophils per 100 leukocytes in blood (LOINC: 26450-7) & $0.38$ \\
Basophil & Basophils per 100 leukocytes in blood (LOINC: 30180-4) & $0.33$\\
Platelets & Platelets [\#/volume] in blood (LOINC: 26515-7) & $0.28$ \\
C-reactive protein & C-reactive protein [mass/volume] in serum or plasma (LOINC: 1988-5) & $0.52$\\
hs. C-reactive protein & C-reactive protein [mass/volume] in serum or plasma by high sensitivity method (LOINC: 30522-7) & $<0.005$ \\
Procalcitonin & Procalcitonin [mass/volume] in serum or plasma (LOINC: 33959-8)& $0.46$ \\
Fibrin D-dimer & Fibrin D-dimer Fibrinogen Equivalent Units (FEU) [mass/volume] in platelet poor plasma & $0.13$ \\
Ferritin & Ferritin [mass/volume] in serum or plasma (LOINC: 2276-4) & $0.12$ \\
Cardiac troponin T & Cardiac troponin T [mass/volume] in serum or plasma (LOINC: 6598-7) & $0.42$ \\
Creatinine & Creatinine [mass/volume] in Serum or Plasma (LOINC: 2160-0) & $0.03$ \\
Lactate dehydrogenase & Lactate dehydrogenase [enzymatic activity/volume] in serum or plasma (LOINC: 2532-0, 14804-9) & $0.01$ \\
Gamma glutamyl transferase & Gamma glutamyl transferase [enzymatic activity/volume] in serum or plasma (LOINC: 2324-2) & $0.69$ \\
Aspartate aminotransferase & Aspartate aminotransferase [enzymatic activity/volume] in serum or plasma (LOINC: 1920-8) & $0.10$ \\
Creatine kinase & Creatine kinase enzymatic activity/volume] in serum or plasma (LOINC: 2157-6) & $0.60$ \\
Bilirubin & Bilirubin [mass/volume] in serum or plasma (LOINC: 1975-2) & $0.28$\\
Albumin & Albumin [mass/volume] in serum or plasma (LOINC: 1751-7) & $<0.005$ \\
Interleukin 6 (IL-6) & Interleukin 6 [mass/volume] in serum or plasma (LOINC: 26881-3) & $0.68$ \\
pH & pH of blood (LOINC: 2744-1, 2746-6) & $<0.005$ \\
PCO$_2$ & Carbon dioxide [partial pressure] in arterial blood (LOINC: 2019-8) & $0.02$ \\
PaO$_2$ & Oxygen [partial pressure] in arterial blood (LOINC: 2703-7) & $0.63$ \\
HCO$_3$ & Bicarbonate [moles/volume] in blood (LOINC: 1959-6, 1960-4) & $0.17$ \\
CO$_2$ & Carbon dioxide, total [moles/volume] in serum or plasma (LOINC: 2028-9) & $0.21$ \\
\bottomrule
\end{longtable}
\clearpage

\begin{figure}[t!]
\centering
\includegraphics[width=1.0\textwidth]{figs/test_frequency2.pdf}
\caption{Distributions of the number of observations per patient (y-axis, log scale, violin plots) for time-varying covariates (x-axis) in the Optum and TriNetX datasets for patients that have at least one observation of a given covariate. The percentage of patients with no observation for each covariate is shown in \Cref{tbl:cohort_description}.}
\label{fig:test_frequency}
\end{figure}

\paragraph{Data Characteristics.} The cohort statistics of the two datasets are presented in \Cref{tbl:cohort_description}, the ICD-9 and ICD-10 codes corresponding to the diagnoses shown in \Cref{tbl:cohort_description} are given in \Cref{tbl:icd_codes}, and the number of observations per patient for time-varying covariates for the two datasets is visualised in \Cref{fig:test_frequency}. As is characteristic for clinical data collected in real-world contexts, missing covariates are common in both datasets. Missingness in real-world EHR data is caused primarily by differences in laboratory testing guidelines, data collection practices, available testing resources and measurement devices between hospitals, and may in some cases depend on patient status and preferences of clinical staff. For example, \Cref{tbl:missingness_hospitalised} compares the missingness between the Optum test set and the non-hospitalised patient subgroup of the Optum test set. In contrast to traditional clinical studies, realistic missingness patterns in both the training and evaluation datasets are a desirable feature in the context of our study as CovEWS{} is designed to be deployed in clinical contexts with similar missingness, and therefore has to be trained and evaluated in the presence of missingness patterns seen across a representative range of heterogeneous hospitals. Covariates were mostly balanced across the Optum and TriNetX datasets. The primary differences were a higher observed mortality rate, and higher ratios of intubations, connective tissue disease, and rheumatoid arthritis in the TriNetX data compared to the Optum data. In addition, we note that the majority of admissions were recorded as being of the "Unknown" type in the TriNetX database. Since the large fraction of unknown admission entries limited potential admission outcome analyses, we reported hospital and ICU admission outcomes as not available for TriNetX (\Cref{tbl:cohort_description}). In compliance with the HIIPA Privacy Rule Section $\S$164.514(a), patients' exact dates of death were not available to protect patient privacy. In our analysis, we therefore imputed the last recorded EHR entry date as the reference date of death for deceased patients. The actual dates of death may have happened at a later point, and our performance estimates are therefore potentially underestimating actual predictive performance, since (i) correct predictions that happened later would mean CovEWS{} predicted sooner than we thought for that patient (which is generally harder, see \Cref{fig:baselines}), and (ii)  incorrect predictions of CovEWS{} may actually have been outside of the prediction time horizon. We believe this approximation of the exact date of death is therefore an acceptable trade-off, since underestimation of performance is not as much a concern as overestimation would be.

\begin{table}
\caption{International Statistical Classification of Diseases and Related Health Problems (ICD) codes corresponding to the disease classifications and symptoms used in our analysis. HIV = Human Immunodeficiency Virus, COPD = Chronic Obstructive Pulmonary Disease.}
\label{tbl:icd_codes}
\centering
\begin{small}
\begin{tabular}{l@{\hskip 4ex}p{35ex}p{15ex}}
\toprule
Disease classification & ICD-9 & ICD-10 \\
\midrule
Kidney disease & 585 - 586 & N18 - N19 \\
Ischemic heart disease & 410 - 414 & I20 - I25 \\
Other heart disease & 390 - 398, 401 - 405, 416 - 417, 420 - 429, 115.03, 115.04, 115.13, 115.14, 115.93, 115.94 & I27 - I52 \\
Cerebrovascular disease & 430 - 434 & I60 - I69\\
Hypertension & 401 - 405 & I10 - I15 \\
Diabetes & 249 - 250, 357.2, 366.41& E10 - E14\\
Hyperlipidemia & 272 & E78 \\
Cancer & 140 - 239 & C \\
Dyspnea & 786.1, 786.2, 786.8, 786.9 & R06 \\
COPD & 496, 491.21, 491.22 & J44 \\
Asthma & 493 & J45 \\
Pulmonary embolism & 415 & I26 \\
Connective tissue disease & 446, 710 - 711, 713, 725, 136.1, 279.8, 517.2, 728.5 & I30 - I36 \\
Inflammatory bowel disease & 555 - 556 & K50 - K51 \\
Osteoarthritis & 715 & M15 - M19 \\
Rheumatoid arthritis & 274, 712 - 714, 716, 719 & M05 - M14 \\
HIV & 42 & B20 - B24 \\
\bottomrule
\end{tabular}
\end{small}
\end{table}


\subsection{Data Normalisation}
\label{sec:data_normalisation}
The EHR data across both data sources used two different, but compatible, underlying data models consisting of recorded diagnoses, demographics, lab tests, procedures, medications and clinical observations. For our risk factors of interest, we converted records from both datasets into a unified data representation. We used ICD-9 and ICD-10 to extract diagnoses (\Cref{tbl:icd_codes}), Logical Observation Identifiers Names and Codes (LOINC) to extract lab tests, Current Procedural Terminology (CPT), and ICD-9 Clinical Modification (ICD-9-CM) and ICD-10 Procedural Coding System (ICD-10-PCS) to extract intubation events from the EHR records. For lab tests, we additionally normalised the unit of each category of lab tests to be the same for each measured record of that category.

\subsection{Stratification}
\label{sec:stratification}

We split the Optum cohort used for model development into training (50\%), validation (20\%) and held-out test folds (30\% of all patients) at random stratified by patient age, gender, presence of mortality events, presence of intubation events, presence of ICU admission and presence of a human immunodeficiency virus (HIV) diagnosis. We added HIV to the set of stratification covariates since its low prevalence could otherwise have led to imbalances in this risk factor across the folds. Stratification produced balanced cohorts across the three folds (\Cref{tbl:cohort_description}). The Optum training fold was used to train CovEWS{}, the validation fold was used to select the optimal hyperparameter configuration for CovEWS{}, and the held-out test fold was used in addition to the external TriNetX test cohort to evaluate the out-of-sample generalisation performance of CovEWS{}.

\begin{table}[t!] 
\setlength{\tabcolsep}{0.91em}
\caption{Comparison of missingness (in \%), i.e. the fraction of patients that did not have any record for a specific covariate in their EHR, in several important clinical covariates between the Optum test set and the non-hospitalised subgroup of the Optum test set. SBP = Systolic Blood Pressure, RR = Respiratory rate, WBC = White blood cells, CRP = C-reactive protein.}
\label{tbl:missingness_hospitalised}
\centering
\begin{small}
\begin{tabular}{l@{\hskip 2.85ex}cccccc}
\toprule
Cohort & SBP & RR & WBC & CRP & Albumin & Platelets\\
\midrule
Optum Test Set & 38.56 & 42.21 & 48.40 & 67.81 & 53.48 & 48.45 \\
Non-hospitalised &  54.00 & 59.57 & 64.67 & 81.10 & 69.02 & 64.72 \\
\bottomrule
\end{tabular}
\end{small}
\end{table}

\subsection{Preprocessing}
\label{sec:preprocessing}

Discrete covariates with $p$ different values were transformed into their one-hot encoded representation with one out of $p$ indicator variables set to 1 to indicate the discrete value for this patient. All continuous features were standardised to have zero mean and unit standard deviation using observed covariate distribution on the Optum training fold. Missing values of continuous covariates were imputed in an iterative fashion using multiple imputation by chained equations (MICE) \cite{white2011multiple}. After the preprocessing stage, continuous input features were standardised and fully imputed, and discrete input covariates were one-hot encoded. All preprocessing operations were derived only from the training fold, and na\"ively applied without adjustment to other folds and datasets in order to avoid information leakage. 

\subsection{Method}
\label{sec:method}

We adopt a variation of the widely used Cox proportional hazard model that is adapted to accommodate nonlinear and time-varying effects of covariates on the log-hazard function. In the following, the basics of time-to-event analysis that is the main subject of this paper is briefly presented. Then we touch upon the Cox proportional model for continuous-time covariates that is followed by the modifications we applied to this model to prepare it for the current work.


\subsubsection{Survival Analysis}
\label{sec:survival_analysis}
Survival analysis which is also known as Time-To-Event (TTE) analysis included a large body of work consisting of mathematical tools to give a statistical analysis of the time duration until a specified event occurs. In the current work, the event is defined to be the time when a patient dies.

An important tool in time-to-event analysis is~\emph{hazard function}. In discrete-time setting, (e.g. if times are given in specified periods) the hazard function is a conditional probability defined in discrete-time as
\begin{equation}
    \label{eq:hazard_function_discrete}
    h(t|\xb) = P(T=t|T\geq t; \xb),\; t=1, 2,\ldots
\end{equation}
that represents the risk of dying at time $t$ if the patient has survived until that time. The relevant covariates of the patients up to time $t$ are encapsulated in the vector $\xb\in\RR^d$. Age, sex, and lab tests are examples of such covariates that can take either binary or standardised real values after preprocessing. Intuitively, the hazard function captures the underlying dynamics of the transition of the condition of the patient from alive to dead. The larger $h$ is at time $t$, the more likely it is for the patient to die at time $t$.

Another useful function is called~\emph{survival function} that is denoted by $S(t)$ and in discrete-time defined as
\begin{equation}
    \label{eq:survival_function_discrete}
    S(t) = P(T>t) = \prod_{s=1}^t (1-h_s).
\end{equation}

\noindent{}Similar functions can be defined in the continuous-time regime. Let $T_c$ be the continuous survival time with the probability density function $f_c$ and cumulative distribution function $F_c$. Similar to~\Cref{eq:survival_function_discrete}, the continuous survival function represents the probability of surviving until time $t$ that is defined as

\begin{equation}
    \label{eq:survival_function_continuous}
    S_c(t|\xb) = P(T_c>t|\xb)=1-F_c(t|\xb).
\end{equation}

\noindent{}Likewise, the continuous hazard function is defined as
\begin{equation}
    \label{eq:hazard_function_continuous}
    h_c(t|\xb) = \lim_{\Delta t\to 0}\frac{P(t\leq T_c\leq t+\Delta t | T_c\geq t;\xb)}{\Delta t}.
\end{equation}
Notice that unlike discrete hazard function~\labelcref{eq:hazard_function_discrete}, the continuous hazard function~\labelcref{eq:hazard_function_continuous} is not a probability distribution and can take values larger than one.

The last useful function in continuous survival analysis is the cumulative hazard function
\begin{equation}
    \label{eq:cumulative_hazard_function}
    H_c (t) = \int_0^t h_c(u) du.
\end{equation}

\noindent{}The connection between these quantities can be simply derived:

\begin{align}
    h_c(t)&=\frac{f_c(t)}{S_c(t)}\label{eq:hazard_density_survival},\\
    S_c(t)&=\exp(-\int_0^t h_c(u) du) = \exp(-H_c(t))\label{survival_cumulative_hazard},\\
    f_c(t)&=h_c(t)\exp(-\int_0^th_c(u) du)=h_c(t) \exp(-H_c(t))\label{density_hazard_cumulative_hazard}.
\end{align}

Before introducing the simple yet flexible Cox model, we discuss a few important issues that must be taken into account in survival analysis.

\paragraph{Censoring.} What makes the survival analysis different from a simple regression from the covariates $\xb$ to $T$ -- observed duration up to the occurrence of the event -- is the concept of~\emph{censoring}. An observation is called censored -- or more precisely~\emph{right-censored} -- if its survival time has not been fully observed. There are several causes for a censored observation. For example, if a patient is not under observation when the event occurs or if the information of the patient is lost for some reason, only a lower bound to the time-to-event $T$ is observed that is the last time the condition of the patient is recorded. 

\paragraph{Discrete vs. Continuous.} Although time is a continuous physical quantity, in practice, it is measured at discrete points. Especially, in medical sciences, the condition of the patient is measured on a regular daily or bi-daily basis. This implies that even though the change of the covariates of a patient occurs at certain points of time, the exact time is not known. The transition point is only known up to the resolution of the measurement. We assume the time at which an event of interest occurs is denoted by $T$. As the resolution of the measurement is hours in the datasets used in the current work, $T=t$ refers to an event that occurs within the $t^{\rm th}$ hour after the patient is admitted to the hospital and its health condition is recorded.

\paragraph{Ties.} In the limited resolution measurement of time, some observations may have the same survival times, e.g., two patients die on the same day even though it is extremely unlikely that both die at the same moment. However, even in continuous time data, ties may occur which is a hint of underlying discrete sampling in time.


A major difference between continuous and discrete-time survival analysis is that the hazard function is a probability distribution in discrete settings while it can take any positive value in continuous settings. However, the traditional continuous-time approach can still be used for discrete event times especially when the measurements are equally spaced.


\subsubsection{Cox Hazards Model}
\label{sec:cox_hazard}
The most widely known model in the analysis of continuous survival time is {Cox's proportional hazard model}~\cite{cox1972regression} that parameterizes the hazard function as
\begin{equation}
    \label{eq:cox_proportional_hazard}
    h_c(t|\xb)=h_0(t)\exp(\betab\tran\xb),
\end{equation}
where $h_0$ is called the {baseline hazard} that is modulated by the effect of covariates via $\exp(\betab\tran\xb)$. Notice that in the traditional Cox model~\labelcref{eq:cox_proportional_hazard} the covariates $\xb$ are assumed constant over time. Consequently, the temporal variation of the hazard function is separated from the influence of the covariates.


\subsubsection{Time-varying Covariates}
\label{sec:cox_proportional_hazard_with_time_varying_covariates}
In many experimental settings, the assumption of time-invariant covariates in~\labelcref{eq:cox_proportional_hazard} does not hold. For example, many entries in the electronic health records such as heart rate, temperature, and blood measurements do not remain constant over the course of the hospitalisation of a patient. Therefore, the traditional Cox model~\labelcref{eq:cox_proportional_hazard} is extended to a time-varying setting by replacing $\xb$ in~\labelcref{eq:cox_proportional_hazard} with $\xb_t$ that is the measured covariates at time $t$. Assume a dataset consists of $N$ patients indexed by $n=1,2,\ldots,N$. As a notational convention, $\xb^{(n)}_t$ denotes the vector of the corresponding covariates to the patient $n$ at time $t$. 

If the Cox model holds and continuous events are observed, the following function called~\emph{partial likelihood} is maximized to estimate $\betab$:
\begin{equation}
    \label{eq:linear_time_varying_cox_partial_likeliood}
    \Lcal(\betab):=\prod_{i=1}^k\frac{\exp (\betab\tran \xb^{(i)}_{t_i})}{\sum_{j\in\Rcal(t_i)}\exp(\betab\tran \xb^{(j)}_{t=t_i})},
\end{equation}
where $t_1<t_2<\ldots<t_k$ are the ordered times at which the events occur and $\xb^{(1)}_{t_1}, \xb^{(2)}_{t_2}, \ldots, \xb^{(k)}_{t_k}$ are the corresponding set of covariates at those times. Notice that the equality of the superscript of the covariate vector $\xb^{(i)}_{t_i}$(patient's index) and the subscript of time $t_i$ emphasizes the continuous event times and the fact that at most one patient experiences the event at each time. For the moment, we assume time is continuous that results in distinct event times. The set $\Rcal(t_i)$ is the set of the patient's indices that are at risk at time $t_i$. Being at risk means they are alive and can potentially experience the event.


\subsubsection{Nonlinear Time-varying Covariates}
\label{sec:nonlinear_time_varying}
One clear limitation of~\labelcref{eq:linear_time_varying_cox_partial_likeliood} that is caused by the definition of the hazard function~\labelcref{eq:cox_proportional_hazard}, is the fact that the exponent of the modulating function $\exp(\betab\tran \xb)$ is a linear function of $\xb$. Hence higher order interactions among different dimensions of the covariate vector cannot be captured by this method. To improve the expressiveness of the model, we replace the linear function $\betab\tran\xb$ with a nonlinear function realised by a neural network. Let $\phi(\cdot;\thetab):\RR^d\to\RR$ be the function implemented by the neural network and parameterised by $\thetab$. Therefore, the hazard model is represented as

\begin{equation}
    h(t|\xb)=h_0(t)\exp(\phi(\xb_t;\thetab)),
    \label{eq:cox_hazard_nonlinear_body}
\end{equation}
where $h_0(t)$ is the baseline population-level hazard that is independent of the associated covariates to each patient. Time-varying covariates are transformed by the function $\phi(\cdot;\thetab)$ to log-hazard. The parameters $\thetab$ are learned via maximising the partial log-likelihood~\cite{cox1972regression}. Despite traditional Cox proportional hazard model where the gradient and Hessian can be computed analytically, here, we use automatic differentiation to compute gradients with respect to $\thetab$. The nonlinear function $\phi(\cdot;\thetab)$ is implemented as a $2-$layer multilayer perceptron -- see \Cref{sec:method} for a detailed description. The hazard function~\labelcref{eq:cox_hazard_nonlinear_body} estimates the instantaneous risk of death at each time for each patient. Integrating with respect to time and exponentiating the result gives the survival function defined as

\begin{equation}
    S(t|\xb_{0:t})=P(T>t|\xb_{0:t})=\exp\left(-\int_0^t h(u|\xb_u) du\right).\label{eq:survival_function_main}
\end{equation}
Notice that $\xb_{0:t}$ denotes the set of covariates until time $t$, meaning that, the probability of survival up to time $t$ depends on the history of the covariates. 

\noindent{}The partial likelihood~\labelcref{eq:linear_time_varying_cox_partial_likeliood} is re-written as
\begin{equation}
    \label{eq:nonlinear_time_varying_cox_partial_likeliood}
    \Lcal(\thetab):=\prod_{i=1}^k\frac{\exp (\phi(\xb^{(i)}_{t_i};\thetab))}{\sum_{j\in\Rcal(t_i)}\exp (\phi(\xb^{(j)}_{t=t_i};\thetab))}.
\end{equation}

\noindent{}To give an intuition of~\labelcref{eq:nonlinear_time_varying_cox_partial_likeliood}, observe that the partial log-likelihood that is computed by taking logarithm of the right-hand side of~\labelcref{eq:nonlinear_time_varying_cox_partial_likeliood} will consist of $k$ terms corresponding to $k$ observed events. The parameter vector $\thetab$ is perturbed such that the hazard increases for the covariates of a patient who dies at time $t_i$ while it decreases for the covariates of the patients who remain alive at $t_i$.

\subsubsection{Resolving Ties}
\label{sec:resolving_ties}
Even though we adopt a continuous-time approach due to the non-normalised parametric form of the hazard function~\labelcref{eq:cox_proportional_hazard} and the resultant partial likelihood~\labelcref{eq:nonlinear_time_varying_cox_partial_likeliood}, the ties can still occur as we work in hourly resolution. Hence, it is possible that two patients die at the same time. When an event occurs for two patients at the same time, the partial likelihood~\labelcref{eq:nonlinear_time_varying_cox_partial_likeliood} is not valid anymore. Several methods exist in the literature to break the ties and remove the ambiguity such as {average partial likelihood} \cite{cox1972regression} and~\emph{Berslow's method}~\cite{breslow1975analysis} that lives on two ends of a spectrum. The former takes average among all possible orders of the events that can break the tie. Hence, it is the most accurate method but computationally prohibitive. The latter gives a partial likelihood almost exactly like the original Cox likelihood by assuming that every ordering of tied events results in the same partial likelihood. This method gives a crude estimate but is easy to implement. A midway approach that we adopted in this work is called Efron's tie-breaker \cite{efron1977efficiency}. In this method, a weighted average likelihood of tied cases is subtracted from the denominator of \Cref{eq:nonlinear_time_varying_cox_partial_likeliood}. Efron's method gives good accuracy and is moderately easy to work with -- see \cite{efron1977efficiency} for details. 

\begin{figure}[t!]
    \centering
  \includegraphics[width=0.50\textwidth]{figs/illustrations/neuralnet.pdf}
    \caption{Schematic illustration of a neural network realisation of the nonlinear exponent of the time-varying Cox hazard function with $L=1$ hidden layers with $N_\text{dim}$ hidden units and one output layer.}
    \label{fig:neuralnet_schematic}
\end{figure}

\subsection{Algorithm Details}
\label{sec:algorithm_details}
Survival analysis by the Cox model is done via  maximum likelihood estimation, where the aim is to maximise the logarithm of~\labelcref{eq:linear_time_varying_cox_partial_likeliood} in the original Cox's proportional hazard model and~\labelcref{eq:nonlinear_time_varying_cox_partial_likeliood} in the nonlinear extension. In the original method with linear exponent both gradient $\partial \log \Lcal / \partial \betab$ and Hessian $\partial^2 \log \Lcal / \partial^2 \betab$ can be computed analytically. This is not the case for our proposed extension~\labelcref{eq:nonlinear_time_varying_cox_partial_likeliood} with nonlinear exponent. Instead of an analytical gradient, we use automatic differentiation to compute the gradient $\partial \log \Lcal/\partial \thetab$. Once the gradient is derived, an appropriate gradient-based method is used to perturb $\thetab$ in the direction that increases the partial likelihood.

As mentioned in~\cref{sec:nonlinear_time_varying}, the linear exponent $\betab\tran \xb$ in~\labelcref{eq:linear_time_varying_cox_partial_likeliood} is replaced with a nonlinear function $\phi(\cdot;\thetab)$. We use a neural network with $L$ hidden layers to realise this function. The employed network linearly transforms the input features to a $N_\text{dim}$-dimensional hidden layer. The transformed features are passed through a leaky rectified linear unit (LeakyReLU) \cite{xu2015empirical} nonlinear activation function. The hidden activations are then transformed by a linear transformation to a single node and finally passes through a tangent hyperbolic (${\rm tanh}$) activation function. In summary the network function can be represented as

\begin{equation}
    \phi(\xb;\thetab)=\rm{tanh}(W_2 (\rm{LeakyReLU}(W_1 \xb))),
\end{equation}
 where $\thetab=\{W_1, W_2\}$ and $W_i,\;i=1,2$ are the trainable weight matrices of the network (\Cref{fig:neuralnet_schematic}). We used Xavier's method~\cite{glorot2010understanding} to initialise the weights $\thetab$ of the model. To prevent overfitting, we additionally applied dropout with a dropout probability of $p_\text{dropout}$. In our PyTorch \cite{paszke2017automatic} implementation of CovEWS{}, we observed stable  convergence of our model using the Adam \cite{kingma2014adam} optimiser with a learning rate of $0.001$ for up to 100 epochs.

\begin{table}[t!] 
\caption{Hyperparameter ranges used for hyperparameter optimisation of CovEWS{} and CovEWS{} (linear). Comma-delimited lists indicate discrete choices with equal selection probability. Hyperparameters selected after hyperparameter optimisation (\Cref{sec:hyperparameter_optimisation}) are highlighted in bold.}
\label{tbl:hyperparameters}
\centering
\begin{small}
\begin{tabular}{l@{\hskip 5.25ex}l@{\hskip 30.5ex}r}
\toprule
& Hyperparameter & Range / Choices\\
\midrule
\parbox[t]{2mm}{\multirow{3}{*}{\rotatebox[origin=c]{90}{Linear}}} & & \\
&Regularisation strength $\lambda$ & 0.01, \textbf{0.1}, 1.0\\
& & \\
\midrule
\parbox[t]{2mm}{\multirow{5}{*}{\rotatebox[origin=c]{90}{CovEWS{}}}} & & \\
& Number of layers $L$ & \textbf{1}, 2\\
& Number of hidden units $N_\text{dim}$ & 64, \textbf{128}, 256\\
& Dropout percentage $p_\text{dropout}$ & 10\%, \textbf{20\%}\\
& & \\
\bottomrule
\end{tabular}
\end{small}
\end{table}

\subsection{Hyperparameter Optimisation}
\label{sec:hyperparameter_optimisation}

For the methods trained on the Optum training fold (CovEWS{} and CovEWS{} [linear]), we used a systematic approach to hyperparameter optimisation where each prediction algorithm was given a maximum of 15 hyperparameter optimisation runs with different hyperparameter configurations chosen at random without duplicates from predefined ranges (see \Cref{tbl:hyperparameters}). Out of the models generated in the hyperparameter optimisation runs, we then selected the model that achieved the highest specificity at greater than 95\% sensitivity on the validation set of the Optum cohort.

\subsection{Postprocessing and Calibration}
\label{sec:postprocessing}

After training CovEWS{} using the Optum training cohort, the predicted hazard for a patient with state $\xb$ is the hazard function $h(t|\xb)$ \cref{eq:cox_proportional_hazard} evaluated at $t=128 \text{ hours}$ ($\approx 5.33$ days) given the current patient covariates $\xb$ and under the assumption that patient covariates stay constant. To produce CovEWS{} scores, we additionally apply post-processing using a percentile transformation that converts $h(t|\xb)$ into the percentile of patient states in the Optum validation set that are assigned a lower $h(t|\xb)$ than the evaluated patient state $\xb$. We chose to output CovEWS{} scores in form of percentiles to aid in the clinical interpretation of CovEWS{} as a risk score relative to a representative set of reference states, and to discourage interpretation as a mortality probability. Interpretation of CovEWS{} scores as a mortality probability is difficult since the mortality risk of a patient depends on their uncertain future trajectory and the prediction horizon, and is influenced by clinical interventions that may be initiated in the future. As shown in \Cref{fig:timeline_1}, patients' states may change rapidly and frequently, and clinical interventions can significantly and positively alter the trajectory of patients. We also verified experimentally that, when interpreted as a probability of mortality, CovEWS{} scores overestimate the actually observed probability of death on the Optum and TriNetX test sets since patients' states may improve, due to intervention or otherwise, between the prediction time and the end of the prediction horizon (\Cref{fig:calibration_curves}; similar results with CovEWS{} [linear] \Cref{fig:calibration_curves_linear}). We, therefore, decided to instead output CovEWS{} scores as relative risk percentiles between 0 and 100 that discourages interpretation as a probability of mortality. To aid in the use of CovEWS{}, the following \Cref{sec:thresholds} outlines calibrated thresholds that can be used to maximise specificity at the desired target level of sensitivity for different prediction horizons.

\subsection{Thresholds}
\label{sec:thresholds}

\begin{table}[t!] 
\setlength{\tabcolsep}{0.162em}
\caption{Optimal thresholds of CovEWS{} scores to maximise specificity at greater than 85\%, 90\% and 95\% sensitivity (Sens.) for each prediction horizon as selected on their respective receiver operator characteristic (ROC) curves on the held-out Optum test set.}
\label{tbl:thresholds}
\centering
\begin{small}
\begin{tabular}{l@{\hskip 0.35ex}|ccccccccc}
\toprule
Sens. & 1 hour & 2 hours & 4 hours & 8 hours & 16 hours & 24 hours & 48 hours & 96 hours & 192 hours \\
\midrule
85\% & 61 & 62 & 56 & 54 & 51 & 48 & 45 & 39 & 36 \\
90\% & 44 & 44 & 42 & 41 & 39 & 38 & 38 & 32 & 27 \\
95\% & 34 & 34 & 33 & 31 & 28 & 27 & 27 & 22 & 19 \\
\bottomrule
\end{tabular}
\end{small}
\end{table}

A key question for clinical decision making is which threshold should be used for CovEWS{} scores to indicate severe risk, and potentially trigger an automated warning. To provide guidance in choosing the appropriate CovEWS{} score depending on the desired trade-off between sensitivity and specificity, we evaluated the optimal observed thresholds of CovEWS{} scores for various target sensitivity levels using their respective receiver operator characteristic (ROC) curves for each prediction horizon (\Cref{tbl:thresholds}). Optimal score thresholds to maximise specificity at high levels of sensitivity were between 61 and 36, 44 and 27, and 34 and 19 depending on the prediction horizon for target sensitivity levels greater than 85\%, 90\% and 95\%, respectively. We note that lower thresholds are necessary to achieve high sensitivity for prediction horizons farther in the future as patients' deterioration has to be identified earlier in its progression.

\subsection{Feature Importance}
\label{sec:feature_importance}

Highlighting the clinical risk factors that positively or negatively influenced CovEWS{} to output a certain score is of high utility as it enables clinical users to contextualise CovEWS{} scores, and, in some cases, these highlights could potentially even point towards opportunities for timely intervention. We utilised the differentiability of our prediction model as outlined in \Cref{sec:algorithm_details} to provide a real-time visualisation of the clinical covariates that are most important for CovEWS{} at any given time point (see \Cref{fig:timeline_1} for an example). To compute the importance scores at each time point, we used the Integrated Gradients (IG) \cite{sundararajan2017axiomatic} method that calculates relative importance scores $a_i \in (-100\%,100\%)$ for each input feature $x_{t,i}$ in the feature vector $\xb_t$ with $i\in[0 \mathrel{{.}\,{.}}\nobreak d-1]$ where $d$ is the number of input features. We used IG with the mean feature vector $\bar{\xb}_t$ across the Optum training set as a reference, calculated $50$ intermediate steps for each explained $\xb_t$, and normalised $a_i$ to the range of $(-100\%,100\%)$ by dividing each $a_i$ by the sum $\Sigma_{i=0}^{d-1}|a_i|$ of all feature attributions for $\xb_t$. To obtain a timeline of attributions as shown in \Cref{fig:timeline_1}, we calculate attributions $a_i$ whenever a change in patient status was recorded in the patient's EHR.

\begin{figure}[t!]
{\centering\textbf{\textsf{CovEWS{} Calibration}}\par\medskip}
\begin{subfigure}{.49\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_calibration.pdf}
\caption{Held-out Optum Test Set} 
\end{subfigure}\quad
\begin{subfigure}{.49\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_calibration_ext.pdf}
\caption{External TriNetX Test Set} 
\end{subfigure}
\caption{Calibration plots \cite{niculescu2005predicting} and counts of observed predictions of CovEWS{} on the held-out Optum test set (left) and the external TriNetX test set (right) when interpreting the predicted risk percentile of CovEWS{} as the probability of a mortality event being observed within the next 24 hours. The reference time point for those patients that did not have an observed mortality event is the date of their respective last observed EHR entry. When interpreted as a patient's mortality probability, CovEWS{} overestimates the risk since it can not account a priori for clinical interventions and potential future patient trajectory changes that may occur rapidly and frequently. Direct interpretation of CovEWS{} scores as a mortality probability is discouraged, and CovEWS{} scores should instead be seen as a relative risk score to stratify patients into risk groups (\Cref{fig:survival_curves}).}
\label{fig:calibration_curves}
\end{figure}

\begin{figure}[t!]
{\centering\textbf{\textsf{CovEWS{} (linear) Calibration}}\par\medskip}
\begin{subfigure}{.49\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_calibration_linear.pdf}
\caption{Held-out Optum Test Set} 
\end{subfigure}\quad
\begin{subfigure}{.49\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_calibration_ext_linear.pdf}
\caption{External TriNetX Test Set} 
\end{subfigure}
\caption{Calibration plots \cite{niculescu2005predicting} and counts of observed predictions of CovEWS{} (linear) on the held-out Optum test set (left) and the external TriNetX test set (right) when interpreting the predicted risk percentile of CovEWS{} (linear) as the probability of a mortality event being observed within the next 24 hours. The reference time point for those patients that did not have an observed mortality event is the date of their respective last observed EHR entry. Like the version of CovEWS{} that models non-linear interactions, CovEWS{} (linear) overestimates mortality risk when interpreted as a patient's mortality probability.}
\label{fig:calibration_curves_linear}
\end{figure}

\subsection{Baselines}
\label{sec:baselines}
In our analysis, we compared the performance of CovEWS{} to the following existing generic and COVID-19 specific clinical risk scores, and  baselines:

\paragraph{CovEWS{} (linear).} A linear time-varying survival Cox model as described in \Cref{sec:cox_proportional_hazard_with_time_varying_covariates} trained using the same Optum training set and using the same pipeline as the non-linear CovEWS{}. We used the implementation provided in version 0.24.8 of the lifelines \citep{cameron2020lifelines} Python package.

\paragraph{COVID-19 Estimated Risk for Fatality (COVER\_F).} The COVER\_F scoring system for COVID-19 as described in \citep{williams2020seek}. Since COVER\_F uses a single flag for any heart disease diagnosis, we aggregated all diagnoses in the diagnosis categories ischemic heart disease, pulmonary embolism, other heart diseases in our dataset into one single joint diagnosis code if any diagnosis in those three categories was present. All other input features used by COVER\_F were direct matches with the input covariates of the same name also used by CovEWS{} (\Cref{tbl:cohort_description}).

\paragraph{Sequential Organ Failure Assessment (SOFA).} SOFA scoring is commonly used in clinical contexts to indicate the risk of organ failure in critical patients. We, therefore, used SOFA as a generic risk scoring baseline that was not specifically designed for COVID-19 to demonstrate the comparative benefits in the predictive performance of a COVID-19 specific risk scoring system. Since we did not have FiO$_2$ values available in our EHR datasets, we assumed a default FiO$_2$ value of $21\%$ (equal to the fraction of oxygen in inhaled air) for patients that were not intubated, and an average of 71\% for patients that are intubated (FiO$_2$ is often set to 100\% initially and then progressively lowered as the patient stabilises, see e.g. \cite{rachmale2012practice} for an example). In addition, we did not have access to Glasgow coma scale (GCS) scores in the EHRs, and potential additional points from a high GCS score (a maximum of +4) were therefore not reflected in our calculated SOFA scores.

\paragraph{Yan et al. 2020.} \citet{yan2020interpretable} derived a simple and interpretable decision rule using three features (Figure 2 in \cite{yan2020interpretable}) for mortality prediction in COVID-19 patients from data collected from \numprint{485} COVID-19 positive patients seen in Wuhan, China. In their validation cohort, the decision rule showed a respectable cross-validated prediction performance of $96.1 \pm 0.03$ (mean $\pm$ standard deviation, 5-fold cross validation) \cite{yan2020interpretable}. All input features used by \citet{yan2020interpretable} were direct matches with the input covariates of the same name also used by CovEWS{} (\Cref{tbl:cohort_description}).

\paragraph{Liang et al. 2020.} \citet{liang2020early} developed a prediction model for critical COVID-19 related illness using data from \numprint{1590} patients seen at \numprint{575} medical centers in China using deep learning and 10 input covariates, including observed X-ray abnormalities. On three external cohorts from different Chinese provinces, they reported a predictive performance in terms of concordance index (c-index) of 0.890, 0.852 and 0.967 for predicting critical illness under the missingness of input covariates, respectively. Since we did not have access to radiologic assessments in our EHR datasets, we evaluated their model with the X-ray abnormality covariate missing for all evaluated patients (i.e. set to zero). To the best of our knowledge, \citet{liang2020early} did not specify which co-morbidities were included in their collected dataset. However, their study reports a maximum of 6 co-morbidities diagnosed in one patient. In our evaluation, we counted existing patient diagnoses of pulmonary embolism, kidney disease, inflammatory bowel disease, asthma, rheumatoid arthritis, and diabetes towards these 6 co-morbidities.

\subsection{Software}
\label{sec:software}

The source code used for developing CovEWS{} and for conducting the presented experiments and analyses was implemented using \texttt{Python} (version 3.7), \texttt{scikit-learn} (version 0.22.2), \texttt{numpy} (version 1.19.1), \texttt{scipy} (version 1.4.1), \texttt{pandas} (version 1.5.0), \texttt{PyTorch} (version 1.5.1), and \texttt{lifelines} (version 0.24.8). All plots shown were generated using the \texttt{ggpot2} R package \cite{wickham2016ggplot2}.

\subsection{Hardware}
\label{sec:hardware}

We used the high-performance computing (HPC) infrastructure provided by the Personalised Healthcare Informatics group at F. Hoffmann-La Roche Ltd to run the presented experiments. The compute nodes used 1st and 2nd generation Intel Xeon Platinum 8000 series processors and had access to  72 GB random access memory (RAM) each.

\subsection{Performance Evaluation}
\label{sec:performance}

In addition to the results presented in the main body of this work, we also present Receiver operator characteristic (ROC) curves for CovEWS{} for various prediction horizons between 1 and 192 hours evaluated on the held-out Optum test set (\Cref{fig:roc_optum}) and the external TriNetX test set (\Cref{fig:roc_trinetx}), the same ROC curves for CovEWS{} (linear) (\Cref{fig:roc_optum_linear} an \Cref{fig:roc_trinetx_linear}) a comparison of CovEWS{}, Time Varying Cox \citep{cameron2020lifelines}, COVER\_F \citep{williams2020seek}, SOFA \citep{vincent1996sofa}, \citet{yan2020interpretable}, and \citet{liang2020early} at various prediction horizons in terms of AUC, AUPR, F$_1$, sensitivity, specificity and specificity at greater than 95\% sensitivity (Spec.@95\%Sens.) for predicting COVID-19 related mortality on the held-out Optum test set (\Cref{tbl:results_all_optum}), on the external TriNetX test set (\Cref{tbl:results_all_trinetx}), on the Optum Future cohort (\Cref{tbl:results_all_optumfuture}), and on the Black or African American (\Cref{tbl:results_all_black}), Hispanic (\Cref{tbl:results_all_hispanic}), Asian (\Cref{tbl:results_all_asian}), Caucasian (\Cref{tbl:results_all_caucasian}) and non-hospitalised (\Cref{tbl:results_all_optumhospitalised}) subgroups of the Optum test set.

\begin{figure}[h]
{\centering\textbf{\textsf{CovEWS{} Receiver Operating Characteristic (Optum Test Set)}}\par\medskip}
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_1.pdf}
\caption{1 hour} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_2.pdf}
\caption{2 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_4.pdf}
\caption{4 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_8.pdf}
\caption{8 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_16.pdf}
\caption{16 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_24.pdf}
\caption{24 hours} 
\end{subfigure} \quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_48.pdf}
\caption{48 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_96.pdf}
\caption{96 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_192.pdf}
\caption{192 hours} 
\end{subfigure}\quad
\caption{Receiver operating characteristic (ROC) curves for CovEWS{} for various prediction horizons between 1 and 192 hours evaluated on the held-out Optum test set. The black dot indicates the optimal decision threshold for each prediction horizon selected on the Optum validation set as the closest point on the ROC curve to the top left coordinate (closest-to-top-left heuristic).}
\label{fig:roc_optum}
\end{figure}

\begin{figure}[h]
{\centering\textbf{\textsf{CovEWS{} Receiver Operating Characteristic (TriNetX)}}\par\medskip}
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_1.pdf}
\caption{1 hour} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_2.pdf}
\caption{2 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_4.pdf}
\caption{4 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_8.pdf}
\caption{8 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_16.pdf}
\caption{16 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_24.pdf}
\caption{24 hours} 
\end{subfigure} \quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_48.pdf}
\caption{48 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_96.pdf}
\caption{96 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_192.pdf}
\caption{192 hours} 
\end{subfigure}\quad
\caption{Receiver operating characteristic (ROC) curves for CovEWS{} for various prediction horizons between 1 and 192 hours evaluated on the external TriNetX test set. The black dot indicates the optimal decision threshold for each prediction horizon selected on the Optum validation set as the closest point on the ROC curve to the top left coordinate (closest-to-top-left heuristic).}
\label{fig:roc_trinetx}
\end{figure}


\begin{figure}[h]
{\centering\textbf{\textsf{CovEWS{} (linear) Receiver Operating Characteristic (Optum Test Set)}}\par\medskip}
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_1_linear.pdf}
\caption{1 hour} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_2_linear.pdf}
\caption{2 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_4_linear.pdf}
\caption{4 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_8_linear.pdf}
\caption{8 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_16_linear.pdf}
\caption{16 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_24_linear.pdf}
\caption{24 hours} 
\end{subfigure} \quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_48_linear.pdf}
\caption{48 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_96_linear.pdf}
\caption{96 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/test_roc_192_linear.pdf}
\caption{192 hours} 
\end{subfigure}\quad
\caption{Receiver operating characteristic (ROC) curves for CovEWS{} (linear) for various prediction horizons between 1 and 192 hours evaluated on the held-out Optum test set. The black dot indicates the optimal decision threshold for each prediction horizon selected on the Optum validation set as the closest point on the ROC curve to the top left coordinate (closest-to-top-left heuristic).}
\label{fig:roc_optum_linear}
\end{figure}


\begin{figure}[h]
{\centering\textbf{\textsf{CovEWS{} (linear) Receiver Operating Characteristic (TriNetX)}}\par\medskip}
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_1_linear.pdf}
\caption{1 hour} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_2_linear.pdf}
\caption{2 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_4_linear.pdf}
\caption{4 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_8_linear.pdf}
\caption{8 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_16_linear.pdf}
\caption{16 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_24_linear.pdf}
\caption{24 hours} 
\end{subfigure} \quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_48_linear.pdf}
\caption{48 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_96_linear.pdf}
\caption{96 hours} 
\end{subfigure}\quad
\begin{subfigure}{.31\textwidth}
\includegraphics[width=\linewidth,page=1]{figs/ext_roc_192_linear.pdf}
\caption{192 hours} 
\end{subfigure}\quad
\caption{Receiver operating characteristic (ROC) curves for CovEWS{} (linear) for various prediction horizons between 1 and 192 hours evaluated on the external TriNetX test set. The black dot indicates the optimal decision threshold for each prediction horizon selected on the Optum validation set as the closest point on the ROC curve to the top left coordinate (closest-to-top-left heuristic).}
\label{fig:roc_trinetx_linear}
\end{figure}











