Bayesian Elastic‑Net Cox Models for Time‑to‑Event Prediction: Application with Breast‑Cancer Cohort

1. Introduction

Accurate time–to–event prediction guides decisions on therapy, follow-up, and trial eligibility. Modern molecular profiling yields wide, correlated predictor sets for comparatively modest cohorts, and in this high-dimensional regime Cox regression can become unstable and risk scores can miscalibrate. Penalised variants such as lasso, ridge, and the elastic–net stabilise estimation and encourage structure, yet they return only point estimates and depend on tuning, leaving uncertainty quantification and calibration largely ad hoc. This motivates a fully Bayesian formulation that delivers calibrated risk and honest uncertainty while respecting correlation among predictors.

The classical model for such time–to–event problems is the Cox proportional–hazards regression [1]. In its original form the model works well when the number of predictors p is small or at most moderate. Modern gene–expression studies, however, typically start from tens of thousands of microarray gene-level features measured on only a few thousand patients. In METABRIC, for example, the Illumina HT-12 v3 platform interrogates roughly $p\approx 24,000$ gene features, but the publicly available cBioPortal file that we analyse here already provides gene–level summaries rather than raw feature intensities. After removing genes with missing values in this file we obtain $p_{\mathrm{raw}}=489$ gene-expression variables, and after the 10th–percentile variance filter we work with $p=440$ features (see Section 3). Thus, although the original assay is ultra–high–dimensional at the feature level, the resulting design matrix is still high–dimensional and strongly correlated relative to n, which makes the partial–likelihood surface relatively flat in many directions and the Cox estimates unstable.

One solution is to add penalties. The Lasso shrinks many coefficients exactly to zero, while ridge regression keeps all coefficients but pulls them toward zero. The elastic–net (EN) combines the two ideas and tends to keep groups of correlated genes together [2]. In addition, penalised Cox models have been widely used for high-dimensional survival and gene-expression data, including lasso and ridge-type penalties [3], as well as comparative studies and cross-validated prognostic pipelines for microarray cohorts [4]. More recently, post-shrinkage strategies have been proposed for high-dimensional Cox models [5]. Related to ridge-type regularisation, weighted ridge ideas and post-selection shrinkage strategies based on weighted ridge have also been developed; see Gao et al. [6] and the follow-up line of work on shrinkage and penalty estimation under censoring (see Hossain and Ahmed [7], Ahmed et al. [8]). Still, the elastic–net needs cross–validation to pick two tuning numbers and only returns point estimates, so we do not know how much uncertainty we have in each hazard ratio.

A key statistical motivation behind shrinkage is the bias–variance trade-off: in moderate-to-high dimensions, shrinking noisy estimates toward a central value can reduce mean squared error (see [9] for a modern post-shrinkage view in high-dimensional modeling). This is not only a heuristic idea: the Stein-type shrinkage approach shows that, for multivariate normal means in dimension $p\ge 3$, the usual unbiased estimator can be inadmissible and can be improved by shrinkage estimators. The same principle underlies many regression shrinkage methods. In particular, ridge regression can be interpreted as a Bayes estimator under a Gaussian prior, and many Stein-type and empirical Bayes procedures choose the shrinkage intensity by estimating hyperparameters from data. In this sense, modern penalised survival models sit very close to empirical Bayes thinking, even when they are presented as optimisation-based methods. The BEN–Cox model developed here follows this idea directly: the ${\ell }_1/{\ell }_2$ shrinkage strengths are treated as unknown and learned from the data via hyperpriors, providing a principled Bayesian analogue of data-driven shrinkage.

Bayesian thinking treats a penalty as the negative log of a prior. A Laplace prior gives the Bayesian Lasso [10], and mixing a Laplace with a Gaussian prior gives the Bayesian elastic–net (BEN) [11,12]. Putting hyper–priors on the mixing scales lets the data decide how much ${\ell }_1$ and ${\ell }_2$ shrinkage to apply, and posterior samples provide full credible intervals. Recent theory also shows that such scale–mixture priors contract toward the truth at near-optimal rates under sparsity [13].

In this paper we fit a fully Bayesian EN Cox model to the METABRIC cohort ($n=1,903$, $p=440$ gene–level features derived from an original $\sim 24,000$–probe array). Our aims are to

• show that Hamiltonian Monte Carlo can draw stable posterior samples in a correlated, high–dimensional survival setting,
• compare predictive accuracy, discrimination, and calibration against a tuned ridge Cox baseline, and
• check that the genes kept by the model contain known biology, for example the PAM50 signature first reported in Curtis et al. [14],
• and, importantly, provide theoretical justification for the BEN–Cox in which we formalize a Bayesian elastic–net grouping effect and a posterior contraction result under sparsity (see Section 2.5 and Appendix A).

The rest of this article is organised as follows. Section 2 introduces the proposed BEN–Cox model, prior specification, and theoretical properties. Section 3 describes the METABRIC data set and the pre-processing procedure. Section 4 presents the application results, including predictive performance metrics and calibration diagnostics based on Demler et al. [15]. Section 5 discusses implications and limitations of the proposed method, and Section 6 provides concluding remarks.

2. Methodology

2.1. Model Specification

For patient $i\in \{1,\dots ,n\}$, let $(y_i,{\delta }_i,{\mathbf{x}}_i)$ denote the observed data, where $y_i$ represents the follow-up time, ${\delta }_i$ is the event indicator (${\delta }_i=1$ if the event is fully observed and ${\delta }_i=0$ if the event time is right-censored), and ${\mathbf{x}}_i$ is a p-dimensional vector of z-standardized ${log}_2$-transformed gene expression feature intensities. Here, ${log}_2$ denotes the base-2 logarithm, corresponding to the conventional ${log}_2$ transformation applied to raw expression intensity measurements. The Cox proportional hazards model [1] assumes

$$h(t\mid {\mathbf{x}}_i)=h_0\left(t\right)exp\left({\mathbf{x}}_i^{\mathsf{T}}\beta \right),$$

(1)

with unspecified baseline hazard $h_0\left(t\right)$. Inference uses the log partial likelihood

$$\ell \left(\beta \right)=\sum _{i=1}^n{\delta }_i\left\{{\mathbf{x}}_i^{\mathsf{T}}\beta -log\sum _{j:y_j\ge y_i}exp\left({\mathbf{x}}_j^{\mathsf{T}}\beta \right)\right\}.$$

