A Bayesian semiparametric mixed effects proportional hazards model for clustered partly interval-censored data

Abstract

Clustered partly interval-censored survival data naturally arise from many medical and epidemiological studies. We propose a Bayesian semiparametric approach for fitting a mixed effects proportional hazards (PH) model to clustered partly interval-censored data. The proposed method allows for not only a random intercept as most frailty models do for clustered survival data, but also random effects of covariates. We assume a normal prior for each random intercept/random effect, seeing the instability of a gamma prior for a frailty in this situation. Simulation studies with data generated from both mixed effects PH model and mixed effects accelerated failure times model are conducted, to evaluate the performance of the proposed method and compare it with the three methods currently available in the literature. The application of the proposed approach is illustrated through analyzing the progression-free survival data derived from a phase III metastatic colorectal cancer clinical trial.

Keywords

Bayesian semiparametric clustered partly interval-censored data spline approximation data augmentation mixed effects proportional hazards model

1 Introduction

Partly interval-censored data often occur in health-related studies that involve periodic examinations for symptoms of interest such as oncology clinical trials, epidemiological studies, dental research, and studies of infectious diseases. With partly interval-censored data, the exact time-to-event is known for some subjects, while is only known to be within certain time intervals for the others. It is a combination of exact failure times and general interval-censored (Huang and Wellner, 1997) failure times; or equivalently, a combination of exact, left-censored, interval-censored, and right-censored failure times. For instance, in cancer clinical trials, progression-free survival, defined as time from randomization to disease progression or death due to any cause, is actually partly interval-censored as disease progression needs to be detected by periodic imaging assessments. Also disease-free survival, defined as the length of time a patient stays free of a disease or cancer after a particular treatment, is also partly interval-censored as the recurrence of cancer needs to be detected by imaging assessments too. Some of the main methodological papers for partly interval-censored data are Huang, (1999), Kim, (2003), Komáek and Lesaffre (2007), Zhao et al. (2008), Gao et al. (2017), Zhou and Hanson (2018a), Pan et al. (2020), and Pan and Cai (2020).

When we have patients treated at different medical centres or recurrent events within different patients, these patients or recurrent events may be dependent, that is clustered by centre or by patient. This requires the inclusion of a random intercept/random effects in our survival models to account for the dependency within clusters. Traditionally survival models with a random intercept/random effects are for clustered right-censored data. Klein (1992) proposed a semiparametric estimation method of the proportional hazards (PH) model with a random intercept based on the expectation-maximization (EM) algorithm. Aslanidou et al. (1998) proposed a Bayesian PH model with a shared gamma frailty. Sargent (1998) proposed a Bayesian PH model with a normal random intercept. Ripatti and Palmgren (2000) proposed a mixed effects PH model with the random effects following a multivariate normal distribution. Only during recent years, has there been growing attention to clustered interval-censored data. Goethals et al. (2009) developed a shared gamma frailty model where they integrated out the frailty analytically. Pan et al. (2015) proposed a Bayesian PH model for clustered interval-censored data with multiple gamma frailties. Yavuz and Lambert (2016) developed a Bayesian semiparametric shared-frailty model for clustered interval-censored data, assuming a normal and a flexible P-splines form for the frailty distribution. Chen et al. (2020) developed a multiple imputation approach to convert clustered interval-censored data to right-censored data and analyze it under the additive hazards model. As we can see, most of the methods for both clustered right-censored data and clustered interval-censored data contain only a shared-frailty, or equivalently, a random intercept.

There are only three publications in the literature that we are aware of for modelling clustered partly interval-censored data. Both Komáek and Lesaffre (2007) and Komáek et al. (2007) developed a Bayesian mixed effects accelerated failure times (AFT) model for clustered partly interval-censored data, with the random effects following a multivariate normal prior. They differ in how the error term (equivalently, baseline distribution) is modelled. The accompanying R functions that fit these two models are bayessurvreg1 and bayessurvreg2 in the R package bayesSurv Komáek (2018). Zhou and Hanson (2018a) developed a Bayesian semiparametric model with a random intercept for clustered partly interval-censored data, under PH, proportional odds, and AFT models. The accompanying R function that fits this model is survregbayes in the R package spBayesSurv Zhou and Hanson (2018b).

In this article, we introduce an efficient and easy-to-implement Bayesian semiparametric approach for fitting a mixed effects PH model to clustered partly interval-censored data. The main differences between the proposed method and the three Bayesian methods mentioned above are: (1) Komáek and Lesaffre (2007) fit a mixed effects AFT model with the error term specified as a normal mixture where the number of components follows a Poisson distribution. Komáek et al. (2007) fit a mixed effects AFT model with the error term specified as a mixture of G-splines (normal densities with equidistant means and constant variances) where the number of components is user-specified. While we fit a mixed effects PH model with the cumulative baseline hazard function specified as a linear combination of basis I-splines. (2) Both Komáek and Lesaffre (2007) and Komáek et al. (2007) assume a multivariate normal prior for the vector of random effects, which forces all random effects to jointly follow a multivariate normal distribution as posterior. This distributional constraint may limit the exploration of sampling space during the Markov chain Monte Carlo (MCMC) updates of the random effects and the associated fixed effects. We have tried a multivariate normal prior for the random effects under the proposed method, and the estimated fixed effect with a corresponding random effect has a relatively large bias. (3) Zhou and Hanson (2018a)’s method only allows for a random intercept which assumes a normal prior. A mixed effects model with a general Q-variate random effects covariate vector is in demand to also account for random variations of treatment effect and risk factors. Such a mixed effects PH model has not yet been considered in the partly interval-censored data literature. (4) We have no way to evaluate how Zhou and Hanson (2018a)’s method would perform if random effects of covariates were included. In their MCMC algorithm, the vector of regression coefficients is drawn from a multivariate normal with mean equal to the current value and variance equal to sample variance with a small adjustment. If we imagine this sampler were also to be used for random effects, then it may face the same distributional constraint as in the Komáek and Lesaffre (2007)’s method and the Komáek et al. (2007)’s method. (5) Adding random effects of covariates to a random intercept/shared-frailty only model poses non-trivial challenges to both method development and MCMC algorithm computation. This will be demonstrated in Section 3, where the Komáek and Lesaffre (2007)’s method and the Komáek et al. (2007)’s method perform well when only a random intercept is included. However, the MCMC chain of a regression coefficient will not converge if a corresponding random effect is added.

The current work is related to our previous work (Pan et al., 2015, 2020; Pan and Cai, 2020), in the sense that they all use a monotone spline to flexibly approximate the cumulative baseline hazard function and use a two-step data augmentation based on a nonhomogeneous Poisson process to facilitate posterior computation. Pan et al. (2015) proposed a multiple-frailty model for clustered interval-censored data with frailty selection. Pan et al. (2020) proposed a Bayesian semiparametric PH model for partly interval-censored data. Pan and Cai (2020) proposed a PH model with spatial frailty for spatially dependent partly interval-censored data. Here our main advancement is to modify the model from Pan et al. (2020) by inclusion of random intercept and random effects, such that heterogeneity of baseline survival and heterogeneity of covariate effects are accounted for.

The article is organized as follows. Section 2 presents the proposed Bayesian semiparametric approach under the mixed effects PH model. Section 3 presents two simulation studies that compare the performance of the proposed method versus current models under two data-generating settings. In Section 4, we illustrate the application of the proposed method using the progression-free survival data derived from a phase III metastatic colorectal cancer clinical trial. Finally, Section 5 provides conclusions and discussions.

2 The model

2.1 Data structure

Suppose there are $i = 1, \dots, I$ clusters, with $n_{i}$ . subjects within cluster $i$ . In cluster $i$ , suppose failure times are exactly known for the first $n_{i 1}$ subjects, denoted as $T_{i j}$ , $j = 1, \dots, n_{i 1}$ . But failure times are only known to be within certain time intervals for the other $n_{i} - n_{i 1}$ subjects, denoted as $(L_{i j}, R_{i j}]$ . $j = n_{i 1} + 1, \dots, n_{i}$ . Here $L_{i j}$ can be $0$ which corresponds to a left-censored observation and $R_{i j}$ can be $\infty$ which corresponds to a right-censored observation. So partly interval-censored data encompass exact, left-censored, interval-censored, and right-censored observations. We assume that failure times and examination times are independent given covariates.

2.2 Model and likelihood

