{'fraction_non_alphanumeric': 0.04006150691952845, 'fraction_numerical': 0.013593029215786777, 'mean_word_length': 4.818442085172372, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SurvLatent ODE : A Neural ODE based time-to-event model with competing risks for longitudinal data improves cancer-associated Venous Thromboembolism (VTE) prediction

Effective learning from electronic health records (EHR) data for prediction of clinical outcomes is often challenging because of features recorded at irregular timesteps and loss to follow-up as well as competing events such as death or disease progression. To that end, we propose a generative time-to-event model, SurvLatent ODE, which adopts an Ordinary Differential Equation-based Recurrent Neural Networks (ODE-RNN) as an encoder to effectively parameterize dynamics of latent states under irregularly sampled input data. Our model then utilizes the resulting latent embedding to flexibly estimate survival times for multiple competing events without specifying shapes of event-specific hazard function. We demonstrate competitive performance of our model on MIMIC-III, a freely-available longitudinal dataset collected from critical care units, on predicting hospital mortality as well as the data from the Dana-Farber Cancer Institute (DFCI) on predicting onset of Venous Thromboembolism (VTE), a life-threatening complication for patients with cancer, with death as a competing event. SurvLatent ODE outperforms the current clinical standard Khorana Risk scores for stratifying VTE risk groups, while providing clinically meaningful and interpretable latent representations.


Introduction
Electronic Health Records (EHR), which consist of longitudinal measurements of patients' measurements often extending over decades, offer unique opportunities to build data-driven clinical decision support tools. Accurate stratification of future disease can allow clinicians to focus clinical resources on the most at-risk patients while minimizing potential side effects and overtreatment for low-risk patients (Myers et al., 2017;Maniruzzaman et al., 2018;Zeiberg et al., 2019). However, most existing patient stratification tools still rely on simple, outcome-specific scoring systems that utilize a small number of static features (Wilson et al., 1998;Khorana et al., 2008). As patients are often treated over the course of multiple interactions, a data-driven clinical decision model should relate the heterogeneous temporal evolution of clinical measurements to the clinical outcome in a flexible form and provide updated outcome prediction as a patient's condition and clinical features change. As EHR data is typically collected incidentally, the model should appropriately handle data missingness due to irregular spacing of visits or loss to follow-up, as well as informative missingness -where the collection of a measurement may itself be predictive of a future outcome. When estimating risk of a clinical outcome, an effective clinical decision model should also consider potential competing events a patient may experience and update the risk for the outcome of interest accordingly. This enables clinicians to jointly consider the risk of cancer-specific death or a treatment-related toxicity when planning a therapeutic intervention.
One widely utilized risk prediction framework is a time-to-event (or survival) analysis, which models the lifetime of a patient until a failure event. Importantly, time-to-event analysis can directly model outcomes that are lost to follow-up (also known as right-censoring), which is very common in EHR due to the lack of formal patient recruitment. One of the most widely utilized time-to-event models is the Cox proportional hazard model, a multivariable, semi-parametric framework proposed by Cox (1972), which relates log of relative hazard to a linear function of baseline features. Recently, there have been efforts to augment the Cox proportional hazard model by learning a non-linear relationship of baseline features to outcomes through deep neural networks (Katzman et al., 2018;Nagpal et al., 2021b). However, these models make time-to-event predictions based on a single snapshot of features and fail to capture evolution of time-varying features, often informative of predicting disease onset. To address this gap, time-to-event methods capable of handling time-varying, longitudinal measurements have been proposed. For example, some models utilized class of Recurrent Neural Network (RNN) models to learn input temporal representation which then renders estimated time-to-event predictions (Ren et al., 2019;Lee et al., 2019;, while other models utilized a deep generative approach to model the event time distribution in a Bayesian framework (Ranganath et al., 2016;Miscouridou et al., 2018). However, presented with missing measurements, these models rely on strategies agnostic to latent time-dependent dynamics between observations and may lead to ill-defined latent representation.
In a multivariate time-series framework, many approaches have been developed to handle informative missingness to improve prediction performance (Choi et al., 2016;Lipton et al., 2016;Futoma et al., 2017;Soleimani et al., 2017;Che et al., 2018). One approach to deal with missing values involves deploying a decay mechanism on both input data and RNN hidden states for describing dynamics under missing observations (Che et al., 2018). Other methods utilized a Gaussian process (GP)-based framework to learn dynamics across observed measurements and provide continuous latent states (Futoma et al., 2017;Soleimani et al., 2017). However, these models rely on parametric assumptions on the latent dynamics, specified by either exponential decay functions or GP kernel functions with a stationary property and may not necessarily describe true latent dynamics. On the other hand, a recently proposed family of Neural Ordinary Differential Equations (Neural ODE) based models (Chen et al., 2018;Rubanova et al., 2019) learns the latent dynamics using highly flexible neural network parameterized functions and is capable of handling input data which may arrive at arbitrary time points, providing an avenue to model longitudinal data with missing values without strong parametric assumptions on the latent dynamics.
In this work, we present a generative, time-to-event model, SurvLatent ODE, which utilizes the Neural ODE framework to effectively learn temporal dynamics of the input representation under irregularly sampled measurements. Adopting a multi-task learning framework, our proposed model allows the underlying mechanism of multiple events to be shared in the latent representation while utilizing a cause-specific decoder module to flexibly learn signals specific to each event from the shared latent representation. As a result, the proposed model provides survival functions for each patient across multiple clinical outcomes they may encounter while incorporating temporal evolution of their features. Using MIMIC-III, a publicly available longitudinal dataset collected at critical care units, we demonstrate that the proposed model significantly outperforms conventional as well as state-of-the-art time-to-event models in predicting time to hospital mortality. Furthermore, we utilize the longitudinal data from the Dana-Farber Cancer Institute (DFCI) for predicting time to VTE with all-cause mortality as a competing event, and show that the proposed model significantly outperforms the current clinical standard, Khorana scores (Khorana et al., 2008). Finally, we demonstrate the learned latent representation offers interpretable clusters of patients with meaningfully different outcomes. The implementation of SurvLatent ODE is available at https://github.com/itmoon7/survlatent_ode.

Generalizable Insights about Machine Learning in the Context of Healthcare
• To the best of our knowledge, our proposed model is the first demonstration of the ODE-based variational autoencoder time-to-event model for longitudinal data where temporal latent dynamics in the input data are explicitly modeled via neural networks.
• Combined with the effective longitudinal modeling, our proposed framework enables flexible estimation of hazard functions for the event of interest as well as competing events via a multi-task learning framework, which gains significant improvements over conventional and recently published deep learning based survival models.
• Applied to the in-house dataset of Venous Thromboembolism (VTE) events, our model significantly outperforms current clinical standards, Khorana scores (Khorana et al., 2008) and provides insights into influential features for elevated VTE risks via interpretable latent representations.

Related Work
Deep learning based approaches have been widely utilized in many published work in the context of time-to-event analysis to learn complex, non-linear relationship between features via neural networks (Lee et al., 2018a;Katzman et al., 2018;Nagpal et al., 2021b), Convolutional Neural Network (CNN) (Jarrett et al., 2018)   incorporates time-varying features from longitudinal data and estimates conditional survival distribution through a fixed mixture of parametric distributions like Weibull or Log-Normal. However, a parametric assumption on the underlying time-toevent process may lead to a model miss-specification and limit the flexibility of relating neural network-learned representations to the conditional survival estimates. Furthermore, the model only handles a single risk, which may lead to an overestimation of the disease risk since in healthcare applications there is often a varying degree of dependence across competing events (Berry et al., 2010;Austin et al., 2016). Lee et al. (2019) proposed Dynamic-Deephit, a discrete time-to-event model, which incorporates time-varying features and utilizes a multi-task learning framework to estimate joint distribution of the first hitting time and competing events. A key assumption of the model is that a patient experiences an event over the predefined finite time horizon with a probability of 1, which often does not hold in real data. When handling missing values, both models rely on population-level statistics (e.g. means) and/or missing data indicators (Lipton et al., 2016), which does not consider dynamics of underlying patient-specific health trajectory. Ranganath et al. (2016) and Miscouridou et al. (2018) demonstrated time-to-event models where missing observations were handled by a shared latent process which models observed measurements as well as event times in a Bayesian framework. However, their models relied on strong exponential parametric assumption on the underlying data generating process and were not evaluated in the longitudinal setting.
In the context of more conventional time-to-event frameworks for competing risks, causespecific Cox regression and Fine-Gray Cox regression (Fine and Gray, 1999) models are widely used for the analysis of continuous event times. These two models have been extended to the analysis involving event times measured on a discrete time scale (Lee et al., 2018b;Berger et al., 2020). While the Fine-Gray model enables estimation of covariates' effects on the cumulative incidence function for the event of interest, the cause-specific Cox model captures effect of covariates on the cause-specific hazard function for each event, which denotes the instantaneous rate of the corresponding event occurrence for patients who are currently event free (Putter et al., 2007;Austin et al., 2016). Both models assume linear relationship between patients' covariates and log of the relative hazard as well as proportional hazard where effects of covariates on the relative hazards remain constant over time. However, these assumptions are limiting because true underlying time-to-event processes are often described by complex nonlinear relationships between patients' biomarkers and effects of biomarkers on the survival may change over time as patients' health status are changing.
Ordinary Differential equations parameterized by neural networks, popularized by Chen et al. (2018), have been utilized in time-to-event analysis in a multi-state setting (Groha et al., 2020) and a single-event continuous-time setting (Tang et al., 2022). Groha et al. (2020) proposes a Neural ODE approach to estimate the Kolmogorov forward equations which then provides transition probabilities in a multi-state framework. Tang et al. (2022) models the distribution of survival time for a single event by learning the dynamics of the cumulative hazard function via neural networks. However, both approaches only considered the baseline data and left out informative signals from time-varying features. To the best of our knowledge, our proposed model is the first demonstration of full ODE-based encoder-decoder architecture for modeling longitudinal data in a time-to-event framework.

Preliminaries
We first introduce three key modeling choices of the proposed framework : Neural Ordinary Differential Equations (Neural ODEs), discrete time-to-event analysis, and competing risks.

Neural Ordinary Differential Equations
First, we model a patient's health trajectory using Neural Ordinary Differential Equations (Neural ODEs). Neural ODEs are a recently proposed family of neural networks-based continuous time models (Chen et al., 2018), which parameterize the dynamics of a hidden state h(t). This hidden state may represent the latent health trajectory of a patient and a function of their time-varying features (i.e. h(t) = z 1 (X t ), where z 1 (·) may be neural networks and X t is a set of patient's features up until time t). Given that the function f θ (parameterized by the neural network) specifies the dynamics of the hidden state, we can define the hidden state h(t) (e.g. a latent health trajectory) at any arbitrary time t as follows: where h(t 0 ) is an initial hidden state. Rubanova et al. (2019) demonstrated that this Neural ODE framework can be incorporated into a Recurrent Neural Network (RNN) model, allowing f θ to learn the dynamics of the RNN hidden state h(t) such that h(t) is well-defined in a continuous time. Such an ODE-RNN model can incorporate irregularly sampled input data without having to explicitly impute missing values. Finally, the hidden state trajectory may be further decoded into, for example, the hazard function trajectory for a time-to-event prediction.

Discrete time-to-event analysis
Second, we model the outcome event times as a discrete time process. A common constraint of healthcare data is that event times such as time of hospital death and time of stroke diagnosis are often recorded on a discrete time scale, for example, in hours or days. In this case, the exact time of an event is not known and all we know is that the event occurred between two consecutive time points. Formally, a continuous time horizon is divided into discrete intervals (i.e. [0, t 1 ), [t 1 , t 2 ), ...), where the time resolution is often determined by a domain knowledge, and an event time is denoted by T , where T = t means that the event of interest has happened in the interval [t − 1, t) (Tutz et al., 2016).

Competing risks
Third, we seek to model multiple competing outcomes (e.g. cancer-specific death, metastasis, and VTE) that a patient is under risk for simultaneously. One common approach to handle competing risks in a discrete time-to-event analysis is to model the cause-specific discrete hazard function (i.e. rate of the event occurrence in a current time interval for event free patients in the previous time interval) (Austin et al., 2016). At a patient level, this framework provides insights on the individual event mechanism by modeling the contribution of a patient's features on the hazard function over time. Formally, given a set of b different events K = {1, ..., b} and a set of features X, the cause-specific discrete hazard function λ k (t) for event k is where k ∈ K. And the resulting overall discrete hazard function is An overall event-free survival probability S(t|X) is Finally, the cause-specific cumulative incidence function (CIF) for event k, F k (t|X), which captures the risk for event k at time t under the other competing risks conditioned on X, is

Notation and setting
In brief, our model takes as input a set of individuals with (possibly irregularly spaced) time-varying features and event times of one or more outcome events, and learns the functions that map the features to each event specific hazard over time. We formalize these variables below.
Survival dataset D with sample size N is a set of tuples where t i and k i are the observed survival time and the event type for sample i, respectively. δ i indicates whether the event occurred for sample i. Note that t i = min(T i , C i ) whereT i is the true event time and C i is the right censoring (i.e. loss to follow-up) time of sample i. Therefore, for right-censored samples (i.e. those with δ i = 0), t i corresponds to C i . Note that throughout this paper, we assume that the censoring mechanism is non-informative.
In the competing risks scenario, there may be multiple events where an occurrence of one of those events prevents the occurrence of the other events (for example, illness cannot occur after death). To that end, given the finite set of b mutually exclusive competing events, X i is the set of features for sample i longitudinally measured at some irregular timestamps denoted by τ i . In other words, given that τ i,0 and τ i,l are the initial measurement time and latest measurement time for sample i, respectively, and M is the number of features. Let T r be a random variable for the remaining time-to-event from the latest measurement or pre-defined landmark such as treatment start date. Our goal is to estimate the discrete cause-specific hazard function denoted by where Z t i is the latent trajectory representing the data X i up until time t. For the rest of our paper, we use λ * i,k (t) to denote λ k (t|X i , Z t i ). Finally, we can relate λ * i,k (t) to the overall event-free survival function, S(t|X i ), as well as cause-specific CIF, F k (t|X i ), to assess individual risk of experiencing event k.  Fig. 1, SurvLatent ODE encodes a patient-specific temporal trajectory of features with a varying time length, denoted by X i , into the latent embedding with a fixed dimension via the ODE-RNN encoder. Then, it decodes the patient specific embedding into the latent trajectory with pre-specified length t m (i.e. Z tm i = (z i,0 , z i,1 , ..., z i,tm )) by evaluating the integral with g φ (·) as the integrand, where g φ (·) is parameterized by the neural networks and describes the latent dynamics. Finally, we utilized cause-specific decoder modules, which consists of fully connected neural networks, to estimate hazard function for each event over time.

Overview of the model architecture
The detailed model architecture of SurvLatent ODE is shown in Fig. 2. On the encoder side, we utilized an ODE-RNN (Rubanova et al., 2019) to learn the latent dynamics of the input temporal data and parameterize approximate posterior over z 0 (i.e. q(z 0 |X )). A function f γ (·), parameterized by neural networks with γ as trainable weights, specifies the dynamics of RNN hidden states, which allows the RNN to have continuously welldefined hidden states and incorporate inputs arriving at arbitrary times. This is the key element in the ODE-RNN encoder, which makes it unnecessary to explicitly impute missing measurements. When feeding the data into the model, we zero-fill missing measurements and concatenate the vector of indicators for missing measurements, m(·), as well as time elapsed since each feature's last observation, ∆(·) (Che et al., 2018), into the input data x(·) over the union of all time points in each batch. For a RNN model, we utilize a Gated Recurrent Unit (GRU) cell (Cho et al., 2014), which resolves the vanishing or exploding gradient problems of conventional RNNs and enables learning long-term temporal dependencies across the data. Given the set of longitudinal data X with the latest measurement time of t , the ODE-RNN encoder runs backwards in time from t to 0 to get the approximate posterior q(z 0 |X ) as follows: Figure 2: Detailed model architecture of SurvLatent ODE. The ODE-RNN encoder takes a patient-specific time-varying features X i along with missing pattern information defined over the union of discrete time points (0, 1, ..., t −1, t ) and parameterizes approximate posterior over the initial latent variable z i,0 . A black-box differential equation solver ODESolve(g φ , z i,0 , t = 0, ..., t m ) is then called over the sampled initial latent variable z i,0 to obtain the latent trajectory, Z tm i = (z i,0 , z i,1 , ..., z i,tm ). Finally, d β (·), which consists of cause-specific decoder modules, subsequent fully connected layer, and softmax layer, maps the corresponding latent trajectory to estimate cause-specific hazard function for each event.
where θ corresponds to a set of trainable weights in the GRU cell and u υ is a neural network with weights υ and relates the final hidden state of the ODE-RNN encoder to the mean and variance over z 0 . As a result of the ODE-RNN encoder, temporal latent dynamics of the hidden states are flexibly modeled by neural networks without resorting to any strong parametric assumptions on the dynamics (Rubanova et al., 2019).
On the decoder side, we utilize a multi-task learning framework (Caruana, 1997;Alaa and van der Schaar, 2017;Lee et al., 2019) and allow for the latent trajectory Z t to be shared across multiple events. For patient i, the latent trajectory Z t i is obtained by sampling the initial state z i,0 from the approximate posterior q(z i,0 ) and evaluating the integral, z i,0 + t 0 g φ (z i,u )du at pre-specified time points (i.e. t = 0, ..., t m , where t m is the end time point of the prediction window and may be arbitrarily far into the future). This integral can be estimated by calling a black-box differential equation solver, ODESolve(g φ , z i,0 , t = 0, ..., t m ) (Chen et al., 2018). Then, each cause-specific decoder module out of total b modules learns a flexible event-specific function which maps the latent trajectory Z t shared across b events to the cause-specific hazard function for each event λ k (t) (see Section 3.3.2). This operationalizes the intuition that multiple clinical outcomes can have partially shared latent mechanisms. For example, the clinical outcomes such as VTE, cancer-specific death, and metastasis can have a shared latent mechanism driven by types of malignant tumors and blood protein tests. Finally, the data decoder module reconstructs a patient-specific data trajectory across the pre-specified time window (i.e.x i (0),x i (1), ...,x i (t m )), which could potentially be used for interpolation and extrapolation tasks.

Survival function estimation
Survival function estimation in our model comes down to learning the function, d β (·), which maps the sampled patient-specific latent trajectory, Z tm i to the cause-specific hazard function for event k, λ * i,k (t) evaluated across the pre-specified time window [0, t m ]. Then, we utilize the hazard function to estimate event free survival probability as well as cause-specific cumulative incidence function over the appropriate time window. β corresponds to a set of trainable weights in cause-specific decoder modules and the subsequent fully-connected layer as shown in Fig. 2. Recall that in our discrete time-to-event analysis framework, where the hazard remains fixed in each time; bin widths are specified based on domain knowledge (e.g. a day or an hour). Note that our model does not specify any functional forms for how the data X relate to the λ * k (t) and allow the neural networks to flexibly parameterize hazard functions. With the assumption that data (X ) is independent of remaining time-to-event (T r ) conditioning on the sampled latent state (Z t ), the event free survival survival probability S(t|X i , Z t i ) for patient i is estimated using the corresponding latent trajectory Z t i (i.e. patient i's latent trajectory up until t) as follows:Ŝ where τ i,l is the latest measurement time for the patient. Setting a lower-bound for τ with τ i,l is necessary since S(t|Z t i ) should always be 1 for t ∈ [0, τ i,l ]. The proof for Equation 7 is shown in Appendix A. Notice that at any given time point and denotes probability of not experiencing any events at time τ given patient i survived up to the prior time point (i.e. τ − 1). We enforce this constraint using the softmax layer as shown in Fig. 2. Finally, the estimated cause-specific cumulative incidence functionF k (t|Z t i ) iŝ

Inference
We devise the loss function which handles time-varying features as well as loss to follow-up (i.e. right-censoring). The total loss consists of two components : the Kullback-Leibler (KL) divergence loss (i.e. D KL [q(z 0 |X ) || p(z 0 |X )]), where q(z 0 ) is the approximate posterior over the initial latent state and X is the set of time-varying features, and the log of the total survival likelihood. In practice, the KL loss cannot be directly evaluated since the posterior distribution p(z 0 |X ) is intractable because computing evidence p(X ) involves marginalizing out high-dimensional latent variable z 0 . Therefore, we instead maximize the evidence lower bound (ELBO) : where p(z 0 ) is assumed to be a standard normal, Φ corresponds to a set of trainable weights {θ, γ, υ, φ} in the ODE-RNN encoder and g φ (·) which specifies dynamics of the latent trajectory on the decoder side, and ζ is a set of trainable weights in the data decoder module, r ζ (·). Maximizing the ELBO in this way is equivalent to minimizing the KL divergence (Bishop and Nasrabadi, 2006). The total survival likelihood L surv , which enables handling of right-censored patients, is estimated as follows: where t r i is the observed remaining follow-up time from the latest measurement, k i is the observed event type of patient i, δ i is an event indicator, and β is a set of trainable weights in cause-specific decoder modules and the subsequent fully connected neural network (i.e. d β (·)). Note that the likelihood contribution from a right-censored patient (i.e. those with δ i = 0) is that they are alive at t r i , which corresponds to the event free survival probability S(t r i |Z i ) (see Equation 7). Therefore, the total loss we want to minimize is We perform an end-to-end training and utilize a standard back-propagation method to update weights in each module. For ODE solvers, we use the fifth order Dormand-Prince method (Dormand and Prince, 1980) from torchdiffeq Python package (Chen et al., 2018). Alternatively, the adjoint method, proposed in Chen et al. (2018), can be used to compute more memory-efficient gradients for ODE solvers, but at a cost of longer computation time.

Setup
To demonstrate clinical utility of our model, we evaluated SurvLatent ODE on two real-world observational datasets: MIMIC-III for predicting time to hospital mortality (i.e. singleoutcome) and the data from the Dana-Farber Cancer Institute (DFCI) for predicting time to VTE with all-cause mortality as a competing event. The datasets include routinely collected longitudinal measurements of patients' biomarkers where some measurements are often informatively missing. For both prediction tasks, patients were aligned at the first measurement time upon their admission. Finally, we randomly split each data into a training set (55%), validation set (15%), and test set (30%), tuned hyperparameters using train and validation sets (see Section 4.5), and obtained model performances on the held-out test set across 25th, 50th, and 75th percentiles of corresponding event times.

Datasets
MIMIC-III MIMIC-III is a large, publicly available longitudinal dataset of patients who were admitted to critical care units at the Beth Israel Deaconess Medical Center in Boston, Massachusetts (Johnson et al., 2016). In critical care scenarios, it is important to accurately identify high-risk patients as well as patient-specific time to mortality based on their timevarying biomarkers. This allows clinicians to determine the right level of care for patients in a timely manner based on their estimated risks and improve health outcomes while reducing unnecessary resource utilization.
For predicting time to hospital mortality, we utilized the cohort which consists of patients who were at risk after the first 36 hours of the admission (n = 21,728) and estimated their remaining time to mortality based on a total of 42 features. The features include static features such as age and sex as well as time-varying features such as vital signs and laboratory tests across the first 36 hours at a time resolution of 1 hour. For the longitudinal features, we included the top 40 most frequent vital signs and laboratory test results (see Appendix B for the full set of features). As an outcome, we used the time to hospital mortality measured from the 36th hour while right-censoring those who did not experience the mortality at their hospital discharge time. Among the total of 21,728 patients in our cohort, 2,281 patients (10.5 %) died at the hospital.
Dana-Farber Cancer Institute (DFCI) dataset Venous Thromboembolism (VTE) is a frequent, yet fatal complication in patients with active cancer, especially while they are receiving chemotherapy. VTE is associated with worse clinical outcomes such as elevated risk of mortality and reduced quality of life, and may lead to a significant disturbance of cancer treatment regimens (BLOM et al., 2006;Khorana et al., 2007;Lloyd et al., 2018). For ambulatory cancer patients, preventive anti-thrombotic measures such as thromboprophylaxis, a mechanical method to promote venous outflow, and anticoagulants drugs are effective at reducing the incidence of VTE (Rutjes et al., 2020;Key et al., 2020;Xiong, 2021), however with an increased risk of bleeding (Key et al., 2020). Therefore, accurate stratification of the VTE risk among patients with cancer may allow clinicians to improve clinical outcome while minimizing side effects due to overtreatment.
The dataset includes patients with active cancer who were admitted to the DFCI from 2/12/2013 to 12/14/2021 and received targeted panel sequencing of their tumor biopsy. We set the sequence date for each patient as the prediction time to avoid immortal time bias. For the experiment cohort, we included those who received the sequencing within a year of their admission and were at the risk of VTE without any recorded prior VTE event (n = 8,734). The total number of features is 64, which includes static features such as age, sex, presence of metastasis, and diagnosed cancer types as well as time-varying features such as laboratory test results and body mass index (BMI) across the first year from admission at a time resolution of one day (see Appendix B for the full set of features). We utilized ICD-10 diagnosis codes, I82.4 and I82.6, to determine the time of the venous thromboembolism (VTE) onset and I26 to determine the time of pulmonary embolism (PE) onset. And, we combined VTE and PE onset as the VTE event. We used the National Death Index (NDI) up until 12/31/2020 and the in-house death registry after the date to determine time of mortality. Finally, we right-censored patients at their last hospital visits, who did not experience VTE and death events. The demographics summary of the chosen cohort as well as causes of death (Table 5 and 6) are shown in Appendix D.
The main goal of this experiment was to develop a risk stratification framework for cancer-related VTE by accurately estimating the time to VTE event in the presence of a competing risk for death as well as missing measurements and loss to follow-ups. We thus estimated a patient specific remaining time to VTE from the sequence date with all-cause mortality as a competing event. Finally, to demonstrate clinical utility of the proposed model, we benchmarked against the categorical Khorana score (Khorana et al., 2008), which is the current clinical standard for cancer-related VTE stratification and utilizes four clinical measurements together with cancer type (see Section 4.3 for details).

Time to hospital mortality prediction (MIMIC-III)
Surv VAE-RNN, modified from (Che et al., 2018), adopts a variational autoencoder framework where both encoder and decoder are recurrent neural networks. This approach is closest to the proposed model except that in our model we utilized deep neural networks to learn the latent state dynamics.
Recurrent Deep Survival Machine (RDSM, Nagpal et al. (2021a)) is a fully parametric time-to-event model where it learns representation of the input longitudinal data using a RNN-based model and derives the conditional survival distribution (i.e. P (T r = t|X)) as a fixed mixture of parametric distributions (e.g. Weibull or Log-Normal).
Dynamic-Deephit (Lee et al., 2019) is a discrete time-to-event model which, like RDSM, learns input representation using a RNN-based model and estimates the conditional joint distribution of the survival time and event type (i.e. P (T = t, K = k|X)).
Cox Proportional Hazard model (Cox PH), proposed by Cox (1972), is a popular semi-parametric, time-to-event model which relates log of the relative hazard to a linear function of baseline features. For the experiment, we forward-filled the missing measurements to obtain the baseline features.
Refer to Table 1 for the summary of the baseline models as well as the proposed model in terms of key strategies in handling longitudinal data.  Cause-specific Cox model (CS Cox) is a Cox PH based model for handling competing risks and captures effects of the static baseline features on the cause-specific hazard function for each event (Putter et al., 2007). Similar to the Cox PH setting above, we forward-filled the missing values at the prediction time.
Fine-Gray Cox model (FG Cox), proposed by Fine and Gray (1999), is another Cox PH based model, widely used for competing risks. The model directly captures effects of the baseline features on the cumulative incidence function to make survival predictions. We adopted the missing measurement imputation strategy exactly identical to CS Cox.

Evaluation Metrics
Time-dependent cumulative/dynamic AUC corresponds to the probability that, given a pair of patients where one experienced the event of interest before t and the other is event-free at t, the model correctly ranks their risks of the event (Kamarudin et al., 2017).
Time-dependent Brier Score is a generalization of the conventional Brier score (Brier et al., 1950) which can account for the right-censored data and measures the mean squared error at a given time point t. This metric evaluates the model's calibration as well as discrimination performance (Xiu et al., 2020).
See Appendix C for details on the evaluation metrics.

Hyperparameters
We used the validation set for optimizing hyperparameters and utilizing early-stopping to avoid over-fitting. We utilized a random search (Bergstra and Bengio, 2012), with the set of hyperparameters including the number of dimensions in the encoder and decoder networks, the number of nodes in RNN hidden layers, the number of hidden layers and nodes in fully-connected neural networks for ODEs. For a fair comparison, we similarly tuned the deep learning baseline models including Surv RNN-VAE, RDSM, and Dynamic-Deephit. See Appendix B for the full list of hyperparameters and their possible values across all models.  (  state dynamics are modeled (i.e. parameterized exponential decay vs. neural networks), showed fairly poor performance across all metrics. This demonstrates that a Neural-ODE based approach allows the variational autoencoder framework to be successfully implemented for time-to-event modeling of longitudinal data, which would have been infeasible otherwise. SurvLatent ODE demonstrated competitive and robust performance over randomly subsampled test data with increasing missing rate (Fig. 3). To investigate each model's generalizability to more sparse patients' data, we evaluated the best model from the earlier experiment (Table 2) on 10 different randomly subsampled versions of the test set at each missing rate. At the baseline, each patient data in the test set already contains substantial missing measurements (∼ 70% on average) across the first 36 hours of admission and 40 time-varying features, driven by real-world mechanisms. RDSM, which relies on mean imputation for handling missing measurements, showed relatively poor performance across the induced missing rates. While Dynamic-Deephit showed consistent performance across the missing rates, SurvLatent ODE, with the substantial performance gain at the baseline real-world condition, continued to outperform all the other models through the artificially induced missing rate of 0.5.  Table 3: Performance summary of the models across the 25th, 50th, and 75th percentiles of event times on the held-out test set of the in-house DFCI data (n = 2,630) for predicting time to VTE event with an all-cause mortality as a competing event. Time-dependent AUC (higher the better) as well as Brier scores (lower the better) are used as evaluation metrics. Standard errors, shown in parenthesis, as well as statistical significance ( * : p < 0.05, * * : p < 0.001) were obtained exactly in the same manner as in the MIMIC-III experiment (see Table 2 caption). Note that Khorana score adopts the integer scoring system and its survival calibration performance cannot be evaluated.  Table 4: Performance summary of the models across the 25th, 50th, and 75th percentiles of event times for predicting time to all-cause mortality as the competing event. Performances for VTE (Table 3) and all-cause mortality (Table 4) were obtained from a single model evaluated on the same held-out test set of the in-house data (n = 2,630). Time-dependent AUC (higher the better) as well as Brier scores (lower the better) are used as evaluation metrics.

Results on VTE
As shown in Table 3 and 4, SurvLatent ODE significantly outperformed all the baseline models including Dynamic-Deephit in the discriminative as well as the calibration metrics, with a statistically significant improvement over the next-best method in two out of three time horizons for the discriminative performance. For predicting time to all-cause mortality (i.e. the competing event), SurvLatent ODE again outperformed all the baseline models in the discriminative metric, with statistically significant improvement for two out of three time horizons, while achieving competitive but weaker performance in the calibration metric. Specifically, Dynamic-Deephit showed superior calibration at the short survival quantile 2. As mentioned in Section 4.3, this cohort includes 1,031 patients with recorded relevant lab test results and BMI to compute Khorana scores at the prediction time (i.e. panel sequencing data).
(25%) while the Fine-Gray Cox model (FG Cox) showed superior calibration for the longer survival quantiles (50% and 75%). SurvLatent ODE significantly outperformed the current clinical standard, Khorana scores (Khorana et al., 2008) for predicting time to VTE events across all percentiles (Table 3). Khorana scores range from 0 to 6 and are meant to stratify risk for VTE among patients with cancer based on 5 different features including diagnosed cancer type and lab test results (see Section 4.3 for details). The performance of Khorana scores was comparatively poor, even after the imputation. Notably, SurvLatent ODE continued to outperform the Khorana score when focusing on the highest risk cancer types for VTE, pancreatic and stomach cancer (Fig. 5 in Appendix D). Furthermore, the proposed model demonstrated a strong performance across various subsets of patients defined by diagnosed cancer types and cancer stage (see Table 7 in Appendix D). The cumulative incidence trajectories were estimated by the Aalen-Johansen Estimator (Aalen and Johansen, 1978). Right : Heatmap visualization of feature-wise coefficient for each latent cluster. See Appendix D for detailed description of the method. SurvLatent ODE provided clinically meaningful and interpretable latent states. Features significantly associated with Cluster 4 include well-established risk factors for VTE. Note that * : p < 0.05, * * : p < 0.001 for statistical significance of feature coefficient.
As shown in Fig. 4, SurvLatent ODE provided clinically meaningful and interpretable latent states with the obtained latent clusters (k = 4) exhibiting significantly different VTE risk trajectories. The latent state for each patient was obtained from the output of the cause-specific decoder for the VTE event (see Appendix D for more detail). We obtained a cluster specific coefficient for each baseline feature by running a logistic regression with cluster membership as the binary outcome. Khorana scores were also included as a categorical feature to highlight the contribution of features beyond the current clinical guidelines. Well-established VTE risk factors were associated with the highest risk cluster (i.e. Cluster 4), including cancer type (e.g. pancreatic cancer) and known lab test results (platelet count : PLT, hemoglobin level : HGB, and white blood cell count : WBC) with consistent directions (Khorana et al., 2008). While some of these factors are included in the Khorana score, the fact that they were still significantly associated after conditioning on the Khorana score indicates that SurvLatent ODE identified additional latent trajectory variations among these factors. In addition to the features for the Khorana score, we observed significant associations with neutrophils (NEUT), albumins (ALB), and SGOT which were supported by the literature (Kapoor et al., 2018;Folsom et al., 2010Folsom et al., , 2014.

Discussion
We propose a full ODE-based variational autoencoder time-to-event framework for modeling longitudinal data, which encodes a patient-specific trajectory of time-varying features and decodes the resulting embedding into the latent trajectory for time-to-event prediction, where dynamics of the trajectories are learned by neural networks. Combined with the flexible estimation of cause-specific hazard functions for multiple events without explicit parametric assumptions on event time distributions, our proposed model outperforms the baselines including recently published deep-learning based models in two real-world datasets. Our proposed method offers a promising deep-learning based time-to-event framework for risk predictions in healthcare scenarios where data irregularities such as missing measurements and loss to follow-up are very common.
Our model significantly outperforms the current clinical standard, Khorana scores, for predicting the VTE event among patients with cancer. In addition, by learning complex nonlinear relationships between time-varying features, our model makes accurate time-to-event predictions for a wider range of patients while effectively capturing heterogeneity of the VTE risk, including identifying low-risk patients in conventionally high-risk cancer types. As future clinically-oriented work, we consider incorporating major fixed-effect variables, such as patient's tumor genetics, as well as more granular outcomes such as stroke, myocardial infarction, and cancer disease progression.
On the technical side, future work may involve estimating well-calibrated uncertainty of time-to-event predictions, reflective of training data availability for a given task as well as sample-wise prediction difficulty. In addition, as demonstrated in Garnelo et al. (2018);Norcliffe et al. (2021), one can explicitly capture uncertainty over underlying dynamics and render time-to-event estimates aware of the uncertainty in the dynamics across latent states. From the time-to-event modeling perspective, it is of interest to additionally adopt the subdistribution hazard framework which directly models effects of covariates on cause-specific cumulative incidence functions and is often used for an individual risk prediction (Austin et al., 2016;Berger et al., 2020).

Appendix A. Discrete time-to-event analysis
In Section 3.3.2, given that T r |= X |Z, we claim the following Proof For a discrete event time T r , we notice that where S(τ |Z t i ) = 1 for τ ≤ τ i,l . And, we have Summing up the hazard function across b events, we have Finally, rearranging Equation 14 and plugging it into Equation 12, we have   Table 6: Summary table for causes of death among those who died (n = 4,070) in the in-house data. We obtained the death causes through National Death Index (NDI) data. We considered all-cause mortality as a competing event for time to VTE event prediction.

D.2. Additional performance results
SurvLatent ODE outperformed the Khorana score for the cohort of patients with pancreatic and stomach cancer, where VTE risks are particularly high. SurvLatent ODE identified a low risk subpopulation that had event rates comparable to the population baseline (shown as a dotted line in Fig. 5), and substantially lower than the corresponding lowest Khorana score. In this set-up, to stratify risk groups at the prediction time, we computed the cause-specific restricted mean failure time (RMFT) from Day 0 (prediction date) to Day 365 (i.e.

365
0F k (t|X i )dt), which is interpreted as the expected number of days lost due to cause k before Day 365 (Mozumder et al., 2021). cancers (n = 239) among the held-out test set. By definition, these patients have at least a Khorana score of 2 due to their high-risk cancer types (Khorana et al., 2008). Cumulative incidence trajectories were estimated by Aalen-Johansen estimator (Aalen and Johansen, 1978), a non-parametric estimator for cumulative incidence, taking into account competing risks. Right : the same cohort of patients was divided into two groups (i.e. above median and below median) based on their cause-specific restricted mean failure time (RMFT) (i.e.

365
0F k (t|X i )dt) estimated by SurvLatent ODE. SurvLatent ODE based score was able to capture heterogeneity within the elevated risk group for VTE.  Table 7: Performance summary of SurvLatent ODE and Khorana scores on various groups of patients defined by diagnosed cancer types and cancer stage for predicting VTE risk. As a performance metric, we utilized the mean of time-dependent AUC from 25th percentile to 75th percentile of VTE event times. Although model performances vary widely across the groups, SurvLatent ODE outperforms imputed Khorana scores for all the groups except GI group. Standard errors, shown in parenthesis, as well as statistical significance ( * : p < 0.05, * * : p < 0.001) were obtained exactly in the same manner as in the MIMIC-III experiment (see Table 2 caption).

D.3. Latent state explanation
To obtain latent clusters, we utilized the temporal latent states provided from the causespecific decoder for VTE just before the fully-connected layer (see Fig. 2) at the inference time for the held-out test cohort (n = 2,630). Then, we summed up the latent state of each patient across one year from their latest measurement (i.e. prediction time) and obtained the latent state matrix with dimension of R Ntest×L , where L is the dimension of the latent trajectory. We then ran k-means clustering on the latent state using sklearn Python package (Pedregosa et al., 2011). Finally, we utilized the Aalen-Johansen estimator (Aalen and Johansen, 1978) to non-parametrically estimate the survival trajectory of patients in Figure 6: Cumulative Incidence trajectories, estimated by Aalen-Johansen Estimator for the all-cause mortality event across the clusters. The clusters were obtained, based on latent states output from the cause-specific decoder for all-cause mortality. each cluster. We used VTE-based latent clusters to obtain cumulative incidence trajectories for Fig. 4 and all-cause mortality based latent clusters (i.e. clusters based on the latent states provided from the cause-specific decoder for all-cause mortality) to obtain cumulative incidence trajectories for Fig. 6. To interpret each cluster, we first obtained the baseline data for the test cohort at the prediction time using a forward-fill imputation method. Then, we used the standardized baseline data to run logistic regression (Seabold and Perktold, 2010) on each feature conditioning on Khorana scores and see how significantly the feature is associated with each cluster encoded as the binary label (i.e. one-versus-rest). Finally, upon investigating the significance of coefficients for each feature, we dropped a set of features without nominal significance (p < 0.05) for any of the clusters. The result is shown in the main text (Fig. 4).