Bayesian inference for Cox proportional hazard models with partial likelihoods,nonlinear covariate effects and correlated observations

Abstract

We propose a flexible and scalable approximate Bayesian inference methodology for the Cox Proportional Hazards model with partial likelihood. The model we consider includes nonlinear covariate effects and correlated survival times. The proposed method is based on nested approximations and adaptive quadrature, and the computational burden of working with the log-partial likelihood is mitigated through automatic differentiation and Laplace approximation. We provide two simulation studies to show the accuracy of the proposed approach, compared with the existing methods. We demonstrate the practical utility of our method and its computational advantages over Markov Chain Monte Carlo methods through the analysis of Kidney infection times, which are paired, and the analysis of Leukemia survival times with a semi-parametric covariate effect and spatial variation.

1. Introduction

For problems involving time-to-event data, the combination of Cox proportional hazard (Cox PH) models and inference via partial likelihood has been the dominant methodology following its development by Cox. ¹ The Cox PH model assumes that any two subjects’ event hazards are proportional as a function of time, with the ratio depending on unknown covariate effects which are inferred from the observed data. Event times may be correlated within the sample, for example when the response is time to kidney failure for the left and right kidneys from the same subject. Inference that is conducted via partial likelihood does not require assumptions to be made about the form of the baseline hazard. Further, the use of Bayesian inference with the Cox PH model is desirable as this (a) allows the use of substantive prior information on the hazard ratios and (b) provides uncertainty quantification for all parameters of interest in the presence of complex models for the hazard, which would be difficult to achieve otherwise. However, existing methods for approximate Bayesian inference based on integrated nested Laplace approximations (INLA) ² cannot be applied to the Cox PH model with partial likelihood because the Hessian matrix of the log partial-likelihood is fully dense while INLA requires this matrix to be diagonal.

Alternative methods of making Bayesian inference for this kind of survival model have been considered in the literature. Dykstra and Laud ³ considered a fully non-parametric approach for Bayesian survival analysis, where the entire hazard function is modeled with an extended gamma process prior and the posterior distribution is derived to be another extended gamma process. Kim and Kim ⁴ considered Bayesian analysis on Cox PH model on partial likelihood and on full likelihood with an extended gamma process prior for the baseline hazard, and carried out inferences using Markov Chain Monte Carlo (MCMC). Martino et al. ⁵ considered application of the INLA methodology to the Cox PH model, using the full likelihood with baseline hazard modeled semi-parametrically. Kalbfleisch ⁶ derived the partial likelihood to be the limiting posterior when baseline hazards are modeled with non-informative priors, and Sinha et al. ⁷ later extended the result to allow the inclusion of grouped survival data, and implemented partial likelihood-based Bayesian inference with a Gibbs sampling algorithm. Henschel et al. ⁸ proposed a Bayesian inference method using MCMC on the full likelihood, with baseline hazard function modeled either as piecewise constant or as linear combination of B-splines, and they accommodated the inclusion of different types of frailties in their method. Hennerfeind et al. ⁹ developed a general geo-additive Cox PH model that allows the inclusion of components such as non-linear covariate effect, spatial effect and group level frailties, with inference carried out using MCMC on full likelihood and baseline hazards modeled using P-splines. Kneib ¹⁰ generalized the method of Hennerfeind et al. ⁹ to accommodate left truncation, left censoring, and interval censoring of the survival times. Most of the existing methods for Bayesian inference of Cox PH model have been utilizing MCMC method to obtain the posterior, and are based on the full likelihood with an explicit form to model the unknown baseline hazard.

Stringer et al. ¹¹ developed an approximate Bayesian inference methodology for case-crossover model, which can be viewed as a special case of Cox PH model, by applying the approximation strategy of INLA to a log-partial likelihood with a non-diagonal Hessian matrix. Their methodology includes nonlinear covariate effects and yields full posterior uncertainty for the corresponding smoothness parameters, an improvement over existing frequentist methods. The partial likelihood they considered corresponds to one of the simplest special case[s] of the general Cox PH model, and the Hessian matrix of their log-partial likelihood is block-diagonal and sparse. In contrast, the Hessian matrix of log-partial likelihood of Cox PH model is generally fully dense, leading to increased computational burden when compared to the model considered by Stringer et al. ¹¹ Further, they use a manual integration strategy which requires the user to supply their own quadrature points, which requires specialist knowledge to do properly. This limits the practical utility of their method. In order to make approximate Bayesian inferences for the Cox PH model with partial likelihood, an alternative computational strategy is needed.

In this paper, we develop an approximate Bayesian inference method for Cox PH models with partial likelihood, that allows linear and nonlinear covariate effects, spatial effects and frailties for modeling correlations between survival times. The proposed inference method utilizes the Laplace approximation-based strategy^2,11 in a novel way that accommodates the use of partial likelihood. Through two simulation studies, we illustrate the circumstances under which the proposed method yields improved results compared to existing methods based on full likelihood, and demonstrate the accuracy of the posterior approximation and its computational advantages compared to partial likelihood method fit with MCMC. We then further demonstrate the practical utility of the proposed method through the analysis of two classical datasets with correlated survival times, nonlinear covariate effects and spatial variations.

The remainder of this paper is organized as follows. In section 2 we describe the Cox PH model and the method of semi-parametric smoothing that will be used in the inference of the nonlinear covariate effect in this paper. In section 3, we describe our proposed methodology and the introduced improvements to solve the computational challenges presented by the complicated partial likelihood. In section 4 we illustrate advantages of the proposed methodology in two simulation studies and through the analysis of the Kidney catheter data analyzed by McGilchrist and Aisbett ¹² and Leukemia survival data analyzed by Martino et al. ⁵ We conclude in section 5 with a discussion.

2. Model

2.1. A general Cox PH model

Suppose we observe $n$ groups indexed by $i$ , each with $n_i$ observations indexed by $j$ . For example, we may observe $n$ subjects with $n_i$ measurements per subject. Denote the random variable representing the $j$ th survival time in the $i$ th group by $T_{ij}$ , and denote its realization by $t_{ij}$ . Let $c_{ij}$ denote the censoring time for observation $T_{ij}$ such that $T_{ij}$ is not directly observable when $c_{ij}<T_{ij}$ . The observed survival time is $y_{ij}=min\{t_{ij},c_{ij}\}$ . Define $d_{ij}=1$ if $y_{ij}=t_{ij}$ (a survival time) and $d_{ij}=0$ if $t_{ij}>y_{ij}$ (a censoring time). The observations for each $i,j$ are hence denoted by pairs $\mathit{y}=\{(y_{ij},d_{ij}):i\in [n];j\in [n_i]\}$ . The total number of rows in the data set is denoted by $N=\sum _{i=1}^n{n}_i$ . By default, we assume the assumption of independent censoring holds to ensure valid inferences from the Cox PH model, which states that the hazard functions conditional on the same covariate values are the same between the censored and uncensored individuals. ¹³

