Sage Journals: Discover world-class research

Abstract

Illness-death models are a class of stochastic models inside the multi-state framework. In those models, individuals are allowed to move over time between different states related to illness and death. They are of special interest when working with non-terminal diseases, as they not only consider the competing risk of death but also allow us to study the progression from illness to death. The intensity of each transition can be modelled including both fixed and random effects of covariates. In particular, spatially structured random effects or their multivariate versions can be used to assess spatial differences between regions and among transitions. We propose a Bayesian methodological framework based on an illness-death model with a multivariate Leroux prior for the random effects. We apply this model to a cohort study regarding progression after an osteoporotic hip fracture in elderly patients. From this spatial illness-death model, we assess the geographical variation in risks, cumulative incidences and transition probabilities related to recurrent hip fracture and death. Bayesian inference is done via the integrated nested Laplace approximation.

Keywords

Bayesian inference integrated nested Laplace approximation multi-state models spatial correlation transition probabilities

1. Introduction

Multi-state models are stochastic models which generalise a wide class of survival scenarios, from unidimensional survival models to multi-event models such as competing risks models or repeated events.^1,2 In the multi-state framework, events are the states of the process, and their respective occurrences are transitions between the state of departure and the state of interest. The uncertainty associated with transitions is modelled via transition probabilities or, equivalently, transition intensities. The latter are analogue to hazard functions in the field of survival analysis. Multi-state models are especially useful in medical research because they provide a natural setting for dealing with the natural history of complex diseases.³

The so-called illness-death model⁴ is one of the simplest and most studied multi-state models. It has three states: an initial state, an illness-related state and a death state. The process starts in the initial state from which it can progress to the illness transient state or to death, which is an absorbent state. Death is also accessible from the illness state. This model is particularly useful to the study of chronic diseases,⁵ cancer progression⁶ or cardiovascular diseases⁷ in which there is a considerable risk of death over time.

The Cox proportional hazards model⁸ is the most popular regression tool in the survival framework to model hazard functions associated with survival times. It expresses hazard functions as the product of a time-dependent baseline hazard function and the exponential of a regression term including covariates and latent elements. The popularity of this model is primarily due to two reasons. First, due to the interpretability of hazard ratios to evaluate differences in the risk of the event of interest among the different covariate levels. Second, because under the frequentist paradigm, that the objective does not require to make any assumptions about the baseline hazard function. From Bayesian reasoning, a model for the baseline risk function needs to be specified,⁹ either parametrically or semi-parametrically.¹⁰ In particular, when we work only with covariates and use a Weibull baseline risk function, the overall risk function will also be Weibull. This property is the basis of the correspondence between the Weibull Cox proportional hazards model and the Weibull accelerated failure time (AFT) regression model.¹¹

The Bayesian paradigm provides a flexible framework for statistical inferences and the generation of knowledge. Under this framework, any measure of interest is subject to uncertainty: not only random variables but also parameters, hypotheses, models, etc. On the other hand, the Bayesian inferential process via the Bayes’ theorem allows for sequentially updated previous knowledge of those measures using new information. The procedure is conceptually simple: the elicitation of a prior distribution for all uncertainties in the model, the computation of the likelihood function for the data obtained, and the estimation of the posterior distribution which updates the relevant knowledge. This posterior is the starting point to approximate the posterior distribution of any measure of interest, such as sojourn times, transition and occupation probabilities, and cumulative incidence functions.

Regression survival models can include not only covariates but also latent effects that account for some non-explained heterogeneity between groups of the target population. In particular, the existence of differences among spatial regions is especially common in epidemiological studies. Uncontrolled risk factors may be relevant to explain high or low risks of disease in some regions, leading to this heterogeneity. We focus in this article on lattice data, that is, data for a finite number of sub-regions of a larger one. Often neighbouring regions can be expected to be similar and thus random effects with a spatial correlation structure can be assumed.

There are a plenty of models for assessing spatial correlation in the statistical literature. Conditional autoregressive (CAR) models¹² and their variants based on a neighbourhood definition of the correlation have been widely used in disease mapping. In particular, the model proposed by Besag, York and Molliè (BYM) in 1991¹³ has been postulated as the main choice over the past decades to deal with counts assuming a Poisson process. It is defined by means of two random effects, the first based on the neighbourhood structure and thus summarizing the spatial correlation between regions, and the second unstructured accounting for heterogeneity among regions. Leroux et al.¹⁴ proposed an alternative specification for the precision matrix of the spatially distributed random effects that better distinguishes between spatial dependence and dispersion effects. Under this model, random effects are defined as a mixture of independent and spatially correlated scenarios. Some authors assessed the behaviour of spatial models inside the survival framework such as Banerjee, Wall and Carlin,¹⁵ comparing different models without random effects (usually referred as frailties in the survival setting), with non-spatial frailties and with a CAR frailty.

Regarding illness-death models, not only spatial correlation can be modelled, but also a correlation between the three transitions, resulting in a multivariate model for random effects. In this regard, Carlin and Banerjee¹⁶ proposed a multivariate CAR model for spatially correlated survival times. However, despite its interest, there are few studies considering spatial components in the illness-death model framework. The most remarkable research work in this direction is Nathoo and Dean¹⁷ in which various structures for region-specific random effects are proposed, with special attention to the comparison of different baseline functions such as Weibull distributions, piecewise-exponential forms and cubic B-splines.

We propose a Bayesian methodological framework to deal with spatially correlated random effects within the illness-death scenario. In particular, a multivariate version of the Leroux model is used to jointly model that spatial correlation as well as the correlation between the transition survival times. The Bayesian procedure involving the approximation of the relevant posterior distribution is done via the integrated nested Laplace approximation (INLA), which in general provides accurate estimations and reduces the computational time compared to Markov chain Monte Carlo (MCMC) methods. In the context of the proposed methodological framework, the computation of posterior outcomes such as sojourn times, transition and occupation probabilities and cumulative incidence functions results in natural outcomes. Moreover, those quantities can be mapped providing rich information about the spatial distribution of illness and death in terms of probabilities that may have important clinical implications as interpreted by clinicians and epidemiologists. We apply this model to a real-world study involving recurrent hip fractures in old people. Data come from the PREV2FO cohort of patients from the Comunitat Valenciana (Spain) aged 65 and over who have been discharged from hospital after a hip fracture. They include individual baseline information of those patients and their progression over time. In addition to these individual characteristics, Health Areas where patients belong to are included to assess geographical differences in their corresponding risks and transition probabilities.

This article is structured as follows: illness-death model is presented in Section 2; the methodological framework regarding the proposed spatial illness-death model is described in Section 3, including the sampling model in Section 3.1, Bayesian inference and prior specification in Section 3.2, and some posterior outcomes in Section 3.3. Section 4 includes the analysis of the study of recurrent hip fracture, in particular, Section 4.1 presents the PREV2FO cohort, Section 4.2 includes some relevant results regarding posterior inference of the parameters of the model, and Section 4.3 includes the results regarding relevant measures of the process such as cumulative incidences and transition probabilities. Finally, we present a discussion in Section 5.

