arXiv:2605.16208v1 [stat.ML] 15 May 2026
A Scalable Nonparametric Continuous-Time Survival Model through Numerical Quadrature
Chaeyeon Lee Department of Statistics Ewha Womans University Seoul, Korea
Sehwan Kim Department of Statistics Ewha Womans University Seoul, Korea
Hyungrok Do Department of Population Health NYU Grossman School of Medicine New York, NY, USA
Abstract Flexible continuous-time survival modeling is critical for capturing complex timevarying hazard dynamics in high-dimensional data; however, training such models remains challenging due to the intractable integral required for likelihood estimation. We introduce QSurv, a scalable deep learning framework that enables nonparametric continuous-time modeling without relying on time discretization or restrictive distributional assumptions. We propose a training objective based on Gauss-Legendre numerical quadrature, which approximates the cumulative hazard with high-order accuracy while facilitating efficient end-to-end training via standard backpropagation. Furthermore, to effectively capture non-stationary hazard dynamics in complex architectures, we introduce time-conditioned low-rank adaptation, a mechanism that conditions general neural backbones on time by dynamically modulating weights via low-rank updates. We provide theoretical analysis establishing approximation error bounds for cumulative-hazard evaluation. Comprehensive experiments across synthetic benchmarks, large-scale real-world tabular datasets, and high-dimensional medical imaging tasks demonstrate that QSurv achieves competitive predictive performance with advantages in instantaneous hazard function estimation, enabling more interpretable characterization of time-varying risk patterns.
1
Introduction
Survival analysis models the time to an event of interest and is fundamental to medical prognosis, engineering reliability, and financial risk stratification. Its primary challenge lies in handling censoring, where the event time is not fully observed for some subjects within the study period. In this work, we focus on right-censoring, which is the most prevalent form of censoring in practice. Classical methods such as the Kaplan-Meier (KM) estimator [29] and the Cox proportional hazards (PH) model [7] are widely used in survival analysis. The KM estimator serves as a nonparametric approach to estimating the population-level survival function, whereas the Cox PH model provides a semiparametric framework linking covariates to the hazard via a log-linear predictor. While these methods remain standard tools owing to their simplicity and interpretability, they are often insufficient to capture the complex, high-dimensional interactions between covariates and time-to-event. Although machine learning approaches like Random Survival Forests (RSF) [27] and classification-based frameworks such as survival stacking [8] offer improved risk prediction over classical methods, they remain limited in their ability to directly learn representations from unstructured and highdimensional data, such as electronic health records and medical imaging, and can impose a substantial computational burden at scale. To address these limitations, the field has increasingly turned to deep learning. One prominent line of work preserves the Cox framework while using a neural network to learn a nonlinear risk function (e.g., DeepSurv [30]). Another common direction adopts a discrete-time formulation, predicting Preprint.
hazards or event probabilities on a fixed time grid to enable flexible training objectives (e.g., DeepHit [33]). Alternatively, parametric deep models support likelihood-based training but constrain the hazard through a chosen distribution family (e.g., Deep Survival Machines [35]). Notably, these approaches typically rely on structural assumptions about the survival function or hazard. Cox-based neural models primarily learn a relative-risk score under proportional hazards, discrete-time models represent risk on a prespecified time grid, and parametric models restrict the hazard through a chosen distributional family. Thus, their limitations are not only computational or predictive, but also concern how directly they characterize the continuous-time instantaneous hazard. Direct characterization of the instantaneous hazard is clinically important because the hazard function is the most natural object for describing time-resolved risk. While the survival function summarizes accumulated risk up to time t, the instantaneous hazard h(t|x) describes how risk changes locally among individuals who remain at risk. Many clinical settings exhibit such time-local structure, including delayed treatment effects, transient early toxicity, crossing hazards, late recurrence, and decaying or accelerating risks over follow-up. These patterns are difficult to summarize by a single proportional-hazards effect and may be obscured when models focus only on survival probabilities or fixed-grid event probabilities. Prior work has therefore emphasized both the limitations of routine hazard-ratio reporting under non-proportional hazards [23, 41] and the practical value of examining hazard functions to understand disease dynamics and recurrence patterns [39, 24]. Despite this need, fully nonparametric deep learning for continuous-time hazard modeling remains underexplored because the survival likelihood requires both pointwise hazard evaluation and cumulativehazard integration. Directly parameterizing λθ (t|x) with a flexible time-conditioned network is therefore computationally challenging. In our setting, the integrand is the modeled hazard λθ (t|x), so the quadrature error is controlled by its temporal smoothness. Existing deep survival models often avoid this difficulty through proportional-hazards, discrete-time, or parametric assumptions. More flexible continuous-time approaches have recently been proposed, but SODEN [42] depends on adaptive ODE solvers, while DeSurv [10] models the cumulative distribution function rather than the hazard as the primitive quantity. Thus, scalable neural survival modeling that directly parameterizes the instantaneous hazard remains underdeveloped. Beyond the integration challenge, direct hazard modeling also raises the question of how time should enter the network. Since h(t|x) is a local-in-time risk function, the model must represent interactions between covariates x and follow-up time t, rather than treating time as a passive index. A standard approach is to concatenate t with x at the input, which is often adequate for tabular data processed by multilayer perceptrons (MLPs); however, for more expressive backbones such as ResNet [22], DenseNet [26] and Transformers [45, 14], simple input concatenation may be insufficient. A more effective strategy is to let temporal information modulate intermediate feature representations, enabling the model to capture time-varying risk more directly. To address the gap between scalable deep learning and direct hazard modeling, we introduce Quadrature-based Nonparametric Continuous-Time Deep Survival Analysis (QSurv). QSurv parameterizes the instantaneous hazard λθ (t|x) as a flexible function of time and covariates, and uses numerical quadrature to evaluate the cumulative hazard required by the continuous-time likelihood. In this way, QSurv makes the hazard function the primary modeled object while avoiding the computational burden of adaptive ODE solvers and the precision loss of time discretization. Our main contributions are summarized as follows: • We propose a nonparametric continuous-time deep survival model that directly parameterizes the instantaneous hazard as a flexible function of time and covariates, enabling time-resolved characterization of risk without fixed time discretization. • We introduce a likelihood objective based on numerical quadrature. This eliminates the need for ODE solvers, enabling efficient end-to-end training with standard backpropagation and minibatch stochastic gradient descent. • We move beyond simple input concatenation by proposing time-conditioned low-rank approximation (Time-LoRA), which conditions general neural network backbones on t by dynamically modulating network weights, resulting in richer temporal representations. • We perform comprehensive experiments across synthetic benchmarks, real-world tabular data, and high-dimensional medical imaging tasks. Our results confirm QSurv’s efficacy and feasibility compared to state-of-the-art baselines. 2
2
Related Works
Integrating deep learning approaches into survival analysis has led to a wide range of methodological advances. While these deep learning-based survival analysis methods can be categorized in several ways, to keep the discussion focused, we organize the literature into two broad settings: (1) discretetime models and (2) continuous-time models. For a comprehensive review, see [46] and the references therein. In discrete-time survival models, time is partitioned into a set of intervals, and the model outputs interval-specific hazards or a discrete probability distribution over discretized time points. Multi-Task Logistic Regression (MTLR) [47] casts survival estimation as a sequence of binary classification tasks, one for each time interval, using logistic regression. This framework has been extended to neural MTLR (N-MTLR) to capture nonlinear effects [16]. DeepHit [33] similarly models discrete-time event probabilities with a neural network and accommodates competing risks by incorporating an additional ranking-based loss. In contrast, Nnet-survival estimates the conditional hazard probabilities in discrete time points [17]. For continuous-time survival models, existing deep learning approaches can be broadly grouped into three directions. Cox models: Many methods start from the Cox proportional hazards (PH) framework and relax parts of the original formulation to improve flexibility. A widely used example is DeepSurv, which replaces the linear predictor in the Cox model with a neural network [30]. Several extensions have been proposed. CoxCC modifies the loss to mimic a case–control sampling scheme, and CoxTime allows the hazard to depend on time through a time-varying effect model [32]. Deep Cox Mixtures further extends DeepSurv by modeling heterogeneous subpopulations via a mixture formulation [36]. More generally, the Deep Extended Hazard model provides a flexible hazard formulation that generalizes the baseline hazard structure [48]. Parametric survival models: Another line of work adopts fully parametric specifications, where the survival time is modeled as a function of covariates and error term, typically assuming distributions such as Weibull or log-normal. Neural networks are then used to learn covariate-dependent parameters [38, 4, 1, 35]. On the other hand, Survival Mixture Density Networks model the event-time density using a Gaussian mixture and use a neural network to map covariates to the mixture parameters [21]. ODE-based models: More recently, ODE-based approaches have been proposed that model the evolution of the cumulative hazard (or related quantities) via an ordinary differential equation [42, 43]. This yields a principled continuous-time formulation capable of capturing complex temporal dynamics. A key limitation, however, is that these models can be computationally very expensive to train, and they also introduce additional hyperparameters related to ODE-solver. A related approach, DeSurv [10], formulates an ODE for the cumulative distribution function and discretizes it via Gauss–Legendre quadrature, but recovers the event density via time differentiation rather than directly modeling the hazard. The literature highlights key limitations that motivate our work. Existing methodologies often impose structural assumptions on the hazard or event-time distribution: discrete-time models depend on a chosen time grid, Cox models inherit restrictions from the PH formulation, and parametric models rely on specific distributional assumptions. While ODE-based formulations offer flexibility in modeling hazard functions, they incur significant computational costs, and the existing quadrature-based alternative bypasses the hazard altogether by modeling the cumulative distribution function. This leaves a clear gap: scalable, fully nonparametric deep learning that directly models the instantaneous hazard remains underexplored. Yet the instantaneous hazard is itself a clinically meaningful quantity, providing a direct measure of risk at any given time that informs time-sensitive decisions such as patient monitoring and intervention timing.
3
Problem Formulation
Let T, C ∈ R+ denote the random variables for the true time-to-event and censoring time, respectively. We consider a setting with covariates X ∈ X ⊆ Rd . Under the assumption of covariate-dependent censoring, we assume that the event time and censoring time are conditionally independent given the covariates, i.e., T ⊥ C|X. In the observational setting, the true event time T is not always observed 3
directly due to censoring. Instead, we observe the pair (O, ∆), where O = min(T, C) denotes the observed time and ∆ = I(T ≤ C) is the binary event indicator. Given a dataset of n i.i.d. realizations D = {(xi , oi , δi )}ni=1 , our primary objective is to learn the parameters θ of the conditional hazard function λθ (t|x). More broadly, however, survival analysis aims to characterize the full time-to-event distribution, which may involve estimating the survival function Sθ (t|x), the cumulative hazard Λθ (t|x), or discrete-time event probabilities. The hazard function λθ (t|x) fully characterizes the time-to-event distribution. It is linked to the cumulative hazard Λθ (t|x) and the survival function Sθ (t|x) through the following equations: Z t Λθ (t|x) = λθ (u|x)du, Sθ (t|x) = exp (−Λθ (t|x)) . (1) 0
To estimate the model parameters θ, we minimize the negative log-likelihood (NLL). Under the standard assumption of covariate-dependent censoring (T ⊥ C|X), the terms governing the censoring distribution are independent of θ and factor out of the optimization. This yields the standard survival analysis objective function: n
L(θ) = −
1X [δi log λθ (oi |xi ) − Λθ (oi |xi )] . n i=1
(2)
Minimizing this objective with respect to θ requires a modeling approach that is both expressive enough to capture complex temporal dynamics and computationally tractable for evaluating the cumulative hazard integral. Existing literature navigates this trade-off through various approximations. This objective exposes the central computational difficulty in direct neural network parameterized hazard modeling: the event contribution requires evaluating the hazard locally at oi , whereas the survival contribution requires integrating the same learned hazard over the full interval [0, oi ]. QSurv is designed to preserve the hazard as the primitive modeled quantity while making this cumulativehazard evaluation efficient and differentiable. 3.1
Quadrature-based Nonparametric Continuous-Time Deep Survival Analysis (QSurv)
We model the log-hazard with a network fθ : X × R+ → R and set λθ (t|x) = exp{fθ (x, t)}. Minimizing the negative log-likelihood requires evaluating the cumulative hazard Λθ (oi |xi ) = R oi λ (s|x ) ds. While the log-hazard requires only a single forward pass, the cumulative hazard θ i 0 is analytically intractable due to the lack of parametric structure. To resolve this, we approximate Λθ (t|x) using K-point Gauss-Legendre quadrature. By applying the change of variables s = tτ to map the integration domain [0, t] to the canonical interval [0, 1], we obtain the approximation: K
Λ̂θ (t|x) =
tX wk λθ (tτk |x), 2
(3)
k=1
where τk = (ξk + 1)/2 are the Gauss-Legendre nodes shifted from the standard interval [−1, 1] to [0, 1] viaPthe affine map, and {wk }K k=1 are the corresponding standard Gauss-Legendre weights satisfying k wk = 2. The factor t/2 accounts for the Jacobian of the change of variables s = tτ combined with the [−1, 1] → [0, 1] affine transformation. Unlike standard fixed-grid approximations, such as Riemann sums or the trapezoidal rule, which are constrained to equidistant points, Gauss-Legendre quadrature optimizes node locations to maximize precision. The standard nodes are derived from the roots of the K-th degree Legendre polynomial PK (ξ) and are mapped to our target domain [0, 1] via the affine transformation τk = (ξk + 1)/2. The resulting rule is exact for polynomials of degree up to 2K − 1, and therefore achieves high accuracy for sufficiently smooth hazards with relatively few evaluations. Consequently, the approximation error is small for smooth hazard functions even with a small K. Crucially, this approximation is implemented as a vectorized set of K network evaluations that remains fully differentiable with respect to θ, enabling end-to-end training via standard backpropagation with a time complexity of O(K) per cumulative hazard evaluation. Importantly, compared to these Riemann or trapezoidal approximations, Gauss-Legendre quadrature is not tied to a dense time discretization. While grid-based methods require a large number of points to minimize discretization error, scaling computational cost linearly with the desired resolution, GaussLegendre leverages optimal node placement to achieve superior accuracy with a fixed, low budget of 4
Algorithm 1 Training procedure for QSurv K 1: Require: Dataset D = {(xi , oi , δi )}n and weights i=1 , Gauss-Legendre nodes τ ∈ [0, 1] K w∈R . 2: Initialize network parameters θ. 3: while not converged do 4: Sample minibatch B ⊂ D of size B. 5: for all (xi , oi , δi ) ∈ B in parallel do 6: Compute log-hazard by feed-forward: fθ (xi , oi ) 7: Estimate cumulative hazard (vectorized): Λ̂i ← o2i · w⊤ exp fθ (xi , oi τ ) . 8: end for h i P 1 9: Compute loss: L(θ; B) ← − |B| (xi ,oi ,δi )∈B δi fθ (xi , oi ) − Λ̂i . 10: Update parameters: θ ← θ − η∇θ L(θ; B) 11: end while
K evaluations. The parameters of QSurv are learned by minimizing the negative log-likelihood of the observed data. Given a dataset the objective function is defined as: n K i 1 Xh oi X wk exp fθ (xi , oi τk ) , L(θ; D) = − δi fθ (xi , oi ) − n i=1 2
(4)
k=1
where the first term captures the instantaneous risk for observed events, and the second term, approximated via the quadrature method described above. Algorithm 1 describes the training procedure for QSurv. 3.2
Approximation Error Bound
The validity of QSurv relies on the convergence of the quadrature approximation to the exact cumulative-hazard integral induced by the modeled hazard. This follows from classical properties of Gauss-Legendre quadrature: the K-point rule is exact for polynomials of degree up to 2K − 1 and admits a high-order error bound for sufficiently smooth integrands [18, 11]. In our setting, the integrand is the neural network parameterized hazard λθ (t|x), so the quadrature error is controlled by the temporal smoothness of the learned hazard. We state the resulting bound below and defer the proof, which is a direct application of the standard Gauss-Legendre error formula after an affine transformation from [−1, 1] to [0, t], to Appendix B. Theorem 3.1 (Approximation Error Bound of Cumulative Hazard Function). Let the hazard function λθ (·|x) be 2K-times continuously differentiable on [0, tmax ]. Let Λθ (t|x) denote the exact cumulative hazard induced by λθ and Λ̂θ,K (t|x) be the approximation computed using K-point Gauss-Legendre quadrature. Then, for any t ∈ [0, tmax ], the approximation error is bounded by: εK (t|x) = Λθ (t|x) − Λ̂θ,K (t|x) ≤
t2K+1 (K!)4 (2K) max λ (τ |x) . (2K + 1)[(2K)!]3 τ ∈[0,t] θ
(5)
Furthermore, if λθ (·|x) is a polynomial of degree at most 2K − 1 with respect to time, the approximation is exact, i.e., Λ̂θ,K (t|x) = Λθ (t|x). Theorem 3.1 shows that the quadrature error is controlled by the number of nodes K and the temporal smoothness of the learned hazard. For sufficiently smooth hazards, the factorial term in the bound yields rapid error decay as K increases, and the Gauss–Legendre rule is exact for hazards that are polynomial in time of degree at most 2K − 1. Thus, a moderate fixed value of K can provide accurate cumulative-hazard estimates in practice, as confirmed empirically in Appendix C. The result also motivates using smooth hazard parameterizations, such as networks with Softplus or GELU activations, since high-order quadrature accuracy relies on the existence of temporal derivatives. Overall, K should be viewed as a numerical precision parameter rather than a model-complexity or regularization parameter. We establish a direct link between the numerical precision of the integration scheme and the training objective. Specifically, we show that the error in the negative log-likelihood is bounded by the error in the cumulative hazard approximation. 5
Corollary 3.2 (Propagation of Approximation Error). Let L(θ; t, δ, x) = −δ log λθ (t|x) + Λθ (t|x) be the true negative log-likelihood for a sample with time t, event indicator δ, and features x. Let L̂K denote the approximate loss computed using the quadrature estimate Λ̂θ,K . Then the loss approximation error equals the cumulative hazard approximation error: L(θ; t, δ, x) − L̂K (θ; t, δ, x) = εK (t|x).
(6)
This result relates the stability of the learning process to the integration accuracy. Because the cumulative hazard enters the negative log-likelihood linearly, the approximation error is additive and does not compound through non-linear transformations. Consequently, the bound εK (t|x) translates directly to the loss function. When combined with the convergence properties established in Theorem 3.1, this suggests that controlling the quadrature error effectively limits the deviation of the training objective from the true likelihood, allowing the approximate loss to serve as a reliable proxy during optimization. 3.3
Choice of the Number of Quadrature Nodes
The number of quadrature nodes K controls the numerical accuracy of the cumulative-hazard approximation and the computational cost of training. Importantly, K is not a time-discretization parameter in the same sense as the number of bins in discrete-time survival models. QSurv still represents the hazard as a continuous function of time; K only determines how accurately the integral of this hazard is approximated inside the likelihood. Very small values of K may provide an inaccurate approximation when the hazard varies rapidly over time, which can distort the training objective. Increasing K improves the fidelity of the quadrature approximation but requires additional evaluations of the neural network parameterized hazard at the quadrature nodes. Thus, K should be interpreted as a numerical precision parameter rather than a regularization hyperparameter: reducing K does not impose a simpler hazard model, but instead coarsens the approximation to the cumulative hazard. In practice, we use a moderate fixed value of K and examine sensitivity to this choice empirically in Section C. 3.4
Time-Conditioned LoRA
For tabular data, a common strategy to model time-varying hazards is to treat time t as an additional input feature, simply concatenating it with the static covariates x, e.g., as in CoxTime [32], DeSurv [10], and SODEN [42]. However, this approach does not trivially extend to high-dimensional data such as medical images processed by backbones like ResNet [22], DenseNet [26] and Transformers [45, 14]. Concatenating a scalar t directly to an image input is not naturally compatible with convolutional or transformer backbones. To address this, we introduce time-conditioned LoRA, which uses a time-dependent weight matrix at the network’s penultimate layer to make the hazard prediction adaptive in t without re-running the backbone. This avoids K redundant forward passes of the heavy backbone, substantially reducing the computational cost compared with fully timeentangled architectures. For this penultimate layer with weight W ∈ Rdout ×din and bias b ∈ Rdout , we parameterize the time-varying weight W (t) as a low-rank perturbation of the static base [25]: W (t) = W + ∆W (t) = W + U · diag(s(t)) · V,
(7)
where U ∈ Rdout ×r and V ∈ Rr×din are learnable low-rank projection matrices (r ≪ min(din , dout )). The temporal dynamics are captured by the modulation coefficients s(t) = gψ (ϕ(t)) ∈ Rr , which are generated by a lightweight MLP gψ acting on a time embedding ϕ(t) [31]. The forward pass can be re-interpreted as a time-gated residual branch operating within a learned subspace: z(t) = W h + U s(t) ⊙ (V h) + b. (8) Here, V projects the input h into a latent subspace, s(t) applies a time-varying activation, and U projects the result back to the output dimension; a final time-independent linear layer then maps z(t) to the scalar log-hazard fθ (x, t). This structural decoupling yields a significant computational advantage during numerical integration, where the hazard must be evaluated at multiple quadrature nodes {τk }K k=1 for a single subject. Because Time-LoRA is applied only at the penultimate layer, 6
Table 1: Comparison of L1 error of instantaneous hazard (smaller the better) across different parametric distributions with standard deviations over 20 random seeds. Best results are bold, second best are underlined. Full results are presented in Table 4 in the appendix. Distribution
CoxCC
CoxTime
NnetSurv
MDN
DeSurv
SODEN
QSurv
Exponential Weibull Gamma Gompertz Log-normal Log-logistic
0.1775 ± 0.0410 0.1576 ± 0.0322 0.9012 ± 0.0921 0.1713 ± 0.0262 0.1235 ± 0.0248 0.1621 ± 0.0210
0.1775 ± 0.0422 0.1182 ± 0.0207 0.5661 ± 0.1769 0.1890 ± 0.0549 0.1114 ± 0.0294 0.1266 ± 0.0242
0.3830 ± 0.0438 0.2383 ± 0.0254 0.7249 ± 0.0370 0.2921 ± 0.0276 0.1938 ± 0.0219 0.3290 ± 0.0525
0.0568 ± 0.0446 0.0644 ± 0.0102 0.5646 ± 0.0459 0.0920 ± 0.0140 0.1882 ± 0.0520 0.1031 ± 0.0090
0.0448 ± 0.0159 0.0295 ± 0.0090 0.2930 ± 0.0907 0.0502 ± 0.0278 0.0473 ± 0.0201 0.0513 ± 0.0299
0.0330 ± 0.0138 0.0272 ± 0.0119 0.1716 ± 0.0901 0.0262 ± 0.0174 0.0462 ± 0.0183 0.0517 ± 0.0295
0.0352 ± 0.0164 0.0282 ± 0.0139 0.1718 ± 0.0859 0.0258 ± 0.0180 0.0457 ± 0.0185 0.0491 ± 0.0291
the embedding h is independent of t, so the base computation W h and the initial projection V h depend solely on the static input x and are computed only once and cached. For each quadrature node τk , only the lightweight temporal branch is evaluated to obtain s(oi τk ), followed by the low-rank recombination. This avoids K redundant forward passes of the heavy backbone, drastically reducing the inference cost compared to fully time-entangled architectures. This caching is only valid because the modulation is confined to the head: applying Time-LoRA to an intermediate backbone layer would make the inputs to all downstream layers time-dependent and require a full backbone pass per quadrature node, which is why we restrict it to the penultimate layer. Unlike FiLM, which applies element-wise affine modulation, Time-LoRA allows time-dependent low-rank mixing across feature channels, providing a more expressive head for modeling time-varying hazards [15].
4
Experiments
We performed a series of experiments across a comprehensive set of benchmark datasets to evaluate the proposed methods. To ensure a rigorous comparison across different modeling paradigms, we evaluated against Cox-style deep learning models CoxCC and CoxTime [32], discrete-time models DeepHit [33] and NnetSurv [17], and continuous-time estimators including Survival MDN [21], DeSurv [10], and SODEN [42]. For a controlled comparison of time-conditioned continuous-time architectures, CoxTime, DeSurv, and SODEN were implemented with the same time-conditioned LoRA mechanism used in QSurv. Details of the deep survival models are provided in Appendix F.6. Code for reproducing the models and experiments is available at https://github.com/hyungrok-do/qsurv. 4.1
Simulation Study
We conducted a simulation study to assess whether QSurv can recover instantaneous hazard functions under nonlinear covariate effects. Details of the data-generating mechanisms and additional results are provided in Appendix D. Table 1 reports L1 error on the test set for the estimated instantaneous hazard functions. QSurv and SODEN consistently produced the smallest hazard-estimation errors, particularly for distributions with more complex hazard shapes. The main difference was computational: QSurv achieved similar accuracy with much shorter training time. Across simulation settings, SODEN required approximately 60 seconds on average, compared with approximately 6 seconds for QSurv. Figures 6–11 show the corresponding test-set estimates of the instantaneous hazard, cumulative hazard, and survival functions. 4.2
Real-world Datasets
We evaluated QSurv on nine benchmark datasets, comprising six tabular datasets (FLCHAIN, FRAMINGHAM, GBSG, METABRIC, NWTCO, SUPPORT2) and three radiology imaging datasets (COVID-19-NY, C4KC-KiTS, and BraTS). These datasets cover diverse medical domains and vary substantially in sample size, feature dimensionality, and censoring rates; their main characteristics are summarized in Appendix F. Model performance was evaluated using the time-dependent concordance index (Ctd ), Integrated Brier Score (IBS), Integrated Binomial Log-Likelihood (IBLL), and Dcalibration. Higher Ctd and IBLL and lower IBS indicate better performance. We define these metrics in Appendix F.3. To assess performance across follow-up durations, we evaluated each method over three horizons: the full horizon defined by Ĝ(τ ) = 0.001, and two restricted quantile-based horizons, TQ1 and TQ2 , with TQ2 corresponding to the median horizon. For robust comparison, we performed random-search hyperparameter optimization with 30 configurations per model for tabular datasets and 20 configurations per model for imaging datasets. The full procedure was repeated 20 times for tabular datasets and 5 times for imaging datasets. Due to the high computational burden of adaptive 7
Table 2: Performance of QSurv vs deep survival baselines across the full and quantile-based horizons. Best results are bold, second best are underlined. The right-most D-cal column is the mean p-value of the discrete-time D-calibration goodness-of-fit test (higher is better; p > 0.05 indicates no calibration violation). Horizon: Full (Ĝ(τ )≈0.001)
Horizon: TQ1
Horizon: TQ2 (median)
Dataset
Model
Ctd ↑
IBS ↓
IBLL ↑
Ctd ↑
IBS ↓
IBLL ↑
Ctd ↑
IBS ↓
IBLL ↑
D-cal ↑
COVID-19-NY
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv SODEN QSurv
0.5693 ± 0.0961 0.5587 ± 0.0847 0.6677 ± 0.0793 0.6660 ± 0.0844 0.6889 ± 0.0626 0.7420 ± 0.0923 0.7143 ± 0.0530 0.7456 ± 0.1046
0.1885 ± 0.0882 0.2079 ± 0.1211 0.0708 ± 0.0209 0.0797 ± 0.0159 0.1892 ± 0.0797 0.0617 ± 0.0269 0.1000 ± 0.0282 0.0716 ± 0.0155
-0.5211 ± 0.2118 -0.6197 ± 0.3588 -0.2426 ± 0.0600 -0.2629 ± 0.0406 -0.5626 ± 0.2426 -0.1981 ± 0.0802 -0.2967 ± 0.0709 -0.2318 ± 0.0489
0.5925 ± 0.1640 0.5901 ± 0.1031 0.4951 ± 0.1224 0.5753 ± 0.0990 0.5669 ± 0.1987 0.7173 ± 0.2020 0.7356 ± 0.1669 0.7260 ± 0.1002
0.0048 ± 0.0052 0.0044 ± 0.0046 0.0045 ± 0.0046 0.0045 ± 0.0046 0.0185 ± 0.0100 0.0052 ± 0.0049 0.0048 ± 0.0049 0.0062 ± 0.0050
-0.0309 ± 0.0349 -0.0306 ± 0.0337 -0.0291 ± 0.0282 -0.0265 ± 0.0259 -0.0954 ± 0.0358 -0.0367 ± 0.0258 -0.0335 ± 0.0275 -0.0405 ± 0.0259
0.5504 ± 0.0633 0.5636 ± 0.0891 0.6516 ± 0.0859 0.6793 ± 0.1175 0.7116 ± 0.0536 0.7944 ± 0.0615 0.7689 ± 0.1099 0.7951 ± 0.0908
0.0366 ± 0.0090 0.0334 ± 0.0100 0.0371 ± 0.0130 0.0344 ± 0.0069 0.0406 ± 0.0107 0.0322 ± 0.0089 0.0324 ± 0.0098 0.0352 ± 0.0157
-0.1656 ± 0.0339 -0.1592 ± 0.0327 -0.1660 ± 0.0491 -0.1463 ± 0.0171 -0.1660 ± 0.0323 -0.1410 ± 0.0334 -0.1360 ± 0.0357 -0.1441 ± 0.0488
0.7264 ± 0.4338 0.6061 ± 0.5335 0.7470 ± 0.3440 0.8787 ± 0.2473 0.7658 ± 0.2906 0.8502 ± 0.1983 0.9221 ± 0.1042 0.7391 ± 0.2436
C4KC-KiTS
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv SODEN QSurv
0.4116 ± 0.1798 0.5808 ± 0.1985 0.5090 ± 0.1089 0.5422 ± 0.1375 0.5230 ± 0.2206 0.5740 ± 0.1630 0.4987 ± 0.1159 0.6381 ± 0.1446
0.1099 ± 0.0050 0.1086 ± 0.0015 0.2179 ± 0.1236 0.1262 ± 0.0240 0.1106 ± 0.0139 0.1156 ± 0.0152 0.1074 ± 0.0045 0.1260 ± 0.0453
-0.3644 ± 0.0240 -0.3571 ± 0.0139 -0.7631 ± 0.3839 -0.5230 ± 0.1907 -0.4742 ± 0.2106 -0.6108 ± 0.3646 -0.3558 ± 0.0307 -0.4002 ± 0.1113
0.4834 ± 0.2542 0.5877 ± 0.2812 0.5290 ± 0.0644 0.7345 ± 0.1797 0.5290 ± 0.1920 0.6274 ± 0.2445 0.5679 ± 0.1812 0.7404 ± 0.1547
0.0307 ± 0.0196 0.0304 ± 0.0194 0.0310 ± 0.0195 0.0340 ± 0.0172 0.0352 ± 0.0186 0.0329 ± 0.0232 0.0307 ± 0.0196 0.0288 ± 0.0178
-0.1523 ± 0.0805 -0.1479 ± 0.0788 -0.1570 ± 0.0838 -0.1964 ± 0.1066 -0.1908 ± 0.1246 -0.2136 ± 0.1838 -0.1401 ± 0.0765 -0.1220 ± 0.0679
0.4030 ± 0.1679 0.5806 ± 0.1851 0.4900 ± 0.1264 0.5371 ± 0.1172 0.5270 ± 0.2244 0.5685 ± 0.1636 0.5150 ± 0.0939 0.6207 ± 0.1596
0.0740 ± 0.0132 0.0734 ± 0.0124 0.0820 ± 0.0192 0.0867 ± 0.0277 0.0773 ± 0.0155 0.0826 ± 0.0228 0.0754 ± 0.0131 0.0776 ± 0.0235
-0.2840 ± 0.0468 -0.2787 ± 0.0408 -0.3271 ± 0.0690 -0.4008 ± 0.1529 -0.3829 ± 0.2099 -0.4910 ± 0.2947 -0.2879 ± 0.0463 -0.2910 ± 0.0644
1.0000 ± 0.0001 0.9999 ± 0.0001 0.6826 ± 0.4021 0.9574 ± 0.0529 0.9693 ± 0.0473 0.9889 ± 0.0225 0.9995 ± 0.0004 0.9845 ± 0.0209
BraTS
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv SODEN QSurv
0.4925 ± 0.0344 0.4541 ± 0.0359 0.5196 ± 0.0383 0.5355 ± 0.0395 0.5333 ± 0.0537 0.5459 ± 0.0297 0.5672 ± 0.0510 0.5811 ± 0.0387
0.2021 ± 0.1203 0.1171 ± 0.0130 0.1164 ± 0.0175 0.1243 ± 0.0394 0.1274 ± 0.0103 0.1038 ± 0.0176 0.1089 ± 0.0121 0.1038 ± 0.0133
-0.6222 ± 0.3545 -0.3764 ± 0.0363 -0.3765 ± 0.0453 -0.5086 ± 0.3472 -0.4131 ± 0.0252 -0.3405 ± 0.0707 -0.3489 ± 0.0442 -0.3327 ± 0.0371
0.4976 ± 0.0676 0.4414 ± 0.0549 0.5020 ± 0.0679 0.5371 ± 0.0533 0.5321 ± 0.0342 0.5570 ± 0.0372 0.5984 ± 0.0894 0.6149 ± 0.0559
0.1414 ± 0.0227 0.1355 ± 0.0128 0.1357 ± 0.0125 0.1806 ± 0.1093 0.1754 ± 0.0082 0.1355 ± 0.0141 0.1511 ± 0.0308 0.1380 ± 0.0254
-0.5501 ± 0.2447 -0.4456 ± 0.0388 -0.4582 ± 0.0378 -0.5546 ± 0.2757 -0.5387 ± 0.0175 -0.4378 ± 0.0363 -0.4757 ± 0.0717 -0.4393 ± 0.0628
0.4949 ± 0.0382 0.4460 ± 0.0406 0.5130 ± 0.0457 0.5471 ± 0.0465 0.5429 ± 0.0479 0.5438 ± 0.0297 0.5763 ± 0.0607 0.5968 ± 0.0537
0.2168 ± 0.0411 0.2021 ± 0.0127 0.2010 ± 0.0160 0.2599 ± 0.1342 0.2120 ± 0.0088 0.2044 ± 0.0195 0.2217 ± 0.0373 0.1993 ± 0.0329
-0.7373 ± 0.3262 -0.5939 ± 0.0319 -0.6005 ± 0.0461 -0.8321 ± 0.5395 -0.6149 ± 0.0186 -0.5930 ± 0.0451 -0.6358 ± 0.0900 -0.5868 ± 0.0868
0.0702 ± 0.1141 0.3090 ± 0.3068 0.1414 ± 0.1793 0.3879 ± 0.5029 0.0000 ± 0.0000 0.1535 ± 0.1623 0.0951 ± 0.1047 0.2161 ± 0.2395
FLCHAIN
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv QSurv
0.7967 ± 0.0101 0.7968 ± 0.0101 0.7883 ± 0.0107 0.7887 ± 0.0113 0.7954 ± 0.0096 0.7968 ± 0.0097 0.7968 ± 0.0096
0.0936 ± 0.0024 0.0936 ± 0.0024 0.1056 ± 0.0026 0.0935 ± 0.0028 0.0946 ± 0.0020 0.0936 ± 0.0021 0.0935 ± 0.0023
-0.3103 ± 0.0075 -0.3100 ± 0.0074 -0.3400 ± 0.0069 -0.3103 ± 0.0086 -0.3142 ± 0.0068 -0.3106 ± 0.0067 -0.3097 ± 0.0070
0.8024 ± 0.0120 0.8027 ± 0.0125 0.7891 ± 0.0145 0.7898 ± 0.0152 0.8012 ± 0.0117 0.8030 ± 0.0116 0.8031 ± 0.0118
0.0680 ± 0.0023 0.0680 ± 0.0023 0.0755 ± 0.0029 0.0681 ± 0.0026 0.0689 ± 0.0021 0.0680 ± 0.0023 0.0678 ± 0.0023
-0.2395 ± 0.0075 -0.2389 ± 0.0076 -0.2587 ± 0.0081 -0.2402 ± 0.0082 -0.2425 ± 0.0069 -0.2398 ± 0.0073 -0.2387 ± 0.0077
0.7988 ± 0.0099 0.7989 ± 0.0102 0.7896 ± 0.0114 0.7903 ± 0.0121 0.7973 ± 0.0097 0.7988 ± 0.0097 0.7989 ± 0.0094
0.0861 ± 0.0023 0.0862 ± 0.0023 0.0974 ± 0.0028 0.0864 ± 0.0026 0.0872 ± 0.0020 0.0862 ± 0.0022 0.0860 ± 0.0022
-0.2914 ± 0.0076 -0.2910 ± 0.0074 -0.3190 ± 0.0076 -0.2922 ± 0.0084 -0.2949 ± 0.0071 -0.2919 ± 0.0070 -0.2907 ± 0.0073
0.9224 ± 0.1189 0.9354 ± 0.1325 0.0192 ± 0.0482 0.7866 ± 0.2870 0.4600 ± 0.2720 0.8706 ± 0.1812 0.8819 ± 0.1852
FRAMINGHAM
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv QSurv
0.7407 ± 0.0142 0.7411 ± 0.0156 0.7262 ± 0.0136 0.7316 ± 0.0165 0.7378 ± 0.0128 0.7400 ± 0.0138 0.7414 ± 0.0146
0.0882 ± 0.0029 0.0881 ± 0.0031 0.1780 ± 0.0162 0.0902 ± 0.0033 0.0890 ± 0.0025 0.0887 ± 0.0030 0.0884 ± 0.0029
-0.2929 ± 0.0088 -0.2925 ± 0.0095 -0.5220 ± 0.0412 -0.2997 ± 0.0110 -0.2960 ± 0.0084 -0.2954 ± 0.0104 -0.2933 ± 0.0087
0.7556 ± 0.0131 0.7561 ± 0.0146 0.7403 ± 0.0136 0.7439 ± 0.0159 0.7541 ± 0.0125 0.7565 ± 0.0137 0.7578 ± 0.0145
0.0607 ± 0.0032 0.0608 ± 0.0034 0.0819 ± 0.0089 0.0615 ± 0.0031 0.0611 ± 0.0031 0.0610 ± 0.0035 0.0607 ± 0.0033
-0.2162 ± 0.0099 -0.2160 ± 0.0109 -0.2948 ± 0.0254 -0.2200 ± 0.0103 -0.2183 ± 0.0108 -0.2178 ± 0.0116 -0.2161 ± 0.0105
0.7407 ± 0.0142 0.7411 ± 0.0156 0.7262 ± 0.0136 0.7316 ± 0.0165 0.7378 ± 0.0128 0.7400 ± 0.0138 0.7414 ± 0.0146
0.0882 ± 0.0029 0.0881 ± 0.0031 0.1780 ± 0.0162 0.0902 ± 0.0033 0.0890 ± 0.0025 0.0887 ± 0.0030 0.0884 ± 0.0029
-0.2929 ± 0.0088 -0.2925 ± 0.0095 -0.5220 ± 0.0412 -0.2997 ± 0.0110 -0.2960 ± 0.0084 -0.2954 ± 0.0104 -0.2933 ± 0.0087
0.9139 ± 0.1385 0.9543 ± 0.0693 0.0000 ± 0.0000 0.6171 ± 0.3546 0.8712 ± 0.2163 0.8984 ± 0.1424 0.9004 ± 0.1399
GBSG
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv QSurv
0.6554 ± 0.0294 0.6495 ± 0.0294 0.6336 ± 0.0280 0.6245 ± 0.0314 0.6529 ± 0.0286 0.6592 ± 0.0245 0.6558 ± 0.0297
0.1759 ± 0.0117 0.1822 ± 0.0176 0.2071 ± 0.0301 0.2001 ± 0.0782 0.1818 ± 0.0108 0.1769 ± 0.0120 0.1780 ± 0.0129
-0.5220 ± 0.0328 -0.5385 ± 0.0464 -0.5975 ± 0.0793 -0.7529 ± 0.9511 -0.5406 ± 0.0321 -0.5248 ± 0.0344 -0.5247 ± 0.0353
0.6697 ± 0.0430 0.6683 ± 0.0425 0.6497 ± 0.0401 0.6401 ± 0.0423 0.6666 ± 0.0484 0.6784 ± 0.0364 0.6765 ± 0.0453
0.0866 ± 0.0088 0.0877 ± 0.0084 0.0903 ± 0.0087 0.1172 ± 0.1302 0.0916 ± 0.0093 0.0873 ± 0.0074 0.0865 ± 0.0076
-0.2991 ± 0.0284 -0.3027 ± 0.0258 -0.3162 ± 0.0270 -0.4358 ± 0.5667 -0.3187 ± 0.0309 -0.3029 ± 0.0203 -0.2983 ± 0.0219
0.6678 ± 0.0291 0.6656 ± 0.0293 0.6531 ± 0.0314 0.6388 ± 0.0301 0.6649 ± 0.0286 0.6727 ± 0.0252 0.6697 ± 0.0299
0.1408 ± 0.0139 0.1433 ± 0.0142 0.1556 ± 0.0185 0.1724 ± 0.1255 0.1442 ± 0.0124 0.1403 ± 0.0113 0.1402 ± 0.0134
-0.4386 ± 0.0366 -0.4473 ± 0.0371 -0.4759 ± 0.0436 -0.6716 ± 0.9790 -0.4537 ± 0.0334 -0.4378 ± 0.0280 -0.4369 ± 0.0343
0.7995 ± 0.2822 0.7181 ± 0.3025 0.4138 ± 0.4239 0.6485 ± 0.3673 0.4328 ± 0.3027 0.3674 ± 0.2050 0.6001 ± 0.3362
METABRIC
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv QSurv
0.6346 ± 0.0200 0.6565 ± 0.0229 0.6626 ± 0.0182 0.6459 ± 0.0232 0.6365 ± 0.0178 0.6551 ± 0.0207 0.6542 ± 0.0198
0.1549 ± 0.0099 0.1578 ± 0.0102 0.1609 ± 0.0102 0.1626 ± 0.0101 0.1708 ± 0.0089 0.1554 ± 0.0110 0.1548 ± 0.0100
-0.4679 ± 0.0294 -0.4739 ± 0.0304 -0.4844 ± 0.0281 -0.4870 ± 0.0263 -0.5143 ± 0.0231 -0.4656 ± 0.0302 -0.4638 ± 0.0290
0.6229 ± 0.0322 0.6739 ± 0.0320 0.6669 ± 0.0300 0.6531 ± 0.0312 0.6309 ± 0.0334 0.6660 ± 0.0343 0.6683 ± 0.0313
0.1055 ± 0.0050 0.1032 ± 0.0063 0.1040 ± 0.0046 0.1044 ± 0.0055 0.1147 ± 0.0052 0.1033 ± 0.0039 0.1031 ± 0.0049
-0.3561 ± 0.0157 -0.3460 ± 0.0230 -0.3494 ± 0.0131 -0.3505 ± 0.0170 -0.3852 ± 0.0142 -0.3460 ± 0.0123 -0.3451 ± 0.0153
0.6368 ± 0.0260 0.6637 ± 0.0246 0.6608 ± 0.0232 0.6505 ± 0.0259 0.6413 ± 0.0223 0.6599 ± 0.0261 0.6602 ± 0.0220
0.1562 ± 0.0055 0.1533 ± 0.0070 0.1571 ± 0.0039 0.1565 ± 0.0058 0.1600 ± 0.0048 0.1533 ± 0.0052 0.1531 ± 0.0048
-0.4798 ± 0.0158 -0.4707 ± 0.0244 -0.4803 ± 0.0100 -0.4787 ± 0.0164 -0.4911 ± 0.0129 -0.4699 ± 0.0140 -0.4687 ± 0.0137
0.6223 ± 0.2984 0.5473 ± 0.3588 0.5653 ± 0.2963 0.3262 ± 0.3765 0.0000 ± 0.0000 0.3543 ± 0.2764 0.3580 ± 0.2688
NWTCO
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv QSurv
0.6422 ± 0.0345 0.6377 ± 0.0347 0.6360 ± 0.0266 0.6301 ± 0.0199 0.6532 ± 0.0273 0.6492 ± 0.0236 0.6480 ± 0.0209
0.1126 ± 0.0039 0.1139 ± 0.0047 0.1147 ± 0.0023 0.1155 ± 0.0027 0.1102 ± 0.0032 0.1120 ± 0.0030 0.1118 ± 0.0027
-0.3688 ± 0.0141 -0.3740 ± 0.0198 -0.3744 ± 0.0108 -0.3771 ± 0.0113 -0.3622 ± 0.0117 -0.3664 ± 0.0115 -0.3667 ± 0.0112
0.6343 ± 0.0347 0.6310 ± 0.0346 0.6282 ± 0.0263 0.6236 ± 0.0209 0.6448 ± 0.0269 0.6398 ± 0.0227 0.6398 ± 0.0204
0.0800 ± 0.0031 0.0810 ± 0.0038 0.0807 ± 0.0025 0.0817 ± 0.0030 0.0808 ± 0.0028 0.0798 ± 0.0030 0.0793 ± 0.0027
-0.2871 ± 0.0107 -0.2904 ± 0.0136 -0.2891 ± 0.0078 -0.2930 ± 0.0080 -0.2938 ± 0.0089 -0.2862 ± 0.0088 -0.2848 ± 0.0075
0.6394 ± 0.0335 0.6350 ± 0.0342 0.6342 ± 0.0259 0.6284 ± 0.0197 0.6498 ± 0.0263 0.6459 ± 0.0231 0.6449 ± 0.0202
0.1023 ± 0.0029 0.1036 ± 0.0038 0.1033 ± 0.0018 0.1046 ± 0.0023 0.1032 ± 0.0019 0.1020 ± 0.0028 0.1015 ± 0.0021
-0.3503 ± 0.0102 -0.3549 ± 0.0156 -0.3527 ± 0.0077 -0.3570 ± 0.0070 -0.3565 ± 0.0080 -0.3489 ± 0.0095 -0.3477 ± 0.0075
0.9777 ± 0.0471 0.9341 ± 0.2133 0.8889 ± 0.2197 0.7297 ± 0.3042 0.9426 ± 0.0962 0.9934 ± 0.0155 0.9401 ± 0.1839
SUPPORT2
CoxCC CoxTime DeepHit NnetSurv MDN DeSurv QSurv
0.6770 ± 0.0061 0.6941 ± 0.0052 0.6847 ± 0.0063 0.6722 ± 0.0058 0.6746 ± 0.0057 0.6950 ± 0.0048 0.6936 ± 0.0063
0.1798 ± 0.0046 0.1778 ± 0.0045 0.1914 ± 0.0041 0.1773 ± 0.0043 0.1886 ± 0.0041 0.1792 ± 0.0053 0.1778 ± 0.0042
-0.5449 ± 0.0201 -0.5336 ± 0.0140 -0.5653 ± 0.0102 -0.5293 ± 0.0117 -0.5564 ± 0.0106 -0.5357 ± 0.0157 -0.5308 ± 0.0110
0.7122 ± 0.0079 0.7508 ± 0.0092 0.7082 ± 0.0083 0.7115 ± 0.0060 0.7017 ± 0.0123 0.7501 ± 0.0094 0.7512 ± 0.0077
0.1202 ± 0.0022 0.1147 ± 0.0023 0.1238 ± 0.0021 0.1137 ± 0.0017 0.1370 ± 0.0023 0.1163 ± 0.0032 0.1154 ± 0.0028
-0.3896 ± 0.0065 -0.3680 ± 0.0071 -0.3954 ± 0.0055 -0.3650 ± 0.0049 -0.4407 ± 0.0060 -0.3723 ± 0.0092 -0.3697 ± 0.0080
0.6897 ± 0.0057 0.7061 ± 0.0052 0.6875 ± 0.0070 0.6827 ± 0.0059 0.6853 ± 0.0075 0.7067 ± 0.0055 0.7067 ± 0.0059
0.1927 ± 0.0031 0.1884 ± 0.0031 0.2090 ± 0.0020 0.1883 ± 0.0030 0.2138 ± 0.0022 0.1897 ± 0.0035 0.1885 ± 0.0032
-0.5714 ± 0.0106 -0.5573 ± 0.0081 -0.6042 ± 0.0045 -0.5565 ± 0.0081 -0.6146 ± 0.0055 -0.5597 ± 0.0083 -0.5562 ± 0.0075
0.0781 ± 0.0903 0.1508 ± 0.2010 0.0657 ± 0.0990 0.3449 ± 0.2590 0.0000 ± 0.0000 0.0160 ± 0.0420 0.0583 ± 0.1195
ODE integration under this evaluation protocol, SODEN was evaluated on the imaging datasets but omitted from the tabular experiments. Table 2 shows that QSurv is consistently competitive across datasets, horizons, and evaluation metrics. Careful hyperparameter optimization makes several baselines strong, including CoxCC, CoxTime, NnetSurv, DeepHit, MDN, and DeSurv. Nevertheless, QSurv frequently achieves the best or second-best performance, particularly for discrimination and survival probability estimation at restricted horizons. These results suggest that directly parameterizing the instantaneous hazard, while evaluating the cumulative hazard through quadrature, provides a robust balance between flexible time-to-event modeling and stable likelihood-based training. The advantage of QSurv is especially clear in the high-dimensional imaging settings, where flexible time-varying hazard modeling is most relevant. On COVID-19-NY, QSurv achieves the best fullhorizon Ctd and the best median-horizon Ctd , while remaining competitive with DeSurv and SODEN on likelihood- and calibration-oriented metrics. On C4KC-KiTS, QSurv obtains the best Ctd across all horizons and the strongest short-horizon IBS and IBLL. On BraTS, QSurv also achieves the best Ctd across all horizons and strong survival probability estimation, including the best median-horizon IBS and IBLL. These findings support the main motivation of QSurv: preserving continuous event times and directly learning instantaneous hazard can be beneficial when risk evolves nonlinearly over follow-up and covariates are high-dimensional. On tabular datasets, performance differences are generally smaller, and several classical or discretetime baselines remain highly competitive. QSurv nevertheless remains stable across datasets and horizons, often matching or improving on the best-performing alternatives. This suggests that the proposed quadrature-based hazard formulation does not sacrifice robustness in lower-dimensional 8
Figure 1: Predicted instantaneous hazard functions on COVID-19-NY across representative risk clusters. Lines show cluster-level predicted hazards, with event times marked by plus signs and censored observations by dots. NnetSurv produces irregular spike-like hazards, MDN produces smooth but constrained trajectories, SODEN produces near-flat covariate-stratified hazards, and DeSurv captures temporal variation but with pronounced oscillations and spikes. QSurv yields smoother and more interpretable cluster-specific temporal patterns, including elevated early risk in the high-risk cluster followed by gradual decline over time. settings, while providing additional flexibility for settings where proportional-hazards, parametric, or fixed-grid assumptions may be restrictive. To further examine model behavior beyond aggregate prediction metrics, we compared predicted instantaneous hazard functions across representative risk clusters on COVID-19-NY (Figure 1). Hazard clusters were formed by applying three-cluster K-means to subject-specific predicted hazard trajectories evaluated over a common time grid. This analysis highlights qualitative differences in how each method characterizes time-resolved risk. NnetSurv produces highly irregular, spike-like hazards. MDN yields smooth but relatively constrained, mostly monotone trajectories. SODEN produces nearly flat hazards that separate clusters mainly by risk level. DeSurv captures temporal variation, but its hazards remain irregular, with sharp oscillations and isolated spikes over follow-up. In contrast, QSurv produces smoother and more interpretable cluster-specific patterns: a high-risk cluster with elevated early hazard followed by gradual decline, an intermediate-risk cluster with persistent moderate risk, and a low-risk cluster with near-flat hazard. Unlike survival curves, which summarize accumulated risk up to each time point, these hazard trajectories directly identify when risk is concentrated over follow-up. The early high-risk pattern estimated by QSurv aligns with COVID-19 studies reporting concentrated early excess mortality after infection and elevated readmission-or-death risk shortly after hospital discharge, with attenuation over time [28, 13]. This comparison suggests that the qualitative difference between QSurv and SODEN may arise partly from how the cumulative hazard is evaluated during training, rather than from the hazard parameterization alone. Compared with SODEN, the fixed Gauss-Legendre quadrature rule evaluates the hazard at predetermined points over each subject’s observed interval, rather than relying on adaptive solver behavior. Compared with DeSurv, QSurv treats the hazard itself as the primitive modeled quantity rather than recovering it indirectly from a learned distribution function. The resulting hazard estimates are flexible enough to capture non-monotone and cluster-specific temporal risk patterns, while remaining substantially smoother and more stable than the discrete-time or indirectly derived alternatives.
5
Conclusion
We introduced QSurv, a quadrature-based nonparametric continuous-time survival model that directly parameterizes the instantaneous hazard as a flexible function of time and covariates. By evaluating the cumulative hazard with Gauss-Legendre quadrature, QSurv enables likelihood-based training without proportional-hazards assumptions, parametric survival distributions, time discretization, or adaptive ODE solvers. The resulting objective is simple, differentiable, and compatible with minibatch SGD, while Time-LoRA extends the framework to high-dimensional backbones by efficiently conditioning hazard predictions on time. Across simulations, tabular benchmarks, and imaging datasets, QSurv achieves competitive and robust performance. More importantly, its learned hazards provide time-resolved risk characterization, revealing clinically interpretable temporal patterns that may be obscured by aggregate survival summaries. A remaining limitation is that QSurv, like other deep survival models, depends on the hazard being learnable from available covariates and may require further validation under severe censoring, distribution shift, or highly irregular hazard dynamics. These results support QSurv as a practical framework for scalable and flexible continuous-time hazard modeling. 9
Acknowledgments and Disclosure of Funding The research of Chaeyeon Lee and Sehwan Kim was supported by the Global-Learning & Academic Research Institution for Master’s and PhD Students, and Postdocs (G-LAMP) Program of the National Research Foundation of Korea (NRF), funded by the Ministry of Education (No. RS-2025-25442252).
References [1] A. Avati, T. Duan, S. Zhou, K. Jung, N. H. Shah, and A. Y. Ng. Countdown regression: Sharp and calibrated survival predictions. In R. P. Adams and V. Gogate, editors, Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, volume 115 of Proceedings of Machine Learning Research, pages 145–155. PMLR, 22–25 Jul 2020. [2] S. Bakas, H. Akbari, A. Sotiras, M. Bilello, M. Rozycki, J. S. Kirby, J. B. Freymann, K. Farahani, and C. Davatzikos. Advancing the cancer genome atlas glioma mri collections with expert segmentation labels and radiomic features. Scientific data, 4(1):170117, 2017. [3] S. Bakas, M. Reyes, A. Jakab, S. Bauer, M. Rempfler, A. Crimi, R. T. Shinohara, C. Berger, S. M. Ha, M. Rozycki, et al. Identifying the best machine learning algorithms for brain tumor segmentation, progression assessment, and overall survival prediction in the brats challenge. arXiv preprint arXiv:1811.02629, 2018. [4] A. Bennis, S. Mouysset, and M. Serrurier. Estimation of conditional mixture weibull distribution with right censored data using neural network for time-to-event analysis. In Advances in Knowledge Discovery and Data Mining: 24th Pacific-Asia Conference, PAKDD 2020, Singapore, May 11–14, 2020, Proceedings, Part I, Berlin, Heidelberg, 2020. Springer-Verlag. ISBN 9783-030-47425-6. doi: 10.1007/978-3-030-47426-3_53. URL https://doi.org/10.1007/ 978-3-030-47426-3_53. [5] N. E. Breslow and N. Chatterjee. Design and analysis of two-phase studies with binary outcome applied to wilms tumour prognosis. Journal of the Royal Statistical Society: Series C (Applied Statistics), 48(4):457–468, 1999. [6] A. F. Connors, N. V. Dawson, N. A. Desbiens, W. J. Fulkerson, L. Goldman, W. A. Knaus, J. Lynn, R. K. Oye, M. Bergner, A. Damiano, et al. A controlled trial to improve care for seriously iii hospitalized patients: The study to understand prognoses and preferences for outcomes and risks of treatments (support). Jama, 274(20):1591–1598, 1995. [7] D. R. Cox. Regression models and life-tables. Journal of the Royal Statistical Society. Series B (Methodological), 34(2):187–220, 1972. ISSN 00359246. URL http://www.jstor.org/ stable/2985181. [8] E. Craig, C. Zhong, and R. Tibshirani. Survival stacking: casting survival analysis as a classification problem, 2021. URL https://arxiv.org/abs/2107.13480. [9] C. Curtis, S. P. Shah, S.-F. Chin, G. Turashvili, O. M. Rueda, M. J. Dunning, D. Speed, A. G. Lynch, S. Samarajiwa, Y. Yuan, et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature, 486(7403):346–352, 2012. [10] D. Danks and C. Yau. Derivative-based neural modelling of cumulative distribution functions for survival analysis. In G. Camps-Valls, F. J. R. Ruiz, and I. Valera, editors, Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pages 7240–7256. PMLR, 28–30 Mar 2022. URL https://proceedings.mlr.press/v151/danks22a.html. [11] P. J. Davis and P. Rabinowitz. Methods of numerical integration. Courier Corporation, 2007. [12] A. Dispenzieri, J. A. Katzmann, R. A. Kyle, D. R. Larson, T. M. Therneau, C. L. Colby, R. J. Clark, G. P. Mead, S. Kumar, L. J. Melton III, et al. Use of nonclonal serum immunoglobulin free light chains to predict overall survival in the general population. In Mayo Clinic Proceedings, volume 87, pages 517–523. Elsevier, 2012. 10
[13] J. P. Donnelly, X. Q. Wang, T. J. Iwashyna, and H. C. Prescott. Readmission and death after initial hospital discharge among patients with covid-19 in a large multihospital system. Jama, 325(3):304–306, 2021. [14] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, et al. An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010.11929, 2020. [15] V. Dumoulin, E. Perez, N. Schucher, F. Strub, H. d. Vries, A. Courville, and Y. Bengio. Featurewise transformations. Distill, 3(7):e11, 2018. [16] S. Fotso. Deep neural networks for survival analysis based on a multi-task framework. arXiv preprint arXiv:1801.05512, 2018. [17] M. F. Gensheimer and B. Narasimhan. A scalable discrete-time survival model for neural networks. PeerJ, 7:e6257, 2019. [18] G. H. Golub and J. H. Welsch. Calculation of gauss quadrature rules. Mathematics of computation, 23(106):221–230, 1969. [19] E. Graf, C. Schmoor, W. Sauerbrei, and M. Schumacher. Assessment and comparison of prognostic classification schemes for survival data. Statistics in medicine, 18(17-18):2529–2545, 1999. [20] H. Haider, B. Hoehn, S. Davis, and R. Greiner. Effective ways to build and evaluate individual survival distributions. Journal of Machine Learning Research, 21(85):1–63, 2020. [21] X. Han, M. Goldstein, and R. Ranganath. Survival mixture density networks. In Machine Learning for Healthcare Conference, pages 224–248. PMLR, 2022. [22] K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770– 778, 2016. [23] M. A. Hernán. The hazards of hazard ratios. Epidemiology, 21(1):13–15, 2010. [24] K. R. Hess and V. A. Levin. Getting more out of survival data by using the hazard function. Clinical Cancer Research, 20(6):1404–1409, 2014. [25] E. J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, W. Chen, et al. Lora: Low-rank adaptation of large language models. ICLR, 1(2):3, 2022. [26] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger. Densely connected convolutional networks. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4700–4708, 2017. [27] H. Ishwaran, U. B. Kogalur, E. H. Blackstone, and M. S. Lauer. Random survival forests. The Annals of Applied Statistics, 2(3):841 – 860, 2008. doi: 10.1214/08-AOAS169. URL https://doi.org/10.1214/08-AOAS169. [28] T. J. Iwashyna, S. Seelye, T. S. Berkowitz, J. Pura, A. S. Bohnert, C. B. Bowling, E. J. Boyko, D. M. Hynes, G. N. Ioannou, M. L. Maciejewski, et al. Late mortality after covid-19 infection among us veterans vs risk-matched comparators: a 2-year cohort analysis. JAMA internal medicine, 183(10):1111–1119, 2023. [29] E. L. Kaplan and P. Meier. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282):457–481, 1958. ISSN 01621459. URL http://www.jstor.org/stable/2281868. [30] J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, and Y. Kluger. Deepsurv: personalized treatment recommender system using a cox proportional hazards deep neural network. BMC medical research methodology, 18(1):24, 2018. 11
[31] S. M. Kazemi, R. Goel, S. Eghbali, J. Ramanan, J. Sahota, S. Thakur, S. Wu, C. Smyth, P. Poupart, and M. Brubaker. Time2vec: Learning a vector representation of time. arXiv preprint arXiv:1907.05321, 2019. [32] H. Kvamme, Ø. Borgan, and I. Scheel. Time-to-event prediction with neural networks and cox regression. Journal of machine learning research, 20(129):1–30, 2019. [33] C. Lee, W. Zame, J. Yoon, and M. Van Der Schaar. Deephit: A deep learning approach to survival analysis with competing risks. In Proceedings of the AAAI conference on artificial intelligence, volume 32, 2018. [34] B. H. Menze, A. Jakab, S. Bauer, J. Kalpathy-Cramer, K. Farahani, J. Kirby, Y. Burren, N. Porz, J. Slotboom, R. Wiest, et al. The multimodal brain tumor image segmentation benchmark (brats). IEEE Transactions on Medical Imaging, 34(10):1993–2024, 2015. doi: 10.1109/TMI.2014.2377694. [35] C. Nagpal, X. Li, and A. Dubrawski. Deep survival machines: Fully parametric survival regression and representation learning for censored data with competing risks. IEEE Journal of Biomedical and Health Informatics, 25(8):3163–3175, 2021. [36] C. Nagpal, S. Yadlowsky, N. Rostamzadeh, and K. Heller. Deep cox mixtures for survival regression. In Machine Learning for Healthcare Conference, pages 674–708. PMLR, 2021. [37] National Heart, Lung, and Blood Institute. About BioLINCC. Online, 2022. URL https: //biolincc.nhlbi.nih.gov/about/. [38] R. Ranganath, A. Perotte, N. Elhadad, and D. Blei. Deep survival analysis. In F. Doshi-Velez, J. Fackler, D. Kale, B. Wallace, and J. Wiens, editors, Proceedings of the 1st Machine Learning for Healthcare Conference, volume 56 of Proceedings of Machine Learning Research, pages 101–114, Northeastern University, Boston, MA, USA, 18–19 Aug 2016. PMLR. [39] P. Royston. Flexible parametric alternatives to the cox model, and more. The Stata Journal, 1 (1):1–28, 2001. [40] M. Schumacher, G. Bastert, H. Bojar, K. Hübner, M. Olschewski, W. Sauerbrei, C. Schmoor, C. Beyerle, R. Neumann, and H. Rauschecker. Randomized 2 x 2 trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients. german breast cancer study group. Journal of Clinical Oncology, 12(10):2086–2093, 1994. [41] M. J. Stensrud and M. A. Hernán. Why test for proportional hazards? Jama, 323(14):1401–1402, 2020. [42] W. Tang, J. Ma, Q. Mei, and J. Zhu. Soden: A scalable continuous-time survival model through ordinary differential equation networks. Journal of Machine Learning Research, 23(34):1–29, 2022. [43] W. Tang, K. He, G. Xu, and J. Zhu. Survival analysis via ordinary differential equations. Journal of the American Statistical Association, 118(544):2406–2421, 2023. [44] H. Uno, T. Cai, M. J. Pencina, R. B. D’Agostino, and L.-J. Wei. On the c-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data. Statistics in medicine, 30(10):1105–1117, 2011. [45] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin. Attention is all you need. In Advances in neural information processing systems, volume 30, 2017. [46] S. Wiegrebe, P. Kopper, R. Sonabend, B. Bischl, and A. Bender. Deep learning for survival analysis: a review. Artificial Intelligence Review, 57(3), Feb. 2024. ISSN 1573-7462. doi: 10. 1007/s10462-023-10681-3. URL http://dx.doi.org/10.1007/s10462-023-10681-3. [47] C.-N. Yu, R. Greiner, H.-C. Lin, and V. Baracos. Learning patient-specific cancer survival distributions as a sequence of dependent regressors. Advances in neural information processing systems, 24, 2011. 12
[48] Q. Zhong, J. W. Mueller, and J.-L. Wang. Deep extended hazard models for survival analysis. Advances in Neural Information Processing Systems, 34:15111–15124, 2021.
13
A
Gauss-Legendre Quadrature
Gauss-Legendre quadrature approximates a definite integral by a weighted sum of function evaluations at carefully chosen nonuniform nodes. Unlike grid-based rules such as Riemann sums or the trapezoidal rule, which evaluate the integrand at equally spaced points, Gauss-Legendre quadrature chooses nodes and weights to maximize algebraic precision. For an integral over the standard interval [−1, 1], the K-point rule is Z 1 K X q(u)du ≈ wk q(ξk ), (9) −1
k=1
where ξk are the roots of the K-th Legendre polynomial PK and 2 (10) wk = 2 ′ (ξ )]2 . (1 − ξk )[PK k The nodes and weights are fixed by the quadrature rule, not learned from data. With K function evaluations, the rule is exact for all polynomials of degree at most 2K − 1. For an integral over a general interval [a, b], the affine transformation b−a b+a s= u+ 2 2 maps u ∈ [−1, 1] to s ∈ [a, b]. Therefore, Z b K b−a X b−a b+a q(s)ds ≈ wk q ξk + . 2 2 2 a
(11)
(12)
k=1
When q is 2K times continuously differentiable, the classical Gauss-Legendre error formula gives (b − a)2K+1 (K!)4 (2K) q (ζ) (13) (2K + 1)[(2K)!]3 for some ζ ∈ (a, b). Thus, the approximation error is controlled by the high-order smoothness of the integrand. For analytic functions, Gauss-Legendre quadrature can converge much faster than low-order grid-based rules, often requiring substantially fewer evaluations to reach a comparable accuracy. EK (q) =
In QSurv, the integrand is the neural network parameterized hazard q(s) = λθ (s|x) and the interval is subject-specific, [0, oi ]. Applying the above transformation with a = 0 and b = oi gives Z oi K oi X Λθ (oi |xi ) = λθ (s|xi )ds ≈ wk λθ (oi τk |xi ), (14) 2 0 k=1
where τk = (ξk + 1)/2. This is the approximation used in the QSurv likelihood: the hazard remains a continuous-time function, while the cumulative hazard is evaluated through a fixed number of differentiable neural network parameterized hazard evaluations.
B
Proofs
Rt Proof of Theorem 3.1. The cumulative hazard is defined as the integral Λθ (t|x) = 0 λθ (s|x) ds. To apply the standard error bound for Gauss-Legendre quadrature, which is defined on the reference interval [−1, 1], we perform an affine transformation of variables. Let s(u) = 2t (u + 1). This maps u ∈ [−1, 1] to s ∈ [0, t]. The differential is ds = 2t du. The standard error term for K-point 2K+1
4
2 (K!) (2K) Gauss-Legendre quadrature on [−1, 1] for a function g is given by EK (g) = (2K+1)[(2K)!] (ξ) 3 g t for some ξ ∈ (−1, 1). Let g(u) = λθ ( 2 (u + 1)|x). By the chain rule, the 2K-th derivative scales as (2K) g (2K) (u) = ( 2t )2K λθ (s(u)|x). Evaluating at the point ξ ∈ (−1, 1) given by the quadrature error (2K) theorem and letting τ = s(ξ) = (t/2)(ξ + 1) ∈ (0, t), we obtain g (2K) (ξ) = (t/2)2K λθ (τ |x). t The total error is the standard error scaled by the integration measure 2 : " # 2K t 22K+1 (K!)4 t (2K) Error = · · λθ (τ |x) . (15) 2 (2K + 1)[(2K)!]3 2
14
Figure 2: Convergence of approximation error and training efficiency on two simulation scenarios. The first three plots display the integrated absolute error (IAE) for survival probability, cumulative hazard, and instantaneous hazard as a function of quadrature nodes (K). All metrics exhibit rapid convergence, saturating beyond K = 5 across both simulation scenarios. The fourth panel confirms that while training time increases linearly with K. Shaded regions denote standard deviation (n = 20). See Appendix D and C for simulation details and estimated curves, respectively. 2K+1
4
(2K)
t (K!) Grouping terms, the factors of 2 cancel, leaving the coefficient (2K+1)[(2K)!] 3 . Bounding |λθ
(τ |x)|
(2K) by maxτ ∈[0,t] |λθ (τ |x)| yields the inequality.
Proof of Corollary 3.2. The negative log-likelihood decomposes into a log-hazard term and a cumulative hazard term. Since the hazard function λθ (t|x) is computed exactly via the neural network forward pass at time t, the term −δ log λθ (t|x) is identical in both the true and approximate losses. Consequently, the difference depends solely on the integral approximation: |L − L̂K | = |(−δ log λ + Λ) − (−δ log λ + Λ̂K )| = |Λ − Λ̂K | = εK (t|x).
C
Number of Quadrature Nodes
We further investigate how the number of quadrature nodes K affects the estimation of the instantaneous hazard, cumulative hazard, and survival function. Although QSurv represents λθ (t|x) as a continuous-time neural function, the training objective evaluates this function only at finitely many time points: the observed event times through the log-hazard contribution and the subject-specific quadrature nodes through the cumulative-hazard contribution. Therefore, if K is too small relative to the temporal complexity of the hazard, the quadrature-based objective may fail to resolve rapid changes in risk over time. This observation clarifies the role of K in QSurv. Unlike the number of bins in discrete-time survival models, K does not define the temporal resolution of the hazard model itself. The hazard remains a 15
Figure 3: Instantaneous hazard, cumulative hazard, and survival functions for simulation scenario 1.
Figure 4: Instantaneous hazard, cumulative hazard, and survival functions for simulation scenario 2. continuous function of time. Instead, K controls the numerical precision with which the cumulative hazard is approximated inside the likelihood. Reducing K should therefore not be interpreted as imposing a simpler or more regularized hazard model. Rather, very small values of K may lead to a poor approximation of the continuous-time likelihood, especially when the hazard has sharp or non-monotone temporal structure. Conversely, once K is large enough to resolve the relevant temporal variation, increasing K further should have little effect on estimation accuracy, aside from increased computation. To empirically assess this behavior, we conducted a simulation study using two synthetic scenarios with complex, non-monotone hazards. Scenario 1 is a crossing-hazards setting with a binary covariate x ∼ Bernoulli(0.5), defined by S(t|x) = exp(−t1+x ),
λ(t|x) = (1 + x)tx .
(16)
Thus, the x = 0 group follows S(t|0) = exp(−t) with constant hazard λ(t|0) = 1, whereas the x = 1 group follows S(t|1) = exp(−t2 ) with linearly increasing hazard λ(t|1) = 2t. The hazards intersect at t = 0.5, while the survival curves cross at t = 1, where S(t|0) = S(t|1) = e−1 . Censoring times are sampled from C ∼ Uniform(0, 2). Scenario 2 introduces higher-frequency temporal variation through anti-phase sinusoidal hazards: λ(t|x) = 1 + 0.8(1 − 2x) sin(4t), Λ(t|x) = t + 0.2(1 − 2x){1 − cos(4t)}, (17) S(t|x) = exp(−Λ(t|x)). For x = 0, the hazard is 1 + 0.8 sin(4t), while for x = 1, the hazard is 1 − 0.8 sin(4t). This scenario creates oscillating risk dynamics with multiple hazard intersections, requiring the model to resolve rapid temporal changes in risk. Survival times are generated by inverse transform sampling, with censoring drawn from an exponential distribution with rate approximately 0.33. 16
We trained QSurv models with K ∈ {1, 2, 3, 5, 7, 10} quadrature nodes across 20 independent runs. The results are summarized in Figure 2. We observe a rapid decrease in integrated absolute error for the instantaneous hazard, cumulative hazard, and survival function as K increases. The error curves stabilize after a moderate number of quadrature nodes, indicating that the quadrature approximation has become sufficiently accurate for these settings. As expected, training time increases approximately linearly with K, reflecting the additional neural network parameterized hazard evaluations required at the quadrature nodes. Figures 3 and 4 show the corresponding estimated trajectories. With very small values of K, the fitted hazard can exhibit visible approximation artifacts and may fail to capture important temporal features of the true hazard. In contrast, moderate values of K recover stable estimates of the instantaneous hazard, cumulative hazard, and survival function. These findings support our practical use of a fixed moderate quadrature order in the main experiments.
D
Simulation Study
To evaluate the capability of the proposed models to capture complex, non-linear relationships between covariates and survival outcomes, we generated synthetic datasets using six parametric distributions: Exponential, Weibull, Gamma, Gompertz, Log-Normal, and Log-Logistic. For all scenarios, we simulated a one-dimensional covariate x drawn from a uniform distribution between -1 and 1. The distribution parameters, such as shape, scale, or rate, were modeled as non-linear functions of x. We employed third-order polynomial transformations followed by an exponential link function to ensure positivity, except for the location parameter in the Log-Normal distribution. The polynomial coefficients were manually tuned to ensure diverse hazard shapes and stable event times within the observation window. The chosen distributions cover a wide range of survival behaviors. The Exponential and Gompertz scenarios represent Proportional Hazards settings, whereas the Weibull, Gamma, Log-Normal, and Log-Logistic distributions induce non-proportional hazards where the shape of the distribution changes conditional on the covariate. Notably, the Log-Normal and Log-Logistic distributions generate unimodal hazard functions that rise and then fall, mimicking patterns often seen in medical contexts. Censoring times were generated independently of the covariates using a uniform distribution, where the upper bound was dynamically calibrated for each dataset to achieve a target censoring rate of approximately 20%. The observed time was defined as the minimum of the true event time and the censoring time. For each simulation run, we generated 2,000 training and testing samples, Table 3: Data generation mechanisms. All covariate-dependent parameters are modeled as θ(x) = exp(P(x; w)) to ensure positivity, except for µ(x) in the Log-Normal distribution, which is modeled directly as µ(x) = P(x; w). The polynomial function is defined as P(x; w) = w0 + w1 x + w2 x2 + w3 x3 . Distribution
Hazard / PDF Definition
Parameter Coefficients
Exponential
h(t|x) = λ(x)
wλ = [−1.0, 0.5, −0.3, 0.15]
k(x)−1
Weibull
k(x) h(t|x) = λ(x)
Gamma
tk(x)−1 e−β(x)t f (t|x) = β(x) Γ(k(x))
wk = [1.8, 0.3, −0.1, 0.05] wβ = [0.3, −0.4, 0.15, −0.05]
Gompertz
h(t|x) = b(x)ect
wb = [−2.0, 0.4, −0.2, 0.1] c = 0.05
Log-Normal
T = exp(µ(x) + σ(x)Z), Z ∼ N (0, 1)
wµ = [1.5, 0.8, −0.4, 0.2] wσ = [−0.1, 0.25, −0.10, 0.03]
Log-Logistic
h(t|x) = (β(x)/α(x))(t/α(x)) 1+(t/α(x))β(x)
t λ(x)
k(x)
β(x)−1
17
wk = [0.3, 0.2, −0.1, 0.05] wλ = [2.0, 0.3, −0.2, 0.1]
wα = [1.2, 0.4, −0.15, 0.08] wβ = [1.0, 0.3, −0.1, 0.05]
Table 4: Performance comparison of survival models across various distributions. Values are Mean ± Std (20 seeds). Best results are bold, second best are underlined. Training times are wall-clock seconds on a single Apple M4 CPU (10-core, 32 GB RAM); relative comparisons across rows are meaningful but absolute values are machine-dependent. Metric
Distribution
CoxCC
CoxTime
NnetSurv
MDN
DeSurv
SODEN
QSurv
L1 Error of Ŝ
Exponential Weibull Gamma Gompertz Log-Normal Log-Logistic
0.0122 ± 0.0038 0.0268 ± 0.0041 0.0197 ± 0.0041 0.0121 ± 0.0040 0.0166 ± 0.0065 0.0356 ± 0.0040
0.0130 ± 0.0042 0.0144 ± 0.0038 0.0138 ± 0.0035 0.0143 ± 0.0054 0.0166 ± 0.0049 0.0157 ± 0.0042
0.0266 ± 0.0084 0.0236 ± 0.0041 0.0213 ± 0.0032 0.0276 ± 0.0058 0.0199 ± 0.0054 0.0209 ± 0.0035
0.0155 ± 0.0055 0.0354 ± 0.0037 0.1056 ± 0.0028 0.0202 ± 0.0043 0.0272 ± 0.0034 0.0567 ± 0.0041
0.0177 ± 0.0057 0.0127 ± 0.0043 0.0121 ± 0.0035 0.0129 ± 0.0043 0.0152 ± 0.0054 0.0144 ± 0.0041
0.0167 ± 0.0069 0.0139 ± 0.0053 0.0107 ± 0.0029 0.0115 ± 0.0034 0.0165 ± 0.0054 0.0150 ± 0.0042
0.0162 ± 0.0058 0.0133 ± 0.0046 0.0110 ± 0.0030 0.0120 ± 0.0032 0.0160 ± 0.0064 0.0154 ± 0.0042
L1 Error of Λ̂
Exponential Weibull Gamma Gompertz Log-Normal Log-Logistic
0.1780 ± 0.0479 0.7486 ± 0.1803 2.4663 ± 0.3089 0.2773 ± 0.1030 0.8253 ± 0.2925 0.6392 ± 0.1234
0.2013 ± 0.0638 0.2372 ± 0.1229 1.1946 ± 0.6717 0.3699 ± 0.2559 0.6984 ± 0.3496 0.2488 ± 0.0619
0.4163 ± 0.1258 0.5065 ± 0.1114 1.3610 ± 0.4712 0.7631 ± 0.1900 0.4052 ± 0.1085 0.3262 ± 0.0438
0.2244 ± 0.0958 0.5137 ± 0.1197 2.8405 ± 0.4260 0.6015 ± 0.1340 1.0280 ± 0.4111 0.4433 ± 0.0416
0.2187 ± 0.0900 0.1974 ± 0.0644 0.6507 ± 0.2066 0.2870 ± 0.1289 0.3427 ± 0.1220 0.2279 ± 0.1152
0.2144 ± 0.0827 0.1887 ± 0.0842 0.5796 ± 0.2705 0.1964 ± 0.1009 0.4213 ± 0.1943 0.2358 ± 0.1280
0.2221 ± 0.0923 0.1894 ± 0.0938 0.5766 ± 0.2472 0.1990 ± 0.1021 0.4205 ± 0.2091 0.2305 ± 0.1158
L1 Error of λ̂
Exponential Weibull Gamma Gompertz Log-Normal Log-Logistic
0.1775 ± 0.0410 0.1576 ± 0.0322 0.9012 ± 0.0921 0.1713 ± 0.0262 0.1235 ± 0.0248 0.1621 ± 0.0210
0.1775 ± 0.0422 0.1182 ± 0.0207 0.5661 ± 0.1769 0.1890 ± 0.0549 0.1114 ± 0.0294 0.1266 ± 0.0242
0.3830 ± 0.0438 0.2383 ± 0.0254 0.7249 ± 0.0370 0.2921 ± 0.0276 0.1938 ± 0.0219 0.3290 ± 0.0525
0.0568 ± 0.0446 0.0644 ± 0.0102 0.5646 ± 0.0459 0.0920 ± 0.0140 0.1882 ± 0.0520 0.1031 ± 0.0090
0.0448 ± 0.0159 0.0295 ± 0.0090 0.2930 ± 0.0907 0.0502 ± 0.0278 0.0473 ± 0.0201 0.0513 ± 0.0299
0.0330 ± 0.0138 0.0272 ± 0.0119 0.1716 ± 0.0901 0.0262 ± 0.0174 0.0462 ± 0.0183 0.0517 ± 0.0295
0.0352 ± 0.0164 0.0282 ± 0.0139 0.1718 ± 0.0859 0.0258 ± 0.0180 0.0457 ± 0.0185 0.0491 ± 0.0291
Training Time (s)
Exponential Weibull Gamma Gompertz Log-Normal Log-Logistic
5.3210 ± 0.7844 5.0815 ± 0.4088 4.6292 ± 0.3867 4.4591 ± 0.2230 5.8264 ± 2.2276 4.3939 ± 0.8087
6.7244 ± 0.8221 6.6068 ± 0.5126 6.3100 ± 0.6772 5.9313 ± 0.3321 7.8568 ± 3.4891 5.8667 ± 1.2868
2.4345 ± 0.2233 2.4946 ± 0.4445 2.2909 ± 0.2931 1.9610 ± 0.1278 2.8560 ± 1.1757 2.3401 ± 1.8595
3.5458 ± 0.2663 3.6100 ± 0.5023 3.4942 ± 0.4543 3.0111 ± 0.0375 4.3230 ± 1.8299 3.3814 ± 0.9135
6.1273 ± 1.2448 5.9509 ± 0.6691 6.0968 ± 0.9452 5.5127 ± 0.4924 8.7879 ± 6.4600 5.0744 ± 0.6543
87.9615 ± 8.5607 124.1380 ± 17.2815 149.9774 ± 24.4947 82.9457 ± 10.5067 243.5298 ± 117.4764 131.3422 ± 13.7360
6.0749 ± 0.8276 5.7566 ± 0.4297 6.2316 ± 1.1584 5.3545 ± 0.3065 8.0038 ± 4.6329 5.7363 ± 0.6092
respectively. The experiments were repeated 20 times with different random seeds to report the mean and standard deviation of the results. All neural network baselines used a standardized architecture with two hidden layers of 32 units and Tanh activations. To strictly assess how well the models recovered the underlying distributions, we calculated the L1 error for the hazard, cumulative hazard, and survival functions on the testing samples by comparing the model estimates against the true analytical functions. To assess the model’s ability to recover the global system dynamics independent of individual covariate variations, we computed the marginalized survival, cumulative hazard, and instantaneous hazard curves. These are obtained by averaging the conditional functions over the distribution of the covariate space X : Z n 1X S̄(t) = S(t|x)p(x)dx ≈ Ŝ(t|xi ), n i=1 X Z n 1X Λ̄(t) = Λ(t|x)p(x)dx ≈ Λ̂(t|xi ), (18) n i=1 X Z n 1X λ̄(t) = λ(t|x)p(x)dx ≈ λ̂(t|xi ). n i=1 X By comparing these marginalized trajectories against the true analytical integrations, we can visualize whether the model correctly captures the aggregate behavior of the population, including the mean hazard shape and the overall survival decay, even when the underlying individual parameters θ vary non-linearly. Figures 6–11 present the results.
E
LoRA Ablation Study
We performed an ablation study to compare different time-conditioning mechanisms within QSurv. Specifically, we compared three variants: direct concatenation of time and covariates, FiLM-based conditioning, and the proposed Time-LoRA conditioning. All variants used the same simulation settings, training protocol, and evaluation metric. Performance was measured by the L1 error between the predicted and true instantaneous hazard functions, averaged over 20 random seeds. Table 5 shows that no single conditioning mechanism uniformly dominates across all data-generating distributions. Direct concatenation performs best for Weibull, Gamma, and log-normal hazards, suggesting that simple conditioning can be sufficient when the hazard structure is relatively easy to represent in the low-dimensional simulation setting. FiLM performs best for the exponential setting and remains competitive for Gompertz and log-logistic hazards. Time-LoRA achieves the 18
Table 5: Comparison of QSurv conditioning variants on L1 error of instantaneous hazard (smaller the better) across different parametric distributions, with standard deviations over 20 random seeds. Best results are bold, second best are underlined. Distribution QSurv (Concat)
QSurv (FiLM)
QSurv (LoRA)
0.0327 ± 0.0159 0.0278 ± 0.0055 0.1199 ± 0.0344 0.0506 ± 0.0157 0.0431 ± 0.0154 0.0939 ± 0.0283
0.0260 ± 0.0104 0.0345 ± 0.0182 0.2033 ± 0.0894 0.0421 ± 0.0403 0.0469 ± 0.0197 0.0590 ± 0.0507
0.0352 ± 0.0164 0.0282 ± 0.0139 0.1718 ± 0.0859 0.0258 ± 0.0180 0.0457 ± 0.0185 0.0491 ± 0.0291
Exponential Weibull Gamma Gompertz Log-normal Log-logistic
best performance for Gompertz and log-logistic hazards, and is second best for Weibull, Gamma, and log-normal distributions. These results suggest that Time-LoRA is a competitive time-conditioning mechanism for hazard modeling, particularly in settings where temporal risk patterns require more flexible interactions between time and covariates. At the same time, the ablation indicates that in simple low-dimensional simulations, direct concatenation can remain a strong baseline. The main motivation for Time-LoRA is therefore not that it uniformly improves all synthetic settings, but that it provides an efficient and flexible mechanism for conditioning high-dimensional neural backbones on time while preserving competitive hazard-estimation accuracy. Figures 6–11 present the visualizations of the fitted curves.
F
Experiment Details
F.1
Datasets
We evaluated all methods on a mixture of tabular and imaging survival datasets. The main benchmark includes six tabular datasets (GBSG, METABRIC, NWTCO, Framingham, FLCHAIN, and SUPPORT2) and three imaging datasets (COVID-19-NY, C4KC-KiTS, and BraTS). GBSG comes from a multicenter randomized clinical trial conducted by the German Breast Cancer Study Group to evaluate adjuvant treatment strategies in patients with lymph node-positive breast cancer [40]. The endpoint is breast cancer recurrence after treatment. The public dataset contains 686 patients, with 299 observed recurrence events (43.6%) and 387 right-censored observations (56.4%). We used eight covariates covering demographic features, tumor characteristics, hormone receptor status, and treatment information. Categorical variables were one-hot encoded. METABRIC comes from the Molecular Taxonomy of Breast Cancer International Consortium and was designed to study molecular subtypes of breast cancer [9]. We used the version distributed with Pycox, following the Immunohistochemical 4 plus Clinical (IHC4+C) covariate set [30]. The input variables include four gene-expression markers (MKI67, EGFR, PGR, and ERBB2) and five clinical variables (hormone therapy, radiotherapy, chemotherapy, estrogen receptor status, and age). The dataset contains 1,904 patients, with 1,103 deaths (57.9%) and 801 censored observations (42.1%). Missing categorical values were imputed with the most frequent category before one-hot encoding. NWTCO is from the third and fourth National Wilms Tumor Study Group trials and includes pediatric patients diagnosed with Wilms tumor between 1979 and 1994 [5]. We analyzed 4,028 patients with complete relapse outcome and covariate information. The covariates include disease stage, age, and subcohort membership. There were 571 observed events (14.1%) and 3,457 censored observations (85.8%). Framingham uses the Framingham Heart Study longitudinal teaching dataset distributed by NHLBI/BioLINCC [37]. This is a Framingham-derived teaching dataset rather than the full Framingham cohort. It contains longitudinal clinic, questionnaire, laboratory, and adjudicated event data on 4,434 participants from three examination periods, yielding 11,627 person-exam observations and 39 variables. Outcomes include angina, myocardial infarction, stroke or cardiovascular disease, hypertension, and death, with corresponding time-to-event or censoring variables. In our benchmark, we used cardiovascular disease (CVD) as the survival endpoint, with CVD as the event indicator and 19
TIMECVD as the observed follow-up time. CVD events were observed in 996 observations (23.3%), and 3,277 observations (76.7%) were right-censored. FLCHAIN comes from a study of serum free light chains and mortality among residents aged 50 years or older in Olmsted County, Minnesota [12]. The final sample contains 7,874 participants selected by stratified random sampling over age and sex. We used six covariates: age, creatinine, kappa and lambda free light chain levels, MGUS status, and sex. There were 2,167 deaths (27.5%) and 5,707 censored observations (72.5%). Missing continuous values were imputed with the mean. SUPPORT2 comes from the SUPPORT Phase II study, a multicenter study of seriously ill hospitalized adults conducted between 1992 and 1994 [6]. We used 9,105 observations with 47 covariates, including demographics, disease severity, physiologic measurements, and clinical status variables. Death was observed in 6,201 observations (68.1%), and 2,904 observations (31.9%) were rightcensored. Observations with follow-up of three days or less were excluded so that baseline covariates preceded the event time. COVID-19-NY consists of patients who tested positive for COVID-19 at Stony Brook University. We used the linked clinical file and AP chest radiographs from the TCIA release. For each patient, the radiograph study closest to the hospital visit date was selected among portable/AP chest views, and the first available DICOM image from that study was converted to a grayscale PNG. The endpoint was in-hospital death, coded from discharge status, and observed time was measured in days from the selected imaging study to discharge or death. Images were padded to a square field of view and resized to 224 × 224 before model input. C4KC-KiTS is a TCIA renal cancer imaging dataset with arterial-phase abdominal CT and kidney/tumor segmentations. Original imaging consisted of CT DICOM series and corresponding segmentation objects. For the 2D slice selection, CT volumes were resampled to 1 mm isotropic spacing, clipped to an abdominal soft-tissue window ([−160, 240] HU), and reduced to the axial slice with the largest kidney tumor cross-section according to the segmentation mask. The selected slice was resized to 224 × 224 and used as a single-channel grayscale image. The endpoint was overall survival after surgery: observed time was vital_days_after_surgery, and death was coded from vital_status; alive or censored patients were treated as right-censored. Only malignant cases with valid follow-up were included. BraTS uses the BraTS 2020 glioma MRI training data [34, 2, 3]. Original imaging consisted of multi-parametric brain MRI volumes with segmentation masks; in this benchmark we used the FLAIR sequence. For the 2D slice selection, we extracted the axial FLAIR slice with the largest tumor cross-section using the segmentation mask, falling back to the center slice only if the mask was unavailable or empty. Intensities were clipped at the within-slice 99th percentile, normalized, and saved as 224 × 224 grayscale images. The endpoint was overall survival time in days from the BraTS survival file. We constructed the event indicator from the available vital-status/alive annotation in the processed survival file; subjects annotated as alive with follow-up time were treated as right-censored, while subjects without an alive annotation were treated as observed events. F.2
Evaluation Horizons
Metrics were reported at three horizons: the full supported time horizon and two shorter horizons. The full horizon was chosen from the training split by estimating the censoring survival function and taking the largest time point with adequate censoring support, operationalized in the code as Ĝ(τ ) ≥ 0.001. The IPCW weights themselves were then capped during metric computation, as described below. The Q1 and Q2 horizons were defined from the observed time distribution of the held-out test split. These horizons were used only for reporting, not for model fitting or hyperparameter selection. In datasets with administrative censoring pile-ups, some upper quantile horizons can coincide with the full horizon; this is reported as part of the dataset behavior rather than manually adjusted. F.3
Evaluation Metrics
Time-dependent C-index. We used an IPCW version of the time-dependent C-index to evaluate discrimination under right censoring [44]. Let Ŝi (t) = Ŝθ (t|xi ) be the predicted survival probability for subject i, and let Ĝ(t) = P (C > t) be the censoring survival function estimated by Kaplan-Meier on the training split. Comparable pairs are formed from subjects with observed events before the 20
evaluation horizon. A pair is concordant when the subject who fails earlier is assigned a lower predicted survival probability at that event time. P Cτtd =
i,δi =1,oi <τ
P
j:oi <oj I{Ŝi (oi ) < Ŝj (oi )}wi
P
i,δi =1,oi <τ
P
j:oi <oj wi
,
wi = min{1/Ĝ2 (oi ), 10}.
The cap at 10 was used for all IPCW-based metrics. It is a numerical stabilization step, not a subjectlevel exclusion rule. Subjects were not removed simply because they occurred late in follow-up. Integrated Brier Score. The Brier score measures the squared error between predicted survival and the observed survival status at a fixed time. We used the IPCW form proposed for right-censored survival data [19]. At time t, n
1X BS(t) = n i=1
"
2 Ŝi (t) I(oi < t, δi = 1) Ĝ(oi )
+
# 2 1 − Ŝi (t) I(oi > t) Ĝ(t)
,
with the same IPCW cap applied to the inverse censoring weights. The integrated Brier score (IBS) was computed by numerical integration over the evaluation interval. Lower IBS indicates better calibrated survival prediction. Integrated binomial log-likelihood. We also report an IPCW binomial log-likelihood, which is sensitive to both calibration and discrimination. At time t, " # n 1 X log 1 − Ŝi (t) I(oi < t, δi = 1) log Ŝi (t) I(oi > t) BLL(t) = + . n i=1 Ĝ(oi ) Ĝ(t) The integrated version, denoted IBLL, was computed over the same time grid as IBS. Because this is a log-likelihood rather than a loss, higher values are better. D-calibration. We also report distributional calibration (D-calibration), following [20]. D-calibration assesses whether the model’s predicted individual survival distributions are calibrated over the full follow-up range, rather than at a single fixed time point. Under a correctly specified continuous survival model, the probability integral transform implies that F (T |X) is uniformly distributed on [0, 1]; equivalently, S(T |X) is also uniformly distributed on [0, 1]. Thus, for uncensored subjects, we evaluate the predicted survival probability at the observed event time, Ŝi (oi ). Right-censoring requires a partial contribution because the true event time is only known to exceed oi . On the survival probability scale, the unobserved value Ŝi (Ti ) must lie in the interval [0, Ŝi (oi )]. Following the standard D-calibration construction, each censored subject’s contribution is therefore distributed uniformly over this compatible interval. We divide [0, 1] into 10 equal-width bins and compare the resulting bin counts with the uniform distribution using a chi-square goodness-of-fit statistic with 9 degrees of freedom. We report the corresponding p-value, where larger values indicate less evidence against distributional calibration. Because this is a hypothesis-test diagnostic rather than an effect-size metric, we use D-calibration as a secondary calibration assessment rather than as the primary model-ranking criterion. F.4
Implementation Details
Tabular Data For tabular datasets, all models used a multilayer perceptron (MLP) backbone to extract fixeddimensional representations from input covariates. The MLP consists of multiple hidden layers with nonlinear activations, dropout, and optional batch normalization. Within each experiment, all models used the same backbone search space, followed by model-specific output heads determined by each survival modeling formulation. We performed hyperparameter optimization using random search. Each candidate configuration was trained for up to 200 epochs, and the final configuration was selected based on the best validation Ctd , with negative validation IBS used as a tie-breaker. We used AdamW as the optimizer with a 21
Figure 5: ResNet-18 backbone architecture used for medical imaging survival modeling. The original fully connected layer is removed, and model-specific heads are applied to the extracted feature vectors to produce survival predictions. cosine learning-rate schedule, setting Tmax to the maximum number of epochs and ηmin = 0. For each tabular dataset and model, we evaluated 30 random hyperparameter configurations. The search space is summarized in Table 6. For NnetSurv and DeepHit, the number of discrete time bins was fixed at 50. For MDN, the number of mixture components was fixed at 5. Table 6: Hyperparameter search space for tabular datasets. Hyperparameter
Search space
Number of dense hidden layers Hidden size Learning rate Weight decay Dropout Batch size Batch normalization
{2, 3, 4} {32, 64, 128, 256} LogUniform(10−4 , 10−2 ) LogUniform(10−8 , 10−3 ) {0.0, 0.1, 0.3, 0.5} {64, 128, 256} {True, False}
Medical Imaging Data For medical imaging datasets, all models used ResNet-18 as the convolutional backbone. The original fully connected layer was removed, and the extracted feature vector was passed to a model-specific survival head. The first convolutional layer was modified when the input channel dimension differed from the standard three-channel RGB setting. For single-channel inputs with pretrained weights, the first-layer filters were initialized by averaging the pretrained RGB filters across channels. As shown in Figure 5, all image-based models shared the same ResNet-18 feature extractor, while the final prediction heads differed according to each model’s output parameterization. We performed random-search hyperparameter optimization for medical imaging models using the search space in Table 7. Each candidate configuration was trained for up to 200 epochs, and the final configuration was selected using validation Ctd , with negative validation IBS used as a tie-breaker. We used AdamW with a cosine learning-rate schedule, again setting Tmax to the maximum number of epochs and ηmin = 0. For each imaging dataset and model, we evaluated 20 random hyperparameter configurations. Hidden-layer dimensions were not searched because the feature dimension was fixed by the ResNet-18 backbone. As in the tabular experiments, the number of time bins was fixed at 50 for NnetSurv and DeepHit, and the number of mixture components was fixed at 5 for MDN. Table 7: Hyperparameter search space for medical imaging datasets. Hidden-layer dimensions are not searched because the feature dimension is fixed by the convolutional backbone. Hyperparameter
Search space
Learning rate Weight decay Dropout Batch size
LogUniform(10−5 , 10−2 ) LogUniform(10−8 , 10−3 ) {0.0, 0.1, 0.3, 0.5} {16, 32, 64}
22
F.5
Computational Resources
All experiments are reproducible on standard CPU or GPU hardware. Tabular experiments were run on CPU only and require a single node with at least 24 GB of system memory. Each hyperparameter optimization cell, corresponding to one method, one dataset, and one seed with 30 trials and up to 200 training epochs, takes approximately 5–30 minutes of wall-clock time depending on dataset size. The full tabular benchmark, consisting of 6 datasets, 9 methods, and 20 seeds, requires approximately 270 CPU-hours. Image experiments require a single GPU with at least 20 GB of memory. We used NVIDIA L40S 44 GB and A100 80 GB nodes interchangeably. Each hyperparameter optimization cell, consisting of 20 trials and up to 200 training epochs, takes approximately 15-90 minutes depending on dataset size. The full image benchmark, consisting of 4 datasets, 9 methods, and 5 seeds, requires approximately 90 GPU-hours. The synthetic simulation study can be run on a laptop-class CPU. In our experiments, the full simulation sweep was run on an Apple M4 machine with 10 CPU cores and 32 GB RAM and completed in approximately 6-8 hours. SODEN with adjoint backpropagation dominated the simulation wallclock time, requiring approximately 130 seconds per fit, whereas all other methods completed in approximately 1-7 seconds per fit. Including preliminary experiments, repeated hyperparameter-configuration runs, and debugging runs that were not included in the final reported results, the total compute used for the project is estimated to be approximately 3–5 times the headline compute reported above. F.6
Deep Survival Analysis Models
CoxCC and CoxTime CoxCC and CoxTime extend the Cox proportional hazards framework by replacing the linear predictor with a neural network [32]. CoxCC retains the proportional hazards assumption, λ(t|x) = λ0 (t) exp{g(x)}, (19) where λ0 (t) is the baseline hazard and g(x) is a time-independent neural risk score. CoxTime relaxes this assumption by allowing the risk score to depend on time, λ(t|x) = λ0 (t) exp{g(x, t)}, (20) thereby permitting time-varying covariate effects and crossing survival curves. We used the PyCox implementation: https://github.com/havakv/pycox. To reduce the computational cost of evaluating full risk sets, these implementations use a case-control approximation to the Cox partial likelihood. For a subject i who experiences an event at time Ti , the model compares the event subject against a sampled control set Ji of individuals who remain at risk at Ti . The resulting loss is n 1X 1 X L= log (1 + exp{g(xj ) − g(xi )}) . (21) n i=1 |Ji | j∈Ji
This softplus form penalizes the model when a control subject j is assigned a higher risk score than the event subject i. After training, a baseline hazard estimate is required to obtain absolute survival probabilities. Since the neural network predicts relative risk, the baseline hazard is recovered nonparametrically using the Breslow estimator. For a unique event time tm with dm events and risk set R(tm ), the CoxCC baseline hazard estimate is dm λ̂0 (tm ) = P . (22) j∈R(tm ) exp{g(xj )} For CoxTime, the risk score depends on the event time, so the network must be evaluated at tm for each subject in the risk set: dm . (23) λ̂0 (tm ) = P exp{g(x j , tm )} j∈R(tm ) The cumulative baseline hazard Λ̂0 (t) is then obtained by summing these estimates over event times tm ≤ t. 23
DeepHit DeepHit formulates survival analysis as a discrete-time prediction problem [33]. The time axis is divided into intervals, and the network outputs a probability mass function p = (p1 , . . . , pK ), where pk is the predicted probability that the event occurs in interval k. The cumulative incidence function is k X F (tk |x) = pj . (24) j=1
We implemeted DeepHit based https://github.com/havakv/pycox.
on
the
PyCox
implementation:
The training objective combines a likelihood term and a ranking loss: L = (1 − α)LNLL + αLRank ,
(25)
where α controls the relative weight of the two components. For an uncensored subject whose event occurs in interval k ∗ , the likelihood contribution is − log pk∗ . For a subject censored in interval k ∗ , the likelihood contribution is − log{1 − F (tk∗ |x)}. The ranking loss encourages subjects who experience earlier events to have higher predicted cumulative incidence at the corresponding event time: X F (TA |xB ) − F (TA |xA ) , (26) LRank = I(TA < TB ) exp σ A,B
where σ is a smoothing parameter. NnetSurv NnetSurv models survival in discrete time through interval-specific conditional event probabilities. Let λj (x) denote the probability of an event in interval [tj−1 , tj ) conditional on survival up to tj−1 : λj (x) = P (T ∈ [tj−1 , tj )|T ≥ tj−1 , x) = σ(fj (x)). (27) The corresponding survival probability at tk is S(tk |x) =
k Y
{1 − λj (x)}.
(28)
j=1
The model is trained with a binary cross-entropy loss over all intervals in which each subject is at risk. For subject i with observed interval index ki and event indicator δi , the target is yij = 0 for j < ki and yiki = δi . The loss is L=−
ki n X X
[yij log λj (xi ) + (1 − yij ) log{1 − λj (xi )}] .
(29)
i=1 j=1
Survival Mixture Density Network The Survival Mixture Density Network (MDN) models the time-to-event distribution using a Gaussian mixture distribution on a latent time scale [21]. We implemented MDN based on the authors’ official codebase: https://github.com/XintianHan/Survival-MDN. To ensure positive event times, the observed time is represented as T = softplus(T ∗ ) = log(1 + exp{T ∗ }), where the latent variable T ∗ follows a conditional Gaussian mixture: K X fT ∗ (u|x) = πk (x)N {u; µk (x), σk (x)}. (30) k=1
The density of the observed event time follows from the change-of-variables formula: fT (t|x) = fT ∗ {softplus−1 (t)|x}
d exp(t) softplus−1 (t) = fT ∗ {softplus−1 (t)|x} . (31) dt exp(t) − 1
The survival function is
ST (t|x) = 1 − FT ∗ {softplus−1 (t)|x}, (32) where FT ∗ is the CDF of the latent Gaussian mixture. The model is trained by minimizing the censored negative log-likelihood, n X L=− [δi log fT (ti |xi ) + (1 − δi ) log ST (ti |xi )] . (33) i=1
24
DeSurv DeSurv parameterizes the conditional distribution function through a monotone neural construction [10]. Rather than modeling the hazard directly, DeSurv defines a non-negative function g(t, x) =
∂u(t, x) ∂t
(34)
using a softplus output, and constructs Z t g(s, x)ds.
u(t, x) =
(35)
0
In our implementation, this integral is evaluated using fixed Gauss–Legendre quadrature with 15 nodes. The conditional CDF is then obtained by mapping u(t, x) to [0, 1]: F (t|x) = tanh{u(t, x)}.
(36)
The event density follows by differentiation: f (t|x) =
∂F (t|x) = 1 − F (t|x)2 g(t, x). ∂t
(37)
Thus, f (t|x) = {1 + F (t|x)}g(t, x). (38) S(t|x) The hazard is therefore a derived quantity rather than the primitive parameterized object. DeSurv is trained by minimizing the standard censored negative log-likelihood, S(t|x) = 1 − F (t|x),
λ(t|x) =
n
L=−
1X [δi log f (ti |xi ) + (1 − δi ) log S(ti |xi )] . n i=1
(39)
SODEN SODEN is a continuous-time survival model based on Neural Ordinary Differential Equations [42]. We implemented SODEN using the authors’ codebase: https://github.com/jiaqima/SODEN. The model represents the instantaneous hazard λ(t|x) with a neural network and computes the cumulative hazard through the ODE integral Z t Λ(t|x) = λ(u|x)du. (40) 0
This integral is evaluated during the forward pass using a differentiable adaptive ODE solver. The model is trained by minimizing the censored negative log-likelihood, n
L=−
1X [δi log λ(ti |xi ) − Λ(ti |xi )] . n i=1
(41)
If the input time scale is transformed during preprocessing, the corresponding Jacobian adjustment is included so that the likelihood remains valid on the original time scale.
25
Figure 6: True vs. predicted survival, cumulative hazard, and instantaneous hazard functions for exponential distribution scenario. Shades are 95% confidence interval.
Figure 7: True vs. predicted survival, cumulative hazard, and instantaneous hazard functions for Weibull distribution scenario. Shades are 95% confidence interval.
Figure 8: True vs. predicted survival, cumulative hazard, and instantaneous hazard functions for gamma distribution scenario. Shades are 95% confidence interval.
Figure 9: True vs. predicted survival, cumulative hazard, and instantaneous hazard functions for Gompertz distribution scenario. Shades are 95% confidence interval.
26
Figure 10: True vs. predicted survival, cumulative hazard, and instantaneous hazard functions for log-normal distribution scenario. Shades are 95% confidence interval.
Figure 11: True vs. predicted survival, cumulative hazard, and instantaneous hazard functions for log-logistic distribution scenario. Shades are 95% confidence interval.
27