tdCoxSNN: Time-Dependent Cox Survival Neural Network for Continuous-time
  Dynamic Prediction
The aim of dynamic prediction is to provide individualized risk predictions
over time, which are updated as new data become available. In pursuit of
constructing a dynamic prediction model for a progressive eye disorder,
age-related macular degeneration (AMD), we propose a time-dependent Cox
survival neural network (tdCoxSNN) to predict its progression using
longitudinal fundus images. tdCoxSNN builds upon the time-dependent Cox model
by utilizing a neural network to capture the non-linear effect of
time-dependent covariates on the survival outcome. Moreover, by concurrently
integrating a convolutional neural network (CNN) with the survival network,
tdCoxSNN can directly take longitudinal images as input. We evaluate and
compare our proposed method with joint modeling and landmarking approaches
through extensive simulations. We applied the proposed approach to two real
datasets. One is a large AMD study, the Age-Related Eye Disease Study (AREDS),
in which more than 50,000 fundus images were captured over a period of 12 years
for more than 4,000 participants. Another is a public dataset of the primary
biliary cirrhosis (PBC) disease, where multiple lab tests were longitudinally
collected to predict the time-to-liver transplant. Our approach demonstrates
commendable predictive performance in both simulation studies and the analysis
of the two real datasets.
Lang Zeng, Jipeng Zhang, Wei Chen, Ying Ding
# Dynamic Prediction using Time-Dependent Cox Survival Neural Network

###### Abstract

The target of dynamic prediction is to provide individualized risk predictions over time which can be updated as new data become available. Motivated by establishing a dynamic prediction model for the progressive eye disease, age-related macular degeneration (AMD), we proposed a time-dependent Cox model-based survival neural network (tdCoxSNN) to predict its progression on a continuous time scale using longitudinal fundus images. tdCoxSNN extends the time-dependent Cox model by utilizing a neural network to model the non-linear effect of the time-dependent covariates on the survival outcome. Additionally, by incorporating the convolutional neural network (CNN), tdCoxSNN can take the longitudinal raw images as input. We evaluate and compare our proposed method with joint modeling and landmarking approaches through comprehensive simulations using two time-dependent accuracy metrics, the Brier Score and dynamic AUC. We applied the proposed approach to two real datasets. One is a large AMD study, the Age-Related Eye Disease Study (AREDS), in which more than 50,000 fundus images were captured over a period of 12 years for more than 4,000 participants. Another is a public dataset of the primary biliary cirrhosis (PBC) disease, in which multiple lab tests were longitudinally collected to predict the time-to-liver transplant. Our approach achieves satisfactory prediction performance in both simulation studies and the two real data analyses. tdCoxSNN was implemented in PyTorch, Tensorflow, and R-Tensorflow.

Cox model; dynamic prediction; neural network; survival analysis; time-dependent covariate +
Footnote †: This paper has been submitted for consideration in _XXX_
## 1 Introduction

For many chronic progressive diseases, the prognosis and severity of the disease change over time. A dynamic prediction model that can forecast the longitudinal disease progression profile is a crucial and unmet need (Jenkins et al. 2018). The unstructured observation times and the varying number of observations across subjects make it challenging to build a dynamic prediction model. The collection of high-dimensional longitudinal data requires the development of novel dynamic prediction models which can handle various inputs such as images.

Joint modeling and landmarking are the two dominating techniques for dynamic prediction. The first prediction approach jointly models the longitudinal and time-to-event data through a longitudinal sub-model and a survival sub-model (Rizopoulos 2011). However, joint modeling is computationally demanding (Rizopoulos et al. 2017) and struggles to directly model the large-scale longitudinal data. In contrast, landmarking is a more pragmatic model which avoids directly modeling the process for the longitudinal covariates. It estimates the effect of predictors through the survival model over all subjects at risk at a given landmark time point (Van Houwelingen 2007). Suresh et al. (2017) and Rizopoulos et al. (2017) compared the prediction accuracy of two models and found that joint modeling performs better than landmarking when the longitudinal process is correctly modeled. In the cases that the longitudinal process is misspecified or difficult to estimate, such as with sparse longitudinal data, the landmarking method provided a good enough prediction (Ferrer et al. 2019). Recently, there have been extensions to joint modeling methods aimed at addressing the nonlinear patterns in longitudinal outcomes. Li et al. (2022) proposed the functional JM model to model the multiple longitudinal outcomes as multivariate sparse functional. Zou et al. (2023) applied the functional JM model to predict the progression of Alzheimer's
disease using pre-specified MRI voxels. However, this approach heavily relies on image registration/pre-processing and disregards the correlation between voxels.

With the development of machine learning and its success in survival analysis (Katzman et al., 2018; Lee et al., 2018; Kvamme et al., 2019), new methods have been developed to integrate the dynamic prediction models with machine learning techniques to expand their application and enhance the prediction accuracy in more complex settings. Lin and Luo (2022) proposed using a neural network to jointly model the survival and longitudinal process. Tanner et al. (2021) combined the landmarking with machine learning ensembles to integrate prediction from standard methods. However, these two approaches can not be directly applied to the situation with high-dimensional longitudinal variables. Lee et al. (2019); Jarrett et al. (2019); Nagpal et al. (2021) proposed different deep-learning models for dynamic prediction under the discrete-time scenario. Although discretizing the time does not necessarily diminish prediction accuracy, the number of time intervals used for discretization significantly impacts accuracy so needs to be carefully tuned in practice (Sloma et al., 2021). In summary, no existing dynamic prediction models could directly handle the high-dimensional, longitudinal variables collected at irregular observational times.