2. Illness-death models

Illness-death models are the most popular multi-state models.⁴ In their simplest version they comprise three states: an initial state (1), an illness-related state (2) and a death state (3). The death state is absorbent and accessible directly from the initial state or through the intermediate and transient state defined by the illness. Figure 1 depicts the model, including transitions between the states.

Figure 1.

Illness-death model with initial state (1), transient illness state (2) and death state (3). Arrows between the states represent the corresponding transitions.

From a probabilistic framework, an illness-death model is defined as a stochastic process ${Z (t), t > 0}$ in continuous time $t$ which takes as values the possible states where individuals can be. In particular, $S = {1, 2, 3}$ is the state space and $C = {1 \to 2, 1 \to 3, 2 \to 3}$ the set which includes all possible transitions of the process. We assume a semi-Markovian¹⁸ structure of the process whose evolution from the initial state to disease or to death only depends on the history of the process through the current state; but the transition from disease to death will depend not only on the present situation of the process but also on how long it has been in the initial state before jumping to the disease state.

The random behaviour of our illness-death model is determined by the initial distribution of the process, $P (Z (0) = i)$ $\forall i \in S$ (usually $P (Z (0) = 1) = 1$ because it is assumed that at $t = 0$ individuals are in state $1$ ), and the so-called transition probabilities between states defined as

\begin{aligned} p_{1 j} (s, t) = P (Z (t) = j ∣ Z (s) = 1), s \leq t, j = 2, 3 \\ p_{23} (s, t ∣ t_{12}) = P (Z (t) = j ∣ Z (s) = 1, T_{12} = t_{12}), s \leq t_{12} \leq t \end{aligned}

(1)

where

T_{12}

is the time the process spends in state

1

before entering into state

2

. Transition probabilities provide intuitive and easily interpretable information about the problem of interest but are difficult to model. For this reason, the statistical modelling usually relies on transition intensities, which are much less intuitive but easier to model. They account for the instantaneous hazard of progression to state

j

conditional on the current state

i

as follows

h_{i j} (t) = lim_{Δ t \to 0} \frac{p_{i j} (t, t + Δ t)}{Δ t}

(2)

Notice that transition

2 \to 3

needs to add the condition

T_{12} = t_{12}

to the subsequent transition intensity.

In the case of our illness-death model, transition probabilities are computed from transition intensities as indicated below.⁶ In particular, note that $p_{13}$ is the probability of total death, regardless of having passed through the state of illness (2) or not

\begin{aligned} p_{11} (s, t) = \exp {- \int_{s}^{t} (h_{12} (u) + h_{13} (u)) d u} \\ p_{22} (s, t ∣ t_{12}) = \exp {- \int_{s}^{t} h_{23} (u - t_{12} | t_{12}) d u} \\ p_{12} (s, t) = \int_{s}^{t} p_{11} (s, u) h_{12} (u) p_{22} (u, t | u) d u \\ p_{13} (s, t) = 1 - p_{11} (s, t) - p_{12} (s, t) \\ p_{23} (s, t ∣ t_{12}) = 1 - p_{22} (s, t) \\ p_{33} (s, t) = 1 \end{aligned}

(3)

The equivalence between some concepts from the world of stochastic processes and of survival analysis is quite natural for the given setting for the transition time from

i

j

can be seen as the survival time between the initiating event

i

and the entrance into the state of interest

j

T_{i j}

. In this regard, transition intensities in the stochastic framework (2) are equivalent to hazard functions in the survival setting as follows

h_{i j} (t) = lim_{Δ t \to 0} \frac{P (t \leq T_{i j} < t + Δ t ∣ T_{i j} \geq t)}{Δ t}

(4)

Transition intensities can be naturally modelled using Cox proportional hazard models,⁸ which expresses hazard functions by means of the product of a baseline hazard function and an exponential regression term with covariates, random effects or any other element that can provide information on the variable of interest.

3. Bayesian spatial illness-death modelling

We propose a Bayesian spatial illness-death model for a finite spatial lattice data with a set of local neighbourhoods defined by geographic vicinity between the sites of the target region. This model includes the joint modelling of the three relevant survival times of the illness-death model associated with each site as well as a spatial structure, based on the Leroux model,¹⁴ which connects the survival process of the different sites of the spatial domain. The notation of the Bayesian content that we will introduce from now on considers all probabilities and derived concepts to be conditional because all the parameters and hyperparameters on which they depend have probability distributions.

3.1. Sampling model

Let $h_{i j}^{(k)} (t ∣ θ)$ denote the conditional hazard function for survival time $T_{i j}$ associated with transition $i \to j \in C$ at time $t$ for an individual from region $k$ , $k = 1, \dots, K$ which we express through the Cox model⁸