(2)

When p is large and predictors are strongly correlated, regularisation is essential; more generally, the same framework covers regimes where p may grow with n.

2.2. Bayesian Elastic–Net Prior

We place an elastic–net prior on $\beta$ that unifies ${\ell }_1$ and ${\ell }_2$ shrinkage:

$$p(\beta \mid {\lambda }_1,{\lambda }_2)\propto exp\left(-{\lambda }_1{\parallel \beta \parallel }_1-\frac{{\lambda }_2}{2}{\parallel \beta \parallel }_2^2\right),$$

(3)

recovering lasso and ridge as ${\lambda }_2\to 0$ or ${\lambda }_1\to 0$ [2]. For computation we use the normal–exponential mixture for the Laplace part with a Gaussian ridge factor:

$$\begin{array}{cc}\hfill z_j\mid {\lambda }_1& \sim \mathrm{Exp}\left(\frac{{\lambda }_1^2}{2}\right),j=1,\dots ,p,\hfill \\ \hfill {\beta }_j\mid z_j,{\lambda }_2& \sim \mathcal{N}\left(0,{\left(\frac{1}{z_j}+{\lambda }_2\right)}^{-1}\right),\hfill \end{array}$$

(4)

which integrates back to (3). We use weakly-informative global hyperpriors,

$${\lambda }_1\sim \mathrm{Half}-\mathrm{Cauchy}(0,5),{\lambda }_2\sim \mathrm{Half}-\mathrm{Cauchy}(0,5),$$

(5)

so the ${\ell }_1/{\ell }_2$ balance is learned from the data. The ${\ell }_1$ part promotes sparsity; the ${\ell }_2$ part stabilises groups of correlated features (useful for gene modules).

2.3. Posterior Distribution

Let $z={(z_1,\dots ,z_p)}^{\mathsf{T}}$. Up to normalising constants, the joint posterior is

$$\pi (\beta ,z,{\lambda }_1,{\lambda }_2\mid \mathbf{y},\delta ,\mathbf{X})\propto exp\left(\ell \left(\beta \right)\right)\prod _{j=1}^p\varphi \left({\beta }_j;0,{\left(\frac{1}{z_j}+{\lambda }_2\right)}^{-1}\right){\displaystyle \frac{{\lambda }_1^2}{2}}{e}^{-{\displaystyle \frac{{\lambda }_1^2}{2}}{z}_j}\pi \left({\lambda }_1\right)\pi \left({\lambda }_2\right).$$

Because the Cox likelihood is non-conjugate, we use Hamiltonian Monte Carlo for full Bayesian inference.

2.4. Hamiltonian Monte Carlo Inference

We fit the model with the No-U-Turn Sampler [16,17]. The parameter block is

$$\mathsf{\Theta }=(\beta ,logz,log{\lambda }_1,log{\lambda }_2),$$

and the target potential $U\left(\mathsf{\Theta }\right)=-log\pi (\mathsf{\Theta }\mid \mathrm{data})$ uses (2) and (4). We run four chains, warm up to adapt step size and diagonal metric, then draw posterior samples. Convergence is assessed with $\widehat{R}<1.01$ and effective sample size thresholds; energy/BFMI diagnostics are checked for geometry pathologies. Posterior predictive survival uses Breslow’s estimator per draw. In this context, let $y_{\left(1\right)}<\dots <y_{\left(K\right)}$ denote the distinct ordered event times in the sample, let $d_k$ be the number of events at time $y_{\left(k\right)}$, and define the risk set $R\left(k\right)=\{j:y_j\ge y_{\left(k\right)}\}$. Then Breslow’s estimator at posterior draw m is given by

$${\widehat{H}}_0^{\left(m\right)}\left(t\right)=\sum _{y_{\left(k\right)}\le t}{\displaystyle \frac{d_k}{{\sum }_{j\in R_{\left(k\right)}}exp\left({\mathbf{x}}_j^{\mathsf{T}}{\beta }^{\left(m\right)}\right)}},{\widehat{S}}^{\left(m\right)}(t\mid \mathbf{x})=exp\left(-{\widehat{H}}_0^{\left(m\right)}\left(t\right)e^{{\mathbf{x}}^{\mathsf{T}}{\beta }^{\left(m\right)}}\right).$$

Before turning to the theoretical properties, we briefly summarise how the proposed procedure works. Figure 1 shows a DAG of the approach, highlighting the parameters used in the model and their probabilistic relationships.

In Figure 1, HC(0,5) is a Half-Cauchy hyperpriors with scale 5, providing weakly informative priors for the global shrinkage parameters ${\lambda }_1$ and ${\lambda }_2$, as specified in (5). ${\lambda }_1$ is a global ${\ell }_1$ penalty parameter controlling the degree of sparsity in the elastic–net prior; learned from the data via the Half-Cauchy(0,5) hyperprior and typically scaling as $\sqrt{(logp)/n}$ in high-dimensional regimes. ${\lambda }_2$ is again a global ${\ell }_2$ penalty parameter controlling ridge regularization strength; keeping ${\lambda }_2>0$ ensures posterior propriety and stabilizes HMC sampling. $z_j$ is latent exponential scale variables with $z_j\sim \mathrm{Exp}({\lambda }_1^2/2)$, enabling the normal–exponential mixture representation of the Laplace component of the elastic–net prior in (4). ${\beta }_j$ is regression coefficient for the j-th gene-expression feature (covariate ${\mathbf{x}}_{\mathbf{j}}$); conditionally Gaussian given $z_j$ and ${\lambda }_2$, with ${\beta }_j\mid z_j,{\lambda }_2\sim \mathcal{N}\left(0,{(1/z_j+{\lambda }_2)}^{-1}\right)$, as in (4). $h_0\left(t\right)$ is unspecified baseline hazard function in the Cox model (1); estimated post hoc via Breslow’s estimator at each posterior draw. ${\eta }_i$ is our linear predictor for subject i and ${\eta }_i={\mathbf{x}}_i^{\mathsf{T}}\beta$, representing the log hazard ratio in (1). ${\mathbf{x}}_i$: p-dimensional vector of z-scored ${log}_2$ gene-expression measurements for subject i, pre-processed as described in Section 3.4. Finally, $y_i$ is a observed survival time for subject i possibly right censored and ${\delta }_i$ is as known event indicator, with ${\delta }_i=1$ if death is observed and ${\delta }_i=0$ if the observation is right-censored at time $y_i$.