We propose a Bayesian semiparametric approach for analyzing clustered partly interval-censored data under the mixed effects PH model. Specifically, it is assumed that the hazard function for the $j$ th subject in the $i$ th cluster (denoted as subject $[i, j]$ ) conditional on the frailty $ϕ_{i}$ of cluster $i$ is:

λ (t_{i j} | x_{i j}, ω_{i j}, ϕ_{i}) = λ_{0} (t_{i j}) exp (β' x_{i j} + ϕ_{i}' ω_{i j}),

(2.1)

where $λ_{0} (t_{i j})$ is the baseline hazard function, $β = \{β_{1}, \dots, β_{P}\}$ is the vector of fixed effects, $ϕ_{i} = \{ϕ_{i 0}, \dots, ϕ_{i Q}\}$ is the vector of random effects of the $i$ th cluster, $x_{i j} = (x_{i j 1}, \dots, x_{i j P})$ is the vector of covariates with fixed effects, and $ω_{i j} = (1, ω_{i j 1}, \dots, ω_{i j Q})$ is the vector of covariates with random effects. Note that the first element of $ω_{i j}$ is 1 which corresponds to a random intercept $ϕ_{i 0}$ and the predictors in $ω_{i j}$ are either equivalent to or a subset of $x_{i j}$ .

For an exact observation $[i, j]$ , failure time $T_{i j}$ is directly observed, and its likelihood function conditional on cluster $i$ ’s random effects $ϕ_{i}$ is

\begin{array}{l} L_{1, i j} (λ_{0} (\cdot), β, ϕ_{i}) = f (t_{i j} | x_{i j}, ω_{i j}, ϕ_{i}) = \\ λ_{0} (t_{i j}) exp (β' x_{i j} + ϕ_{i}' ω_{i j}) exp \{- Λ_{0} (t_{i j}) exp (β' x_{i j} + ϕ_{i}' ω_{i j})\}, \end{array}

(2.2)

where $Λ_{0} (t) = \int_{0}^{t} λ_{0} (s) d s$ is the cumulative baseline hazard function.

For a general interval-censored observation $[i, j]$ , a time interval $(L_{i j}, R_{i j}]$ which contains the failure time $T_{i j}$ is observed, and its likelihood function conditional on cluster $i$ ’s random effects $ϕ_{i}$ is

\begin{array}{l} L_{2, i j} (λ_{0} (\cdot), β, ϕ_{i}) = {\{F (R_{i j} | x_{i j}, ω_{i j})\}}^{δ_{i j 1}} \\ {\{F (R_{i j} | x_{i j}, ω_{i j}) - F (L_{i j} | x_{i j}, ω_{i j})\}}^{δ_{i j 2}} {\{1 - F (L_{i j} | x_{i j}, ω_{i j})\}}^{δ_{i j 3}}, \end{array}

(2.3)

where $F (\cdot | x, ω)$ is the cumulative distribution function of failure time $T$ given $x$ and $ω$ , and $δ_{1}$ , $δ_{2}$ , $δ_{3}$ are the left-, interval-, and right-censoring indicators.

Assuming that subjects are independent conditional on the random effects $ϕ_{i}$ ’s, the overall conditional likelihood function is:

L (λ_{0} (\cdot), β, ϕ_{i}) = \prod_{i = 1}^{I} \{\prod_{j = 1}^{n_{i 1}} L_{1, i j} (λ_{0} (\cdot), β, ϕ_{i}) \prod_{j = n_{i 1} + 1}^{n_{i}} L_{2, i j} (λ_{0} (\cdot), β, ϕ_{i})\} .

(2.4)

2.3 Estimation of

Λ_{0} (t)

and

λ_{0} (t)

Following our previous work (Pan et al., 2015, 2020; Pan and Cai, 2020), we approximate the unspecified cumulative baseline hazard function $Λ_{0} (t)$ semiparametrically with a linear combination of basis I-splines defined in Ramsay (1988):

Λ_{0} (t) = \sum_{l = 1}^{K} γ_{l} I_{l} (t),

(2.5)

where $\{γ_{l}\}$ is a set of non-negative coefficients and $\{I_{l} (t)\}$ is a set of basis I-splines. To construct the set of basis I-splines, we need to specify the degree of the basis I-splines (1 $=$ Linear, 2 $=$ Quadratic, etc.) and the knots. Then the number of basis I-splines $K$ $=$ Degree of basis I-splines $+$ Number of knots $-$ 2. The set of basis I-splines are a set of monotonic polynomial splines that increase from 0 to 1. So Equation (2.5) gives a non-negative and non-decreasing spline that is consistent with the properties of $Λ_{0} (t)$ . As a general rule of thumb, we recommend setting the degree to be 2 or 3 for adequate smoothness and choosing 5 to 15 knots that are equally spaced or at the percentiles for adequate modelling flexibility, based on our extensive experiences.

Based on the fact that the baseline hazard function $λ_{0} (t)$ is the derivative of $Λ_{0} (t)$ , we take the derivative of Equation (2.5) and get:

λ_{0} (t) = \sum_{l = 1}^{K} γ_{l} M_{l} (t) .

(2.6)

Per definition in Ramsay (1988), $\{M_{l} (t)\}$ is called a set of basis M-splines which are the derivatives of their corresponding set of basis I-splines, that is, $M_{l} (t) = \frac{d I_{l} (t)}{d t}$ .

2.4 Data augmentation

It would be difficult to draw MCMC samples from the posteriors derived directly from the conditional likelihood function in Equation (2.4). To facilitate posterior computation, we construct the following data augmentations in order to obtain more posterior distributions of standard forms.

Plugging the I-splines approximation of $Λ_{0} (t)$ (Equation (2.5)) and the M-splines approximation of $λ_{0} (t)$ (Equation 2.6) into Equation (2.2), the likelihood function of an exact observation $[i, j]$ becomes:

\begin{array}{l} L_{1, i j} (λ_{0} (\cdot), β, ϕ_{i}) = \\ \{\sum_{l = 1}^{K} γ_{l} M_{l} (t)\} exp (β' x_{i j} + ϕ_{i}' ω_{i j}) exp \{- \{\sum_{l = 1}^{K} γ_{l} I_{l} (t)\} exp (β' x_{i j} + ϕ_{i}' ω_{i j})\}, \end{array}

It is impossible to extract a specific $γ_{l}$ from the first term for its posterior computation as it is involved in a summation form. We introduce multinomial latent variables $u_{i j} = (u_{i j 1}, u_{i j 2}, \dots, u_{i j K}) \sim Multinomial (1; \frac{1}{K}, \frac{1}{K}, \dots, \frac{1}{K})$ , then we would have the augmented data likelihood of an exact observation $[i, j]$ :

\begin{array}{l} L_{1 a u g, i j} (λ_{0} (\cdot), β, ϕ_{i} | u_{i j}) = \\ \{K \prod_{l = 1}^{K} {(γ_{l} M_{l} (t_{i j}))}^{u_{i j l}}\} exp (β' x_{i j} + ϕ_{i}' ω_{i j}) exp \{- \{\sum_{l = 1}^{K} γ_{l} I_{l} (t)\} exp (β' x_{i j} + ϕ_{i}' ω_{i j})\}, \end{array}

from which we can extract $γ_{l}$ for its posterior computation. Integrating out the latent variables ${\{u_{i j l}\}}_{l = 1}^{K}$ in $K \prod_{l = 1}^{K} {(γ_{l} M_{l} (t_{i j}))}^{u_{i j l}}$ will return $\sum_{l = 1}^{K} γ_{l} M_{l} (t_{i j})$ .

For the general interval-censored observations part, we introduce a two-step data augmentation. Let ${N (t) : t > 0}$ be a nonhomogeneous Poisson process with cumulative intensity function $Λ_{0} (t) exp (β' x + ϕ' ω)$ . Define $T = inf {t : N (t) > 0}$ , time of first occurrence in the Poisson process. Then we have $P (T > t | ϕ) = P (N (t) = 0 | ϕ) = exp \{- Λ_{0} (t) exp (β' x + ϕ' ω)\}$ , indicating that $T$ indeed follows the mixed-effects PH model (Equation (2.1)). Now define two time points $t_{1}$ and $t_{2}$ such that $0 < t_{1} < t_{2}$ . Then the latent variables $Z \equiv N (t_{1}) \sim Poi (Λ_{0} (t_{1}) exp (β' x + ϕ' ω))$ , $W \equiv N (t_{2}) - N (t_{1}) \sim Poi ((Λ_{0} (t_{2}) - Λ_{0} (t_{1})) exp (β' x + ϕ' ω))$ , and they are independent. For a left-censored observation, we set $t_{1} = R < t_{2}$ . For an interval-censored observation, we set $t_{1} = L$ and $t_{2} = R$ . For a right-censored observation, we set $t_{1} < t_{2} = L$ . Then we have for left-censored data, $Z > 0$ and $W$ can take any value. For interval-censored data, $Z = 0$ and $W > 0$ . For right-censored data, $Z = W = 0$ . The augmented data likelihood for a general interval-censored subject $[i, j]$ is:

\begin{array}{l} L_{2 a u g 1, i j} (λ_{0} (\cdot), β, ϕ_{i} | Z_{i j}, W_{i j}) = {\{p (Z_{i j} | μ_{i j 1}) 1 (Z_{i j} > 0)\}}^{δ_{i j 1}} \\ {\{p (Z_{i j} | μ_{i j 1}) p (W_{i j} | μ_{i j 2}) 1 (Z_{i j} = 0) 1 (W_{i j} > 0)\}}^{δ_{i j 2}} \\ {\{p (Z_{i j} | μ_{i j 1}) p (W_{i j} | μ_{i j 2}) 1 (Z_{i j} = 0) 1 (W_{i j} = 0)\}}^{δ_{i j 3}}, \end{array}

(2.7)

where $p (\cdot | μ)$ denotes the Poisson probability mass function with mean $μ$ , $μ_{i j 1} = Λ_{0} (t_{i j 1}) exp (β' x_{i j} + ϕ_{i}' ω_{i j})$ , and $μ_{i j 2} = (Λ_{0} (t_{i j 2}) - Λ_{0} (t_{i j 1})) exp (β' x_{i j} + ϕ_{i}' ω_{i j})$ . Integrating out the latent variables $Z_{i j}$ and $W_{i j}$ will return $L_{2, i j}$ in Equation (2.3).

Furthermore, based on the additive property of Poisson distribution and the additive form of $Λ_{0} (t)$ in Equation (2.5), decompose $Z$ and $W$ respectively into $K$ mutually independent Poisson latent variables ${\{Z_{l}\}}_{l = 1}^{K}$ and ${\{W_{l}\}}_{l = 1}^{K}$ , such that $Z = \sum_{l = 1}^{K} Z_{l}$ with $Z_{l} \sim Poi (γ_{l} I_{l} (t_{1}) exp (β' x + ϕ' ω))$ and $W = \sum_{l = 1}^{K} W_{l}$ with $W_{l} \sim Poi (γ_{l} (I_{l} (t_{2}) - I_{l} (t_{1})) exp (β' x + ϕ' ω))$ . Then the further augmented data likelihood for a general interval-censored subject $[i, j]$ is:

\begin{array}{l} L_{2 a u g 2, i j} (λ_{0} (\cdot), β, ϕ_{i} | Z_{i j l}' s, W_{i j l}' s) = {[\{\prod_{l = 1}^{K} p (Z_{i j l} | μ_{i j l 1})\} 1 (\sum_{l = 1}^{K} Z_{i j l} > 0)]}^{δ_{i j 1}} \\ {[\{\prod_{l = 1}^{K} p (Z_{i j l} | μ_{i j l 1}) p (W_{i j l} | μ_{i j l 2})\} 1 (\sum_{l = 1}^{K} Z_{i j l} = 0) 1 (\sum_{l = 1}^{K} W_{i j l} > 0)]}^{δ_{i j 2}} \\ {[\{\prod_{l = 1}^{K} p (Z_{i j l} | μ_{i j l 1}) p (W_{i j l} | μ_{i j l 2})\} 1 (\sum_{l = 1}^{K} Z_{i j l} = 0) 1 (\sum_{l = 1}^{K} W_{i j l} = 0)]}^{δ_{i j 3}}, \end{array}

(2.8)

where $μ_{i j l 1} = γ_{l} I_{l} (t_{i j 1}) exp (β' x_{i j} + ϕ_{i}' ω_{i j})$ and $μ_{i j l 2} = (γ_{l} I_{l} (t_{i j 2}) - γ_{l} I_{l} (t_{i j 1})) exp (β' x_{i j} + ϕ_{i}' ω_{i j})$ . Integrating out the latent variables ${\{Z_{i j l}\}}_{l = 1}^{K}$ and ${\{W_{i j l}\}}_{l = 1}^{K}$ will return $L_{2 a u g 1, i j}$ in Equation (2.7). As we can see, the likelihood function for the part of general interval-censored observations now becomes a product of Poisson probability mass functions, which will enable more straightforward posterior computation.

The overall augmented data likelihood then becomes:

\begin{array}{l} L_{a u g} (λ_{0} (\cdot), β, ϕ_{i} | u_{i j}, Z_{i j l}^{'} s, W_{i j l}^{'} s) = \\ \prod_{i = 1}^{I} \{\prod_{j = 1}^{n_{i 1}} L_{1 a u g, i j} (λ_{0} (\cdot), β, ϕ_{i} | u_{i j}) \prod_{j = n_{i 1} + 1}^{n_{i}} L_{2 a u g 2, i j} (λ_{0} (\cdot), β, ϕ_{i} | Z_{i j l}^{'} s, W_{i j l}^{'} s)\}, \end{array}

which will enable efficient sampling in the posterior computation to be presented in the Appendix. The latent variables, spline coefficients, and hyperparameters can all be sampled from standard distributions. Only the posteriors of regression coefficients and random effects are of non-standard form, which will be sampled using the Metropolis-Hastings algorithm Hastings (1970).

2.5 Prior specification and posterior computation

The detailed MCMC algorithm and formulas of the posterior distributions are presented in the Appendix. The latent variables $Z_{i j}$ and $W_{i j}$ are sampled from zero-truncated Poisson distributions. The latent variables $(Z_{i j 1}, \dots, Z_{i j K})$ and $(W_{i j 1}, \dots, W_{i j K})$ are sampled from multinomial distributions. The latent variables $(u_{i j 1}, \dots, u_{i j K})$ are sampled from a multinomial distribution too. For spline coefficients, an exponential prior $γ_{l} \sim Exp (η)$ together with a gamma hyperprior $η \sim Ga (a_{η}, b_{η})$ are assumed. This leads to conjugate gamma posteriors for both $γ_{l}$ and $η$ . For the regression coefficient of a numeric covariate, a normal prior $β_{p} \sim N (0, σ_{0}^{2})$ is assumed. The corresponding posterior is not conjugate and the Metropolis-Hastings algorithm is used for sampling from the posterior. For the regression coefficient of a categorical covariate with $c$ levels, we denote it using $c - 1$ dummy variables and sample the exponentiated parameter $ζ_{p} = exp (β_{p})$ for each dummy variable. A gamma prior $ζ_{p} \sim Ga (a_{ζ}, b_{ζ})$ leads to a conjugate posterior. Then we transform $ζ_{p}$ back to $β_{p}$ . For each random effect, we assume a normal prior $ϕ_{i q} \sim N (0, τ_{q}^{- 1})$ , $i = 1, \dots, I$ together with a gamma hyperprior $τ_{q} \sim Ga (a_{τ}, b_{τ})$ . This leads to a conjugate posterior for $τ_{q}$ . The posterior of $ϕ_{i q}$ is not conjugate and the Metropolis-Hastings algorithm is used for sampling from the posterior.

2.6 Model comparison criteria

Spiegelhalter et al. (2002) proposed a Bayesian criterion for comparing complex hierarchical models: deviance information criterion (DIC). Assume a model for observed data $y$ postulates a density $p (y | θ)$ . Then $D (θ) = - 2 log {p (y | θ)}$ is the deviance. The posterior mean deviance $\bar{D} = E [D (θ)]$ was suggested as a measure of fit by Spiegelhalter et al. (1996). The quantity $p_{D} = \bar{D} - D (\bar{θ})$ was derived by Spiegelhalter et al. (2002) as a measure of the effective number of parameters in a model. Then we have $DIC = \bar{D} + p_{D}$ . So DIC, like AIC or BIC, is an aggregate of a measure of fit and a penalty for complexity (Spiegelhalter et al., 2002, 2014). It aims to identify a succinct model that best explains the observed data. For the three Bayesian models, DIC was calculated within each MCMC algorithm. For coxph, DIC is not available.

