Saddle-Point Deep Generative Copula Survival Model: Differentiable Likelihood Surrogates for High-Dimensional Censored Data

Full text

Text extracted from the article PDF for reading and search. Tables, figures and equations are in the PDF, which remains the authoritative version.

Introduction

Survival analysis is concerned with the time until an event of interest occurs (which is extremely relevant in many scenes of research, ranging from biomedical research to reliability engineering and econometrics) (Kleinbaum & Klein, 1996) [21]. Censoring, or the event of interest not being observed, can significantly complicate the analysis (Kalbfleisch & Prentice, 2002) [18]. In its most basic form, censoring occurs when the parties of interest are still under observation at the end of a study period. One common example is that the event of interest occurs after the study has concluded; this is called right-censoring (Kalbfleisch & Prentice, 2002) [18]. Well-known classical methods for modeling survival analysis include the Cox Proportional Hazards model (Moolgavkar et al., 2018) [29] and the Accelerated Failure Time model (Saikia & Barman, 2017) [39], both of which rely on a set of assumptions (such as proportional hazards), which may be extremely constraining, negating one of the main advantages of high-dimensional data, namely, the ability to capture complex and non-linear relationships. Recently, deep learning has been increasingly popular and incorporated into survival modeling due to various limitations with traditional survival approaches being surmounted by the flexible hazard functions learnt through Deep Survival Machines (Nagpal et al., 2021) [31] and DeepHit (Lee et al., 2018) from neural nets. Neural ODE (Chen et al., 2018) [9] moves these advances in survival modeling further towards continuous-time dynamics. Yet, we would like to highlight a major problem that has been almost neglected in the context of these approaches, which is the case of multiple, dependent events. Such can include cases of multiple disease progressions with competing risks or in system reliability with more than one component at a time. As such, we would like to conclude the presentation by discussing copula models from the perspective of multivariate dependence.

Copula theory (Nelsen, 2006) [32] solves the problems of dependence through the more tractable separation of a multivariate dependence structure from its marginal distributions. Aas et al. (2009) [1] proposed a convenient framework for constructing high-dimensional copulas from bivariate pairs, known as vine copulas. The difficulty of estimating the copula depends on the collection of associated copula pairs and their corresponding models, which seems difficult to learn. The idea of constructing high-dimensional copulas by sequentially stacking bivariate copulas intuitively corresponds to a special case of a Deep Generative Copula (Janke et al. (2021)) [16]. These deep generative copula models combine the advantages of vine copula modeling with those of normalizing flows (Rezende and Mohamed, 2015) [37] and hence RealNVP (Dinh et al., 2016) to learn highly flexible marginal distributions. One easy avenue of improvement is to take account of censoring by transforming the intractable contribution to the likelihood from integrating unseen failure times into the surrogate likelihood. The only existing approaches to copula-based models for censoring (Glasserman, 2003) [14] rely on Monte Carlo techniques or variational approximations to correspondingly estimate these integrals (Blei et al., 2017) [6], both of which suffer from stochastic estimates of the gradients that can severely destabilize training with an increasing number of dimensions. The variance of these stochastic estimators is another concern because they can slow convergence or even lead to suboptimal solutions. In this work, we sidestep these problems and drastically improve the scalability of the copula model by replacing an approximation of stochastic integrals with a differentiable, deterministic surrogate. This method corresponds to the Sliced-and-Pruned Deep Generative Copula Survival Model (SP-DGCSM).

We will next reduce the integral in the censored likelihood using the saddle-point approximation methods previously developed in Rue et al. (2009) [38] for Laplace approximations and Cam (2012) [8] for steepest descent. These methods are well-known and widely applied in statistical physics as well as Bayesian statistics to produce asymptotic approximations of otherwise intractable integrals. Our goal will be to again construct the full censored log-likelihood as a continuous path integral over the latent failure times, identify the dominant saddle point configuration of the latent failure times under the censoring observations, and apply a second-order Taylor expansion around this saddle point to reduce the integral to a closed-form Gaussian integral. The result will be a fully differentiable surrogate likelihood that depends smoothly on the model parameters and our computed saddle-point value. Regardless of the model parameters, the surrogate likelihood will be continuous and have a well-defined gradient via automatic differentiation, yielding exact gradients for our copula parameters and for the neural-network weights via automatic differentiation rather than stochastic estimators. Prior work on deep generative copulas is limited to either a Monte Carlo or variational approach for obtaining the gradients of the censored copula likelihood. In this work, we propose a new method to compute the censored copula likelihood gradients that avoids sampling entirely, making it possible to bypass both known approaches. To do this, we integrate the gradient surrogate into a new deep generative copula model with a learnable regular vine dependence structure (Joe et al., 2010) defined on top of the copula. In this manner, we make the marginal distributions learnable using a neural hazard network, finally enabling the full end-to-end optimization of the model. To achieve this, we employ the preconditioned gradient descent optimizer from Anil et al. (2020) [3] in high-dimensional parameter spaces, greatly stabilizing the optimization routine and enabling convergence. A significant advantage of our method is that it preserves probabilistic coherence and copula identifiability while demonstrating impressive scalability of hundreds of dimensions and accurate uncertainty quantification and predictive performance.

Literature Review