Define $h_{ij}(t)$ to be the hazard function for the random variable $T_{ij}$ . The Cox PH model assumes $h_{ij}(t)=h_0(t)\text{exp}({\eta }_{ij})$ where $h_0(t)$ is an unknown baseline hazard function that does not depend on the covariates. Kim and Kim ⁴ considered the inference on linear fixed effects with linear predictor defined as ${\eta }_{ij}={{\mathit{x}}_{\mathit{i}\mathit{j}}}^T\mathit{\beta }$ , and they briefly discussed the possibility of generalizing their method to accommodate survival times correlated within groups. Sinha et al. ⁷ proposed a MCMC method for inference with both fixed effects and group level frailties using partial likelihood, but their method does not accommodate nonlinear covariate effect. Dykstra and Laud ³ on the other hand considered to model the entire hazard function $h_{ij}(t)$ nonparametrically, but their method cannot be directly used to quantify the association of a particular covariate with the survival times.

Stringer et al. ¹¹ considered a general linear predictor that accommodates both linear fixed and nonlinear semi-parametric covariate effects, but the type of likelihood they considered is one of the simplest special cases of the general partial likelihood of Cox PH model, and does not allow the estimation of group-level correlation of survival times.

To accommodate nonlinear covariate effects and correlated survival times, we define an additive predictor ${\eta }_{ij}$ which links the covariates for the $ij$ th observation to the survival time $T_{ij}$ :

$$\begin{array}{cc}{\eta }_{ij}& ={{\mathit{x}}_{\mathit{i}\mathit{j}}}^T\mathit{\beta }+{\displaystyle \sum _{q=1}^r{\gamma }_q(u_{qij})+{\xi }_i,i\in [n],j\in [n_i],}\\ {\xi }_i|{\sigma }_{\xi }& \stackrel{iid}{\sim }\mathcal{N}(0,{\sigma }_{\xi }),i\in [n],\\ {\gamma }_q(\cdot )|{\sigma }_q& \stackrel{ind}{\sim }\mathcal{G}\mathcal{P}(0,{\mathcal{C}}_{{\sigma }_q}),q\in [r].\end{array}$$

(1)

Let $\mathit{\eta }=\{{\eta }_{ij}:i\in [n];j\in [n_i]\}$ be the vector of all the additive linear predictors. Here ${\mathit{x}}_{\mathit{i}\mathit{j}}$ is a $p$ -dimensional vector of covariates that are modeled as having linear associations with the log-hazard, and $\mathit{\beta }=({\beta }_1,\dots ,{\beta }_p)$ are regression coefficients. The ${\mathit{u}}_{\mathit{q}}=\{u_{qij}:i\in [n];j\in [n_i]\},q\in [r]$ are covariates whose associations with the log-hazard are modeled semi-parametrically through unknown smooth functions $\{{\gamma }_i,i\in [r]\}$ . The vector of group intercepts $\mathit{\xi }=\{{\xi }_i:i\in [n]\}$ —referred to as frailties coefficients in the context of survival analysis ¹⁴ —are included to model correlation between survival times coming from the same group $i$ . There is no global intercept ${\beta }_0$ as this would be absorbed by $h_0(t)$ . However, in contrast to the case-crossover model considered by Stringer et al., ¹¹ the group-specific intercepts ${\xi }_i$ are estimable in this general Cox PH model, since the case-crossover model only compares the survival times within the same group but the Cox PH model compares all the survival times across all groups.

2.2. Modeling nonlinear covariate effect

The nonlinear covariate effects ${\gamma }_q$ , $q\in [r]$ are modeled semi-parametrically as $r\in \mathbb{N}$ independent zero-mean Gaussian processes, each defined by its covariance function $C_{{\sigma }_q}$ . The covariance functions are each parametrized by a single parameter ${\sigma }_q>0$ . A typical choice of covariance function is the covariance function of the second fold Integrated Wiener process, ¹⁵ which has a connection to cubic smoothing splines. ¹⁶

To infer the infinite-dimensional parameters ${\gamma }_q$ , $q\in [r]$ , Lindgren and Rue ¹⁷ proposed the use of second order random walk model (RW2) to discretize the Integrated Wiener process prior. Miller et al. ¹⁸ showed that the RW2 model proposed by Lindgren and Rue ¹⁷ can be understood as a penalizing basis expansion of the form $\gamma (u)=\sum _{j=1}^d{\varphi }_j(u){\mathrm{\Gamma }}_j$ for each $\gamma (\cdot )$ (dropping the subscript $q$ ), where the random weights $\mathbf{\Gamma }=\{{\mathrm{\Gamma }}_j,j\in [d]\}$ are parameters to be inferred and ${\varphi }_j(\cdot ),j\in [d]$ are fixed, known basis functions which must be chosen. Yue et al. ¹⁹ note that a similar discretization technique to the one used by Lindgren and Rue ¹⁷ yields the B-spline smoothing with integrated derivative penalty of general order proposed by OSullivan, ²⁰ and Wood ²¹ provide an explicit construction of the corresponding precision matrix. In the method of Lindgren and Rue, ¹⁷ the basis function ${\varphi }_j(\cdot )$ is chosen to be the linear B-spline function, and the random weights are given zero-mean Gaussian prior with a banded precision matrix ${\mathrm{\Sigma }}^{-1}(\sigma )$ controlled by a single variance parameter $\sigma$ . In the proposed approach, we use cubic B-splines for the ${\varphi }_j(\cdot )$ and choose the precision matrix ${\mathrm{\Sigma }}^{-1}(\sigma )$ that is obtained by using an integrated second derivative penalty of Wood. ²¹

We put a sum-to-zero constraint such that $\sum _{i=1}^n\gamma (u_i)=0$ in all the following example as default. This sum-to-zero constraint is needed since (a) the overall intercept parameter cannot be identified in partial likelihood and (b) the semi-parametric model we introduced above is invariant to addition of any constant. The constraint will not impact the shape of the nonlinear effect estimate, but only shift its location, which gives it an interpretation as an effect relative to the overall level. This linear constraint will be handled using the approach of conditioning by Kriging, ²² by simulating from the unconditional distributions and applying the corrections to the unconditional samples.

Finally, define the variance parameter vector $\mathit{\theta }=({\theta }_0,\dots ,{\theta }_r)$ where ${\theta }_q=-2\log {\sigma }_q$ , $q=1,\dots ,r$ , and ${\theta }_0=-2\log {\sigma }_{\xi }$ . The variance parameters are given prior distribution $\theta \sim \pi (\theta )$ .

The proposed method can also accommodate the inference of spatial variations, by replacing the covariance function of ${\gamma }_q(u)$ with a proper spatial covariance function. In the following example, we will take the Matern covariance function $M_{\nu }(\Vert .\Vert ;{\sigma }_q,{\rho }_q)$ as the default choice, such that $\mathrm{Cov}({\gamma }_q({\mathit{u}}_{\mathit{i}}\mathit{+}\mathit{h}\mathit{)}\mathit{,}{\gamma }_{\mathit{q}}\mathit{(}{\mathit{u}}_{\mathit{i}}\mathit{)}\mathit{)}\mathit{=}{\mathit{M}}_{\nu }\mathit{(}\Vert \mathit{h}\Vert \mathit{;}{\sigma }_{\mathit{q}}\mathit{,}{\rho }_{\mathit{q}}\mathit{)}$ for any $u_i,h\in {\mathbb{R}}^2$ . For the Matern covariance function $M_{\nu }(\Vert .\Vert ;{\sigma }_q,{\rho }_q)$ , we take the same parametrization as in Brown et al., ²³ with fixed shape parameter $\nu =1$ , and ${\sigma }_q$ and ${\rho }_q$ , respectively represent the marginal standard deviation parameter of ${\gamma }_q$ and its practical correlation range. The corresponding variance parameter ${\mathit{\theta }}_{\mathit{q}}$ is defined as ${\mathit{\theta }}_{\mathit{q}}=(\log ({\sigma }_q),\log ({\rho }_q))$ in this case.