Geisser and Eddy (1979) proposed a Bayesian predictive approach to model selection. Assume there are several possible models $\{M_{k}, k = 1, \dots, m\}$ for the observed data $y$ , the model $M_{k^{*}}$ that produces the largest likelihood $L_{k}$ is the most appropriate of the models being considered, where $L_{k} = \prod f (y_{r} | y_{(r)})$ Geisser and Eddy, (1979). The product $L_{k}$ is called the pseudo-marginal likelihood by Gelfand et al. (1992) and can be used as a surrogate of the joint marginal likelihood $f (y)$ of the data (Gelfand et al., 1992; Dey et al., 1997). The component $f (y_{r} | y_{(r)})$ is called the conditional predictive ordinate (CPO) by Dey et al. (1997). If $\{Y_{r} : r = 1, \dots, n\}$ are conditionally independent given the set of model parameters $θ$ , then $f (y_{r} | y_{(r)}, θ^{s})$ simplifies to $f (y_{r} | θ^{(s)})$ and the Monte Carlo estimate of ${CPO}_{r}$ becomes: ${CPO}_{r} = {\{\frac{1}{B} \sum_{s = 1}^{B} (1 / f (y_{r} | θ^{(s)}))\}}^{- 1},$ where $f (y_{r} | θ^{(s)})$ is the $r$ th individual likelihood evaluated at iteration $s$ Dey et al. (1997). For the clustered partly interval-censored data, the cross-validation is with respect to individual subject $[i, j]$ rather than cluster $i$ and the individual subjects are conditionally independent given the random effects $ϕ_{i}$ . This is the formula used to calculate the CPOs for the three Bayesian models. We compute the negative log pseudo-marginal likelihood (NLLK) for the three Bayesian models. For coxph, we present the negative of the log-likelihood.

3 Simulation studies

3.1 Simulation I

We evaluate the performance of the proposed method through simulation studies. We fit the mixed effects PH model, the Zhou and Hanson (2018a)’s PH model with a random intercept with the survregbayes function in the spBayesSurv package, the Komáek and Lesaffre (2007)’s mixed effects AFT model with the bayessurvreg1 function in the bayesSurv package, the Komáek et al. (2007)’s mixed effects AFT model with the bayessurvreg2 function in the bayesSurv package, and the traditional Cox PH model (Cox, 1972) for right-censored data with the coxph function in the survival package (Therneau et al., 2021). For coxph, we ignore the clustering and convert the partly interval-censored data to right-censored data as conventionally done by practitioners, that is, take the right endpoints of finite time intervals (left-censored and interval-censored observations) as the observed event times. The purpose of including coxph is to see how the conventional approach can introduce bias in the estimation of fixed effects and survival function.

Failure times are generated from a PH model with a random intercept (Data 1) and a PH model with a random intercept and a random treatment effect (Data 2):

Data 1 : λ (t | x_{i j 1}, x_{i j 2}, ϕ_{i 0}) = λ_{0} (t) exp (β_{1} x_{i j 1} + β_{2} x_{i j 2} + ϕ_{i 0})

Data 2 : λ (t | x_{i j 1}, x_{i j 2}, ϕ_{i 0}, ϕ_{i 1}) = λ_{0} (t) exp (β_{1} x_{i j 1} + β_{2} x_{i j 2} + ϕ_{i 0} + ϕ_{i 1} x_{i j 2})

where $Λ_{0} (t) = log (1 + t)$ , $β_{1} = β_{2} = 1$ , $x_{i j 1}$ ’s $\sim N (0, {0.5}^{2})$ , $x_{i j 2}$ ’s $\sim Bernoulli (0.5)$ , $ϕ_{i 0} \sim N (0, {0.5}^{2})$ , and $ϕ_{i 1} \sim N (0, {0.5}^{2})$ . The number of medical examinations performed for each person equals 1 plus a random number generated from a Poisson distribution with mean 2. The gap times between adjacent examinations are generated from an Exponential distribution with mean 1. The observed time intervals are the ones that contain the true failure times, and then around $20 %$ of which are set to have exact event times observed. It can be deemed as a multi-centre randomized clinical trial with $I$ $=$ 25 centres, each with $n_{i}$ $=$ 20 patients. In each data-set, there are $N$ $=$ 500 subjects. A total of 100 datasets are generated in Data 1 and Data 2 respectively.

In fitting the proposed method, we set the degree $=$ 3 and the knots to be 0, 2, 6, and the maximum of the observed time points plus 1 for the basis I-splines, based on the choice of Pan and Cai (2020). For hyperparameters, we set $σ_{0}^{2} = 100$ , $a_{η} = b_{η} = 1$ , $a_{ζ} = b_{ζ} = 1$ , and $a_{τ} = b_{τ} = 1$ . Fast convergence and good mixing were observed for all key parameters. For each MCMC chain, we set the total number of iterations $=$ 6000, burn-in $=$ 1000, and thin $=$ 1. As sensitivity analysis, we have also tried the knots as 10 equally spaced points from 0 to the maximum of observed time points plus 1, $σ_{0}^{2} = 10, 1000$ , $a_{η} = b_{η} = 0.1, 0.01$ , $a_{ζ} = b_{ζ} = 0.1, 0.01$ , and $a_{τ} = b_{τ} = 0.1, 0.01$ . The results are very similar and are not presented for space concern.

Table 1 presents the simulation results of fitting the proposed method, survregbayes, bayessurvreg1, bayessurvreg2, and coxph to the 100 datasets generated under Data 1. For each parameter, the point estimate is the average of the 100 posterior means, SSD is the sample standard deviation of the 100 posterior means, ESE is the average of the 100 empirical standard errors, 95CP is the coverage probability of the 100 $95 %$ credible intervals, and ESS is the effective sample size computed using the coda package Plummer et al. (2019). The two model selection criteria introduced in Section 2.6 are computed: DIC and NLLK.

Table 1

Simulation Ia: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph under Data 1.

R function	True	Estimate	SSD	ESE	95CP	ESS	DIC	NLLK
Proposed method	$β_{1} = 1$	1.020	0.142	0.132	0.92	230	734	368
	$β_{2} = 1$	0.996	0.126	0.127	0.98	353
	$τ = 4$	3.206	0.822	1.082	0.89	580
survregbayes	$β_{1} = 1$	1.021	0.143	0.133	0.93	419	732	367
	$β_{2} = 1$	1.020	0.134	0.125	0.96	233
	$τ = 4$	3.190	0.803	1.085	0.90	1201
bayessurvreg1	$β_{1} = 1$	−1.352	0.198	0.183		874	754	378
	$β_{2} = 1$	−1.373	0.180	0.178		485
bayessurvreg2	$β_{1} = 1$	−1.352	0.197	0.183		757	789	394
	$β_{2} = 1$	−1.377	0.179	0.177		367
coxph	$β_{1} = 1$	0.618	0.128	0.106	0.11			1957
	$β_{2} = 1$	0.646	0.104	0.106	0.08

As seen in Table 1, for Data 1, the proposed method and survregbayes provide very good estimation for the regression coefficients, with small bias, sample standard deviation close to empirical standard error, and coverage probability close to nominal level. Their overall model goodness-of-fits are very similar and both outperform bayessurvreg1 and bayessurvreg2, as indicated by the model selection criteria: DIC and NLLK. Both also provide reasonably close estimates of the random intercept precision $τ$ . The coxph function performs poorly with large bias, low coverage probability, and large negative log-likelihood.

Table 2 summarizes the simulation results of fitting the proposed method, survregbayes, bayessurvreg1, bayessurvreg2, and coxph to the 100 datasets generated under Data 2. The proposed method performs very well with small bias, sample standard deviation close to empirical standard error, and coverage probability close to nominal level, and the smallest DIC and NLLK. The Zhou and Hanson (2018a)’s method, which is implemented by the survregbayes function, only allows for a random intercept. We still included it here to see how ignoring the random treatment effect might affect estimation and inference. The point estimate of $β_{2}$ from survregbayes turns out to be robust, which is understandable as the mean of the random treatment effect is zero. However, the estimated standard error of $β_{2}$ is reduced (0.126 vs 0.158), which leads to a lower coverage probability (0.88) and could also lead to wrongly significant inference. This finding is consistent with Yavuz and Lambert (2016) who conclude that ‘the main consequence of ignoring the frailty is a reduction in the standard errors’. The result from Komáek and Lesaffre (2007)’s method (implemented by the bayessurvreg1 function) is interesting. Most importantly, as shown in Figure 1, the MCMC chain of $β_{2}$ does not converge. Of course, the ESS are super small. The result from Komáek et al. (2007)’s method (bayessurvreg2) shares the similar issue and furthermore, the MCMC chain of $β_{1}$ does not converge either (traceplot not shown). We will further investigate in Simulation II to see whether the non-convergence of $β_{2}$ (even $β_{1}$ ) is caused by model mis-specification (fitting an AFT model to data generated from a PH model), or it is an internal shortcoming of the method itself. In addition, under Data 2, the DIC and NLLK values from bayessurvreg1 and bayessurvreg2 are significantly larger than these from the proposed method and survregbayes due to model mis-specification.