Copula-Based Dependence Modeling in Survival Analysis The field of multivariate survival analysis, which deals with correlated (possibly multivariate) lifetimes, is quite common in biomedical sciences, reliability engineering, and economics. Most traditional approaches for handling this type of analysis assume that failures of individuals are uncorrelated or rely on frailty models (Wienke 2010) [42] in which correlated failures are seen as a result of a shared random effect. However, the latter modeling imposes one single global structure of dependence across all types of events, which might be quite unrealistic when there are complicated and non-parametric associations between types of events. Copula models (Nelsen 2006) [32] have the advantage of separately modeling the joint association among a multitude of types of events and the marginal survival functions, offering more flexibility by allowing different shapes of associations and, in particular, providing the possibility of asymmetric and varying tail-dependence properties. The Clayton, Gumbel, and Frank copulas are very simple graphical tools often used in typical multivariate survival analysis problems (Othus and Li 2010) [34], but with growing popularity, a more flexible methodology called vine copulas (Aas et al. 2009) [1] can construct high-dimensional dependence structures in a modular sense from a cascade of bivariate copulas. With the advent of vine copulas, conditional dependencies in high-dimensional settings where the number of heterogeneous pairwise interactions grow quadratically can be accounted for. The Dynamic Vine Copula (Nagler et al. 2022) [30] and the Conditional Vine Copula (Wu et al. 2021) extend vine copulas into the time series and censored survival domains, respectively. These models build on prior work but, like other copula models, require pre-specified parametric marginal distributions. Because marginal distributions can be mis specified, the fundamental assumption of parametric marginals may not hold in practice. We overcome this limitation by the proposed form of learning the marginal hazard functions through deep neural networks, circumventing the ill-posed problem of intractable integral likelihoods with a deterministic surrogate.

Deep Generative Models for Survival Data At present, most deep learning-based survival models, such as DeepSurv (Wiegrebe et al., 2024) [41], which implements the Cox proportional hazards model via deep neural network, and DeepHit, directly learn the discrete event time distribution. This relies on implicit rather than explicit dependence between competing events, as information is passed across the network through weight sharing rather than a model of the probability of failure time instants. Generative models step closer to this explicit modeling of the processes that generate failure times. Deep generative MCMC models (Norcliffe et al., 2023) [33] learn a Markov model that generates counterfactual failure times for censored individuals using latent variables. The dependence on expensive MCMC inference along with the computational burden of training are drawbacks of such models in both training and inference stages. The Variational Autoencoders (VAEs)-based survival models (Kingma & Welling, 2013) [20] employ a latent variable that influences the hazard over time, therefore allowing generative modeling in the context of survival analysis (Apellániz et al., 2024) [4]. Existing work leveraging VAEs requires estimating a Kullback–Leibler divergence term from samples, introducing stochastic noise into the gradients. In contrast, our saddle-point approximation provides a deterministic surrogate and hence avoids sampling noise by ensuring differentiability throughout. The proposed method thus enjoys much improved convergence properties. Two recent works propose methods for obtaining survival models that account for censoring using GANs (Goodfellow et al., 2014) [15], one by Zhong et al. (2024) [46] and the other by Norcliffe et al. (2023) [33]. However, the GAN training procedure presents challenges to convergence and requires considerable fine-tuning of hyperparameters. In contrast to these existing methods, we propose estimators that directly optimize a principled surrogate for the marginal likelihood without incurring all the problems associated with min-max optimization.

Integral Approximation Techniques in Deep Learning In probabilistic machine learning, latents are often troublesome; namely, when integrals over them arise, they usually cannot be computed in closed form. Variational inference (Blei et al., 2017) [6] techniques work around the problem by parametrizing the approximate posterior using a form amenable to variational inference, with the ELBO (Qi, 2005) [36] optimized with stochastic gradient descent. However, stochastic optimization can require careful variance-reduction techniques; otherwise, gradients can be too noisy to converge at an effective rate due to the variance induced by sampling. When approximating integrals with Monte Carlo estimation (Glasserman, 2003) [14], the variance of the estimator can scale exponentially with the dimensionality of the latent space—an obstacle when working with high-dimensional censored data. Next, consider the familiar Laplace approximation that is widely used in Bayesian inference (Rue et al., 2009) [38]. In Laplace approximation we approximate an intractable integral by a Gaussian approximation that is determined from the second-order Taylor expansion about the mode. Using this idea in the context of survival analysis with censored latent failure times leads to steepest descent approximations of the appropriate path integral and differentiable approximations of the Laplace-type integrals by forming the determinant of the Hessian of the log likelihood function. This yields a differentiable surrogate likelihood that permits end-to-end gradient-based learning.

Comparison with the Proposed Work Existing deep generative survival models focus on approximating the censored copula likelihood stochastically (based on the Monte Carlo method), with VAEs or GANs used as prominent examples of this approach. Alternatively, frailty models impose simple parametric structures on dependence. In contrast, we develop a deterministic deep generative copula survival model (SP-DGCSM model), which provides a computationally tractable approximation of the censored copula likelihood. This deterministic integral approximation removes sampling noise and allows for more stable gradient estimates, which potentially leads to faster convergence than Monte Carlo sampling or variational approaches. The component responsible for the structural dependence is constructed from a learnable vine copula, which enables flexible modelling of high-dimensional dependence. We thereby show that the proposed SP-DGCSM constitutes a differentiable proxy for the intractable censored copula likelihood. The proposed model can therefore serve as a golden standard for other deep generative survival models. Moreover, our proposed model differs from prior work by Janke et al. (2021) [16] as it exploits vine copulas that are fully learnable from data.

Preliminaries This section provides a concise review of fundamental concepts that will be relevant for what follows. First, the theory of copulas is introduced to facilitate modeling the correlation structures that exist in high-dimensional settings through vine copulas. Second, the survival analysis formalism that is employed is presented.