The time-dependent Cox model (Fisher et al., 1999; Thomas and Reyes, 2014) is a straightforward continuous-time method used to incorporate the relationship between the longitudinal and time-to-event processes. This approach has received numerous criticisms as it may not accurately reflect the longitudinal process (Sweeting and Thompson, 2011). Nonetheless, we discovered that it can be easily combined with neural network techniques to create a dynamic prediction method capable of handling complex longitudinal markers (e.g. images) and their nonlinear relationship with survival outcomes. This paper proposes a dynamic prediction method on a continuous time scale when the longitudinal covariates could be high-dimensional and measured at unstructured time points. Specifically, we combined the
time-dependent Cox model with the Cox survival neural network. The proposed method can incorporate structured inputs (e.g. images, texts) through an additional neural network architecture. The rest of the article is organized as follows: Section 2 describes the AREDS dataset which highly motivates this study. Section 3 introduces notation and the standard dynamic prediction techniques. Section 4 presents the proposed model. Section 5 introduces the accuracy metrics for the evaluation of prediction performance. Section 6 and 7 presents the simulation and two real data analysis results. Finally, we conclude with a discussion in Section 8.

## 2 AMD progression prediction and existing works

This research is highly motivated by developing a dynamic prediction model for progressive eye disease, Age-related Macular Degeneration (AMD), using longitudinal fundus images. AMD is a polygenic and progressive neurodegenerative eye disease, which is a leading cause of blindness in the older population, especially in developed countries. It has been reported that by 2040, AMD is going to affect about 288 million people worldwide Peng et al. (2020). Once the disease is in the late stage, it is typically not curable. Therefore, accurate predictive models for predicting the risk of progressing to late-AMD at an early stage are needed. This will allow clinicians to identify high-risk individuals for late-AMD at their subclinical stage so that they can initiate preventive interventions for those individuals.

Colored fundus photographs have been routinely used to examine and document the presence and severity of AMD in clinical practice and trials. Our motivating study is the Age-related Eye Disease Study (AREDS), which is a large multi-center, controlled, randomized clinical trial of AMD and age-related cataract (Group 1999). It was designed to assess the clinical course and risk factors related to the development and progression of AMD and cataract. The participants were followed up for 12 years and the fundus photographs were performed every six months during the first six years, and annually thereafter.
Many models have been proposed in recent years to characterize and predict AMD progression. Sun et al. (2020) built a survival neural network prediction model to predict the progression risk with the baseline demographic and genotype data. Yan et al. (2020) used both genotype and fundus image data to predict the risk of progression to late-AMD at discretized given time points, where images were processed through the convolutional neural network (CNN). Peng et al. (2020) used the deep features of baseline fundus images obtained from DeepSeeNet (Peng et al., 2019) along with demographic and genotype data to generate individualized progression curves. Ghahramani et al. (2021) used the deep features of the fundus images from baseline, year 2, and 3 through a recurrent neural network (RNN). Ganjdanesh et al. (2022)trained a image generation model on all consecutive time-points data to predict the fundus image at next visit. These prediction models either rely solely on baseline predictors or predictors from specific years.

## 3 Notation and existing approaches for dynamic prediction

### Notations

Let \(\{T_{i},\delta_{i},\{\mathcal{Y}_{i}(t),0\leqslant t\leqslant T_{i}\},X_{i};i =1,\ldots,n\}\) denote \(n\) samples. \(T_{i}=\min(T_{i}^{*},C_{i}^{*})\) denotes the observed event time for subject \(i\), with \(T_{i}^{*}\) and \(C_{i}^{*}\) denoting the underlying true event time and censoring time, and \(\delta_{i}=I(T_{i}^{*}\leqslant C_{i}^{*})\) is the event indicator. \(X_{i}\) is time-inarying measurement for subject \(i\). \(\{\mathcal{Y}_{i}(t),0\leqslant t\leqslant s\}\) is the measurements in the interval \([0,s]\). It's often the case that the history can not be fully measured and we only observed \(\{\mathcal{Y}_{i}(0),\mathcal{Y}_{i}(t_{i1}),\ldots,\mathcal{Y}_{i}(t_{in_{i} }),t_{in_{i}}\leqslant s\}\). We focus on the scenario where \(t\) is continuous and allow the number of longitudinal observations \(n_{i}+1\) to differ across individuals.

For dynamic prediction, we are interested in predicting the probability that a new patient \(j\), with time-varying measurements up to time \(s\), will survive up to time \(u\) for \(u>s\), denoted as \(\pi_{j}(u|s)=Pr(T_{j}^{*}>u|T_{j}^{*}>s,\{\mathcal{Y}j(t),0\leqslant t \leqslant s\},X_{j})\). In contrast to static predictions,
dynamic models allow predictions to be updated and obtain \(\hat{\pi}_{j}(u|s^{\prime})\) when new information is available at a new time \(s^{\prime}>s\).

### Existing approaches

Joint modeling is a popular approach for dynamic prediction. It consists of a sub-model for the longitudinal process and a sub-model for the survival process linked through the shared random effects (Tsiatis and Davidian 2004) or some functional forms of the longitudinal outcomes Li and Luo (2019); Mauff et al. (2020). The estimation of the joint modeling is performed by the Bayesian approach (Rizopoulos 2011). It is usually computationally expensive and especially challenging with high-dimensional longitudinal features such as images. After estimating the parameters, the prediction of the survival probability \(\hat{\pi}_{j}(u|s)\) is obtained from the posterior predictive distribution of the survival process (Rizopoulos 2014).

Another popular dynamic prediction approach is the landmarking method (LM). Different from joint modeling, LM does not model the longitudinal process. It predicts \(\pi_{j}(u|s)\) through a model fitted based on subjects still at risk at time \(s\), which is called the landmark time. And the administrative censoring at \(u\) will be applied, hence the estimated effect (\(\hat{\beta}_{LM(s,u)}\)) of the predictors at landmark time can approximate the effect of the time-varying covariates on survival outcome within the time window \((s,u]\) (Van Houwelingen 2007). The predictors at landmark time \(s\) is usually a summary of the history. For example, one may use the mean or maximum of \(\{\mathcal{Y}_{j}(t),0\leqslant t\leqslant s\}\) as a summary. The choice of the time-dependent covariate depends on the research problem and is discussed in Fisher et al. (1999). The idea of generating a summarized time-dependent covariate can also be applied to JM. For simplicity, we used the instantaneous measurements \(\mathcal{Y}_{i}(s)\) as time-varying covariates in this work.
### Using time-dependent Cox model for dynamic prediction