Table 2

Simulation Ib: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, bayessurvreg1, bayessurvreg2 and coxph under Data 2.

$R$ function	True	Estimate	SSD	ESE	95CP	ESS	DIC	NLLK
Proposed method	$β_{1} = 1$	1.065	0.140	0.137	0.92	194	730	367
	$β_{2} = 1$	0.990	0.149	0.190	0.97	127
survregbayes	$β_{1} = 1$	1.010	0.135	0.132	0.94	415	738	370
	$β_{2} = 1$	0.991	0.158	0.126	0.88	238
bayessurvreg1	$β_{1} = 1$	−1.363	0.211	0.185		836	859	436
	$β_{2} = 1$	−1.313	0.236	0.174		7
bayessurvreg2	$β_{1} = 1$	−1.375	0.202	0.182		53	865	441
	$β_{2} = 1$	−1.366	0.224	0.182		10
coxph	$β_{1} = 1$	0.642	0.116	0.107	0.11			1943
	$β_{2} = 1$	0.590	0.134	0.106	0.06

Figure 1

Traceplot of MCMC chain of $β_{2}$ – Left panel: Proposed method; middle panel: bayessurvreg1; right panel: bayessurvreg2.

Figure 2 plots the estimated baseline survival curves versus the true under Data 1 and Data 2 based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2, and coxph. The proposed method and survregbayes very closely approximate the baseline survival curve. The estimation from bayessurvreg1 and bayessurvreg2 diverges off starting from $T = 2$ , due to the different underlying distribution assumption. The estimated baseline survival curve by the traditional Cox PH model is significantly different from the true curve.

Figure 2

Plot of ${\hat{S}}_{0} (t)$ – Left panel: Data 1 with a normal random intercept; right panel: Data 2 with a normal random intercept and a random treatment effect.

3.2 Simulation II

To give an advantage to the mixed effects AFT model by Komáek and Lesaffre (2007) and Komáek et al. (2007) (implemented by the R functions bayessurvreg1 and bayessurvreg2), we conducted a second simulation study by generating datasets under mixed effects AFT models. Specifically, 100 datasets are generated from an AFT model with a random intercept (Data 3) and another 100 datasets are generated from an AFT model with a random intercept and a random treatment effect (Data 4):

Data 3 : log (T_{i j}) = β_{0} + β_{1} x_{i j 1} + β_{2} x_{i j 2} + ϕ_{i 0} + σ ϵ_{i j}

Data 4 : log (T_{i j}) = β_{0} + β_{1} x_{i j 1} + β_{2} x_{i j 2} + ϕ_{i 0} + ϕ_{i 1} x_{i j 2} + σ ϵ_{i j}

where we set $σ = 1$ and $ϵ_{i j} \sim N (0, 1)$ . So both data models are log-normal AFT model with random effects. We can derive the baseline survival function as $S_{0} (t) = 1 - Φ (\frac{log (t)}{σ})$ , where $Φ (\cdot)$ is the cumulative distribution function of standard normal distribution. Other settings are exactly the same as in Simulation I.

We compare the proposed method, Zhou and Hanson (2018a)’s method (implemented by survregbayes), Komáek and Lesaffre (2007)’s method (implemented by bayessurvreg1), Komáek et al. (2007)’s method (implemented by bayessurvreg2), and the conventional Cox PH model (implemented by coxph). Similar to Simulation I, Zhou and Hanson (2018a)’s method only allows for a random intercept but is fitted to both the 100 datasets generated under Data 3 and the 100 datasets generated under Data 4. Prior specifications are exactly the same as in Simulation I.

As we can see in Table 3, bayessurvreg1 and bayessurvreg2 provide the best fit under Data 3, as would be expected because the datasets are generated under an AFT model. However, the proposed method and survregbayes perform almost as well, as indicated by the very close model fitting criteria values (DIC and NLLK) to those of bayessurvreg1 and bayessurvreg2. The coxph function performs poorly with large negative log-likelihood.

Table 3

Simulation IIa: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph under Data 3.

$R$ function	True	Estimate	SSD	ESE	95CP	ESS	DIC	NLLK
Proposed method	$β_{1} = 1$	−1.138	0.156	0.142		321	870	437
	$β_{2} = 1$	−1.101	0.146	0.143		793
survregbayes	$β_{1} = 1$	−1.119	0.153	0.142		472	867	435
	$β_{2} = 1$	−1.067	0.143	0.134		327
bayessurvreg1	$β_{1} = 1$	1.023	0.131	0.124	0.93	1330	858	429
	$β_{2} = 1$	1.001	0.123	0.119	0.94	710
bayessurvreg2	$β_{1} = 1$	1.024	0.130	0.122	0.94	824	862	433
	$β_{2} = 1$	1.003	0.120	0.118	0.97	406
coxph	$β_{1} = 1$	−0.809	0.151	0.121				1461
	$β_{2} = 1$	−0.826	0.115	0.123

The results under Data 4 are presented in Table 4. Similar to what we have seen under Data 2, by ignoring the random treatment effect, the estimated standard error of treatment effect ( $β_{2}$ ) from survregbayes is reduced (0.135 vs 0.183). The MCMC chains of $β_{2}$ from bayessurvreg1 and bayessurvreg2 still do not converge (Figure 3). Also the MCMC chain of $β_{1}$ from bayessurvreg2 fails to converge, for which the traceplot is not shown but can still be revealed from the small ESS. So now we can conclude that the inclusion of a random effect of a covariate would adversely impact the convergence of the MCMC chain of at least the same covariate’s fixed effect. This is an internal shortcoming of the methods of Komáek and Lesaffre (2007) and Komáek et al. (2007). The proposed method can provide better sample quality and efficiency. The proposed method turns out to perform better even for the data generated from a mixed effects AFT model. It is indicated by the smallest DIC and NLLK values. This is not expected but explainable. The reason, as we have just described, is because: (1) The MCMC chains of $β_{2}$ from the methods of Komáek and Lesaffre (2007) and Komáek et al. (2007) do not converge. (2) It seems that the semiparametric distribution of the proposed method is able to approximate the AFT distribution reasonably well.

Table 4

Simulation IIb: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph under Data 4.

$R$ function	True	Estimate	SSD	ESE	95CP	ESS	DIC	NLLK
Proposed method	$β_{1} = 1$	−1.146	0.151	0.149		277	858	433
	$β_{2} = 1$	−1.133	0.204	0.219		243
survregbayes	$β_{1} = 1$	−1.072	0.144	0.142		467	866	435
	$β_{2} = 1$	−1.029	0.183	0.135		321
bayessurvreg1	$β_{1} = 1$	0.991	0.107	0.126	0.97	1215	991	508
	$β_{2} = 1$	0.996	0.151	0.130	0.89	11
bayessurvreg2	$β_{1} = 1$	0.996	0.110	0.124	0.98	47	992	510
	$β_{2} = 1$	1.006	0.154	0.135	0.88	16
coxph	$β_{1} = 1$	−0.765	0.128	0.121				1467
	$β_{2} = 1$	−0.818	0.166	0.124

Figure 3

Traceplot of MCMC chain of $β_{2}$ – Left panel: Proposed method; middle panel: bayessurvreg1; right panel: bayessurvreg2.

Figure 4 plots the estimated baseline survival curves averaged over the 100 datasets generated under Data 3 and Data 4 respectively, based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2, and coxph, versus the true. For Data 3 and Data 4, the proposed method, survregbayes, bayessurvreg1, and bayessurvreg2 all provide very close estimation of the baseline survival curve. We recall that for datasets generated from PH models (Data 1 and Data 2), the estimated baseline survival curves from bayessurvreg1 and bayessurvreg2 diverge off starting from $T = 2$ . The fact that our model can still closely estimate the baseline survival curve of an AFT model shows its flexibility. The estimation from the traditional Cox PH model still deviates significantly from the true curve.