Copula Theory A copula is a multivariate cumulative distribution function (CDF) with uniform marginal distributions on [0,1] (Nelsen, 2006) [32]. The central result in copula theory is Sklar’s theorem, which states that for any multivariate distribution ?????? with marginal CDFs ??????1,…,??????�, there exists a copula � such that:

The copula C is unique due to the preservation of continuity of the margins, and it thus underlines the dependence structure between the random variables without the influence of the margins. The copula density c is obtained as the derivative of the copula CDF.

Estimation of copula parameters. Once we have chosen a parametric version of the copula in this paper, either the Clayton, Gumbel, or Gaussian copula, we proceed to estimate its parameter(s) by maximum likelihood estimation (MLE). If no censoring took place, the joint density can be expressed in terms of the copula density and the marginals as illustrated in the example.

Where �� is the marginal density of the �-th variable. This factorization is a cornerstone of our model, as it allows us to separate the learning of marginal survival functions from the learning of the dependence structure.

Vine Copulas For high-dimensional data (�≫2), directly parameterizing a full �-dimensional copula is challenging due to the curse of dimensionality. Vine copulas (Aas et al., 2009) [1] circumvent this issue by decomposing a �-dimensional copula density into a product of �(�−1)/2 bivariate copula densities, organized in a hierarchical tree structure. Specifically, a regular vine (R-vine) (Joe et al., 2010) is specified by a sequence of �−1 nested trees �1,…,��−1. Each tree �� has nodes corresponding to a set of variables (or conditional variables) and edges representing bivariate copulas. The density of a regular vine can be expressed as:

where for an edge � connecting nodes �� and �� with conditioning set ��, ��|�� is the bivariate copula density, and the arguments are conditional distribution functions (h-functions). While in practice vine copula structures are often rather intuitively selected, employing measures of dependence such as Kendall’s τ, the present framework directly takes into account the structure learning as well as estimating the bivariate copula parameters. By taking advantage of complete end-to-end gradient descent, the entire model becomes fully adaptable to the data.

Survival Analysis and Censoring The variable of interest in a right-censored survival experiment is the failure time, denoted by T in the following. However, as we will see, not all subjects experience the event before dropping out, and thus we will not observe all failure times. A surviving subject is therefore described by the pair (Y, δ), where Y = min(T, C), C is the censoring time, and δ = I(T ≤ C) is a variable indicating whether the event is observed. The hazard function (Kalbfleisch and Prentice, 2002) [18] is defined as follows. One considers the phenomenon of interest as the outcome of a counting process, and the hazard function, denoted λ(t), is an instantaneous failure rate interpreted conditionally given that the event has not occurred before time t.

0. The

Laplace-Approximated Censored Copula Gradients When faced with right-censored observations, rather than exactly giving the copula density evaluated at the observed censor time, as would be the case if there were no censoring, we instead know some lower bound on the failure time (�� > ??????�). Thus, the contribution of the joint-likelihood of (��) will consist of an integral from the censoring observation (??????�) to infinity, and this integral may not have a tractable closed form. In particular, for more complex copula densities like those that arise from neural networks in higher dimensions, finding the contribution exactly is impossible. Instead, we will treat the censoring integral as a path integral over the latent variables, we will solve it using the method of steepest descent by finding the most probable saddle point, and finally we will use a Gaussian integral expansion of that saddle point. The result is a deterministic differentiable approximation.

Saddle-Point Location via Steepest Descent The exponent in Equation 10 is:

The Newton update for the censored variables, holding the uncensored variables fixed, is:

Gaussian Integral Approximation and Differentiable Surrogate Once the saddle point �∗ is found, we perform a second-order Taylor expansion of the exponent around this point. The Hessian of � is:

The Hessian is negative-definite at the saddle point, as required. The integral in Equation 10 is then approximated by a Gaussian integral:

Taking the logarithm, we obtain the surrogate log-likelihood for the observation:

The term ??????�/??????�∗ is zero by the definition of the saddle point. The gradient of the log-determinant can be computed using standard matrix calculus (Magnus & Neudecker, 2019) [28]. This yields a deterministic, unbiased gradient estimator.

Recursive Saddle-Point Approximation for Vine Copulas Another and more structured way of deriving the saddle-point approximation is to note that the density of a vine copula has a hierarchical structure that can likewise be derived recursively. For example, to arrive at a learnable vine copula like that just described, the tree-structured architecture could be derived starting from the joint density based on a product of bivariate copulas ���|� which can be computed conditional on the information in the observed variables D. More precisely, two latent failure times can be linked to the observed information component of the density in order to allow each integral in the saddle-point factorization. Consider a bivariate pair (�,�)∈� where both components are censored. The integral becomes:

where � is the conditional correlation parameter. This closed form for the one-dimensional case is computationally efficient and avoids iterative solvers for each pair.

Integration into The Deep Generative Vine Copula Module As illustrated in Figure 1, each component's marginal distribution is modeled by a neural hazard network. Here, a covariate � (which may be multivariate) is input into the neural network, producing a parametric baseline hazard for flexible distributions, such as a Weibull, or a piecewise constant baseline. Dependencies between components can be flexibly represented by a learnable vine copula structure. Each bivariate copula can be parameterized with a neural network, which receives input � and produces a scalar correlation parameter �. Finally, the parameters of the model are trained through a negative surrogate log-likelihood.