Before we dive into the technical details, we want to highlight the difference between the time-dependent Cox method and the standard LM. The standard LM directly models the bilateral relationship between predictor variables at landmarking time \(s\) and the corresponding residual survival time with administrative censoring time \(u\) through a survival model (e.g. time-independent Cox model), with all observations after \(s\) disregarded. Both \(s\) and \(u\) need to be prespecified which makes the model less flexible. Moreover, using \(\hat{\beta}_{LM(s,u)}\) to approximate the true effect of time-dependent covariates can be inaccurate when \(u\) is large. With \(u\) more away from \(s\), the effect estimated in LM attenuates (Van Houwelingen, 2007; Putter and van Houwelingen, 2017) compared to the estimation from the time-dependent Cox model which considers the observations after landmarking time. For those reasons, we will not consider LM + the time-independent models such as (Van Houwelingen, 2007) and Pickett et al. (2021) in our simulations and analyses. Instead, following the spirit of LM, we fitted the time-dependent Cox model on the representative subjects with a given landmark time to evaluate the performance of the LM + time-dependent Cox model when we do have a landmark time of interest in simulation 1 (section 6.2.1).

Using the time-dependent Cox model is a trade-off between JM and the standard LM. Similar to JM, it fully uses the available longitudinal variables and is flexible (without pre-specifying landmark time \(s\) and administrative censoring time \(u\)), while keeping the simplicity of LM (without modeling the longitudinal process). We propose to construct a survival neural network under the time-dependent Cox model to incorporate the non-linear effect of the time-varying predictors, and can further use high-dimensional longitudinal predictors (e.g. image) as input.
## 4 Dynamic prediction using time-dependent Cox survival neural network

### Cox model with time-dependent covariates

The time-dependent Cox model takes the form

\[h_{i}(t)=h_{0}(t)\exp[g_{\beta,\gamma}(X_{i},\mathcal{Y}_{i}(t))]=h_{0}(t)\exp[ \beta^{T}X_{i}+\gamma^{T}\mathcal{Y}_{i}(t)]\]

and estimate \((\beta,\gamma)\) through the partial likelihood \(pL=\prod_{i}^{n}[\frac{\exp[g_{\beta,\gamma}(X_{i},\mathcal{Y}_{i}(T_{i}))]}{ \sum_{j:T_{j}\geqslant T_{i}}\exp[g_{\beta,\gamma}(X_{j},\mathcal{Y}_{j}(T_{i}) )]-E_{i}(g_{\beta,\gamma})}]^{\delta_{i}}\). With Efron's approximation for handling the tied events (Efron, 1977), we have \(E_{i}(g)=\frac{\sum_{j}1(j>i,T_{j}=T_{i})}{\sum_{j}1(T_{j}=T_{i})}\sum_{j:T_{j} =T_{i}}\exp[g(X_{j},\mathcal{Y}_{j}(T_{i}))]\).

In the partial likelihood for the time-dependent Cox model, at each event time \(T_{i}\), \(\mathcal{Y}(T_{i})\) is required for all the subjects still at risk at \(T_{i}\), which is not always available. Therefore, interpolation between longitudinal measurements is required. The last observation carried forward (LOCF) method is commonly used (Thomas and Reyes, 2014; Therneau et al., 2017), which assumes the values of the longitudinal predictors stay constant until the next measurement is available.

To predict \(\pi_{j}(u|s)\), similarly, we assume \(\mathcal{Y}_{j}(t)=\mathcal{Y}_{j}(s)\) for all \(s<t\leqslant u\). Therefore, the predicted survival probability is

\[\hat{\pi}_{j}(u|s)=\exp\{-(\hat{H}_{0}(u)-\hat{H}_{0}(s))\exp[g_{\hat{\beta}, \hat{\gamma}}(X_{j},\mathcal{Y}_{j}(s))]\}, \tag{4.1}\]

where \(\hat{H}_{0}(t)=\sum_{i=1}^{n}\frac{I(T_{i}\leqslant t)\delta_{i}}{\sum_{j:T_{ j}\geqslant T_{i}}\exp[g_{\beta,\hat{\gamma}}(X_{j},\mathcal{Y}_{j}(T_{i}))]-E_{i}(g_{ \beta,\hat{\gamma}})}\) is the Breslow estimator of the cumulative baseline hazard function (Breslow, 1972; Lin, 2007) with Efron's approximation.

### Cox survival neural network with time-dependent covariates

To establish a dynamic prediction model using high-dimensional image data, we consider the use of the neural network to augment the time-dependent Cox model. A neural network is an architecture that models the relationship between the input \(\mathbf{x}\in\mathbb{R}^{p_{0}}\) and output \(f(\mathbf{x})\in\mathbb{R}^{p_{L+1}}\) through the recursive layer structures \(f(\mathbf{x})=W_{L}\times\sigma_{\mathbf{V}_{L}}\left(W_{L-1}\times\sigma_{ \mathbf{V}_{L-1}}(\ldots W_{1}\times\sigma_{\mathbf{V}_{1}}(W_{0}\times \mathbf{x}))\right).\)\(L\) is the total number of hidden layers (depth of neural network) with \(p_{k}\) nodes (width) in each
layer \(k\) (\(k=1,\ldots,L\)). The activation function \(\sigma\) with the shift vector \(\mathbf{V}_{k}\) is a nonlinear transformation that operates componentwise \(\sigma_{\mathbf{V}_{i}}((y_{1},\ldots,y_{p_{k}})^{T})=(\sigma(y_{1}-v_{1}), \ldots,\sigma(y_{p_{k}}-v_{p_{k}}))^{T}\). The depth \(L\), width \(\mathbf{p}=(p_{1},\ldots,p_{L})\), and activation function \(\sigma\) should be pre-specified before the model fitting. The weight matrix \(W_{k}\in p_{k+1}\times p_{k}\) and the shift vector \(\mathbf{V}_{k}\in\mathbb{R}^{p_{k}}\) are the parameters that will be estimated by minimizing the loss function.

