The Cox-Polya-Gamma Algorithm for Flexible Bayesian Inference of Multilevel Survival Models

Semiparametric Cox proportional hazards (PH) regression is an important tool used for time to event data (Cox, 1972, 1975). Many extensions of the Cox model have been proposed, such as cure models (or zero-inflation) (Sy & Taylor, 2000), frailty models (or random effects) (Wienke, 2010), smoothing splines (Therneau et al., 2000; Cui et al., 2021), convex set constraints (McGregor et al., 2020), variable selection (Tibshirani, 1997), and other methods. However, nonparametric inference for semiparametric Cox regression involves difficult monotonicity constrained inference of the cumulative hazard (Hothorn et al., 2018). The Bayesian paradigm provides a convenient framework to naturally address many of these variations of the canonical Cox model (Ibrahim et al., 2001; Chen et al., 2006). Bayesian methods can naturally incorporate convex set constraints (Valeriano et al., 2023; Pakman & Paninski, 2014; Maatouk & Bay, 2016) that can be applied to estimation and inference of monotonic cumulative hazard functions. In addition, Bayesian hierarchical modeling can be extended to different types of regression, such as random effects (Wang & Roy, 2018), zero-inflation (Neelon, 2019), smoothing splines (Wand & Ormerod, 2008), and variable selection (Carvalho et al., 2009; Bhadra et al., 2019; Makalic & Schmidt, 2015) with straightforward Gibbs samplers for multilevel Gaussian models. In order to leverage this rich resource, we propose a novel Markov chain Monte Carlo (MCMC) algorithm that allows Cox model parameters to be modeled with hierarchical Gaussian models

The Cox model can be described as locally a Poisson process that is constrained to binary outcomes of recording death or right censoring (Ibrahim et al., 2001; Aalen et al., 2008). In recent years, Polson et al. 2013 proposed a Polya-gamma (PG) Gibbs sampler for logistic and negative binomial (NB) regression that is based in Gaussian sampling. Consequently, the negative binomial process is a useful analog to the Poisson process of Cox regression. We use a standard gamma process Bayesian modification to Cox regression that results in a NB kernel and enables Gaussian-based modeling to incorporate the necessary structure and constraints. By specifying a univariate gamma frailty or gamma prior on the cumulative hazard, we establish an analogous frailty model to locally negative binomial relationship (Aalen et al., 2008; Winkelmann, 2008). Typically, a vague gamma frailty is specified to bridge Poisson and NB models (Kalbfleisch, 1978; Sinha et al., 2003; Duan et al., 2018). The practicality of a vague gamma prior allows access to the PG sampler, benefiting from the use of Gaussian Gibbs samplers. In addition, maintaining Gibbs transitions is crucial as we can take advantage of the marginal transitions to derive succinct Metropolis-Hastings (MH) acceptance probabilities that can be be used to remove bias due to gamma frailty augmentation (Duan et al., 2018).

However, the NB form alleviates much of the algebra associated with the cumulative hazard but does not address the entirety of the baseline hazard rate which has positivity constraints. There are a broad class nonparametric estimators that rely on partitioning the time interval (Nieto-Barajas & Walker, 2002; Nelson, 1972; Aalen, 1978; Kaplan & Meier, 1958; Ibrahim et al., 2001). Evaluating these estimators without considering the computational algorithm often results in difficult MCMC procedures due to monotonic constraints on the baseline cumulative hazard with MH steps that require tuning. There are many tools available to enforce convex set constraints on Gaussian sampling (Valeriano et al., 2023; Pakman & Paninski, 2014; Maatouk & Bay, 2016; Cong et al., 2017; Li & Ghosh, 2015) which we will demonstrate can be used to enforce monotonicity and positivity constraints of Cox regression. Inspired by previous work in nonparametric analysis, we proposed the use of a sufficient statistic partition-based estimator with closed forms for Gibbs sampler conditional distributions. In order to do so, we introduce a novel monotonic spline system with a half-space constraint (Ramsay, 1988; Meyer et al., 2018) to model the baseline log cumulative hazard function that factorizes into beta sufficient statistics in the likelihood. Next, we introduce slice sampler auxiliary variables to remove the sufficient statistics associated with the baseline hazard rate from the algebra (Damien et al., 1999; Mira & Tierney, 2002; Neal, 2003). Our spline system collapses the linear operator of our Markov transitions on to beta sufficient statistics to improve the efficiency of Bayesian computation while retaining the Gaussian structure. This results in a flexible Bayesian procedure we refer to as the Cox-Polya-Gamma (Cox-PG) algorithm that is appealing for being able to adapt multilevel Gaussian inference to Cox regression.