\begin{aligned} h_{i j}^{(k)} (t ∣ θ, ψ) = h_{i j, 0} (t ∣ θ) \exp {η_{i j}^{(k)}} \\ η_{i j}^{(k)} = x^{'} β_{i j} + b_{i j}^{(k)} \end{aligned}

(5)

where

θ

is the vector of parameters and hyperparameters of the model,

ψ

the vector of all random effects,

h_{i j, 0} (t ∣ θ)

a baseline hazard function and

η_{i j}

a regression term defined in term of covariates

x = (x_{1}, \dots, x_{L})^{'}

, a vector of regression coefficients

β_{i j} = (β_{i j, 1}, \dots, β_{i j, L})^{'}

and a random effect

b_{i j}^{(k)}

associated with transition

i \to j

in region

k

. Baseline hazard functions can be approached in several ways. From parametric models such as the Weibull, undoubtedly the traditional and most widely used model in biometric applications, to more flexible modelling such as piecewise constant functions or B-splines.¹⁹ In our case, we propose a full parametric approach via Weibull baseline hazard functions defined as

h_{i j, 0} (t ∣ θ) = α_{i j} λ_{i j} t^{α_{i j} - 1}

Random effects $b_{i j}^{(k)}$ in (5) depend on the different model transitions as well as of the different sites in the target region. Let $B$ be a matrix which comprises all random effects

B = (\begin{array}{ccc} b_{12}^{(1)} & b_{13}^{(1)} & b_{23}^{(1)} \\ ⋮ & ⋮ & ⋮ \\ b_{12}^{(k)} & b_{13}^{(k)} & b_{23}^{(k)} \\ ⋮ & ⋮ & ⋮ \\ b_{12}^{(K)} & b_{13}^{(K)} & b_{23}^{(K)} \end{array})

(6)

B (, i)

the

i

-th column of the matrix

B

, and

v e c (B) = [B (, 1)^{'}, B (, 2)^{'}, B (, 3)^{'}]^{'}

a column vector

(3 K \times 1)

including each of the columns of matrix

B

. We assume a conditional multivariate Gaussian Markov random field for the random effects

v e c (B)

with a mean vector whose elements are all zero and a matrix of variances-covariances

Σ

(v e c (B) ∣ Σ) \sim N_{3 K} (0, Σ)

(7)

It is worth mentioning that

v e c (B)

is precisely the set of random effects that we have generically represented as

ψ

above.

The structure of $Σ$ includes both multivariate dependence between the illness-death transitions of a given site (between columns of $B$ ) and spatial dependence for each transition of the model (within columns of $B$ ) in the form

Σ = Σ_{b e t w e e n} \otimes Σ_{w i t h i n}

(8)

where

\otimes

represents the Kronecker product.

Models for the spatial variability $Σ_{w i t h i n}$ have a wider range of options: from the simplest independent scenario to the conditional autoregressive (CAR) models¹² and their variants (intrinsic CAR, proper CAR, Besag York & Molliè, Leroux model). The proposal by BYM¹³ has been widely used for disease mapping in the epidemiological literature over the past decades. Part of its popularity remains in its interpretability, as it consists of a random effect, which considers the spatial correlation between regions according to a neighbourhood structure, and an unstructured random effect accounting for heterogeneity among regions. With the BYM model, however, only the sum of both sets of random effects is identifiable, failing thus to identify both random effects separately.²⁰ The Leroux model¹⁴ circumvents this problem since there is only one set of random effects defined in terms of a mixture of independent and spatially-dependent elements that allow to assessing the intensity of them as follows

Σ_{w i t h i n} = (τ [(1 - γ) I + γ (D - W)])^{- 1}

(9)

where

τ

is a dispersion hyperparameter,

I

the identity matrix,

D

a diagonal matrix whose non-zero elements on the diagonal are the number of neighbours in the corresponding site,

W

an adjacency matrix, that is,

W_{k l} = 1

if sites

k

and

l

are neighbours,

k \neq l

, and 0 otherwise, and hyperparameter

γ \in [0, 1]

determines how matrices

I

and

D - W

are combined. A value of

γ = 0

simplifies to an independent random effects model without spatial patterns, whilst

γ = 1

corresponds to an intrinsic CAR model.

We model the variance-covariance matrix $Σ_{b e t w e e n}$ of the times between the three transitions as

Σ_{b e t w e e n} = (\begin{array}{ccc} \frac{1}{τ_{12}} & \frac{ρ_{(12) (13)}}{\sqrt{τ_{12} τ_{13}}} & \frac{ρ_{(12) (23)}}{\sqrt{τ_{12} τ_{23}}} \\ \frac{ρ_{(12) (13)}}{\sqrt{τ_{12} τ_{13}}} & \frac{1}{τ_{13}} & \frac{ρ_{(13) (23)}}{\sqrt{τ_{13} τ_{23}}} \\ \frac{ρ_{(12) (23)}}{\sqrt{τ_{12} τ_{23}}} & \frac{ρ_{(13) (23)}}{\sqrt{τ_{13} τ_{23}}} & \frac{1}{τ_{23}} \end{array})

(10)

This matrix includes two types of hyperparameters:

τ_{i j}

, the marginal precision of the random effects associated to transition

i \to j

, and

ρ_{(i j) (i^{'} j^{'})}

, the correlation between the random effects on transitions

i \to j

and

i^{'} \to j^{'}

. Note that for identifiability reasons, we fix

τ = 1

for the multivariate version of the Leroux model, so that precision parameters

τ_{i j}

from the covariance matrix

Σ_{b e t w e e n}

become the ones capturing dispersion.

3.2. Bayesian inference and prior specification

A Bayesian approach based on the integrated nested Laplace approximation (INLA)²¹ has been considered to estimate the posterior distribution of all the quantities of interest of the model. Bayesian inference combines prior knowledge of all unknown parameters and hyperparameters $θ$ of the model in probabilistic terms throughout the prior distribution with the likelihood function obtained from data $D$ by means of the Bayes’ theorem to derive the joint posterior distribution $π (θ, ψ ∣ D)$ of the parameters and hyperparameters $θ$ and random effects $ψ$ . As models get more complex, it is harder to find an analytic expression for those posterior distributions and computational methods are required to approach them. The most popular procedures are MCMC methods²² which, in most cases, imply large computational times to ensure convergence of the estimations. Alternatively, the INLA is a fast and accurate option. It uses Laplace approximations to obtain the approximated marginal posterior distribution of the parameters, hyperparameters and latent terms of the sampling model. Survival models, including Cox proportional hazards models, can be adapted and implemented in INLA because they can be expressed in terms of Gaussian Markov random field models.²³ In particular, competing risks models,²⁴ and illness-death models as an extension of them, can be approached using INLA. It also allows the inclusion of Gaussian random effects in the regression term of the Cox proportional hazards model, and thus the proposed spatial illness-death model is naturally approachable with INLA.

As just discussed, we need to complete the Bayesian model with a prior distribution for all parameters and hyperparameters $θ$ of the sampling model. We have considered a framework of prior independence between the different elements in $θ$ . The shape parameters $α_{i j}$ of the baseline hazard functions in (5) were assumed to follow a penalised complexity prior (PC prior) as described in the INLA documentation (see inla.doc(‘pc.alphaw’) for a detailed definition). Those PC priors consider an exponential as the base model, that is, a Weibull model with $α_{i j} = 1$ , and penalise the departure from this exponential model. The more general Weibull model would be preferred only if there exists enough evidence supporting it Van Niekerk et al.²⁵ Meanwhile, the scale parameter for the baseline transition intensities, $λ_{i j}$ , has not a prior by itself, but through considering an intercept $β_{i j, 0} = l o g (λ_{i j})$ which follows, as well as the regression coefficients $β_{i j, l}$ , a Gaussian distribution with mean 0 and precision $0.001$ .

The covariance matrix including the correlation between transitions, $Σ_{b e t w e e n}$ , was assumed to follow an inverse Wishart distribution, or equivalently, a Wishart distribution for the precision matrix $Σ_{b e t w e e n}^{- 1}$ . The Wishart distribution is a multivariate generalization of the gamma distribution.²⁶ It is specially relevant when modelling correlated normal random effects as it is a conjugated prior distribution for the precision matrix from multivariate normal distributions, being the most common choice when inferring covariance matrixes. In our case, the prior values for the Wishart parameters were those provided by default in the INLA specification of the model for correlated random effects, that is, $Σ_{b e t w e n}^{- 1} \sim {Wishart}_{3} (ν = 7, R = I)$ , where $ν = 7$ are the degrees of freedom. Wishart distribution with the identity as the scale matrix is typically set as a relatively uninformative prior. Note however that several authors have discussed its appropriateness and some alternatives have been proposed. For instance, this prior specification might not be appropriate in the presence of parameters with small variances, resulting in a strongly informative prior distribution.²⁷ On the other hand, separation-strategies decompose the covariance matrix into variance and correlation components, being it possible to specify separate priors for each component. The correlation matrix obtained after this separation may be assumed to follow an inverse Wishart distribution,²⁸ or a Lewandowski-Kurowicka-Joe,²⁹ for instance. Meanwhile, many positive priors can be set for variances such as truncated-normal, half-normal, half-Cauchy or uniform distributions.²⁶

A non-informative uniform prior, $U (0, 1)$ was assumed for the mixture parameter $γ$ of the Leroux modelling. To define the multivariate Leroux model for random effects we used the rgeneric latent effect. Using this mechanism, latent effects can be implemented in INLA via R.³⁰ Despite they are not specifically applied to survival models, some multivariate versions of random effects have already been defined using this method, such as intrinsic multivariate CAR latent effects, and collected in the INLAMSM package for R.³¹

Note that when working with INLA, the number of hyperparameters is an important element to take into consideration. INLA can deal with models with a large number of components in the latent field, including fixed effects for covariates and random effects. However, a small number of hyperparameters is a critical assumption required to ensure accuracy in the approximations and for computational reasons. Rue et al.²¹ pointed out that a moderate amount of hyperparameters would be in the range of 6–12, while Rue et al.³² set a maximum of 20 hyperparameters. In our analysis, we have hyperparameters from two different sources: the illness-death model and the spatial structure. As we assumed baseline risk Weibull distributions in the multistate model, we would have three $α$ -shape hyperparameters. In the spatial model, the correlation matrix of the random effects has a $γ$ hyperparameter from the Leroux prior, three between-transition correlations $ρ$ , and three marginal precisions $τ$ . We have thus a total number of 10 hyperparameters, which is in the desired range.

3.3. Posterior outcomes

The posterior distribution $π (θ, ψ ∣ D)$ contains all updated information on the random behaviour of the illness-death population. Nevertheless, it provides unclear practical evidence about the prognostic clinical status of a patient over time. Compound measures such as sojourn time distributions, transition and occupation probabilities and cumulative incidence functions^33,34 mix the time-evolution information from the illness-death setting as well as the risk-variation among regions. They are specially relevant in order to gain insight into the clinical setting. From a statistical point of view, posterior inferences for those compound quantities is straightforward as they are indeed defined as functions of $(θ, ψ)$ . The procedure for doing so is as follows: first, we sample from the joint posterior distribution, obtaining samples of parameters, hyperparameters and random effects; second, we apply the target function to all elements of the sample and obtain a sample from the output of interest, which indeed approximates its posterior distribution; finally, we can summarise this distribution by means of the sample mean, median or credible intervals, etc, which are approximations of the posterior mean, posterior median and posterior credible intervals of the corresponding outcome.

It is worth mentioning that, by default, INLA provides approximate samples from the marginal posterior distributions and not from the joint posterior. In that sense, sampling from marginal information to infer about the outcomes of interest is not appropriate because it can produce misleading results that do not take into account the possible relationships between the components of the latent field. Fortunately, the INLA is also able to sample from the joint posterior distribution.^30,35 That is what we have done in this article. In the following, we will introduce some of those outcomes discussed above and discuss their posterior estimation.

3.3.1. Sojourn times

Sojourn time in state $i$ for an individual living in the site $k$ refers to the time an individual remains in that state without leaving. Possibly, the most interesting sojourn time in illness-death models corresponds to the initial state. It is defined in terms of the conditional survival function as follows

\begin{aligned} S_{1}^{(k)} (t ∣ θ, ψ) & = P (T^{(k)} > t ∣ θ, ψ) \\ = exp (- \int_{0}^{t} (h_{12}^{(k)} (u ∣ θ, ψ) + h_{13}^{(k)} (u ∣ θ, ψ)) d u) \end{aligned}

(11)

where

T^{(k)} = min {T_{12}^{(k)}, T_{13}^{(k)}}

. Because sojourn time in state

1

depends on

(θ, ψ)

, the subsequent posterior distribution

π (S_{1}^{(k)} (t ∣ θ, ψ) ∣ D)

\forall t

, can be easily approximate from a simulated sample of the posterior

π (θ, ψ ∣ D)

3.3.2. Transition and occupation probabilities

Transition probabilities depend on the parameters through the subsequent hazard functions according to (3). Therefore, their posterior distribution $π (p_{i j}^{(k)} (s, t ∣ θ, ψ ∣ D)$ associated with an individual in region $k$ will be also computed from a simulated sample of $π (θ, ψ ∣ D)$ . Occupation probabilities refer to the probabilities associated with the presence of the process in each of the different states at a given time $t$ . They can be expressed as transition probabilities $π (p_{i j}^{(k)} (0, t ∣ θ, ψ ∣ D)$ , and consequently its posterior distribution could be also approximated from an approximate sample from $π (θ, ψ ∣ D)$ .

3.3.3. Cumulative incidence functions

Cumulative incidence functions are more frequently used in competing risks environments, but they are also useful for illness-death models, specially when the illness is relevant by itself and not only as an intermediate state between the initial state and death. They can be defined equivalently to the competing scenario for survival times $T_{12}^{(k)}$ and $T_{13}^{(k)}$ as follows

\begin{aligned} F_{12}^{(k)} (t ∣ θ, ψ) = & P (T^{(k)} \leq t, η^{(k)} = 1 ∣ θ, ψ) \end{aligned}

(12)

\begin{aligned} F_{13}^{(k)} (t ∣ θ, ψ) = & P (T^{(k)} \leq t, η^{(k)} = 0 ∣ θ, ψ) \end{aligned}

(13)

where

η^{(k)}

is the indicator function with value 1 if

T_{12}^{(k)} < T_{13}^{(k)}

and 0 otherwise. They can be interpreted as the probability at time

t

of having moved directly from the initial state 1 to state

j

j = 2, 3

, keeping this sense of accumulation as its name suggests. Cumulative incidence regarding the illness state 2 is highly informative because it indicates how many individuals are expected to suffer the illness. It can also be directly compared with the transition probability from state 1 to 2, which indicates the expected rate of patients who experienced illness and are still alive. Cumulative incidence functions are also expressed in terms of

(θ, ψ)

F_{1 j}^{(k)} (t ∣ θ, ψ) = \int_{0}^{t} h_{1 j}^{(k)} (s) \exp {- \int_{0}^{s} (h_{12}^{(k)} (u) + h_{13}^{(k)} (u) d u} d s, j = 1, 2

(14)

Consequently, the posterior distribution of each of these cumulative incidences,

π (F_{1 j}^{(k)} (t ∣ θ, ψ) ∣ D)

, can also be approached by simulated samples of the posterior distribution

π (θ, ψ ∣ D)

4. A study of recurrent hip fractures in elderly patients

Clinical settings involving the progression of non-terminal diseases, repeated events and populations with a considerable competing risk of death are the main scenarios where multi-state models can be applied. We illustrate here the application of the previous model on a study of recurrent hip fracture.

4.1. The PREV2FO cohort

We analyse the PREV2FO cohort, a population-based cohort comprising patients aged 65 years and older discharged after hospitalization for an osteoporotic hip fracture in the Valencia Region (Spain) from 1 January 2008 to 31 December 2015.³⁶ The Valencia Region is an autonomous community of Spain, with a population of roughly 5 million people (10% of the Spanish population). The region provides universal healthcare services through the Valencia Health System (VHS) which is an extensive network of public hospitals, primary care centres and other public resources managed autonomously by the regional government. It is divided into 24 Health Areas, each one corresponding to the administrative area of influence of a public hospital from the VHS.

Patients were followed after the index fracture until death or end of study (31 December 2016), accounting for recurrent hip fractures during the follow-up period. Figure 2 shows a diagram of this process as and illness-death model with an initial state of discharge after a first hip fracture (F), an intermediate state that accounts for discharge after a refracture (R) and the state of death (D). Transition times between states were right-censored only due to end of study or death (see Supplemental Table 1 for specific information about the number of patients that progresses between states by sex and Health Area).

Figure 2.

Illness-death model with an initial state of hip fracture, a recurrent hip fracture state and a death state.

From a clinical point of view, there is a possibility of more than one refracture. We have dismissed this possibility because in our study only a reduced number of the patients suffered from them. In particular, we have 2532 patients with a refracture, only 26 of them with a second refracture and only one with a third. Due to the complexity of the model and the limited information regarding more than one refracture, we decided not to include this in the model and therefore, the exact definition of our refracture state would be thus ‘having at least one refracture’. Moreover, note that our model is easily generalizable and in case we had a higher number of patients with second or third refractures, we could add the corresponding states to represent those additional transitions.

In order to define a basic patient profile we have considered sex, age at the discharge and the Health Area in which patients were hospitalised as covariates. The study involved 34491 patients discharged alive after hip fracture, 25807 (74.8 $%$ ) were women and 8684 (25.2 $%$ ) men. Age was included as a continuous and mean-centered predictor in the model. The mean age at the first fracture was 83.4 years (IQR: 79.0–88.3). Patients were followed a median time of 5.0 years (IQR: 3.0-7.0 years). By age group, 12.4 $%$ of patients were under 75 years old, 43.6 $%$ between 75 and 85 years old, 40.6 $%$ between 85 and 94 years old and 3.4 $%$ were over 95 years old.

Survival times from state $F$ to $D$ , from $F$ to $R$ and from $R$ to $D$ in the Health Area $k$ , $k = 1, \dots, 24$ , $T_{F D}^{(k)}$ , $T_{F R}^{(k)}$ and $T_{R D}^{(k)}$ , respectively, are modelled by means of a Bayesian spatial Cox proportional hazards model as proposed in the previous section.

4.2. Posterior distribution

We present the approximated posterior distribution sequentially, first the parameters, then the hyperparameters and finally, the random effects. Table 1 summarises the approximate posterior marginal distribution of all parameters of the model. Estimations of the shape parameters of the baseline risk functions, $α_{F D}$ , $α_{F R}$ and $α_{R D}$ , indicate decreasing hazards over time, specially for the risks of death without and after refracture. $α_{F R}$ is closer to 1 which is the threshold which changes the behaviour of the Weibull hazard functions, from increasing to decreasing. Women and men showed no relevant differences in the risk of recurrent hip fracture (E $(β_{F R, W o m a n} ∣ D) = 0.021$ ), whereas women showed lower mortality risks as compared to men (E $(β_{F D, W o m a n} ∣ D) < 0$ , E $(β_{R D, W o m a n} ∣ D) < 0$ ). Age was found as a risk factor for refracture and for death without and after refracture.

Table 1.
Summary of the approximate posterior distribution of the parameters from an illness-death model with a multivariate-Leroux model for random effects. Transition-related parameters: shape of Weibull distribution and regression coefficients.

Time Parameter Mean Median SD 2.5% 97.5%

From $F$ to $R$ $α_{F R}$ 0.921 0.921 0.016 0.891 0.953

$λ_{F R}$ 0.028 0.027 0.005 0.018 0.040

$β_{F R, W o m a n}$ 0.021 0.021 0.050 −0.076 0.119

$β_{F R, A g e}$ 0.024 0.024 0.003 0.018 0.030

From $F$ to $D$ $α_{F D}$ 0.776 0.776 0.005 0.766 0.786

$λ_{F D}$ 0.335 0.331 0.054 0.238 0.460

$β_{F D, W o m a n}$ −0.510 −0.510 0.017 −0.543 −0.477

$β_{F D, A g e}$ 0.070 0.070 0.001 0.068 0.073

From $R$ to $D$ $α_{R D}$ 0.628 0.628 0.016 0.597 0.659

$λ_{R D}$ 0.593 0.579 0.131 0.374 0.897

$β_{R D, W o m a n}$ −0.634 −0.634 0.065 −0.761 −0.505

$β_{R D, A g e}$ 0.049 0.049 0.005 0.040 0.059

Time	Parameter	Mean	Median	SD	2.5%	97.5%
From $F$ to $R$	$α_{F R}$	0.921	0.921	0.016	0.891	0.953
	$λ_{F R}$	0.028	0.027	0.005	0.018	0.040
	$β_{F R, W o m a n}$	0.021	0.021	0.050	−0.076	0.119
	$β_{F R, A g e}$	0.024	0.024	0.003	0.018	0.030
From $F$ to $D$	$α_{F D}$	0.776	0.776	0.005	0.766	0.786
	$λ_{F D}$	0.335	0.331	0.054	0.238	0.460
	$β_{F D, W o m a n}$	−0.510	−0.510	0.017	−0.543	−0.477
	$β_{F D, A g e}$	0.070	0.070	0.001	0.068	0.073
From $R$ to $D$	$α_{R D}$	0.628	0.628	0.016	0.597	0.659
	$λ_{R D}$	0.593	0.579	0.131	0.374	0.897
	$β_{R D, W o m a n}$	−0.634	−0.634	0.065	−0.761	−0.505
	$β_{R D, A g e}$	0.049	0.049	0.005	0.040	0.059

Figure 3 shows the posterior expectation of the baseline hazard function associated to each of the three survival times. Note that the risks of death after recurrent hip fracture are higher than those of death without refracture. Transition intensity from fracture to refracture is notably lower than transitions to death. Note that baseline hazard functions are indeed the hazard functions for the reference values of predictors: average-aged men from a Health Area with a random effect equal to 0. Baseline functions suggest higher hazards during the first year, including the hazard of refracture, despite it cannot be appreciated graphically. It results in a sharper increase in the cumulative incidence of those events, as well as greater increases or decreases in the transition probabilities during the initial follow-up.

Figure 3.

Posterior mean of the baseline hazard function for each survival transition: from $F$ to $R$ , from $F$ to $D$ and from $R$ to $D$ . Horizontal axis indicates years from initial fracture, for transitions $F \to R$ and $F \to D$ , whilst years from refracture for the transition $R \to D$ .

Table 2 presents a summary of the approximate posterior marginal distribution of the hyperparameters of the spatial illness-death model, all them associated with the variability between the transition survival times and within the different Health Areas in the Valencia Region. The estimation of the parameter $γ$ from the Leroux model is $0.841$ thus indicating that the mixture of an independent scenario and an intrinsic CAR model lends toward the second (Figure 4). A 95% credible interval excludes lower values suggesting a relevant spatial correlation between areas. Correlation parameters between transitions showed posterior distributions not only including 0 but also zero-centred, which indicates irrelevant correlation parameters thus indicating an uncorrelated scenario. The highest value however was estimated for the correlation between death without refracture and death after refracture, $ρ_{(F D) (R D)}$ , showing a slight correlation between both types of mortality. Uncertainty about random effects is given by precision parameters $τ$ . Higher precision estimations indicate lower variability among random effects. Although the magnitude of the three is very similar, ordered from least to most uncertainty, we have random effects with transitions of death without refracture, refracture and death after refracture.

Figure 4.

Approximated posterior distribution of the $γ$ parameter from an illness-death model with multivariate-Leroux random effects.

Table 2.

Summary of the approximate posterior distribution of the hyperparameters from the illness-death model estimated with a multivariate-Leroux model for the spatial random effects. $γ$ parameter from the Leroux model, precision of the random effects and correlation between transitions times.

Parameter	Mean	Median	SD	2.5%	97.5%
$γ$	0.841	0.862	0.101	0.591	0.973
$τ_{F R}$	14.257	13.581	4.620	7.197	25.185
$τ_{F D}$	19.896	19.228	5.595	10.915	32.737
$τ_{R D}$	11.743	11.137	4.181	5.386	21.625
$ρ_{(F R) (F D)}$	−0.044	−0.047	0.181	−0.388	0.315
$ρ_{(F R) (R D)}$	−0.076	−0.078	0.178	−0.415	0.275
$ρ_{(F D) (R D)}$	0.109	0.111	0.164	−0.217	0.423

Figure 5 displays the posterior mean of the random effects associated with each transition time and Health Area of the Valencia Region. Health Areas coloured red indicate a higher risk of experiencing the event of interest compared to the overall average for all areas. Areas shaded in yellow indicate the opposite. The random effects associated with the three survival times of the illness model from the same Health Area do not always behave the same. We can observe some areas with positive random effects in the three survival times considered, but also some cases where the effects show negative relationships. There are some particular areas with a particular spatial pattern. This is the case of Requena-Utiel (the most western Health Area) and Denia (located at the cape in the east of the Valencia Region). The first shows a lower risk of recurrent hip fracture and a higher risks of death without and after refracture. The latter shows the opposite scenario, higher risk of refracture and lower mortality. Both cases illustrate a negative association between the risk of refracture and mortality, whilst a positive association between both risks of death.

Figure 5.

Posterior mean of the region-specific random effects by Health Area of the Valencia Region, from an illness-death model with multivariate-Leroux random effects.

4.3. Outcome measures of the hip fracture process

The examination of the raw estimations provided by the posterior distribution $π (θ, ψ ∣ D)$ contains full information about the differences in the risk of each outcome. Nevertheless, it provides unclear evidence about which will be the prognostic for a patient with a hip fracture in each particular Health Area or what would be the general evolution of the survival time transitions in the target population. Information regarding time-evolution from the illness-death setting and the variation in the risk among Areas are indeed combined in posterior distributions for cumulative incidences and for transition probabilities.

The cumulative incidence of a hip refracture at time $t$ can be interpreted as the probability of having a hip fracture at a time before $t$ without having died before that time, as death plays a major role in censoring refracture. Note that a higher risk of death leads to the observation of fewer refractures. Therefore, two Areas with the same risk of refracture could show different incidences of refracture depending on the risk of death. Figure 6 shows the posterior mean of the cumulative incidence of refracture for 80-year-old women and men in the different Health Areas of the Valencia Region at $t = 1, 2, \dots, 5$ years after a first hip fracture. In broad terms, higher incidences of recurrent hip fracture are estimated for those regions with a higher risk of refracture as it is expected. Those differences become more visible after some years from the initial fracture. The Health Area of Requena-Utiel (the most western region) shows a particular low incidence despite its not so low risk, which can be related to being the region with the highest risk of death without refracture. Men show lower incidences of refracture due to their increased risk of death, as we found no differences in the risk of refracture compared with women. Men reach the same incidence values than women with a delay of 1–2 years approximately.

Figure 6.

Posterior mean of the cumulative incidence of refracture in 80-year-old women and men at $t = 1, 2, \dots, 5$ years after a first fracture by Health Area.

Regarding transition probabilities from fracture to refracture (Figure 7), they are also higher for women, as they are more likely to experience a refracture than men. It shows an increasing trend during 2 and 4 years from the initial fracture for men and women, respectively. After this time, the probability of being refractured and alive remains stable, as the number of patients at risk of refracture decreases and the mortality after refracture offsets the number of new refractures.

Figure 7.

Posterior mean of the transition probability from fracture to refracture ( $p_{F R}$ ) in 80-year-old women and men at $t = 1, 2, \dots, 5$ years after a first fracture by Health Area.

Mortality is higher in men for both, total mortality and after recurrent hip fracture, only. Women approximately reach at 4 years the same mortality than men at 2 years (Figure 8). This difference is even higher for death after a refracture. Women reach at 5 years after refracture the same mortality rates than men only 2 years after a refracture (Figure 9).

Figure 8.

Posterior mean of the total-death probability ( $p_{F D}$ ) in 80-year-old women and men at $t = 1, 2, \dots, 5$ years after a first fracture by Health Area.

Figure 9.

Posterior mean of the probability of death after refracture ( $p_{R D}$ ) in 80-year-old women and men at $t = 1, 2, \dots, 5$ years after a refracture by Health Area.

The number of patients who die after a refracture represents a low fraction of the total mortality. Cumulative incidence of refracture indicates that less than 10% of women experience a refracture 5 years after the initial fracture (even lower in men). Thus, the spatial pattern of the total-death probability (Figure 8) is similar to that showed by the random effects on the risk of death without refracture (Figure 5).

The probability of death is higher for those patients with a recurrent hip fracture. Mortality one year after refracture is similar to that expected two years after the initial fracture. Its spatial pattern is also different with respect to that shown by total-mortality, and is identical to that shown by the respective random effects on the transition from refracture to death in Figure 5. This is due to the fact that the probability of death after refracture is the only one which depends exclusively on one transition intensity, in particular, the transition intensity from refracture to death.

4.4. Stability of the inference

The stability of posterior inferences to the specification of the different elements of the Bayesian model is an important issue in Bayesian analysis.³⁹ In this section, we only intend to take a brief look at the subject by analysing, in a non-exhaustive way, the sensitivity of the posterior results to the prior specification of some of the elements of the model. Obviously, these analyses do not demonstrate the robustness of the model, but they do show evidence in its favour. Likewise, and as the comparison between MCMC and INLA methods is always interesting, we have worked with a simpler model and compared the results provided by the INLA and those obtained through MCMC using JAGS.

We start focusing on the posterior inferences carried out by modifying the prior information on some hyperparameters of the spatial precision matrix. We assessed both differences in the estimations of the random effects as well as parameters and hyperparameters. Firstly, we considered the $γ$ parameter, which is a relative measure of spatial dependence. Note that in order to work with a proper prior distribution for $γ$ , it takes values in the range $[0, 1)$ , and thus we cannot reach the upper bound. We ran our model with $γ$ values of 0, 0.5 and 0.99. Random effects slightly differ despite constraining our model to a specific amount of spatial correlation (Supplemental Figure 3), and estimations of the parameters remained stable (Supplemental Table 4). The scenario with a value of $γ = 0$ , implying no spatial correlation at all, differed the most. We also set the degrees of freedom for the between-transition precision matrix to 6, 8, 9 and 10 (7 is the reference value that we used in our work). Again, posterior estimates showed few differences (Supplemental Figure 4; supplemental Table 5). Of course, the transition from refracture to death showed more variability as the sample size is quite low (Supplemental Table 1).

For the comparison between MCMC methods, via JAGS, and INLA we worked with a less complex Bayesian illness-death model and reduced the space to a group of five neighbouring Health Areas. The model (see Supplemental File 1) includes Gaussian random effects with between-transition correlation, but without spatial correlation. A post-sweeping has been applied to random effects obtained from JAGS, leading to identifiable intercepts and random effects.^37,38 Convergence was assessed graphically and via Gelman and Rubin’s $\hat{R}$ (Supplemental Table 2). Posterior estimations of the parameters and random effects from INLA and MCMC via JAGS were similar (Supplemental Figures 1 and 2) and since both procedures sample from the joint posterior distribution, we expect similar posterior inferences about the posterior outcomes. In summary, we found that INLA inferences were reasonable for those types of models, allowing us to benefit from INLA faster inference.

On the other hand, we estimated Monte Carlo Standard Errors (MCSEs)⁴⁰ for the samples of the posterior outcomes (Figures 6 to 9). Note that we first simulated values of the parameters from the joint posterior distribution with INLA, and then we obtained samples of cumulative incidences and transition probabilities. They are thus susceptible to Monte Carlo uncertainty, and we should provide the MCSE of those samples, analogously to what is done when simulating from MCMC or other methods which work with random draws. We calculated those errors for each Health Area, $k = 1, \dots, 24$ , for each $t = 1, \dots, 5$ , for each quantity of interest, and for both women and men (Supplemental Table 3). In broad terms, MCSEs were low compared to posterior standard deviations.

5. Discussion

The potential and usefulness of illness-death models are highly increased after combining them with spatial information modelled by multivariate random effects associated with a set of spatial units. Differences among regions can be studied under this joint framework, in addition to the progression of individuals over time. Despite the correlation between transitions was not observed to be relevant in our real-world study, the possibility of modelling it jointly with the spatial correlation from a model based on a neighbourhood structure, such as the Leroux model, defines an extensive number of options depending on the needs of the clinical framework. Moreover, one might be interested in using more complex multi-state models with additional states and transitions. It would be analogous as the concept throughout our work regarding how random effects are included in transition intensities is quite general and, thus, it is not restricted to illness-death models.

Some extensions of our work regarding censoring and truncation⁴¹ could be implemented without much difficulty. Patients with a hip fracture require hospitalization and, therefore, we exactly know when events occur. As a result, we had no other censoring than administrative censoring (i.e. due to end of study) or death. Note, however, that interval-censoring could arise in many survival health applications in which visits are the only way to know patients’ status.⁴² A patient could reach an illness state (state 2) between two consecutive visits and thus we would only know the time interval in which the event of interest occurred. On the other hand, we have considered time since fracture as the timescale of our analysis. Patients enter the process at that time and are followed up thereafter. An alternative study of the process according to a time scale marked by the age of the patient,⁴³ with time zero determined by the age of 65 years, would generate left-truncated data. They would correspond to those patients who made their debut in the study, leaving the hospital already healed from the first fracture, at ages over 65 years. INLA functions for Weibull survival models can deal with right censored data, as in our case, but also interval-censored and left-truncated.

The assumption of Weibull hazard functions in the study appears to be quite consistent. Weibull hazards can only be either increasing or decreasing (as in our study). For instance, regarding recurrent hip fracture and death, decreasing hazards seem reasonable as it is known that the occurrence of those events is particularly high during the first year after a fracture. However, after some time, the risk of death might be expected to increase, showing a smooth U-shaped curve, which is not possible assuming Weibull distributions. Flexible specifications, such as those based on piecewise functions or cubic splines functions,¹⁰ could be a good alternative to explore in future work.

The usage of the INLA is another strong point of our work. Computational time is highly reduced compared to MCMC methods, and the introduction of Gaussian random effects in the regression terms is something natural in the INLA. However, it is not a popular choice when assessing multi-state models yet, although this in the near future might change given its benefits. Multi-state modelling in INLA remains unexplored and further work regarding this issue would be of high interest.

Finally, assessing transition probabilities and cumulative incidences instead of only analysing the estimations of random effects on transition intensities provides a deeper understanding of the clinical or epidemiological problem. In fact, the risk assessment alone seldom provides predictive information regarding the prognostic of a patient given some specific characteristics. Meanwhile, trajectories expressed in terms of probabilities are dynamic and interpretable outcomes which are essential to reach an individualised care and that can be determinant to clinicians and policy makers.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802231172034 - Supplemental material for A Bayesian multivariate spatial approach for illness-death survival models

Supplemental material, sj-pdf-1-smm-10.1177_09622802231172034 for A Bayesian multivariate spatial approach for illness-death survival models by Fran Llopis-Cardona, Carmen Armero and Gabriel Sanfélix-Gimeno in Statistical Methods in Medical Research

Footnotes

Acknowledgements

We thank the editor and reviewers for their thorough comments and suggestions which have substantially improved our manuscript.

Declaration of conflicting interests

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

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been funded by Instituto de Salud Carlos III (ISCIII) through the projects [PI14/00993, PI18/01675], ‘RD16/0001/0011 – Red Temática de Servicios de Salud Orientados a Enfermedades Crónicas (REDISSEC)’, and ‘RD21/0016/0006 – Red de Investigación en Cronicidad, Atención Primaria y Promoción de la Salud (RICAPPS)’, and co-funded by the European Union. FLC was funded by Instituto de Salud Carlos III (ISCIII) [grant number FI19/00190], and co-funded by the European Union. CA was partially funded by Ministerio de Ciencia e Innovación (MCI, Spain) [grant number PID2019-106341GB-I00].

ORCID iD

Fran Llopis-Cardona

Supplemental material

All R code is available at the GitHub repository: . Supplementary file 1 contains the comparison between INLA and JAGS estimations and the assessment of the stability of the inference.

References

Andersen

Abildstrom

Rosthøj

. Competing risks as a multi-state model. Stat Methods Med Res 2002; 11: 203–215.

Kneib

Hennerfeind

. Bayesian semi parametric multi-state models. Stat Model 2008; 8: 169–198.

Le-Rademacher

Therneau

. The utility of multistate models: a flexible framework for time-to-event data. Curr Epidemiol Rep 2022; 9: 183–189.

Andersen

Keiding

. Multi-state models for event history analysis. Stat Methods Med Res 2002; 11: 91–115.

Vejakama

Ingsathit

McEvoy

, et al. Progression of chronic kidney disease: an illness-death model approach. BMC Nephrol 2017; 18: 205.

Armero

Cabras

Castellanos

, et al. Bayesian analysis of a disability model for lung cancer survival. Stat Methods Med Res 2016; 25: 336–351.

Kuhn

Olié

Grave

, et al. Estimating the future burden of myocardial infarction in France until 2035: an illness-death model-based approach. Clin Epidemiol 2022; 14: 255–264.

Cox

. Regression models and life-tables. J R Stat Soc Series B Stat Methodol 1972; 34: 87–220.

Christensen

Johnson

Branscum

, et al. Bayesian ideas and data analysis: an introduction for scientists and statisticians. 1st ed. CRC Press, 2010.

10.

Lázaro

Armero

Alvares

. Bayesian regularization for flexible baseline hazard functions in Cox survival models. Biom J 2021; 63: 7–26.

11.

Collett

. Modelling survival data in medical research. 3rd ed. New York: Chapman and Hall/CRC, 2014.

12.

Besag

. Spatial interaction and the statistical analysis of lattice systems. J R Stat Soc Series B Stat Methodol 1974; 36: 192–236.

13.

Besag

York

Molliè

. Bayesian image restoration with two applications in spatial statistics. Ann Inst Stat Math 1991; 43: 1–59.

14.

Leroux

Lei

Breslow

. Estimation of Disease Rates in Small Areas: A new Mixed Model for Spatial Dependence. In: Halloran ME and Berry D (eds) Statistical Models in Epidemiology, the Environment, and Clinical Trials. New York, NY: Springer, 2000, pp.179–191.

15.

Banerjee

Wall

Carlin

. Frailty modeling for spatially correlated survival data, with application to infant mortality in Minnesota. Biostatistics 2003; 4: 123–142.

16.

Carlin

Banerjee

. Hierarchical multivariate CAR models for spatio-temporally correlated survival data. In: Bernardo JM, Berger JO, Dawid AP, Smith AFM (eds). Bayesian Statistics 7. Oxford: Oxford University Press, 2003, pp.45–63.

17.

Nathoo

Dean

. Spatial multistate transitional models for longitudinal event data. Biometrics 2008; 64: 271–279.

18.

Meira-Machado

de Uña-Alvarez

Cadarso-Suárez

, et al. Multi-state models for the analysis of time-to-event data. Stat Methods Med Res 2009; 18: 195–222.

19.

Ibrahim

Chen

Sinha

. Bayesian survival analysis. 1th ed. New York, NY: Springer, 2001.

20.

Eberly

Carlin

. Identifiability and convergence issues for Markov chain Monte Carlo fitting of spatial models. Stat Med 2000; 19: 2279–2294.

21.

Rue

Martino

Chopin

. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J R Stat Soc Series B Stat Methodol 2009; 71: 319–392.

22.

Gilks

Richardson

Spiegelhalter

, (eds) Markov Chain Monte Carlo in Practice. 1st ed. New York: Chapman and Hall/CRC, 1995.

23.

Martino

Akerkar

Rue

. Approximate bayesian inference for survival models. Scand J Stat 2011; 38: 514–528.

24.

Van Niekerk

Bakka

Rue

. Competing risks joint models using R-INLA. Stat Model 2021; 21: 56–71.

25.

Van Niekerk

Bakka

Rue

. A principled distance-based prior for the shape of the Weibull model. Stat Probab Lett 2021; 174: 109098.

26.

Gelman

Carlin

Stern

, et al. Bayesian data analysis. 3rd ed. Boca Raton: Chapman and Hall/CRC, 2013.

27.

Schuurman

Grasman

RPPP

Hamaker

. A comparison of inverse-Wishart prior specifications for covariance matrices in multilevel autoregressive models. Multivar Behav Res 2016; 51: 185–206.

28.

O’Malley

Zaslavsky

. Domain-level covariance analysis for multilevel survey data with structured nonresponse. J Am Stat Assoc 2008; 103: 1405–1418.

29.

Lewandowski

Kurowicka

Joe

. Generating random correlation matrices based on vines and extended onion method. J Multivar Anal 2009; 100: 1989–2001.

30.

Gómez-Rubio

. Bayesian Inference with INLA. Boca Raton, FL: Chapman and Hall/CRC Press, 2020.

31.

Palmí-Perales

Gómez-Rubio

Martinez-Beneito

. Bayesian multivariate spatial models for lattice data with INLA. J Stat Softw 2020; 98: 1–29.

32.

Rue

Riebler

Sørbye

, et al. Bayesian computing with INLA: a review. Annu Rev Stat Appl 2017; 4: 395–421.

33.

Meira-Machado

Sestelo

. Estimation in the progressive illness-death model: a nonexhaustive review. Biom J 2019; 61: 245–263.

34.

Touraine

Helmer

Joly

. Predictions in an illness-death model. Stat Methods Med Res 2016; 25: 1452–1470.

35.

Chiuchiolo

van Niekerk

Rue

. Joint posterior inference for latent Gaussian models with R-INLA. J Stat Comput Simul 2022; 93: 723–752.

36.

Llopis-Cardona

Armero

Hurtado

, et al. Incidence of subsequent hip fracture and mortality in elderly patients: a multistate population-based cohort study in eastern Spain. J Bone Miner Res 2022; 37: 1200–1208.

37.

Ogle

Barber

. Ensuring identifiability in hierarchical mixed effects Bayesian models. Ecol Appl 2020; 30: 1–19.

38.

Gelman

Hill

. Data analysis using regression and Multilevel/Hierarchical models (Analytical Methods for Social Research). Cambridge: Cambridge University Press, 2006, p. 423.

39.

Berger

. An overview of robust Bayesian analysis. TEST 1994; 3: 5–124.

40.

Kruschke

. Doing Bayesian data analysis: a tutorial with R, JAGS, and Stan. London: Academic Press, 2014, p. 187.

41.

Klein

Moeschberger

. Survival analysis: techniques for censored and truncated data. New York: Springer, 2003, pp. 63–90.

42.

Zhang

Sun

. Interval censoring. Stat Methods Med Res 2010; 19: 53–70.

43.

Lamarca

Alonso

Gómez

, et al. Left-truncated data with age as time scale: an alternative for survival analysis in the elderly population. J Gerontol A Biol Sci Med Sci 1998; 53: M337–M343.

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

3.20 MB