We apply the neural network to model the nonlinear effect of time-dependent covariates and incorporate the high-dimensional longitudinal data. Instead of assuming the unknown risk score function \(g_{0}(X_{i},\mathcal{Y}_{i}(t))\) takes the linear form \(g_{\beta,\gamma}(X_{i},\mathcal{Y}_{i}(t))=\beta^{T}X_{j}+\gamma^{T}\mathcal{ Y}_{j}(t)\), we leave the form of \(g_{0}(X_{i},\mathcal{Y}_{i}(t))\) unspecified and model it through a neural network \(g_{\theta}(X_{i},\mathcal{Y}_{i}(t))\) parameterized by the weight matrixes and shift vectors of the neural network \(\theta=(\mathbf{W},\mathbf{V})\).

The time-dependent Cox survival neural network (Cox SNN) is a feed-forward neural network that models the effect of time-dependent covariates on their hazard function. The input is the covariates \((X_{i},\mathcal{Y}_{i}(t))\) and the output \(g_{\theta}(X_{i},\mathcal{Y}_{i}(t))\) is a single node with linear activation function so that \(g_{\theta}(X_{i},\mathcal{Y}_{i}(t))\in\mathbb{R}\) (Figure 1a). \(\theta\) is estimated by minimizing the negative log partial likelihood \(\hat{\theta}=\arg\min[l(\theta|g_{\theta})]\), where

\[l(\theta|g_{\theta})=-\frac{1}{n}\sum_{i=1}^{n}\delta_{i}\left[g_{\theta}(X_{ i},\mathcal{Y}_{i}(T_{i}))-\log\left(\sum_{j:T_{j}\geqslant T_{i}}\exp\{g_{ \theta}(X_{j},\mathcal{Y}_{j}(T_{i}))\}-E_{i}(g_{\theta})\right)\right]. \tag{4.2}\]

Following the prediction formula (4.1) for the time-dependent Cox model, the predicted probability under the Cox SNN is given by

\[\hat{\pi}_{j}(u|s)=\exp\{-(\hat{H}_{0}(u)-\hat{H}_{0}(s))\exp[g_{\tilde{\theta }}(X_{j},\mathcal{Y}_{j}(s))]\}, \tag{4.3}\]

where \(\hat{H}_{0}(t)=\sum_{i=1}^{n}\frac{I(T_{i}\leqslant t)\delta_{i}}{\sum_{j:T_{ j}\geqslant T_{i}}\exp[g_{\tilde{\theta}}(X_{j},\mathcal{Y}_{j}(T_{i}))]-E_{i}(g_{ \tilde{\theta}})}\). \(\hat{\pi}_{j}\) can be updated when new information of subject \(j\) is available (Figure 1b).

The estimator of \(\theta\) which can minimize the loss (4.2) is not unique. For a given minimizer \(\hat{\theta}\), it is able to find \(\tilde{\theta}\) such that \(g_{\tilde{\theta}}:=g_{\tilde{\theta}}+c\). Note that \(\tilde{\theta}\) is also a minimizer since \(l(\tilde{\theta}|g_{\theta})=l(\hat{\theta}|g_{\theta})\). However, the constant shift \(c\) of the estimated hazard function will not change the predicted
probability in (4.3) (as it appears both in the numerator and the denominator). Therefore, \(\hat{\pi}_{j}(u|s)\) in (4.3) is robust to a constant shift of the neural network output.

Besides the ability to model the nonlinear effects, another benefit of using the neural network structure is that a pre-trained neural network, which can process a specific data type such as image or text, can be readily combined with SNN (transfer learning). For example, ResNet50 (He et al., 2016) is a CNN for image classification whose weights have been trained over one million images. Adding the pre-trained CNN on the top of SNN allows SNN to take the raw images as input, which is the approach we take for this study.

## 5 Prospective accuracy metrics

Methods for assessing the predictive performance of survival models concentrates either on calibration, i.e., how well the model predicts the observed data, or on discrimination, i.e., how well the model discriminates between subjects with the event from subjects without. We consider both calibration and discrimination metrics for model performance evaluation under the time-dependent setting.

Similar to previous studies, we compare dynamic prediction models on a time-window \((t,t+\Delta t]\) where the landmark time \(t\) and the length of time window \(\Delta t\) are pre-specified (Rizopoulos, 2011; Rizopoulos et al., 2017; Tanner et al., 2021). This is to evaluate, for subjects survived to time \(t\) on a separate test data with time-dependent covariates \(\mathcal{Y}_{i}(t)\) up to time \(t\) collected, how well the predicted \(\hat{\pi}_{i}(s|t)\) (\(t<s<t+\Delta t\)) agrees with the observed data. The censoring-free probability \(G(t)=P(C>t)\) is used to account for the censoring as a weight in the metric.
### Calibration metric: time-dependent Brier Score

The time-dependent Brier Score measures the mean squared error between the observed survival status and the predicted survival probability weighted by the inverse probability of censoring (IPCW) (Gerds and Schumacher, 2006). A lower Brier Score indicates a higher prediction accuracy. For a given landmark time \(t\), the estimated Brier Score at time \(t+\Delta t\) is \(\hat{B}S(t,\Delta t;\hat{\pi})=\frac{1}{\sum_{i}1(T_{i}>t)}\sum_{i=1}^{n} \left(\mathbbm{1}(T_{i}>t)\hat{W}_{i}(t,\Delta t)\{1(T_{i}>t+\Delta t)-\hat{ \pi}_{i}(t+\Delta t|t)\}^{2}\right),\) where \(\hat{W}_{i}(t,\Delta t)=\{\frac{1(T_{i}>t+\Delta t)}{\hat{G}(t+\Delta t|t)}+ \frac{1(T_{i}\leqslant t+\Delta t)\delta_{i}}{\hat{G}(T_{i}^{-}|t)}\}\) is the IPCW weight and \(\hat{G}(s|t)=\frac{\hat{G}(s)}{\hat{G}(t)}\) is the Kaplan-Meier estimate of the conditional censoring distribution.

