Case-Base Neural Networks: survival analysis with time-varying,
  higher-order interactions
In the context of survival analysis, data-driven neural network-based methods
have been developed to model complex covariate effects. While these methods may
provide better predictive performance than regression-based approaches, not all
can model time-varying interactions and complex baseline hazards. To address
this, we propose Case-Base Neural Networks (CBNNs) as a new approach that
combines the case-base sampling framework with flexible neural network
architectures. Using a novel sampling scheme and data augmentation to naturally
account for censoring, we construct a feed-forward neural network that includes
time as an input. CBNNs predict the probability of an event occurring at a
given moment to estimate the full hazard function. We compare the performance
of CBNNs to regression and neural network-based survival methods in a
simulation and three case studies using two time-dependent metrics. First, we
examine performance on a simulation involving a complex baseline hazard and
time-varying interactions to assess all methods, with CBNN outperforming
competitors. Then, we apply all methods to three real data applications, with
CBNNs outperforming the competing models in two studies and showing similar
performance in the third. Our results highlight the benefit of combining
case-base sampling with deep learning to provide a simple and flexible
framework for data-driven modeling of single event survival outcomes that
estimates time-varying effects and a complex baseline hazard by design. An R
package is available at https://github.com/Jesse-Islam/cbnn.
Jesse Islam, Maxime Turgeon, Robert Sladek, Sahir Bhatnagar
# Case-Base Neural Networks: survival analysis with time-varying, higher-order interactions

###### Abstract

[Summary of Case-Base Neural Networks: survival analysis, machine learning, case-base, neural network

## 1 Introduction

Smooth-in-time accelerated failure time (AFT) models can estimate absolute risks by modeling the hazard directly through a user-specified baseline hazard distribution (Kleinbaum & Klein, 2012). Cox proportional hazards models are used more often than AFT models, causing analyses to be based on hazard ratios and relative risks rather than on survival curves and absolute risks (Hanley & Miettinen, 2009). The identification of an appropriate distribution for the baseline hazard in an AFT model may be difficult for common diseases that have many interacting risk factors, or a Cox model where the disease pathogenesis may change with age, making it difficult to maintain the proportional hazards assumption. For example, previous studies of breast cancer incidence have discovered time-varying interactions with covariates of interest, such as tumor size (Coradini et al.,
2000). One approach to provide flexibility in the baseline hazard involves using the basis of splines on time in our model (Royston and Parmar, 2002). However, regression-based models are limited in that they require prior knowledge about potential time-varying interactions and their quantitative effects.

Neural networks provide a data-driven approach to approximating interaction terms. For example, DeepSur is a neural network-based proportional hazards model that implements the Cox partial log-likelihood as a custom loss function (Katzman et al., 2018), resulting in a stepwise absolute risk curve that cannot accommodate time-varying interactions. Compared to Cox regression, DeepSurv shows better performance on the Study to Understand Prognoses Preferences and Risk Treatments (SUPPORT) dataset (Knaus et al., 1995). To handle non-proportional hazards, a modification to the loss function was proposed (Faraggi and Simon, 1995). To remove the need for tuning the number of layers and nodes, the concept of extreme learning machines has been applied to a Cox neural network model (Wang and Li, 2019).

As an alternative method that assumes a baseline hazard distribution, DeepHit specifies each survival time of interest in the model and directly estimates survival curves, rather than deriving a hazard function (Lee et al., 2018). It assumes an inverse Gaussian distribution as the baseline hazard and it outperforms DeepSur on the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) dataset (Curtis et al., 2012).

Providing further flexibility, Deep Survival Machines (DSM) is a parametric survival model using neural networks with a mixture of distributions as the baseline hazard (Nagpal et al., 2021). On both the SUPPORT and METABRIC datasets, DSM outperforms DeepSur and DeepHit (Nagpal et al., 2021). However, like DeepHit and DeepSur, DSM cannot model time-varying interactions.

We note that these alternative neural network approaches all require custom loss functions (Katzman et al., 2018) (Lee et al., 2018) (Nagpal et al., 2021). DeepHit introduces a hyperparameter weighing its two loss functions (negative log-likelihood and ranking losses) while DSM requires a two-phase learning process and user implementations for distributions beyond Log-Normal or Weibull (Lee et al., 2018) (Nagpal et al., 2021). Regression-based approaches require prior specification of all interaction terms, which makes it challenging to model covariate effects that change over time. The current neural network models provide flexibility at the cost of opacity, while regression models provide clarity at the cost of flexibility.

In this article, we propose Case-Base Neural Networks (CBNN) as a method that models time-varying interactions and a flexible baseline hazard using commonly available neural network components. Our approach to modeling the full hazard uses case-base sampling (Hanley and Miettinen, 2009). This sampling technique allows probabilistic models to predict survival outcomes. As part of the case-base framework, we use transformations of time as a feature (covariate) to specify different baseline hazards. For example, by including splines of time as covariates, we can approximate the Royston-Parmar flexible baseline hazard model (Royston and Parmar, 2002), however, this still requires explicit use of time-varying interactions. CBNN can model both without extra tuning parameters.

