Bayesian Autoregressive Frailty Models for Inference in Recurrent Events

1.1 Bayesian nonparametric framework

This paper lies within the gap times approach and develops a Bayesian semiparametric model for gap times between events. We assume that the joint distribution of the finite sequence of gap times for each individual is the product of the conditional distributions of each gap time, given the previous ones. Moreover, we specify a regression model for each of these conditional distributions to link the realization of each gap time to possibly time-varying covariates and previous gap times. To account for inter-subject variability, we introduce individual specific frailty parameters which we model flexibly using a Dirichlet process mixture prior as random effects distribution. Dirichlet process mixture (DPM) models [20, 21] are arguably the most common nonparametric Bayesian prior and have proved successful in many applications due to their flexibility and ease of computation. DPM models are mixtures of a parametric distribution where the mixing measure is the Dirichlet process (DP) introduced by Ferguson [22]. It is well known that the DP is almost surely discrete, and that if G is a DP(M, G₀) with total mass parameter M and baseline distribution G₀, then G admits the following representation [23]:

where δ_θ is a point-mass at θ, the weights follow a stick-breaking process, wh=Vh∏j<h(1−Vj), with Vh∼iidBeta(1,M), and the atoms {θh}h≥1 are such that θh∼iidG0. As the discreteness of G is inappropriate in many applications, it is common to convolve a parametric kernel k(y∣θ) with respect to G, obtaining a DPM:

Due to the discreteness of G, the DPM prior induces clustering of the subjects in the sample based on the trajectory of the recurrent events over time, where the number K of clusters is unknown and learned from the data. Earliest examples of the use of a Bayesian nonparametric distribution for the random effects can be found in Müller and Rosner [24] and Kleinman and Ibrahim [25]. More recently, Pennell and Dunson [26] employ a Dirichlet process prior to build semiparametric dynamic frailty models for multiple event time data, allowing also the frailty parameter to change over time. Moreover, Paulon et al. [27] develop a joint model for gap times between two consecutive hospitalizations and survival time using the DP as random effect distribution, but do not explicitly account for the temporal dependence between gap times observed on the same individual.

In this work, we investigate different strategies to link gap times at time t with previous gap times. We start by assuming a standard Markov model where also the order of dependence p is unknown and is object of inference. We explore two different strategies to specify a prior distribution on p: one involves eliciting a prior directly on the space of all possible Markov models for p ∈ {0, 1, …, P}, while the other approach employs spike and slab base measures and it is in the spirit of stochastic search variable selection [28]. In conclusion, the main contribution of this work is to provide, within a Bayesian nonparametric framework, a model for gap times which accommodates for temporal dependence, performs variable selection on the order of dependence with two different approaches, allows for clustering of individuals and inter-subject variability, is able to handle different numbers of gap times per individual and missingness.

In Section 2 we introduce the model, while in Section 3 we explain how to perform inference on the order of dependence in the Markov structure. In Section 4 we investigate the performance of the proposed approach in simulations and compare the different strategies to model time dependency and to select the order p. Sections 5 and 6 present two medical applications involving recurrent hospitalizations and urinary tract infections, respectively. We conclude the paper in Section 7.

2.1 Nonparametric AR(1)-type models

Following a similar modelling strategy to the one described in Di Lucca et al. [29], a straightforward way to introduce dependence among random effects at different times is to allow the distribution of α_ij to depend on some summary of the observations up to time j – 1:

When j = 1, the distribution of the random effect α_i1 simplifies as the autoregressive term in (2) disappears and it reduces to the Normal distribution with mean m_i0.

We assume conditional independence among subjects, given the parameters, and that (mi0,mi1) are independent under the base measure G₀, which becomes the product of a Normal density for m_i0 and a Normal density for the autoregressive coefficient m_i1. The prior specification is completed as follows:

assuming a priori independence among the different parameters.

When implementing the model in JAGS, we opt for marginal uniform priors on bounded interval for σ² and τ² as suggested by. Gelman [30] as inverse-gamma distributions give numerical instability problems when hyperparameters are small. Moreover, Gelman [30] recommends uniform distributions (on bounded intervals) for scale parameters in hierarchical models as being weakly-informative, opposite to inverse-gamma priors for variance parameters, which are very sensitive to the choice of the hyperparameters themselves.

The choice of g in (2) is obviously crucial and depends on the context and the goals of the inference problem. Common choices in the literature are:

1. g(Yi1,…,Yij−1)=Yij−1, i. e. the random effect at time j depends only on the observation at time j – 1, conditionally on the remaining parameters;