### Discrimination metric: time-dependent AUC

The area under the receiver operating characteristic curve (AUC) measures the discrimination of the prediction model. It ranges from 0 to 1 with 0.5 indicating that the discriminability is no better than random guessing and AUC further away from 0.5 suggesting better discrimination. We use cumulative sensitivity and dynamic specificity AUC (cdAUC) (Kamarudin et al., 2017) to evaluate the discrimination performance of the models at different time points. Given a threshold \(b\) and predictor \(X\), the cumulative sensitivity is defined as \(Se^{C}(b,\Delta t)=P(X_{i}>b|T_{i}\leqslant\Delta t)\) while the dynamic specificity is \(Sp^{D}(b,\Delta t)=P(X_{i}\leqslant b|T_{i}>\Delta t)\). The term "cumulative" is to differentiate this sensitivity from the incident sensitivity \(Se^{I}(b,\Delta t)=P(X_{i}>b|T_{i}=\Delta t)\) which assesses the sensitivity for the population whose survival time exactly equals \(\Delta t\). With the cumulative sensitivity and the dynamic specificity, the corresponding cdAUC\((t,\Delta t)=P(X_{i}>X_{j}|t<T_{i}\leqslant\Delta t,T_{j}>\Delta t)\), \(i\neq j\). Specifically, for a given time interval \((t,t+\Delta t]\), the IPCW estimator of cdAUC is

\[\text{cd}A\hat{U}C(t,\Delta t)=\frac{\sum_{i=1}^{n}\sum_{j=1}^{n} \mathbbm{1}_{(\hat{\pi}_{i}(t+\Delta t|t)<\hat{\pi}_{j}(t+\Delta t|t))}\delta _{i}\mathbbm{1}_{(t<T_{i}\leqslant t+\Delta t)}\mathbbm{1}_{(T_{j}>t+\Delta t )}W_{i}(t,\Delta t)W_{j}(t,\Delta t)}{\sum_{i=1}^{n}\sum_{j=1}^{n}\delta_{i} \mathbbm{1}_{(t<T_{i}\leqslant t+\Delta t)}\mathbbm{1}_{(T_{j}>t+\Delta t)}W_ {i}(t,\Delta t)W_{j}(t,\Delta t)}.\]

It computes the IPCW weighted percentage of the comparable subject pairs \((i,j)\) where their predicted survival probabilities are consistent with their observed data for the given
time interval \((t,t+\Delta t]\). The comparable pair \((i,j)\) is two subjects in which the subject \(i\) experiences the event within the time interval \((t,t+\Delta t]\) and subject \(j\) is event-free by \(t+\Delta t\).

## 6 Numerical implementation and simulations

### Numerical implementation