2.3. Partial likelihood

Our inference is carried out via a partial likelihood function. Define the risk set $R_{ij}=\{k,l:y_{kl}\ge y_{ij}\}$ , which contains all the alive and uncensored observations across all the groups at time $y_{ij}$ . Assuming $y_{ij}\ne y_{kl}$ when $(i,j)\ne (k,l)$ , the partial likelihood can be written as follows:

$$\begin{array}{rl}\pi (\mathit{y}|\mathit{\eta })& =\prod _{i=1}^n\prod _{j=1}^{n_i}\{\frac{\exp [{\eta }_{ij}]}{{\displaystyle \sum _{l,k\in R_{ij}}^{}\exp [{\eta }_{lk}]}}{\}}^{d_{ij}},\\ & =\prod _{i=1}^n\prod _{j=1}^{n_i}\{\frac{1}{1+{\displaystyle \sum _{l,k\in R_{ij},(l,k)\ne (i,j)}\exp [{\mathrm{\Delta }}_{lk,ij}]}}{\}}^{d_{ij}},\end{array}$$

(2)

where ${\mathrm{\Delta }}_{lk,ij}={\eta }_{lk}-{\eta }_{ij}$ . Note that $h_0(t)$ does not appear in the partial likelihood, and hence inference may be carried out in the absence of assumptions about $h_0(t)$ .

The partial likelihood (2) can be written in the following form:

$$\begin{array}{cc}\pi (\mathit{y}|\mathit{\eta })& =\prod _{i=1}^n\prod _{j=1}^{n_i}\pi (y_{ij}|\mathit{\eta }),\end{array}$$

(3)

while in order for a model to be compatible with INLA, its likelihood must have the form:

$$\begin{array}{cc}\pi (\mathit{y}|\mathit{\eta })& =\prod _{i=1}^n\prod _{j=1}^{n_i}\pi (y_{ij}|{\eta }_{ij}).\end{array}$$

(4)

Stringer et al. ¹¹ extend this to permit partial likelihoods of the form:

$$\begin{array}{cc}\pi (\mathit{y}|\mathit{\eta })& =\prod _{i=1}^n\prod _{j=1}^{n_i}\pi (y_{ij}|{\mathit{\eta }}_{\mathit{i}}),\end{array}$$

(5)

with ${\mathit{\eta }}_{\mathit{i}}:=\{{\eta }_{ij},j\in [n_i]\}$ , which still does not include (2). Martino et al. ⁵ are able to write the likelihood for their Cox PH model in the form (4) using the full, not partial likelihood (2). Because of this, they require assumptions to be made about the baseline hazard.

Further define ${\mathrm{\Delta }}_{lk,ij}={\eta }_{lk}-{\eta }_{ij}$ in terms of the additive predictors (1). Note that ${\mathrm{\Delta }}_{lk,ij}={\mathrm{\Delta }}_{11,ij}-{\mathrm{\Delta }}_{11,lk}$ for every $(i,j,l,k)$ . To simplify notation, define ${\mathrm{\Delta }}_{ij}={\mathrm{\Delta }}_{11,ij}$ , and note that ${\mathrm{\Delta }}_{11}=0$ . The entire partial likelihood (2) depends on $\eta$ only through $\mathbf{\Delta }=\{{\mathrm{\Delta }}_{ij}:i\in [n];j\in [n_i]\}$ . For the remainder of the paper we reflect this in our notation, writing $\pi (\mathit{y}|\mathrm{\Delta })\equiv \pi (\mathit{y}|\mathit{\eta })$ and defining the log-likelihood $\ell (\mathbf{\Delta };\mathit{y})=\log \pi (\mathit{y}|\mathbf{\Delta })$ .

In existing Laplace approximations for posterior distributions,^2,5,11 the latent parameters $\mathit{W}$ are defined as $\mathit{W}=(\mathbf{\Delta },\mathbf{\Gamma },\mathit{\beta },\mathit{\xi })$ , where the (differenced) linear predictors $\mathbf{\Delta }$ are included as part of the latent parameter vector. Approximate Bayesian inference of this type requires the precision matrix of $\mathit{W}$ to be non-singular, ²⁴ and hence a small noise term ${ϵ}_{ij}\stackrel{iid}{\sim }\text{N}(\mathbf{0},{\tau }^{-1})$ (for some large, fixed $\tau$ ) is added into the model to make the required matrices non-singular. Define the noised linear predictors as

$${\tilde{\mathrm{\Delta }}}_{ij}={\eta }_{11}-{\eta }_{ij}+{ϵ}_{ij},$$

then the resulting precision matrix of (noised) latent parameters $\tilde{\mathit{W}}:=(\tilde{\mathbf{\Delta }},\mathbf{\Gamma },\mathit{\beta },\mathit{\xi })$ is non-singular even when improper prior such as the RW2 prior is used.

Such posterior approximation methods have the advantage that, when the likelihood can be factored out in the form of (4), the resulting log-likelihood Hessian matrix is diagonal and hence efficient to be computed and stored. ² Alternatively, if the likelihood is in the form of (5), the Hessian matrix is still sparse even it is no longer diagonal. ¹¹ However, if one considers applying such approximate Bayesian inference on Cox PH model with partial likelihood, the resulting Hessian matrix will be completely dense and with the number of elements growing quadratically with sample size $N$ . The novel inference method that we introduced in this paper on the other hand does not use the noised predictor, and hence will have a Hessian matrix with a fixed dimension that will not grow with the sample size.

3. Methods

3.1. Approximate Bayesian inference

To deal with the problem of dense Hessian matrix, we proposed a new way to utilize the Laplace approximation, by defining the latent parameter vector to only include the parameters of interest, $\mathit{W}=(\mathbf{\Gamma },\mathit{\beta },\mathit{\xi })$ . As before, the differenced additive linear predictors $\mathbf{\Delta }$ can be computed from $\mathit{W}$ , but $\mathbf{\Delta }$ will not be included as part of $\mathit{W}$ . Note that in our definition, the dimension of this latent parameter vector will be constant, and hence the size of the dense Hessian matrix will be small regardless of the sample size $N$ . This will significantly reduce the memory requirement and computational challenge introduced by computing, storing and factorizing the Hessian matrix, which is necessary for the inferential procedures.

When nonlinear semi-parametric covariate effect is included in the model, $\mathit{W}$ will have a singular precision matrix as the precision matrix of $\mathbf{\Gamma }$ is rank deficient, and hence direct application of the Laplace approximation of Tierney and Kadane ²⁴ will be problematic. This problem of the singular precision matrix is typically fixed by introducing a small Gaussian noise into the additive linear predictors, which makes the precision matrix full rank.^11,2 In our proposed approach, no noises will be added into the linear predictors; instead we will follow the approach of Wood, ²⁵ to fix this problem by adding a very small constant term into the diagonal terms of the precision matrix of $\mathbf{\Gamma }$ , ${\mathrm{\Sigma }}_{\mathrm{\Gamma }}^{-1}$ , which will also result in a full rank precision matrix for $\mathit{W}$ . The proposed modification only shifts the diagonal terms of ${\mathrm{\Sigma }}_{\mathrm{\Gamma }}^{-1}$ by a very small constant, hence will not change any conditional independence structure in the original prior.