Aside from being able to incorporate prior information into the Cox model, our Cox-PG approach is desirable for situations where inference on the baseline hazard is necessary such as Kaplan-Meier models (Kaplan & Meier, 1958). We also show that case weights such as inverse probability weighting and power priors, can be incorporated into the Cox-PG algorithm through fractal beta kernels (Shu et al., 2021; Ibrahim et al., 2015). In addition, we explore conditions under which the Cox-PG sampler is uniformly ergodic. Consequently, we are able to derive ergodic theory results for our Gibbs algorithm under mild assumptions. Hierarchical models from Gaussian Gibbs samplers are natural extensions of our approach with straightforward modifications to the Cox-PG algorithm. This includes cure models or zero-inflation (Sy & Taylor, 2000; Neelon, 2019), joint models (Wulfsohn & Tsiatis, 1997; Wang et al., 2022), copulas (Wienke, 2010), competing risk (Kalbfleisch & Prentice, 2011) and variable selection (Carvalho et al., 2009; Bhadra et al., 2019; Makalic & Schmidt, 2015) and other multilevel Gaussian models. As a result, we propose a unified Bayesian framework that envelopes many of the models found throughout survival analysis Ibrahim et al. 2001. In contrast to influence function based semiparametric theory (Tsiatis 2006 Chapter 5), the Cox-PG sampler provides a flexible inference procedure for finite sample Cox modeling, addressing a critical gap in many survival analysis settings such as clinical trials (Ildstad et al., 2001).

The rest of the article is organized as follows. In Section 2, we outline the Cox-PG algorithm. In Section 3, we present uniform ergodicity theoretical results. We present additional methods for accelerating MCMC algorithms, removing bias, and variations of Cox semiparametric regression in Section 4 and 5. Simulation studies and real data analysis involving multilevel survival models are presented in Section 6, 7 and 8.

Semiparametric Cox PH regression is closely related to the exponential family with an additional hazard rate term that needs to be accounted for with constrained inference. Exponential family Bayesian models are well studied and have appealing theoretical properties (Choi & Hobert, 2013). In order to emphasize Cox model constraints and its Poisson kernel, we show the derivation of the transformation model for Cox PH regression (Hothorn et al., 2018). The extreme value distribution function of PH models is given as $F_{\mathrm{MEV}}\{A(t)\}=1-\exp[-\exp\{A(t)\}]$ with survival function $S_{\mathrm{MEV}}\{A(t)\}=\exp[-\exp\{A(t)\}]$. In this case, study time $T_{i}$ denotes the censoring or event time for subject $i\in\{1,\dots,N\}$. Let $A(t_{i})=\log\{\Lambda(t_{i})\}+\mathbf{x}^{\top}_{i}\bm{\beta}=\alpha(t_{i})+\mathbf{x}^{\top}_{i}\bm{\beta}=\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}+\mathbf{x}^{\top}_{i}\bm{\beta}=\mathbf{m}_{i}^{\top}\bm{\eta}$. Monotonic function $\alpha(t_{i})=\sum^{J}_{j=1}z_{\alpha,j}(t_{i})u_{\alpha,j}=\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}$ is the baseline log cumulative hazard approximated using additive monotonic splines with the constraints $\mathbf{u}_{\alpha}>\mathbf{0}$ represented through $\mathcal{C}_{0}=\left\{\bm{\eta}:\mathbf{u}_{\alpha}>\mathbf{0}\right\}$. The baseline hazard is $D_{t}\Lambda(t_{i})=D_{t}\exp\left\{\alpha_{0}+\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}=\left\{D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}\exp\left\{\alpha_{0}+\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}$ where the differential operator is represented as $d/dt=D_{t}$ and $\alpha_{0}$ is an intercept incorporated into $\mathbf{x}^{\top}_{i}\bm{\beta}$. The marginal likelihood contribution of the model (Kalbfleisch & Prentice 2011 Chapter 2) is

$$\begin{array}[]{c}\left[D_{t}F_{\mathrm{MEV}}\{A(t_{i})\}\right]^{y_{i}}\left[S_{\mathrm{MEV}}\{A(t_{i})\}\right]^{1-y_{i}}=\left\{D_{t}\alpha(t_{i})\right\}^{y_{i}}\exp\{y_{i}A(t_{i})\}S_{\mathrm{MEV}}\{A(t_{i})\}\end{array}$$