Figure 4

Plot of ${\hat{S}}_{0} (t)$ – Left panel: Data 3 with a normal random intercept; right panel: Data 4 with a normal random intercept and a random treatment effect.

4 Metastatic colorectal cancer data

The proposed method can be implemented with the main function PICBayes in the PICBayes package (Pan, 2021). The dataset is called mCRC in the package. An R script to reproduce the analysis presented here is in the supplementary material of this article. The mCRC dataset was derived by the first author from the raw datasets of a randomized phase III study that compares the efficacy of panitumumab $+$ FOLFIRI versus FOLFIRI in patients with previously treated metastatic colorectal cancer. FOLFIRI is a combination of chemotherapy drugs: fluorouracil, leucovorin, and irinotecan. Panitumumab is a fully human monoclonal antibody directly against the epidermal growth factor receptor (EGFR). The efficacy endpoint is progression-free survival (PFS). Two binary covariates are of interest: treatment arm (FOLFIRI vs panitumumab $+$ FOLFIRI) and patient tumour KRAS mutation status (wild-type vs mutant). KRAS stands for the gene Kirsten rat sarcoma viral oncogene homolog. The dataset contains $I = 185$ sites with number of patients within sites ranging from 1 to 23. There are a total of $N = 855$ patients with 52 deaths on-study, 168 left-censored, 329 interval-censored, and 336 right-censored.

We analyze the data with a random intercept using the proposed method, survregbayes, bayessurvreg1, bayessurvreg2, and coxph. Since retrospective analysis of phase II and III studies have already demonstrated that mCRC patients with tumour KRAS mutation do not derive clinical benefit from panitumumab (Peeters et al., 2010, 2014), we expect there would be a significant treatment by mutation interaction. Hence, we include this interaction term in each of the five models. For the proposed method, we set degree = 3. Given that the maximum time point is 974 days, which is large, we set the knots to be 79 equally spaced points from 0 to 975. The results are summarized in Table 5. The proposed method and survregbayes have similar DIC and NLLK values, which indicates that they fit the data comparably well. However, only the proposed method and bayessurvreg2 were able to detect a significant interaction. We should also note that the Zhou and Hanson (2018a)’s method is not applicable if one wants to account for potential random variation of treatment effect across sites. If random treatment effect does exist and the Zhou and Hanson (2018a)’s method is fitted, then it may give false positive inference about the treatment effect. Then panitumumab would be added to tumour KRAS mutant patients, which would not stop cancer cells from growing but instead cause severe toxicities. So under some situations, such as clinical trials, it is very important that we account for random treatment effect if it exists, so as not to make wrong even life-threatening medical decisions. The bayessurvreg1 and bayessurvreg2 functions have larger DIC and NLLK values, which indicates that they do not fit the data as well as the two Bayesian methods under the PH model. The $95 %$ credible interval for treatment effect among wild-type patients is $(- 0.687, - 0.215)$ . The $95 %$ credible interval for treatment effect among mutant patients is $(- 0.334, 0.199)$ , which contains zero. We conclude that adding panitumumab monoclonal antibody to the FOLFIRI regimen can significantly increase PFS for KRAS wild-type patients. However, it would not be beneficial for KRAS mutant patients. This is consistent with current common clinical understanding.

Table 5

Metastatic colorectal cancer trial ( $N = 855$ ): Estimation result based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph.

$R$ function		Estimate	SE	95% CI	ESS	DIC	NLLK
Proposed method	treatment	−0.453	0.120	(−0.687, −0.215)	572	3089	1552
	mutation	−0.019	0.125	(−0.256, 0.234)	457
	interaction	0.388	0.172	(0.044, 0.717)	581
survregbayes	treatment	−0.387	0.132	(−0.637, −0.143)	370	3088	1547
	mutation	0.023	0.131	(−0.242, 0.284)	466
	interaction	0.335	0.191	(−0.050, 0.712)	398
bayessurvreg1	treatment	0.417	0.127	(0.157, 0.660)	278	3154	1578
	mutation	0.023	0.148	(−0.287, 0.302)	193
	interaction	−0.352	0.195	(−0.721, 0.048)	405
bayessurvreg2	treatment	0.445	0.109	(0.230, 0.664)	374	3153	1581
	mutation	0.027	0.129	(−0.225, 0.281)	270
	interaction	−0.381	0.173	(−0.720, −0.036)	557
coxph	treatment	−0.340	0.116	(−0.567, −0.112)			3229
	mutation	0.025	0.123	(−0.215, 0.265)
	interaction	0.274	0.172	(−0.063, 0.610)

Figure 5 presents the estimated survival curves for the four subgroups formed by treatment arm and mutation status, based on the proposed method. The survival probability is obviously the highest for wild-type patients receiving panitumumab $+$ FOLFIRI, indicating the clinical effectiveness of panitumumab for KRAS wild-type patients. The other three curves (mutant patients receiving panitumumab $+$ FOLFIRI, wild-type patients receiving FOLFIRI, and mutant patients receiving FOLFIRI) are pretty much close to each other, indicating the futility of panitumumab for KRAS mutant patients and the equivalent efficacy of FOLFIRI to either KRAS wild-type patients or KRAS mutant patients.

Figure 5

Metastatic colorectal cancer trial ( $N = 855$ ): Estimated survival curves for four subgroups based on the proposed method.

5 Conclusions and discussions

There are currently only three methods in the literature for analyzing clustered partly interval-censored data: Komáek and Lesaffre (2007), Komáek et al. (2007), and Zhou and Hanson (2018a). We developed a Bayesian approach under the PH model which is flexible by its semiparametric nature. The proposed method outperforms the Zhou and Hanson (2018a)’s method by allowing for not only a random intercept but also random effects of covariates. As we have shown in Simulation I and Simulation II, if the random effect of a covariate exists, fitting the Zhou and Hanson (2018a)’s method would lead to a reduced estimated standard error for the covariate’s fixed effect. That could give a false positive inferential conclusion. The Komáek and Lesaffre (2007)’s method and the Komáek et al. (2007)’s method include both a random intercept and random effects of covariates. However, their methods may not quite work when the random effect of a covariate is included in the model, as we have shown in both Simulation I and Simulation II. For the regression coefficient $β_{p}$ with a corresponding random effect, its MCMC chain does not converge anymore. In other words, their method may not give a good estimation for the regression coefficient of a covariate if one also includes a random effect for this covariate in the model. This finding gives us reasons to believe that adding random effects to a model with a random intercept only is not a trivial task from both method development and MCMC computation perspectives, as we would have assumed. This may also explain why most of the models for clustered survival data only include a random intercept/shared frailty, as we have seen in the literature review. This also further confirms the meaningfulness of the proposed method that includes random effects of covariates, given the existence of the Zhou and Hanson (2018a)’s method which fits a PH model with a random intercept only.

As we have seen in the literature review, more than half of the current publications assume a gamma frailty, for the benefit of a conjugate gamma posterior distribution. We also used gamma frailties for clustered interval-censored data (Pan et al., 2015). However, due to increased complexity of data structure, we have found that gamma frailties do not work well for our specific data structure and model here. The algorithm could get stuck due to some very large values generated from gamma posterior distributions for frailties. It makes sense if we think of the skewed nature of a gamma distribution. While with a normal distribution which is symmetric and hence less likely to give extreme values, we did not have this problem. Curious about a Bayesian nonparametric prior, we have also tried a Dirichlet process normal mixture prior (Görür and Rasmussen, 2010) for each random effect. Specifically, we assume each random effect to follow a normal distribution with mean zero and precision parameter following a Dirichlet process. The results are very similar to our current method and the algorithm is significantly more complicated due to the involvement of the blocked Gibbs sampler (Ishwaran and Zarepour, 2000; Ishwaran and James, 2001). After a certain number of MCMC iterations, all clusters converge to sharing the same precision parameter. So the Dirichlet process normal mixture prior reduces to a normal prior. We have also explored a multivariate normal prior for the vector of random effects. It turns out the estimation for fixed effects are not as good as the current method, which assumes a normal prior for each random intercept/random effect separately. We think the reason is that more constraints are actually being imposed if we force the joint distribution of the random effects to be of a specific form. Even though it has been a common practice for linear models, it may not be as good a choice for clustered partly interval-censored data. Komáek and Lesaffre (2007) and Komáek et al. (2007) assumed a multivariate normal prior for the random effects, and we have seen that the MCMC chain for the fixed effect with an associated random effect does not converge. Of course, it is unknown whether assuming a normal prior for each random intercept/random treatment effect separately under their model framework would solve this non-convergence problem. Yavuz and Lambert (2016) employed a flexible P-splines specification to approximate the distribution of a frailty, this nonparametric form could be something that is worthy of exploring in the future.