In Section 2, we describe how case-base sampling and neural networks are combined both conceptually and algebraically, along with our hyperparameter choices and software implementation. Section 3, describes our metrics and compares the performance of CBNN, DeepSur, DeepHit, DSM, Cox regression and case-base using logistic regression (CBLR) on simulated data. Section 4 describes the real-data analysis, while Section 5 explores the implications of our results and contextualizes them within neural network survival analysis in a single event setting.

## 2 Case-base neural network methodology, metrics and software

In this section, we define case-base sampling, which converts the total survival time into discrete person-specific moments (person-moments). Then, we detail how neural networks can be used within this framework, explicitly incorporating time as a feature while adjusting for the sampling bias. Finally, we report on the software versions used. An R package is available for use at [https://github.com/Jesse-Islam/cbnn](https://github.com/Jesse-Islam/cbnn). The entire code base to reproduce the figures and empirical results in this paper is available at [https://github.com/Jesse-Islam/cbnnManuscript](https://github.com/Jesse-Islam/cbnnManuscript).

### Case-base sampling

Case-base sampling is an alternative framework for survival analysis (Hanley and Miettinen, 2009). In case-base sampling, we sample from the continuous survival time of each person in our dataset to create a _base series_ of _person-moments_. This _base series_ complements the _case series_, which contains all person-moments at which the event of interest occurs.

For each person-moment sampled, let \(X_{i}\) be the corresponding covariate profile \((x_{i1},x_{i2},...,x_{ip})\), \(T_{i}\) be the time of the person-moment and \(Y_{i}\) be the indicator variable for whether the event of interest occurred at time \(T_{i}\). We estimate the hazard function \(h(t\mid X_{i})\) using the sampled person-moments. Recall that \(h(t\mid X_{i})\) is the instantaneous potential of experiencing the event at time \(t\) for a given set of covariates \(X_{i}\), assuming \(T_{i}\geq t\).
Now, let \(b\) be the (user-defined) size of the _base series_ and let \(B\) be the sum of all follow-up times for the individuals in the study. If we sample the base series uniformly across the study base, then the hazard function of the sampling process is equal to \(b/B\). Therefore, we have the following equality 1:

Footnote 1: We are abusing notation here, conflating hazards with probabilities. For a rigorous treatment, see Saarela and Hanley (2015) section 3 Saarela and Hanley (2015).

\[\frac{P\left(Y_{i}=1\mid X_{i},T_{i}\right)}{P\left(Y_{i}=0\mid X_{i},T_{i} \right)}=\frac{h\left(T_{i}\mid X_{i}\right)}{b/B}. \tag{1}\]

To provide some intuition for this equation, the hazard is the rate at which the event (\(Y_{i}=1\)) occurs at a given time, conditional on your survival up to the follow-up-time of interest. As we are looking at person-moments, we may retrieve the rate at which an event does not occur (\(Y_{i}=0\)) at a given time by taking the number of samples in \(b\) and dividing it by the total number of moments in the study base \(B\) at the current follow-up-time. by taking the ratio. of these two metrics, we end up with our odds on the left hand side. The odds of a person-moment being a part of the _case series_ is the ratio of the hazard \(h(T_{i}\mid X_{i})\) and the uniform rate \(b/B\). Though an intuitive approach is used to describe this equation, a complete derivation of this equation can be found at (Saarela and Hanley, 2015).

Sampling techniques are known to result in a loss of information. To address this concern, the original case-base methodology paper uses a measure for relative information when comparing two averages, based on the size of our base-series and case-series. If we let the size of \(b=100c\), then the variances and covariances are expected to be proportional to \(\frac{1}{c}+\frac{1}{100c}\) rather than \(\frac{1}{c}+\frac{1}{\ln t}\), about 1 percent larger. If there is a concern with a 1 percent inflation of these estimates, a larger ratio than 100 may be used. As there are concerns about dependency due multiple person-moments coming from each individual, the information loss was studied in Saarela 2015 ((Saarela and Hanley, 2015)), where they found a ratio of 100:1 base-series to case series gives an efficiency comparable to Mantel-Haenszel and conditional logistic regression methods (Saarela and Hanley, 2015).

\[\log\left(h\left(t\mid X_{i}\right)\right)=\log\left(\frac{P\left(Y_{i}=1 \mid X_{i},t\right)}{P\left(Y_{i}=0\mid X_{i},t\right)}\right)+\log\left( \frac{b}{B}\right). \tag{2}\]

To estimate the correct hazard function, we adjusting for the bias introduced when sampling a fraction of the study base \(B\) by adding an offset term \(\log\left(\frac{B}{b}\right)\) (2). Next, we propose using neural networks to model the odds.

### Neural networks to model the hazard function

After case-base sampling, we pass all features, including time, into any user-defined feed-forward neural network, to which an offset term is added, then passed through a sigmoid activation function (Figure 1). We use the sigmoid activation function as it is the inverse of the log-odds and commonly available in popular neural network packages (keras, pytorch, etc) (Allaire and Chollet, 2021)(Paszke et al., 2019). which we can use to calculate the hazard. The general form for the neural network using CBNN is:

\[P\left(Y=1|X,T\right)=\mathrm{sigmoid}\left(f_{\theta}(X,T)+\log\left(\frac{B }{b}\right)\right), \tag{3}\]

where \(T\) is a random variable representing the event time, \(X\) is the random variable for a covariate profile, \(f_{\theta}(X,T)\) represents any feed-forward neural network architecture, \(\log\left(\frac{B}{b}\right)\) is bias term set by case-base sampling, \(\theta\) is the set of parameters learned by the neural network and \(\mathrm{sigmoid}(x)=\frac{1}{1+e^{-x}}\). By approximating a higher-order polynomial of time using a neural network, the baseline hazard specification is now data-driven, where user-defined hyperparameters such as regularization, number of layers and nodes control the flexibility of the hazard function. We provide a detailed description of the choices we made in the next sub-section.
The following derivation shows how our probability estimate is converted to odds:

\[\log\left(h(t\mid X)\right) =\log\left(\frac{\operatorname{sigmoid}\left(f_{\theta}(X,T)+\log \left(\frac{B}{b}\right)\right)}{1-\operatorname{sigmoid}\left(f_{\theta}(X,T) +\log\left(\frac{B}{b}\right)\right)}\right)+\log\left(\frac{b}{B}\right)\] \[=\log\left(\frac{\frac{\exp\left(f_{\theta}(X,T)+\log\left(\frac{ B}{b}\right)\right)}{\exp\left(f_{\theta}(X,T)+\log\left(\frac{B}{b}\right) \right)+1}}{1-\frac{\exp\left(f_{\theta}(X,T)+\log\left(\frac{B}{b}\right) \right)}{\exp\left(f_{\theta}(X,T)+\log\left(\frac{B}{b}\right)\right)+1}}+ \log\left(\frac{b}{B}\right)\right.\] \[=\log\left(\exp\left(f_{\theta}(X,T)+\log\left(\frac{B}{b}\right) \right)\right)+\log\left(\frac{b}{B}\right)\] \[=f_{\theta}(X,T)+\log\left(\frac{B}{b}\right)+\log\left(\frac{b} {B}\right)\] \[=f_{\theta}(X,T).\]

We use binary cross-entropy as our loss function (Gulli and Pal, 2017):

\[L(\theta)=-\frac{1}{N}\sum_{i=1}^{N}y_{i}\cdot\log(\hat{f}_{\theta}(x_{i},t_{ i}))+(1-y_{i})\cdot\log(1-\hat{f}_{\theta}(x_{i},t_{i})),\]

where \(\hat{f}_{\theta}(x_{i},t_{i})\) is our estimate for a given covariate profile and time, \(y_{i}\) is our target value specifying whether an event occurred and \(N\) represents the number of individuals in our training set.

Backpropagation with an appropriate minimization algorithm (e.g. Adam, RMSPropagation, stochastic gradient descent) is used to optimize the parameters in the model (Gulli and Pal, 2017). For our analysis, we use Adam as implemented in Keras (Gulli and Pal, 2017). Note that the size of the _case series_ is fixed as the number of events, but we can make the _base series_ as large as we want. A ratio of 100:1 _base series_ to _case series_ is sufficient (Hanley and Miettinen, 2009). We pass our feed-forward neural network through a sigmoid activation function (Figure 1). Finally, we can convert this model output to a hazard. When using our model for predictions, we manually set the offset term to 0 in the new data, as we account for the bias during the fitting process.

Since we are directly modeling the hazard, we can readily estimate the risk function (\(F\)) at time \(t\) for a covariate profile \(X\), viz.

\[F\left(t\mid X\right)=1-\exp\left(-\int\limits_{0}^{t}h(u|X)\,\mathrm{d}u \right). \tag{4}\]

We use a finite Riemann sum (Hughes-Hallett, Gleason, & McCallum, 2020) to approximate the integral in (4).

Figure 1: Steps involved in CBNN from case-base sampling to the model framework we use for training. The first step is case-base sampling, completed before training begins. Next, we pass this sampled data through a feed-forward neural network. We add the offset and pass that through a sigmoid activation function, whose output is a probability. Once the neural network model completes its training, we can convert the probability output to a hazard, using it for our survival outcomes of interest.
### Hyperparameter selection

Neural networks are flexible when defining the architecture and optimization parameters. These hyperparameter decisions can affect the estimated parameters and were chosen during a set of initial simulations to determine if CBNN can learn complex interactions in practice. We set a max number of epochs to be \(2000\), batch size as \(512\), learning rate as \(10^{-3}\), decay as \(10^{-7}\), patience as \(10\) epochs, \(\{50,50,25,25\}\) nodes in each hidden layer with \(50\%\) dropout at each layer, a minimum delta loss on the validation set of \(10^{-7}\) over 10 epochs and Adam [12] as our optimizer. These choices may permit the model to approximate higher-order interactions while preventing over-fitting [13]. We fix a train-validation-test split that allows us to update the weights with a subset of the data (training), assess performance at each epoch (validation) and gauge performance of the final model (test) for each method. We select the best weights after training based on validation loss [12].

### Software implementation

R [14] and python [15] are used to evaluate methods from both languages. We fit the Cox model using the **survival** package [13], the CBLR model using the **casebase** package [16], the DeepSurv model using the **survivalmodels** package [17], the DeepHit model using **pycox** [18] and the DSM model using **DeepSurvivalMachines** [19]. We made the components of CBNN using the **casebase** package for the sampling step and the **keras**[1] package for our neural network architecture. The **simsurv** package [16] is used for our simulation studies, while **flexsurv**[1] is used to fit a flexible baseline hazard using splines for our complex simulation. We use the implementation of \(C_{IPCW}\) from the python package **sksurv**[12]. The **riskRegression** package [16] is used to get the Index of Prediction Accuracy (IPA metric). Both metrics are described in detail in the following section. We modify the **riskRegression** package to be used with any user supplied risk function \(F\). To ensure that both R and Python-based models are running in union on the same data through our simulations and bootstrap, we use the **reticulate** package [14].

## 3 Simulation studies

In this section, we use simulated data to evaluate the performance of CBNN and compare our approach with existing regression-based (Cox, CBLR) and neural network-based (DeepHit, DeepSurv, DSM) methods. We specify a linear combination of each covariate as the linear predictor in regression-based approaches (Cox, CBLR), which contrasts with neural network approaches that allow for non-linear interactions. We simulate data under a simple exponential model and a complex baseline hazard with time-varying interactions, each with 2000 individuals and 10% random censoring. For both settings, we simulate three covariates:

\[z_{1}\sim\mathrm{Bernoulli}(0.5) z_{2}\sim\begin{cases}N(0,0.5)&\text{if }z_{1}=0\\ N(10,0.5)&\text{if }z_{1}=1\end{cases} z_{3}\sim\begin{cases}N(8,0.5)&\text{if }z_{1}=0\\ N(-3,0.5)&\text{if }z_{1}=1.\end{cases}\]

The DeepHit-specific hyperparameter alpha is set to \(0.5\) (equal weight between its negative log-likelihood and ranking losses [18]). We modify the **DeepSurvivalMachines**[19] package to include dropout and a minimum delta loss during the fitting process. For DSM, we define a mixture model of six Weibull distributions for the baseline hazard. All other hyperparameters are held constant across all neural network methods in both the simulation studies and real data applications.

Besides the methods mentioned above, we include the Optimal model in our comparisons using CBLR. That is, we include the exact functional form of the covariates in a CBLR model (referred to as Optimal for simplicity). We calculate \(t\)-based 95% confidence intervals using 100 replications of the simulated data. For all analyses, we use 80% for training and 20% for the test set. 20% of the training set is kept for validation at each epoch. We predict risk functions \(F\) using (4) for individuals in the test set, which are used to calculate our \(C_{IPCW}\) and IPA scores.

### Performance metrics

We use two metrics to assess the performance of the different methods of interest: 1) (IPA) [16] and 2) inverse probabililty censoring weights-adjusted concordance index (\(C_{IPCW}\)) [10], which we define below. We may be interested in an absolute risk a specific number of years in the future. For the sake of comparison, we plot each metric over time \(t\) so that the performance is transparent. As we go further along follow-up-time, we expect the performance of all models to drop as there are fewer individuals left in the study by then. For \(C_{IPCW}\), the score demonstrates performance up to time \(t\)[10], while IPA is performance at time \(t\).
#### 3.1.1 Index of prediction accuracy (IPA)