with death (or event) variable $y_{i}$ and right censoring $1-y_{i}$. The hazard function is given by $\left\{D_{t}\alpha(t_{i})\right\}^{y_{i}}\exp\{y_{i}A(t_{i})\}$ following the canonical factorization of hazard-survival functions. The flexibility of this framework is seen by noting that if $\alpha(t)=\alpha_{0}+\log(t)$, the likelihood is an exponential PH regression while if $\alpha(t)=\alpha_{0}+\gamma\log(t)$ we have the Weibull PH model and the parameterization of $\alpha(t)$ determines the hazard family. Now, we set $A(t_{i})=\mathbf{m}^{\top}_{i}\bm{\eta}$, where $\bm{\eta}=\left[\mathbf{u}^{\top}_{\alpha},\bm{\beta}^{\top}\right]^{\top}$ and $\mathbf{m}_{i}=\left[\mathbf{z}^{\top}_{\alpha}(t_{i}),\mathbf{x}^{\top}_{i}\right]^{\top}$, and the resulting generalized likelihood is

$$\begin{array}[]{c}L_{\mathrm{PH}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})=\left[\prod_{i=1}^{N}\left\{D_{t}\mathbf{z}_{\alpha}^{\top}(t_{i})\mathbf{u}_{\alpha}\right\}^{y_{i}}\right]\mathbb{I}(\bm{\eta}\in\mathcal{C}_{0})\exp\left(\mathbf{y}^{\top}\mathbf{M}\bm{\eta}\right)\exp\left\{-\sum_{i=1}^{N}\exp\left(\mathbf{m}_{i}^{\top}\bm{\eta}\right)\right\}.\end{array}$$

The derivative of monotonic splines, $D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})$, is positive and the basis coefficients are constrained to be positive. The design matrix $\mathbf{M}\in\mathbb{R}^{N\times(J+P)}$ is composed of rows $\mathbf{m}_{i}$. Under this construction, slice sampler auxiliary variables are suitable to maintain positive $D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\in\mathbb{R}_{+}\coloneqq(0,\infty)$ (Damien et al., 1999; Mira & Tierney, 2002).

Slice samplers naturally remove positive terms in the likelihood by enforcing constraints in conditional distributions. In order to effectively incorporate a slice sampler, we start by placing a gamma mixture variable $z_{\gamma,i}\sim\mathrm{G}(\epsilon,\epsilon)$ on the Poisson kernel of PH models to obtain a frailty model analogous to negative binomial regression. Note that the first two moments are given as $\mathbb{E}(z_{\gamma,i})=1$ and $\mathrm{Var}(z_{\gamma,i})=\epsilon^{-1}$ and a vague mixture $\epsilon\gg 0$ can be used. As noted in Duan et al. (2018) this limit approximation works well for Poisson processes. In addition, we show how to remove this bias later on in the article through a simple MH step. A quick Laplace transform reveals the modification to the baseline hazard under frailty: $D_{t}\Lambda_{\text{frailty}}(t)=D_{t}\Lambda(t)/\left\{1+\epsilon^{-1}\Lambda(t)\right\}$ (Aalen et al. 2008 Chapter 6) and by sending $\epsilon\rightarrow\infty$ we have the Cox model hazard. We remove the hazard component of the likelihood with multiplication of slice sampling auxiliary uniform variable $\nu_{i}\mid\bm{\eta}\sim\mathrm{U}\{0,D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\}$ for uncensored events, $y_{i}=1$. The resulting likelihood is given as

$$\begin{array}[]{rl}L_{\text{PH}}(\bm{\eta},\mathbf{z}_{\gamma}\mid\mathbf{y},\mathbf{M})\pi(\bm{\nu}\mid\bm{\eta})\pi(\mathbf{z}_{\gamma})=&\left[\prod_{i=1}^{N}\mathbb{I}\left\{\nu_{i}\leq D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}^{y_{i}}\right]\mathbb{I}(\bm{\eta}\in\mathcal{C}_{0})\\ &\times\left[\prod_{i=1}^{N}\left\{\exp\left(\mathbf{m}_{i}^{\top}\bm{\eta}\right)z_{\gamma,i}\right\}^{y_{i}}\right]\pi(\mathbf{z}_{\gamma})\\ &\times\exp\left\{-\sum_{i=1}^{N}\exp\left(\mathbf{m}_{i}^{\top}\bm{\eta}\right)z_{\gamma,i}\right\}\end{array}$$

