Bayesian variable and hazard structure selection in the General Hazard model

The proportional hazards (PH) and accelerated failure time (AFT) models are the most widely used hazard structures for analysing time-to-event data. When the goal is to identify variables associated with event times, variable selection is typically performed within a single hazard structure, imposing strong assumptions on how covariates affect the hazard function. To allow simultaneous selection of relevant variables and the hazard structure itself, we develop a Bayesian variable selection approach within the general hazard (GH) model, which includes the PH, AFT, and other structures as special cases. We propose two types of g-priors for the regression coefficients that enable tractable computation and show that both lead to consistent model selection. We also introduce a hierarchical prior on the model space that accounts for multiplicity and penalises model complexity. To efficiently explore the GH model space, we extend the Add-Delete-Swap algorithm to jointly sample variable inclusion indicators and hazard structures. Simulation studies show accurate recovery of both the true hazard structure and active variables across different sample sizes and censoring levels. Two real-data applications are presented to illustrate the use of the proposed methodology and to compare it with existing variable selection methods.

💡 Research Summary

This paper introduces a Bayesian framework for simultaneous variable selection and hazard‑structure selection within the General Hazard (GH) model, which subsumes the proportional hazards (PH), accelerated failure time (AFT), and accelerated hazards (AH) models as special cases. The authors begin by formalising the GH model: the hazard function depends on two sets of covariates, one acting on the hazard scale (β) and another on the time scale (α). By re‑parameterising the baseline log‑location‑scale distribution, the likelihood becomes analytically tractable, allowing the inclusion of right‑censored observations without numerical integration.

To encode the dual role of each covariate, a categorical inclusion indicator γj∈{0,…,4} is defined, distinguishing “no effect”, “time‑scale effect only”, “hazard‑scale effect only”, and two versions of “both effects” (different or identical magnitudes). This yields a model space of size 4p + 2p − 1, where p is the number of candidate covariates.

The methodological core consists of two novel g‑prior constructions for the model‑specific coefficients (θ, η). The first, a likelihood‑curvature‑matching prior, replaces the intractable expected Fisher information with the sample average of the observed Fisher information and evaluates the covariance at the model‑specific maximum‑likelihood estimates. The second, a product‑g‑prior, treats the time‑scale and hazard‑scale coefficients independently, greatly simplifying computation. Both priors involve a single hyper‑parameter g, calibrated via the unit‑information principle, and are shown to lead to consistent Bayes factors.

A hierarchical prior on the model space extends the Beta‑Binomial multiplicity correction to accommodate covariates that may appear in multiple roles. This prior penalises model complexity while automatically adjusting for the enlarged combinatorial space, thereby mitigating over‑fitting.

For posterior computation, the authors adapt the Add‑Delete‑Swap (ADS) algorithm to jointly update inclusion indicators and hazard‑structure assignments. Three move types are employed: (i) addition or deletion of a covariate, (ii) swapping a covariate’s role between time and hazard scales, and (iii) switching the overall hazard structure (e.g., PH ↔ GH). Acceptance probabilities are derived from Metropolis‑Hastings ratios that incorporate marginal likelihood approximations obtained via Laplace and integrated‑Laplace methods. These approximations are fast and accurate enough to drive efficient exploration of the high‑dimensional model space.

Theoretical results establish model‑selection consistency: under mild regularity conditions, the posterior probability concentrates on the true γ vector and the correct hazard structure as the sample size n→∞, provided the g‑prior’s scaling satisfies the usual asymptotic requirements. The hierarchical model prior plays a crucial role in achieving this consistency for both variable and structure selection.

Simulation studies explore a range of scenarios (sample sizes 100–500, censoring rates 10%–50%, and true models drawn from PH, AFT, AH, and GH families). The proposed method consistently outperforms existing Bayesian variable‑selection approaches that are restricted to a single hazard structure (e.g., spike‑and‑slab for PH or AFT) and also beats frequentist penalised‑likelihood GH selection. Performance metrics include correct‑structure recovery, variable‑selection precision/recall, and predictive discrimination (C‑index), all of which show substantial gains.

Two real‑world applications illustrate practical utility. In an oncology cohort with ~300 patients and 20 clinical/genomic covariates, the GH model attains the highest posterior probability, revealing that some genes influence survival through the time scale while clinical factors act on the hazard scale. In a reliability dataset of ~400 machine components with 15 sensor readings, the accelerated‑hazards (AH) model is selected, indicating that only a few sensors affect the hazard directly. In both cases, model‑averaged predictions improve over PH‑ or AFT‑only analyses.

The paper concludes by acknowledging limitations such as reliance on parametric log‑location‑scale baselines and the need for extensions to semiparametric baselines, high‑dimensional settings, and non‑local priors. Nonetheless, it delivers the first comprehensive Bayesian solution for joint variable and hazard‑structure selection, offering theoretical guarantees, computational efficiency, and demonstrable empirical advantages.