The IPA is a function of the Brier score (\(BS(t)\)) [11], which is defined as

\[BS(t)=\frac{1}{N}\sum_{i=1}^{N}\left(\frac{\left(1-\widehat{F}(t\mid X_{i}) \right)^{2}\cdot I(T_{i}\leq t,\delta_{i}=1)}{\widehat{G}(T_{i})}+\frac{ \widehat{F}(t\mid X_{i})^{2}\cdot I(T_{i}>t)}{\widehat{G}(t)}\right), \tag{5}\]

where \(\delta_{i}=1\) shows individuals who have experienced the event, \(N\) represents the number of samples in our dataset over which we calculate \(BS(t)\), \(\widehat{G}(t)=P[c>t]\) is a non-parametric estimate of the censoring distribution, \(c\) is censoring time and \(T_{i}\) is an individual's survival or censoring time. The Brier score provides a score that accounts for the information loss because of censoring. There are three categories of individuals that may appear within the dataset once we fix our \(t\) of interest. Individuals who experienced the event before \(t\) are present in the first term of the equation. The second term of the equation includes individuals who experience the event or are censored after \(t\). Those censored before \(t\) are the third category of people. The inverse probability censoring weights (IPCW) adjustment (\(G(\cdot)\)) is to account for these category three individuals whose information is missing. The IPA score as a function of time is given by