2.5. Theoretical Properties

This section summarizes the theoretical properties for the introduced BEN–Cox approach. In particular, Theorem 1 gives a Bayesian analogue of the elastic-net grouping idea at the posterior mode, and Theorem 2 establishes posterior contraction under sparsity. A complete proof of Theorem 2 including test construction, prior mass, and entropy bounds in the Cox partial-likelihood setting is provided in Appendix A. In this context we use the following notation: ${\parallel ·\parallel }_q$ is the ${\ell }_q$ norm; $I_p$ is the $p\times p$ identity matrix; $R_i=\{j:y_j\ge y_i\}$ is the Cox risk set. In addition, $I\left(\beta \right)=-{\nabla }^2\ell \left(\beta \right)$ are the observed partial-likelihood information. Posterior is $\mathrm{\Pi }(·\mid \mathrm{data})$ under (3)–(4). Also the following assumptions are needed to ensure the theoretical soundness.

Assumptions

(A1)

Columns of $\mathbf{X}$ are z-scored on the training split; $\parallel {\mathbf{x}}_{·j}{\parallel }_2^2/n=1$ and entries are bounded.

(A2)

As a Cox model regularity condition, event times lie in a compact interval and the baseline hazard is locally bounded; $\ell \left(\beta \right)$ is twice continuously differentiable (see [1]).

(A3)

There exists $\kappa >0$ such that for all v with support S, $\left|S\right|\le s_0$, $v^{\mathsf{T}}I\left(\beta \right)v\ge \kappa {\parallel v\parallel }_2^2$ in a neighbourhood of ${\beta }^{★}$.

(A4)

As a sparsity assumption, the truth ${\beta }^{★}$ is $s_0$–sparse with $s_{0}logp=o\left(n\right)$.

(A5)

${\lambda }_1,{\lambda }_2$ follow the half–Cauchy priors in (5) with positive density near $\{{\lambda }_1\asymp \sqrt{(logp)/n},0<{\lambda }_2\le C\}$ (see [18]).

We can explain the details of assumptions individually as follows: (A1) standardizes columns so correlation and curvature arguments are interpretable and avoids scale pathologies in HMC, which is crucial for Bayesian inference. When it comes to (A2), it is the usual smoothness condition that guarantees a well-behaved partial likelihood. (A3) is the high-dimensional analogue of identifiability: restricted curvature prevents flat directions on sparse supports, which can be interpreted as the Cox counterpart of restricted eigenvalues. (A4) encodes the high-dimensional regime with effective dimensionality $s_0$; it is the minimal condition for $\sqrt{(s_{0}logp)/n}$ rates and finally (A5) ensures the prior places enough mass near the optimal ${\lambda }_1$ scale while keeping ${\lambda }_2>0$ so the geometry is strongly convex and the posterior proper.

In the following two basic geometric properties are introduced that we will use repeatedly: (i) the concavity of the Cox partial likelihood, which controls the curvature of the log-posterior, and (ii) the mixture representation of the elastic–net prior, which links our hierarchical formulation (4) back to the penalty in (3). These results are standard, but they are provided for completeness.

Lemma 1

(Concavity of the Cox partial likelihood.). Under (A2), $\ell \left(\beta \right)$ in (2) is concave in β and $I\left(\beta \right)=-{\nabla }^2\ell \left(\beta \right)$ is positive semi-definite (see Cox [1] for standard log-sum-exp concavity).

Lemma 2.

The latent parameterisation (4) integrates to the elastic–net prior (3):

$${\int }_0^{\infty }\varphi \left({\beta }_j;0,{\left(\frac{1}{z_j}+{\lambda }_2\right)}^{-1}\right)·{\displaystyle \frac{{\lambda }_1^2}{2}}{e}^{-({\lambda }_1^2/2)z_j}d{z}_j\propto exp\left(-{\lambda }_1\left|{\beta }_j\right|-\frac{{\lambda }_2}{2}{\beta }_j^2\right).$$

See Li and Lin [11] for normal–exponential representation of Laplace.

With ${\lambda }_2>0$ the Gaussian ridge factor controls the tails in all directions. Since $\ell \left(\beta \right)$ is concave (see Lemma 1) and the hyperpriors are proper according to assumption A5, the posterior normalizes. Moreover, $-\ell \left(\beta \right)+{\displaystyle \frac{{\lambda }_2}{2}}{\parallel \beta \parallel }_2^2$ is strongly convex with curvature $\ge {\lambda }_2$, so the posterior mode in $\beta$ is unique. Therefore, this stabilizes the optimization and HMC geometry.

With these building blocks in place, we now formalise two behaviours that are particularly relevant in our high-dimensional gene-expression setting. First, the grouping effect shows that highly correlated gene features tend to receive similar coefficients, so the model naturally picks gene modules rather than arbitrary single gene features. Second, a contraction result shows that, under sparsity, the BEN–Cox posterior concentrates around the true coefficient vector at the usual high-dimensional rate, providing reassurance that our calibrated risk scores are statistically well behaved as n grows. We state a grouping theorem for the posterior mode, followed by a contraction theorem for the full posterior.

Theorem 1

(Bayesian elastic–net grouping at the posterior mode). Assume (A1)–(A2). Let ${\widehat{\beta }}_{\mathrm{MAP}}$ be the posterior mode under (3) with fixed $({\lambda }_1,{\lambda }_2)$ and ${\lambda }_2>0$. For standardized columns $j,k$ with sample correlation ${\rho }_{jk}$,

$$\left|{\widehat{\beta }}_{\mathrm{MAP},j}-{\widehat{\beta }}_{\mathrm{MAP},k}\right|\le {\displaystyle \frac{\parallel \nabla \ell \left({\widehat{\beta }}_{\mathrm{MAP}}\right){\parallel }_{\infty }+2{\lambda }_1}{{\lambda }_2}}\sqrt{2(1-{\rho }_{jk})}.$$