Define $\mathit{W}|\mathit{\theta }\sim \text{N}[\mathbf{0},Q_{\theta }^{-1}]$ , where $Q_{\theta }$ is the precision matrix for $\mathit{W}$ . We are interested in estimating and sampling from the joint posterior distribution of the latent parameters:

$$\begin{array}{c}\pi (\mathit{W}|\mathit{y})={\displaystyle \int \pi (\mathit{W}|\mathit{y},\mathit{\theta })\pi (\mathit{\theta }|\mathit{y})d\mathit{\theta }.}\end{array}$$

(6)

We are also interested in the joint posterior distributions of the variance parameters:

$$\begin{array}{c}\pi (\mathit{\theta }|\mathit{y})=\frac{{\displaystyle \int \pi (\mathit{W},\mathit{y},\mathit{\theta })d\mathit{W}}}{{\displaystyle {\int }_{}{\displaystyle {\int }_{}\pi (\mathit{W},\mathit{y},\mathit{\theta })d\mathit{W}d\mathit{\theta }}}}.\end{array}$$

(7)

These are used for point estimates and uncertainty quantification of the variance parameter $\mathit{\theta }$ , and appear as integration weights in (6).

For the posterior of variance parameter (7), we follow the procedure of Stringer et al. ¹¹ to approximate it with its corresponding Laplace approximation ${\tilde{\pi }}_{LA}(\mathit{\theta }|\mathit{y})$ . The posterior of the latent parameter vector (6) is approximated by $\tilde{\pi }(\mathit{W}|\mathit{y})$ defined as:

$$\begin{array}{cc}\tilde{\pi }(\mathit{W}|\mathit{y})& ={\displaystyle \sum _{k=1}^K{\tilde{\pi }}_G(\mathit{W}|\mathit{y},{\mathit{\theta }}^{\mathit{k}}){\tilde{\pi }}_{LA}({\mathit{\theta }}^{\mathit{k}}|\mathit{y}){\delta }_k,}\end{array}$$

(8)

where $\{{\mathit{\theta }}^{\mathit{k}},{\delta }_k{\}}_{k=1}^K$ is a set of nodes and weights corresponding to an adaptive Gauss-Hermite quadrature rule. The ${\tilde{\pi }}_G(\mathit{W}|\mathit{y},{\mathit{\theta }}^{\mathit{k}})$ is a Gaussian approximation for $\pi (\mathit{W}|\mathit{y},{\mathit{\theta }}^{\mathit{k}})$ and the ${\tilde{\pi }}_{LA}({\mathit{\theta }}^{\mathit{k}}|\mathit{y})$ is a Laplace approximation for $\pi ({\mathit{\theta }}^{\mathit{k}}|\mathit{y})$ , which we describe at below.

For any fixed $\mathit{\theta }$ , define

$$\begin{array}{rl}{\hat{\mathit{W}}}_{\theta }=({\hat{\mathbf{\Gamma }}}_{\theta },\hat{\mathit{\beta }},{\hat{\mathit{\xi }}}_{\theta })& ={\text{argmax}}_{\mathit{W}}\log \pi (\mathit{W}|\mathit{\theta },\mathit{y}),\\ H_{\theta }(\mathit{W})& =-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }\mathit{W}\mathrm{\partial }{\mathit{W}}^T}\log \pi (\mathit{W}|\mathit{\theta },\mathit{y}).\end{array}$$

(9)

For the conditional posterior

$$\begin{array}{c}\pi (\mathit{W}|\mathit{\theta },\mathit{y})\propto \exp \{-\frac{1}{2}{\mathit{W}}^T{Q}_{\theta }\mathit{W}+\ell (\mathbf{\Delta };\mathit{y})\},\end{array}$$

(10)

a second-order Taylor expansion of $\log \pi (\mathit{W}|\mathit{\theta },\mathit{y})$ about $\mathit{W}={\hat{\mathit{W}}}_{\theta }$ yields a Gaussian approximation:

$$\begin{array}{cc}\pi (\mathit{W}|\mathit{\theta },\mathit{y})& \approx {\tilde{\pi }}_G(\mathit{W}|\mathit{y},\mathit{\theta })\\ & \propto \text{exp}\{-\frac{1}{2}{(\mathit{W}-{\hat{\mathit{W}}}_{\theta })}^T{H}_{\theta }\left({\hat{\mathit{W}}}_{\theta }\right)(\mathit{W}-{\hat{\mathit{W}}}_{\theta })\}.\end{array}$$

(11)

Define ${\hat{\mathbf{\Delta }}}_{\theta }$ be the differenced additive linear predictors computed at ${\hat{\mathit{W}}}_{\theta }$ . For the joint posterior of the variance parameters, the method of Tierney and Kadane ²⁴ yields a Laplace approximation:

$$\pi (\mathit{\theta }|\mathit{y})\approx {\tilde{\pi }}_{LA}(\mathit{\theta }|\mathit{y})\propto \pi (\mathit{\theta })\{\frac{|Q_{\theta }|}{|H_{\theta }({\hat{\mathit{W}}}_{\theta })|}{\}}^{1/2}\exp \{-\frac{1}{2}{{\hat{\mathit{W}}}_{\theta }}^T{Q}_{\theta }{\hat{\mathit{W}}}_{\theta }+\ell ({\hat{\mathbf{\Delta }}}_{\theta };\mathit{y})\}.$$

(12)

With these approximations available, inference for $\mathit{\theta }$ can be directly obtained using the analytical form of ${\tilde{\pi }}_{LA}(\mathit{\theta }|\mathit{y})$ as in equation (12). Inference for $\mathit{W}$ or its marginal component can be easily obtained if it is possible to get independent samples from $\tilde{\pi }(\mathit{W}|\mathit{y})$ .

To sample from $\tilde{\pi }(\mathit{W}|\mathit{y})$ , note that by equation(6), $\tilde{\pi }(\mathit{W}|\mathit{y})$ is Gaussian mixture distribution with $K$ mixture components each with mixture weight being

$${\varphi }_k={\tilde{\pi }}_{LA}({\mathit{\theta }}^{\mathit{k}}|\mathit{y}){\delta }_k,k\in [K].$$

For a large integer $B$ , we sample independent $\{Z_i{\}}_{i=1}^B$ from $\text{Multinomial}({\varphi }_1,\dots ,{\varphi }_K)$ , and then for each $Z_i\in [K]$ , sample ${\mathit{W}}_{\mathit{i}}$ from ${\tilde{\pi }}_G(\mathit{W}|\mathit{y},{\mathit{\theta }}^{Z_i})$ . The resulting sample $\{{\mathit{W}}_{\mathit{i}}{\}}_{i=1}^B$ then contains $B$ independent draws from $\tilde{\pi }(\mathit{W}|\mathit{y})$ , and all the posterior summaries for $\mathit{W}$ can be obtained using this independent sample.