\[\mathrm{IPA}(t)=1-\frac{BS_{model}(t)}{BS_{null}(t)},\]

where \(BS_{model}(t)\) represents the Brier score over time \(t\) for the model of interest and \(BS_{null}(t)\) represents the Brier score if we use an unadjusted Kaplan-Meier (KM) curve as the prediction for all observations [10]. Note that IPA has an upper bound of one, where positive values show an increase in performance over the null model and negative values show that the null model performs better. These scores show how performance changes over follow-up time.

A potential artifact of IPA is that the score is unstable at earlier and later survival times. This is because of a near equivalent Brier score among each model and the null model. At small values, a difference of \(0.1\) creates a more significant fold change than at larger values. As the Brier score is potentially small at the start and at the end of their respective curves, The IPA score may be unstable at the same locations.

#### 3.1.2 Inverse probability censoring weights-adjusted concordance index

The \(C_{IPCW}\) is a non-proper, rank-based metric that does not depend on the censoring times in the test data [12]. The \(C_{IPCW}\) is given by

\[C_{IPCW}(t)=\frac{\sum_{i=1}^{N}\sum_{j=1}^{N}\delta_{i}\left\{\widehat{G}(T_ {i})\right\}^{-2}I(T_{i}<T_{j},T_{i}<t)I\left(\widehat{F}(t|X_{i})>\widehat{F }(t|X_{j})\right)}{\sum_{i=1}^{N}\sum_{j=1}^{N}\delta_{i}\left\{\widehat{G}(T _{i})\right\}^{-2}I(T_{i}<T_{j},T_{i}<t)}. \tag{6}\]