2. g(Yi1,…,Yij−1)=(Yi1+⋯+Yij−1)/(j−1), i. e. the conditional expected value of each α_ij depends on the sample mean of the observations up to time j – 1;

3. g(Yi1,…,Yij−1)=Yi1×⋯×Yij−11/(j−1), this is equivalent to the geometric mean of the observations up to time j – 1.

Note that, when g(Yi1,…,Yij−1)=Yij−1, then (2)–(3) imply that the random effects distribution at time j is a DPM of AR(1) processes, with dependence only on the gap time at time j – 1.

2.2 Nonparametric AR(p) models

The model in Section 2.1 can be extended to include higher order dependence, by modifying (2) – (3) as follows:

The distribution of α_ij for j ≤ p, depends only on the available past observations as in any AR(p) model. Prior specification for the remaining parameters is the same as described in Section 2.1.

3.1 Spike and slab base measure

Kim et al. [33] introduce the Spiked Dirichlet process prior in the context of regression. A key feature of their method is to employ the spike and slab distribution, i. e. a mixture of a point mass at 0 and a continuous distribution as centering distribution of the DP. This implies that, in a regression context, some coefficients have a positive probability of being equal to 0 and therefore they are not influential on the response of interest. Their technique is easily accommodated in our context by simply modifying G₀ in (7) as

where the introduction of hyperpriors on the weights of the mixture induces sparsity.

3.2 Prior on the Order of Dependence

Following Quintana and Müller [34], we specify a prior directly on the order p of the autoregressive process and then, conditioning on p, we set a Dirichlet Process prior of appropriate dimension for the parameters of the AR(p), i. e. the vector (mi0,mi1,…,mip). Let P be the maximum possible order. Then we can specify the following hierarchy:

When p = 0, the process simplifies as the autoregressive term in (9) disappears and the base measure of the DP reduces to the Normal distribution for the intercept term.

An MCMC algorithm to perform posterior inference on the order of dependence is described in [34], who sample from the posterior distribution of p given the data. For the JAGS implementation we have opted for a pseudo-prior approach, see [35].

4.1 Simulation scenario 1

We consider a simulated dataset of N = 300 subjects, with n_i randomly generated from a Poisson distribution. In particular we select n_i in such a way that all subjects have at least three gap times: ni∼P(2.5)+3, for all i. One third of the observations are generated from

while another third is generated from

and the last 100 observations are generated from

In simulating the data, we assume independence across subjects. In this example, for ease of explanation, we do not include covariates. We introduce censoring according to the following procedure. Firstly, we set the percentage of the sample for which we have censoring. In this simulation we have used three different censoring rates: 10 %, 25 % and 50 %. We consider at least two gap times for each individual and so for all patients i with censoring we generate a censoring time C_i from a Uniform distribution defined on the interval (Ti2,Tini). The results obtained without censored data and with the different censoring rates are compared.

We fit the model (1), (5)–(6) with p = 3:

where G₀ is given by the product of spike and slab distributions as defined in (8):

In fitting the model we set:

Hyperparameters are chosen in order to specify vague marginal prior distributions. Figure 1 shows the predictive distributions of m_i0, m_i1, m_i2 and m_i3 obtained with the spike and slab base measures. Red, blue, green and grey lines correspond to estimates obtained using the dataset with 0 %, 10 %, 25 % and 50 % censoring rates, respectively. By visual inspection, it is clear that the marginal posterior predictive distributions of the m_ij’s agree with the true values used to generate the dataset. In fact, the predictive distribution of m_i0 is concentrated around 0, while the predictive distributions of m_i1, m_i2 and of m_i3 are bimodal. Furthermore, the variance of the predictive distributions of each m_ij becomes larger with increasing censoring rate, because we have less information about the population parameters. The marginal posterior distributions of η₁, η₂, and η₃ in the prior (8), not reported here, concentrate most mass on 1.

Figure 1:

Simulation scenario 1: predictive marginal distributions of m_i0 (a), m_i1 (b), m_i2 (c) and m_i3 (d) under the spike and slab prior. Red, blue, green and grey lines correspond to estimates obtained using the dataset with 0 %, 10 %, 25 % and 50 % censoring rates, respectively.

Moreover, there is a clear difference between the posterior distribution of K, the number of distinct components in the mixture (5)–(6), using the data with and without censoring (results shown in Figure 2 of Supplementary Material). Using data without censoring, the configuration involving K = 3 clusters has the highest posterior probability: posterior inference on K is in agreement with the 3 components used to generate the data. Instead, using data with censoring, the number of clusters increases as the model needs to accommodate for the lack of information and the greater uncertainty in the distribution of the gap times. We also fit model (1), (9) to this dataset and similar results are obtained using the prior on the order of dependence. For further details, see Section 1 and Figure 1 of Supplementary Material.