3.2. Adaptive quadrature and automatic differentiation

Computing the approximations (8) requires choosing a quadrature rule consisting of nodes $\{{\mathit{\theta }}^{\mathit{k}}{\}}_{k=1}^K$ and weights $\{{\delta }_k{\}}_{k=1}^K$ for some chosen $K\in \mathbb{N}$ . Stringer et al. ¹¹ lay a user-chosen grid over a range of $\mathit{\theta }$ that is thought to be plausible, and then compute the Gaussian (11) and Laplace (12) approximations at each point on this grid. This requires the user to choose the location and spread of the grid points, as well as a number $K$ of points that is large enough such that the structure of the resulting posterior approximations is captured. The function $\pi (\mathit{W}|\mathit{y},\mathit{\theta })$ must be optimized, and the Hessian matrix must be stored, for each of these $K$ points. In addition to this strategy requiring the user to have specialist knowledge to implement, it is potentially computationally wasteful since $K$ has to be chosen large enough such that the quadrature points densely cover the range where the majority of mass in $\pi (\mathit{\theta }|\mathit{y})$ lies. In our case, this problem is made more severe by the presence of a dense Hessian. Martino et al. ⁵ use the INLA software and its adaptive quadrature rule which avoids the need for the user to choose points but may still result in a large number of points being used for this same reason.

To mitigate the computational challenges associated with applying a manual quadrature rule for (8), we implement Adaptive Gauss-Hermite Quadrature (AGHQ). This technique has been motivated as a useful tool for Bayesian inference ²⁶ and work has been done to show that it is very accurate when using only a very small number of quadrature points,^28,27 for example attaining $O(N^{-1})$ asymptotic accuracy with $K=3$ and $O(N^{-2})$ with $K=5$ . The use of a small number of quadrature points means only a small number of dense Hessian matrices need to be stored in memory, an improvement over Stringer et al. ¹¹ that is necessary to extend their method to work with the partial likelihood of the Cox PH model.

Computing the AGHQ rule requires computation of the mode of the Laplace approximation:

$$\begin{array}{cc}\hat{\mathit{\theta }}& =\text{argmax}\log {\tilde{\pi }}_{LA}(\mathit{\theta }|\mathit{y}),\end{array}$$

(13)

as well as the Hessian matrix of $\log {\tilde{\pi }}_{LA}(\mathit{\theta }|\mathit{y})$ and its Cholesky decomposition. For the optimization, we use the TMB package, which implements automatic computation of the Laplace approximation and its gradient, which avoids repeated inner optimizations to find ${\hat{\mathit{W}}}_{\theta }$ at different values of $\mathit{\theta }$ . ²⁹ Because of the use of automatic differentiation algorithm, computing the AGHQ rule in the proposed method will not introduce significant computational difficulty.

4. Examples

In this section, we present two simulation studies and two data analysis examples. All the codes are available in the online supplementary materials.

4.1. Simulation studies

We will provide two simulation studies to demonstrate the accuracy of our proposed method and under which situations the accuracy is improved over the existing full likelihood method INLA. Also, we will show that the proposed method provides posterior approximations that are comparable to results of MCMC, with much shorter runtime.

4.1.1. Simulation with sparse frailties

In the first simulation study, we considered the Bayesian inference problem for models with sparse frailties. In other words, survival times were correlated within groups while the number of observations in each group is small. We randomly generated $n=60$ groups, each group with $n_i\equiv m$ observations. The group-level frailties $\{{\xi }_i,i\in [n]\}$ were simulated independently from $\mathcal{N}(0,{\sigma }_{\xi }^2)$ , with varying levels of ${\sigma }_{\xi }$ . Besides the independent frailties, we also assumed there is a covariate $x$ generated from $\mathcal{N}(0,1)$ , with covariate effect $\beta =0.2$ . Among all the survival times generated in this study, $10\text{\%}$ of observations were randomly selected to be right-censored as default. This corresponds to the random censoring mechanism, which directly implies the assumption of independent censoring holds in the simulation. ¹³ In this simulation study, we consider the baseline function to be a simple step function. This choice of piece-wise constant baseline function corresponds to the piece-wise exponential model, which is a type of Cox PH model that is frequently used in the literature. ³⁰ The baseline hazard function in this simulation study is shown in Figure 1(a). We consider six different levels of frailty sparsity in this simulation study by respectively setting $m$ to $1$ , $2$ , $3$ , $4$ , $5$ and $10$ . The parameters of primary interest in this study are the group-level frailties.

Figure 1.

True Baseline Hazards in the two examples in section 4.1. (a) Simple stepwise baseline. (b) Oscillating stepwise baseline. (c) Complicated baseline.

The fixed effect $\beta$ was given a prior $\mathcal{N}(0,1000)$ . The variance parameter ${\sigma }_{\xi }$ was given an Exponential prior with median of 1, which is a penalized complexity prior of Simpson et al. ³¹ The same priors were used for implementations of both our proposed method and INLA. For the adaptive quadrature we used in our inference, the number of grid points for variance parameter was set to be $K=15$ . For the implementation of INLA, we used its first-order random walk model for the baseline hazard run under its default settings. To compare the accuracies between the two methods, we used the metrics of posterior mean square error (MSE) and coverage rates of the $95\text{\%}$ posterior credible intervals, for both the fixed effect parameter and the frailties. All the metrics were computed by averaging through 5000 independent replications.

The comparison metrics when ${\sigma }_{\xi }=1$ are shown in Figure (2). Based on Figure (2)(a) and (b), it can be noticed that our proposed method in general gives more accurate inferential results than INLA for frailty effects, both in terms of smaller MSE and coverage rates closer to the nominal level (i.e. $95\text{\%}$ ), and these differences get larger as the frailties get sparser (smaller $m$ ). The comparison metrics for fixed effect are similar between the two approaches, as shown by Figure (2)(c) and (d). When $m=1$ , INLA has slightly better MSE for the fixed effect parameter $\beta$ than the proposed method, because of the additional correction to the marginal posterior of $\beta$ implemented in its software. ²

Figure 2.

Results for the first simulation with sparse frailty in section 4.1. Left: Plots of MSE of frailties (a) and fixed effect (c) from 5000 independent replications with different $m$ settings, using INLA (red) and the proposed method (blue). Right: Bar-plots of Coverage Rate of 95% posterior credible intervals of frailties (b) and of fixed effect (d) from 5000 independent replications with different $m$ settings, using INLA (red) and the proposed method (blue). The red horizontal line is the nominal rate of 95%. INLA: integrated nested laplace approximation; MSE: mean square error.

For completeness, we implemented the same simulation setting for ${\sigma }_{\xi }=0.4,0.8$ , and $1.3$ . As ${\sigma }_{\xi }$ grows, we found that as a contrast to the other method, inferences from the proposed method are more robust to the sparse frailties, especially in terms of coverage probability for $\mathit{\xi }$ . The detailed comparison metrics for each setting of ${\sigma }_{\xi }$ can be found in Supplemental Appendix B. These results show that unless the group level frailties have trivial variation, the proposed method yields more reliable result than INLA, especially when frailties are sparse.

