The Cox-Pólya-Gamma algorithm for flexible Bayesian inference of multilevel survival models

Bayesian Cox semiparametric regression is an important problem in many clinical settings. The elliptical information geometry of Cox models is underutilized in Bayesian inference but can effectively bridge survival analysis and hierarchical Gaussian models. Survival models should be able to incorporate multilevel modeling such as case weights, frailties, and smoothing splines, in a straightforward manner similar to Gaussian models. To tackle these challenges, we propose the Cox-Pólya-Gamma algorithm for Bayesian multilevel Cox semiparametric regression and survival functions. Our novel computational procedure succinctly addresses the difficult problem of monotonicity-constrained modeling of the nonparametric baseline cumulative hazard along with multilevel regression. We develop two key strategies based on the elliptical geometry of Cox models that allows computation to be implemented in a few lines of code. First, we exploit an approximation between Cox models and negative binomial processes through the Poisson process to reduce Bayesian computation to iterative Gaussian sampling. Next, we appeal to sufficient dimension reduction to address the difficult computation of nonparametric baseline cumulative hazards, allowing for the collapse of the Markov transition within the Gibbs sampler based on beta sufficient statistics. We explore conditions for uniform ergodicity of the Cox-Pólya-Gamma algorithm. We provide software and demonstrate our multilevel modeling approach using open-source data and simulations.