where the follow-up period of interest is (0,\(t\)), \(I(\cdot)\) is an indicator function, \(\widehat{F}(X_{i},t)\) is the risk function estimated for everyone in the study at time \(t\) and \(C_{IPCW}\) can compare the performance of different models, where a higher score is better. Note that the \(C_{IPCW}\) may produce misleading performance, as it ranks based on survival times, not event status [1]. This metric is considered an unbiased population concordance measure because of the IPCW adjustment [12].

### Simple simulation: constant baseline hazard

We simulate data from a simple model that primarily depends on a constant baseline hazard:

\[\log h(t\mid X_{i})=\beta_{1}z_{1}+\beta_{2}z_{2}+\beta_{3}z_{3},\]

The covariate effects are given by \(\beta_{1}=0.1,\beta_{2}=0.1,\beta_{3}=0.1\). Once we simulate survival times, we introduce 10% random censoring.

#### 3.2.1 Performance comparison in simple simulation

Figure 2 A, B and Table 1 A show the results for the simple simulation. The regression-based methods (CBLR, Cox, Optimal) outperform the neural network ones in the simple simulation setting. Among the neural network approaches, CBNN outperforms all other methods in terms of both IPA and \(C_{IPCW}\) (Figure 2 A, B). Specifically, we see CBNN is consistent across time with smaller confidence intervals compared to DeepHit, DeepSurv and DSM.

In this simple setting, the regression models are much closer to the Optimal model, while the neural network models perform worse than the KM null model. The wide confidence bands in \(C_{IPCW}\) suggest the neural network models may be over-parameterized.
### Complex simulation: flexible baseline hazard, time-varying interactions

This simulation demonstrates performance with the presence of a complex baseline hazard and a time-varying interaction. Originally used to show the spline-based hazard model proposed by Royston and Parmar (Royston and Parmar 2002), the breast cancer dataset provides a complex hazard from which we simulate, available in the **flexsurv** R package (Jackson 2016). To increase the complexity of our data-generating mechanism for this simulation, we design the model as follows:

\[\log h(t\mid X_{i})=\sum_{i=1}^{5}(\gamma_{i}\cdot\psi_{i})+\beta_{1}(z_{1})+ \beta_{2}(z_{2})+\beta_{3}(z_{3})+\tau_{1}(z_{1}\cdot z_{2}\cdot t)+\tau_{2}(z_ {1}\cdot z_{3})+\tau_{3}(z_{2}\cdot z_{3}),\]

where \(\gamma_{1}=3.9,\gamma_{2}=3,\gamma_{3}=-0.43,\gamma_{4}=1.33,\gamma_{5}=-0.86, \beta_{1}=1,\beta_{2}=1,\beta_{3}=1,\tau_{1}=10,\tau_{2}=2,\tau_{3}=2\) and \(\psi\) are basis splines. The \(gamma\) coefficients are obtained from an intercept-only cubic splines model with three knots using the _flexsurv spline_ function from the **flexsurv** package (Jackson 2016). Note that we fix these values for the analysis. The \(\beta\) coefficients represent direct effects, \(\tau_{2}\) and \(\tau_{3}\) represent interactions and \(\tau_{1}\) is a time-varying interaction.

#### 3.3.1 Performance comparison in complex simulation

Figure 2 C, D and Table 1 B show the performance over time on a test set in the complex simulation. Apart from the Optimal regression model, CBNN outperforms the competing models when examining IPA and up to the 75-percentile of follow-up time for \(C_{IPCW}\). We expected the Optimal model to perform best in both metrics. However, this was not the case for \(C_{IPCW}\) and may be due to an artifact of concordance-based metrics, where a misspecified model may perform better than a correctly specified one (Blanche et al. 2019). We attribute the performance of CBNN to its flexibility in modeling time-varying interactions and baseline hazard, flexibility the other neural network models do not have.