We then study the performance of the two inference methods when the group sizes $n_i$ ’s are varying across groups instead being fixed at $m$ for all groups. Again the results suggest that as long as the majority of the groups have small group sizes, the proposed approach based on partial likelihood will still be more advantageous in the inference of frailties. The detailed results can be found in Supplemental Appendix A.

To compare the performance at different levels of censoring rates, we use the setting of ${\sigma }_{\xi }=1$ and $m=2$ at censoring rates of 10%, 20% and 40%, and compare the proposed method with INLA as above. The detailed comparison results can be found in Supplemental Appendix C, which shows that the above conclusion still holds for sparse frailty across different levels of censoring rates. As the censoring rates get higher, both methods have worse results for the inference of $\xi$ measured by MSE and coverage rates, but INLA’s full-likelihood approach is more severely affected.

The simulation result above seems to be related to the type of full likelihood INLA utilized in its inference. ⁵ In Cox, ³² Cox pointed out one problem associated with the use of the type of full likelihood considered by Martino et al., ⁵ that is the large number of parameters introduced in order to model the unknown baseline hazard function. In this simulation study, the sample size is only $60m$ , but the latent parameter contains more than $60$ parameters not counting the additional parameters INLA introduced for modeling the baseline hazard, which is likely to become a problem if $m$ is small.

To demonstrate the accuracy of our proposed approximation and the computational advantage compared to existing method, we also fitted the same partial likelihood model using MCMC method, through STAN’s No U-turn Sampler (NUTS), ³³ using two replications where ${\sigma }_{\xi }=1$ and $m=10$ and $2$ , respectively. With four chains each with $35,000$ total iterations with $25,000$ warmups running in parallel, MCMC took a total time of $40.89$ minutes when $m=10$ and $3.56$ minutes when $m=2$ . As a comparison, the proposed method only took $29.64$ seconds when $m=10$ and $1.31$ seconds when $m=2$ , to obtain $10,000$ independent samples of $\mathit{W}$ from the approximate posterior. On average, our method can be run in the time that takes to perform $423$ iterations of MCMC when $m=10$ , and $215$ iterations when $m=2$ .

The difference between posterior distributions yielded by the MCMC method and the approximate posterior obtained by the proposed method is quantified using Kolmogorv–Smirov (KS) statistic, which denotes the maximal absolute difference between the two cumulative posterior distributions, with a larger value indicating less similarity between two distributions. The KS statistics have been computed for both ${\sigma }_{\xi }$ and $\beta$ , as well as the mean and maximal KS statistic for the $60$ frailties. These results are summarized in Table (1). As illustrated in the table, the proposed method yielded approximate posteriors that are very similar to those yielded by MCMC method, but took a significantly shorter time for computation than MCMC. As $m$ decreases from $10$ to $2$ , the KS statistic for ${\sigma }_{\xi }$ increases by 0.088. That is because of the use of Laplace approximation for $\pi ({\sigma }_{\xi }|\mathit{y})$ , which has been known to be less accurate when random effects are sparse. ³⁴ However, the result from our proposed method is still more accurate than the existing Laplace approximation-based method as shown in Figure (2).

Table 1.

KS statistic for each parameter in the first simulation study with sparse frailty in the “Examples section”, to compare the proposed approach with MCMC.

Parameters:	${\sigma }_{\xi }$	$\beta$	${\xi }_i$
Number of Measurements	KS	KS	max KS	mean KS
m = 2	0.104	0.016	0.059	0.039
m = 10	0.016	0.010	0.033	0.019

KS: Kolmogorv-Smirov; MCMC: Markov Chain Monte Carlo.

4.1.2. Simulation with non-smooth baseline

To compare the accuracy of our method with INLA when the smoothness assumption for baseline hazard function is violated, we performed our second simulation study. We generated $n=1000$ uncorrelated data points from a distribution with known hazard function. For the baseline hazard functions, we consider three different settings corresponding to three different levels of wiggliness. Specifically, baseline hazard function is respectively set to simple step function, oscillating step function and an extremely complicated function that switches between linear and constant. Again, the first two piece-wise constant baseline hazards correspond to piecewise Exponential models with different complexity. ³⁰ All of the three baseline hazards $h_0(t)$ are shown in Figure (1). The additive predictor is ${\eta }_i=\gamma (u_i)$ with $\gamma (u)=1.5[\text{sin}(0.8u)+1]$ in all the three simulation settings. We generated the covariates $u$ as $u_i\stackrel{ind}{\sim }\text{Unif}(-6,6)$ , $i\in [n]$ and randomly censored $10\text{\%}$ of all the survival times.

To infer the unknown risk function $\gamma$ , we used the Bayesian cubic B-spline smoothing method mentioned in the section 2 in our proposed method, with $50$ equally spaced knots. For the smoothing method in INLA, we placed the values of $u$ into $50$ discrete bins, and fitted its second-order random walk model for $\gamma$ . ¹⁷ As shown in Stringer et al., ¹¹ the Bayesian semi-parametric smoothing methods we considered here are not sensitive to the choice of number and placement of knots. As before, we implemented INLA under its default setting, with a first-order random walk model for the baseline hazard. This implicitly assumes that $h_0(t)$ is smooth. In contrast, our procedure does not infer $h_0(t)$ , and does not make assumptions about its smoothness. In both of the smoothing methods, the single variance parameter $\sigma$ that controls the smoothness of $\gamma$ , was modeled with an Exponential( $\lambda$ ) prior with $\lambda$ chosen such that $\mathbb{P}(\sigma >2)=0.5$ . For the adaptive quadrature we used in our inference, the number of grid points for variance parameter to is set to be $K=7$ . As in the first simulation study, we compared the accuracy of our proposed method with INLA, using the metrics of MSE computed using posterior mean and coverage rate computed using the $95\text{\%}$ posterior credible interval. These metrics were computed from 1000 independent replications.

The comparison metrics under the three settings of baseline hazard are shown in Figure (3). Figure (3)(a) shows boxplots of the mean squared errors provided by the two methods for $\{\gamma (u_i){\}}_{i=1}^n$ in each replication, for the three different settings of baseline hazard function. Based on the boxplots, the proposed method provides inference results that are at least as accurate as INLA in terms of MSE in all the settings. At the same time, as shown in Figure (3)(b) the proposed method provides consistent coverage rate that is close to the nominal level $95\text{\%}$ in all the settings, but INLA has coverage rate lower than the nominal level when baseline hazard is complicated.

Figure 3.

Results for the second simulation with non-smooth baseline in the section 4.1.2. (a): Plots of MSE from 1000 replications with different baseline settings, using INLA (red) and the proposed method (blue). (b): Bar-plot of Coverage Rate of 95% posterior credible intervals from 1000 replications with different baseline settings, using INLA (red) and the proposed method (blue). The red horizontal line is the nominal rate of 95%. INLA: integrated nested laplace approximation; MSE: mean square error.

While the performance of the proposed method is not affected by the choice of true baseline hazard function, the performance of INLA is sensitive to the true baseline hazard function, as shown in the corresponding boxplots and barplots. This is not unexpected as the full-likelihood used in INLA’s inference implicitly requires that the baseline hazard is smooth enough to be approximated well by its first-order random walk, which will not hold under setting such as Figure (1)(c) where the baseline hazard is varying rapidly as time changes. On the other hand, the inference of our proposed method relies on the partial likelihood, which makes no assumption on the form of the baseline hazard, and hence unaffected by the wiggliness of the baseline hazard in this study.