Where � is the number of observations. Fig 1: Architecture of the Saddle-Point Deep Generative Copula Survival Model Lastly, at every iteration, we compute gradients both from the saddle-point solver and directly from the log-determinant term using automatic differentiation. We then feed them into a preconditioned optimizer like Shampoo (Anil et al., 2020) [3] with an explicit Hessian structure from the saddle-point approximation, i.e., the term H in equation (16) acts as a natural preconditioner for copula parameters. In each batch, we solve for z* using Newton’s method, which typically requires 5–10 iterations for convergence, then evaluate the surrogate loss, backpropagate its gradients, and finally update the weights of the neural network using the preconditioned gradient.

Deterministic Uncertainty Quantification during Inference At inference time, we now choose the new subject’s covariates, �∗, and observed times, �∗, to compute actual probabilities. The goal at inference, then, is simply to compute specific conditional survival probabilities, namely, the joint probability of survival at future times, �∗, with the most likely latent configuration �∗ chosen for the latent variables, and conditional on whatever values we observed for the censored component (component j in this case).

where �� is the vector of correlations between component � and the observed components, and ��∗ are the saddle-point values for the observed components at their observed times. This expression provides a closed-form, sampling-free estimate of the conditional survival probability, which is essential for fast and reliable uncertainty quantification in clinical decision support and reliability analysis.

Experimental Results We validate the SP-DGCSM on both synthetic and real high-dimensional survival datasets, focusing on a few objectives. The first is to empirically assess the speedup achieved by our method over the standard maximum likelihood estimation in terms of both stability and scalability by comparing the deterministic saddle-point surrogate objective with the Monte Carlo-based approximation of the maximum likelihood estimate. The second is to compare the predictive performance and calibration of our approach to existing state-of-the-art survival analysis baselines. Our experiments are performed in PyTorch on a single NVIDIA A100 GPU with 40 GB of memory, and we set the Newton solver for finding the saddle-point location to at most 10 iterations with a tolerance of 10−6 for �∗ with respect to relative change.

Experimental Setup We consider four datasets with gradually increasing dimension, censoring rate, and complexity of dependence.

Synthetic High-Dimensional Vine Copula Data The first simulation study identifies the dependence structure and hazard function with censored data. First, we want to see whether we can identify the correct dependence structure of the data. To this end, we simulate a data set of size 10,000 with the following characteristics. For simulation of the multivariate random variable, we first simulate a dendrogram from a regular vine copula, for which we randomly assign the edges one of the three pair-copulas, Clayton, Gumbel, and Gaussian. We then simulate the marginal distributions by Weibull distributions with a shape parameter uniformly drawn from the interval (0.5, 2.0) and a scale parameter fixed to 1.0. These marginal distributions are then combined with the random variable to realize the data set. For the introduction of censoring effects, we simulate an independent censoring random variable from an exponential distribution, controlling the rate parameter so that approximately 40 percent of the data gets censored. Finally, we randomly split the resulting data set into two training sets (7,000 subjects) and one evaluation/test set (1,500 subjects).

SEER Cancer Incidence Data The data set used in this application is the Surveillance, Epidemiology, and End Results (SEER) dataset from 2011. This dataset contains over 200,000 cancer patients that span 12 cancer sites and various demographic subgroups. The survival times from these cancer patients are modeled as a joint distribution over cancer sites. All combinations of clinical and demographic covariates are included in the marginal hazard networks. The survival times are subject to right censoring and have a censoring rate of approximately 55%.

MIMIC-III Clinical Database For an application, we used data from the MIMIC-III ICU database (Li & Li, 2022). We modeled each of 20 common laboratory measurements and vital signs as separate survival endpoints. Events were defined as either discharge or death because of an unrelated cause, and static covariates included age, gender, and admission type. We considered a total of approximately 40,000 ICU stays, with a mean censoring proportion of 63% across all 20 endpoints.

High-Dimensional Synthetic Stress Test Finally, we perform a scalability experiment of PAC with respect to increasing dimensions of the survival problem at hand. To achieve that, we construct synthetic survival data sets of increasing dimensionalities �∈ {50,100,200,300,500}. We accomplish this using a Gaussian copula of dimension d with an autoregressive correlation structure for the dimension d components, that is ���=0.7|�−�|, together with log-normal marginals calibrated to have 30% censoring over the training set. Each of the five configurations has 5,000 observed training units, and they allow isolating the effect of increasing the dimensionality on both saddle-point convergence and quality of the gradients.

Baselines We compare SP-DGCSM against the following methods: ● DeepSurv (Wiegrebe et al., 2024) [41]: A feedforward neural network, trained using the Cox partial likelihood independently for each dimension, can serve as a non-multivariate benchmark for the performance of univariate applications. ● DeepHit (Lee et al., 2018): We applied this model to obtain the joint distribution of the competing events. To extend the discrete-time model to the multivariate setting, we treat each dimension as a separate competing risk, using a single representation across dimensions. ● Copula-GLM (Othus & Li, 2010) [34]: In line with existing literature, we use a parametric copula model anchored by Weibull marginals, which allows us to fit the model via two-stage maximum likelihood; it is a multivariate Gaussian copula that fixes the correlation matrix at the training-data estimates. ● Deep Generative Copula (DGC) (Janke et al., 2021) [16]: The baseline method corresponds to (12) but with neural marginals and the copula defined using normalizing flows. The censored likelihood is approximated using importance sampling, which means 500 Monte Carlo samples are used per observation. This corresponds to the direct baseline from section 3. ● VAE-Surv (Apellániz et al., 2024) [4]: The proposed encoder-decoder architecture is referred to as variational autoencoder (VAE) since the latent codes are used to represent shared information across all dimensions, whereas the marginal hazards are decoded by different neural networks. The evidence lower bound (ELBO) is maximized by stochastic gradient descent using 100 Monte Carlo samples per observation.