For the time-dependent Cox SNN, we implemented the log partial likelihood function with Efron tie approximation (4.2) using Tensorflow (Abadi et al., 2016). It is implemented through matrix operations so the calculation is fast (details can be found in the Appendix). The PyTorch (Paszke et al., 2019) and R-Tensorflow (Allaire and Tang, 2019) version can also be found at [https://github.com/langzeng/tdCoxSNN](https://github.com/langzeng/tdCoxSNN). The neural network was optimized through the Adam optimizer (Kingma and Ba, 2014). We used the following survival neural network structure in all simulations and real data analysis: input layer \(\rightarrow\) hidden layer \(\rightarrow\) batch normalization layer \(\rightarrow\) dropout layer \(\rightarrow\) output layer. The batch normalization layer (Ioffe and Szegedy, 2015) accelerates the neural network training and the dropout layer (Srivastava et al., 2014) protects the neural network model from over-fitting. Hyper-parameters were also fixed in all analyses: 30 nodes in the hidden layer, Scaled Exponential Linear Unit (SeLU) as the hidden layer activation function, batch size 50, epoch size 20, learning rate 0.01, and dropout rate 0.2.

### Simulation

#### 6.2.1 Simulation 1: Low dimensional predictors

We carried out simulation studies to empirically compare the dynamic prediction performance of the proposed time-dependent Cox SNN with the time-dependent Cox model and joint modeling in a low-dimensional setting. Data were generated through the joint models. For all simulations, for sample \(i\), one-dimensional longitudinal covariate \(\mathcal{Y}_{i}(t)\) was generated through \(\mathcal{Y}_{i}(t)=y_{i}(t)+\epsilon\) where \(y_{i}(t)\) is the true longitudinal trajectory over time and \(\epsilon\sim N(0,0.3^{2})\) represents a measurement
error. We considered the true value of time-varying covariate to be given by

\[y_{i}(t)=\beta_{0}+\beta_{1}t+\beta_{2}t^{2}+b_{i0}+b_{i1}t+b_{i2}t^{2}\text{ and }\mathbf{b}_{i}=(b_{i0},b_{i1},b_{i2})^{T}\sim N(0,\Sigma_{3\times 3})\]

where \(\beta_{0}=3.2,\beta_{1}=-0.07\). \(\Sigma\) denotes a 3 by 3 inter-subject variance matrix with \(\Sigma_{11}=1.44,\Sigma_{22}=0.6\) and we assume the covariances \(\Sigma_{ij}\) are zero in all simulations. \((\beta_{2},\Sigma_{33})\) capture the non-linearity of the trajectory and we considered the trajectory of the longitudinal measurement to be linear (\(\beta_{2}=0,\Sigma_{22}=0\)) or nonlinear (\(\beta_{2}=0.004,\Sigma_{22}=0.09\)). \(y_{i}(t)\) was measured regularly per time unit at \(t=0,1,2,\ldots,14\). The survival time \(T_{i}^{*}\) was obtained through a Weibull model \(h_{i}(t)=\lambda\rho t^{\rho-1}\exp\{g(X_{i},y_{i}(t))\}\) with \(\rho=1.4\) and \(\lambda=0.1\). The censoring time \(C_{i}^{*}\) was generated through an exponential distribution \(\exp(\frac{2}{14})\). The observed survival time \(T_{i}=\min(T_{i}^{*},C_{i}^{*})\) and the event indicator \(\delta_{i}=1\left(T_{i}^{*}\leqslant C_{i}^{*}\right)\) were then calculated. Longitudinal values \(\mathcal{Y}_{i}(t)\) with \(t\geqslant T_{i}\) were disregarded, only \(\mathcal{Y}_{i}(t)\) measured before \(T_{i}\) were kept as the observed longitudinal measurements for subject \(i\). The baseline covariates \(x_{k}\) (\(k=1,2,3,4\)) were independently generated from the continuous uniform distribution on \([-0.5,1.5]\). Models were fitted using variables \((x_{1},x_{2},x_{3},x_{4},\mathcal{Y}_{i}(t))\). We considered four different cases (see below) with the linear or nonlinear risk function \(g(X,y(t))\) and the linear or quadratic trend (in time) of longitudinal measurement \(y(t)\). We added the intercept term -10 to make the censoring rate close to that of AREDS data (around 80%) in each simulation.

\(\bullet\): Case 1: \(\begin{cases}g(X,y(t))=x_{1}+2x_{2}+3x_{3}+4x_{4}+0.3y(t)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 2: \(\begin{cases}g(X,y(t))=x_{1}+2x_{2}+3x_{3}+4x_{4}+0.3y(t)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t+(0.004+b_{2})t^{2}\end{cases}\)\(\bullet\): Case 3: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{\frac{1} {3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 4: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 5: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 6: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 7: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 8: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 9: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 10: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 11: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 12: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 13: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 14: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 15: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 16: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\end{cases}\)\(\bullet\): Case 17: \(\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{ \frac{1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(
* Case 4: \[\begin{cases}g(X,y(t))=(\{x_{1}^{2}x_{2}^{3}+\log(x_{3}+1)+(0.3y(t)x_{4}+1)^{\frac{ 1}{3}}+\exp(\frac{x_{4}}{2})+0.3y(t)\}^{2}/3)-10\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t+(0.004+b_{2})t^{2}\end{cases}\]

The time-dependent Cox model was fitted through the R function {survival::coxph} (Therneau and Lumley, 2015). For joint modeling, we used the R package {JMBayes} (Rizopoulos, 2014; Rizopoulos et al., 2021) and included the linear effect of time \(t\) in the longitudinal sub-model. The time-dependent Brier Score and time-dependent AUC were calculated through the R packages {pec} (Mogensen et al., 2012) and the R package {timeROC} (Blanche et al., 2013), respectively.

In each setting, we performed 100 simulation runs. Models were fitted on the training set and the prediction metrics were evaluated over the separate test samples. We compared \(n_{train}=500\) and 1000 to evaluate the sample size effect on the performance of SNN, since a large sample size is usually required to train a neural network well. Separate test samples with \(n_{test}=200\) were generated in each run to evaluate the fitted models.

We set landmark time \(t=1\) in all simulations. Comparisons were made across seven models including the landmark time-dependent Cox model (LM-CoxPH), the landmark time-dependent Cox SNN (LM-CoxSNN), the time-dependent Cox model (CoxPH), the time-dependent Cox SNN (CoxSNN), the Joint modeling (JM), the Kaplan-Meier (KM), and the true conditional survival curve (Truth). The LM-CoxPH and LM-CoxSNN were fitted among training samples still at risk at the landmark time. CoxPH and CoxSNN were fitted over all training subjects. For the longitudinal sub-model in JM, we included the main effect of time in the fixed-effects part and an intercept and a time term in the random-effect design matrix. Predictions using the KM and the true model represented the worst and the best prediction we could obtain. Calibration metric \(\hat{BS}(t,\Delta t)\) and discrimination metric \(A\hat{UC}(t,\Delta t)\) for the seven methods were calculated at \(\Delta t=1,2,3,4\) from the landmark time.

[Figure 2 about here.]
#### Simulation 2: High-dimensional predictors

We also evaluated the performance of time-dependent Cox SNN with high-dimensional predictors. The simulation mechanism is the same as simulation 1 in section 6.2.1. After data was generated, at each visit time, we mapped the true risk score \(g(X,y(t))\) to handwritten digit images from the MNIST (LeCun 1998). Specifically, we standardized \(g(X,y(t))\) to \([0,0.99]\) and rounded it to 2 decimal places. Then the two handwritten digit images representing the tenth and hundredth numbers were randomly sampled in the corresponding digit class of MNIST. The two \(28\times 28\times 1\) images are treated as the observed predictor at time \(t\) to train the model as well as make the prediction.

To deal with the images, the convolutional layers and max pooling layers were added on top of the SNN structure introduced in section 6. Details can be found in Appendix. The time-dependent Cox SNN (td-CoxSNN) was fitted directly using the longitudinal images. As a comparison, the baseline Cox SNN (Base-CoxSNN) was fitted on the baseline images only. The baseline Cox model (Oracle-BaseCoxPH) and time-dependent Cox model (Oracle-tdCoxPH) were fitted on the true \(g(X,y(t))\) to represent the best performance that CoxSNN can achieve with image predictors.

* Case 5: \(\begin{cases}g(X,y(t))=x_{1}+2x_{2}+3x_{3}+4x_{4}+0.3y(t)-5\\ y(t)=(3.2+b_{0})-(0.07+b_{1})t\\ g(X,y(t))=x_{1}+2x_{2}+3x_{3}+4x_{4}+0.3y(t)-5\\ y(t)=3.2+b_{0}\end{cases}\)

Two simulation cases were considered. The setting of Case 5 is to evaluate the performance of the proposed method with time-varying longitudinal high-dimensional predictors (images). Case 6 is a scenario without time-varying covariates but the longitudinal images vary since the mapping of risk scores to images was performed separately on each visit. Therefore, the
only difference between Base-CoxSNN and td-CoxSNN is that td-CoxSNN would see more images. We performed 100 simulation runs for each case \(n_{train}=2,000\) and \(10,000\). The prediction accuracy \(\hat{BS}(t,\Delta t)\) and \(\hat{AUC}(t,\Delta t)\) was evaluated on a separate \(n_{test}=200\) test samples at \(t=1\) and \(\Delta t=1,2,3,4\).

### Simulation results

#### 6.3.1 Simulation 1: Low dimensional predictors

Figure 2 and Figure 3 display the boxplots of cdaUC and BS at times 1,2,3,4 after the landmark time over 100 simulation runs for each case. The plots indicate that the model with better discrimination ability (the higher cdaUC) tends to have better calibration ability (the lower BS) as well. In general, the CoxSNN is strongly competitive at the different time points after the landmark time. Under the complex setting of case 3 and case 4, CoxSNN outperformed CoxPH, and LM-CoxSNN outperformed LM-CoxPH. In simpler case 1 and case 2 where the effects of predictors are linear, CoxSNN and LM-CoxSNN performed similarly to the CoxPH and LM-CoxPH. This demonstrates that the time-dependent Cox SNN is able to learn the nonlinear effect in complex settings while maintaining a similar performance as the time-dependent Cox model when the effect is linear. Besides, JM performed the best when both the effect of longitudinal predictors and the sub-model for the longitudinal process were correctly specified (case 1). In case 3, JM correctly modeled the longitudinal process but didn't outperform neural network models as it failed to reflect the nonlinear effect of the predictors. Moreover, JM performed the worst in cases 2 and 4, when the longitudinal sub-model was incorrectly specified. This implies the performance of joint modeling highly depends on the correctness of the longitudinal sub-model.

We further evaluate the choice of sample for the neural network fitting. In case 1 and case 2, the prediction accuracy of models fitted with all subjects (CoxPH and CoxSNN) is close to the accuracy of landmarking models fitted with those still
time (LM-CoxPH and LM-CoxSNN). In settings where the effect of predictors is complex (case 3 and case 4), the landmarking methods demonstrated slightly better performance, suggesting that landmarking can potentially enhance prediction accuracy, particularly when there is a landmark time of interest. The landmark individuals can serve as an appropriate representation of the target population who survived the landmark time.

Lastly, when \(n_{train}\) decreased from 1000 to 500, the performance of SNN is relatively stable. Note that when \(n_{train}=500\) with 80% censoring rate, there is only 100 observed event in the training data, suggesting that the SNN can perform well with moderate sample sizes.

Overall, our simulation shows that time-dependent Cox survival neural network can achieve satisfactory prediction performance across diverse settings.

#### 6.3.2 Simulation 2: High-dimensional predictors

Figure 4 presents the results of the simulation with longitudinal images as predictors in case 6. In this scenario, the risk score \(g(X,\mathcal{Y}(t))\) is constant over time while the longitudinal images mapped to the risk score are time-dependent (the information embedded in images is time-independent). For oracle-CoxPH models, since they used the true risk score as the predictor and it is time-independent, they had the same prediction accuracy as the prediction from the real survival curve. When using the images as predictors, td-CoxSNN outperformed baseline-CoxSNN, especially when the training sample size (\(n_{train}=2,000\)) is small. This demonstrates the priority of td-CoxSNN over baseline-CoxSNN even the longitudinal images are uninformative.

When the image embeddings are time-varying (case 5), the baseline models were less accurate than the time-dependent models, as the baseline models ignored the longitudinal measurements and hence failed to evaluate the effects of time-varying predictors. Additionally, the baseline-CoxSNN may also be suffered from seeing fewer images as demonstrated above. The boxplots of cdaUC and BS for case 5 are provided in the appendix.

[FIGURE:S6.F4][ENDFIG
## 7 Real data analysis

### Application to AMD progression prediction

We applied the proposed time-dependent SNN on AREDS data and built a dynamic prediction model using longitudinal fundus images. In AREDS data, at each follow-up for each eye, there is a fundus image along with multiple manually graded image features performed by a medical center, for example, the size of the abnormal area (drusen) of a fundus image. After removing the eyes with late-AMD at baseline and images with low quality (i.e., image features not gradable or missing), our working data included 53,076 eye-level observations from 7,865 eyes of 4,335 subjects. The median follow-up time is 9.9 years, and 20.5% of the eyes progressed to late-AMD by the end of the study (Table 1). For the prediction perspective, we excluded observations that were measured after the diagnosis of late-AMD and conducted the eye-level analysis without considering the correlation between the two eyes from the same individual.

We performed a 5-fold cross-validation as follows. Data were split into five folds where models were trained on each 4-folds (training set) and evaluated on the remaining 1-fold (test set). The split was done on the subject level to ensure the two eyes from the same subject were included in the same fold. To fit the time-dependent Cox model and the time-dependent Cox SNN, longitudinal data were formatted into multiple intervals [tstart, tstop) where each interval represents the time window between one visit and the next visit (Figure 1a). To evaluate the prediction accuracy, because the fundus images were taken at different time points across subjects, on the test set, we chose individualized landmark time \(t_{i}\) as follows. For each subject \(i\) in the test set, \(t_{i}\) was chosen randomly from their longitudinal measurement time points prior to the last observed time \(T_{i}\). This is to mimic the real-world scenario when a new subject \(i\) comes at time \(t_{i}\), we use their measurements up to \(t_{i}\) to predict their survival
at \(t_{i}+\Delta t\). Prediction metrics (BS and cdauc) were evaluated on the \(\Delta t=1,\ldots,7\) windows across all subjects (Figure 5a).

We compared three dynamic prediction approaches: joint modeling (JM), time-dependent Cox model (CoxPH), and time-dependent Cox SNN (CoxSNN). The predictors include baseline demographic variables (age at enrollment, educational level, and smoking status) and seven longitudinal manually graded image features, which were found significant from a multi-variable baseline Cox regression model (Table 1). We also fitted the time-dependent Cox SNN directly on the longitudinal fundus images where the images were modeled by a CNN. Specifically, on top of the survival neural network, we added DeepSeeNet, which is a well-trained CNN for grading late-AMD using fundus images (Peng et al., 2019, 2020). The 256 nodes from the DeepSeeNet hidden layer and the baseline demographic nodes together formed the input layer of time-dependent Cox SNN. The details of the survival neural network structure and hyper-parameters can be found in section 6.1.

Figure 5a shows that the Brier Scores are lowest for the time-dependent Cox SNN either using the seven manually graded image features or directly using the fundus images across all prediction time points. As for the cdauc, the time-dependent Cox SNN using the longitudinal images is comparable to the three dynamic prediction models fitted on the seven manually graded image features. In practice, grading the image features is labor intensive and requires special medical image expertise. This prediction model (time-dependent Cox SNN with fundus images) can directly handle the raw colored fundus images without any further input from clinicians. Overall, joint modeling performed the worst in terms of the two accuracy metrics except for the long-term prediction (\(\Delta t=7\)).

For the interpretation of the SNN model fitted with longitudinal images, we generated the saliency map to visualize the regions with the most significant impact on the risk score for
a given subject. Figure 5b displays a fundus image from a participant's left eye at year 2.3 (since being on the AREDS study) where the eye showed multiple large drusens (i.e., the yellow spots). Our model predicts this eye will develop late-AMD with a high probability (60.1% since \(\Delta t=2\) year and 70.5% since \(\Delta t=2.8\) year) and the truth is this eye developed the disease 2.7 years later. We can see the saliency map successfully detects the pathological areas which are predictive of the disease progression.

After completing the model training, it is also possible to identify individuals with a higher risk of developing the disease through the estimated risk score \(\hat{g}\), which is the output of tdCoxSNN. For illustrative purposes, we used the training and test data from the first cross-validation split. The baseline risk score \(\hat{g}\) for 1,582 eyes in test data was estimated using their baseline fundus images and demographics. We identified two subgroups from the gaussian mixture model and compared their survival curves (Figure 6). During the follow-up period, a majority of the individuals in high-risk groups developed the disease, whereas those in low-risk groups maintained better health. We further compared the seven baseline manually extracted image features between the two groups (Table 2). The baseline fundus images in the high-risk group exhibited a greater number of higher-risk image features (\(p<0.001\), for each feature), which are significantly associated with the late-AMD (Table 1). The difference between the two groups for subject-level characteristics (excluding those with two eyes of differing risk) was found to be small. This suggests that the identification of subgroups was primarily based on the baseline fundus image.
### Application to PBC2 data

We also applied the proposed method (continuous-time model) to a publicly-available dataset with low-dimensional longitudinal predictors and compared with three discrete-time deep learning models. This data was collected between 1974-1984 for the research of primary biliary cirrhosis (PBC) disease (Fleming and Harrington, 2011). The dataset consists of 312 subjects (1,912 visits) with the event of interest being time-to-liver transplant and a censoring rate of 55%. The predictors include 12 longitudinal variables (7 continuous lab tests, such as albumin, and 5 categorical, such as liverenlargement) and 3 baseline variables (e.g. gender, age at start-of-study, and treatment indicator).

We followed the data processing steps from Putter and van Houwelingen (2017) for a fair comparison. The LM-tdCoxSNN (fitted for each landmark time) and tdCoxSNN (using all subjects) were fitted on the discrete-time scale. The prediction accuracies for the three discrete-time deep learning models were obtained directly from (Putzel et al., 2021). Each entry represents an average of accuracy calculated at 2,4,6, and 8 months after the landmark times.

Table 3 presents the prediction accuracy of the five methods across different landmark times. There is no universally optimal deep learning method, and our method is comparable to the existing discrete-time methods when applied to the discrete-time data. We observed that tdCoxSNN outperformed LM-tdCoxSNN in terms of both BS and dynamic C-Index, suggesting that the proposed method could benefit from retaining more samples during training, especially given the relatively small total sample size (1,912 visits in total).

[Table 3 about here.]
## 8 Discussion

We combined the time-dependent Cox model with the survival neural network to establish a dynamic prediction model in a continuous time scale. The proposed approach not only provides a powerful tool to model the non-linear effect of predictors on the risk but also allows users to directly incorporate the longitudinal high-dimensional features such as images without modeling the longitudinal process. Due to the neural network nature of SNN, existent neural networks can be added on top of the time-dependent Cox SNN to take advantage of well-developed deep learning structures for complex data, for example, RNN for sequential data (Lee and Dernoncourt 2016) and CNN for images (He et al. 2016). With the availability of more and more high-dimensional biomarkers with non-linear effects and unstructured longitudinal trends, our approach makes it possible to build a dynamic prediction model using complex longitudinal biomarkers (e.g., MRI, metabolomics data) for future research.

A limitation of the proposed method is the Last Observation Carried Forward (LOCF) assumption between visits, which may not accurately reflect reality. One potential direction for future research involves integrating joint modeling with machine learning techniques, although this could be computationally demanding. Additionally, our method was only compared to a limited number of discrete-time machine learning dynamic prediction approaches using a single dataset (section 7.2). Further efforts are necessary to thoroughly benchmark these methods across various settings.

Compared with discrete time dynamic prediction models (Lee et al. 2019; Tanner et al. 2021; Lin and Luo 2022), our model does not require the selection of time intervals for discretization, which makes it easier to process the longitudinal data for model fitting. Our model achieves satisfactory prediction performance in both simulation and real data analysis. It is worth noting that the same SNN structure and hyper-parameters introduced in section 6.1 worked well through all analyses. The training procedure with 20 epochs makes
fitting the SNN very fast. Additional tuning of the hyper-parameters may further improve the prediction accuracy of the model.

In this work, we used the saliency map to help interpret the fitted SNN model. One alternative solution for model interpretation is to use a partially linear Cox model where the risk score consists of a parametric component for predictors of interpretation interest and a nonparametric component modeled through the neural network. Zhong et al. (2022) proved the semiparametric efficiency of the parametric estimator which allows the method to make inferences along with using the SNN model for nuisance covariates. For future research, one may consider the partially linear structure in the time-dependent Cox SNN to improve model interpretability while maintaining the flexibility for modeling complex non-linear effects.