Appendix

After initializing values for the parameters, the proposed MCMC algorithm proceeds in the following steps.

Let $Z_{i j} = 0$ and $W_{i j} = 0$ for all $i$ and $j$ , $Z_{i j l} = 0$ and $W_{i j l} = 0$ for all $i$ , $j$ , and $l$ . If $δ_{i j 1} = 1$ , then sample

Z_{i j} \sim Poi (Λ_{0} (R_{i j}) e^{β' x_{i j} + ϕ_{i}' ω_{i j}}) 1 (Z_{i j} > 0),

(Z_{i j 1}, \dots, Z_{i j K}) \sim Multinomial (Z_{i j}; γ_{1} I_{1} (R_{i j}), \dots, γ_{K} I_{K} (R_{i j})) .

If $δ_{i j 2} = 1$ , then sample

W_{i j} \sim Poi (\{Λ_{0} (R_{i j}) - Λ_{0} (L_{i j})\} e^{β' x_{i j} + ϕ_{i}' ω_{i j}}) 1 (W_{i j} > 0),

(W_{i j 1}, \dots, W_{i j K}) \sim Multinomial (W_{i j}; γ_{1} \{I_{1} (R_{i j}) - I_{1} (L_{i j})\}, \dots, γ_{K} \{I_{K} (R_{i j}) - I_{K} (L_{i j})\}) .

Sample $(u_{i j 1}, \dots, u_{i j K}) \sim Multinomial (1; γ_{1} M_{1} (t_{i j}), \dots, γ_{K} M_{K} (t_{i j}))$ .

For $β_{p}$ corresponding to a numeric covariate, use the Metropolis-Hastings algorithm to sample from its full conditional distribution

p (β_{p} | \cdot) \propto exp [\sum_{i = 1}^{I} \sum_{j = 1}^{n_{i 1}} \{x_{i j p} β_{p} - Λ_{0} (t_{i j}) e^{β' x_{i j} + ϕ_{i}' ω_{i j}}\}

+ \sum_{i = 1}^{I} \sum_{j = n_{i 1} + 1}^{n_{i}} \{x_{i j p} β_{p} (Z_{i j} δ_{i j 1} + W_{i j} δ_{i j 2}) - e^{β' x_{i j} + ϕ_{i}' ω_{i j}} (Λ_{0} (R_{i j}) (δ_{i j 1} + δ_{i j 2}) + Λ_{0} (L_{i j}) δ_{i j 3})\}] e^{- \frac{β_{p}^{2}}{2 σ_{0}^{2}}} .

For $β_{p}$ corresponding to a dummy variable, let $ζ_{p} = exp (β_{p})$ , sample $ζ_{p}$ from

Ga (a_{ζ} + \sum_{i = 1}^{I} \sum_{j = 1}^{n_{i 1}} x_{i j p} + \sum_{i = 1}^{I} \sum_{j = n_{i 1} + 1}^{n_{i}} x_{i j p} (Z_{i j} δ_{i j 1} + W_{i j} δ_{i j 2}),

b_{ζ} + \sum_{i = 1}^{I} \sum_{j = 1}^{n_{i 1}} Λ_{0} (t_{i j}) e^{β_{- p}' x_{i j, - p} + ϕ_{i}' ω_{i j}} x_{i j p}

+ \sum_{i = 1}^{I} \sum_{j = n_{i 1} + 1}^{n_{i}} e^{β_{- p}' x_{i j, - p} + ϕ_{i}' ω_{i j}} \{Λ_{0} (R_{i j}) (δ_{i j 1} + δ_{i j 2}) + Λ_{0} (L_{i j}) δ_{i j 3}\} x_{i j p}),

where $β_{- p} = \{β_{k} : k \neq p\}$ and $x_{i j, - p} = \{x_{i j k} : k \neq p\}$ .

Sample $γ_{l}$ , $l = 1, \dots, K$ , from

Ga (1 + \sum_{i = 1}^{I} \sum_{j = 1}^{n_{i 1}} u_{i j l} + \sum_{i = 1}^{I} \sum_{j = n_{i 1} + 1}^{n_{i}} (Z_{i j l} δ_{i j 1} + W_{i j l} δ_{i j 2}),

η + \sum_{i = 1}^{I} \sum_{j = 1}^{n_{i 1}} e^{β' x_{i j} + ϕ_{i}' ω_{i j}} I_{l} (t_{i j}) + \sum_{i = 1}^{I} \sum_{j = n_{i 1} + 1}^{n_{i}} e^{β' x_{i j} + ϕ_{i}' ω_{i j}} \{I_{l} (R_{i j}) (δ_{i j 1} + δ_{i j 2}) + I_{l} (L_{i j}) δ_{i j 3}\}) .

Sample $η$ from $Ga (a_{η} + K, b_{η} + \sum_{l = 1}^{K} γ_{l})$ .

For a specific $q$ , sample $ϕ_{i q}$ , $i = 1, \dots, I$ , using the Metropolis-Hastings algorithm from its full conditional distribution