Figure 2: Summarizes the simple simulation (A, B), complex simulation (C, D), SUPPORT case study (E, F) and METABRIC case study (G, H) results. The first row shows the IPA score for each model in each study over follow-up time. Negative values mean our model performs worse than the null model and positive values mean the model performs better. The second row demonstrates the \(C_{IPCW}\) score for each model in each study over follow-up time. A score of 1 is the maximum performance for either metric. Each model-specific metric in each study shows a 95-percent confidence interval over 100 iterations. The models of interest are case-base with logistic regression (CBLR), Case-Base Neural Networks (CBNN), Cox, DeepHit, DeepSurv, Deep Survival Machines (DSM), Optimal (a CBLR model with the exact interaction terms and baseline hazard specified) and Kaplan-Meier (to serve as a baseline, predicting the average for all individuals). Note than Cox and CBLR are superimposed, due to very similar prediction outcomes.
\begin{table}
\begin{tabular}{c c c c c c c c c} \hline \hline \multicolumn{1}{c}{**A: Simple**} & \multicolumn{3}{c}{IPA} & \multicolumn{3}{c}{c} \multicolumn{3}{c}{c-ipcw} \\ \hline Method & 25 & 50 & 75 & 100 & 25 & 50 & 75 & 100 \\ \hline Cox & **0.05 (0.05,0.06)** & **0.07 (0.07,0.08)** & **0.08 (0.07,0.08)** & **0.07 (0.06,0.08)** & **0.0 (0.39,0.6)** & **0.0 (0.39,0.6)** & **0.0 (0.39,0.6)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** \\ CBLR & **0.05 (0.05,0.06)** & **0.07 (0.07,0.08)** & **0.08 (0.07,0.08)** & **0.07 (0.06,0.08)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** \\ DeepSurv & -0.26 (0.27,0.25) & -0.36 (0.38,0.34) & -0.22 (0.24,0.21) & 0 (0.01,0) & 0.51 (0.50,52) & 0.51 (0.50,5.52) & 0.51 (0.50,5.52) & 0.51 (0.50,5.52) \\ DeepHt & -0.03 (0.04,0.02) & -0.05 (0.07,0.04) & -0.03 (0.04,0.02) & -0.01 (0.01,0) & 0.52 (0.51,0.52) & 0.52 (0.51,0.52) & 0.52 (0.51,0.53) & 0.50 (0.49,0.5) \\ CBNN & -0.03 (0.04,0.03) & -0.04 (0.05,0.04) & -0.04 (0.05,0.03) & -0.03 (0.04,0.03) & 0.54 (0.54,0.55) & 0.54 (0.54,0.55) & 0.54 (0.54,0.55) & 0.54 (0.54,0.55) & 0.54 (0.54,0.55) \\ Optimal & **0.05 (0.05,0.06)** & **0.07 (0.07,0.08)** & **0.08 (0.07,0.08)** & **0.07 (0.06,0.08)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** & **0.0 (0.59,0.6)** \\ DSM & -0.28 (0.3,0.27) & -0.4 (0.43,0.37) & -0.28 (0.3,0.25) & 0.01 (0.01,0.02) & 0.55 (0.54,0.56) & 0.55 (0.54,0.56) & 0.55 (0.55,0.56) \\ \hline \hline \multicolumn{1}{c}{**B: Complex**} & \multicolumn{3}{c}{IPA} & \multicolumn{3}{c}{c-ipcw} \\ \hline Method & 25 & 50 & 75 & 100 & 25 & 50 & 75 & 100 \\ \hline Cox & 0.55 (0.54,0.56) & 0.45 (0.44,0.46) & -0.19 (0.21,0.17) & -0.51 (0.54,0.48) & 0.77 (0.76,0.79) & 0.74 (0.72,0.76) & 0.56 (0.54,0.58) & 0.54 (0.52,0.56) \\ CBLR & 0.53 (0.52,0.54) & 0.45 (0.44,0.46) & -0.15 (0.16,0.13) & 0.46 (0.48,0.43) & 0.79 (0.78,0.79) & 0.79 (0.78,0.79) & 0.79 (0.78,0.79) & **0.78 (0.76,0.79)** \\ DeepSurv & 0.03 (0.01,0.05) & 0.06 (0.05,0.07) & -0.27 (0.29,0.25) & -0.04 (0.07,0.02) & 0.68 (0.56,0.7) & 0.68 (0.56,0.7) & 0.68 (0.66,0.7) & 0.68 (0.66,0.7) \\ DeepHt & 0.26 (0.22,0.31) & 0.07 (0.06,0.00) & -0.06 (0.08,0.04) & 0.05 (0.04,0.06) & 0.72 (0.70,7.04) & 0.65 (0.62,0.67) & 0.37 (0.33,0.4) & 0.35 (0.33,0.30) \\ CBNN & 0.83 (0.81,0.84) & 0.69 (0.67,0.71) & 0.49 (0.45,0.52) & 0.56 (0.53,0.59) & **0.86 (0.85,0.86)** & **0.87 (0.87,0.88)** & **0.8 (0.79,0.81)** & 0.67 (0.65,0.68) \\ Optimal & **0.94 (0.94,0.95)** & **0.9 (0.9,0.91)** & **0.84 (0.84,0.85)** & **0.85 (0.84,0.85)** & 0.74 (0.74,0.75) & 0.82 (0.82,0.82) & 0.63 (0.63,0.63) & 0.57 (0.56,0.57) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Four tables representing performance at certain follow-up times for the simple simulation, complex simulation, SUPPORT and METABRIC. Each table shows performance for each method in each study at \(25\%\), \(50\%\). \(75\%\) and \(100\%\) of follow-up time. The bold elements show the best model for each study, at each follow-up time of interest. These tables are included to provide exact measures at certain intervals. The models of interest are: Cox, case-base with logistic regression (CBLR), DeepSurv, DeepHt, Case-Base Neural Network (CBNN), Optimal and Deep Survival Machines (DSM).
## 4 Application to Support and Metabric Data