In particular, if ${\rho }_{jk}\to 1$ and $(\parallel \nabla \ell \left({\widehat{\beta }}_{\mathrm{MAP}}\right){\parallel }_{\infty }+2{\lambda }_1)/{\lambda }_2$ remains bounded, then ${\widehat{\beta }}_{\mathrm{MAP},j}-{\widehat{\beta }}_{\mathrm{MAP},k}\to 0$. Hence the difference tends to zero (Bayesian analogue of the ElasticNet’s grouping effect; see Zou and Hastie [2], Hans [12]).

Proof. 

Consider $Q\left(\beta \right)=-\ell \left(\beta \right)+{\lambda }_1{\parallel \beta \parallel }_1+{\displaystyle \frac{{\lambda }_2}{2}}{\parallel \beta \parallel }_2^2$. Karush–Kuhn–Tucker (KKT) conditions give $\mathbf{0}=-\nabla \ell \left({\widehat{\beta }}_{MAP}\right)+{\lambda }_{1}s+{\lambda }_2{\widehat{\beta }}_{MAP}$ with $s_j\in [-1,1]$. Taking the j–k difference,

$${\lambda }_2({\widehat{\beta }}_{MAP,j}-{\widehat{\beta }}_{MAP,k})={[\nabla \ell \left({\widehat{\beta }}_{MAP}\right)]}_j-{[\nabla \ell \left({\widehat{\beta }}_{MAP}\right)]}_k-{\lambda }_1(s_j-s_k),$$

so $|{\widehat{\beta }}_{MAP,j}-{\widehat{\beta }}_{MAP,k}|\le {\displaystyle \frac{{|[\nabla \ell ]}_j-{[\nabla \ell ]}_k|}{{\lambda }_2}}+{\displaystyle \frac{2{\lambda }_1}{{\lambda }_2}}.$ For the score difference, express the Cox score using risk-set weights and apply Cauchy–Schwarz to the difference column $\mathbf{d}={\mathbf{x}}_{·j}-{\mathbf{x}}_{·k}$. With standardisation (A1), ${\parallel \mathbf{d}\parallel }_2^2=2n(1-{\rho }_{jk})$, yielding ${|[\nabla \ell ]}_j-{[\nabla \ell ]}_k|\le \left({\sum }_i{\delta }_i\right)\sqrt{2n(1-{\rho }_{jk})}.$ Scaling per observation and absorbing ${\sum }_i{\delta }_i/n=O\left(1\right)$ gives the stated bound with $\parallel \nabla \ell \left({\widehat{\beta }}_{MAP}\right){\parallel }_{\infty }/{\lambda }_2·\sqrt{2(1-{\rho }_{jk})}$. If ${\rho }_{jk}\to 1$ and ${\lambda }_1/{\lambda }_2$ is not diverging, the difference vanishes. Since the posterior typically concentrates in a neighborhood of ${\widehat{\beta }}_{\mathrm{MAP}}$ under standard regularity, the grouping behavior at the mode is often reflected in posterior summaries as well; we observe this empirically in the METABRIC analysis.    □

Since the posterior distribution is (locally) log-concave and concentrates around ${\widehat{\beta }}_{MAP}$, the grouping effect at the mode also appears in posterior means. In other words, highly correlated gene features tend to have similar posterior mean coefficients as well; see Li and Lin [11], Hans [12] for related results on Bayesian elastic–net concentration.

Theorem 2.

(Posterior contraction). Assume (A1)–(A5) and let ${ϵ}_n\asymp \sqrt{(s_{0}logp)/n}$. Then there exists $M>0$ such that

$$\mathrm{\Pi }\left(\parallel \beta -{\beta }^{★}{\parallel }_2<M{ϵ}_n|\mathbf{y},\delta ,\mathbf{X}\right)\stackrel{P}{\to }1.$$

This theorem is a posterior contraction result: as n grows, the BEN–Cox posterior for $\beta$ concentrates inside an ${\ell }_2$–ball of radius ${ϵ}_n\asymp \sqrt{(s_{0}logp)/n}$ around the sparse truth ${\beta }^{★}$. The proof follows the standard Bayesian contraction arguments of Ghosal and van der Vaart [13], adapted to the Cox partial-likelihood setting, and combines three ingredients. (i) Tests: Assumption (A3) provides local quadratic lower bounds for $\ell \left(\beta \right)$, which yield exponential tests that separate ${\beta }^{★}$ from the complement of an ${\ell }_2$–ball of radius $M{ϵ}_n$. (ii) Prior mass: By Lemma 2 and (A5), the prior allocates at least $exp(-C{s}_{0}logp)$ mass to such balls (taking ${\lambda }_1\asymp \sqrt{(logp)/n}$ and using ${\lambda }_2>0$ to stabilize curvature). (iii) Complexity: The class of $s_0$–sparse vectors within radius $M{ϵ}_n$ has logarithmic complexity of order $s_{0}logp$, which is dominated by $n{ϵ}_n^2$. A complete, proof specialised to our Cox partial-likelihood setting (including the testing argument and the hyperprior treatment) is given in Appendix A.

Under Assumptions (A1)–(A4) and the contraction result in Theorem 2, the posterior for $(\beta ,H_0)$ concentrates around the true pair $({\beta }^{★},H_0^{★})$. The survival function $S(t\mid \mathbf{x};\beta ,H_0)$ is continuous in $(\beta ,H_0)$ on compact time intervals, and Breslow’s estimator is uniformly consistent for the baseline cumulative hazard in the Cox model [19]. Combining these facts, one can obtain:

$$\underset{t\in [0,\tau ]}{sup}\left|\mathbb{E}\left[{\widehat{S}}^{\left(m\right)}(t\mid \mathbf{x})\mid \mathrm{data}\right]-S_0(t\mid \mathbf{x})\right|\stackrel{P}{\to }0,$$

for any fixed covariate vector $\mathbf{x}$ and finite horizon $\tau$. In other words, the posterior predictive survival curves converge to the true survival function uniformly on compact time intervals.

From a computational point of view, the ridge component with ${\lambda }_2>0$ induces strong convexity in $\beta$, while the z-scaling in (A1) keeps the design entries bounded. [16,17].

For practical implementation, it is helpful to summarize the full BEN–Cox workflow from raw gene-expression data to posterior summaries and predictive metrics. The following step-by-step algorithm makes explicit how the model, prior, MCMC inference, and evaluation pieces fit together in our analysis.

