NeuralSurv: Deep Survival Analysis with Bayesian Uncertainty Quantification
Authors: Mélodie Monod
Alessandro Micheli
Samir Bhatt
Venue: NeurIPS 2025 poster ()
Keywords: Survival analysis
Variational Inference
Deep Learning
Bayesian Neural Networks
Local Linearization
Polya-Gamma
Primary area: probabilistic_methods

Abstract: We introduce *NeuralSurv*, the first deep survival model to incorporate Bayesian uncertainty quantification. Our non‑parametric, architecture‑agnostic framework flexibly captures time‑varying covariate–risk relationships in continuous time via a novel two‑stage data‑augmentation scheme, for which we establish theoretical guarantees. For efficient posterior inference, we introduce a mean‑field variational algorithm with coordinate‑ascent updates that scale linearly in model size. By locally linearizing the Bayesian neural network, we obtain full conjugacy and derive all coordinate updates in closed form. In experiments, *NeuralSurv* delivers superior calibration compared to state-of-the-art deep survival models, while matching or exceeding their discriminative performance across both synthetic benchmarks and real-world datasets. Our results demonstrate the value of Bayesian principles in data‑scarce regimes by enhancing model calibration and providing robust, well‑calibrated uncertainty estimates for the survival function.

NeuralSurv: Deep Survival Analysis with Bayesian Uncertainty Quantification
Mélodie Monod"
Alessandro Micheli' Imperial College London London, United Kingdom
Imperial College London
London, United Kingdom
melodie monod Gemnod@@@imeeriee ac uk
a. mi mnicheei19@imperaall ac uk
Samir Bhatt
Imperial College London; University of Copenhagen London, United Kingdom; Copenhagen, Denmark s. bhatt@imperial.. acc.a.lll............. uk
Abstract
We introduce NeuralSurv, the first deep survival model to incorporate Bayesian uncertainty quantification. Our non-parametric, rrchttetture-gnootiic framework captures time-varying covariate-risk relationships in continuous time via a novel two-stage Sata-augmentation scheme, for which we establish theoretical guarantees. For efficient posterior inference, we introduce mean-field variational algorithm with coordinate- ing the Bayesian neural network, we obtain full conjugacy and derive all coordinate updates in closed form. In experiments, NeuralSurv delivers superior calibration compared to state-of-the- their discriminative performance across both synthetic benchmarks and real-world datasets. Our results demonstrate the value of Bayesian principles in data-scarce regimes by enhancing model calibration and providing robust, well-calibrated uncertainty estimates for the survival function.
updates that scale linearly in model size. By locally lineariz-
deep survival models, while matching or exceeding
Introduction
Survival analysis branch of statistics focused the study of time-to-event data, usually called event times. This type of data appears in a wide range of applications such as medicine [33], engineering [35], and social sciences [42]. A key objective of survival analysis is estimate the hazard function and the survival function that govern the distribution of event times.
Traditional survival models like the Cox proportional hazards model [10] and accelerated failure time models [8] have long delivered reliable inference under strong parametric assumptions. However, such assumptions may fail adequately capture complex and evolving baseline hazards, especially when risk relationships vary over time. To overcome these limitations, recent work has begun incorporating modern machine-learning techniques [48], and in partieular deep architectures [49, 24, 32]. which can learn rich, hierarchical representations directly from data. Yet most deep-survival approaches remain purely frequentist, optimizing point-estimate losses and offering coherent uncertainty quantification. high-stakes settings like medicine, this lack of reliable uncertainty estimates can undermine trust and impede adoption.
Bayesian statistics, by contrast, inherently quantifies uncertainty: prior beliefs are combined with ob- served data yield posterior distribution model parameters [15]. In survival analysis, Bayesian
'Equal contribution.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