where $\mathbf{z}_{\gamma}$ is a vector of the $z_{\gamma,i}$ (Wienke 2010 Chapter 3) and $\nu_{i}\geq 0$. We can marginalize out $\mathbf{z}_{\gamma}$ by carrying out the Mellin transform to obtain a negative binomial kernel

$$\begin{array}[]{c}\frac{1}{y_{i}!}\int_{0}^{\infty}\exp(-\lambda_{i}{z_{\gamma}})\lambda_{i}^{y_{i}}z_{\gamma}^{y_{i}}\frac{\epsilon^{\epsilon}}{\Gamma(\epsilon)}z_{\gamma}^{\epsilon-1}e^{-\epsilon z_{\gamma}}dz_{\gamma}=\frac{\Gamma(y_{i}+\epsilon)}{\Gamma(y_{i}+1)\Gamma(\epsilon)}\left(\frac{\epsilon}{\epsilon+\lambda_{i}}\right)^{\epsilon}\left(\frac{\lambda_{i}}{\epsilon+\lambda_{i}}\right)^{y_{i}}.\end{array}$$

Let $\psi_{i}=\mathbf{m}_{i}^{\top}\bm{\eta}-\log(\epsilon)=\mathbf{m}_{i}^{\top}\bm{\eta}+\eta_{\epsilon}$ and $\lambda_{i}=\exp(\mathbf{m}_{i}^{\top}\bm{\eta})$, then the NB likelihood is