Again, to assess the accuracy of the proposed posterior approximation, we fitted the same partial likelihood model using MCMC using NUTS with four chains each with 10,000 iterations and 8000 warmups, for an arbitrarily chosen replication. The KS statistic for variance parameter $\sigma$ that controls the smoothness of the inferred $\gamma$ is 0.035, and the maximum and mean KS statistics for the risk function evaluations $\{\gamma (u_i){\}}_{i=1}^n$ are 0.027 and 0.021. We also computed the mean and maximum absolute difference between the posterior mean for $\{\gamma (u_i){\}}_{i=1}^n$ given by the proposed method and by MCMC, which are respectively 0.002 and 0.007. The small KS statistics and absolute difference between posterior mean demonstrate that the proposed posterior approximation is comparable to MCMC. On the other hand, the average runtime of each Markov chain was around 18.7 hours for the MCMC method to get 10,000 iterations, but the proposed method with 10,000 independent samples only requires 28.3 seconds. On average, our method can be run in the time that takes to perform $8$ iterations of MCMC.

4.2. Kidney catheter data

Therneau et al. ³⁵ analyzed a Kidney Catheter dataset using their proposed penalized partial likelihood method. The Kidney Catheter dataset contains 76 times to infection at the point of insertion of a catheter, for $n=38$ patients. An observation for the survival time of a kidney is censored if the catheter is removed for reasons other than an infection. Each patient $i=1,\dots ,n$ forms a group, and the survival times are the time to infection of each patient’s $n_i=2$ kidneys. Because of the small group sizes in the data, this is therefore a practical example of the type of sparse frailty model on which the partial likelihood approach performed better than the full likelihood approach in the simulations of the section 4.1.

We first analyzed this dataset on partial likelihood using the proposed method, with $\mathcal{N}(0,1000)$ priors on the linear covariate effects for age, sex and pre-existing disease types. Subject-specific intercepts ${\xi }_i\stackrel{\stackrel{iid}{\sim }}{\mathcal{N}}(0,{\sigma }_{\xi }^2)$ were included to account for correlation between kidneys from the same subject. We used an $\text{Exponential}$ prior distribution for ${\sigma }_{\xi }$ with median 2. For the adaptive quadrature we used in our inference, the number of grid points for variance parameter to is set to be 15. Ties are handled using the method of Breslow. ³⁶ As a comparison, we also used INLA to fit a Cox PH model to these grouped data on its full likelihood, with the same set of priors as above. The setting for baseline hazard is set to default in INLA’s implementation.

Table (2) shows the results of our procedure compared to the results obtained using INLA and MCMC. Based on the table, our procedure gave different posterior means and reported larger posterior standard deviations compared to INLA, especially for the effects of different disease types. These results are also summarized in Figure (4), where the posteriors obtained from the proposed method are shown to be much closer to MCMC than posteriors obtained from INLA. As we have shown in the section 4.1, when sparse frailties exist, Bayesian inference on partial likelihood tends to be more stable than on full likelihood.

Figure 4.

Posteriors of parameters for the kidney data in the section 4.2. Posterior distribution was obtained using the proposed method(solid line), using MCMC(gray histogram) and using INLA (dashed). (a) Posteriors of σ_ξ. (b) Posteriors of β_Age. (c) Posteriors of β_Sex. (d) Posteriors of β_GN. (e) Posteriors of β_AN. (f) Posteriors of β_PKD. MCMC: Markov Chain Monte Carlo; INLA: integrated nested Laplace approximation.

Table 2.

Estimated means and standard deviations of linear effects by proposed method, INLA and MCMC for the kidney data in the section 4.2.

Methods:		Proposed		INLA		MCMC
Variables/Reference	Levels	Mean	SD	Mean	SD	Mean	SD
Age		0.00484	0.0147	0.00235	0.0130	0.00516	0.0158
Sex/Male	Female	$-$ 1.65	0.473	$-$ 1.64	0.385	$-$ 1.72	0.507
Disease Type/Other	GN	0.172	0.538	0.111	0.474	0.172	0.576
	AN	0.394	0.542	0.519	0.467	0.415	0.573
	PKD	$-$ 1.19	0.826	$-$ 1.06	0.708	$-$ 1.26	0.859

INLA: integrated nested Laplace approximation; MCMC: Markov Chain Monte Carlo.

Again, to quantify the difference between the proposed posterior approximation and the posterior yielded by MCMC, we computed the KS statistics for the five fixed covariate effects, the $38$ patient level frailties and the variance parameter ${\sigma }_{\xi }$ . The KS statistic for ${\sigma }_{\xi }$ is 0.083, and the mean/maximal KS statistics for fixed effects and frailties are 0.034/0.069 and 0.036/0.051, respectively. As a comparison, when using INLA, the KS statistics for ${\sigma }_{\xi }$ is 0.357, and the mean/maximal KS statistics for fixed effects and frailties are respectively 0.102/0.127 and 0.143/0.278.

The KS statistics show that the proposed method give similar posteriors to MCMC method, for most of the parameters. As shown by the KS statistics, the posterior for the variance parameter ${\sigma }_{\xi }$ is approximated less accurately compared to other parameters, because of the problem with Laplace approximation for sparse frailties. ³⁴ This deviation can also be noticed from Figure (4)(a). However, the inference for ${\sigma }_{\xi }$ obtained from the proposed method is still closer to MCMC result than that from INLA for all the parameters. The runtimes are respectively 0.64 seconds for our approach with 10,000 independent samples and 1.98 minutes for MCMC with 35,000 iterations with 25,000 warmups. The number of MCMC iterations that can be obtained with the runtime of the proposed method is only 157.

4.3. Leukemia data

In this example, we implemented our proposed procedure to fit a Cox PH model to the Leukemia data set analyzed by Martino et al. ⁵ as well as previously by Lindgren et al. ³⁷ ; Henderson et al. ³⁸ The dataset contains information from $n=1043$ adult Leukemia patients, with $16\text{\%}$ of observations right-censored. We are interested in quantifying the relationship between survival rate of Leukemia patients with the Townsend deprivation index (tpi), and with their residence locations, controlling effect of the age of the patient, the count of white blood cells at diagnosis (wbc), and sex of the patients.

Based on the model comparison result from Kneib and Fahrmeir, ³⁹ we follow the same analysis in Martino et al. ⁵ to model the effects of age, sex and white blood cell count linearly, and model the effect of the deprivation index (tpi) nonlinearly with the semi-parametric model described in the section 2 using 50 equally spaced knots. Prior distributions $\beta \stackrel{iid}{\sim }\mathcal{N}(0,1000)$ , were used for the linear regression coefficients. The semi-parametric nonlinear effects ${\gamma }_1(\mathbf{t}\mathbf{p}\mathbf{i}):=\{{\gamma }_1({\mathrm{tpi}}_i),i\in [n]\}$ were modeled with the reference constraint $\sum _{i=1}^n{\gamma }_1({\mathrm{tpi}}_i)=0$ . The single variance parameter ${\sigma }_1$ was given an $\text{Exponential}$ prior with a prior median of 2. Furthermore, let $\mathcal{M}$ denote the studied region in this dataset, $\mathit{s}:=\{s_i,i\in [n]\}\subset \mathcal{M}$ denote the residence locations and ${\gamma }_2(\mathit{s}):=\{{\gamma }_2(s_i),i\in [n]\}$ be their corresponding spatial effects. To model the spatially varying Gaussian process ${\gamma }_2(\mathit{s})$ , a Matern covariance function $M_{\nu }(\Vert .\Vert ;{\sigma }_2,\rho )$ is used. For this Matern covariance function, we follow the procedure of Lindgren et al. ³⁷ to fix its shape parameter $\nu$ at 1, and use independent Exponential priors on the standard deviation parameter ${\sigma }_2$ and correlation parameter $\rho$ such that $\mathcal{P}({\sigma }_2>1)=\mathcal{P}(\rho <20\text{km})=0.5$ . For the adaptive quadrature we used in our inference, the number of grid points is set to be $K=4$ .