4.2 Simulation Scenario 2

In this section we simulate a dataset of N = 200, with n_i randomly generated from a Poisson distribution. Similarly as in the previous example, we select ni∼P(2.5)+3 for all i, so that all subjects have at least three gap times. Half observations are generated independently from

while the other half is independently generated from

As in the previous example, covariates are not included and we introduce censoring using the same strategy. Also in this case we set the censoring rate equal to 10 %, 25 % and 50 %. We consider at least two gap times for each individual and for all the patients with censoring we generate a censoring time C_i from a Uniform distribution defined on the interval (Ti2,Ti10). The results obtained without censoring and with the different censoring rates are compared.

Also for this scenario, we fit both model (1), (5)–(6), (8) and model (1), (9). Here we report posterior inference of the latter model, with the prior on the order of dependence, where the maximum order of dependence is P = 3 and the prior hyperparameters (corresponding to a vague prior) are set as follows:

The mode of the marginal posterior distribution of p is 2 for all datasets without and with censoring, with corresponding posterior probability reported in the Table 1. In particular, the mass placed on p = 2 decreases with increasing censoring rate.

Conditional on p = 2, Figure 2 reports the predictive distributions of m_i0, m_i1 and m_i2. Red, blue, green and grey lines correspond to estimates obtained using the dataset with censoring rates equal to 0 %, 10 %, 25 % and 50 %, respectively.

More in detail, the predictive distribution of m_i0 is concentrated around 0, while the predictive distributions of m_i1 and of m_i2 are bimodal with mode around {–0.9, 0.9} and {–0.7, 0.7}, respectively. Also in this scenario, the variance of the predictive distributions of m_ij becomes larger with increasing censoring rate.

Figure 2:

Simulation scenario 2: predictive marginal distributions of m_i0 (a), m_i1 (b), m_i2 (c) and m_i3 (d) under the prior on the order of dependence, conditioning on p = 2. Red, blue, green and grey lines correspond to estimates obtained using the dataset with 0 %, 10 %, 25 % and 50 % censoring rates, respectively.

Finally, in the dataset without censoring time, conditioning on p = 2, the posterior mode for the number K of clusters is 2 (results not shown), with associated posterior probability equal to 0.5, while for the datasets with censoring, the number of clusters increases because of the greater uncertainty on the distribution of gap times. See Section 2 of Supplementary Material for the results obtained by fitting the model with the spike and slab base measure.

To evaluate the computational performance of the two approaches for model selection, we report in Table 2 summary statistics for the Effective Sample Size (ESS). ESS measures the amount by which autocorrelation within a MCMC chain increases uncertainty in estimates. The ESS of a parameter aims to estimate how many iid samples have the same estimation performance as the correlated samples found in the MCMC chain and it will be typically smaller than the final sample size. Hence large values of the ESS denote smaller autocorrelation in the chains; for further details, see [39]. Summary statistics are evaluated across the ESS of all estimated parameters. We report results for both simulations scenarios and four different censoring rates. The empirical distributions of the ESS in Table 2 seem comparable across different prior specifications and scenarios.

4.3 Comparison with the time-dependent frailty model

In this subsection, we compare our model with the time-dependent frailty model of [40] on the same simulated data described in the previous section. This strategy not only accounts for dependence of recurrent failure times within observations for the same subject by introducing a subject-specific random frailty, but also for the order of events, which is a significant feature of dependence. In this setup, time independence is introduced in the model for the intensity function of the recurrent event process at time t:

where the variable c_i(t) is the censoring indicator equal to 1 if patient i is observed at time t and 0 otherwise. The function λ₀(t) is an unspecified baseline intensity, while V_i(t) is a subject-specific random frailty capturing the risk not accounted for by potential risk variables included in the analysis. The model can be extended to include time-dependent covariates. Following [40] we assume λ0(t)=exp(β0) and a first-order autocorrelated process prior for the frailties:

where Vi(0)∼N(0,σW2) and e_i(t) are i.i.d. random variables having a N(0,σW2) distribution. The parameter ϕ is constrained to lie between – 1 and 1 (a priori ϕ∼ Translated Beta(3, 3), i. e. with density proportional to (ϕ+1)2(1−ϕ)2I(−1,1)(y)), and it measures the degree of serial correlation in the subject-specific frailty. The prior specification is completed as follows:

Also for this model posterior inference can be performed through a standard Gibbs sampler algorithm, which we implement in JAGS [36]. For each simulated scenarios, we run the algorithm for 251,000 iterations, discarding the first 1,000 iterations as burn-in and thinning every 50 iterations to reduce the autocorrelation of the Markov chain. The final sample size is 5,000 as before. Table 3 reports the posterior mean and standard deviation of the frailty correlation parameter ϕ, for the first simulation scenario, described in Section 4.1; Table 4 reports similar summaries of ϕ for the second simulation scenario (see Section 4.2).

From Table 3 it is clear that Manda and Meyer’s model is able to detect a significant positive correlation between recurrent events (as evident from the posterior means and standard deviations). In fact, in the first scenario, the data are generated so that observations at time j depend on observations at time t – j in two clusters, while in the third one the observations are i.i.d.

On the other hand, the posterior mean of ϕ in the second scenario is centred around 0. This shows the limitations of a parametric approach as observations are generated from two equal size groups, with the same order of dependence, but opposite effect. Therefore, it is difficult to capture the effect of the correlation between the recurrent event times. Moreover, from Figure 1 and Figure 2, it is clear that our model is more flexible. Both Figures show multimodal distributions for the autoregressive coefficients; this implies that our approach captures the clustering of observations and the cluster specific dependence structure. For example in Figure 1, the predictive distribution for m_i3 has a mode in 0, i. e. in one cluster there is no dependence of third order, while the other mode is centred over 0.4, i. e. in the other cluster there is a positive dependence of the third order.

5.1 Testing for the order of dependence

When testing the order of dependence, we first fit model (1), (5)–(6) and (8) with p = 3 (G₀ being the spike and slab distribution) and then model (1), but specifying a prior directly on model size as in (9) with P = 3.

Figure 3 reports the marginal posterior predictive marginal distributions of m_{i, l}, for l = 0, 1, 2, 3, obtained with the spike and slab base measure. Since the predictive distributions of mi0,mi1,mi2 are not concentrated around 0, unlike that of m_i3, we can conclude that the process best describing the readmission dataset has a dependency of the second order. This result is confirmed also using the approach described in Section 3.2. Indeed, the posterior distribution of p, displayed in Figure 9 of Supplementary Material, places most of its mass on 2.

Figure 3:

Readmission dataset: marginal posterior predictive distributions of m_i0 (a), m_i1 (b), m_i2 (c) and m_i3 (d) under the spike and slab prior.

5.2 Posterior analysis

We compare now the results of the nonparametric AR(2) model for the random effects α_ij’s as in (5)–(7), selected in the previous section, with models (2)–(3) built using different choices of g. In particular we consider two summary statistics: g(Yi1,…,Yij−1)=(Yi1+⋯+Yij−1)/(j−1) and g(Yi1,…,Yij−1)=Yi1×⋯×Yij−11/(j−1). The goal is to understand if higher order temporal dependency can be approximated by an AR(1)-type process built on some appropriate function of past observations as described in (2). We analyze the posterior distribution of K, the number of components in the mixture (5)–(6) under different alternatives. In particular, the results show that the posterior modes of K are 2 or 3 with a probability of around 30 % for the AR(1)-type models. On the other hand, the AR(2) model, suggests the existence of 3, 4 or 5 groups (See Figure 10 of Supplementary Material for results).

In Figure 4 we present the marginal posterior predictive distributions of Yjnew for a hypothetical new subject, for each time j,j=1,..., 6, setting the values of the covariates to the empirical mode. From the figure, it is evident that the two AR(1)-type models produce very similar results. Obviously, for j = 1 and j = 2 the three distributions are almost identical, as the models are closer. For j > 2, it is clear that the posterior predictive distributions of Yjnew have a larger variance and are more skewed under the AR(2) model. This experiment shows that it is not straightforward to approximate higher order dependency using summary statistics.

Figure 4:

Readmission dataset: posterior predictive distributions of Yjnew, j=1,…,6 for a hypothetical new subject, with covariates equal to the empirical mode of all covariates in the sample. The green and blue lines represent AR(1)-type models, with g(Yi1,…,Yij−1)=(Yi1+⋯+Yij−1)/(j−1) and g(Yi1,…,Yij−1)=Yi1×⋯×Yij−11/(j−1), respectively, and the red distribution indicates AR(2) model.

5.3 Posterior inference on the regression parameters