Our complex simulation demonstrates the superior performance of CBNN in ideal conditions with clean data. To obtain a more realistic performance assessment, we compared models using two real datasets with a time-to-event outcome. The first case study examines the SUPPORT dataset (Knaus et al., 1995). The second case study examines the METABRIC dataset (Curtis et al., 2012). We use the same hyperparameters as in the simulation studies. As we do not know the true model for the real data, we exclude the Optimal model. We split the datasets keeping 20% of the observations as a test set. 20% of the training set is kept aside for validation at each epoch. We predict risk functions for everyone in the test set, which is used to calculate our metrics. We conduct 100 bootstrap re-samples for the real data applications to obtain confidence intervals.

### Performance evaluation using the SUPPORT dataset

The SUPPORT dataset tracks the time until death for seriously ill patients at five American hospitals (Knaus et al., 1995). We use a pre-processed version of the dataset made available in the DeepSurv package (Katzman et al., 2018). This dataset contains 9104 samples and 14 covariates (age, sex, race, number of comorbidities, presence of diabetes, presence of dementia, presence of cancer, mean arterial blood pressure, heart rate, respiration rate, temperature, white blood cell count, serum's sodium and serum's creatinine) (Katzman et al., 2018). Patients with missing features were excluded and 68.10% of the patients died during the 5.56-year study period (Katzman et al., 2018).

Figure 2 E, F and Table 1 C demonstrates the performance over time on a test set. The regression models (CBLR, Cox) perform best considering IPA, followed by CBNN from the \(25^{th}\) to \(100^{th}\) percentile of follow-up time. We note a drop in performance for CBNN before the \(25^{th}\) percentile of follow-up. For \(C_{IPCW}\), CBNN outperforms the competing models consistently over follow-up time. Note that performance is similar for all models, aside from DeepSurv whose \(C_{IPCW}\) is lower than the rest (Figure 2 E, F and Table 1 C).

### Performance evaluation using the METABRIC dataset

METABRIC is a 30-yearlong study aiming to discover the molecular drivers of breast tumors, following 2000 individuals with breast cancer until death (Curtis et al., 2012). They described these growths as primary invasive breast carcinomas, with the goal of discovering both genetic and clinical risk factors for breast cancer survival (Curtis et al., 2012). We used the processed dataset made available through DeepSurv (Katzman et al., 2018), which includes 1980 patients of which 57.72% die due to breast cancer within a median 10 years of follow-up (Katzman et al., 2018). There are 10 covariates in total: time, 4 RNA-Seq gene expressions (MK167, EGFR, PGR and ERBB2) and 5 clinical features (hormone treatment indicator, radiotherapy indicator, chemotherapy indicator, ER-positive indicator and age at diagnosis) (Katzman et al., 2018).

Figure 2 G, H and Table 1 D shows the performance on a test set over time. The IPA scores suggest that regression models outperform competing models on this dataset, as all the neural network models are equal to or perform worse than KM over follow-up time. Our CBNN model is comparable to KM until around the 50-percentile of follow-up time, after which CBNN, DeepSurv and DSM drop in performance. The \(C_{IPCW}\) produces a different ranking. Our CBNN model outperforms the other models up to the \(25^{th}\) percentile of follow-up time, whereas DeepHit performs best from the \(50^{th}\) to \(75^{th}\) percentile of follow-up. With this metric, CBNN and DeepHit outperform the regression models. This disagreement between IPA and \(C_{IPCW}\) may be due to the misspecification issue of concordance-based metrics (Blanche et al., 2019). The neural network models may be over-parameterized, as shown by the wide confidence bands in \(C_{IPCW}\).

## 5 Discussion

CBNN models survival outcomes by using neural networks on case-base sampled data. We incorporate follow-up time as a feature, providing a data-driven estimate of a flexible baseline hazard and time-varying interactions in our hazard function. The three competing neural network models we evaluated cannot model time-varying interactions by design (Nagpal et al., 2021)(Katzman et al., 2018)(Lee et al., 2018). DSM requires a mixture component distributions for a flexible baseline hazard to be fit (Nagpal et al., 2021). As our goal is to fix the design and compare performance, we did not change any shared hyperparameters. With our choice of shared hyperparameters and model design, DSM did not converge in the complex simulation and SUPPORT case study. Despite this limitation, we include this method when it did function as it can fit flexible baseline hazards. Compared to CBNN, the remaining two models also have limitations. DeepSurv is a proportional hazards model and does not estimate the baseline hazard (Katzman et al., 2018). DeepHit requires an alpha hyperparameter, is restricted to a single distribution for the baseline hazard and models the survival function directly (Lee et al., 2018). The alternative neural network methods match on time, while CBNN models time directly.