$$\begin{array}[]{rl}L_{\text{NB}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\propto&\left[\prod_{i=1}^{N}\left\{D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}^{y_{i}}\right]\mathbb{I}(\bm{\eta}\in\mathcal{C}_{0})\left[\prod_{i=1}^{N}(p_{i})^{y_{i}}(1-p_{i})^{\epsilon}\right]\\ L_{\text{NB}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\pi(\bm{\nu}\mid\bm{\eta})\propto&\left[\prod_{i=1}^{N}\mathbb{I}\left\{\nu_{i}\leq D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}^{y_{i}}\right]\mathbb{I}(\bm{\eta}\in\mathcal{C}_{0})\left[\prod_{i=1}^{N}(p_{i})^{y_{i}}(1-p_{i})^{\epsilon}\right]\end{array}$$

(1)

where $p_{i}=\mathrm{expit}(\psi_{i})=\exp(\psi_{i})/\{1+\exp(\psi_{i})\}$ and we condition on $\epsilon$. We can integrate out $\nu_{i}\geq 0$ to obtain the proportional frailty likelihood. The frailty likelihood is multiplied by constant $\prod_{i=1}^{N}\epsilon^{-y_{i}}$ to get the NB likelihood. After accounting for covariates and conditioning on $\epsilon$, through the Laplace transform, the hazard function is proportional to $D_{t}\Lambda_{\text{frailty}}(t_{i}\mid-)\propto\left\{D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}^{y_{i}}p_{i}^{y_{i}}$, which allows for data augmentation (Wienke 2010 Chapter 3). The survival function under frailty is $(1-p_{i})^{\epsilon}$ resulting in the canonical hazard-survival factorization.

To account for the logit link, we can use a PG auxiliary variable $\omega\mid b,c\sim\mathrm{PG}\left(b,c\right)$, where $b>0$ (Polson et al., 2013). The PG density is

$$\begin{array}[]{c}\pi(\omega\mid b,c)=\frac{1}{2\pi^{2}}\sum_{k=1}^{\infty}\frac{g_{k}}{(k-1/2)^{2}+c^{2}/\left(4\pi^{2}\right)}\end{array}$$

where $g_{k}$ are independently distributed according to $\mathrm{G}(b,1)$. We have the following exponential tilting or Esscher transform (Esscher, 1932) relationship for PG distributions

$$\begin{array}[]{c}\frac{\left(e^{\psi}\right)^{a}}{\left(1+e^{\psi}\right)^{b}}=2^{-b}e^{\kappa\psi}\int_{0}^{\infty}e^{-\omega\psi^{2}/2}\pi(\omega\mid b,0)d\omega,\quad\pi(\omega\mid b,c)=\frac{\exp\left(-c^{2}\omega/2\right)\pi(\omega\mid b,0)}{\mathbb{E}_{\omega\mid b,0}\left[\exp\left(-c^{2}\omega/2\right)\right]}\end{array}$$

(2)

where $\kappa=a-b/2$ and $\pi(\omega\mid b,0)$ denotes a $\mathrm{PG}(b,0)$ density. We use $\omega_{i}\sim\mathrm{PG}\left(y_{i}+\epsilon,0\right)$ and get

$$\begin{array}[]{rl}\frac{\left(e^{\psi_{i}}\right)^{y_{i}}}{\left(1+e^{\psi_{i}}\right)^{y_{i}+\epsilon}}\propto\int_{0}^{\infty}e^{\kappa_{i}\psi_{i}}e^{-\omega_{i}\psi_{i}^{2}/2}\pi\left(\omega_{i}\mid y_{i}+\epsilon,0\right)d\omega_{i}\end{array}$$

to tease out a Gaussian kernel. Note that

$$\begin{array}[]{c}\pi(\omega_{i}\mid y_{i}+\epsilon,\psi_{i})=\cosh^{y_{i}+\epsilon}(|\psi_{i}/2|)\exp\left(-\omega_{i}\psi_{i}^{2}/2\right)\pi(\omega_{i}\mid y_{i}+\epsilon,0),\end{array}$$

can be shown with a Laplace transform (Polson et al., 2013; Choi & Hobert, 2013). We get the NB posterior

$$\begin{array}[]{rl}L_{\text{NB}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\pi(\bm{\nu},\bm{\omega}\mid\bm{\eta})\propto&\left[\prod_{i=1}^{N}\mathbb{I}\left\{\nu_{i}\leq D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}^{y_{i}}\right]\mathbb{I}(\bm{\eta}\in\mathcal{C}_{0})\\ &\times\prod_{i=1}^{N}e^{\kappa_{i}\psi_{i}}e^{-\omega_{i}\psi_{i}^{2}/2}\pi\left(\omega_{i}\mid y_{i}+\epsilon,0\right)\end{array}$$

(3)

where $\kappa_{i}=(y_{i}-\epsilon)/2$. Integrating out $\bm{\nu}$ and $\bm{\omega}$ results in the original NB likelihood. We can write our Gaussian density conditioned on $\eta_{\epsilon}$ through $\begin{array}[]{c}e^{\kappa_{i}\psi_{i}}e^{-\omega_{i}\psi_{i}^{2}/2}\propto e^{\kappa_{i}\varphi_{i}}e^{-\eta_{\epsilon}\omega_{i}\varphi_{i}}e^{-\omega_{i}\varphi_{i}^{2}/2}\end{array}$ with $\varphi_{i}=\mathbf{m}_{i}^{\top}\bm{\eta}$. The conditional distributions, after applying the Esscher transform to both the PG and Gaussian variables, are given as

$$\begin{array}[]{rl}\left(\omega_{i}\mid-\right)&\sim\mathrm{PG}\left(y_{i}+\epsilon,\psi_{i}\right)\\ (\nu_{i}\mid-)&\sim\mathrm{U}\{0,D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\}\\ (\bm{\eta}\mid-)&\sim\mathrm{TN}\left[\left(\mathbf{M}^{\top}\bm{\Omega}\mathbf{M}\right)^{-1}\mathbf{M}^{\top}(\bm{\kappa}-\eta_{\epsilon}\bm{\omega}),\left(\mathbf{M}^{\top}\bm{\Omega}\mathbf{M}\right)^{-1},\mathcal{C}_{\nu}\bigcap\mathcal{C}_{0}\right]\end{array}$$

(4)

where the truncated normal is constrained to set: $\mathcal{C}_{\nu}\bigcap\mathcal{C}_{0}$ and $\mathcal{C}_{\nu}=\bigcap_{\{i:y_{i}=1\}}^{N}\left\{\bm{\eta}:\nu_{i}\leq D_{t}\mathbf{z}^{\top}_{\alpha}(t_{i})\mathbf{u}_{\alpha}\right\}$. The truncated normal density is proportional to the Gaussian kernel multiplied by an indicator for set constraints: $\bm{\eta}\sim\mathrm{TN}(\bm{\mu},\bm{\Sigma},\mathcal{C})$,

$$\begin{array}[]{c}\pi(\bm{\eta}\mid\bm{\mu},\bm{\Sigma},\mathcal{C})=\frac{\exp\left\{-\frac{1}{2}(\bm{\eta}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{\eta}-\bm{\mu})\right\}}{c_{1}(\bm{\mu},\bm{\Sigma},\mathcal{C})}\mathbb{I}(\bm{\eta}\in\mathcal{C})\propto{\exp\left\{-\frac{1}{2}(\bm{\eta}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{\eta}-\bm{\mu})\right\}}\mathbb{I}(\bm{\eta}\in\mathcal{C})\end{array}$$

where $c_{1}(\bm{\mu},\bm{\Sigma},\mathcal{C})=\oint_{\bm{\eta}\in\mathcal{C}}\exp\left\{-\frac{1}{2}(\bm{\eta}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{\eta}-\bm{\mu})\right\}d\bm{\eta}$. Here $\bm{\Omega}=\operatorname{diag}\left(\omega_{1},\cdots,\omega_{N}\right)$, $\omega_{i}$ make up the vector $\bm{\omega}$ and $\kappa_{i}=(y_{i}-\epsilon)/2$ and make up the vector $\bm{\kappa}$. Note that convex set intersection preserves convexity (Boyd & Vandenberghe 2004 Chapter 2), allowing for the exploitation of convex set constrained Gaussian sampling (Maatouk & Bay, 2016). At this point, we have shown how to incorporate hazard rate constraints into the general Gibbs sampler. The priors for $\bm{\eta}$ are flat and we outline hierarchical Gaussian structures later in the article.

We establish uniform ergodicity for our Gibbs sampler and show that the moments of posterior samples for $\bm{\eta}$ exists which guarantees central limit theorem (CLT) results for posterior MCMC averages and consistent estimators of the associated asymptotic variance (Jones, 2004). Our posterior inference of $\alpha(t)$ involves matrix multiplication of coefficients and basis functions, $\bm{\mathcal{Z}}_{\alpha}\mathbf{u}_{\alpha}$ or a decomposition of additive functions and have the same CLT properties. Here $\bm{\mathcal{Z}}_{\alpha}$ is the matrix representation of the monotonic spline bases over the range of the observed covariate data. Suppose $\bm{\eta}^{(n)}$ is a draw from the posterior distribution, then $\alpha^{(n)}(t)=\alpha_{0}^{(n)}+\bm{\mathcal{Z}}_{\alpha}\mathbf{u}^{(n)}_{\alpha}$ are posterior samples of nonparametric effects with MCMC averages of $\widehat{\alpha}(t)=\left(N_{\text{MC}}\right)^{-1}\sum^{N_{\text{MC}}}_{n=1}\alpha^{(n)}(t)$ with an example provided in Figure 1c.

First, we show uniform ergodicity of $\bm{\eta}$ for our Gibbs sampler by taking advantage of the truncated gamma distribution (Choi & Hobert, 2013; Wang & Roy, 2018) and slice sampler properties (Mira & Tierney, 2002). A key feature of this strategy is truncating the gamma distribution at a small $\tau_{0}$ results in useful inequalities related to ergodicity while allowing for a vague prior to be specified for $\tau_{\mathcal{B}}$. We denote the $\bm{\eta}$-marginal Markov chain as $\Psi\equiv\{\bm{\eta}(n)\}^{\infty}_{n=0}$ and Markov transition density (Mtd) of $\Psi$ as

$$\begin{array}[]{c}k\left(\bm{\eta}\mid\bm{\eta}^{\prime}\right)=\int_{\mathbb{R}_{+}^{J}}\int_{\mathbb{R}_{+}}\int_{\mathbb{R}_{+}^{N}}\pi(\bm{\eta}\mid\bm{\omega},{\tau}_{\mathcal{B}},\mathbf{v},\mathbf{y})\pi\left(\bm{\omega},{\tau}_{\mathcal{B}},\mathbf{v}\mid\bm{\eta}^{\prime},\mathbf{y}\right)d\bm{\omega}d{\tau}_{\mathcal{B}}d\mathbf{v}\end{array}$$

where $\bm{\eta}^{\prime}$ is the current state, $\bm{\eta}$ is the next state. Space $\mathbb{R}_{+}^{N}$ is the PG support that contains $\bm{\omega}$, $\tau_{\mathcal{B}}\in\mathbb{R}_{+}$ and $\mathbf{v}\in\mathbb{R}_{+}^{J}$. In addition, placing upper bounds of the baseline hazard ensures finite integrals in the Markov transitions while accommodating vague priors. Under initial constraints $\mathcal{C}_{0}$, we see that $k\left(\bm{\eta}\mid\bm{\eta}^{\prime}\right)$ is strictly positive and Harris ergodic (Tierney, 1994; Hobert, 2011). We show the Mtd of satisfies the following minorization condition: $k\left(\bm{\eta}\mid\bm{\eta}^{\prime}\right)\geq\delta h(\bm{\eta})$, where there exist a $\delta>0$ and density function $h$, to prove uniform ergodicity (Roberts & Rosenthal 2004, Meyn & Tweedie 2012 Chapter 16). Uniform ergodicity is defined as bounded and geometrically decreasing bounds for total variation distance to the stationary distribution in number of Markov transitions $n$, $\left\|K^{n}(\bm{\eta},\cdot)-\Pi(\cdot\mid\mathbf{y})\right\|_{\text{TV}}:=\sup_{A\in\mathscr{B}}\left|K^{n}(\bm{\eta},A)-\Pi(A\mid\mathbf{y})\right|\leq Vr^{n}$ leading to CLT results (Van der Vaart 2000 Chapter 2). Here $\mathscr{B}$ denotes the Borel $\sigma$-algebra of $\mathbb{R}^{P+M}$, $K(\cdot,\cdot)$ is the Markov transition function for the Mtd $k(\cdot,\cdot)$

$$K\left(\bm{\eta}^{\prime},A\right)=\mathbb{P}\left(\bm{\eta}^{(r+1)}\in A\mid\bm{\eta}^{(r)}=\bm{\eta}^{\prime}\right)=\int_{A}k\left(\bm{\eta}\mid\bm{\eta}^{\prime}\right)d\bm{\eta},$$

and $K^{n}\left(\bm{\eta}^{\prime},A\right)=\mathbb{P}\left(\bm{\eta}^{(n+r)}\in A\mid\bm{\eta}^{(r)}=\bm{\eta}^{\prime}\right)$. We denote $\Pi(\cdot\mid\mathbf{y})$ as the probability measure with density $\pi(\bm{\eta}\mid\mathbf{y})$, $V$ is bounded above and $r\in(0,1)$.

1. (C1)

   The gamma prior parameters are constrained as $a_{0}+M/2\geq 1$, $b_{0}>0$ and gamma support is constrained to be $\tau_{\mathcal{B}}\geq\tau_{0}$.

2. (C2)

   The monotonic splines system are comprised of Lipschitz continuous basis functions with coefficients constrained to $\mathbf{u}_{\alpha}\in\left(\mathbf{0},\mathbf{u}_{\alpha}^{+}\right)$, such that $\mathcal{C}_{0}\subseteq\{\bm{\eta}:\mathbf{u}_{\alpha}\in\left(\mathbf{0},\mathbf{u}_{\alpha}^{+}\right)\}$ and $D_{t}\mathbf{z}_{\alpha}^{\top}(t_{i})\mathbf{u}_{\alpha}<D_{t}\mathbf{z}_{\alpha}^{\top}(t_{i})\mathbf{u}_{\alpha}^{+}<\infty$.

We leave the proofs to the supplement. Condition (C1) can be disregarded when considering a Gibbs sampler without mixed effects; the Gibbs sampler is still uniformly ergodic under (C2). Condition (C1) bounds the second order variation of the random effects while allowing for vague priors through the truncated gamma support (Wang & Roy, 2018). Condition (C2) enforces distributional robustness (Blanchet et al., 2022) on how much the baseline hazard can vary, with a broad class of vague and informative priors that can be accommodated. Specifically, it is a common condition that bounds the slopes of the log cumulative hazard function. In large sample theory, condition (C2) or bounding the baseline cumulative hazard is also used to facilitate the dominated convergence theorem and achieve strong consistency (Wang, 1985; McLain & Ghosh, 2013). This proof also guarantees uniform ergodicity for Poisson regression through NB approximations and PG samplers (Zhou et al., 2012; Duan et al., 2018). The following results for coupling time $\left\|K^{n}(\bm{\eta},\cdot)-\Pi(\cdot\mid\mathbf{y})\right\|_{\text{TV}}\leq(1-\delta)^{n}$ is derived from the minorization condition (Jones, 2004). Numerical integration can be used to calculate $\delta$ prior to MCMC sampling to evaluate convergence rates and determine number of MCMC samples $n$ (Jones & Hobert, 2001; Geweke, 1989). More importantly, these results establish conditions for convergence in total variation distance and posterior expectations for constrained nonparametric objects under Bayesian inference.

From Gibbs samplers, we have the useful property of reversible marginal transitions. This allows us to construct a simple Metropolis-Hasting correction to debias our MCMC estimate. The pure Poisson process version of this problem was studied in Duan et al. (2018) with slow mixing but accurate inference using $\epsilon=1000$. A large $\epsilon$ artificially induces an imbalanced logistic regression problem (Zens et al., 2023). Data calibration can be used once a Gibbs sampler for approximate PH regression has been established (Duan et al., 2018). We can use Gibbs sampler (7) as a proposal distribution to sample from the target posterior constructed with the beta kernel Cox semiparametric regression posterior, $L_{\mathrm{PH}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta})$ where $\pi_{0}(\bm{\eta})$ is the prior. We accept the new Gibbs draw from proposal distribution $q(\bm{\eta}\mid\bm{\eta}^{\prime})$ using a MH step with probability

$$\begin{array}[]{c}\min\left\{1,\frac{L_{\mathrm{PH}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta})q\left(\bm{\eta}^{\prime}\mid\bm{\eta}\right)}{L_{\mathrm{PH}}(\bm{\eta}^{\prime}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta}^{\prime})q\left(\bm{\eta}\mid\bm{\eta}^{\prime}\right)}\right\}=\min\left\{1,\prod_{i=1}^{N}\frac{\exp\left\{\lambda_{i}^{\prime}\right\}}{\exp\left\{\lambda_{i}\right\}}\frac{\left\{1+\exp\left(\psi_{i}\right)\right\}^{y_{i}+\epsilon}}{\left\{1+\exp\left(\psi_{i}^{\prime}\right)\right\}^{y_{i}+\epsilon}}\right\}\end{array}.$$