We now discuss the inference on the regression parameters in order to understand how covariates influence the recurrent event process. Although some covariates are fixed and do not vary over time, we still assume that their effect can be different in time and therefore we estimate a different vector of regression coefficients for all covariates in the model for each gap time j, j = 1, …,6. Covariates chemo and sex are binary variables, while dukes and charlson are 3 levels categorical variables and we need to introduce 2 dummy variables for each of them in the model, with baseline set to A–B for dukes and to 0 for charlson. Therefore, the final covariate vector for individual i is given by xi=(xi1,xi2,xi3,xi4,xi5,xi6)=(indicator for chemotherapy, indicator for female, indicator for dukes equal to C, indicator for dukes equal to D, indicator for charlson in 1–2, indicator for charlson equal to 3). The vector of regression parameters βj=(β1j,β2j,β3j,β4j,β5j,β6j) for each gap time j, j = 1, …, J (= 6), is therefore 6-dimensional.

Figure 5 displays 95 % credible intervals for the marginal posteriors of the regression parameter of each covariate for the first 5 gap times, i. e. of β˜r:=(βr1,βr2,βr3,βr4,βr5). In general there is no evident effect of chemotherapy on the outcome on any gap time. However, β˜2, which measures the effect of sex on the gap times, indicates that women have mainly larger gap times.

Figure 5:

Readmission dataset: 95 % posterior credible intervals for the regression parameters of each covariate across the first five gap times. (a) Regression coefficients of chemo. (b) Regression coefficents of sex. (c) Regression coefficients of dukes = C. (d) Regression coefficients of dukes = D. (e) Regression coefficents of charlson = 1–2. (f) Regression coefficents of charlson = 3.

The dummy variable relative to Dukes stage C of the tumour and the regression coefficients corresponding to Dukes stage D are never significant. Finally, posterior CI’s of the parameters corresponding to the Charlson index show that, in general, patients with index 0 will have larger gap times with respect to the one with index 3, in fact the medians of the posterior distribution are concentrated on negative values. Moreover, credible intervals are larger for the last gap times, as expected, since few individuals have a large number of events. The regression coefficients at time j = 6 are not shown as the credible intervals are not comparable with those of the previous times.

5.4 Comparison with existing methods

In this subsection, we compare our model with the frequentist shared frailty model of [41] as implemented in the R package frailtypack [42]. In this setup the dependence between recurrences is modelled through a frailty variable, such that events which occur several times within the same subject during the period of observation share an unobserved random effect, i. e. frailty. The hazard function for the j-th event (j=1,…,ni+1) and i-th individual (i=1,…,N), conditional on the frailty parameter ω_i, is

where h₀(t) is the baseline hazard function, β is the vector of the regression coefficients associated to the covariate vector xij. Table 6 reports the estimates of the regression coefficients and the associated p-value.

The estimate of the variance θ of the frailty term is 0.67 with standard error 0.14, so θ is significantly different from 0, meaning that there is heterogeneity between subjects. Differently from our approach, in this model the regression coefficients do no vary across time. The most significant variables are the indicator for female, indicator for dukes equal to D and indicator for charlson equal to 3. A negative coefficient implies a lower risk of rehospitalization, while a positive a higher risk. In our analysis we obtain that (i) the coefficient for indicator of female (sex=1) remains in general positive over time, implying longer times between rehospitalizations; (ii) the coefficient for duke equal to D is initially negative but becomes positive for the last gap times; (iii) the regression parameters corresponding to charlson index equal to 3 are mostly negative over time. These results show a good agreement between the two models. Finally, Section 6 of Supplementary Material reports a comparison of our model with the dynamic frailty models of [26], who also adopt a Bayesian nonparametric approach based on the DP.

6.1 Posterior inference

We perform inference on the order of dependence using both approaches introduced in Section 3. The marginal posterior predictive distributions of m_{i, l}, for l = 0, 1, 2, 3, obtained with the spike and slab base measure, are displayed in red in Figure 7. Panels (b), (c) and (d) show that the predictive distributions of m_i,1, m_i,2, m_i,3 are all concentrated around 0. In addition, specifying directly a prior on p, with P = 3, leads to a posterior distribution for the order of temporal dependence with mode in 0 (result shown in blue). These results indicate that there is no evidence of dependence between gap times, although the presence of a small cluster whose distribution is concentrated slight away from zero is detectable. This result is confirmed also by running the model for the three choices of function g described in Section 2.1, assuming an AR(1) structure. We obtain similar posterior predictive marginal distributions for m_i0 and m_i1, as well as the same posterior inference for K. In particular, a posteriori, the marginal distribution of m_i1 is concentrated around 0, supporting the hypothesis of independence between gap times.

Figure 7:

UTI dataset: marginal posterior predictive distributions of m_i0 (a), m_i1 (b), m_i2 (c) and m_i3 (d) obtained with the spike and slab base measure (red) and with the prior on the order of dependence (blue).