Evaluation Metrics We report the following metrics across all experiments: ● Concordance Index (C-index) (Krishnapuram et al., 2007) [23]: Rank correlation coefficients between the predicted and observed survival times are obtained, evaluated separately for each dimension, averaged, and used as a measure of discriminant validity. ● Integrated Brier Score (IBS) (Fernandez et al., 2024) [12]: A measure of prediction accuracy and calibration over a time horizon, with lower values indicating better performance. ● Negative Log-Likelihood (NLL) on the held-out test set: The copula-based models (SP-DGCSM, DGC and Copula-GLM) all make use of the joint likelihood, which assesses the fit for both the marginal and dependence structure, whereas the non-joint models use the summed marginal negative log-likelihoods.

Kendall’s ?????? Estimation Error To evaluate the accuracy of the recovery of the dependence structures, we used synthetic datasets with nested hierarchies of dependencies. The mean absolute error between the true values and those learned as Kendall’s τ are reported. Main Results Table 1 contains comparative results of SP-DGCSM across all datasets. SP-DGCSM consistently ranks first or second across the different measures and datasets, and its best results are observed on the simulated high-dimensional stress test, as well as the SEER dataset with severe censoring.

Fig 2: Convergence of the deterministic Newton solver to the saddle point overlaid on the exact log-likelihood contours for a bivariate censored case, showing the elliptical approximation of the Hessian at convergence.

We also observe that the Taylor expansion is very accurate for estimating the log-likelihood near the saddle point, while the elliptical Gaussian resembles the true local log-likelihood shape quite closely, which is congruous and indicates that the approximation is accurate. We apply SP-DGCSM on the SEER dataset, achieving a C-index of 0.729 and an integrated Brier score of 0.171 over 12 kinds of cancers. This dataset has heterogeneous hazards associated with good correlation across cancer types. Algorithms like DeepSurv and DeepHit that independently fit a baseline for each cancer type are unable to take advantage of the shared information across the different types of cancer. DGC also gives good performance but suffers from high variance in the gradient estimates, which necessitates a smaller learning rate and results in around 2.3× more wall-clock training time than SP-DGCSM. VAE-Surv also gave a satisfactory prediction but suffers from optimizing the stochastic ELBO due to highly censored tails. In our implementation, we found the Newton iteration usually converges within 5-7 steps.

Despite the potentially extreme sparsity of sample observations for the target events, the DGC’s tendency for gradients to remain stable despite heavily censored observations stubbornly persists, as shown in the MIMIC-III example for prediction of rare and heavily censored endpoints. SP-DGCSM achieves a C-index of 0.694 and an IBS of 0.187 when trained with the physiological measurement dependencies across the 20 endpoints. These features of stability proved essential for the MIMIC-III dataset, but what about the synthetic stress test of 500 dimensions? Here the Monte Carlo variance in DGC proved large, and the number of samples needed was beyond what could be reasonably handled in standard training. Ultimately, we found that more than 2,000 samples per observation were necessary to stabilize from high variance at the gradients through Monte Carlo. In SP-DGCSM this problem is elegantly solved through the recursive saddle-point approximation of the vine copula, where utilizing the closed-form solution illustrated in Equation 22 of Section 3 for the univariate latent integrals significantly reduces the complexity. The resulting training times were similar to the saddle-point inverting 6.3 (average across observations) Newton iterations. Finally, Figure 3 gives an example of the final 500-dimensional dependence structure with pairwise conditional Kendall’s τ errors shown in the heatmap plot. Fig 3: Absolute estimation errors of pairwise conditional Kendall’s tau parameters between the learned vine copula structure and the known ground truth in a 500-dimensional simulation.

Figure 3 suggests that this is indeed the case, as the autoregressive correlation across these dimensions is recovered exactly. The mean absolute error across all (500 2) pairs is only 0.042. This is primarily a result of the characteristic difficulty of separating strong signal from noise in the case of weak dependencies, since all the errors occur on pairs with small strength of association, i.e., true |??????|<

0.1. It should be noted, however, that the vine structure does

successfully identify the hierarchy in the importance of variables, since edges with maximum pair-copulas (the edges of the selected maximum spanning tree) revolve around pairs with equally high pairwise dependencies.