3. Data

3.1. Cohort Description

We analyse the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) cohort [14,20]. The study combines tumours from two centres (Cambridge, UK and Vancouver, Canada) and is publicly available through the cBioPortal interface [21]. The dataset contains clinical variables, overall survival, and a gene–expression matrix. From this file we extract overall survival and gene–expression information and remove subjects with missing survival time or status. After this quality control step we retain $n=1903$ primary invasive breast cancers for our analysis.

3.2. Outcome

The time–to–event end-point is overall survival (OS). In the cBioPortal file OS is reported in months; we convert this to days for the modelling step and define the event indicator ${\delta }_i=1$ for deaths and ${\delta }_i=0$ for right-censored observations. Among the 1903 subjects, 800 experience the event (42.0 %), and the median follow-up is approximately 16.4 years, estimated using the reverse Kaplan–Meier method.

3.3. Predictors

• Gene expression : The METABRIC dataset provides Illumina HT-12 v3 micro-array measurements summarised at the gene level. In the version used here this block contains $p_{\mathrm{raw}}=489$ gene–expression features, indexed from the BRCA1 gene onwards. While the original METABRIC platform measures many more micro-array gene-level features, we restrict attention to this gene-level expression matrix for a more stable and reproducible design.
• Clinical covariates : age at diagnosis, tumour size, histological grade, lymph-node status and treatment indicators (hormone/chemotherapy) are also available in the downloaded dataset. These variables are not used in the present study but are retained for potential future work and for possible extensions of the BEN–Cox model.

3.4. Preprocessing

We first remove genes whose variance lies below the 10th percentile across patients, leaving $p=440$ gene–expression features. This variance filter removes almost-constant gene features and improves the conditioning of the design matrix. The data are then split once into an 80 % training set and a 20 % test set, stratified by the event indicator so that the censoring fraction is similar in both splits. Each gene is z-scored in the training set to have mean 0 and variance 1, and the same centring and scaling are applied to the test set. Finally, the survival tuples $(y_i,{\delta }_i)$ are encoded as described in Section 2.1, with $y_i$ measured in days and ${\delta }_i$ the event indicator.

4. Application

4.1. Experimental Design

We treat METABRIC as a single large cohort and mimic a typical prognostic modelling workflow: a model is fitted once on a development set and then evaluated on a separate test set. The data are therefore split once into 80 % training and 20 % testing, stratified on the event indicator so that the proportion of deaths and censored observations remains comparable across the two splits.

All tuning and model fitting are performed using only the training data. For the ridge-penalised Cox baseline we use five-fold cross-validation to select the penalty parameter ${\lambda }_{\mathrm{ridge}}$, with the concordance index as the optimisation criterion. In contrast, the Bayesian elastic-net Cox (BEN–Cox) model treats the elastic–net penalties as random: the global shrinkage parameters $({\lambda }_1,{\lambda }_2)$ follow Half-Cauchy hyperpriors and are learned jointly with the regression coefficients. This avoids an additional tuning loop and lets the data determine the balance between ${\ell }_1$ and ${\ell }_2$ shrinkage.

Posterior inference for BEN–Cox relies on Hamiltonian Monte Carlo as implemented in Stan, using the No-U-Turn Sampler (NUTS) with adaptively tuned step size and diagonal mass matrix. Starting from zero coefficients and diffuse initial values for $({\lambda }_1,{\lambda }_2)$, the chain is run until trace plots stabilise and effective sample sizes for the regression coefficients are comfortably above standard thresholds. The initial part of the chain is discarded as warm-up, and all posterior summaries are based on the retained draws. In particular, we monitor mixing and autocorrelation for representative coefficients and for the global shrinkage parameters to ensure that the posterior exploration is adequate for downstream prediction and uncertainty quantification.

4.2. Baseline Models

To show the BEN–Cox model’s advantages and disadvantages and how well it performs, we compare it against two baseline Cox models:

• Null Cox: a Cox model with no gene-expression covariates, i.e., all subjects share the same baseline hazard $h_0\left(t\right)$ (no covariate effects). This provides a no-information baseline for discrimination and absolute-risk calibration.
• Ridge Cox : an ${\ell }_2$-penalised Cox model fitted with glmnet, where the penalty ${\lambda }_{\mathrm{ridge}}$ is chosen by five-fold cross-validation on the training set. This is a natural frequentist benchmark for high-dimensional survival modelling when sparsity is not enforced explicitly.

Also note that, all models are fitted on exactly the same training data and are evaluated on the same held-out test set.

4.3. Evaluation Metrics

We assess predictive performance along three ways: time-dependent prediction error, discrimination between subjects with different risk, and calibration of predicted survival probabilities.

Let $T_i$ denote the true event time for patient i, $C_i$ the censoring time, and $Y_i=min(T_i,C_i)$ the observed time, with ${\delta }_i=\mathbf{1}\{T_i\le C_i\}$ the event indicator. For a given model, let ${\widehat{S}}_i\left(t\right)$ be the predicted survival probability at time t for patient i, and let ${\eta }_i$ be its linear predictor or risk score.

As mentioned in the Algorithm 1, we have 3 main performance metrics that are IBS, C-index and calibration slope and metrics. These metrics are defined as follows:

Algorithm 1: BEN–Cox pipeline

Input: gene matrix $\mathbf{X}$, times $\mathbf{y}$, events $\delta$.

Output: posterior draws ${\left\{{\beta }^{\left(m\right)}\right\}}_{m=1}^M$, test-set performance metrics.

1. Pre-processing

1.1 Remove gene-level features whose variance is below the 10th percentile [14].

1.2 z-score each column of $\mathbf{X}$ using training-set moments.

2. Train–test split

2.1 Randomly split subjects 80/20 into training and test sets, stratified by $\delta$.

2.2 Fix the random seed and record the split for reproducibility.

3. Model and prior

3.1 Specify the Cox model (1) with elastic–net prior (3).

3.2 Use the normal–exponential mixture representation (4) and

Half-Cauchy hyperpriors for $({\lambda }_1,{\lambda }_2)$ as in (5).

4. HMC sampling with Stan (training set only)

4.1 Fit the model using 4 chains, warm-up + sampling.