Figure 5(a) shows the posterior mean and posterior 95 % credible interval of the exponentiated covariate effect of tpi using the proposed approach. Figure 5(b) shows the posterior and prior for the spatial correlation parameter $\rho$ in unit of km. To better understand the spatial variation of risk, we simulate ${\gamma }_2({\mathit{s}}^{\mathbf{\ast }}):=\{{\gamma }_2(s_i^{*}),i\in [L]\}$ for a large $L$ from its corresponding posterior predictive distribution, with $s_i^{*}$ being an arbitrarily chosen point from the region of interest $\mathcal{M}$ . This is done by first simulating ${\gamma }_2(\mathit{s})$ from ${\gamma }_2(\mathit{s})|\mathit{y}$ using the method described in the section 3, and then simulating ${\gamma }_2({\mathit{s}}^{*})$ from ${\gamma }_2({\mathit{s}}^{*})|{\gamma }_2(\mathit{s})$ using the method of Schlather et al. ⁴⁰ For each $s_i^{*}$ , we simulate 100 independent samples, and their corresponding (exponentiated) posterior means and posterior probabilities of having effects exceeding 1.5 are shown in Figure 5(c) and (d). Based on the figure, it can be observed that the exponentiated spatial effect will likely be upper bounded by 1.5 in most of the areas and the isolated areas of high exceeding probability stand out clearly. Posterior exceeding probability at other levels can be found at Supplemental Appendix D.

Figure 5.

Results for the leukemia data in the section 4.3, (a): Posterior mean of (exponentiated) tpi effect (solid) and its 95% credible interval (dashed). (b) Posterior (solid) and prior (dashed) for the spatial $\rho$ parameter. (c): Posterior mean of (exponentiated) effects of residence locations (d): Posterior probability of the (exponentiated) effect of residence locations being larger than 1.5.

There are two major differences between our analysis and the analysis in Martino et al. ⁵ First, our analysis is carried out on partial likelihood using the proposed method, but the analysis in Martino et al. ⁵ is carried out on full likelihood and hence requires assumption on the form of true baseline hazard function. Second, the proposed approach is able to achieve higher resolution estimate for the spatial variation than the approaches taken by Martino et al. ⁵ because we include the full continuously sampled locations ${\gamma }_2(\mathit{s})$ in the latent parameter vector of our analysis without using any approximation to the continuous process. Since ${\gamma }_2(\mathit{s})$ has a dense precision matrix, the method of Martino et al. ⁵ approximates ${\gamma }_2(\mathit{s})$ with a piece-wise constant function with sparse precision matrix. Therefore, Martino et al. ⁵ considers an approximation of the original spatial variation process in order to retain sparseness of the log posterior Hessian, while our model considers the exact Matern process which yields a dense log posterior Hessian, a computational burden that we have already overcome in our procedure.

To contrast the proposed inference method with the method of Stringer et al., ¹¹ we compared the runtimes of Cholesky factorization as well as the memory requirements for storage, of $H_{\theta }({\hat{\mathit{W}}}_{\theta })$ obtained by the two methods at a given $\theta$ . The $H_{\theta }({\hat{\mathit{W}}}_{\theta })$ obtained by the proposed approach takes 0.02 seconds for a single Cholesky factorization, while the $H_{\theta }({\hat{\mathit{W}}}_{\theta })$ obtained by directly applying the approach of Stringer et al. ¹¹ takes 1.08 seconds. As for the memory requirements, $H_{\theta }({\hat{\mathit{W}}}_{\theta })$ obtained by the proposed approach requires 13.7 Mb of memory for its storage, while the one obtained by directly applying the approach of Stringer et al. ¹¹ requires 45.4 Mb to store a single matrix. Since $H_{\theta }({\hat{\mathit{W}}}_{\theta })$ needs to be repeatedly computed and stored during the optimization in equation (9), the computational differences between the two approaches will accumulate during the inference procedures. Therefore, these computational differences demonstrate the necessity of the methodological improvements in our proposed approach as described in the section 3, in order to extend the approach of Stringer et al. ¹¹ to general Cox PH model, especially for complex datasets as above.

5. Discussion

The methodology we proposed in this paper provides a flexible way to carry out Bayesian inference for Cox PH models with partial likelihood, that accommodates the inference for nonlinear covariate effects, spatial variations and correlated survival times. The use of partial likelihood does not require any assumption on the baseline hazard function, which is an advantage over existing approaches for Bayesian inference in this model. We have demonstrated the accuracy and the computational efficiency of our new approach through simulation studies and analysis of two classical datasets in survival analysis. Our proposed method is an appealing option to adopt for the analysis of time-to-event data, when the inference of baseline hazard is not of primary interest.

As shown in our first simulation example, when frailties are sparse, Laplace approximation-based methods tend to yield less accurate approximation for $\pi ({\sigma }_{\xi }|\mathit{y})$ , a problem that is also pointed out in Ogden. ³⁴ For such an application, if the inference of the variance parameter ${\sigma }_{\xi }$ is of primary interest, MCMC type method might be preferred for higher inference accuracies, at the cost of longer runtime. However, we have shown that the proposed approach under such setting yields more accurate inferences than the existing Laplace approximation-based method. Therefore for the analysis with sparse frailties, if the primary interest is the inference of $\pi (\mathit{W}|\mathit{y})$ instead of $\pi ({\sigma }_{\xi }|\mathit{y})$ , or if computational efficiency is of major concern, the proposed method then should be considered as a competent alternative to MCMC method.

The framework of this proposed methodology can be extended to fit more complex models, by modifying the covariance structure of the covariate with nonlinear semi-parametric effect or the covariance structure of the spatial variations. Because we accommodate the dense Hessian matrix of the log-likelihood, our approach could be extended to approximate Bayesian inference for other models with a dense Hessian matrix. We leave such extensions to future work. A R package that implements our proposed method called abcoxp is available at https://github.com/AgueroZZ/abcoxpGitHub.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802221134172 - Supplemental material for Bayesian inference for Cox proportional hazard models with partial likelihoods, nonlinear covariate effects and correlated observations

Supplemental material, sj-pdf-1-smm-10.1177_09622802221134172 for Bayesian inference for Cox proportional hazard models with partial likelihoods, nonlinear covariate effects and correlated observations by Ziang Zhang, Alex Stringer, Patrick Brown and Jamie Stafford in Statistical Methods in Medical Research