methods can produce full posterior distributions for individual survival functions summarizable via credible intervals that communicate model confidence [21]. Traditional Bayesian survival tools, such as Gaussian processes (GPs) [13, 25], offer nonparametric flexibility and built-in uncertainty but often falter in high-dimensional settings due to scalability issues. To date, no method has combined the representational power of deep learning with full Bayesian uncertainty quantification in scalable survival framework. Such a synthesis would hold the potential to learn complex, high-dimensional survival dynamics while retaining principled probabiliiticc interpretations.
In this work, we introduce NeuralSurv, an architecture aggnostic. Bayesian deep-learning framework for survival analysis which integrates with modern deep learning architectures. NeuralSurv leverages deep Neural Networks (NNs) to learn hierarchical representations from covariates and uses a princi- pled variational inference framework to provide rigorous uncertainty quantification over the survival function. We develop two-stage data-auggentatiin strategy using latent marked Poisson processes and Pólya-Gamma variables to enable exact continuous-time likelihood computation, and provide novel theoretical guarantees for this approach. By locally linearizing the Bayesian Neural Network (BNN), we achieve conjugacy and derive closed-form coordinate-ascent updates that scale linearly with network size.
Through extensive experiments on synthetic and real survival datasets, in data-scarce settings, Neu- ralSurv consistently delivers superior calibration compared to of-e-of-theear deep survival models, and matches or exceeds their discriminative performance. Its Bayesian formulation captures epis- temic uncertainty to prevent overfitting, while informative priors induce a soft regularization that yields smooth, plausible survival functions. The code to reproduce our experiments is available on the GitHub repository https:/ //gp::://iiii http://gtcuu.c.cm........... www/MLGGoballaeeee wateeealltteeeussuurrruu under the MIT License.
2 NeuralSurv
In this section, we outline the main assumptions underlying NeuralSurv. We begin by briefly reviewing key concepts in survival analysis. Survival analysis focuses on modeling time-to-event data. Let T be a continuous nonnegative random variable with probability density f and cumulative distribution function F, representing the time until particular event occurs. Its survival function S(t)= P(T > t) = 1 - F(t) gives the probability of not experiencing the event by time t, while the functionn
hazard function X(I) = [(L)/((L) represents the instantaneous risk of the event at time conditional
on having survived up to time L. In practice, the event time may not be observed for all individuals, because some observations are subject to right-censoring, where the event has not yel occurred by the end of the observation period. For each observation i == l,.... N, denote the event time by T, and the censoring time by C.. We observe y min(Ti, C.) and ó; (T, where y represents the observed time (which may correspond either to the event or to censoring). and ò; indicates whether the event time was observed (6, 1) or the observation period was censored (6, 0). We assemble the dataset as D = (1,, 5,) i = 1, .N. Each observation also carries a covariate vector x; E R", gathered into X = {x,
1, 1.... N}. Throughout this paper, we assume that the censoring time C is independent of the event time T; given x; (known as non- non-informative censoring). Further details on survival analysis theory are provided in Appendix A.
2.1 Sigmoidal Hazard Function
Our goal is to model the hazard function X, i.e. the instantaneous event rate at time conditional on survival to and covariates x. We employ the sigmoid function az)) 1/(1 + exp(-2)), which
maps real-valued inputs to the interval (0, 1). The sigmoidal hazard model is constructed as the product of normalized baseline hazard function (Ap) and a modulation function (ar):
A(t x; c. g(: 0)) == Ao x; Φ) (g(t, x: 0)),
(I)
where the normalized baseline hazard is given by
Ao x; 0) = Ao Φ)
(2)
2(1,
for the baseline hazard Ag RR + R4, parametrized by & E RR,, and a normalization factor 7(1, x) that depends on both time and covariates. The term (t, x; Φ) encodes prior "best-guess" hazard

---

profile over time. The flexible function g R4 R' + R. parametrized by 0 E R provides data-driven adjustment: once passed through the sigmoid, multiplicatively attenuates the baseline hazard, continuously scaling it between 7ero and λg. The normalization factor (1,, x) ensures that the overall hazard remains properly sealed after modulation by the sigmoidal function (see Section 2.3 for details). Modeling hazard and intensity functions using sigmoidal transformations is common in applications such as survival analysis [13] and point process models [52, 12, 45, 1, 2]. This approach popular due to the balance it offers between modeling flexibility and analytical tractability.
2.2 Likelihood Distribution
Given the hazard function in (1), the likelihood density for the observation corresponding to the 7th observation is given by:
p(yi, 8 xx, 6, (::
(do(u), x;: (g(n,, x;; 0)) exp (-/," 2o(1, x;; ( x;; Ox;d)) (3)
Assuming (yi, 6,) are i.i.d. conditional on (xi, s( 0)), the full-sample likelihood simply the product
p(D X. x.0,((:0)
II p x,, 9(;;0),
2.3 Prior Distributions
Prior Distribution on 0. We assume that g(: 0) is BNN parameterized by 0. Furthermore, denote by I the m x m identity matrix. We place the following isotropic Gaussian prior with zero mean and identity covariance over the NN weights
(0 00
(5)
This common choice [6] assumes weights are independently distributed and centered around zero, acting as uninformative yet regularizing prior that discourages large weights and helps prevent overfitting via shrinkage.
Prior Distribution on é. We adopt a Weibull-type baseline hazard
(t Φ) = dtP-1 (0( Gamma(au, Bo), p > fixed,
(6)
where ag is the shape and 3o is the rate of the Gamma distribution. The Weibull-type baseline hazard (6) is the hazard of a Weibull distribution, a common choice in survival analysis [13]. When 1, Ao(t; o) becomes constant and the baseline hazard reduces to the hazard of the Exponential distribution.
Normalization Factor Z. We define the normalization factor introduced in (2)
(t, x) := Eg (o(g(t, 0))]
and refer the reader Appendix D for further details how it is computed. Introducing this normalization factor ensures that the prior mean of the sigmoidal hazard in (1) coincides with the prior mean of the baseline hazard, i.e.
E
[A(t x; 90: = Egpps [oo ø)]
This approach, similar to the technique used [13], centers distribution around the baseline hazard Ao (), favouring hazard trajectories that remain close to this prior "best-guess" profile while still permitting data-driven deviations, Notice that if g(; 0) has fully connected architecture, then Z(t, x) = for all (t, x), resulting in the same normalization factor value as in [13].

---

2.4 Posterior Distribution
Let p(ó, 0 D, X) denote the posterior density over the parameters and 0, defined with respect to the product measure do x d0. By Bayes' rule, this posterior is proportional (up to normalization) to
(0. θ X) x p(D X, é, g(: 0)) po( po (o
(7)
The posterior in (7) is generally intractable to compute for three reasons. First, its normalization constant unavailable in closed form. Second, the likelihood from (3) requires evaluating N integrals, none of which admits an analytic solution. Finally, the sigmoid in (1) introduces an extra nonlinearity, rendering inference even more analytically challenging.
3 Data Augmentation Strategy
In this section, we present data augmentation scheme that leverages the properties of Poisson processes and Pólya-Gamma random variables. Specifically, Poisson processes help overcome the challenges associated with computing the integrals of the continuous-ting time function to evaluate the likelihood, while the Pólya-Gamma random variables allow for exact handling of the sigmoid nonlincarity without relying on analytic approximations. This combined approach allows us to efficiently perform posterior inference from the model without resorting discretization.
This approach builds on analogous strategies previously applied other settings, including Bayesian inference for Sigmoid Gaussian Cox Processes [12], nonparametric Hawkes processes [52], and, in the case of Pólya-Gamma augmentation alone, mutually regressive point processes [4]. To the best of our knowledge, this is the first application of such data augmentation strategy in the context of survival analysis. Furthermore, we are the first to provide rigorous theoretical framework that establishes the validity of a method belonging to this broader class of augmentation- approaches (see Theorem 3.1).
Detailed reviews of Pólya-Gam Gamma random variables and Poisson processes are provided in Appen- dices B and C, respectively.
3.1 Polya-Gammn Augmentation Scheme
A primary challenge in model arises from the sigmoid function, whose inherent nonlinearity complicates the posterior inference. To overcome this, we adopt the scheme introduced in [38]. The key insight of this approach is that the sigmoid function can be represented in terms of Pólya-Gamma random variables. Define the function
data augmentation
x²
f(w,)) = 2 2
log(2).
(8)
Then, the following identity holds:
(()) Lo (((,,)) pra(w 0),,,,
(9)
where Pw 1, denotes the density of a Polya-Gamma Gamma random variable with parameters (1,
Since our model considers N observations, we apply this augmentation scheme to each data point. Accordingly, we introduce N independent Pólya-Gamma random variables, denoted by (w: [.v,, each distributed according to (w:))) = (w; 1, 0) and with a joint density
Pe (w) = I
Pa (wx)
II PrG (we 1,,0)
(10)
3.2 Poisson Process Augmentation Scheme
Evaluating the likelihood in (3) requires computing N integrals involving sample function drawn from the BNN prior. This integral generally analytically intractable, due to the nonparametric and highly non- non-linear nature of BNN sample paths. To address this, we leverage a Poisson process-based based

---

data augmentation scheme, drawing inspiration from methodologies proposed in [12, 52]. By substituting the sigmoid identity from (9), the intractable integral for the 7th data point becomes
5 L Ao(t, (g((t, 0))
(1- (t, x;; p)pra 0ddddt (11)
where (1, 0) is the density of a Pólya- GaaaaG random variable. The key insight here is that this double integral can be expressed as an expectation over a marked Poisson process.
Before proceeding further, we briefly review the concept of marked Poisson process. A marked Poisson process extends the standard Poisson process by associating each event (or location) with an additional random variable known accompanied by positive mark w. With this in mind, consider the space [0, ys] x R4 which consists of points (t, w) where [0, yal and R+. We then denote by V marked Poisson process [0, yel x R2 with intensity
mark. In our case, each event occurs at time / and is
(1 wi 6) == Ao x;i 6) pra(w (w] 11,0), (1,w) [0, y/ x R+
(12)
Under suitable assumptions on the BNN g(: 9), Campbell's theorem allows us to express the integral in (1 as
(1 10 1 [- -1- ollxxx))) (t, wi )))wdt
000
exp
Bvg B
II
(13)
(
where Pv, is the path measure of the process V In (13), we take the convention that an empty product equals 1. Equation (13) corresponds the term with the intractable integral on the right-hand side of (3). This representation enables us to avoid time discretization, allowing an exact and efficient evaluation of the integral. Since our model involves N observations, we apply this augmentation scheme to each data point by introducing N independent marked Poisson processes, denoted by y -
3.3 Augmented Likelihood
Leveraging both the -aaya- and the marked Poisson process augmentation schemes, we can reformulate the likelihood given in (3) in a tractable way. With these auxiliary variables, we define the augmented likelihood density for the 7h observation as
(ys, 8, xx, 0= we, V.))
elss, ) (14)
II (t,w),
where f(w, 2) was defined in (8). The following proposition formalizes the data augmentation scheme.
Theorem 3.1 (Data Augmentation). Assume for each i = L...,,N N that the function g(- xxxx;)) ) C((0, yi] x Rm). Let p(ys, 5g x;, Φ. g(:;0) be the likelihood density given in (3). Additionally, let p (yi, x, & 5. 0 ) Wi, V.) be the augmented likelihood density defined (14). Then,
P(y, 8, x. d, g(: 0)) Europp Bu B. Ip (yis 8 x4, g(: 0), wi U.)
The proof Theorem 3.1 postponed to Appendix N.I. Existing augmentation schemes ap- proaches [12, 52, 4] do not offer any theoretical guarantees regarding the validity of the methodology. In contrast, Theorem 3.1 provides the first rigorous theoretical framework that establishes the sound- ness of a method within this class of data augmentation techniques.
5

---

Using the assumption from Section 2.2 that (3). 6,) are i..dd conditional on (x, o, g(= 0)), and given the structure of the data augmentation, we observe that (y/, 6,) are conditionally independent of w and v for all + i. As result, the full-sample augmented likelihood factorizes as simple product:
N
p(D X, x.o
= II plu. 6g x, o, (0.g((;)))))))))))))))))())))( V.).
(15)
Variational Inference in the Augmented Space
In this section, develop novel variational inference algorithm based on this augmentation scheme.
4.1 Variational Mean-Field Approximation
Computing the posterior distribution IP (ó. 0, w. vp X) is analytically intractable because its normalization constant is unavailable in closed form. We consider a variational inference algorithm that aims to find an approximating variational distribution @ 0 w, v that minimizes the KI. divergence from the true posterior distribution.
To make the optimization tractable, we restrict our search to distributions that satisfy the following mean-field factorization:
(66 w V = Qo(d) Qo( x (w) (
Here, we take (o). (e( and Qw (w) to admit densities To (Ø))), go go(() and qw (w) with respect to the Lebesgue measures dø, d0 and dw. The remaining factor @ (V)) is a measure on the space of marked point-process paths, which does not admit density with respect to the Lebesgue measures (see, e.g., a similar discussion GPs [34]).
To handle this within the variational inference framework, we must introduce a reference measure Fw... which plays the role of "Lebesgue- base measure on path space (see Definition E.1 for details). We then assume our variational law Qq is absolutely continuous with respect to Pw... that it admits strictly positive Radon-Nikodym derivative dGv, which satisfies the normalization
al'...
d0 [T)] 1. These conditions ensure that Qo is a valid probability measure on the Ep-P lar ... space of marked point-process paths (see Appendix for further technical details).
This formulation enables us to express the KL divergence between the variational distribution and the true posterior in terms of the ELBO:
DKL (Q 0 w. 1 || (6, 6,w, V D,X)) (ELBO(() + const,
(16)
where the ELBO is defined
(gLBO(g) =
(D de ((: 0), w. v po(#)ro (0)p dFwie (T) ae.
B 0-a0 wmm Bn un@
log
(17)
( ge m (w) (
dG TPo.-
is the Radon-Nykodim derivative of the true conditional law P with respect to Pv... From (16), it follows that minimizing the KL divergence is equivalent to maximizing the ELBO.
dl and where aFs.-
4.2 Local Linearization of the Bayesian Neural Network
crucial insight that the data augmentation strategy transforms the intractable likelihood density in (3) into form that is conditionally Gaussian, shown below:
(vs. 8, x, (: 0), wi V.))
exp 8, [V:. 8,
( x 2) Wi
Σ
gity,
x;: 0)
g(t,, x;;
0)
exp
Wj
2
2
(t,w), ev,

---

This transformation is particularly advantageous when placing a GP prior on g(:: 9), as it induces conjugacy in the model. Conjugacy is crucial for variational inference because enables efficient computation of the BLBO (17), which involves taking expectations over the distribution of O. However, when {(: 0) is BNN, these expectations generally lack closed-form solutions, making exact Bayesian updates intractable. As a result, we seek to approximate g(: 0) in way that retains the expressive power NNs while preserving Gaussian conjugacy to enable tractable inference.
We adopt the local linearization approximation introduced in [22]. This approach approximates the BNN g(: 0) using a first-order Taylor expansion around reference point 0*:
g(t, x; ①) gdt 0) g(t, x; 0*) Jo (1, x) (0 0*),
(18)
where x)}; (g(t,x:0) is the Jacobian of the BNN with respect to the parameters 0. Fol- 20,
lowing [22]. select 0* = MAAP the maximum posteriori (MAP) estimate, which is defined as:
(OMAP; (MAP) arg max p(0, db D,X),
(19)
0.6
where p(0, Φ D) is the posterior density defined in (7). By centering the lincarization at MAPP- we ensure maximal approximation accuracy precisely where Bayesian inference is most sensitive: in the high-probability region of the posterior that dominates both parameter uncertainty quantification and predictive distributions. The procedure used to obtain the MAP estimates of (19) is detailed in Appendix F. Under the assumption of Gaussian prior on the BNN parameters (5), the local linearization induces the GP prior
g GP((,, Re)
with mean function μ and and covariance function K. given by:
p(t, x) == gt, OAPP) + J Joss (t, x) (E [0] M <((t, x), x') == Jem (t, (t', x²)
Incorporating this approximation into our variational framework allows to exploit Gaussian conjugacy for fast, flosed-form updates, while still preserving the flexibility of NNs. Concretely, we take Taylor expansion of the ELBO around git and, by truncating at the lowest order term, obtain the simple approximation
(ELBO(G) th))
Our approach is analogous to the method introduced in [47, Section 3.2], where the authors apply Delta Method Variational Inference by approximating the ELBO around fixed point in parameter space. In contrast, we extend this idea by approximating the ELBO around a reference function شار rather than fixed point.
4.3 Coordinate Ascent Variational Inference
We adopt a Coordinate Ascent Variational Inference (CAVI) approach, allowing us to draw on stan- dard results from variational inference (see, e.g., [5, Chapter 10.1]). In this framework, the optimal variational distributions are derived by maximizing the linearized ELBO, (gLBO(g ")) with each distribution depending the current state of the others. The algorithm proceeds by cyclically updat- ing each variational distribution while keeping the others fixed. This iterative process progressively refines the optimal variational distributions, ultimately leading the best possible approximation of the posterior distribution. A complete derivation of each optimal variational distribution is provided in Appendix G while the complete CAVI algorithm is presented in Appendix H.
At the 4t iteration the optimal variational distributions for the parameters & and 0 are given by
(k) (x)
Gamma (a(A),
(k) (0)
N (1 (R) S (k)
where (a) 3) and (et (k), Σ are given in Appendix G.3 and G.4, respectively. At the htth iteration, the optimal update for the auxiliary parameters w is given by
N
(w) (w)
II
pee (ws 1, x,)(
7

---

where (k) is given in Appendix G.1. Finally, at the h iteration, the optimal variational law Q is the probability measure under which each W, (i = 1..... .......)..).. N) is a marked Poisson process on [0, yel x R, with intensity function 2 Q2(0)) as given in Appendix G.2.
It is important to emphasize that we did not impose a specifie form on the variational distributions; for example, we did not assume (o( to be Gaussian. Instead, we derived our results by minimizing the KL divergence over the full space of distributions. This contrasts with methods that fix a parametric form and use the reparameterization trick with Monte Carlo gradient estimates.
Finally, in Appendix I, we demonstrate that, by exploiting the Woodbury matrix identity, our inference updates require only O (m) time complexity (m is the number of weights in the NN architecture). This linear scaling renders Bayesian framework feasible for contemporary large-scale deep neural architectures, which are well suited to model sighddimensionl data.
Experiments
Details on the experimental setup, including dataset descriptions, benchmark methods specifications and evaluation metrics definitions are provided in Appendix J. Moreover, the implementation details for NeuralSurv provided in Appendix K.
To comprehensively evaluate NeuralSurv, we compare its performance against the following set of
benchmark models: MTLR [51], DeepHit [32], DeepSurv [24], Logistic Hazard [16], CoxTime [29], CoxCC [29], PMF [28], PCHazard [28], BCESurv [30], and DySurv [36), Sumo-Ner [43] and DQS [50]. A detailed overview of these models is provided in Appendix I. and summarized in Table A2. Except for DySurv, which employs an autoencoder framework, we adopt the same NN architecture across all benchmark models and NeuralSurv parameterize the hazard function. For DySurv, we use the original autoencoder architecture specified in its implementation.
We assess discriminative performance using the Antolini's concordance index (C-index) [3], and evaluate model calibration with the inverse probability of censoring weighting (IPCW) integrated Brier score (IBS) [17]. the Distribution Calibration (D-Calibration) [18], and the Kaplan- Calibration (KM Calibration) [9]. The C-index evaluates how well a model performs by measuring the concordance between the rankings of the predicted event times and the true event times. The C-index ranges from 0 to 1, where higher values indicate better discriminative performance; a value of 0.5 corresponds random guessing. Similar to the mean squared error, the Brier score (BS) assesses the accuracy of an estimated survival function at some time t. The IPCW are observation-specific weights that account for censoring in survival data, ensuring that the BS remains unbiased. The IPCW IBS is the integral of the IPCW BS over the observational period. The D-Calibration test bins each individual's 's predicted survival probability at their observed event time into equal-width width bins over and applies x2 test to assess whether those predicted probabilities are uniform across bins. A well- calibrated model should yield non-significant p-value. The KMM-alibration procedure compares the average predicted survival curve with the Kaplan-Meier estimate. For the KM-Calibration the closer the two curves align, the better calibrated the model, where 0 indicates perfect calibration and 1 indicates maximal miscalibration and a random prediction yields 0.25. The C-index and the IPCW IBS metrics are computed using the TorchSurv package [37]. The D-Calibration and KM-Calibration metrics are computed using the SurvivalEVAL package [41].
5.1 Synthetic Data Experiment
In this section, we present experiments conducted on synthetic data. The experimental setup was inspired by [13] and constitutes a broadly applicable evaluation benchmark. We simulate the training sets with increasing sizes N 25, 50, 100 and 150 samples where the event time was drawn from two distributions: (T") = LogNormal(3,, 0,8²) and p1 (T) = LogNormal 1²). Each observation includes a covariate indicating whether the event time is sampled from Pa or Pl. along with three additional noisy covariates generated from a standard normal distribution. The censoring times are drawn from exponential distribution with rate of 0.025 yielding an average censoring rate of 54% across the four synthetic datasets. The test set is generated using the same data- -atteranng process, fixed to 100 observations, and held constant across all experiments.
Figure presents the true survival function alongside the predicted functions from NeuralSurv and the two tor-performing benchmark models, selected based on IPCW IBS, Each panel represents
8

---

50
1,0
Surviva Punecion
0.8-
0.6-
0.4
0.2
Sunchee 0:221
StamoNens 6L:22
Sammehee 01.12
(voCc 0.363
Cooci:: 222
0,3
50
100
0
56
100
50
100
50
100
Time
Figure 1: Comparison of the true survival function (black) with the estimated survival functions from NeuralSurv and the two por-rerforning benchmark models (colored) synthetic data. The time axis truncated at the maximum observed event time in the training data. Each panel represents a different training size. The IPCW IBS score is reported for each method in each panel, with lower values indicating better predictive accuracy. NeuralSury estimates the full posterior over survival functions, and the 90% credible interval is shown as a ribbon around its estimate.
a different training set size. As the number of training samples increases, the predicted survival functions match more closely the true survival function. The results show that NeuralSurv consistently ranks as the best method according to IPCW IBS, and its predictive accuracy improves with larger sample sizes. Beyond its competitive performance, NeuralSurv also provides Bayesian credible intervals, offering uncertainty estimates for survival probabilities, an important feature absent deep learning benchmark models. Notably, these credible intervals appropriately narrow as more data becomes available, demonstrating well-calibrated uncertainty quantification. Corresponding C-index, IPCW IBS score, D-Calibration p-values, and KM-Calibration scores for all methods are reported in Tables A3-A4.
5.2 Real Survival Data Experiments
To comprchensively evaluate NeuralSurv, we conduct experiments on eight real survival datasets: the chemotherapy for colon cancer (COLON), the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC), the Rotterdam and German Breast Cancer Study Group (GBSG), the National Wilm's Tumor Study (NWTCO), the Worcester Heart Attack Study (WHAS), the Study to Understand Prognoses and Preferences for Outcomes and Risks of Treatment (SUPPORT), the Veterans administration Lung Cancer trial (VLC) and the Sac 3 simulation study. Each dataset is subsampled to 125 observations to highlight the advantages of a Bayesian approach in data-scarce regimes. The data is randomly partitioned into five equally sized folds, with each fold serving as a distinct train/test split, comprising 100 training samples and 25 test samples per fold.
Table I presents the C-index and IPCW IBS on the held-out sets for three representative datasets, while results for the remaining datasets, as well as D-Calibration and K.M-Calibration, are shown in Tables A5-A6. Across eight datasets, NeuralSury achieves the best IPCW IBS score on seven, highlighting superior overall calibration compared to benchmarks. It consistently passes D D-calibration, together with Sumo-Ner as the only benchmark achieving this result, while ranks fifth in KM- calibration. This strength arises from its Bayesian framework, which naturally models uncertainty and provides effective regularization data-scarce settings. Beyond calibration, NeuralSurv also demonstrates strong discriminative performance, achieving the best C-index in four datasets and the second best in threc.
An ablation study using larger training set of 250 observations is presented in Tables A7-A8. NeuralSurv continues to outperform benchmark methods under this setting terms of calibration performance demonstrating the robustness of the method to training size. Furthermore, we also include results from traditional survival models, such as the Cox Proportional Hazards model [10], the Weibull Accelerated Failure Time model [8], the Random Survival Forest [23], and the Survival Support Vector Machine [40] in Tables A9-A10. These models often achieve strong performance in data-scarce regimes. However, they are not designed to leverage high-dimensional or complex feature representations, which limits their applicability modern deep learning contexts. Our focus remains evaluating deep survival methods that can scale with data complexity, but include these

---

COLON
METABRIC
Method
C-index f
IPCW IBS
C-index 1
IPCW TBS
C-index 7
IPCW TBS
MTLR [51]
0.562
0.548
0.279
0.602
0.273
DcepHit
0,478
0.511
0.243
0.578
0.309
DuepSurv [24]
0.572
0.326
0.523
0.289
0.618
0.252
Logistic Hazard [16]
0,490
0.321
0.541
0.317
0.618
0.296
Cox [29]
0,578
0.277
0.533
0.307
0.599
0.285
CoxCC [29]
0.584
0.289
0.575
0.257
0.646
0.240
PMF [28]
0,509
0.324
0.440
0.336
0.655
0.250
PCHazard [28]
0.338
0.297
0.541
0.291
0.609
0.249
BCFSary
0.491
0.302
0.616
0.277
0.581
0.273
DySurv [36]
0.488
0.536
0.561
0.465
0.572
0.485
Sumo-Net [43]
0.485
0.241
0.447
0.223
0.476
0.250
DQS [50]
0.635
0.246
0.564
0.261
0.611
0.229
NeuralSury (Ours)
0.584
0.657
Table I:
Performance comparison of deep survival models over five different train/test splits of each dataset. The best results for each metric shown in bold, and the second-best results underlined. f indicates higher is better; + indicates lower is better.
classical baselines for reference and completeness. A prior sensitivity analysis for the parameter é, using priors with double and half the original variance, is presented in Table A11. While the posterior distributions under different priors largely overlapped, their central tendencies occasionally differed, indicating mild to moderate sensitivity to the choice of prior. Incorporating prior information remains important, it helps strike principled balance between model flexibility and regularization.
6 Conclusion
We propose the first fully Bayesian framework for deep survival analysis that models time-varying relationships between covariates and risk. On both synthetic and real-world datasets, in data-scarce regimes, our method consistently achieves better calibration than state-of-the-ee deep survival models and matches or surpasses their discriminative performance. In contrast to previous approaches deep survival analysis, which are either constrained to discrete-time settings [51, 32. 16, 28, 30, 36] or lack the ability to provide Bayesian uncertainty quantification [51, 32, 24, 16, 29, 30, 36], NeuralSurv introduces continuous- time modeling framework that naturally incorporates Bayesian inference, enabling both accurate survival predictions and well-calibrated uncertainty estimates.
Pólya-Gamma Gamma and Poisson data aittaaaugmetaatonn schemes (Section 3.3) have been extensively employed with standard Gaussian process models [12, 52]. Likewise, the local linearization of Bayesian neural networks, which yields a Gaussian process-based approximation, (Section 4.2), is well established technique [22]. To knowledge, this work is the first to integrate these approaches into unified framework that capitalizes on Gaussian process conjugacy. By combining these methods, we contribute novel inference strategy at the intersection of Bayesian deep learning and Gaussian process modeling.
Despite its strengths, NeuralSury relies three key simplifying assumptions. First, assume sigmoidal hazard function, a choice shared by prior work (e.g.. [13, 25]), which may not capture all risk dynamics. Second, mean-field variational inference treats parameters and independent, ignoring posterior correlations. Third, we linearize the network around the MAP estimate to enforce conjugacy. In real-world settings, however, the true posterior be multimodal and strongly correlated, so this local, factorized approximation may overlook secondary modes or misestimate joint uncertainty.
Concerning the computational efficiency of our method, the coordinate-ascent updates scale linearly with network size but still require full-dataset passes each iteration. For very large cohorts, this becomes a bottleneck. Extending the algorithm to use stochastic or mini-batch updates would preserve conjugacy benefits while improving scalability.
believe that NeuralSurv has the potential to make a positive societal impact. For instance, as healthcare data becomes increasingly diverse, there is a growing need for models that can handle multimodal data within time-to-event analyses effectively. NeuralSurv represents an important first step toward accommodating such data within a Bayesian deep learning framework.
10

---

Acknowledgments and Disclosure of Funding
SB acknowledges funding from the MRC Centre for Global Infectious Disease Analysis (reference MR/X020258/1), funded by the UK Medical Research Council (MRC). This UK funded award is car- ried out in the frame of the Global Health EDCTP3 Joint Undertaking. SB acknowledges support from the Novo Nordisk Foundation via The Novo Nordisk Young Investigator Award (NNF200C0059309) acknowledges the Danish National Research Foundation (DNRF160) through the chair grant. SB acknowledges support from The Eric and Wendy Schmidt Fund For Strategic Innovation via the Schmidt Polymath Award (G-22-63345) which also supports AM and MM.
References
[1] Ryan Prescott Adams, lain Murray, and David C. MacKay. Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities. Proceedings of the 26th Annual International Conference on Machine Learning, pages 9-16, 2009.
[2] Virginia Aglietti, Edwin V Bonilla, Theodoros Damoulas, and Sally Cripps. Structured varia- tional inference in continuous cox process models. Advances in Neural Information Processing Systems, 32, 2019.
[3] Laura Antolini. Patrizia Boracchi, and Elia Biganzoli. A time-dependent discrimination index for survival data. Statistics in Medicine, 24(24))3927-3944, 2005.
[4] Ifigeneia Apostolopoulou, Scott Linderman, Kyle Miller, and Artur Dubrawski. Mutually Regressive Point Processes. Advances in Neural Information Processing Systems, 32, 2019.
[5] Christopher Bishop. Pattern Recognition and Machine Learning. Information Science and Statistics. Springer, New York, NY, 2016.
[6] Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra. Weight Uncertainty in Neural Network. Proceedings of the 32nd International Conference on Machine Learning. 37:1613-1622, 2015.
[7] Pierre Bremaud. Point Processes and Queues. Springer Series in Statisties. Springer, 1981.
[8] Kevin J. Carroll. On the use and utility of the Weibull model the analysis of survival data. Controlled Clinical Trials, 2466)7622-701, 2003.
[9] Paidamoyo Chapfuwa, Chenyang Tao, Chunyuan Li, Irfan Khan, Karen J. Chandross, Michael J. Pencina, Lawrence Carin, and Ricardo Henao. Calibration and uncertainty in neural time-to- event modeling. IEEE Transactions on Neural Networks and Learning Systems, 34(4):1666 1666- 1680, 2023.
[10] D. R. Cox. Regression Models and Life-Tables. Journal of the Royal Statistical Society Series B: Statistical Methodology, 34(2):187-202, 1972.
[11] Cameron Davidson-Pilon. lifelines: survival analysis in Python. Journal of Open Source Software, 4(40):1317, 2019. (version 0.30.0).
[12) Christian Donner and Manfred Opper.
Efficient Bayesian Inference of Sigmoidal Gaussian Cox Processes. Journal of Machine Learning Research, 19(67):1-34, 2018.
[13] Tamara Fernandez, Nicolas Rivera, and Yee Whye Teh. Gaussian Processes for Survival Analysis. Advances in Neural Information Processing Systems, 29, 2016.
[14] Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPy- Torch: Blackbox Matrix-Man Matrix Gaussian Process Inference with GPU Acceleration. Advances in Neural Information Processing Systems, 31, 2018.
[15] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin. Bayesian Data Analysis. Chapman Hall/CCC Texts in Statistical Science. Chapman & Hall/CRC, 3 edition, 2013.

---

[16] Michael F. Gensheimer and Balasubramanian Narasimhan, A scalable discrete-ting survival model for neural networks. PeerJ, 7:c6257, 2019.
[17] Erika Graf, Claudia Schmoor, Willi Sauerbrei, and Martin Schumacher. Assessment and comparison of prognostic classification schemes for survival data. Statisticy in Medicine, 18(17-18)2529-5444 1999.
[18] Humza Haider, Bret Hochn, Sarah Davis, and Russell Greiner. Effective ways to build and evaluate individual survival distributions. Journal of Machine Learning Research, 21(1), January 2020.
[19] Frank E Harrell, Kerry L Lee, and Daniel B Mark. Multivariable Prognostic Models: Issues in Developing Models, Evaluating Assumptions and Adequacy, and Measuring and Reducing Errors. Statistics in Medicine, 1514)331-3387, 1996.
[20] Nicholas Higham. Functions of matrices: theory and computation. Cambridge University
Press, 2008
[21] Joseph G Ibrahim, Ming-Hui Chen, and Debajyoti Sinha. Bayesian Survival Analysis. Springer series in statistics. Springer, 2010.
[22] Alexander Immer, Maciej Korzepa, and Matthias Bauer. Improving predictions of bayesian neural nets via local linearization. Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, 130:703-711, 2021.
[23] Hemant Ishwaran, Udaya B. Kogalur, Eugene H. Blackstone, and Michael S. Lauer. Random survival forests. The Annals of Applied Statistics, 2(3), 2008
[24] Jared L. Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval Kluger. DeepSurv: personalized treatment recommender system using a cox proportional hazards deep neural network. BMC Medical Research Methodology, 18(1), 2018.
[25] Minyoung Kim and Vladimir Pavlovic. Variational Inference for Gaussian Process Models for Survival Analysis. UAI, pages 435 2018.
[26] J. F. C. Kingman. Poisson Processes. Oxford Studies in Probability. Clarendon Press, 1992.
[27] Håvard Kvamme. pyсox: Survival analysis with PyTorch. https://pypi.org/prog pycox/, 2024. (version 0.3.0).
[28] Hávard Kvamme and Ørnulf Borgan. Continuous and discrete-ting time survival prediction with neural networks. Lifetime Data Analysis, 27(4):710-736 2021.
[29] Hávard Kvamme, Ørnulf Borgan, and Ida Scheel. Time-to-Event Prediction with Neural Networks and Cox Regression. Journal of Machine Learning Research, 20 20(129):1-30, 2019.
[30] Hävard Kvamme and Ørnulf Borgan. The Brier Score under Administrative Censoring: Prob- lems and Solution. Journal of Machine Learning Research, 24(2):1-26, 2023.
[31] Jerald F Lawless. Statistical models and methods for lifetime data. Wiley Series in Probability and Statistics. John Wiley & Sons, edition, 2002.
[32] Changhee Lee, William Zame, Jinsung Yoon, and Mihaela Van der Schaar. DeepHit: A Deep Learning Approach Survival Analysis With Competing Risks. Proceedings of the AAAI Conference on Artificial Intelligence, 32(1), 2018.
[33] Jialiang Li and Shuangge Ma. Survival Analysis in Medicine and Genetics. Chapman & Hall/CRC Biostatistics Series. Chapman Hall/CRC, 2023.
[34] Alexander G de G Matthews, James Hensman, Richard Turner, and Zoubin Ghahramani. On sparse Variational Methods and the Kullback-Leine Divergence between Stochastic Processes. In Artificial Intelligence and Statistics, pages 231-239, 2016.
[35] J W McPherson. Reliability Physics and Engineering: Time-To-Fallu Modeling. Springer International Publishing, 3 edition, 2019.
12

---

[36] Munib Mesinovic, Peter Watkinson, and Tingting Zhu. DySurv: dynamic deep learning model for survival analysis with conditional variational inference. Journal of the American Medical Informatics Association, page ocae271, 2024.
[37] Mélodie Monod, Peter Krusche, Qian Cao, Berkman Sahiner, Nicholas Petrick, David Ohlssen, and Thibaud Coroller. TorchSurv: A Lightweight Package for Deep Survival Analysis. Journal of Open Source Software, 9(104)7341,, 2024. (version 0.1.4).
[38] James G. Scott Nicholas G. Polson and Jesse Windle. Bayesian Inference for Logistic Mod- els Using Pólya-Gamma Latent Variables. Journal of the American Statistical Association, (085504)1339--11 2013.
[39] Sebastian Pölsterl. scikit-survival A Library for Time-to-Event Analysis Built on Top of scikit-leam. Journal of Machine Learning Research, 21(212):1 2020. (version 0.24.0).
[40] Sebastian Pölsterl, Nassir Navab, and Amin Katouzian. Fast Training of Support Vector Machines for Survival Analysis, page 243-259. Springer International Publishing. 2015.
[41] Shi-ang Qi. Weijie Sun, and Russell Greiner. SurvivalEVAL: A comprehensive open-sourge python package for evaluating individual survival distributions. Proceedings of the AAAI Symposium Series, 2:453-4577 01 2024.
[42] Alejandro Quiroz Flores. Survival Analysis: A New Guide for Social Scientists. Elements in Quantitative and Computational Methods for the Social Sciences. Cambridge University Press, 2022.
[43] David Rindt, Robert Hu, David Steinsaltz, and Dino Sejdinovic. Survival regression with proper scoring rules and monotonic neural networks. Proceedings the 25th International Conference on Artificial Intelligence and Statistics (AISTATS) 2022, Vulencia, Spain, 2022.
[44] Hugh Salimbeni and Marc Deisenroth. Doubly stochastic variational inference for deep gaussian processes, Advances in Neural Information Processing Systems, 30, 2017.
[45] Yee Teh and Vinayak Rao. Gaussian process modulated renewal processes. Advances in Neural Information Processing Systems, 2011.
[46] Terry M Therneau. survival: A package for survival analysis in R. htttps///CAAN... R.-AAN..--rojeccc org/packkag--surivva.. 2024. (version 3.7.0).
[47] Chong Wang and David M. Blei. Variational Inference in Nonconjugate Models. Journal of Machine Learning Research, 14(1):1005-103: 2013.
[48] Ping Wang, Yan Li, and Chandan K. Reddy. Machine Learning for Survival Analysis: Survey. ACM Computing Surveys, 51(6):1-36, 2019,
[49] Simon Wiegrebe, Philipp Kopper, Raphael Sonabend, Bernd Bischl, and Andreas Bender. Deep learning for survival analysis: a review. Artificial Intelligence Review, 57(3), 2024.
[50] Hiroki Yanagisawa. Proper scoring rules for survival analysis. Proceedings of the 40th International Conference on Machine Learning. 2023.
[51] Chun-Nam Yu, Russell Greiner, Hsiu-Chin Lin, and Vickie Baracos. Learning Satient-Speciffe Cancer Survival Distributions as a Sequence of Dependent Regressors. Advances in Neural Information Processing Systems, 24, 2011.
[52] Feng Zhou, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen. Efficient Inference for Nonparametric Hawkes Processes Using Auxiliary Latent Variables. Journal of Machine Learning Research, 21(241):1-31, 2020.
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction clearly state the claims made, including the contributions made the paper and important assumptions and limitations. The claims made match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The conclusion points out all the assumptions made in the paper and reflects on how these assumptions might be violated in practice and what the implications would be. Moreover, the conclusion summarizes the computational efficiency of the proposed algorithms and how they scale with dataset size and Appendix gives more details.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model vell-specifiiattin, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
14

---

3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Theorem 3.1, as well as all the formulas, are numbered and cross-referenced. The assumptions of Theorem 3.1 are stated in the statement and the proof is presented in Appendix N.I. The theorems and Lemmas that the proof relies upon are properly referenced.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: The paper fully discloses all information necessary reproduce the main experimental results. Specifically:
Appendix J details how the benchmark results were obtained, the definitions of the evaluation metrics, and the specific software packages (along with their versions) used for model evaluation and metric computation.
Pseudocode in Algorithms and 2 for NeuralSurv is provided clarify the algorithmic steps. Additional implementation details for NeuralSury are included Appendix K to facilitate replication.
The source code is made publicly available via a GitHub repository, enabling direct access to the implementation used to produce the reported results.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often good way to accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
15

---

(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily
model architecture, the paper should descrihe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (c.g., with an open- open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the partieular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The paper provides open access to both the data and code, along with sufficient instructions to faithfully reproduce the main experimental results:
The source code is made publicly available via GitHub repository, enabling direct access to the implementation used to produce the reported results.
The generation process for the synthetic dataset is described in Section 1. and the resulting data are available for download from the accompanying GitHub repository.
The real-world data used our experiments are publicly available open-source datasets. The source of these datasets, including the specific package and version used to obtain them, detailed in Section J.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cr cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: //nips. http://c/wpccguuieeeeeee/e/////e wwwggguiiss//e/ httpg/uuiddee/eeeeeeaeeeeeeeeeee Cww.deauaiiiiiiioooooioiioooeoee
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question:
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
The paper specifies all training and test details necessary to understand the results:
Justification:
16

---

The train-test data split for the synthetic experiment is described in Section 5.1. The train-tes data split for the real data experiment described in Section 5.2 and further elaborated in Appendix J.1.
The model architecture is detailed in Appendix K.
Hyperparameters for the benchmark models, including how they were selected, well as the optimizer, batch size, and number of training epochs for the benchmark methods, are provided in Appendix J.2.
For the NeuralSurv model, the prior hy hyperparameters and the value of the regularization parameter p specified in Appendix K.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [NA]
Justification: The paper does not report statistical significance tests. It focuses on settings where statistical testing not applicable for supporting the main claims.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative eтror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The information on the computer resources is provided in Appendix K.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
17

---

The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https:/ /ttp:::///nurips... cc
Answer: [Yes]
Justification: To the best of our knowledge and judgment, the research complies fully with the NeurlPS Code of Ethics. our view, the work does not pose potential harms caused by the research process, does it present any adverse societal impact potentially harmful consequences.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g., if there is special consid- cration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: The potential positive societal impacts of the paper are summarized in the conclusion. No negative societal impacts have been identified.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairesss considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g.. gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
18

---

Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model. for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
The benchmark methods used in the paper are properly credited, with citations to the original papers (where applicable) and detailed references to the software packages and their versions Appendix J.2. For the datasets, any original studies associated with the data cited, and the process for obtaining the data, including the relevant package names and versions, is explained in Appendix J.1. All licenses and terms of use have been respected.
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The source code is made publicly available via GitHub repository under the MIT License, enabling direct access to the implementation used to produce the reported results. This is explained in Appendix K.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/code/mmodel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
19

---

At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing. editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non standard components.
Please refer to our LLM policy (https / ///ereurip...// cttpps.cc/////o fcc//oneerencee/222////////////e for what should or should not be described.
20

---

A Review of Survival Analysis
This appendix offers concise summary of the survival-analysis framework on which our approach is built. For an in-depth review, the reader is referred to [31].
Survival data for each observation consist of three components:
Feature vector: A covariate vector x E I' capturing baseline characteristics;
Event time: nonnegative random variable T measuring the time from baseline to the occurrence of the event of interest;
Event indicator: A binary variable 8, which takes the value 1 if the event is observed, and if the event is not observed within the observational period. In the latter case, the observation's 's data is said to be right-censored, meaning that the only available information is the time of the last follow-up before the event could occur.
To handle censoring uniformly, introduce censoring time C' and record the observed time y = min(T, C). The event indicator can then be written succinetly as 1 [Cg Throughout, assume noninformative right-censoring. i.c. conditional the covariates, the censoring time is independent of the event time: C' 1 x. Conditional on x, we let the event time T have cumulative distribution function F(] x) and probability density function f(()x) such that
S.'
x) = P(T s 1
f(s x)
fort E [0, oo). The survival function gives the probability of remaining event-free beyond time
S(t x) = PT PT= F(t x) f.a f(s x)ds,
for t [0, An important modeling quantity is the hazard function, which represents the instantaneous event rate at time given survival up to t:
P(t -
< t At T txx)
f(t x)
A(t x) =
lim Ar=0
At
5(t x)
Equivalently,
d log S(t x),
X(t x)
so that the survival function can be written in terms of the hazard:
S(t x) =xp A(s x) ds)
21

---

B Review of Pólya-Gamma Random Variables
We follow [38] in defining the family of Pólya-Gamma distributions and their properties.
Definition B.1 (Pólya- Gamma Distribution), A random variable is said to follow a Polya- distribution with parameters > 0 and e E IR, denoted by w ~ PG(b, c),
9k
i.d. Gamma(b, 1).
with gk
(A1)
o²
4m
The following result expresses the reciprocal of the hyperbolic cosine function raised to the power b as an infinite Gaussian mixture. This representation central to connecting the Pólya-Gamma density with parameter с + to the case when c= O.
Proposition B.2. The reciprocal of the hyperbolic cosine raised to the power b can be represented as an infinite Gaussian mixture:
cosh 9
exp (-Zw) Pro(w 0) dw.
(-Ew)
Notice that Proposition B.2 can also be read as providing closed-form expression for the expectation (u-=po([=(b_0) e Building on this representation, we can relate the density function of a Pólya-Gamma Gamma random variable with non-zero parameter through exponential tilting of the Galya-Gamma random density with = 0. This connection is summarized in the next proposition.
Proposition B.3, The Pólya-Gamma density (A1) can re-wrilten the form
c² (cosh(c/2)) PrG(w 0).
Pro(w b,c) exp
(A2)
2
The previous propositions not only establish key representations of the also facilitate the derivation of its moment properties. In particular, one can derive the moment generating function, from which the first moment follows directly. This is captured in the next result. Proposition B.4. Let ppo(w] b, c) denote the density function of the random variable w ~ PG(b, c). with > 0 and E R. Using Propositions B.2 and B.3, the moment generating function is given by
density but
10 &" @rol( b, dw
cosh (osh((c/22 (c/2)
(A3)
cosh ( (e 26
In particular, the first moment is obtained by differentiating (A3) with respect to & at 0:
Bwmpnotwlbooouu tanh
2
(A4)
2c
Finally, the following theorem illustrates how the Pólya Gamma distribution can be used derive useful integral identities.
Theorem B.5. Let pro(w b, 0) denote the density function of the random variable w PG 0)
with b > 0. Then, for all E IR, the following integral identity holds:
(1 e =2 be 1
e tha
w
exp
PPG(w b, 0) dw,
where 6 =
The following corollary is direct application of Theorem B.5, Corollary B.6. Let f(w, 2) =
2 12)(22. Then,
e ( ((w..))ppole 1, dww..
(2)
(A5)
(
22

---

C Review of Poisson Processes
This appendix briefly summarizes the properties of Poisson process that most relevant to our analysis. For more comprehensive treatment, see Chapters 3 and 5 of [26].
Definition C.1 (Poisson Process). Let 2 be measurable space. A random countable subset
= (a 2}
said to be Poisson process 2 if it satisfies the following properties:
1. Independence: For any sequence of disjoint subsets {2x C ) the counts
((x))= = n 2xl
are mutually independent.
2. Poisson Counts: For each measurable subset Z₁ 2, the count (2z)) is Poisson distributed with mean
/
2x X(=) dz,
where A 2 - R is the intensity function.
Given point process V, denote its path measure that is, the probability measure induced on its sample-path space - by Pu. If the intensity function A(z) is constant, A(z) = X, then Vi is called homogeneous; otherwise, it is inhomogeneous. We now extend the concept of Poisson process by incorporating additional random attributes, known as marks.
Definition C.2 (Marked Poisson Process). Let = {z E 2} be Poisson process on 2 with intensity function À 2 R2.. Suppose that for each point z, associate random variable w, such that w ~ Puls (w|z). taking values in some space M. Then the collection
M = {(z, E Z x M}
defines a Poisson process on the product space 2 x M. The resulting process is known as a marked Poisson process with intensity
A(z,w) = A(2) (w]=)
Next, we present Campbell's Theorem, which describes the law of sums taken over the points of Poisson process (see [26, Sec. 3.2]).
Theorem C.3 (Campbell's Theorem). Let M be a marked Poisson process on M with intensity function X(2, w) and let f Z x M IR be measurable. Then the sum
H(WM) =
Σ
f(eg,wj)
(M.W))EW
absolutely convergent with probability one if and only if
Jax (2)([[((,,)),,)))33 widzdl < co.
this condition holds, then
ENM [eeH(VN)) exp (
EXM (e<s(ow) - ())zzw))dddd
for any E for which the integral the right converges. Moreover
Ey EuurPay [H)MM)]= = -
f(2,W))(,,
in the sense that the expectation exists if and only if the integral converges.
23

---

D Obtaining the Normalization Factor [tt, x)
In this appendix we derive an efficient approximation for the normalization factor
[1, x) = EB-pe (1,
x; x))]
(A6)
which is needed when computing the CAVI optimal updates (see Appendix G), Recall from (5) that has following prior distribution
(0, Lm.:
where I. is the identity matrix. Moreover, recall from Section 4.2 that approximate the network output (t, x: 0) around some reference 0 by its first-order linearization
(t, x: 0) g(t, 0*) + Jo (t, x) (0 -0"),
where J (1, x) denotes the Jacobian of g(, x: 0) with respect to 0. Because is Gaussian, the linearized output is also Gaussian:
Jm (1, x: 0) N (9(1,x;07) - -Jo (1, xx0 0", | (t, x)]
In order to approximate xt,,)) we wish to leverage well-known asymptotic approximation. Specifically, for normal random variable X NN it holds that
[x-N(a,>)0((X]] a
(
(A7)
VI+
We can apply the result in (A7) to the normal random variable g (t, x; 0) and approximate Z(t, x) as
g(t, x: 0-) Jo- (t. x) T
0*
X(t,x) =
V1+ | (t. x
Since in (A6) we are taking the expectation under the prior pa (0), it is natural to linearize around the prior mean, therefore, we 0 0.
24

---

E Combining Variational Inference with Poisson Processes
In this appendix, outline how our Iriiational-inferemn framework integrates marked Poisson processes an essential part in the mean-field variational approximation of Section 4.1. For a fully rigorous, measure-theorenn treatment, the reader referred Brémaud's text [7]. Our development relies in particular on Theorem T10 in Chapter VIIII of that book, which shows how the law of a marked Poisson process arises via change of measure using the appropriate Radon-Nikodym derivative.
We begin by fixing reference measure path space:
Definition E.1 (Reference measure Pv,:). Let T ( 1..... WN) be N independent marked
Poisson processes, where each V is defined on the product space [0, yol x R.,. We define Pv,- to be their joint law where each V, has intensity
(t, w) = t0- "ro(w 1, 0) for all (t,w) E [0, vsl X R+.
(A8)
Next, let 7/ (t) be deterministic function [0, yil and let n (f,ww) be deterministic density [0, wi] x R4, satisfying
th h (t, w proa 0)dw = and Jt=, (1) <0
(A9)
for all t E [0, 3i] and = ,,...,N. .N. It convenient to introduce the function
(1,w) == 1 (1)) for all (1,w) E [0. Rt,
as well as the functional
N II
o (())))())))))))) exp (1%/ (Nox(l.w)) X3 (1,, dwall).
II
==1
(t.))Ev;
By Theorem T10.b [7, Chapter VIII], whenever Ev-PP,. [L(V]] = 1, the measure Qu (V) defined by dGw ( L( is exactly the law under which each V is a marked Poisson process on dP
[0. yil x R, with intensity (t, w). The above result underpins the analysis in Appendix G.2, where we show that the optimal variational measure Q coincides with the law of collection of independent marked Poisson processes.
Finally, the measure P also admits a Radon-Nykoding derivative with respect to Fs,, which is given by:
dF の dl v.-
exp (f% ()) Aell, u)) dwil)
V II i=1
II (t,)),VV;
Z(1,, xs)
(A10)
Ae(tywoyid)) )
Notice that xx)
i.c. the ratio of the intensities of Pvis and Pv...
25

---

F
Obtaining the Maximum a Posteriori MAP
We seek the maximum posteriori (MAP) estimates
(fMAP, (MAP) arg max (0, X).
Applying Bayes' rule gives the following expression for the posterior density
logp(0, X) log X. g(: Φ)
log po + ps())
where the likelihood density p(D X. () o) and the prior densities po( and p specified in Equations (4), (5), and (6), respectively. Since the log likelihood distribution is intractable, direct optimization of the posterior distribution is infeasible.
F1 Approximating the Log Likelihood distribution
Variational Mean-Field Approximation. Our aim is approximate the log likelihood density
log X, g(: 0). Φ). order do so, we introduce variational distribution @(w, 0. Φ) to approximate the true distribution PP(w, V DD X, ((: 0), ②). Such variational distribution differs from the one used for full-model inference Section 4.3 because it is conditioned on the values and Φb. Hence, we adopt the notation Q (instead of Q) to highlight this difference. We restrict our search to distributions that satisfy the following mean-field factorization:
(ww, y 0, Φ) Qul 0, Φ) x Quio, V 0)
Here, we take (w) 0, o) to admit the density qw ( 0, () with respect to the Lebesgue measure dw.
For the marked point process component, we assume that the variational law Q e is absolutely continuous with respect to ...... so that it admits strictly positive Radon-Nikodyn derivative
ddpi0.0 which satisfies the normalization Ey-.o dFe.-
Td
dP v. (T) 1. These two conditions
guarantee that @ a is indeed a probability measure on the space of marked point-process paths.
We decompose the log-lite ligellkelioodd as follows:
log X, (),()()))))
DRL. Qw. ve 10, V D. g(=,)((;)),)) + LMBO (A11)
where the ELBO is given by:
p(D X, ((: o.w, w 1 (w) dFw (T) aP..
İELBO E Eu- 2010.
log
(A12)
(w 0, devie. V 0, +b)
and where di is the Radon-Nykodim derivative of the true conditional law Peej with respect to aFo.. Pw... cf. (A10).
Minimizing the KL Divergence. When the variational distribution (ww, vp 0, σ) matches the true posterior P(w, V (D, g(; 0), o). the KL divergence term in (A11) vanishes. Consequently, the ELBO becomes equal to the marginal log-likelihood, and maximizing the ELBO is equivalent to maximizing log-likelihood directly. In practice, we minimize the KL divergence so that our ELBO provides the closest possible lower bound the true log-likelihood. Therefore, in order to obtain the closest lower bound to to the log-likelihook log P(D X, g(: 0), Φ) we must find the distribution (ww, V 0, Φ) which minimizes the KI. divergence in (A11).
Using standard meen-field variational inference techniques (see, e.g., Chapter 10.1 of [5]). the optimal distribution for the latent variables w given (0()) (()) is obtained by computing the expectation of the joint lon-density with respect to the other variational factors, that is
log (w))
dP (o(s)
E S Co
logr(D) X, (x,((;0)), ((0),, w, V + log (w)) + log
(T)
+ const.
dP v..
26

---

A similar update applies for y given (0 (0)),
dQ
log
dP dPv..
cP (o!!!! dlP v.+
logp(D) X, 0;; (()))()))))))))))))))))))))))) ,00 v + log Pw Pw (w) + log
(w)
+ const.
Following the same derivation as in Appendix G.I, we find the optimal variational distribution of w given (0, ()
N
N
lojo scc (w ((0) (0()) = II
Uae 0i ,000 ,000
(w, ( (00) = II
PPG (w) 1, 1,
(A13)
where
(1) 6; (9(s:: x;; (()))
(A14)
By mirroring the derivation Appendix G.2, one shows that the optimal measure Q ( 4 ()) exactly the law under which each F, for I,..., N, is marked Poisson process on [0, y/ R. with intensity
(()) " PrG (w 1, [; X (t,w 000) (00) = 2
where we set
1P-1 7(1,x)
gt, x;; ((0) + (g x,; ((0)
( (
(0))))
1 (lg(1, x;; (0)) exp
(A15)
F.2
EM Algorithm for MAP Estimation
Maximization (EM) algorithm provides an efficient framework iteratively maximize the Q-function, At each iteration 0, 1, 2.... we perform the following three steps:
EM Algorithm.
The ExpectatoonMaxiimiztionn
1. Latent Variables Update. Given the current estimates (0, ()) update Twjmnn 1ere (w gl0) (l0)) and Quioin ,stes (v ((0), e(0)) according to (A13) and (A15), respectively.
2. E-Step. Given current estimates (0, (9) compute the Q-function:
(0)))))
log (P(D) X, w, vt) Pw (w)
dIF 2
E B x 00000 201)
+ log po (0) log (o))
(A16)
Note that the entropy term of the ELBO (i.e., the denominator) is not included as it does not depend on the parameters (0, σ) but on the current estimates (0, Φ)) (+), hence it is inrelevant the parameters' optimization.
3. M-Step. Update the parameters by maximizing the Q-function:
(0, ((+1) arg max Q( Φ) (0,
0,0
Steps 1-3 repeated until a given convergence criterion is met. We provide algorithmic description of our EM algorithm in Algorithm 1.
27

---

Algorithm Expectation-Maximizaatee (EM) for maximum a posteriori (MAP) Estimation
1: Initialize: Set initial value for (g(t) ())).
2: Set: iteration counter l <-
3: repeat
4: f 2++1
5: Latent Variables Update:
6: Update gl.)
N
7: Update:
given ( following (A14).
i=1
8: Update
(
9: Update:
(.( ())
i=1 given ( and (((+) following (A22).
, 10: E-step: Evaluate the Q-function Q((0, Φ) (0, ()) given f0 ()) N
} o² and o
following (A16)
11:
M-step: Update parameters by
(0, (0,0) arg max Q((0,)) (0, (0)
0.6
12:
until Convergence criterion is met 13: return (0!t) ()(1)
28

---

Computing the Q-function. The optimal distributions which minimize the KL divergence can now be plugged in the ELBO of (A12) obtain the closest lower bound the log- oogliikeihood. We now recast the MAP optimization problem in term of this lower bound. Specifically, define the following Q-function
(0)(0)))))))
= E w 000) alo) alo)
w. V) [log X, g(: 0), w.
c /PPv dP v,-
+ R Bund woll)) alo) p- 1000
log(pw(w) + log
(V)
+ Po (0) log (1)).
We now wish to derive closed-form expression for the Q-function which can be used in the MAP optimization. Specifically, using the augmented likelihood factorization in (15), obtain
N
Q((0,))(( ()
E
[log x.g(::)) V.]
wi
10000000000
0. (0())
dIP 1i
+ +Ewnn a(e) log Pw (w)) + E
log
Hog Pa (00) const
const.
d'
Next, by substituting the expression for the augmented likelihood in (14), for the priors (0( in (5) and (d) in (2) and for the Radon-Nikood derivative of Paja with respect to Pv,, from (A10),
we obtain
(0. 2m1
6, logo + {1::x;;0) I(Vi, x;; 0)2
=>6
E Burn 10000 .010) [wsl]
+ E Wy B
Σ
f(wy, o(tg/.xx: xxx 01)
S/a 9
(tww),
(t,w)EE A(1y xxi )
+ Byg w ,000)
log
1 oTo 0 + log(0) 1) 33 + const.
2
We apply Camphell's theorem (see Theorem C.3), we substitute the expression for the baseline hazard A1(: σ) from (12) and we substitute the expectation using the optimal variational distribution of w; from (A13), to obtain
N
glyi, x;;
g(vi, x;; 0)2
Q( (0. ()))
Σ
tanh
2
40
2
i=1
g(,, Xs; 0) 0 (.0))
gtt, x;; 0)2
lo(t,x;;0(0)
1 (
o(0))t
tanh
lg(t, x;; (())
N
dx + Io" 5
TTo + log(()
Σ
29(1 000) ol0))dt -1 -1
ar +
1=1
N
1h
1P-1
o 30
+
dt
+ const.,
Jo
Z(t, x:)
i=1
where 2 (t 0(0) ((00 is shown in (A15).
29

---

G Coordinate Ascent Variational Inference Optimal Updates
In this Appendix present heuristic derivation of the CAVI optimal updates presented in Sec- tion 4.3. Before presenting the next results, we define here for convenience
(k) m (t) = E a~ (2) [g ndt (t,x;;0) 0)]
(t) = V E 0~ (t,
for k 2 O.
G.1
Optimal Update for w
Using standard mean-field variational inference techniques (see, e.g., Chapter 10.1 of [5]), the optimal update for the latent variables w obtained by computing the expectation of the joint log-density with respect to the other variational factors. In particular, have
log (
(w) E (
(1-1) 00 1)-1) 0 log (D Φ, (: 0), w. w) + log 1o (w) + const.
Using the augmented likelihood factorization in (15), the expression decomposes as
log (1 (w) Σ Bang 10-1) 0~40 (0-1) N1-06 "
(1-11 log p (yi, 8a xi, "i, (:;0,, W:)
+ log (w) const.
Next, by substituting the expression for the prior Pw from (10) and the augmented likelihood from (14), we obtain
log de (
(w)
Σ
wio:
(5( (s/t-1 (vx))
log PPG (wi(1,))) const.
2
Finally, by applying the identity in (A2), we deduce that the optimal variational distribution factorizes as
(w)
(w,)
II PrG (ws 11 11
where
al 8,
8
(yA)
(A17)
Optimal Variational Expectations for w. for updating the other variational factors with
From Proposition B.4, we obtain the required expectation
(k)
E
[w)]
tanh
(A A11
22
for i = 1,..., N. Notably, since this expectation is always multiplied by 6; when updating other variational factors, it remains well-defined in all cases.
G.2 Optimal Update for T
Using standard mean-field variational inference techniques (see, e.g., Chapter 10.1 of [5]), we obtain the optimal Radon-Nykodim derivative aPe.
d0. by taking the expectation of the joint log-density with
respect to the other variational factors. In particular, we have
dQ
(R) d
log dIP vp
(
E g @m (0 (k-1) w [logp(D) Φ. (:: 0), w. 4)
log dP P 4T)
E A~g (K-1)
+ const., (A19)
dIP di w.
30

---

where the constant tern absorbs all terms irrelevant to the optimisation. Using the augmented likelihood factorization (15), the expression in (A19) decomposes
(R) dQ
log dlF v.
E o (1-1) 1), @
(1- wir
(8) [log di 0), Wi, V:)]
wi,
dP l
E S
(2-1)
log
+ const.
clP
Next, by substituting the augmented likelihood from (14) and the Radon-Nikodym derivative of Puré with respect to Pv,, from (A10), we arrive at the unnormalised form
dQ dIP +p
N
log
B B [f(w), -9 lin (t,, x;;
i=1 (L,wo)y EVV,
Z(tz. x:) )
)
E G (
1k-1)
log
const. (A20)
(+,))EV;
Plugging in the definition f(:,)) from (8) simplifies (A20) to
dQ
(k)
N
m (t,)
(t,))² wj + log(2)
log
(1)
Σ
dlP 4
2
2
(L,wel EVe,
N
ding 1- flog Φ) - llog 2(4 x:) + const.
1=1 (tW
To express this closed form, define for each i 1,
N and (t,w) E [0. R, the functions
[k (1)
8 (t)
h, (200k (t,w) = exp
cosh
2
2
24 (t)
1
o
exp
mi
(1-1)
(t)
8
(k-
(t)
+ E &
(- (log ol
Z(t, xx)
2 20.00 (t,w)= = (k) ht)
2.00
(t, (t,w),,,
where (t, is the intensity defined in (A8). Furthermore, we define for convenience, (k)kk (t) to t0-1
G.(k) ( (t). (t).
(A21)
Notice that by using expression (A2), the function X (k(k)) (t, w) can be written as
e. (t.,)) = X G.(k) (t) PPG (w 1, 8 ()((
(A22)
Finally, enforeing the normalisation condition
do dP v.
Exp-Pa..
V
=1
together with Campbell's theorem (Theorem C.3) yields the normalized derivative
dQ
(V)
dPP v..
2,000 (t,) (t.000))) exp (foh) ( (drslttw) - -20 x2.10)(t,,)) (t,w)) dwat))
II
II
G.(K)
74
i=1
(t.)),
31

---

7i (t,)/; Notice that the products Q.(k) (.( (t) are all strictly positive 2, hence dēg.. is also strictly
dQla
positive. Under suitable regularity conditions on g. one can show that ( (l, w) and 7; (1) satisfy the integrability criteria of (A9), so that 266 is the probability measure under which each V. (i= 1,....N) a marked Poisson Process on [0. 5/) x R with intensity function X ( (1, w).
Optimal Variational Expectations for V. From Proposition B.4, we obtain the required integrals for updating the other variational factors
Jar 0.1k) (l,w)dw 2, (1) (L);
s(a-1((1)
(k-
I
(1, = =1,- (2.(4)(1)
(1) (k
tanh 8
28:
2
G.3 Optimal Update for
Using standard mean-field variational inference techniques (see, e.g., Chapter 10.1 of [5]), the optimal variational factor for the parameter @ obtained by computing the expectation of the joint log-density with respect to the other variational factors. In particular, we have
dF (T)
(k) (d) E. Bo (6- - ( v logp (D "h (::00) V log log 9
dF +,+
+ log (Ø)) const.
Using the augmented likelihood factorization in (15), the expression decomposes as
N
(k) log 9
(0)
E Bmgl (k-1) wigggh) o
log p (y:, 6/ o, O, wi, U.))
dlP v
E
log
v
+ log (o( + const.
dP dP
Next, by substituting the expression for the augmented likelihood from (14), the Radon-Nikodym derivative of IP ws with respect Pv.. from (A10), and the prior of Φ from (2), we obtain,
n 1 xt;(
(k) (0) log 96
6, log Au(u,,
))
E B
log Z(tz,x))
+ (av 1) log(d) Buo const.
(L.))y
We apply Camphell's Theorem (Theorem C.3) and substitute the expression for the baseline hazard )- from (2), to obtain
log
(R) (Ø)
log((Ø) 00
Σ
+
f Vi x (a))))
-
Bo
+
Σ i=1
(
10-1 x11,x,)
dl
const.,
where A (.(k) (t) shown in (A21). We deduce that
(Ø) = Gamma( (k),
where with shape (A) and rate 3 given by
(A)) ao Σ 0, [(1++)
X (t)at)
3 9o +
+
21 10-1
dt. Z(t,x,,)
(A23)
"ee Lemma N.1 for proof of the strict positivity of the normalization factor [tt,
32

---

Optimal Variational Expectation for Φ. We obtain the required expectation for updating the other variational factors with
(4) Dog Φb) = (a(t)) -log (
E
(A24)
where (-) is the digamma function.
G.4
Optimal Update for
Using standard mean-field variational inference techniques (see, e.g. Chapter 10.1 of [5]), the optimal variational factor for the parameters 0 is obtained by computing the expectation of the joint log-density with respect the other variational factors. In particular, we have
log do
(0)
E Lr (2) w 12 Q 14
log (D) ( (=: 0), v) + lo (0) + const.
Using the augmented likelihood factorization in (15), we obtain
(k) (0) E d (k) wirds (2) v. [logr (yi, 6i xi, 0, (=;0) ws, v:)]
log do
+ log (0) + const.
Next, by substituting the expression for the augmented likelihood (14) and for the prior for 0 from (5), we obtain,
N
8
log 9o
(1) (0)
2 9
lin
(yt,
X0; 0) EE wing (ws] I "h (2,, 0)2)
E W
Σ g
0 TT 0 const.
(BW),ER;
We apply Camphell's Theorem (Theorem C.3) obtain,
N
log 90
(k)
(0) =
Σ
8: lin (Bii X;; E w [eslg (2): 0)²
9
(
2
i=1
+ 2 2.
/
(f, (-glin
0) g (t; X (k) (t, w))dd
0 o + const.,
where A (G,0k) (t) is shown in Equation (A22). Next, we recall the expression for (:: 0) from (18)
and we notice that
g (:: 0)
Jo (.) + const.
g Jin 0)² 0" Jo (.) (2g(: 0MAP) 2J exee (-) ( + TJJeww
Josar (-)Jeme (-) T const.,
where the constant term represents terms that do not depend on We substitute the expression for Lin (:; and (:; 0)2 and we obtain,
log dp (0)
(k)(0) = OT ATAU (K) ( + const.,
where
N
A (R)
o,J
j.Jorno (yi, x:) (1- 2E,
mangl (1) [wi] (g(ye, x;; 0MAP) -Jnn (U:, x;) T (MAAP)
2
1=1
(k)
(k) I((2) (MAPP)
1,
+
2
(A25)
N Σ
((k)
2
1
8;
E
wimg
(2)
[wi]J (Vi. X, x;)Jor Jove Xi) T + I + Im
(A26)
1=1
33

---

and
I
J Buw (t, x;) (.(kk) (t)dt (t)
Je Vuar (t, xt)g(t, x,; MMAP) ) (0.00) tanh 8 ( (1)/2)
(()/2())/2)
(K)
2,
1t) dt
2 Io
(k)
U) J Je (t, (t.x;))Bsauu Brea x,) X (2((k)
tanh
(1)/22
dt.
28; kk 11 (t)
A. Σ1,: and 2,, vectors of the same length of 0. B and Z3 are square matrices for which each dimension is the length of 0, and Im is the identity matrix of length of O. We deduce that
(2) 40 (0) =N (it D (M))
where
(b)
(
-1
A (k)
S = 2 (B()
(A27)
Optimal Variational Expectation for 0. We obtain the required expectation for updating the other variational factors,
(6) [0 Td x;; 0)] g(l, x;; 0MAP) Jome (1,, x;) (iilh) (MAP))
E 0~
(g 1n xx; 0)² (g(1, x;; 0MAP) Jowar (1, xx) (nelh) 0MAA)) ))
(A28)
+ J Jonr (1, x;) x;) Ts)
Jouso (1, x:).
34

---

H
Coordinate Ascent Variational Inference Algorithm
Algorithm 2 Coordinate Ascent Variational Inference (CAVI)
1: Compute: Compute 3 following (A23).
2: Initialize:: Set initial values for ( and (i, (0)
(rruge 3: Compute: E (o) [log Øl given (a 3) following (A24).
4: Compute: (m,, (0))) j given ii, 2 (0)
following (A28).
5: Set: iteration counter l +-
6: repeat
7:
A + k +
8: Update two
9: Update:
c i-- given (3:(-)(--1))
) | following (A17).
{ R
mey wil] 1 given (k)
10:
Compute:
J=1 following (A18).
Update (
11:
(A8:00())) 1 given (((m:((, E Bt log [lg ob) follow-
12:
Update:
ing (A22).
13:
Update 4x
(R),
Update: (k) given { {2G: (x((10))))
()) 1= following (A23).
14:
15:
Compute: E gung go [log Φ) given (x) B)) following (A24).
Update (R)
16:
Update: (p, 5 (1)
given بيبي [ 2 (.)00)
17:
following (A27).
im1
(k)
18:
Compute: {-) (+))
i== given
(i, 5
following (A28).
19:
until Convergence criterion met
Return: Optimized variational distributions
90 (k' (k))(0)
( ())) and o (Φ)
Gamma (a(h), 3,). where is the final iteration after convergence.
35

---

I
if k K;, k >> K, Computational Speed-Ups
Survival-analysis cohorts often comprise only few hundred to few thousand observations, yet modern deep learning models may involve millions of parameters, putting us in the N < m regime. complementary strategies that avoid any expensive m- dimensional inversions or factorizations by leveraging the fact that the nontrivial part of our key matrix is low-rank relative to the full parameter dimension m. We also show how heavy censoring further reduces the computational burden.
To exploit this disparity, we develop
To streamline what follows, let introduce the shorthand
@MAP (3/s+ x,) RR=x
for I....,,A N. With this notation (and dropping the CAVI-iteration index for clarity), the matrix B E IR defined in (A25) becomes
N
B
6, Euy [. + Is.s) + I
2
Here, each Iz,; is the integral
Ia 1 Jomae (1, xs) Jona x.) (1) tanh (5:(t)/2) di 28;(1)
and in general admits no closed-form solution. We therefore approximate by any standard quadrature rule (e.g. trapezoid, Simpson's ss or Gauss-Leggndre)) In what follows, we will illustrate the argument with the trapezoid rule, though the same steps apply to any other quadrature method.
We begin by introducing uniform grid of points along the time axis:
t1 11.t2....tt. tk,
where 1 and tx max(v4 N We associate set of quadrature weights {ua | to the time grid points, tailored for each observation These weights correspond to the trapezoidal rule for numerical integration the interval [0, yıl, and are defined
if 1: = and t Vi: Ki and la < yis
1
46-1
if
Usk
tK
2
0,
where K, = max(k E (1...., . K} : ta < yi Further we denote by V, the collection of quadrature weights for observation such that
V = (va...., v.k) E RRK
We collect the Jacobian evaluations into the matrices
Q [JoMAP (41, (;) Jo (12,x;) x:)
JOMAP (IK, x:)] e R=*k
With these definitions in hand, any K K-point quadrature rule yields the approximation
Ia R3
Uik Jomar (tx, xi) Jo (1x, x;) Q. Q.
1=1
Likewise, each term
B, [
can be written in the form J. C.C, where the scalar C; 8, Bu, [ws].
We collect all contributions into single matrix U E mR R where R = N + NK. This matrix is constructed by horizontally concatenating the vectors J, and Q; for i 1...., N, follows:
U = Ji J2.
JN. Q1 .202..... .2....,Q.. Q2
(mx1)
(mxX)
36

---

Further, we define the block diagona weight matrix
C := diag 6 (Eww man [wi].... 8N Ow
[NN V
Vv) ERRRR
(N)
(K)
It is straightforward to verify that
B = 2
(Im + UCUT).
Applying the Woodbury identity (see [20, Appendix B.10J) then reduces the inversion of B to that of an R. R matrix:
B B- = 2( + -
Forming the Gram matrix U requires O(mR²) operations (each of its R² entries an inner product of two length-m vectors) while inverting the resulting dense R x It matrix costs O Therefore, assembling and solving the small system costs
OmmR²) + O(R") O(mR² + +R)
instead of (m)²) for full m X m inversion. Whenever R < this yields dramatic speed-up. By replacing the direct (RR²) factorization with Conjugate-Gradeen Gradient (CG) solver is commonly done in Gaussian Gaussian-process toolkits such as GPyTorch [14] we reduce the cost to O(R²)
Finally, many survival datasets exhibit censoring, i.e. 0 for fraction of observations. Since censored observations contribute only through the integral term, we may further partition the low- rank factor U into blocks for uncensored and censored cases. The effective rank becomes RL' Nuncensored + N where Nuncensorod is the number of uncensored observations, so that any Cholesky or CG solve scales with (Nuuncensored + K) rather than (N -N).. When N Nunceusored < N, this yields an additional, potentially large reduction in computational cost.
37

---

Experiment Set-Up
J.1 Real Survival Data
The real survival data used in Section 5.2 are presented below. In the central experiment, each dataset was subsampled to contain 125 observations in total. In an ablation experiment, each dataset was subsampled contain 250 observations total. Then, we performed 5-fold coss-validation, where the dataset was randomly divided into five equal parts. In each fold, one part (20%) was used as the (central experiment: 25 samples, ablation experiment: 50 samples), while the remaining four parts (80%) formed the training set (central experiment: 100 samples, ablation experiment: 200 samples). From the training set, 20% (central experiment: 20 samples, ablation experiment: 40 samples) was further attributed to the validation set.
Colon. The first successful trials of adjuvant chemotherapy for colon cancer dataset was obtained from the survival package [46]. The dataset contains records of 1,822 observations with 15 covariates among which 49.23% are censored. All rows with missing values were excluded from the dataset.
NWTCO. The National Wilm's Tumor Study (NWTCO) was obtained from the pyсox pack- age [27]. The dataset contains records of 4,028 observations with 7 covariates among which 14.18% are censored.
GBSG. The Rotterdam and German Breast Cancer Study Group (GBSG) was obtained from the Pyсox package [27]. The dataset contains records of 2,232 observations with 7 covariates among which 43.23% are censored.
METABRIC.
The Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) dataset obtained from the pycox package [27]. The dataset contains records of 1,904 observations with covariates among which 42.07% are censored.
WHAS. The Worcester Heart Attack Study (WHAS) dataset was obtained from the sksurv pack- age [39]. The dataset contains records of 500 observations with 14 covariates among which 43.00% are censored.
SUPPORT. The Study to Understand Prognoses and Preferences for Outcomes and Risks of Treatment (SUPPORT) dataset obtained from the pyсoх package [27]. The dataset contains records of 8,873 observations with 14 covariates among which 31.97% are censored.
VLC. The Veterans administration Lung Cancer trial (VLC) dataset was obtained from the sksurv package [39]. The dataset contains records of 137 observations with 8 covariates among which 6.57% % are censored.
SAC The Sac 3 dataset from the simulation study in [28, Appendix A. I] was obtained from the Pyсox package [27]. The dataset contains records of 100,000 observations with 45 covariates among which 37.20% are censored.
J.2
Benchmark Methods
J.2.1
Benchmark Deep Survival Methods
All deep learning methods share the same neural network architecture, which is detailed in Section K. The benchmark deep survival models were trained using the Adam optimizer with learning rate selected via grid search. Batch normalization was applied, and dropout rate of 0.1 was used. Training was conducted for 1,000 epochs with batch size of 256.
MTLR. The Multi-Task Logistic Regression [51] was implemented using the MTLR class from the pyсox package [27].
38

---

DeepHit. The DeepHit method [32] was implemented using the DeepHitSSngll class from the pycox package [27]. The hyperparameters and were set 0.2 and 0.1, respectively. Those are the default values.
DeepSurv. The DeepSurv model [24] was implemented using the CoxPH class from the pycox package [27].
Logistic Hazard. The Logistic Hazard method [51] was implemented using the LogisticHazani class from the pycox package [27].
CoxTime. The CoxTime method [29] was implemented using the CoxTime class from the pycox package [27].
CoxCC. age [27].
The CoxCC method [29] was implemented using the CoxCC class from the pycox pack-
PME
The PMF method [28] was implemented using the PMF class from the pycox package [27].
PCHazard. The PCHazard method [28] was implemented using the PCHazard class from the pycox package [27].
BCESurv. package [27].
The BCESurv method [28] was implemented using the BCESurv class from the pycox
DySurv. The DySurv method [36] was implemented using the official code provided by the authors, com/munt wwweessiooic/DySSrrr hbtp:///lbb////b/ htt/b/lbbbmmmmm/// Maain/Models/Resee available https://github. http://gggthuu...mm..m..ttt
wxwapplle_eeleeee..eeeeea ipynb (Accessed on May 13 2025).
Sumo-Net. The Sumo-Net method [43] was implemented using the official code provided by the
authors, available at attpp:///iiie com/MrHufffSSmmo- -et (Accessed on July 25 2025).
DQS. The DQS method [50] implemented using the official code provided by the authors, available at https://gitie cwwccco/////// dqs (Accessed on July 25 2025).
J.2.2 Traditional Survival Methods
CoxPH. The Cox Proportional Hazards model [10] was implemented using the CoxPHFitted class from the lifelines package [11]. The Breslow method was used to compute the survival function.
Weibull AFT. The Weibull Accelerated Failure Time model [8] was implemented using the WeibullAFTFitEen class from the lifelines package [11].
RSE.
The Random Survival Forest [23] was implemented using the Bandom@urrvivalForeee class
from the sksurv package [39]. The number of trees in the forest is set to 1,000. The minimum number of samples required to split an internal node is 10, and the minimum number of samples required to be at a leaf node is 15. Those were the same byperparameters used in [36].
SSVM. The Survival Support Vector Machine [40] was implemented using the FastSurvivalse
class from the sksurv package [39]. The optimal regularization hyperparameter a was selected grid search by evaluating model performance on the training set using the C-index. This method does not allow for estimation of the survival function. Predicted ranks were used as risk scores for computing the C-index.
39

---

J.3 Evaluation metrics
C-index. Let (t) be the predicted risk score of observation with covariates x, time The C-index estimate [19] is given by
C-index E di (misge (R(a()))>))))) + : (0(0000000)000
Ei=l Ezst 6e (,,r,))
Let (t( be the predicted survival function of observation with covariates x, at time When the predicted risk score is taken to be the negative of the survival function, i.e., (:(() = (t), the C-index is referred as the Antolini's C-index [3] and is found with
E Ents 6g Inis (1 (1 ((,((0)8,)))0)) + #1 (S....).,
C-index
N Ezze 6; (m:s))
The C-index is obtained using the CencordanceInder class from the TorchSurv package [37].
IPCW Integrated Brier Score. Let (t))) be the predicted survival function of observation with covariates X; at time f. Let the inverse probability censoring weight (IPCW) at time be defined as the inverse of the probability of being uncensored, E(I) - (t), where C()) denotes the Kaplan-Meier estimate of the censoring survival function. Under right censorship. the IPCW Brier score (BS) [17] at time is given by
N Σ (1 N
IPCW BS(()=
+ E(t) (1 -(()))
(A29)
The IBS is the integral of the Brier Score in (A29). The IPCW weights and the IPCW IBS are com- puted using the get_ipew function and the BrierScore class from the TorchSurv package [37].
Distribution Calibration. D-Calibration [18] evaluates whether predicted survival probabilities at observed times are uniformly distributed. For an individual i, let S;(y;) be the predicted survival probability their event or censoring time y:- Under perfect calibration, we expect:
(y:)) Uniform(0, 1).
(A30)
The predicted probabilities are binned into B quantiles, and a histogram is constructed over both event and censored observations. For censored data, the probability distributed proportionally across bins beyond the censoring time. A chi-squared test compares the resulting histogram to the expected uniform distribution, and the p-value reflects how well the survival model is calibrated. The D-Calibration is obtained using the LifelinesDy weeiinnneeeaaaa d_callllllll calibration () class from the SurvivalEVAl package [41].
Kaplan-Meie (aM)-Calirration, as introduced by [9]. evaluates how well the average predicted survival curve from model aligns with the empirical KM survival curve. Let Sivg() (t) denote the model's 's average survival probability at time f, and (tM(t) the KM estimate. The KM calibration score is defined as the normalized integrated mean squared error (MSE);
Kaplan-Meier Calibration.
Timax
Tiax (Smg(()- Sem(()) dl.
KM-Calibration
(A31)
This score lies in 1), where indicates perfect calibration, and values near 0.25 represent uninformative predictions. The KM-Calibration is obtained using the LifellinesEvaluator km_ -tm.alibratiinn(( class from the SurvivalEVAL package [41].
40

---

K Implementation Details
Code availability.
The code is available on the GitHub repository https://githubb.co com/ under the MIT License.
Architecture.
We employed feedforward neural network with two hidden layers, each containing 16 neurons and using ReLu activations. The input of the network for observation I...., N is the pair (t, x,).
Time normalization. time value by the maximum observed time in the training set.
The observation period is normalized to the interval [o, 1] by dividing each
EM algorithm. The parameters are initialized so that they match their prior expected values. Specifically, set O 0 and gi / Bo. The maximization step of the EM algorithm is performed using the L-BFGS-B algorithm. The EM algorithm considered to have converged when the relative change in the Q-function between consecutive iterations falls below tolerance threshold of 10-6 for two successive iterations.
CAVI algorithm. The hyperparameters are initialized that the expected values of the model parameters match the MAP estimates. Specifically, we set ( ΦMAP x B, and (eL, y) (MMA.Im). The CAVI algorithm is considered to have converged when the relative change between successive parameter estimates falls below tolerance threshold of 10-6
Integral approximation. The integrals required to compute the Q-function in the EM algorithm, well as those involved in the optimal variational updates of ø and in the CAVI algorithm, approximated using the trapezoidal rule.
Prior and p. For all experiments, we fix the hyperparameters of the prior distribution over o, given in (6), to be ogp, Bo = 1. Furthermore, we fix 1.
Machine. The experiments were conducted NVIDIA RTX A6000 GPUs with 48GB of memory.
Table Al reports the running time for a single fold on the Colon dataset at varying sample sizes. All folds and datasets were processed in parallel across multiple GPUs to ensure consistent timing.
Running time
N 25
N 125
N 250
Method
21-166
21-1661
21.-60
MTLR [51]
0.104
0.102
0.100
0.075
0.093
0.118
0.075
0.099
0.091
DoopHit [32]
0.164
0.139
0.140
0.120
0.126
0.142
0.122
0.145
0,135
DeepSurv [24)
0.082
0.098
0.072
0.063
0.079
0.072
0.071
0.080
0.101
Logistic Hazard [16]
0.077
0.079
0.076
0.070
0.090
0.078
CoxTime [29]
0.139
0.167
0.114
0.113
0.175
0.134
0.136
0.135
0,186
CocCC [29]
0,130
0.118
0.094
0.100
0.119
0.114
0.107
0.132
0,149
PMF [28]
0.129
0.119
0.094
0.069
0.104
0.095
0.072
0.087
0,125
PCHazard [28]
0.083
0.102
0.083
0.087
0.092
0.084
0.091
0.091
0.091
BCESury [30]
0.081
0.082
0.092
0.060
0.078
0.092
0.062
0.087
0.079
DySurv [36]
0.049
Sumo-Net [43]
0.058
0.058
0.060
0.062
0.071
0.070
0.073
0.073
0.078
DQS
0.019
0.020
0.024
0.023
0.025
0.025
0.029
0.030
0,030
NeuralSary (Ours)
0.621
1.165
0.566
3.926
16.454
22,673
7,373
87.098
141.389
Table A1: Inference runtime for the Colon dataset (in minutes). The central analysis presented in
Table A5-A6 is for N = 125 and a MNP with layers (2I.) and 16 units (16U). The ablation study presented in Table A7-A8 for N 250 and MNP with 2L and 16U.
41

---

L
Related Work
Survival analysis methodologies have evolved significantly over the past decades, encompassing parametric, semi- parameeree,r,c pa-paaametric, and more recently, deep learning- -ased approaches. We review these developments, focusing on their applicability to high-dimensional data and uncertainty quantification capabilities.
Parametric and Semi-parametric Traditional Models. Traditional survival models often impose parametric or parimerraettic assumptions on the hazard function. The Accelerated Failure Time (AFT) model [8] assumes linear relationship between covariates and the logarithm of survival time, with parametric baseline distributions (e.g., Weibull). While interpretable, such models struggle with high-dimensional data and nonlinear covariate effects. The Cox Proportional Hazards (CoxPH) model [10], a semi- paniiparametric approach, avoids specifying the baseline hazard but assumes propor- tional hazards. Though widely adopted, CoxPH's linear predictor and proportionality constraints limit its flexibility in complex data regimes.
Non-parn -or-maraett Traditional Models. To mitigate parametric assumptions, non- por-mmametrric meth- ods like Random Survival Forests (RSF) [23] and Survival Support Vector Machines (SSVM) [40] emerged. RSF leverages ensemble learning for risk stratification but faces challenges in high- dimensional settings due to greedy tree induction. GP survival models [13] offer flexibility by modeling the hazard function nonparametrically. with inherent uncertainty quantification. Existing work has sought to address the cubic complexity in sample size of GPs by introducing variational inference techniques [25]. However, GPs remain fundamentally limited in scalability, particularly struggling with nigh-dimensional inputs and lacking the capacity to learn hierarchical representations, such as those required in image-based tasks [44].
Deep Survival Models. The advent of deep learning revolutionized survival analysis by enabling automatic feature learning from high-dimensional inputs. DeepSurv [24] extended CoxPH with neural networks, while DeepHit [32] employed multi-task learning for competing risks via discrete-time hazards. Discrete-time methods, including MTLR [51] and PCHazard [28), discretize the time axis
to simplify likelihood computation, with recent advances like DySurv [36] incorporating conditional variational inference for dynamic prediction. Sumo-Net [43] introduces partially monotonic NN that directly optimizes the right-censored log-likelihood, which is proven be a strictly proper scoring rule achieviii strong log- aigliikeld performance. DQS [50] formulates survival prediction using extensions of strictly proper scoring rules that remain proper under discrete-tine survival settings. Despite their predictive prowess, these models rely on frequentist training. yielding point estimates without uncertainty quantification, a significant shortcoming in safety-critie applications. Comprehensive reviews [49] highlight the rapid growth of deep survival methods but underscore their neglect of probabilistic uncertainty.
Bayesian and Uncertainty-Aware Approaches. Bayesian methods provide a natural framework for uncertainty quantification but have seen limited integration with deep survival models. GP-based approaches [13, 25] inherit GP limitations in scalability and nigh-dimensional processing. Recent works like BCESurv [30] explore bootstrap confidence intervals, yet these post-hoc approximations lack the coherence of Bayesian posteriors. Consequently, existing Bayesian survival models either sacrifice scalability for uncertainty quantification or compromise on model flexibility, leaving a critical gap in high-dimensional, uncertainty-awee -aware survival analysis.
Summary. While parametric and semi-parametric models provide interpretability. they falter in high-dimensional, nonlinear regimes. parameatetri methods like RSF and GP improve flexibility but face scalability challenges. Deep learning approaches excel at feature extraction yet lack princi- pled uncertainty quantification. Bayesian methods, though theoretically sound, remain confined to traditional architectures or partial approximations first scalable, deep Bayesian survival model that harmonizes neural networks with full probabilistic uncertainty, addressing a critical need in modern applications.
Our work bridges this divide by proposing the
42

---

Method
Uncertainty (Bayesian)
Continuous Time
Deep Learning
CoxPH [10]
AFT [8]
RSF [23]
X
SSVM [40]
X
GP survival models [13, 25]
MTLR [51]
X
DeepHit [32]
X
DeepSurv [24]
Logistic Hazard [16]
X
CoxTime [29]
CoxCC [29]
PMF [28]
PCHazard [28]
BCESurv [30]
DySurv [36]
X
X
Sumo-Net [43]
DQS [50]
X
NeuralSury (Ours)
Table A2: Summary of Survival Analysis methods: Bayesian Uncertainty Quantification, Time Domain, and Deep-Learning Status.
43

---

M Further Results
M. Synthetic Data Experiment
N = 25
N
50
N = 100
150
Method
C-index f
IPCW
C-index 1
IPCW
C-index
IPCW IBS
C-index
IPCW IBS
MTLR [51]
0.560
0.284
0,505
0.239
0,491
0.171
0,542
0.17
DeepHfit [32]
0.473
0.239
0.469
0.214
0.502
0.171
0.574
0.114
DeepSury [24]
0.492
0.313
0.471
0.241
0.507
0.517
0.169
Logistic Hozard [16]0.477
0.297
0,498
0.256
0.507
0,499
0.176
ConTine [29]
0.424
0.284
0,532
0.273
0.52
0.184
0.575
0.118
CoxCC [29]
0,421
0,268
0.497
0.229
0,526
0.128
0,513
0.109
PMF [28]
0,573
0.261
0,551
0,334
0,523
0.168
0.607
0.184
PCHazard [28)
0.477
0.337
0,501
0.249
0.467
0.174
0.486
0.193
BCESary [30]
0,545
0.287
0.16
DySurv [36]
0,399
0.237
0,491
0.174
Sumo.Net Net [43]
0.473
0.223
0,503
0.179
0,588
0.127
0.495
0.113
DQS [50]
0.435
0.326
0.48
0.252
0.525
0.131
0,556
0.124
NeuralSury (Ours)
0,378
0.196
0.554
0.160
0.589
0.126
0.589
0.106
Table A3: Performance comparison of survival models over synthetic data. The best results for each metric are shown in bold, and the second-best results are underlined. indicates higher is better; indicates lower is better.
N 25
50
Method
D-Calibration (p-value)
KM-Calibratin +
D-Calibration (p-value)
KM-Calibration J.
MTLR [51]
0.000 (x)
0.032
0,000 (X)
0,056
Deepffit [32]
(x000(()
0.029
0,000 (x)
0,168
DeepSury [24]
0.000 (x)
0.009
0.000 (x)
0.036
Logistic Hazard [16]
(x0)
0.037
0,000 (X)
0.062
CoxTime [29]
(x)
0,000 (x)
CoxCC [29]
(X0)
0.008
0.000 (X)
PMF [28]
0.000 (X)
0.032
0,000 (X)
0,071
PCHaaad [28]
0.000 (x)
0.037
0,000 (x)
0.073
BCESurv [30]
0.000
0.022
0.000 (x)
0,014
DySurv [36]
0.000 (x)
0.201
0.000 (X)
0.241
Sumo-Net [43]
0.132 (V)
0.051 (V)
0.036
DQS
0.000(())
0.051
0.000 (X)
NeuralSurv (Ours)
0.833 (v)
0.014
0.539
0.012
N 100
N = 150
Method
Calibration (p- (p-value)
IM-Calibaation J
Caliibration (p value)
KM CM-Caaibrttion J
MTLR [51]
0.000
0.034
0.000
0.035
DeepHit [32]
(x0)
0.030
0.000 (x)
0.040
DeepSurv [24)
0.000(X)
0.003
0.007 (x)
0.003
Logistic Hazard [16]
0.000 (x)
CoxTime [29]
0.000 (x)
0.011
0.346
0.005
CocCC [29]
0.000 (X)
0.001
0.001
0.001
PMF [28]
0.000
0.033
0.000 (x)
0.034
PCHazard [28]
(x)
0.056
0.000 (x)
0.066
BCESury [30]
(x))
0.038
0.000 (x)
0.034
DySurv [36]
0.000 (x)
Sumo-Net [43]
0.683 (V)
0.004
0.345
0.003
DQS [50]
0.005 (X)
0.032
0.000
0.031
NeuralSary (Ours)
(v419(())
0.004
0.639
0.004
Table A4: Performance comparison of survival models synthetic data (part 2). A checkmark (V) indicates that the null hypothesis of perfect D-Calibration was not rejected at a = 0.05 (model considered well-calibrated); cross (x) indicates rejection of D- D-Calibration (model considered not well-calibrated). The best results for the KM-Calibration are shown in bold, and the second-best results are underlined. J. indicates lower is better.
44

---

M.2
Real Data Experiment
M.2.1
Central Analysis
COLON
GBSG
Method
C-index
IPCW IBS 4
C-index
IPCW IBS L
C-indox
IPCW IBS J
[51]
0.298
DeepHfit [32]
0.28
0.243
0.309
DeepSurv [24)
0.572
0.326
0.523
0.289
0.618
0.252
Logistic Hazard [16]
0.490
0.321
0.541
0.317
0.618
0.296
TimTime [29]
0,578
0.277
0.533
0.307
0.599
0.285
CoaCC [29]
0.584
0,289
0.575
0.257
0.646
0.240
[28]
0.250
PCHazard [28)
0.538
0.297
0.541
0.291
0.249
BCESury [30]
0.491
0,302
0,616
0.277
0
0.273
DySurv [36]
0.488
0.536
0.561
0.465
0.572
0.485
Sumo-Net [43]
0.485
0.24!
0.447
0.223
0.476
0.250
DQS [50]
0.246
0.564
0.261
0.229
NeuralSurv
NWTCO
WHAS
SUPPORT
Method
C-index
IPCW IBS L
C-index
IPCW IBS L
C-index 1
IPCW IBS L
MTLR [51]
0.592
0.315
0.357
Deeptiit [32]
0.516
0.296
0.510
0.303
0.452
0.341
DeepSurv [24]
0.527
0,248
0.654
0.281
0.505
0.354
Logistic Hazard [16]
0,512
0.298
0.545
0.315
0.536
0.378
CoxTime (29]
0.550
0.199
0.678
0.250
0.547
0.327
CoaCC (29)
PMF [28]
0.312
0.299
0.399
PCHazard (28]
0,551
0.209
0.527
0.291
0.514
0.335
BCESSur [30]
0.530
0,272
0.548
0.292
0.446
0.398
DySurv [36]
0.402
0.683
0.424
0.523
0.525
0.342
Sumo-Net [43]
0.595
0.170
0.556
0.260
0.444
0.289
NeuralSurv (Ours)
0.712
0.166
0.602
0.233
0.599
0.333
VL.
SAC3
Method
C-index
IPCW IBS 4
C-index
IPCW IBS
MTLR [51]
0.432
0.299
0.471
0.276
DocpHit [32]
0.409
0,236
0.456
0.289
DeepSurv [24)
0,186
0,530
0.264
Logistic Hazard [16]
0.413
0,272
0,480
0.348
CoxTime [29]
0.671
0,212
0.485
0.276
CoxCC (29)
0.645
0.169
0.533
0,261
PMF [28]
0.445
0,284
0.472
0.270
PCHazard [28)
0.502
0.294
0,527
0,276
BCESury [30]
DySurv [36]
0.436
0.162
0.476
0.303
Sumo-Net [43]
0.527
0.157
0.457
0.237
DQS [50]
0.568
0,218
0.481
0.293
NeuralSary (Ours)
0.667
0.142
0.532
0.204
Table A5: Performance comparison deep survival models over five different train/test splits each dataset. The best results for each metric shown in bold, and the second-best results are underlined. t indicates higher is better; J indicates lower is better.
45

---

COLON
METABRIC
Method
D-Calibration (p-value)
KM-Calibration 4
D-Calibration (p-value)
KM-Calibratio |
MTLR [51]
(x)
0.016
0.000 (x)
0.025
DcepHit (32)
0.001 (x)
0.089
0.101
0.064
DuepSurv [24]
0.047 (x)
0.024
0.006
0.012
Logistic Hazard [16]
0.002 (x)
0.019
0.000 (x)
0.026
CoxTime [29]
(x)
0.011
0.012 (x)
0.012
CoxCC [29]
0.011 (x)
0.011
0.013 (x)
0.019
PMF [28]
0.001 (x)
0.017
0.002 (x)
0.017
PCHazard [28]
0.007
0.029
0.008
0.036
BCFSary [30]
0.000 (x)
0.012
0.000
0.021
DySurv [36]
0.000 (x)
0.362
0.000 (x)
0.306
Sumo-Net [43]
0.741 (V)
0.014
0.600 (v)
0.008
DQS [50]
0.381 (V)
0.017
0.135 (v)
0.020
NeuralSury (Ours)
0.594 (V)
0.020
0.661 (V)
0.012
GBSG
NWTCO
Method
Dalibbrationn (p value)
KM XM-Calibration 4
D. Calibration value)
KM CMlCbbrbttti +
MTLR [51]
(x)
0.011
0.570
0.011
DeepHfit [32]
(x)
0.149
0.575
0.014
DcepSurv [24]
0,084 (V)
0.009
0,887 (v)
0.003
Logistic Hazard [16]
0.008 (x)
0.017
0,397 (v)
0.012
CoxTime [29]
0.247 (V)
0.009
0.883 (v)
0.003
CoxCC [29]
0.256 (V)
0.005
0.954 (V)
0.003
PMF [28]
0.003 (x)
(V)
0.013
PCHazard [28]
0.063 (v)
0.024
0.697 (V)
0.008
BCESury [30]
0.001 (x)
0.015
0,312
0.010
DySurv [36]
0.000 (x)
0.369
0.000 (x)
0.260
Sumo-Ner Net [43]
0.524 (V)
0.010
0.991 (V)
0.004
DQS [50]
0.234 (V)
0.012
0.916
0.005
NeuralSury (Ours)
(v735(())
0,920 (V)
0.007
WHAS
SUPPORT
Method
Callbbrationn (p- value)
KM-Calibration J
D.Calibration (p-value)
KM-Calibration, 1
MTLR [51]
0.001 (x)
0.035
0.000
0.099
DoopHit [32]
0.195 (V)
0.067
0.000 (x)
0.124
DoepSurv [24)
0.076 (V)
0.025
0.000 (x)
0.015
Logistic Hazard [16]
0.013 (x)
0.033
0.000 (x)
0.095
CoxTime [29]
0.035 (x)
0.017
0.000 (x)
0.016
CoxCC [29]
0.106 (V)
0.019
0.000 (x)
0.020
PMF [28]
0.003 (x)
0.031
0.084
PCHazard [28]
(v)
0.022
0.000 (x)
0.070
BCESSrr [30]
0.000 (x)
0.025
0.000 (x)
0.088
DySurv [36]
(x000
0.281
0.000 (x)
0.099
Sumo-Net [43]
(v)
0.143 (v)
0.017
DQS [50]
0.081 (V)
(x)
0.038
NeuralSary (Ours)
(v3
0.031
0.063 (v)
0.083
VLC
SAC3
Method
D- Calibration (a-ue)
KM-Calibration 4
D- D.Calibration (p-value)
KM AM-Calibration +
MTLR [51]
0.000 (x)
0.072
0.000 (x)
0.020
DeepHHit [32]
0.000 (x)
(x)
0.094
DeepSurv [24)
0.004
(x)
Logistic Hazard [16]
0.000(x)
0.078
0.000
0.039
CoxTime [29]
0.021 (x)
0.010
0.000 (x)
0.014
CoxCC [29]
0.004 (x)
0.011
0.000
0.016
PMF [28]
0.000 (x)
0.077
0.000 (x)
0.016
PCHazard [28]
(x)
0.080
0.000 (x)
0.029
BCESurv [30]
(x)
0.074
(x)
0.021
DySurv [36]
0.077 (V)
0.028
0.000
0.146
Sumo-Net Net [43]
0.131 (V)
0.011
0.478 (v)
0.016
DQS [50]
0.000(X)
0.054
0.021 (X)
0.034
NeuralSury (Ours)
0.436 (V)
0.013
0.624 (V)
0.016
Table A6: Performance comparison of deep survival models over five different train/test splits of each dataset (part 2). A checkmark (v) indicates that the null hypothesis of perfect D-Calibration was not rejected at a = 0.05 (model considered well-calibrated); cross (x) indicates rejection of D-Calibration (model considered not well-calibrated). The best results for the RM-Calibration are shown in bold, and the second-best results are underlined. J. indicates lower is better.
46

---

M.2.2 Ablation Study with N 250
COLON
METABRIC
GBSG
Method
C-index f
IPCW
C-index
IPCW
C-index
IPCW
MTLR [51]
0,545
0.291
0.572
0.290
0.567
0.312
DoopHlit [32]
0.564
0.284
0.545
0.301
0.563
0.272
DeepSurv (24)
0.600
0.295
0.6U5
0.265
0.531
0.277
Logistic Hazard [16]
0.501
0.289
0.553
0.252
0.562
0.287
CoxTime [29]
0.621
0,259
0.621
0.264
0,578
0.255
CoxCC [29]
0.640
0.277
0.610
0.254
0.565
0.244
PMF [28]
0,541
0.291
0.554
0.300
0.537
0,304
PCHazard [28]
0,549
0.280
0.561
0.246
0.524
0.295
BCESury [30]
0.537
0.289
0.565
0.289
0.554
0.301
DySurv [36]
0.516
0.506
0.508
Sumo-Net [43]
0.529
0.273
0,473
0.27
0.471
0.255
[5Q(
0.593
0.267
0.600
0.228
0.562
0.237
NeuralSary (Ours)
0.601
0.215
0.543
0.198
0.546
0.212
Table A7: Performance comparison of deep survival models the ablation study with 250 observa- tions, over five different train/test splits of each dataset. The best results for each metric are shown in bold, and the second-best results underlined. 1 indicates higher is better; 4 indicates lower is better.
COLON
METABRIC
Method
D-Calibration (p-value)
KM-Calitration,
D-Calibration (p-value)
KM-Calibration
MILR [51]
0.000(X)
0.007
0.000(
0.017
Deepl [32]
0.000 (x)
0.068
0.044 (x)
0.072
DeepSurv [24]
(x)
0.001 (x)
Logistic Hazord [16]
Cox Tone [29]
0,022
0.005
0.001 (X)
0.008
CoxCC [29]
0.001 (x)
0.003
0.000 (x)
0.006
PMF [28]
0.000 (x)
0.008
0.000 (X)
0.008
PCHazard (28)
0.001 (x)
0.018
0.003 (x)
0.027
BCESurr [30]
(x000(())
DySurv [36]
0.000 (X)
0.273
Sumo-NNt [43]
0.323 (v)
0.007
0.192 (V)
0.007
DQS [50]
0.086 (v)
0.009
0.137 (V)
0.013
NeuralSurv (Ours)
0.404 (v)
0.011
0.708 (V)
0.013
GBSG
Method
D- DaCabibration (p-value)
KM KM-Calibraion +
MTLR [51]
0.000 (x)
0.010
DeepHfit [32]
DeepSurv [24]
Logistic Hazurd [16]
0.000
0.016
CoxTime [29]
0.171 (V)
0.004
CoxCC [29]
0.059 (V)
0.006
PMF [28]
0.000 (x)
0.010
PCHazard [28]
0.004(())
0.015
BCESurv
DySurv [36]
0.000
0.336
Sumo-NNt Net [43]
0.337 (V)
0.004
DQS [50]
0.084 (V)
0.008
NeuralSury (Ours)
0.617 (v)
0.003
Table A8: Performance comparison of deep survival models on the ablation study with 250 obser- vations, over five different train/test splits of each dataset (part 2). A checkmark (v) indicates that the null hypothesis of perfect D D-Calibration was not rejected at a = 0.05 (model considered well- calibrated); cross (x) indicates rejection of D-Calibration (model considered not well-calibrated), The best results for the KM-Calibration are shown in bold, and the second-best results are underlined. I indicates lower is better.
47

---

M.2.3
Comparison to Traditional Survival Models
COLON
NWTCO
GBSG
Method
IPCW IRS 4
C-index
IPCW TRS !
C-index
IPCW IRS !
CoxPH [10]
0.669
0.192
0.710
0.136
0.694
0.171
Weibull AFT [8]
0.681
0.198
0.697
0.134
0.673
0.179
RSF 23)
0.590
0.210
0.156
0.193
[40)
METABRIC
WHAS
SUPPORT
Method
C-index
IPCW IBS 1
C-index
IPCW IBS L
C-index
IPCW IBS L
CoxPH
0.207
Weibull AFT [8]
0.658
0.172
0.622
0.224
0.650
0.239
RSF [23)
0.587
0.189
0.683
0.209
0.601
0.225
SSVM [40]
0.649
0.653
0.636
VLC
SAC3
Method
C-index
IPCW IBS L
C-index 1
IPCW IBS
CoxPH
0.697
0.125
0.569
0.190
Weibull AFT [8]
0.690
0.127
0.607
0.287
RSF [23]
0.687
0,139
0.487
0.182
SSVM [40]
0.698
0.504
Table A9: Performance comparison traditional survival models over five different train/test splits of each dataset. 1 indicates higher is better; 4 indicates lower is better. The SSVM method does not provide estimates of the survival function; the predicted ranks are used for the corresponding C-index evaluations while the IPCW-IBS metric cannot be computed.
COLON
NWTCO
GBSG
Method
D-Calibration
KM-Calibration
D-Calibeation
KM-Calibration J
D-Calibration
KM-Calibeation
(p-value)
(p-value)
(p-value)
CoxPH [10]
0.913 (V)
0.006
0.979 (V)
0.003
0,950 (v)
0.004
Weibull AFT 18]
0.788 (
0.010
0.986 (V)
0.004
0.767 (V)
0.008
[S3]
0.009
0.999 (V)
0.854
METABRIC
WHAS
SUPPORT
Method
D-Calibration
KM-Calibeation 4
D-Calibeation (p-value)
KM-Calibration J
D-Calibration
KM-Calibeation
(p-value)
(p-value)
CoaPH
0,846 (v)
0.008
0.730 (V)
0.017
0.354 (V)
0.010
Weibull AFT [8]
0.759 (V)
0.006
0.650 (V)
0.018
0.420 (V)
0.012
RSF (23)
0,746 (V)
0.009
0.625 (V)
0.022
0,593
0.010
VLC
SAC3
Method
D Calibration
KM SM4-Calibrtionn J
Calibration
KM-Calibration
(p (p-value)
(p valule
[10]
0.597 (<))
(V)
0.005
Weiboll AFT [8]
0.759 (X)
0.007
0.038 (V)
0.018
RSF [23)
0.414 (X)
0.013
0.706 (V)
0.012
Table A10: Performance comparison traditional survival models over five different train/test splits of each dataset (part 2). A checkmark (v) indicates that the null hypothesis of perfect D-Calibration was not rejected
0.05 (model considered well-calibrated); cross (x) indicates rejection of D-Calibration (model considered not well-calibrated), The best results for the EMM-alibration are shown in bold, and the second- best results are underlined. indicates lower better.
48

---

M.3 Prior Sensitivity Analysis
COLON
Gamma Prior of
Gamma Posterior of
Posterior Median
C-index f
IPCW IBS
D. D-Calibration (p-value)
KM-Calihratie +
and 95% CI
(okg 1,
8 78.789)
0.768 [0.591, 0.978)
0.671
0.594 (V)
2. 20
45. 0= 79.789) 38.333, 3= 78.289)
0.571 [0.421. 0.753]
0.593
0.237
0.601 (
0.5, 9o- 0.5)
0.485 [0.347.
0.512
0.229
0.715 (V)
METABRIC
Gamma Prior of
Gamma Posterior of
Posterior Median
C-index f
IPWW IBS
D.Calibration (p-value)
KM- KM-Calibration J
95%
(cg
1
(6 48.406, 8= 70.176)
0.685 [0.509, 0.897]
0.584
0.212
0.661 (V)
0.012
(Og 2, 80
( 48.034,
71.166
0.670 [0,498, 0.879]
0.536
0.201
0.819 (V)
0.009
(ara
0.5, 90
0.5)
(6 47.347, 8 69.676)
0.675 [0.500, 0.886]
0.553
0.200
0.802 (V)
0.009
GBSG
Gamma Prior
Median and 95% CI
f
Calibration (p-value)
J
(cig 1, 1, 1)
(á = 58.650, 8 85.462)
0.682 [0.522, 0.8731
0.657 0.602
0,188
0.735 (V)
0.009
(org 2. 80 2)
(& = 62,386 = 86-462)
0.718 [0.5 (0.554, 0.91 I]
0.195
0.808 (V)
0.010
(Okg 0.5, Bg 0.5)
(& = 57,54, 3 84.922)
0.673 (0.514, 0.863)
0.665
0.189
(V)
0.010
Table All: Prior sensitivity analysis using priors with double and half the original variance..
49

---

N Proofs
N.1 Proof of Theorem 3.1
Before proving Theorem 3.1 we must show some intermediate results.
Lemma N.1. Assume that for each i = . the function g(-x;; E C(O, yil x R.) Then, it follows that
J Xo(1, <
for every = I....,N..
Proof. Fix an arbitrary index E (1...., N). From Section 2.3, recall that po is the probability density function of multivariate normal distribution with zero mean and identity covariance matrix Im. Our goal is to show that the normalization factor (t, x;) admits a strictly positive lower bound on [0, yel. from which the integrability of (o xx; Φ) will follow.
Step 1: Continuity of (t, x,) on [0, vil. Fix any to E [0, v.l and let (tx)x21 be a sequence in [0, yil such that in to as n + oc. Define, for each n, the functions
hn
(g(t,, x,; 0 (O), n21, h(0) == a (g(1a, x,; (O)
Since x): E C([O, di] Rm) and the sigmoid a (-) is continuous function, it follows that lim hn (0) = h()) 100000
pointwise for all 0 E R Furthermore, observe that
hn (0) S pe (0)
since 0 < a()) < 1. Because po (0) integrates to 1 over IR" we may apply the Dominated
Convergence Theorem (DCT) to conclude that:
lim Z(b, X.) = lim
Ir, (0)do
h(0)d0 xi).
Rmm
Rm
Since lo was arbitrary in [0, W:]. Z is continuous everywhere on that interval.
Step 2: Strict positivity of Z(t, x,) on [0, vil- For each fixed t [0, v.l since (g(t,x;;0) >
and pe (0) > 0 for all 0 E 12mh, we have:
/ a x,; 0))o > 0. pmm
2(4,
Since x(),x;) is continuous and strictly positive function the compact interval [0, vil. the Weierstrass Extreme Value Theorem ensures that attains minimum on this interval. Define:
min Ztt x;)
0
AE 10
Step Integrability of (t, d). Note that for all E y:). we have
\o(t:
>o(t:
to x):
Z(t,x;)
2²
It is straightforward to verify that Ao(t: Φ) is integrable on [0, vil. therefore it follows that
y
No Ao(t,x;))de < 1 1 (o(t () co.
This completes the proof.
50

---

Our next result verifies condition needed for applying Campbell's 's Theorem in the proof of Theo- rem 3.1. To establish this, will use the following Pólya-Gamma Gamma identity:
E Bew (wl 1,0) [
(A32)
4
which follows by taking the limit c - 0 in equation (A4). Alternatively, to prove (A32), one can start from the representation in equation (A apply Tonelli's theorem to interchange expectation and infinite summation, and then invoke the series identity
x²
(k 2)
2
We are now ready to present our next result.
Lemma N.2. Assume that for each I,...,A the function g(,x;:-)) C(OO, vil x R.). Then, with probability 1 the sum
H(V,)
Σ
f(w,, -(t;,,;;;0)))
((,w))
absolutely convergent for every I....,N.. N.
Proof. Fix an arbitrary index i E {L...., N}. Recall the definition of f(w, 2) from (8). From Theorem C.3, it suffices to show
I,
min([s(w, -0(, x,; 0))1, 1) o)dddl < co.
(A33)
Since R4, then it follows from the triangle inequality that
min(][(w, - (1, x,; 1) < [/(w, -g(,, x,; 0))1
lg(t. x;: 0)
x1;)x;;0²²
+ log(2).
2
Hence remains prove finiteness of three integrals:
V1
(g(t, x5:
0)1
I =
As(t, wi @)dtdw
2
yi
000
g(t, x;; 0)²
In =
(l wi Adddll,
2
Ue
f >i wi wdwwd..
I3 = log(2)
I, is finite. Since gt, x;; 0) is continuous on the compact interval [0, yil, it is bounded by some
M > 0. Then,
1 = ==
(fo Pra(w)1, pro(w|1.0)dw) f n 1g(1, x,; 0)1
do x x)dddt s M l
co,
where the last inequality is Lemma N.1.
Ig is finite. Likewise g(l, x;; 0)² is bounded by some C over [0, yi] and Euwmpe(w/11000/)
(see (A32)), so
Ig = [Bunpee(w1..)) (w]) 1 g(1,x;;; 0)2
C (f, soll, (x;(,);;))a)) < 00.
x x1)
2
where the last inequality is Lemma N.I.
51

---

Ig is finite.
Finally,
log(2) Sh //" xo(lxx; co,
again by Lemma N.I.
Since I,, I2, Iz, all finite, the condition in (A33) is satisfied and the sum (v;) converges absolutely with probability 1.
The next result presents an integral identity which is key to proving the data augmentation scheme of Theorem 3.1.
Lemma N.3. Assume that for each L...., N the function g(,xx;)) C(O,, vol x RT). Then the double integral
J. (1 -i--elww-o (l-e pro(w[1, ()) x;; x;: dddwdd
converges, and in fact
y/ J.¹ (1 (1- (1-ellu, {,xh
ddwdt
>o
x:: xo(1;x;
x,; 0;) x;; (A34)
for every i = I,...,N. N.
Proof. Fix arbitrary index E (1,...,N), NN. By Lemma N.1
Ao x;; Ø) < oc.
Since 0 < <-)+) < 1, we have
()( (t(,xxxx < (o(1,x;;))
and therefore
xo x;; () x,; 0))dt < oc.
(A35)
This shows the finiteness of the right-hand side (A34). By combining o(2) (-2)) with (A5) we obtain that
V1 [
At(t,x;;0 (f((,;))))
=
|," U1 1 (1-e(wu-ottxxo))) O)) x;; @)dwdl. (A36)
Putting together the finiteness from (A35) with the equality of (A36) completes the proof.
We are now ready to prove Theorem 3.1.
Proof of Theorem 3.1. Fix an arbitrary index i E {1....,N} N). The joint expectation factors into two independent pieces:
Expectation over wi: This term recovers Ao(Us, x;; 0 x;; x;;
2. Expectation over v, This term recovers exp (- f
0))dt)
52

---

Step (1): Expectation over wiy Since d; E {0, 1},
(e (clonso()):0)) & f 1
d.
Hence,
(clouistuna:o))) pra(w:|1,0)du: = (f/ cluur. (wrlluil [1, 1)ddss
By the Pólya- Völya-Gamma Gamma identity (Eq. (A5)), the bracketed integral equals (g x;; 0)). Multiplying by \o x;; 6)th gives exactly
x.: ") (( x(((v;
Step (2): Expectation over V. By Lemma N.2 the random sum
(VV))=
Σ f(w),- g((,,, (t,w),W,
absolutely convergent with probability 1, and by Lemma N.3 the corresponding integral converges. Therefore, may apply Camphell's Theorem (Theorem C.3) together with the PG-sigmoid identity from (A34) to conclude
eep (-/-- t
II EV,,
Putting Steps (I) and (2) together reproduces precisely the two factors of the original likelihood p(y, 6,,,,, o, g). This completes the proof.
53


Reviews:
[{"review_id": "Jhx6S3yDKg", "reviewer": "Reviewer_QbSK", "summary": "This paper introduces a novel survival deep learning method within a Bayesian framework. Survival analysis is a critical tool for generating personalized time-to-event predictions while effectively handling censoring distributions.\n\nThe authors propose a new variational inference algorithm that leverages an augmented data scheme. This scheme utilizes Pólya-Gamma random variables to approximate the sigmoid function and Poisson variables to tractably manage the posterior distribution of their likelihood.\n\nThe proposed method, a new Bayesian Neural Network (BNN), demonstrates superior performance on several datasets when evaluated using the Integrated Brier Score (IBS) and the C-index (a discrimination metric). It consistently outperforms other deep learning methods while maintaining scalability.", "questions": "Theoretical part: \n- Why is a sigmoidal function exclusively assumed to map between 0 and 1? If \\lambda_0 is considered the baseline hazard, does the \"personalized\" prediction only decrease? How then is \\lambda_0 interpreted? Why was a Weibull distribution chosen for \\lambda_0 instead of adopting the more common approach seen in the Cox model?\n\nExperimental part: \nI am willing to change my grade if you could give me an answer on those points: \n- Can you do one experiment to show that your method outperforms the standard ML methods on a real-life dataset? Maybe on the KKbox dataset that contains 2.5M rows and is available on Kaggle or via pycox? \n- Can you add some calibration measures to your benchmark on the real datasets? \n- Can you add SumoNet, Yanamisawa, to the benchmark to understand their calibration power, and their abilities to predict the probabilities and their discriminative power?\n- In Table A.3, some C-index are lower than 0.5, which means that those models are worse than randomness (especially DeepHit on different datasets). Could the authors provide an explanation for these results?\n- Also, on support, all the IBS are worse than it is usually seen in the literature (e.g. Table 12 in Rindt). Could the authors clarify the methodology used to compute these scores?", "limitations": "yes", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 4, "strengths_and_weaknesses": "To enhance the accessibility of the paper, it would be beneficial to include a concise introduction to key concepts in survival analysis e.g. the hazard function (\\lambda(t)), the survival function (S(t), the censoring distribution)\n\nStrengths: \n- Using Bayesian methods with a prior is not really used in the literature, and using BNN is quite a novelty. \n- Integrating this method into survival analysis offers practitioners the valuable ability to incorporate domain-specific constraints or prior beliefs, which can significantly enhance model interpretability and robustness.\n- A lot of details are provided in the manuscript on the experimental part and on the theoretical part. \n- the proofs presented in the appendix are easy to follow, contributing to the overall rigor of the paper.\n\nWeaknesses: \n- My main concern pertains to Table A5, where standard machine learning methods consistently outperform the authors' method in terms of both IBS and C-index across all datasets. While the authors claim their method may excel in high-dimensional feature settings, no experimental evidence is provided.\n- You may need to clearly write your assumptions and the different key points in your reasonning.  \n- The Integrated Brier Score (IBS) can be biased due to the unknown censoring distribution inherent in real-life datasets. Other calibration metrics, such as Kaplan-Meier (KM) calibration [1, 2] and D-calibration [1, 2], have been developed in survival analysis to provide a more accurate assessment of calibration and should be considered.\n- The authors did not compare their method against approaches that train neural networks using proper scoring rules, such as SumoNet [3] or the method proposed by Yanagisawa [4]. This comparison would provide valuable insights into the relative performance of the proposed BNN.\n\n[1]: Haider, Effective Ways to Build and Evaluate Individual Survival Distributions, 2018\n\n[2]: one implementation in Python of those calibration measures: https://github.com/shi-ang/SurvivalEVAL\n\n[3]: Rindt, Survival Regression with Proper Scoring Rules and Monotonic Neural Networks, 2021 \n\n[4]: Yanagisawa, Proper Scoring Rules for Survival Analysis, 2024", "was_revised": true, "final_justification": "I have raised my grade thanks to their answer. \nThe theoretical weaknesses have been answered, as well as more experiments have been made to justify their method on the calibration viewpoint, and more baselines have been added.\n\nNevertheless, due to their low capacity to scale, other ML methods that still outperforms their solution, and a weird result for one method (Deephit performance), I am grading the paper with \"weak reject\".", "rebuttal": "> Q1: Concise introduction to key survival‐analysis concepts\n\nWe appreciate this suggestion and have added a brief overview of primary definitions and quantities at the start of Section 2; Appendix A still provides a more detailed treatment.\n \n> Q2: Concern regarding Table A.5, where standard ML methods outperform NeuralSurv in both IBS and C-index.\n\nWe apologize and recognise that interpreting the performance in Table A5 may not be straightforward and requires a nuanced perspective. It is important to note that the methods listed in Table A5, including the widely used Cox proportional hazards model, are designed for structured tabular data and typically do not require regularisation to mitigate overfitting. In our examples with NeuralSurv, we used a MLP with without regularization techniques such as L2 penalties, dropout, or batch normalization. Our primary reason to omit regularisation is to ensure a fair and transparent comparison to other competing neural network survival models. Indeed, you will note that it is not just our model that performs poorly compared to those listed in Table A5, but all the NN approaches, although ours is better performing than the others. To further reassure the reviewer of this, we have conducted an additional analysis in which NeuralSurv is implemented using a single linear layer instead of a full MLP. In this revised experiment, NeuralSurv performs much more comparably to the other methods in Table A5 and, in one case, achieves the best performance. Finally, it is critical to note that no method should work best on all data in these comparisons - some methods are based on proportional hazards, ours on a logistic hazard, and yet others on varying subtle assumption choices. These choices will match some data sets better than others. Our main point is NeuralSurv is, on average, highly competitive with all methods, and even those in Table A5 depending on the neural network choice. \n\n> Q3: While the authors claim their method may excel in high-dimensional feature settings, no experimental evidence is provided.\n\nWe thank the reviewer for raising this point. Our intention was not to claim demonstrated superiority in high-dimensional settings, but rather to highlight a theoretical strength of our framework. Since NeuralSurv employs neural networks as its predictive backbone, it inherits their established capacity to handle high-dimensional inputs, as evidenced in fields like computer vision and NLP. We acknowledge that our current experiments do not include a high-dimensional benchmark. This was due to practical constraints on computational resources, not a limitation of the method itself. Evaluating NeuralSurv in such settings is an important direction for future work. \n\nTo avoid overstatement, we have revised the sentence in the introduction from:\n\n“we address this critical need by introducing NeuralSurv, an architecture-agnostic, Bayesian deep-learning framework for survival analysis.”\n\nto the more measured:\n\n“we introduce NeuralSurv, an architecture-agnostic, Bayesian deep - learning framework for survival analysis which integrates with modern deep learning architectures.“\n \n > Q4: Clear statement of assumptions \n\nWe have added the following opening sentence to Section 2 to make all model assumptions explicit:\n\n“In this section, we outline the main assumptions underlying our model, NeuralSurv.”\n\nDetails of the Poisson–Pólya–Gamma augmentation appear in Section 3, and those of the variational inference algorithm in Section 4.\n\n> Q5: Why is a sigmoidal function exclusively assumed to map between 0 and 1? If \\lambda_0 is considered the baseline hazard, does the \"personalized\" prediction only decrease? How then is \\lambda_0 interpreted?\n\nWe thank the reviewer for this insightful question. In models where a sigmoidal function modulates the hazard, the baseline hazard $\\lambda_0(t)$ is typically interpreted as a reference or average hazard, reflecting the expected event rate without covariate adjustment. The sigmoidal transformation acts as a gating function that attenuates the baseline hazard based on individual covariates, thus modeling how personalized risk decreases relative to a worst-case or population-average scenario. This formulation is well established and has been widely adopted in both survival analysis [10] and point process models [9, 44, A, B, C], where it provides a balance between interpretability, flexibility, and computational tractability.\n\nTo clarify this modeling choice, we have revised Section 2 of the main text and added the following explanatory sentence and references:\n\n“Modeling hazard and intensity functions using sigmoidal transformations is common in applications such as survival analysis [10] and point process models [ 44 , 9, A, B, C]. This approach is popular due to the balance it offers between modeling flexibility and analytical tractability.”\n\n[A] Teh & Rao, Gaussian process modulated renewal processes. NeurIPS 2011\n\n[B] Adams et al., Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities. ICML 2009.\n\n[C] Aglietti et al., Structured variational inference in continuous cox process models. NeurIPS 2019.\n \n> Q6: Why use a Weibull distribution for \\lambda_0 instead of the Cox model?\n\nThe semi-parametric Cox model lacks a parametric baseline hazard, relying on partial likelihood methods that hinder Bayesian inference. In contrast, the Weibull models specify a parametric baseline, making them more suitable for Bayesian analysis. We chose the Weibull model for its flexibility, as also used in [10].\n \n> Q7: Can you do one experiment to show that your method outperforms standard ML methods on a real-life dataset?\n\nWe thank the reviewer for this suggestion. Our method has already been evaluated on eight real-world datasets (Section 5.2), covering a broad range of survival analysis settings.\nRegarding the KKBox dataset (~2.5M samples), we agree it's a valuable benchmark. However, NeuralSurv is a Bayesian model, which entails higher computational complexity compared to standard frequentist ML methods. As with Gaussian processes, Bayesian inference typically scales poorly with large datasets due to cubic time complexity. While recent advances [1] enable exact GP inference on datasets with up to a million points, they require extensive engineering optimizations that are not yet broadly applicable to general Bayesian deep learning models like ours. Applying NeuralSurv to KKBox would thus require substantial approximations or architectural changes beyond the scope of this work.\n\n[1] Wang et al., Exact Gaussian Processes on a Million Data Points. Advances in Neural Information Processing Systems. NeurIPS 2019.\n \n> Q8: Can you add SumoNet and Yanamisawa to the benchmark to evaluate calibration and predictive performance?\n\nWe thank the reviewer for suggesting SumoNet and DQS. We have added both methods to our benchmark comparisons across all experiments. The updated results are included in the revised table in the appendix. We are pleased to report that the first and second positions of NeuralSurv in term of C-index and IPCW IBS remain unchanged after adding the new benchmark methods.\n \n> Q9: Can you add some calibration measures to your benchmark on the real datasets?\n \nWe thank the reviewer for this suggestion. To address censoring bias in the IBS metric, we already applied IPCW adjustment as noted in Section 5 (l.249) and Appendix J.3.\n \nIn addition, we have now included two calibration metrics, KM calibration and D-calibration, in updated tables consistent with Tables A2–A4 in the Appendix. We are pleased to report that NeuralSurv shows proper calibration under the D-calibration (p-value) metric across all experiments, matching the performance of SumoNet as the only other method with consistent calibration.\nRegarding K-calibration, our method ranks fifth overall across all experiments. The best-performing method in terms of K-calibration is CoxTime.\n \n> Q10: In Table A.3, some C-index values are below 0.5. Why are these models performing worse than random?\n\nC-index values below 0.5 can occur when a model learns spurious or incorrect patterns, leading to risk predictions that are systematically misaligned with actual outcomes. This issue is especially pronounced in small-sample settings, where few comparable pairs make the metric highly sensitive to misrankings. As a result, even random noise or overfitting can yield C-index estimates worse than chance.\n \n> Q11: Also, on SUPPORT, all the IBS are worse than what is typically reported in the literature (e.g., Table 12 in Rindt). Could the authors clarify the methodology used to compute these scores?\n\nWe thank the reviewer for this observation. We note a methodological difference: Rindt et al. report the IBS without inverse probability of censoring weighting (IPCW), whereas our evaluation uses IPCW-adjusted IBS, which is more conservative and typically results in lower scores. This metric difference likely contributes to the discrepancy in reported IBS values.\n \nMoreover, in our experiments, we worked with subsampled versions of the datasets for computational reasons, using N=125 in the main analysis (Table A3) and N=250 in the ablation study (Table A4). These are considerably smaller than the full SUPPORT(N=8,873)  dataset, which may explain the lower C-index values. To investigate this further, we reran the benchmark models used by Rindt et al. on the full SUPPORT dataset. As shown below, we found that the Antolini’s C-index values closely matched those in our subset-based analysis, suggesting our results are qualitatively consistent. \nNo changes were made to the manuscript.\n\n|Support dataset | Our result | Rindt et al. |\n   |--------|------|------|\n|DeepSurv |   0.611  | 0.609  |\n|CoxTime | 0.609 | 0.613 |\n|CoxCC |0.614 | 0.608 |"}, {"review_id": "r8VvWtOAen", "reviewer": "Reviewer_iZcc", "summary": "The authors introduce a Bayesian neural survival analysis model. They use a neural network with sigmoid activation to modulate a parametric baseline hazard model based on observed covariates, which makes their model highly flexible and able to accommodate settings where the relationship between covariates and risk varies over time. After discussing the likelihood and placing suitable priors on neural network baseline hazard model parameters, respectively, they present a data augmentation strategy in which they represent the sigmoid function using a set of Polya-gamma random variables, then express the intractable double integral in the resulting likelihood as an expectation over a marked Poisson process. They use mean field variational inference to approximate the posterior distribution, noting as a limitation that this assumes neural network parameters are independent of the parameters of their hazard model. They conduct experiments on synthetic datasets, with the number of samples ranging from 25 to 150, followed by 8 real-world datasets, each subsampled to 125 samples. The evaluation shows that their method tends to improve on the concordance index and integrated Brier score compared to a comprehensive set of alternatives, and that the true survival function falls within their 90% credible interval in the synthetic data examples.", "questions": "My questions are related to the weaknesses described above. Namely:\n- How well (or poorly) does the method scale to larger datasets? Please be clear about this and show runtimes.\n- How exactly does the current data augmentation strategy relate to those already proposed in [9, 44, 2]? Please be up front and detailed about this.", "limitations": "The authors have described weaknesses but more could be done particularly around scalability & runtime, as I discuss in the \"weaknesses\" section.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- Addresses a true need for a fully Bayesian yet highly flexible survival analysis method\n- Method is flexible and elegant\n- Paper is well organized and generally well written\n- Experiments are compelling and illustrate the method's effectiveness and advantages for small datasets\n\nWeaknesses:\n- The method appears to be limited to very small samples. This is discussed briefly, notably in the conclusions section, and the authors note that extending e.g. to a mini-batch setup is a possible direction for future work. Still, they should expand on the details here and present runtime results for each experiment. If the method is not scalable in its current form, this is OK in my view but should be clearly stated.\n- A key benefit of the method is its flexibility in modeling the hazard function, but typically this kind of flexibility is becomes advantageous only for larger datasets, whereas their method is currently constrained to smaller datasets.\n- A substantial portion of the paper concerns the data augmentation strategy. The authors note in lines 149-150 that \"analagous data augmentation schemes have been proposed in [9, 44, 2]\". However, I think they should note this at the beginning of the section rather than the end, and they should also do more to explain how their approach is similar versus different from those already proposed. I think the authors' overall contribution is substantial even if this particular aspect is not very novel, so there is no need to obscure what is established versus new here.\n- The initial presentation of \\lambda_0 in section 2 is confusing, because it's not clear why or how \\lambda_0 depends on x. This becomes clear in 2.2, but I think it would be better to briefly comment on this when initially introducing \\lambda_0.", "was_revised": true, "final_justification": "The authors have addressed my concerns and clarified the contribution in the data augmentation.", "rebuttal": "> Q1: How well (or poorly) does the method scale to larger datasets? Please be clear about this and show runtimes.\n\nWe thank the reviewer for this helpful suggestion. In response, we have added a runtime analysis. We refer the reviewer to Question 2 of Reviewer 52KX. \n\n> Q2: How exactly does the current data augmentation strategy relate to those already proposed in [9, 44, 2]? Please be up front and detailed about this.\n\nWe thank the reviewer for requesting a more detailed explanation of the connection between our data augmentation strategy and existing work. \n\n- We added this text at the beginning of Section 3:\n\n“This approach builds on analogous strategies previously applied in other settings, including Bayesian inference for Sigmoid Gaussian Cox Processes [9], nonparametric Hawkes processes [44], and, in the case of Pólya–Gamma augmentation alone, mutually regressive point processes [2]. To the best of our knowledge, this is the first application of such a data augmentation strategy in the context of survival analysis. Furthermore, we are the first to provide a rigorous theoretical framework that establishes the validity of a method belonging to this broader class of augmentation-based approaches (see Theorem 3.1).”\n\n\n- Moreover, we have modified the text after Theorem 3.1:\n\n“The proof of Theorem 3.1 is postponed to Appendix N.1. Existing augmentation schemes approaches [2, 9, 44] do not offer any theoretical guarantees regarding the validity of the methodology. In contrast, Theorem 3.1 provides the first rigorous theoretical framework that establishes the soundness of a method within this class of data augmentation techniques”\n\n> Q3: The initial presentation of \\lambda_0 in section 2 is confusing, because it's not clear why or how \\lambda_0 depends on x. This becomes clear in 2.2, but I think it would be better to briefly comment on this when initially introducing \\lambda_0.\n\nWe thank the reviewer for this helpful suggestion. To address it, we have modified the text in Section 2:\n\n“Our goal is to model the hazard function $\\lambda$, i.e. the instantaneous event rate at time $t$ conditional on survival to $t$ and covariates $\\mathbf{x}$. We employ the sigmoid function $\\sigma(z) = 1 / (1 + \\exp(-z))$, which maps real-valued inputs to the interval $(0,1)$. The sigmoidal hazard model is constructed as the product of a normalized baseline hazard function ($\\lambda_0$) and a modulation function ($\\sigma$):\n$$ \\lambda(t\\mid\\mathbf{x}; \\phi, g(\\cdot;\\boldsymbol{\\theta})) := \\lambda_0(t, \\mathbf{x};\\phi)\\  \\sigma(g(t, \\mathbf{x};\\boldsymbol{\\theta})),$$\nwhere the normalized baseline hazard is given by\n$$ \\lambda_0(t, \\mathbf{x};\\phi) := \\frac{\\lambda_0(t;\\phi)}{Z(t, \\mathbf{x})},$$\n\nfor the baseline hazard $\\lambda_0: R_+\\to R_+$,  parametrized by $\\phi\\in R_{+}$, \nand a normalization factor $Z(t,\\mathbf{x})$ that depends on both time and covariates. The term $\\lambda_0(t, \\mathbf{x};\\phi)$ encodes our prior ``best‐guess'' hazard profile over time. The flexible function $g:\\mathbb{R}_+\\times\\mathbb{R}^p\\to \\mathbb{R}$, parametrized by $\\boldsymbol{\\theta}\\in\\mathbb{R}^{m}$, provides a data‐driven adjustment: once passed through the sigmoid, it multiplicatively *attenuates* the baseline hazard, continuously scaling it between zero and $\\lambda_0$. The normalization factor $Z(t,\\mathbf{x})$ ensures that the overall hazard remains properly scaled after modulation by the sigmoidal function (see Section 2.2 for details).\n”"}, {"review_id": "AlcYjaIV5G", "reviewer": "Reviewer_T13r", "summary": "The paper introduces NeuralSurv, a Bayesian deep learning framework for survival analysis that incorporates uncertainty quantification. The method leverages a two-stage data augmentation scheme involving Polya-Gamma variables and marked Poisson processes to enable efficient posterior estimation. The authors demonstrate superior calibration and competitive discriminative performance compared to state-of-the-art deep survival models on both synthetic and real-world datasets.", "questions": "1.\tThe authors should expand their discussion by referencing more literature on deep survival models, particularly focusing on key results and their limitations.  \n2.\tThe sigmoidal hazard model in (1) depends on the distribution of \\(\\theta\\). Since the random variable \\(Z(t,x)\\) in (5) is also influenced by the distribution of \\(\\theta\\), the same applies to model (1). Two critical questions arise: (a). How sensitive is the model to the prior distribution of \\(\\theta\\) when its dimensionality is high?  (b). Does the model (1), (5), and (6) converge to a fixed model as the dimension of \\(\\theta\\) approaches infinity?  \n3.\tWhat is the connection between the proposed model (1), (5), (6) and classical survival models such as the Cox proportional hazards model or the Accelerated Failure Time (AFT) model? Specifically, can the Cox or AFT model be derived as a special case of the proposed framework?  \n4.\tThe integral in the likelihood function involves only one variable, which is computationally tractable. However, the data augmentation strategy employed by the authors appears to be more computationally intensive. Could the authors provide a comparative analysis of the computational times for these two implementations?  \n5.\tCould the authors clarify which classes of models are not encompassed by their proposed framework (1), (5), (6)?", "limitations": "Yes.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "The authors provide a Bayesian uncertainty estimate for survival functions. However, it is unknown how the proposed method relies on the prior distributions.", "was_revised": true, "final_justification": "I'd like to keep the current score.", "rebuttal": "> Q1: The authors should expand their discussion by referencing more literature on deep survival models, particularly focusing on key results and their limitations. \n\nWe thank the reviewer for highlighting the need to expand our discussion of related work. We have added the following text to the conclusion to better position our model within the broader landscape of deep survival models. We have focused on the primary areas where our approach provides novelty.\n\n“In contrast to previous approaches in deep survival analysis, which are either constrained to discrete-time settings [43, 28, 13, 24, 26, 32] or lack the ability to provide Bayesian uncertainty quantification [43, 28, 13, 25, 24, 26, 32], \\emph{NeuralSurv} introduces a continuous-time modeling framework that naturally incorporates Bayesian inference, enabling both accurate survival predictions and well-calibrated uncertainty estimates.”\n\nAppendix  L remains available as a more detailed overview of existing literature on  machine learning methods in survival analysis.\n\n> Q2: Two critical questions arise: (a). How sensitive is the model to the prior distribution of (\\theta) when its dimensionality is high? (b). Does the model (1), (5), and (6) converge to a fixed model as the dimension of (\\theta) approaches infinity?\n\n\nWe thank the reviewer for these insightful questions. (a) NeuralSurv’s inference algorithm is built on the local linearization approach for Bayesian neural networks [18], which imposes an isotropic unit‐variance Gaussian prior on all network weights, this modeling choice is fixed and not subject to alternative specifications. (b) To date, no theoretical results have been published on the sensitivity of this method to high‐dimensional priors or its convergence as the number of parameters grows without bound. Accordingly, we cannot assess how models (1), (5), and (6) behave in the limit $dim(\\theta)\\to\\infty$, since any such asymptotic guarantees would require convergence proofs for the underlying BNN framework, which are not yet available.\n\n> Q3: What is the connection between the proposed model (1), (5), (6) and classical survival models such as the Cox proportional hazards model or the Accelerated Failure Time (AFT) model? Specifically, can the Cox or AFT model be derived as a special case of the proposed framework?\n\nWe thank the reviewer for requesting this clarification. Although the sigmoidal hazard function in equation (1) may appear structurally similar to both the Cox proportional hazards (PH) and AFT models, neither the Cox PH nor the AFT hazard functions arise as special cases of equation (1). In particular, no mathematical transformation renders the sigmoidal hazard function equivalent to either the Cox PH or AFT formulations. Nevertheless, we emphasize that the sigmoidal hazard function falls within a broader class of models that have been widely applied in practice (see answer to Question 5).\n\n> Q4: The integral in the likelihood function involves only one variable, which is computationally tractable. However, the data augmentation strategy employed by the authors appears to be more computationally intensive. Could the authors provide a comparative analysis of the computational times for these two implementations?\n\nWe thank the reviewer for this suggestion. While employing a numeric approximation (such as the trapezoidal or quadrature rule) to approximate the cumulative hazard integral does indeed make the likelihood evaluable, this approach presents two key limitations. First, the numerical error depends critically on the grid spacing and placement. If the hazard function exhibits rapid variation across time, a coarse grid may smooth or entirely miss these sharp changes, introducing systematic bias into parameter estimates and undermining both statistical efficiency and validity [A-D]. Second, although numerical integration makes the likelihood explicitly computable, directly using it as a loss in a neural network framework inherently leads to frequentist (MLE) estimation. To perform Bayesian inference via variational inference would require deriving an evidence lower bound (ELBO) that involves integrating over the variational parameter distributions [E]. This would result in a double integral problem. In contrast, our continuous-time framework avoids any double integrals and avoids any integrals over the parameters, enabling closed-form variational updates that are both computationally efficient and statistically principled. We are willing to include further clarification in the text in the data-augmentation section if the reviewer deems it necessary.\n\n[A] Michael Sloma, Fayeq Syed, Mohammedreza Nemati, Kevin S. Xu Proceedings of AAAI Spring Symposium on Survival Prediction - Algorithms, Challenges, and Applications 2021, PMLR 146:118-131, 2021.\n\n[B] Mert Ketenci, Shreyas Bhave, Noemie Elhadad, Adler Perotte Proceedings of the 8th Machine Learning for Healthcare Conference, PMLR 219:360-380, 2023.\n\n[C] Kvamme, H., Borgan, Ø. Continuous and discrete-time survival prediction with neural networks. Lifetime Data Anal 27, 710–736 (2021).\n\n[D] Voelkle, M. C., & Oud, J. H. L. (2012). Continuous time modelling with individually varying time intervals for oscillating and non‐oscillating processes. British Journal of Mathematical and Statistical Psychology, 66(1), 103–126. \n\n[E] Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. (2017). Variational Inference: A Review for Statisticians. Journal of the American Statistical Association, 112(518), 859–877. \n\n> Q5: Could the authors clarify which classes of models are not encompassed by their proposed framework (1), (5), (6)?\n\nSigmoidal hazard functions were first introduced by [B] and represent a distinct class of models. They have been used in the context of point process modeling and survival analysis. To help orient the reader and highlight existing literature on sigmoidal hazard functions, we have added the following text at the end of the introductory paragraph in Section 2:\n\n“Modeling hazard and intensity functions using sigmoidal transformations is common in applications such as survival analysis [10] and point process models [ 44 , 9, A, B, C]. This approach is popular due to the balance it offers between modeling flexibility and analytical tractability.”\n\nHowever, we would like to emphasize that, to the best of our knowledge, our work is the first to introduce a sigmoidal hazard parameterization in the context of deep survival analysis, where the flexible function g in equation (1) is modeled as the output of a Bayesian neural network (BNN). This formulation required a novel augmentation scheme to enable tractable inference under the Bayesian framework, which constitutes a key technical contribution of our work (see answer 1 to reviewer Pg1A).\n\n[A] Yee Teh and Vinayak Rao. Gaussian process modulated renewal processes. In Advances in Neural Information Processing Systems, 2011\n\n[B] Ryan Prescott Adams, Iain Murray, and David J. C. MacKay. Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities. ICML ’09, page 9–16. ACM, 2009.\n\n[C] Virginia Aglietti, Edwin V Bonilla, Theodoros Damoulas, and Sally Cripps. Structured variational inference in continuous cox process models. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019."}, {"review_id": "pwDCQ2Hexz", "reviewer": "Reviewer_52KX", "summary": "The paper introduces NeuralSurv, a deep survival model that offers an architecture-agnostic framework for continuous-time survival analysis that can capture complex time-varying covariate-risk relationships. The method includes a two-stage data-augmentation scheme utilizing latent marked Poisson processes and Pólya-Gamma variables. The BNN s approximated with local linearization and trained with a variational mean-field objective. The authors demonstrate good calibration and discriminative performance on synthetic and real-world datasets, particularly in data-scarce scenarios.", "questions": "- question about hyperparameters and initializations as elaborated in the weakness section. \n\n- is it possible to have a table to compare runtime with baseline methods as well? i understand that runtime is already discussed in discussion and appendix K, but quantitative comparison will be informative to researchers trying to using your method. (from my understanding NeuralSurv is already quite fast.)", "limitations": "yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths \n\n- The contribution is original and significant. this paper provides a first framework for neural network-based survival analysis, with very reasonable choices and construction that are well-placed in literature. This allows following works to study survival analysis leveraging NN's expressivity to model complex time-varying covariate-risk relationships, while having relatively reliable uncertainty quantification.\n\n- Strong empirical performance on the synthetic and real world datasets. \n\n- overall very well-written paper - all the sections are clear and coherent, and the method is well-explained in a step-by-step manner. although I am not an expert on variational inference, I am still able to follow the methodology and the math in the main text. I appreciate the primers in the appendix as well. \n\nWeaknesses \n\n- Lack of discussion on sensitivity to hyperparameters and initialization. the method relies on multiple steps of data augmentation and approximation for VI. my experience with BNNs is that that they are brittle to hyperparameter choices and initializations. From the experiments section of the paper, it seems like the experiments were performed once with a single architecture. No error bars on metrics from multiple runs (which i would expect from statisticians? maybe i'm wrong). It would be great to have a discussion on the robustness of the method.", "was_revised": true, "final_justification": "The authors have answered my questions and I maintain my positive recommendation.", "rebuttal": "> Q1: question about hyperparameters and initializations as elaborated in the weakness section\n\nWe thank the reviewer for their question. \n\n- Initialization: As described in Appendix K, in the EM algorithm, the parameters are initialized so that they match their prior expected values. Specifically, we set $\\boldsymbol{\\theta}^{(0)} = \\mathbf{0}$ and $\\phi^{(0)} = \\alpha_0 / \\beta_0$. In the CAVI algorithm, the hyperparameters are initialized so that the expected values of the model parameters match the MAP estimates. Specifically, we set $\\tilde{\\alpha}^{(0)} = \\phi_{\\text{MAP}} \\times \\tilde{\\beta}, \\quad \\text{and} \\quad (\\tilde{\\boldsymbol{\\mu}}, \\tilde{\\boldsymbol{\\Sigma}})^{(0)} = (\\boldsymbol{\\theta}_{\\text{MAP}}, \\mathbf{I}_m).$\n\n- Architecture: We chose not to optimize NeuralSurv (via hyperparameter searches) to every data set as we felt this would be an unfair representation of performance in comparison to the benchmarks.\n\n- Error bars: While we chose not to include the metrics variances in the paper in order to align with common practice in the survival analysis literature [24, 25, 32], where only point estimates are typically reported, we are happy to include them if the reviewer believes it would enhance clarity or interpretability.\n\n- Prior sensitivity: We have conducted a prior sensitivity analysis for the hyperparameters of the prior on phi. It is presented in Question 2 of reviewer Pg1A.\n\n> Q2: is it possible to have a table to compare runtime with baseline methods as well? i understand that runtime is already discussed in discussion and appendix K, but quantitative comparison will be informative to researchers trying to using your method. (from my understanding NeuralSurv is already quite fast.)\n\nWe thank the reviewer for this helpful suggestion.\nIn response, we have added a runtime analysis table reporting the inference times for NeuralSurv and several benchmark methods in the Appendix. For this experiment, we used the Colon dataset and explored two axes of variability:\n\n1/ Sample size: $N \\in \\{25, 125, 250\\}$\n\n2/ Neural network architecture:\n- MLP with 1 hidden layer (1L)\n- MLP with 2 hidden layers and 6 units (2L-6U)\n- MLP with 2 hidden layers 16 units (2L-16U)\n\nA summary of the results (runtime in minutes) is provided below for NeuralSurv, the fastest (DQS) and the slowest (DeepHit) benchmark methods. Please note that DQS is a new benchmark method required by reviewer  QbSK.\nAs expected, we observe that runtime increases with both the number of samples and the complexity of the model architecture. Notably, NeuralSurv, being a Bayesian method, is slower than frequentist baselines. This is consistent with prior literature, as Bayesian inference, particularly when combined with neural networks, is computationally more intensive. As illustrated by our experiments, this computational cost currently limits the scalability of NeuralSurv to very large datasets or highly complex models without additional approximation techniques. We hope this additional analysis clarifies the runtime characteristics of our method and provides a more complete picture of its computational trade-offs.\n\n | Method | N=25 ||| N=125 ||| N=250 |||\n   |--------|------|------|------|------|------|------|------|------|------|\n   |        | 1L   | 2L-6U | 2L-16U | 1L | 2L-6U | 2L-16U | 1L | 2L-6U | 2L-16U |\n| DeepHit | 0.164  | 0.139 | 0.140 | 0.120 | 0.126 | 0.142 | 0.122  | 0.145  | 0.135  |\n| DQS | 0.019 | 0.020 | 0.024 | 0.023 | 0.025 | 0.025 | 0.029 | 0.030 | 0.030 |\n| NeuralSurv | 0.621 | 1.165 | 0.566 | 3.926 | 16.454 | 22.673 | 7.373 | 87.098 |141.389 |"}, {"review_id": "rqm3Ur1Cyf", "reviewer": "Reviewer_Pg1A", "summary": "This work introduces NeuralSurv, a deep survival analysis framework with Bayesian uncertainty quantification. To do so, it proposes a two-stage data augmentation strategy using Pólya-Gamma variables and marked Poisson processes to enable exact continuous-time inference. A local linearization technique makes Bayesian neural networks tractable via closed-form coordinate-ascent updates. Experimental results highlight the performance of the method in discrimination and calibration in data-scarce regimes.", "questions": "1. Could the authors more clearly clarify which elements of the work are their own novel contribution versus drawn from other insights in the field of Bayesian inference, as well as adding a related work section providing background Bayesian inference and survival analysis methods? [I would raise my clarity score from \"2: fair\" to \"3: good\"; depending on the result of this discussion I may raise my Originality score from \"2: fair\" to \"3: good\"]\n2. Could the authors implement the following experiments, as discussed in Weaknesses? [I would raise my overall Rating from \"4: Borderline accept\" to \"5: Accept\" with these additional evaluations]\n    * Demonstrate the effect of different priors on concordance / IBS.\n    * Highlight survival curve uncertainty bounds in data-rich vs. data-scarce settings.\n    * Comparison against baseline methods of computing uncertainty intervals on survival curves.\n    * Include additional experimental results that use the whole dataset (rather than the subsampled component).\n3. Could the authors include confidence intervals on their evaluation statistics (c-index, Brier score), along with appropriate significance testing and multiple-comparisons correction? [this is another element that would contribute to me raising my overall Rating to a 5]", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 4, "originality": 2, "strengths_and_weaknesses": "I congratulate the authors on their work, and I am grateful for the opportunity to have reviewed this paper. The paper is well-written and a pleasure to read.\n\nStrengths:\n* The work is **very timely**: it takes steps to tackle an important problem in the in the survival analysis literature that I've noticed over recent years. As we obtain the means to statistically learn ever-more-granular survival curves (e.g., via expressive methods like SODEN [1] or mixture density networks [2]), there may or may not be support in the data for the granularity of survival curves that are being learned. As such, developing rigorous Bayesian mechanisms to credibly imbue individual survival curve estimates is an important problem. This work represents strong first steps in this direction -- e.g., the confidence intervals in Figure 1 are _exactly_ what a lot of practitioners are looking for in applied settings.\n* The experiments, as presented, are strong -- the method seems to outperform standard baselines on both concordance and calibration. I recommend the authors highlight these results not as evidence of utility (discussed more in Weaknesses), but rather as evidence that the model's design decisions (e.g., Weibull-like baseline hazard; isotropic Gaussian prior on $\\theta$; etc.) are not too restrictive for practical settings.\n* The NeuralSurv framework seems fast and scalable. By performing local linearization around the MAP and mean-field variational inference, the authors implement a closed-form coordinate ascent update with linear time complexity in model size.\n\nWeaknesses:\n* Related work and novelty: I admit that I do not have the deepest knowledge of Bayesian inference. With this in mind, the paper makes it rather difficult for me to tease apart which elements of NeuralSurv represent novel contributions, and which represent established methods in Bayesian inference. I would appreciate some clarity on this. My understanding is that the Polya-Gamma scheme in 3.1 is an established contribution; the Poisson Process Augmentation Scheme in 3.2 is novel; the Variational Mean-Field Approximation in 4.1 is an established method that is novel when applied to the survival setting; the local linearization of the BNN in 4.2 is established; and CAVI in 4.3 is established.\n  * I would appreciate if the authors added a Related Work section so that elements like this can be discussed, and so that their framework can be placed into conversation with other survival analysis frameworks.\n* There is a degree of misalignment between the experiments and the proposed contribution of the work. Whereas the experiments highlight concordance and calibration (via integrated Brier score), my understanding is that the core contribution of the work is not to implement the best ranking metric (or calibrated survival curve), but rather to (a) support the integration of prior knowledge into estimation, and (b) permit easy computation of a credible confidence interval. With that in mind, I would prefer that the authors include three additional experiments of the following form and to reframe the discussion of contributions accordingly:\n  1. Demonstrate how inference under different priors influences model performance as evaluated using classical statistics (c-index, IBS). This would inform the reader of the extent to which prior selection influences model performance on held-out data, as well as providing some guidance as to what may make for a \"good\" vs. \"bad\" prior in the survival setting.\n  2. Show several examples of survival curve estimation, both from the data-rich regime of the distribution (e.g., a patient who is similar to many other patients), and from the data-scarce regime (e.g., a patient who is dissimilar to many other patients). The aim of this experiment would be to demonstrate how the uncertainty in survival curves varies as a function of distributional support.\n  3. How do the credible intervals using this Bayesian setting compare to those computed via other methods? The authors may know the baselines here better than I do, but I'd imagine some baseline of bootstrapped confidence intervals and Monte Carlo dropout [3] would be important points of comparison.\n* I understand that the authors chose their experimental setup to highlight the performance of NeuralSurv in data-scarce regimes; but would it be possible to include additional experimental results on the performance of each method using the entire dataset (not subsampled to 100 train / 25 validation samples)? It would be good to get a sense of how the method performs as the dataset size increases. This would also allow for a more standardized comparison against the baseline methods, which are often evaluated on the entirety of COLON, METABRIC, and GBSG.\n4. The authors claim that their method supports \"time-varying covariate–risk relationships\" (Line 300); however, I don't see anything in their experimental framework supporting this. I see that $X=\\{x_i \\,:\\, i=1, ..., N\\}$ (Line 56); if the covariates were time-varying (e.g., as in [4]) I would have expected a different setup (e.g., $X=\\{x_i(t) \\,:\\, i=1, ..., N, t=0, ..., T\\}$). Could this be clarified?\n\n\n[1] Tang, Weijing, Jiaqi Ma, Qiaozhu Mei, and Ji Zhu. \"Soden: A scalable continuous-time survival model through ordinary differential equation networks.\" Journal of Machine Learning Research 23, no. 34 (2022): 1-29.\n\n[2] Han, Xintian, Mark Goldstein, and Rajesh Ranganath. \"Survival mixture density networks.\" In Machine Learning for Healthcare Conference, pp. 224-248. PMLR, 2022.\n\n[3] Gal, Yarin, and Zoubin Ghahramani. \"Dropout as a bayesian approximation: Representing model uncertainty in deep learning.\" In international conference on machine learning, pp. 1050-1059. PMLR, 2016.\n\n[4] Lee, Changhee, Jinsung Yoon, and Mihaela Van Der Schaar. \"Dynamic-deephit: A deep learning approach for dynamic survival analysis with competing risks based on longitudinal data.\" IEEE Transactions on Biomedical Engineering 67, no. 1 (2019): 122-133.", "was_revised": false, "final_justification": "", "rebuttal": "> Q1: the paper makes it rather difficult for me to tease apart which elements of NeuralSurv represent novel contributions, and which represent established methods in Bayesian inference. I would appreciate some clarity on this. \n\nWe thank the reviewer for highlighting the need to distinguish our novel contributions from established methods. All components of our algorithm have precedents in the literature: the Pólya–Gamma augmentation scheme was first introduced by Polson [A], and its combination with Poisson process augmentation appears in  [9,44]. The variational mean‑field approximation and CAVI are well‑known techniques that have previously been applied in survival settings (e.g., [B]). The local linearization method for BNNs, which yields a Gaussian process–based approximation, is also a recent but established approach [18]. Our key innovation lies in unifying these components into a single inference framework that fully exploits Gaussian process conjugacy alongside Poisson–Pólya–Gamma augmentation. To clarify this in the manuscript, we have added the following text in Section 6:\n\n“Pólya–Gamma and Poisson data-augmentation schemes (Section 3.3)  have been extensively employed with standard Gaussian process models [9,44]. Likewise, the local linearization of Bayesian neural networks, which yields a Gaussian process–based approximation, (Section 4.2), is a well established technique [18]. To our knowledge, this work is the first to integrate these two approaches into a unified framework that capitalizes on Gaussian process conjugacy. By combining these methods, we contribute a novel inference strategy at the intersection of Bayesian deep learning and Gaussian process modeling. ”\n\n[A] Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian Inference for Logistic Models Using Pólya–Gamma Latent Variables. Journal of the American Statistical Association, 108(504), 1339–1349.\n\n[B] Michael Komodromos, Eric O Aboagye, Marina Evangelou, Sarah Filippi, Kolyan Ray, Variational Bayes for high-dimensional proportional hazards models with applications within gene expression, Bioinformatics, Volume 38, Issue 16, August 2022, Pages 3918–3926.\n\n> Q2: Demonstrate how inference under different priors influences model performance as evaluated using classical statistics (c-index, IBS). \n\nWe thank the reviewer for their suggestion. \n\na/ Prior on theta: NeuralSurv’s inference algorithm is built on the local linearization approach for Bayesian neural networks [18], which imposes an isotropic unit‐variance Gaussian prior on all network weights, this modeling choice is fixed and not subject to alternative specifications. \n\nb/ Prior on phi: We conducted a prior sensitivity analysis by varying the hyperparameters of the Gamma prior over the precision parameter $\\phi$, where $\\phi \\sim \\text{Gamma}(\\alpha_0, \\beta_0)$. The main results reported in the paper use the default setting $\\alpha_0 = \\beta_0 = 1$, as specified in Appendix K. To assess sensitivity, we explored two additional configurations: (i) Reduced prior variance: $\\alpha_0 = \\beta_0 = 2$, (ii) Increased prior variance: $\\alpha_0 = \\beta_0 = 0.5$. These correspond to priors with half and double the prior variance, respectively, while maintaining the same mean, reflecting a non-informative prior structure.\n\nWe ran this analysis on all datasets. While the posterior distributions under different priors generally overlapped, their central tendencies varied, for example in the GBSG dataset, indicating mild to moderate sensitivity to prior choice. In contrast, the METABRIC results were nearly identical across priors, suggesting robustness in that setting. Overall, these findings suggest that while prior choice can influence results in certain datasets. When prior information is available, it should always be used to guide the model toward a principled balance between flexibility and regularization.\n\nFor each dataset, we will include in the Appendix tables showing the evaluation metrics, as well as the posterior median and 95% credible intervals.\n\nWe reported below the evaluation metrics, as well as the posterior median and 95% credible intervals for phi (only for METABRIC and GBSG datasets for space). \n\nMETABRIC dataset \nPrior| Posterior | C-index | IPCW IBS | \n|-|-|-|-|\nGamma(1, 1) | 0.685 [0.509, 0.897]| 0.584 | 0.212 | \nGamma(2, 2) | 0.670 [0.498, 0.879] | 0.536 | 0.201 | \nGamma(0.5, 0.5) | 0.675 [0.500, 0.886] | 0.553 | 0.200 |\n\nGBSG dataset \nPrior| Posterior | C-index | IPCW IBS | \n|-|-|-|-|\nGamma(1, 1) | 0.682 [0.522, 0.873] | 0.657  | 0.188 | \nGamma(2, 2) | 0.718 [0.554, 0.911] |  0.602 | 0.195 | \nGamma(0.5, 0.5) | 0.673 [0.514, 0.863] | 0.665 | 0.189 |\n\n> Q3: Show several examples of survival curve estimation, both from the data-rich regime of the distribution (e.g., a patient who is similar to many other patients), and from the data-scarce regime (e.g., a patient who is dissimilar to many other patients).\n\nTo identify similar and dissimilar patients, we computed the average pairwise distance of each patient's covariates to all other observations in the dataset. Patients with smaller average distances were those with the larger distributional support (in a data-rich regime), while those with larger distances were considered those with the smaller distributional support (in a data-scarce regime). We focused on two representative individuals: Patient A, who had the largest distributional support across the dataset, is the most similar and Patient B, who had the smallest distributional support across the dataset, is the most dissimilar.\n\nAs expected in a well-calibrated Bayesian framework, we observe consistently greater uncertainty in the data-scarce regime, reflecting limited covariate support for dissimilar patients.\n\nFor each dataset, we will include in the Appendix plots showing the posterior survival functions (median and 95% credible intervals) of these two patients side by side. \n\nTo summarize the results here, we report below the median and 95% interquartile range (IQR, defined as the difference between the 97.5th and 2.5th percentiles) of the posterior survival probabilities at four key timepoints: early (25% of the observational period), midpoint (50%), late (75%), and final (end of follow-up). Additionally, we report the mean IQR across the entire survival curve in the last column. We reported the results only for METABRIC and GBSG datasets for space. \n\n| GBSG dataset | Early | Mid-point | Late | Final | Mean of IQR | \n|-|-|-|-|-|-|\nMost dissimilar | 0.77 (0.14) |  0.62 (0.22)| 0.58 (0.24) | 0.56 (0.27) | 0.18 |\nMost similar | 0.70 (0.13) |  0.49 (0.19) | 0.34 (0.20) | 0.26 (0.20)| 0.15|\n\n| METABRIC dataset | Early | Mid-point | Late | Final | Mean of IQR | \n|-|-|-|-|-|-|\nMost dissimilar |  0.72 (0.17) |  0.53 (0.26) | 0.40 (0.29) |  0.32 (0.30)|  0.21|\nMost similar | 0.70 (0.15) | 0.49 (0.21) | 0.34 (0.23) | 0.26 (0.21)  | 0.17 |\n\n>  Q4: How do the credible intervals using this Bayesian setting compare to those computed via other methods? \n\nWe thank the reviewer for this insightful suggestion. We agree that directly comparing our Bayesian credible intervals with alternative uncertainty estimates, such as those obtained via Monte Carlo dropout [3], would meaningfully enhance our discussion of uncertainty quantification. Unfortunately, we were not able to include these experiments in the current submission due to limited computational resources (we prioritized the survival‐curve estimation experiments described above, as well as the runtime analyses presented below).\nNonetheless, we believe a rigorous comparison remains valuable. Such an analysis would involve implementing Monte Carlo dropout versions of each deep survival benchmark model, then evaluating and contrasting their uncertainty estimates with our credible intervals across both synthetic and real‐world datasets. However, carrying out this comprehensive study falls outside the present paper’s scope and exceeds the time available for revisions.\n\n> Q5:  include additional experimental results on the performance of each method using the entire dataset \n\nWe thank the reviewer for their thoughtful comment. We agree that evaluating model performance on the full datasets would provide a more standardized basis for comparison and insight into scalability. However, at this stage, our method is not yet scalable to full dataset sizes due to the computational demands of our fully Bayesian modeling approach. As noted in our response to Q2 to Reviewer 52KX, the running time becomes prohibitive as sample size increases, for reference, COLON (N = 1822), METABRIC (N = 1904), and GBSG (N = 2232) currently exceed our feasible limits.\nTo partially address this concern, we conducted an additional experiment on the full SUPPORT dataset in response to Reviewer QbSK’s Question 11. This experiment included comparisons with published benchmarks and demonstrated that our inference pipeline yields comparable C-index values, thereby providing some evidence for the validity.\n\n> Q6: The authors claim that their method supports \"time-varying covariate–risk relationships\" (Line 300); however, I don't see anything in their experimental framework supporting this. I see that  (Line 56); if the covariates were time-varying (e.g., as in [4]) I would have expected a different setup (e.g., ). Could this be clarified?\n\nWe thank the reviewer for pointing out the need for clarification. In our manuscript, the phrase “time‑varying covariate–risk relationships” was intended to mean “time‑varying relationships between covariates and risk,” not “relationships involving time‑varying covariates.” To address this, we have revised the conclusion to read:\n “We propose the first fully Bayesian framework for deep survival analysis that models time-varying relationships between covariates and risk.”."}]

Author rebuttal:

Meta-review: This work introduces a Bayesian neural network to model survival data, where the hazard function is expressed in a multiplicative form. It consists of a parametric baseline hazard multiplied by a risk function that varies with time, the covariate (X), and an (m)-dimensional parameter (\theta). The parametric baseline hazard function (\lambda_0(t, \phi)) does not involve the covariate (X) (see equations (5) and (1)) and is assumed to follow a Weibull distribution with a one-dimensional parameter (\phi) that has a Gamma prior. This results in a proportional hazards model, where the flexible function ($g(t, x, \theta)$) in the risk component is modeled using a Bayesian neural network (BNN) with an isotropic Gaussian prior on ($\theta$). As the authors claim, this appears to be the first deep survival model to incorporate Bayesian uncertainty quantification.

The paper proposes a novel algorithm that demonstrates superior performance on several synthetic and real-world datasets, consistently outperforming other deep learning methods. The paper is well-written and organized. However, one weakness is the scalability of the method, which seems limited to very small samples. This limitation contradicts the advantage of flexible modeling, which is particularly beneficial for larger datasets. Another weakness is the lack of sensitivity analysis regarding hyperparameters and initializations.

Overall, this is a significant contribution that addresses an important gap in Bayesian machine learning for survival data. I recommend its publication.
Decision: Accept (poster)