p (ϕ_{i q} | \cdot) \propto exp [\sum_{j = 1}^{n_{i 1}} (ϕ_{i}^{'} z_{i j}) - \sum_{j = 1}^{n_{i 1}} Λ_{0} (t_{i j}) e^{β' x_{i j} + ϕ_{i}' ω_{i j}}]

exp [\sum_{j = n_{i 1} + 1}^{n_{i}} (ϕ_{i}^{'} ω_{i j}) (Z_{i j} δ_{i j 1} + W_{i j} δ_{i j 2}) - \sum_{j = n_{i 1} + 1}^{n_{i}} e^{β' x_{i j} + ϕ_{i}' ω_{i j}} \{Λ_{0} (R_{i j}) (δ_{i j 1} + δ_{i j 2}) + Λ_{0} (L_{i j}) δ_{i j 3}\}]

exp (- 0.5 τ ϕ_{i q}^{2}) .

For a specific $q$ , sample $τ_{q}$ from $Ga (a_{τ} + 0.5 I, b_{τ} + \frac{1}{2} ϕ_{q}' ϕ_{q})$ .

Supplementary Materials

Supplementary files are available online.

Supplemental Material for Bayesian semiparametric mixed effects proportional hazards model for clustered partly interval-censored data by Chun Pan and Bo Cai, in Statistical Modelling

Footnotes

Acknowledgements

The authors thank the reviewers and the associate editor for their insightful comments which greatly improved this article.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

CP was funded by the National Institute of General Medical Sciences of the National Institutes of Health under Award Number SC2GM135078.

References

Aslanidou

, Dey

and Sinha

(1998) Bayesian analysis of multivariate survival data using Monte Carlo methods. The Canadian Journal of Statistics , 26, 33–48.

Bellamy

, Li

, Ryan

, Lipsitz

, Canner

and Wright

(2004) Analysis of clustered and interval censored data from a community-based study in asthma. Statistics in Medicine , 23, 3607–3621

Chen

, Shen

and Lee

(2020) A gamma-frailty model for interval-censored data with dependent examination times: a computationally efficient approach. Communications in Statistics-Simulation and Computation . doi: 10.1080/03610918.2020.1790600.

Cox

(1972) Regression models and life-tables (with discussion). Journal of the Royal Statistical Society: Series B , 34, 187–220.

Dey

, Chen

and Chang

(1997) Bayesian approach for nonlinear random effects models. Biometrics , 53, 1239–1252.

Gao

, Zeng

and Lin

(2017) Semiparametric estimation of the accelerated failure time model with partly interval-censored data. Biometrics , 73, 1161–1168.

Geisser

and Eddy

(1979) A predictive approach to model selection. Journal of the American Statistical Association , 74, 153–160.

Gelfand

(1992) Model determination using predictive distributions with implementation via sampling-based methods (with discussion). In Bayesian Statistics 4 , edited by Bernardo

, Berger

, Dawid

and Smith

AFM

, pages 147–167. Oxford: Oxford University Press.

Glidden

and Vittinghoff

(2004) Modelling clustered survival data from multicentre clinical trials. Statistics in Medicine , 23, 369–388.

10.

Goethals

, Ampe

, Berkvens

, Laevens

, Janssen

and Duchateau

(2009) Modeling interval-censored, clustered cow udder quarter infection times through the shared gamma frailty model. Journal of Agricultural, Biological, and Environmental Statistics , 14, 1–14.

11.

Görür

and Rasmussen

(2010) Dirichlet process Gaussian mixture models: Choice of the base distribution. Journal of Computer Science and Technology , 25, 615–626.

12.

Hastings

(1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika , 57, 97–109.

13.

Huang

and Wellner

(1997) Interval censored survival data: A review of recent progress. In Proceedings of the First Seattle Symposium in Biostatistics , edited by Lin

and Fleming

, pages 123–169. New York: Springer.

14.

Huang

(1999) Asymptotic properties of nonparametric estimation based on partly interval-censored data. Statistica Sinica , 9, 501–519.

15.

Ishwaran

and Zarepour

(2000) Markov chain Monte Carlo in approximate Dirichlet and beta two-parameter process hierarchical models. Biometrika , 87, 371–390.

16.

Ishwaran

and James

(2001) Gibbs sampling methods for stick-breaking priors. Journal of the American Statistical Association , 96, 161–173.

17.

Kim

(2003) Maximum likelihood estimation for the proportional hazards model with partly interval-censored data. Journal of the Royal Statistical Society: Series B , 65, 489–502.

18.

Klein

(1992) Semiparametric estimation of random effects using the Cox model based on the EM algorithm. Biometrics , 48, 795–806.

19.

Komárek

and Lesaffre

(2007) Bayesian accelerated failure time model for correlated interval-censored data with a normal mixture as an error distribution. Statistica Sinica , 17, 549–569.

20.

Komárek

, Lesaffre

and Legrand

(2007) Baseline and treatment effect heterogeneity for survival times between centers using a random effects accelerated failure time model with flexible error distribution. Statistics in Medicine , 26, 5457–5472.

21.

Komárek

(2018) bayesSurv: Bayesian survival regression with flexible error and random effects distributions. https://cran.r-project.org/package=bayesSurv. R package version 3.2.

22.

Pan

, Cai

and Wang

(2015) Multiple frailty model for clustered interval-censored data with frailty selection. Statistical Methods in Medical Research , 26, 1308–1322.

23.

Pan

, Cai

and Wang

(2020) A Bayesian approach for analyzing partly interval-censored data under the proportional hazards model. Statistical Methods in Medical Research . doi: 10.1177/0962280220921552.

24.

Pan

and Cai

(2020) A Bayesian model for spatial partly interval-censored data. Communication in Statistics-Simulation and Computation . doi: 10.1080/03610918.2020.1839497.

25.

Pan

(2021) PICBayes: Bayesian models for partly interval-censored data. https://cran.r-project.org/package=PICBayes. R package version 1.0.

26.

Peeters

, Price

, Cervantes

, Sobrero

, Ducreux

, André

, Chan

, Lordick

, Punt

CJA

, Strickland

, Wilson

, Ciuleanu

, Roman

, Cutsem

, Tzekova

, Collins

, Oliner

, Rong

and Gansert

(2010) Randomized phase III study of panitumumab with fluorouracil, leucovorin, and irinotecan (FOLFIRI) compared with FOLFIRI alone as second-line treatment in patients with metastatic colorectal cancer. Journal of Clinical Oncology , 28, 4706–4713.

27.

Peeters

, Price

, Cervantes

, Sobrero

, Ducreux

, Hotko

, André

, Chan

, Lordick

, Punt

CJA

, Strickland

, Wilson

, Ciuleanu

, Roman

, Tian

and Sidhu

(2014) Final results from a randomized phase 3 study of FOLFIRI ± panitumumab for second-line treatment of metastatic colorectal cancer. Annals of Oncology , 25, 107–116.

28.

Plummer

, Best

, Cowles

, Vines

, Sarkar

, Bates

, Almond

and Magnusson

(2019) coda: Output analysis and diagnostics for MCMC . https://cran.r-project.org/package=coda. R package version 0.19-3.

29.

Ramsay

(1988) Monotone regression splines in action. Statistical Science , 3, 425–441.

30.

Ripatti

and Palmgren

(2000) Estimation of multivariate frailty models using penalized partial likelihood. Biometrics , 56, 1016–1022.

31.

Sargent

(1998) A general framework for random effects survival analysis in the Cox proportional hazards setting. Biometrics , 54, 1486–1497.

32.

Spiegelhalter

, Thomas

, Best

and Gilks

(1996) BUGS Examples Volume 1, Version 0.5 (Version ii) . Cambridge: Medical Research Council Biostatistics Unit.

33.

Spiegelhalter

, Best

, Carlin

and van der Linde

(2002) Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B , 64, 583–639.

34.

Spiegelhalter

, Best

, Carlin

and van der Linde

(2014) The deviance information criterion: 12 years on. Journal of the Royal Statistical Society: Series B , 76, 458–493.

35.

Therneau

, Lumley

, Atkinson

and Crowson

(2021) survival: Survival analysis . https://cran.r-project.org/package=survival. R package version 3.2-13.

36.

Yavuz

and Lambert

(2016) Semi-parametric frailty model for clustered interval-censored data. Statistical Modelling , 16, 360–391.

37.

Zhao

, Zhao

, Sun

and Kim

JS.

(2008) Generalized log-rank tests for partly interval-censored failure time data. Biometrical Journal , 50, 375–385.

38.

Zhou

and Hanson

(2018a) A unified framework for fitting Bayesian semiparametric models to arbitrarily censored survival data, including spatially-referenced data. Journal of the American Statistical Association , 113, 571–581.

39.

Zhou

and Hanson

(2018b) spBayesSurv: Bayesian modeling and analysis of spatially correlated survival data . https://cran.r-project.org/package=spBayesSurv. R package version 1.1.3.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

A Bayesian semiparametric mixed effects proportional hazards model for clustered partly interval-censored data

Abstract

Keywords

1 Introduction

2 The model

2.1 Data structure

2.2 Model and likelihood

2.6 Model comparison criteria

3 Simulation studies

3.1 Simulation I

Table 1

Simulation Ia: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph under Data 1.

Simulation Ib: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, bayessurvreg1, bayessurvreg2 and coxph under Data 2.

Traceplot of MCMC chain of β 2 – Left panel: Proposed method; middle panel: bayessurvreg1; right panel: bayessurvreg2.

Plot of S ^ 0 t – Left panel: Data 1 with a normal random intercept; right panel: Data 2 with a normal random intercept and a random treatment effect.

Table 3

Simulation IIa: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph under Data 3.

Simulation IIb: Estimation of regression coefficients, ESS, NLLK, and DIC based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph under Data 4.

Traceplot of MCMC chain of β 2 – Left panel: Proposed method; middle panel: bayessurvreg1; right panel: bayessurvreg2.

Plot of S ^ 0 t – Left panel: Data 3 with a normal random intercept; right panel: Data 4 with a normal random intercept and a random treatment effect.

Table 5

Metastatic colorectal cancer trial ( N = 855 ): Estimation result based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph.

Metastatic colorectal cancer trial ( N = 855 ): Estimated survival curves for four subgroups based on the proposed method.

Appendix

Supplementary Materials

Supplementary files are available online.

Footnotes

Acknowledgements

Declaration of Conflicting Interests

Funding

References

Supplementary Material

Traceplot of MCMC chain of $β_{2}$ – Left panel: Proposed method; middle panel: bayessurvreg1; right panel: bayessurvreg2.

Plot of ${\hat{S}}_{0} (t)$ – Left panel: Data 1 with a normal random intercept; right panel: Data 2 with a normal random intercept and a random treatment effect.

Traceplot of MCMC chain of $β_{2}$ – Left panel: Proposed method; middle panel: bayessurvreg1; right panel: bayessurvreg2.

Plot of ${\hat{S}}_{0} (t)$ – Left panel: Data 3 with a normal random intercept; right panel: Data 4 with a normal random intercept and a random treatment effect.

Metastatic colorectal cancer trial ( $N = 855$ ): Estimation result based on the proposed method, survregbayes, bayessurvreg1, bayessurvreg2 and coxph.

Metastatic colorectal cancer trial ( $N = 855$ ): Estimated survival curves for four subgroups based on the proposed method.