Ablation Study To isolate the contributions of individual components in SP-DGCSM, we conduct an ablation study on the synthetic vine copula dataset (&#3627408465;=100). Table 2 reports the results of five model variants:

Table 2: An ablation study on the synthetic vine copula dataset (&#3627408465;=100) evaluating the impact of key design choices: saddle-point surrogate, learnable vine, and preconditioned optimization. Model Variant C-index (↑) IBS (↓) NLL (↓) Training Time (s/epoch) Full SP-DGCSM 0.782 0.143 2.645 48.2 SP-DGCSM w/ Fixed R-Vine 0.769 0.151 2.718 45.1 SP-DGCSM w/ MC Integration (1000 samples) 0.774 0.149 2.683 72.4 SP-DGCSM w/o Preconditioner 0.771 0.154 2.712 68.9 SP-DGCSM w/ Full Gaussian Copula 0.765 0.156 2.746 41.7 Fixed R-Vine Structure While the vine structure can be learned from the data at hand, i.e., our learnable vine, for demonstration we also fit the data with the vine structure determined at a pilot study using training data. This results in a fixed R-vine determined with regard to marginal Kendall’s τ. We denote this vine of structure R-set. The evaluations for this reduced model yield a C-index of 0.769 and an NLL of 2.718. Thus, there is a substantial gain in jointly learning the vine structure and the parameters of the pair-copula factorization as demonstrated by model comparison through a two-stage fitting approach where the vine structure is fixed. There are dependencies that do not admit of monotonic relationships. A learnable vine structure is therefore attractive because it can adapt to the data-generating process.

Monte Carlo Integration We now replace the saddle-point approximation with Monte Carlo integration, taking 1000 importance samples per observation. The performance degrades to C-index 0.774 and NLL 2.683. Training takes 50% longer, from 48.2 to 72.4 seconds per epoch. The problem is that the variance of the gradients is much greater, giving an oscillating behavior of the loss curve (see Figure 4). The deterministic surrogate is much better. It delivers better accuracy and trains faster and more smoothly.

Fig 4: Training curves of the negative log-likelihood for the saddle-point surrogate versus Monte Carlo integration with 1,000 samples, highlighting the variance reduction and faster convergence of the deterministic method.

Preconditioned Optimizer To substantiate this claim, Table 2 compares Shampoo with Adam (Kingma & Ba, 2014) [19]. Adam converges more slowly, yielding a C-index of 0.771 and an NLL of 2.712. The reason lies in the block-diagonal structure of the Hessian, which is exploited by Shampoo owing to the saddle-point approximation jumping to conclusions regarding the parameters of copula correlation. Without preconditioning, optimization requires many more epochs and runs a greater risk of being trapped in inferior local minima.

Full Gaussian Copula The vine model is not only better at estimating joint likelihoods but also at estimating conditional dependencies. This is not just a nice property but rather the main reason for the large C-index that the vine model achieves compared with the significantly poorer C-index of 0.765 that a single 100- dimensional Gaussian copula achieves 41.7 seconds per epoch, yielding a significant speed advantage because no recursive saddle-point calculations are required.

Calibration and Uncertainty Quantification The SP-DGCSM is deterministic and, as such, can be employed for solving uncertainty quantification problems. Calibration plots obtained from the method proposed by Arrieta-Ibarra et al. (2022) [5] are used here to evaluate the joint survival predictions of the SP-DGCSM. In these plots, the predicted values of the joint survival probabilities are compared to the empirical ones provided via the Kaplan– Meier approach on the test set. Figure 5 shows these plots with DGC and SP-DGCSM regarding MIMIC-III. Focus here is placed in the tails of the distribution, especially for survival probabilities less than 0.3.

Fig 5: Joint survival calibration curves comparing predicted probabilities from the saddle-point surrogate against empirical estimates, emphasizing performance in the distributional tails where Monte Carlo methods exhibit high variance.

The calibration of the SP-DGCSM is excellent over the entire range of survival probabilities, its maximum absolute deviation from the empirical curve equal to 0.028, while that of the DGC model based on Monte Carlo is impaired at the tails by deviations larger than 0.05 for survival probabilities less than 0.2. This loss of calibration occurs in extreme quantiles because of noisy estimates of the Monte Carlo estimator associated with high variance in the importance sampling weights. On the other hand, the calibration for rare events predictions of the SP-DGCSM is perfect, as illustrated by the smoothness of the curve since the surrogate saddle-point construction is entirely deterministic. In all of the datasets we have investigated, the empirical coverage of the 95 percent uncertainty intervals for the conditional probabilities of surviving that we computed via equation (24) at the median time point is greater than 93 percent, which gives empirical support for the Gaussian approximation via the saddle-point Hessian that we relied upon to obtain the uncertainty intervals.

Discussion

The experimental evidence provided in the previous section enriches our understanding of the SP-DGCSM framework in several ways. Besides addressing the question of non-uniqueness when the vine copula relies on neural-parameterized marginal inputs, it allows us to investigate the saddle-point. approximation’s theoretical validity under various censoring schemes, as well as to speculate upon its extensibility to increasingly complex data structures.

Validity Boundaries and Failure Modes of the Saddle-Point Approximation While the surrogate performs well empirically, its validity should be assessed under theoretical conditions. Rue et al. (2009) [38] discusses the conditions under which the Laplace approximation is valid. When the integrand in Equation 10 is sharply concentrated around a single dominant saddle point, the approximation is valid. Thus, the surrogate is likely to perform better when the number of observed, uncensored components is large compared to the number of censored components and when the correlation matrix R is well-conditioned (i.e., not too ill-conditioned). In turn, we should consider cases where the aforementioned conditions are not met, leading to the surrogate failing: first in the scenario of weak concentration, then the case of two competing saddle points, and finally when censoring is relatively extensive.

Nevertheless, we identify three specific scenarios where the approximation may degrade:

1. Extreme Censoring Regimes: Highly censored data

with early failure observations (i.e., &#3627408486;&#3627408471; is small) leads to high-dimensional and unbounded (W) with a slowly converging integrand. In this situation, the Hessian (H), like the Gaussian approximation itself, may become numerically ill-conditioned, which would make the approximation even less reliable. A simple diagnostic for problems like these is the condition number, which becomes problematic when it exceeds about 104, because it leads to numerical instability in the computation of the log-determinant term. Keeping track of these condition numbers during training is advised, and if numerical instability is observed, the safest course is a more conservative Laplace approximation, using the tail correction for censored components whose observation times are early.

2. Multimodal Posterior Distributions. Because the

saddle-point approximation is based on the assumption that the posterior distribution is unimodal, we have sketched some methods that can be used when highly negative correlations among the components are expected so that a multihumped posterior distribution is likely. This is particularly problematic in competing risks models, where one of the events is censored. It is known that in these circumstances there can be an interpretation of the estimates of the latent time points as being bimodal; one mode (or peak) corresponds to the events occurring quite soon after the observation corresponding to the censored event, while the other mode has the feature of the events never occurring at all up until quite a late date (with the censoring occurring then). It is argued that the saddle-point approximation presumes only one of these peaks to exist in order to be useful, so that it cannot approximate accurately the total mass of the probability distribution. To cope with this situation, a mixture of saddle-point approximations (Crawford, 1994) [10] is proposed, which corresponds to making a sum of Gaussian integrals at all the stationary points where the mass of the distribution is located.

3. Non-Gaussian Copula Families. Density equations (7)

and (9) are derived within Section 4 under the assumption of Gaussian bivariate pair-copulas. The conceptual hurdle involved in deriving saddle-point equations in the presence of non-Gaussian pair-copulas like the Clayton or Gumbel copulas, is that one would have to deal with the more complex function g(z) in equation (11) and with the fact that no unique saddle-point is guaranteed. In that regard, one alternative is to closely approximate the Archimedean copulas with Gaussian copulas after rank transformations (see Ling et al. (2020)) [26]. The other is to work out the saddle-point approximations in uniform or observable space rather than in Gaussian latent space.

Solving the Newton updates is the most taxing part of the optimizations, due to the need to invert the |C| matrix present in Equation 15. Because the cost scales as &#3627408450;(|&#3627408438;|3) this is okay for the hundreds of parameters we have encountered already, but it becomes infeasible for even thousands of parameters. A scalable solution is to approximate &#3627408453;−1 via the vine structure: all updates to a latent variable depend on a block-diagonal structure for which individual updates can be solved independently as bivariate or trivariate regression problems. A caveat of this structure is that it only holds if the pairs in our copula are conditionally independent, which is quite an assumption. A better alternative is to run quasi-Newton optimizations, such as L-BFGS (Liu & Nocedal 1989) [27]. The key advantage here is that no Hessian needs to be computed or inverted, meaning that the cost scales linearly with the number of latent variables.

Interpretability Challenges in Neural-Parameterized Vine Structures The learnable vine copula in SP-DGCSM is parameterized by neural networks that output the conditional correlation parameter &#3627409164;&#3627408462;&#3627408463;|&#3627408439;(&#3627408485;) for each pair-copula edge as a function of the covariates &#3627408485;. While this provides tremendous flexibility, it also introduces significant challenges for model interpretation. In classical vine copula models, each pair-copula has a fixed, parametric form (e.g., Clayton with a single parameter), and the resulting dependence structure is directly interpretable in terms of bivariate conditional dependencies. In our model, each &#3627409164;&#3627408462;&#3627408463;|&#3627408439; is the output of a multi-layer perceptron, making it difficult to trace how a change in a specific covariate affects the dependence between two failure times.

Previously, our neural network that predicts the correlation parameter of the pair-copula depended on all of our covariates. However, it is often the case where we want to understand how certain covariates modulate dependence between certain endpoints. For instance, we may expect the covariate age or primary site to account for the apparent association between breast cancer survival and ovarian cancer survival, while other covariates do not play a role.

Thus, it is desirable to have a model where each pair-copula depends on a subset of covariates. This can be accomplished using a structured prior on the neural network weights, thus inducing a type of sparsity, or by replacing the neural network with a generalized additive model (GAM). With a GAM each covariate’s effect is decomposed into smooth and interpretable function components. A drawback of using the GAM approach in place of a neural network, however, may be that substantially non-linear modulations of dependence across covariates are less well represented; nonetheless, this “interpretability vs. capacity bias” trade-off may be worthwhile in some applications such as clinical decision support. Currently, the implementation learns the vine structure from scratch in each epoch via maximum spanning trees on pairs of dependence measures and dynamically updates the tree with each training iteration stymies effective analysis of changes in the conditional dependence graph across training and patient subpopulations. One potential extension of these ideas is to adopt a fully Bayesian approach to the vine structure that would marginalize over its uncertainties during prediction, resulting in a posterior distribution over edges that would amount to a principled quantity measuring the strength of the dependence. However, this raises more difficult questions of model fitting, since a Bayesian treatment would require computation of integrals over a discrete space of regular vines—a combinatorial optimization problem that is not trivial to solve. Efficient methods for sampling from the space of regular vines remain a challenge worth pursuing. While the advantage of the deterministic saddle-point approximations is computational efficiency, the drawback is that these formulations obfuscate the stochastic nature of the relevant epistemic uncertainty via the omission of its random component. In a full Bayesian treatment, the posterior variance across latent failure times would be captured only through sampling, while the deterministic implementation discussed in this chapter approximates this variation through the Hessian-based Gaussian approximation (Equation 18). The approximated uncertainty interval was well-calibrated according to the examples shown; however, more work is warranted to measure the interval coverage in relation to samples from the full posterior distribution produced via a Markov chain Monte Carlo simulator, especially in the case of model misspecification, which may lead to systematic under coverage.

Extensions to Complex Censoring Regimes and Domain-Specific Applications The SP-DGCSM we have discussed here assumes, for simplicity, that observations are right-censored with censoring times distributed independently of event times. However, many real-life situations lead to left-censoring, interval-censoring, and possibly competing terminal events. These patterns are summarized by Pintilie (2006) [35], and our next challenge is to extend the SP-DGCSM to deal with them through the rich structure provided by the saddle-point approximation. The constrained saddle-point equations can be solved iteratively, yielding a computationally feasible way to evaluate the Laplace approximation. When failure times are interval-censored, the variables for which this occurs have lower and upper bounds, i.e., [&#3627408484;&#3627408471;&#3627408474;&#3627408470;&#3627408475;,&#3627408484;&#3627408471;&#3627408474;&#3627408462;&#3627408485;] and [&#3627408473;&#3627408471;,&#3627408479;&#3627408471;], and simple changes of variables yield to hyper-rectangular integration domains. This is not a problem for a truncated Gaussian copula, where one can simply integrate over the truncated multi-normal distribution. But the saddle-point equations must be solved subject to the box constraints imposed by the interval censoring. Fortunately, the saddle-point equations can be solved without the constraints, and the approximation is simply clipped to remain feasible (possibly in the saddle-point equations themselves). Wongkamthong and Akande (2023) [43] employ variations on this clipping technique to handle imputation of ordinal data, and the same techniques may be effective here. The saddle point equations begin as in Eq 12 and then start clipping the approximation to the feasible box. When clipping leads to a boundary solution, the constrained saddle point will lie on the boundary of the box constraint, yielding a lower-dimensioned integral over the unconstrained components only.

The class of problems one might have in mind when seeking the type of augmented saddle-point approximations considered in this paper would be situations arising in competing risks. In such problems only the minimum observed time of failure is known, together with its cause; the failure times associated with the other causes remain hidden. The likelihood in such cases must integrate over all the associated times for the different causes. Fortunately, a multivariate model for survival times provides a natural framework since one can view the dependent failure times in this situation. The interested reader should follow recently developed models based on vine copula approaches, which enable one to model the joint distribution of such dependent failure times without needing to assume independence of causes. The saddle-point approximation can then be applied to the latent times associated with the hidden causes. Such a model may well outperform models based on independent causes that dominate much of the literature (Pintilie, 2006) [35]. The classic situation considered so far is that covariates are fixed throughout an individual’s lifetime. The SP-DGCSM can be extended to incorporate time-varying covariates and models of longitudinal survival data. For instance, the hazards associated with time-dependent biomarkers can be modeled through a neural hazard net (Kleinbaum & Klein, 1996) [21], or a piecewise constant covariate that is measured only at clinic visits can be modeled with a piecewise linear cumulative hazard, and the saddle point solver can then be used independently within each segment.

When multiple biomarkers are observed, one can model them as a multivariate stochastic process with copulas capturing contemporaneous dependencies. One popular strategy for modeling such multivariate time-dependent survival patterns is through dynamic vine copulas (Nagler et al., 2022) [30], where correlation parameters associated with each pair copula will evolve according to observed histories. Again, the saddle-point approximation could be employed at each time step to support differentiable learning of the resulting dynamic survival dependencies. The SP-DGCSM can be applied to numerous survival-analysis settings, including in oncology and critical-care medicine (such as was done for SEER and MIMIC-III), reliability engineering, finance, and demography. The reliability-engineering variant of the model estimates the joint failure times of an entire system’s components based on partial observations of system failures. In a finance context, the model specifies the joint default time of a portfolio of firms, a classic problem originating in the field of credit risk (Altman et al., 1998) [2]. In demography, versions of the model may be used in the study of the joint survival patterns of families and/or populations, accommodating genetic and environmental correlations. In any of these high-dimensional, consequential applications, the SP-DGCSM’s scalable and deterministic saddle-point approximation ensures computational feasibility and prediction accuracy.

Conclusion

Lastly, we propose the Saddle-Point Deep Generative Copula Survival Model (SP-DGCSM), a novel approach that addresses the critical challenge in learning high-dimensional dependence structures, particularly for right-censored survival data. This principled procedure, based on SP-DGCSM, relies on a deterministic saddle-point approximation of the surrogate likelihood, instead of a standard stochastic integral approximation. The censored likelihood can be rewritten as a path integral (Cameron-Martin theorem) over latent variables uniformly distributed in the unit sphere with unit exponential distributions. We then identify the saddle point, yielding the dominant contribution through the steepest descent method. A second-order Taylor expansion at the saddle point leads to a closed-form expression of the Gaussian integral. Automatic differentiation is then applied to compute the resulting exact gradient, yielding superior stability for gradient-based optimization, leading to efficient learning even in high-dimensional spaces. Remarkably, the recursive saddle-point approximations serve to greatly mitigate the bivariate pair-copula problem by facilitating the computation of the latent integrals over one dimension with closed-from expressions. This architecture jointly learns both the neural marginal hazards and the vine correlation parameters in an end-to-end manner, with a preconditioned optimizer making use of the Hessians to further enhance convergence speed. The performance on synthetic and real-world examples (including a 500-dimensional stress test) typically exceeds the state-of-the-art in terms of concordance, calibration, and joint likelihood, while being suitable for uncertainty quantifications by virtue of the determinism in the inference—two-dimensional calibration curves are almost indistinguishable from the empirical estimates in the extreme quantile regions. The SP-DCGM provides an elegant solution to common survival analysis challenges, including handling right-censoring and discovering latent structures for predictive variables, all while avoiding ad hoc sampling from distributions. These models elegantly combine rich statistical capabilities for fitting arbitrary dependence structures with the computational tractability often desired in high-dimensional applications. Future work could focus on implementing the SP-DCGM with other common scenarios in survival analysis, such as interval-censoring, competing risks, and dynamic covariates. We also envision future research aimed at exploring interpretable neural-parameterized dependence structures.