4.2 Check convergence: $\widehat{R}<1.01$ and ESS$>400$ for all coefficients.

4.3 Retain posterior draws $\left\{{\beta }^{\left(m\right)}\right\}$ and, for each draw, compute

the Breslow baseline cumulative hazard ${\widehat{H}}_0^{\left(m\right)}\left(t\right)$.

5. Prediction on the test set

5.1 For each test subject i and draw m, compute

linear predictor ${\eta }_i^{\left(m\right)}={\mathbf{x}}_i^{\mathsf{T}}{\beta }^{\left(m\right)}$ and

survival curve ${\widehat{S}}^{\left(m\right)}(t\mid {\mathbf{x}}_i)=exp\left(-{\widehat{H}}_0^{\left(m\right)}\left(t\right)e^{{\eta }_i^{\left(m\right)}}\right)$.

5.2 Form posterior summaries (e.g., mean or median survival probabilities).

6. Performance metrics (test set)

6.1 Compute the integrated Brier score (IBS).

6.2 Compute the concordance index (C-index) for discrimination.

6.3 Compute calibration slope and intercept, following Demler et al. [15].

Integrated Brier Score:

The Brier score at time t is a squared prediction error adapted to the survival setting. Following the inverse-probability-of-censoring weighting (IPCW) approach, we write