(8)

We denote current state as $\psi_{i}^{\prime}=\log(\lambda_{i}^{\prime})-\log(\epsilon)=\mathbf{m}_{i}^{\top}\bm{\eta}^{\prime}-\log(\epsilon)$ and next state as $\psi_{i}=\log(\lambda_{i})-\log(\epsilon)=\mathbf{m}_{i}^{\top}\bm{\eta}-\log(\epsilon)$. The derivation of this probability can be done quickly by noting that Gibbs samplers are MH algorithms with acceptance probability one. Therefore under the Cox-PG Gibbs sampler for target posterior $L_{\mathrm{NB}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta})$ we have

$$\begin{array}[]{c}\frac{L_{\mathrm{NB}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta})q\left(\bm{\eta}^{\prime}\mid\bm{\eta}\right)}{L_{\mathrm{NB}}(\bm{\eta}^{\prime}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta}^{\prime})q\left(\bm{\eta}\mid\bm{\eta}^{\prime}\right)}=1\implies\frac{q\left(\bm{\eta}^{\prime}\mid\bm{\eta}\right)}{q\left(\bm{\eta}\mid\bm{\eta}^{\prime}\right)}=\frac{L_{\mathrm{NB}}(\bm{\eta}^{\prime}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta}^{\prime})}{L_{\mathrm{NB}}(\bm{\eta}\mid\mathbf{y},\mathbf{M})\pi_{0}(\bm{\eta})}\end{array}.$$

(9)

Substituting (9) into (8), we have a ratio of likelihoods which simplifies to approximately a ratio of survival functions or a difference of cumulants on the log scale. This nested detailed balance structure removes the bias of $\text{G}(\epsilon,\epsilon)$ convolution and allows moderate $\epsilon=100$ to be used to accelerate mixing. The frailty model with $\epsilon=100$ modifies the original Cox model to have more heterogeneity, making it an ideal proposal distribution because proposal distributions should have slightly more variation than their target distribution. Through our experiments, we found $\epsilon=100$ and the MH step to work well with over 90% acceptance rates in various settings. The bias due to frailty in Gibbs samplers (7) can be removed with MH acceptance probability defined in (8). For example, draws of $\bm{\eta}$ are sampled sequentially from (7), each new draw of $\bm{\eta}$ is accepted with probability (8). After removing the bias, we have draws from the PH likelihood based posterior and we can use the log–log link (i.e. $\exp[-\exp\{A(t_{i})\}]$) to map to the survival probabilities.