To assess performance among these models, we use both IPA and \(C_{IPCW}\) metrics. Concordance-based measures are commonly used in survival analysis to compare models and we opt to keep them in our analyses. However, \(C_{IPCW}\) is a non-proper metric and may cause misspecified models
to appear better than they should (Blanche et al., 2019). Therefore, we contextualize our \(C_{IPCW}\) results in relation to IPA, a proper scoring rule. The model rankings between \(C_{IPCW}\) and IPA differed for both the complex simulation and METABRIC application. Wider \(C_{IPCW}\) confidence intervals for DeepHit and DeepSurv show potential misspecification of the models (Figure 2 D, H).

With this interpretation of \(C_{IPCW}\) in mind, we assess two simulations with minimal noise and two case studies in real scenarios. The simple simulation demonstrates potential pitfalls associated with neural network models, particularly overparameterization. All neural network-based approaches performed worse than null KM model (IPA score), while the regression-based performed better. We attribute this to potential overparameterization in the neural network models, as the wide confidence intervals suggest over-fitting, even with dropout. Both DeepSurv and DSM are affected, while DeepHit and CBNN are less so. This is a limitation to our strategy for evaluating the methods, which uses a fixed study design across all assessments.

From this baseline benchmark, we move on to a complex simulation that requires a method that can learn both time-varying interactions and have a flexible baseline hazard. Here, CBNN demonstrates a distinct advantage over all other methods. The regression models show improved performance over the null KM model, while the competing neural network models performed worse. Based on our complex simulation results (Figure 2 C, D and Table 1 B), CBNN outperforms the competitors when time-varying interactions and a complex baseline hazard are present. This simulation shows how CBNN can perform under ideal conditions, while the following two analyses on real data serve to assess its performance in realistic conditions.

In the SUPPORT and METABRIC case studies, flexibility in both interaction modeling and baseline hazard did not improve CBNNs relative performance, suggesting that this flexibility does not aid prediction in either case study. From a biological perspective, the 30-year follow-up time in the METABRIC study may contain competing causes of death. The causes at the start of the study may not match the causes towards the end, potentially explaining the drop in performance as we reach later survival times with fewer individuals. The baseline hazard is unlikely to cause this drop as DSM is the most flexible competing model in our comparison. Over-fitting is also unlikely given the tight confidence intervals. Further research is required to determine the cause of the drop in performance seen in the METABRIC case study. In both case studies, CBNN outperforms the competing neural network methods.

The way neural networks are used in this paper are not directly interpretable. This is by design, as we only wish to compare to competitors based on predictive ability. All feed-forward neural network models can make use of different architectures to improve performance or interpretability. While a deep neural network as used in this paper is not directly interpretable, techniques such as Local Interpretable Model-Agnostic Explanations (LIME) (Ribeiro, Singh, & Guestrin, 2016) may be used to interpret the relevance of each feature to a given prediction. An example of architectural interpretability is a convolutional neural network for images, which provide positional information about the relevant pixels in an image (Allaire & Chollet, 2021). Note that neither of these are specific to CBNN. As CBNN is concerned with learning time-varying interactions, a concern may be interpreting which time-varying interactions have been learned. Other secondary methods after the fitting process may be used to retrieve said information [https://ui.adsabs.harvard.edu/abs/2017arXiv1705049777/abstract](https://ui.adsabs.harvard.edu/abs/2017arXiv1705049777/abstract).

## 6 Conclusions

CBNN outperforms all competitors in the complex simulation, demonstrating its value in survival settings that may involve time-varying interactions and a complex baseline hazard. Once we perform case-base sampling and adjust for the sampling bias, we can use a sigmoid activation function to predict our hazard function. Our approach simplifies the incorporation of censored individuals, allowing survival outcomes to be treated as binary ones. Forgoing the requirement of custom loss functions, CBNN only requires the use of standard components in machine learning libraries (specifically, the add layer to adjust for sampling bias and the sigmoid activation function). Due to the simplicity in its implementation and by extension user experience, CBNN is both a user-friendly approach to data-driven survival analysis and is easily extendable to any feed-forward neural network framework.

### Data and code availability statement

The pre-processed data for the SUPPORT case study can be found at [https://github.com/jaredleekatzman/DeepSurv/tree/master/experiments/data/support](https://github.com/jaredleekatzman/DeepSurv/tree/master/experiments/data/support). The pre-processed data for the METABRIC case study can be found at [https://github.com/jaredleekatzman/DeepSurv/tree/master/experiments/data/metabric](https://github.com/jaredleekatzman/DeepSurv/tree/master/experiments/data/metabric). The data is accessed using [https://github.com/havakv/pycox](https://github.com/havakv/pycox). The code for this manuscript and its analyses can be found at [https://github.com/Jesse-Islam/cbnnManuscript](https://github.com/Jesse-Islam/cbnnManuscript). The software package making CBNN easier to use can be found at [https://github.com/Jesse-Islam/cbnn](https://github.com/Jesse-Islam/cbnn).
## Acknowledgements

We would like to thank Dr. James Meigs, the project leader of these awards, for his support and helpful discussions. The work was also supported as part of the Congressionally Directed Medical Research Programs (CDMRP) award W81XWH-17-1-0347. We would also like to thank Dr. James Hanley for his support and discussions while extending the case-base methodology.

## Financial disclosure

This work was supported by subcontracts from UM1DK078616 and R01HL151855 to Dr. Rob Sladek.

## Conflict of interest

The authors declare no potential conflict of interests.