$$BS\left(t\right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left[{\displaystyle \frac{\mathbf{1}\{Y_i\le t,{\delta }_i=1\}}{\widehat{G}(Y_i-)}}{\left(0-{\widehat{S}}_i\left(t\right)\right)}^2+{\displaystyle \frac{\mathbf{1}\{Y_i>t\}}{\widehat{G}\left(t\right)}}{\left(1-{\widehat{S}}_i\left(t\right)\right)}^2\right],$$

where $\widehat{G}$ is an estimate of the censoring survival function (e.g. Kaplan–Meier for $C_i$). Intuitively, patients who are still under follow-up at time t contribute squared error for predicting survival, and patients who die before t contribute squared error for predicting death, with both contributions reweighted to account for censoring.

To summarize performance over a time window $[0,\tau ]$ we use the integrated Brier score (IBS),

$$IBS\left(\tau \right)={\displaystyle \frac{1}{\tau }}{\int }_0^{\tau }BS\left(t\right)dt,$$

which can be interpreted as a time-averaged mean squared prediction error for survival probabilities up to horizon $\tau$. In practice we approximate the integral numerically on a fine grid of time points.

Concordance (C) Index:

Discrimination is quantified using Harrell’s concordance index (C-index). The C-index estimates the probability that, for a randomly chosen pair of comparable patients $(i,j)$, the one who dies first also has the higher predicted risk. Formally,

$$\widehat{C}={\displaystyle \frac{{\sum }_{i<j}\mathbf{1}\{Y_i<Y_j,{\delta }_i=1\}\mathbf{1}\{{\eta }_i>{\eta }_j\}+\mathbf{1}\{Y_j<Y_i,{\delta }_j=1\}\mathbf{1}\{{\eta }_j>{\eta }_i\}}{{\sum }_{i<j}\mathbf{1}\{Y_i<Y_j,{\delta }_i=1\}+\mathbf{1}\{Y_j<Y_i,{\delta }_j=1\}}},$$

where the denominator counts all comparable pairs (those with an observed event before the other’s follow-up has finished) and the numerator counts how often the model orders their risk correctly.

Grouped Calibration (GND ${\chi }^2$):

In addition to the calibration slope and intercept, we also report a Greenwood–Nam–D’Agostino (GND) type grouped-calibration chi-square type statistic [15] as a global calibration measure. To achieve that we consider, $\tau$ which means a fixed, clinically relevant follow-up time (in this study, horizon $\tau =5$ years), and we only look at whether each patient has experienced the event before this time or not.

For each patient we define a horizon event indicator

$$Z_i\left(\tau \right)=\mathbf{1}\{T_i\le \tau ,{\delta }_i=1\},$$

which is 1 if the event occurs on or before time $\tau$ and 0 otherwise. We then group patients into K risk groups (for example, deciles, $K=10$) according to their predicted event probabilities ${\widehat{\pi }}_i\left(\tau \right)$. Let $G_k$ denote the index set of patients in group k, $n_k=\left|G_k\right|$ be the group size, and

$$o_k=\sum _{i\in G_k}{Z}_i\left(\tau \right),\overline{Z}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{Z}_i\left(\tau \right),e_k=n_k\overline{Z}$$

be the observed and expected numbers of events in group k. The GND statistic is then given by

$$\mathrm{GND}=\sum _{k=1}^K{\displaystyle \frac{{\left(o_k-e_k\right)}^2}{e_k}},$$

which summarize observed and expected event counts across the risk groups; smaller GND values indicate better global calibration of absolute risk at horizon $\tau$. Please note that this is a simplified Greenwood–Nam–D’Agostino–type statistic, using the overall event rate as the reference expected rate.

Calibration:

Calibration asks whether predicted risks match observed event frequencies. We focus on a clinically relevant time horizon $\tau$ (e.g. 5 years) and consider the predicted event probability ${\widehat{\pi }}_i\left(\tau \right)=1-{\widehat{S}}_i\left(\tau \right)$ for each patient. Following Demler et al. [15], we fit a recalibration model of the form

$$logit\left\{Pr\left(T_i\le \tau \mid {\widehat{\pi }}_i\left(\tau \right)\right)\right\}={\alpha }_0+{\alpha }_{1}logit\left({\widehat{\pi }}_i\left(\tau \right)\right),$$

with suitable weighting to account for censoring at $\tau$. A calibration slope ${\alpha }_1$ close to one and intercept ${\alpha }_0$ close to zero indicate that the model’s predicted risks agree well with the observed data; values ${\alpha }_1<1$, for example, suggest that the model overstates differences between low- and high-risk patients.

Bootstrap Uncertainty:

To quantify uncertainty in these performance estimates we use paired bootstrap resampling of the test set. In each of B bootstrap replicates we resample test patients with replacement, recompute the IBS, C-index, and calibration parameters $({\alpha }_0,{\alpha }_1)$ for each model, and then summarise the resulting distribution by its mean and standard error. This provides a simple way to assess the sampling variability of the performance metrics and to compare BEN–Cox with the ridge and null baselines.

4.4. Results

This section reports the METABRIC data analysis results for the Bayesian elastic–net Cox (BEN–Cox) model and the two baseline Cox models on the held-out test set. A brief summary of the cleaned dataset is given in Table 1 to keep the analysis context clear.

The final BEN–Cox model is fitted with Stan using Hamiltonian Monte Carlo (HMC). The global shrinkage parameters concentrate around ${\lambda }_1\approx 3.08$ (L1 component) and ${\lambda }_2\approx 0.84$ (L2 component), indicating a non-trivial combination of sparsity-inducing and ridge-type shrinkage. Convergence diagnostics for the fitted model are provided in Table 2. The split-$\widehat{R}$ values are essentially 1.00 for the median regression coefficient and for both global penalty parameters, and the effective sample sizes (ESS) are in the low thousands, which indicates stable posterior exploration. Therefore, it can be said that the model parameters have converged well. Finally, although we estimate the recalibration intercept and slope $({\alpha }_0,{\alpha }_1)$ for each model using the scheme in Section 4.3, we do not tabulate them here; instead, we summarize calibration through the GND ${\chi }^2$ statistic in Table 3 and the calibration plot in Figure 4.

Trace plots for the penalty parameters ${\lambda }_1$ and ${\lambda }_2$ and randomly selected regression coefficients are given in Figure 2. In these trace plots, the chains mix well with no visible trends or drifts, which is consistent with the numerical diagnostics. Based on posterior credible intervals, the BEN–Cox model retains 48 genes with 95% credible intervals that exclude zero, yielding a relatively sparse signature compared with the 440 pre-filtered gene-level features.

Regarding the predictive performance of the introduced BEN-Cox model, Table 3 is given below which summarizes predictive performance on the 20% held-out test set for three models: a null Cox model (no covariates), a ridge-penalized Cox model, and the BEN–Cox model. Performance is evaluated using the integrated Brier score (IBS), Harrell’s C-index, and a grouped-calibration which is simplified version of Greenwood–Nam–D’Agostino (GND) goodness-of-fit statistic (GND ${\chi }^2$). All values are reported as mean ± standard error over 100 bootstrap resamples of the test cohort.

As expected, the null model has no discrimination (C-index 0.50) and serves as a calibration baseline. Both penalized models clearly improve discrimination: the ridge Cox model reaches a C-index of 0.647, and BEN–Cox achieves a slightly higher value of 0.655. The integrated Brier scores follow the same pattern: BEN–Cox gives a smaller IBS (0.216) than ridge (0.224), with the difference on the order of one standard error. Overall, both penalized models extract meaningful prognostic signal from the gene-expression data, but the fully Bayesian BEN–Cox model is marginally more accurate and more discriminative on the held-out test set with much less number of covariates.

The GND statistic points in the same direction. BEN–Cox shows a moderate lack of fit (GND $\approx 19$), whereas the ridge Cox model exhibits a substantially larger value (around 89 with considerable variability), indicating more severe global miscalibration of absolute event probabilities. The null model, which predicts a constant hazard for all patients, naturally achieves a small GND but at the cost of zero discrimination. In this sense, BEN–Cox achieves a more favourable balance between discrimination and calibration than the ridge baseline.

In order to visualise how the BEN–Cox risk score stratifies patients, we divided the test set into quintiles of the linear predictor and plotted Kaplan–Meier curves for each group (Figure 3). It can be said that there is a clear gradient across quintiles: higher-risk groups show consistently worse survival, and the lowest-risk quintile has visibly better outcomes. This pattern is consistent with the C-index values and illustrates that BEN–Cox provides clinically meaningful risk stratification on the held-out data.

Calibration of predicted event probabilities is examined in Figure 4, which shows observed versus predicted 5-year event probabilities for BEN–Cox. The overall trend follows the identity line reasonably well, with some visible deviations, especially towards the extremes of the predicted risk distribution. These departures are in line with the moderate GND statistic and suggest that, while BEN–Cox is not perfectly calibrated, its absolute risk estimates are substantially better behaved than those of the ridge model and adequate for exploratory prognostic modelling.

The posterior forest plot for the 30 largest BEN–Cox coefficients (Figure 5) highlights a sparse set of genes with clearly non-zero effects and reasonably tight credible intervals. Among the 48 genes selected by BEN–Cox, we recover ERBB2, a key component of the PAM50 breast-cancer panel, which supports the biological plausibility of the learned signature with less covariates.

Figure 4. Calibration of BEN–Cox 5-year event probabilities on the test set. Points show observed event proportions across groups of patients with similar predicted risk; the dashed line is the ideal ${45}^{\circ }$ line. Deviations from the line reflect residual miscalibration, consistent with the moderate GND statistic.

Figure 4. Calibration of BEN–Cox 5-year event probabilities on the test set. Points show observed event proportions across groups of patients with similar predicted risk; the dashed line is the ideal ${45}^{\circ }$ line. Deviations from the line reflect residual miscalibration, consistent with the moderate GND statistic.

Figure 5. Posterior means and 95% credible intervals for the 30 largest BEN–Cox regression coefficients. The selected panel is sparse (48 genes overall) and includes ERBB2 from the PAM50 signature.

Figure 5. Posterior means and 95% credible intervals for the 30 largest BEN–Cox regression coefficients. The selected panel is sparse (48 genes overall) and includes ERBB2 from the PAM50 signature.

Table 4 summarizes coverage of the PAM50 gene panel for the ridge and BEN–Cox models. The ridge Cox model, which retains all 440 pre-filtered features with non-zero coefficients, effectively achieves full recall of the PAM50 genes present in the expression matrix. In contrast, BEN–Cox retains 48 genes in total and recovers about one third of the PAM50 genes available in this reduced panel (33.3% recall). This behaviour is consistent with the design of the elastic–net prior: it aggressively shrinks weak or redundant signals, aiming to identify a compact subset of strongly prognostic genes rather than reproduce the full PAM50 signature.

In summary, the BEN–Cox model on METABRIC satisfies standard convergence diagnostics, yields a sparse and interpretable panel of 48 prognostic genes, and achieves *slightly better* discrimination and prediction error on the held-out test set than a tuned ridge-penalised Cox model. Both penalised models substantially outperform the null Cox model in terms of discrimination, but BEN–Cox attains a smaller IBS and a higher C-index than ridge, with differences that are small but consistent across bootstrap resamples.

From a calibration perspective, the GND statistic indicates that BEN–Cox has a moderate lack of fit, whereas the ridge model shows a much larger global miscalibration signal. Taken together, these results suggest that BEN–Cox offers a more favourable trade-off: it delivers marginally better discrimination and lower prediction error than ridge while also exhibiting substantially better calibration and a much more compact gene signature. The fact that the selected panel still recovers ERBB2 from the PAM50 panel aligns with known biology and supports the interpretability of the learned effects.

These findings support the idea that a fully Bayesian elastic–net prior can deliver a sparse, biologically plausible high-dimensional survival model that slightly improves upon ridge in discrimination, prediction error, and calibration, while offering full posterior uncertainty quantification for hazard ratios and survival curves.

5. Discussion

The Bayesian elastic–net Cox model achieved discrimination on METABRIC that is slightly better than a tuned ridge-penalized Cox baseline, while producing a substantially sparser set of prognostic genes and providing coherent posterior uncertainty estimates. At the same time, the BEN–Cox model showed more favourable calibration according to the GND statistic, suggesting that a carefully constructed global–local prior can improve both discrimination and global fit relative to a purely ${\ell }_2$-penalized approach. The grouping result in Theorem 1 helps explain the consistent retention of correlated gene-level feature clusters, in line with known gene modules and signatures such as ERBB2-related pathways.

Beyond the empirical METABRIC analysis, the methodological contribution is the theoretical support for BEN–Cox in the high-dimensional Cox setting. Especially, Theorem 1 formalizes a Bayesian analogue of the elastic-net grouping effect at the posterior mode, which is practically important in gene-expression studies because strongly correlated genes in the dataset tend to move together and form biologically meaningful modules. Moreover, Theorem 2 establishes posterior contraction under sparsity for the Cox partial-likelihood geometry (with a ridge component ensuring stability), which provides reassurance that the inferred risk score and the associated uncertainty quantification are statistically well-behaved as n grows. In this sense, the paper is not only an applied Bayesian survival analysis, but also a theoretical shrinkage construction specialized to the Cox partial-likelihood framework.

Limitations of the current study include computational cost (HMC is more expensive than deterministic optimisation), and the fact that all results are derived from a single cohort. External validation on independent datasets, and extensions to multi-cohort or hierarchical BEN–Cox formulations, would provide a stronger assessment of transportability. Another limitation is that we focused solely on gene-expression features; integrating clinical covariates and mutation data may further improve both calibration and interpretability.

Regarding the assumption justification of the introduced estimator in this specific METABRIC data analysis, by considering (A1–A5), we can infer the following: Our pre-processing enforces (A1) by z-scoring each feature on the training split and carrying scalings to the test set, which also bounds design entries after low-variance filtering. Regarding (A2), overall-survival times are observed on a finite horizon and the Cox partial likelihood is smooth; Breslow’s baseline estimator is standard in this setting [19]. For (A3), the compatibility/restricted-eigenvalue condition is plausible after variance filtering and with the ${\ell }_2$ stabiliser in the prior; empirically, the observed partial-likelihood information has well-behaved restricted spectra on supports of the selected size. The sparsity/growth condition (A4) is consistent with our use of elastic–net shrinkage to target a low-to-moderate number of effective gene features relative to n. Finally, (A5) holds by construction: half–Cauchy hyperpriors put sufficient mass near the canonical ${\lambda }_1\asymp \sqrt{(logp)/n}$ scale while keeping ${\lambda }_2>0$ with high probability [18]. Overall, these checks support the applicability of the theoretical guarantees to the METABRIC data analysis.

Future work will explore multi-study hierarchical BEN models, integration of additional omics layers (e.g. somatic mutation data), and explicit calibration-improving extensions, for example, by combining BEN–Cox with Bayesian isotonic regression or flexible baseline hazard modelling.

6. Conclusion

In this study, we introduced a fully Bayesian elastic–net Cox (BEN–Cox) model for high-dimensional time–to–event analysis and applied it to the METABRIC breast cancer cohort using only gene–expression features. The model combines an elastic–net type prior with HMC to produce sparse coefficient estimates together with full posterior uncertainty for hazard ratios and survival curves. A fully Bayesian elastic–net Cox model offers an interpretable and sparse risk score using only gene–expression features from the METABRIC dataset. On this moderate level high–dimensional dataset, with 440 gene–level predictors after quality control, the BEN–Cox model achieves slightly lower integrated Brier score, slightly higher C–index, and a clearly smaller GND ${\chi }^2$ statistic than a tuned ridge–penalised Cox model on the held–out test set (Table 3). In other words, within this setting BEN–Cox provides a modest but consistent improvement in prediction error, discrimination, and global calibration while still yielding a compact number of genes and full posterior uncertainty quantification.

The main conclusions of this study can be summarized as follows:

• BEN–Cox selects a small set of prognostic genes (48 in our analysis), with posterior credible intervals that help quantify the strength and direction of each effect, and it recovers biologically meaningful markers such as ERBB2.
• In terms of predictive performance, compared with a tuned ridge Cox model, BEN–Cox shows slightly better integrated Brier score, slightly higher C–index, and noticeably better global calibration as measured by the GND statistic, while both models clearly outperform the null Cox baseline.
• The Stan HMC implementation delivers full posterior distributions for regression coefficients and survival curves, offering coherent uncertainty quantification at both the gene and patient level.
• Despite the Bayesian formulation and MCMC, the complete pipeline (fitting, bootstrap evaluation, and plotting) remains computationally feasible on standard hardware for a cohort of this size.
• Here it can be said that the introduced theoretical results emphasize the theoretical contribution of the paper which is the Bayesian grouping property explains stable behavior under strong gene correlations, and the posterior contraction result provides an asymptotic guarantee that the BEN–Cox posterior concentrates around sparse truth at the usual high-dimensional rate in the Cox partial-likelihood setting.

At the same time, the remaining imperfections in calibration and the reliance on a single cohort highlight the need for future work on Bayesian calibration strategies, multi–cohort or hierarchical extensions of BEN–Cox, and integration of additional clinical and molecular covariates in high–dimensional time–to–event prediction.