Competing risks joint models using R-INLA

Abstract

The methodological advancements made in the field of joint models are numerous. None the less, the case of competing risks joint models has largely been neglected, especially from a practitioner's point of view. In the relevant works on competing risks joint models, the assumptions of a Gaussian linear longitudinal series and proportional cause-specific hazard functions, amongst others, have remained unchallenged. In this article, we provide a framework based on R-INLA to apply competing risks joint models in a unifying way such that non-Gaussian longitudinal data, spatial structures, times-dependent splines and various latent association structures, to mention a few, are all embraced in our approach. Our motivation stems from the SANAD trial which exhibits non-linear longitudinal trajectories and competing risks for failure of treatment. We also present a discrete competing risks joint model for longitudinal count data as well as a spatial competing risks joint model as specific examples.

Keywords

accelerated failure time model competing risks integrated nested Laplace approximation joint model non-Gaussian spline

1 Introduction

Time-to-event data is observed in its simplest form as the time when a particular event happens or when the monitoring process halts (censoring time). A competing risks model (Hoel, 1972; Moeschberger and David, 1971; Altshuler, 1970) arises in the case where patients are at risk of one of multiple events, that is, multiple mutually exclusive outcomes. The occurrence of an event precludes all other events so that each observational unit is at risk of each of the $C$ events until one event occurs. After this, the risk of all other events is zero. In classical survival analysis, the hazard and survival functions are estimated from the data to infer the risk of one specific event. In competing risks, the risk of each event is dependent on the assumption of surviving all other events and these classical methods cannot be used for each of the multiple events without the marginalization of all other events. There are mainly two modeling approaches in this respect, the subdistribution (Fine and Gray, 1999) and cause-specific hazards approaches (Prentice et al., 1978). In the subdistribution approach, the data enters the model through a predictor term in each of the cumulative incidence functions (CIF), the marginal cumulative probability of experiencing the specific event. On the other hand, the cause-specific hazards approach incorporates the data through a predictor term for each of the cause-specific hazard functions. Debates surrounding the choice between the two approaches are still ongoing with the consensus being that the cause-specific hazards approach should be chosen for causal inference of covariate effects. The CIF and the cause-specific hazard functions can be deduced from each other through certain functional relationships, see Crowder (2001) for more details.

In most situations, additional data is available in the form of repeated measurements over time constituting a longitudinal data set. The joint modeling of the competing risks data and the longitudinal data (Hu et al., 2009; Elashoff et al., 2007) are often necessary since the missing longitudinal observations are not missing at random but rather due to the occurrence of one of the competing events. Ignoring this phenomenon, will bias the longitudinal estimates. If the treatment efficacy is being evaluated, then joint modeling enables us to evaluate the treatment effect in both the longitudinal and time-event endpoints simultaneously. Also, unobserved heterogeneity in the individuals due to some latent effect over time can be successfully accommodated through the joint modeling.

Computational frameworks for competing risks have been developed in recent years such as the CFC package (Sharabiani and Mahani, 2017) and the cmprsk package (Gray, 2014). The existence of these frameworks is essential for catalysing the use of competing risks models in practice.

For joint modeling of competing risks and longitudinal data, the R package joineR (Philipson et al., 2018) can be used. This package implements the model developed by (Williamson et al., 2008). This model is a cause-specific Cox proportional hazards regression model (for two competing events) and a linear mixed effects model with an assumed latent Gaussian association process. An accelerated failure time (AFT) model, a non-Gaussian distribution for the longitudinal model component or spatial data thus excludes the use of joineR. Generic Bayesian software platforms such as BayesX, STAN and JAGS can also be used to make an inference on competing risks joint models through simulation, like Markov chain Monte Carlo (MCMC). In this article we will not use simulation-based tools but rather approximate methods which have considerable computational advantage.

Practitioners depend largely on the availability of wieldy computational frameworks to ease the computational burden of implementing such joint models. On the flip side, too many different developments could confuse the potential users and exacerbate the hesitation to apply competing risks models. This is especially true when more complex models are needed and new packages are developed for each additional complexity. This is cumbersome for users since many different packages should be familiarized for common use. The exigency for a more cohesive approach to competing risks joint models is clear. We believe that we address this deficiency in this article.

To achieve our aim of presenting a cohesive computational framework for competing risks joint models, we are spurred to characterize competing risks models within the class of latent Gaussian models. This characterization enables us to use the integrated nested Laplace approximation (INLA) method proposed by (Rue et al., 2009). Our rationale for this endeavour is the computational efficiency and ability to incorporate various modeling elements available in the accompanying R-INLA package (Rue et al., 2011). The INLA framework has been used in various fields of application since its inception. In survival analysis, (Martino et al., 2011) provided the details for basic survival models and (Álvaro-Meca et al., 2013) presented an approach for competing risks models using R-INLA, while (Van Niekerk et al., 2019) presented the case of joint survival longitudinal modeling. Competing risks joint models with frailties (Martino et al., 2011), spatial structures (Bakka et al., 2018) and non-linear covariate effects (Lindgren and Rue, 2008) with non-Gaussian longitudinal submodels are all natural possibilities within our framework.

In the SANAD trial (Marson et al., 2007), the aim was to evaluate the efficacy of various treatments for seizure control in patients with epilepsy. Two treatment failure types, unacceptable adverse effects (UAE) and inadequate seizure control (ISC) were considered as competing events. From exploratory data analysis, the non-linear mean trajectory of the titration dosage over time for the two different treatments is clear and also noted by (Williamson et al., 2008). A piecewise linear model with two different slopes (with the changepoint at time $500$ ) was used by (Williamson et al., 2008) as an attempt to properly model this phenomenon. We, however, propose a fully non-parametric spline component to capture this non-linear behavior over time in the form of a random walk order two model. This type of non-linear competing risks model can be formulated as a latent Gaussian model which can then be fitted using the established R-INLA framework for accurate and efficient computation. This endeavour resulted in better results for the SANAD trial as well as the development of a class of latent Gaussian competing risks joint models for various other realistic data situations such as longitudinal count data, spatial data and correlated competing risks to name a few.

In Section 2 we present an introduction and necessary details of competing risks joint models in general which are then formulated as latent Gaussian models in Section 3. We also provide two specific examples of competing risks joint models which might arise in practice to give the reader an intuition for the proposed methodology. We apply the proposed method to the SANAD trial in Section 4 and a simulation study for a discrete competing risks joint model (longitudinal count data) is presented in Section 5. Finally, we conclude with a discussion in Section 6 on the proposed method and some exciting possibilities thereof.

2 Joint models for competing risks and longitudinal data

The applications of joint models have increased over the last few years. Essentially, a joint model is a statistical method to model two or more dependent endpoints/responses simultaneously. An example of joint models in spatial statistics are coregionalization models (Schmidt and Gelfand, 2003; Krainski et al., 2018) where there is dependency between multiple spatial fields. In survival analysis, the use of joint models to infer joint survival longitudinal models is very popular. In these cases, joint modeling is fundamentally correct since the time-event process and the longitudinal data per individual is inherently dependent since it is governed by the same biological process. Although joint models are popular in survival analysis with one event of interest, the potential for their application in more complicated survival models is yet to be explored. In this article, we focus on competing risks joint models with various realistic features such as multiple events of interest, spatial random effects and non-linear longitudinal trajectories. To this end, we present cause-specific competing risks joint models with proportional or non-proportional cause-specific hazard functions and non-Gaussian or Gaussian longitudinal series as latent Gaussian models. This venture allows the use of the fast and accurate INLA methodology Rue et al. (2009) for Bayesian inference of the model.

2.1 Cause-specific competing risks joint model

The competing risks joint model is used to make an inference for the biological process that underpins the longitudinal and competing risks data. For simplicity we assume $C$ competing risks and one longitudinal marker in this section, even though our approach can be used for various longitudinal markers.

For the longitudinal marker $y$ , we formulate a linear predictor $η^{L}$ and vector of hyperparameters, $θ^{L}$ for $N$ individuals with a total of $N^{L} = N_{1} + N_{2} + . . . + N_{N}$ observations. Note that we do not make any distributional assumptions, hence denote the likelihood of $y$ as follows:

π_{L} (y | η^{L}, θ^{L}) = \prod_{l = 1}^{N^{L}} π_{L} (y_{l} | η_{l}^{L}, θ^{L}) .

(2.1)

Denote the cause-specific hazard function for cause $j, j = 1, . . ., C$ at time $t$ as follows:

h_{j} (t | η^{j}, θ^{j}),

where $η^{j}$ is the linear predictor. The vector of hyperparameters $θ^{j}$ is associated with the cause-specific hazard function, for example the baseline hazard parameters for the Cox proportional hazard model or the shape parameter of the Weibull accelerated failure time (AFT) model. The likelihood function for the observed event times $t = {t_{1}, . . ., t_{N}}$ with event indicators $d = {d_{ij}}, i = 1, . . ., N, j = 1, . . ., C$ is as follows:

π_{S} (t | η, θ) = \prod_{i = 1}^{N} \prod_{j = 1}^{C} h_{ij} (t_{i} | η_{i}^{j}, θ^{j})^{d_{ij}} exp (- \int_{0}^{t_{i}} h_{ij} (u | η_{i}^{j}, θ^{j}) du),

(2.2)

where $d_{ij}$ is an indicator if individual $i$ died from cause $j$ . The model is illustrated in Figure 1. The linear predictors, $η^{j}$ and $η^{L}$ have an additive structure of fixed and random effects as follows:

Figure 1:

Graphical illustration of the setup of the competing risks joint model

η = β Z + \sum_{k = 1}^{K} f_{k} (u_{ik}),

where $f_{k} (u_{ik})$ is the $k^{th}$ random effect of covariate $u_{ik}$ for individual $i$ . This formulation allows us to include spatial fields, spline components and frailties within our approach to joint competing risks models. The form of the linear predictor facilitates the contextualization of this joint competing risks model as a latent Gaussian model.

3 Competing risks joint models as latent Gaussian models

In this section we briefly present the details of latent Gaussian models (LGMs) and the integrated nested Laplace approximations (INLA) method for approximate Bayesian inference. If we can formulate competing risks joint models as LGMs and consequently use INLA for the inference, then we can fit these types of models very efficiently and with complex modeling components.

3.1 Latent Gaussian models and the integrated nested Laplace approximations method

Latent Gaussian models are a specific subset of hierarchical Bayesian additive models. This class comprises of well-known models such as mixed models, temporal and spatial models amongst many others. An LGM is defined as a model having a specific hierarchical structure, as follows:

The observed responses, ${y_{i}}$ , are conditionally independent based on the likelihood hyperparameters $θ_{1}$ and the linear predictors ${η_{i}}$ .

The linear predictor is formulated as follows:

η_{i} = β_{0} + β X + \sum_{k = 1}^{K} f_{k} (u_{ik}) + ε_{i},

(3.1)

where $β$ represent the linear fixed effects of the covariates $X$ , $ε$ is the unstructured random effects and the ${f_{k} (.)}$ ’s represent unknown (non-linear) functions of the covariates $u$ . The unknown functions ${f_{k} (.)}$ ’s, also known as structured random effects, include spline models, spatial effects, temporal effects, non-separable spatio-temporal effects, frailties, subject or group-specific intercepts and slopes etc.

The latent field $X$ is formed from the structured additive predictor as ${β_{0}, β, {f_{k} (.)}, η}$ with a multivariate Gaussian prior such that the latent field forms a Gaussian Markov random field with sparse precision matrix $Q (θ_{2})$ , that is, $X \sim N (0, Q^{- 1} (θ_{2}))$ .

A prior (possibly non-Gaussian), $π (θ)$ can then be formulated for the set of hyperparameters $θ = (θ_{1}, θ_{2})$ .

From this hierarchical Bayesian formulation the joint posterior distribution is then given as follows:

π (X, θ | y) \propto π (θ) π (X | θ) \prod_{i} π (y_{i} | X, θ) .

(3.2)

Within this framework the joint posterior density (3.2) and subsequently the marginal posterior densities, $π (X_{i} | y), i = 1 . . . n$ and $π (θ | y)$ can be efficiently and accurately calculated using the INLA methodology developed by Rue et al. (2009), without using Gibbs sampling.

3.2 Competing risks joint models as LGMs

In this section we show that most competing risks joint models are indeed LGMs and can therefore be defined as one model. This enables the use of the INLA methodology for many different competing risks joint models, for example discrete competing risks joint models, spatial competing risks joint models, non-linear competing risks joint models and correlated competing risks joint models.

Consider the case with $C$ competing risks, $N$ individuals with time to event observations $t$ and event variable $d$ , and $M$ series’ of longitudinal observations each of size $N_{l}, l = 1, . . ., M$ . Let $d_{i}, i = 1, . . ., C$ be a vector of indicators if the event type is $i$ .

Now define the stacked response list $y^{stacked} = {y_{1}, y_{2}, . . ., y_{M}, (t_{1}, d_{1}), (t_{2}, d_{2}), . . ., (t_{C}, d_{C})}$ . Also define the stacked linear predictor $η^{stacked} = {η_{1}^{L}, η_{2}^{L}, . . ., η_{M}^{L}, η^{1}, η^{2}, . . ., η^{C},}$ as in (2.1) and (2.2) where each component is of the form (3.1). For each component we define two vectors of hyperparameters, $θ_{1}$ for the hyperparameters from the likelihood (like a precision parameter if a Gaussian distribution is assumed or the shape parameter for a Weibull model etc.) and $θ_{2}$ for the hyperparameters from the random effects in the linear predictor (like a range parameter for a spatial effect or a precision parameter for a spline effect etc.), such that the stacked vector of hyperparameters is $θ^{stacked} = {θ_{1, 1}^{L}, θ_{2, 1}^{L}, θ_{1, 2}^{L}, θ_{2, 2}^{L}, . . ., θ_{1, M}^{L}, θ_{2, M}^{L}, θ_{1, 1}, θ_{2, 1}, θ_{1, 2}, θ_{2, 2}, . . ., θ_{1, C}, θ_{2, C}}$ . The latent field is then defined as $X = {β^{stacked}, η^{stacked}, {f (.)}^{stacked}}$ and we assume a multivariate Gaussian prior for $X$ .

From this construction the cause-specific competing risks joint model (in general, with various random effects such as spatial fields, spline effects and correlated models.) is an LGM and the joint posterior distribution of the latent field and the hyperparameters is as follows:

π (X, θ^{stacked}) \propto π (θ^{stacked}) π (X | θ^{stacked}) \prod_{i = 1}^{length (y^{stacked})} π (y_{i} | X, θ^{stacked}) .

Since this proposition is abstract, we provide some specific examples of relevant competing risks joint models in the next section.

3.3 Specific examples

Next, we present some specific (non-exhaustive) examples of competing risks joint models that fit in this framework.

3.3.1 Example 1. Poisson competing risks random-effects joint model

Suppose that the longitudinal observations, $y_{i} (t)$ follow a Poisson distribution with mean $λ_{i} (t)$ , i.e., $y_{i} (t) \sim Poisson (λ_{i} (t))$ with $λ_{i} (t) = exp (η_{i}^{L} (t))$ . Also assume three competing causes and a Weibull model for all the cause-specific hazard functions with shape parameters $α_{1}, α_{2}, α_{3}$ , that is, for observed time $t_{i}$ with event indicators $d_{i} = {d_{i 1}, d_{i 2}, d_{i 3}}$ , $h_{ij} (t) = {[α_{j} (γ_{ij})^{α_{j}} t^{α_{j} - 1}]}^{d_{ij}}$ with $γ_{ij} = exp (η_{i}^{(j)}), j = 1, 2, 3$ .

Now define the associated linear predictors as follows:

\begin{matrix} η^{L} (t) & = & β^{L} X^{L} (t) + v + w t \\ η^{(1)} (t) & = & β^{(1)} X^{(1)} (t) + γ^{(1)} v + κ^{(1)} w t \\ η^{(2)} (t) & = & β^{(2)} X^{(2)} (t) + γ^{(2)} v + κ^{(2)} w t \\ η^{(3)} (t) & = & β^{(3)} X^{(3)} (t) + γ^{(3)} v + κ^{(3)} w t, \end{matrix}

(3.3)

where $[\begin{matrix} v_{i} \\ w_{i} \end{matrix}] \sim N ([\begin{matrix} 0 \\ 0 \end{matrix}], [\begin{matrix} τ_{v}^{- 1} & ρ \\ ρ & τ_{w}^{- 1} \end{matrix}]) .$ Note that we do not enforce the same scaling parameter for the two random effects but accommodate for different scalings of the random intercept and the random slope. Now define $θ = {α_{1}, α_{2}, α_{3}, τ_{v}^{- 1}, τ_{w}^{- 1}, ρ}$ as the vector of hyperparameters. Define the latent field as follows

X = {β^{L}, β^{(1)}, β^{(2)}, β^{(3)}, η^{L}, η^{(1)}, η^{(2)}, η^{(3)}, v, w},

then this competing risks joint model is an LGM following Section 3.

3.3.2 Example 2. Spatial competing risks joint model with smoothed longitudinal predictor as the association term and unknown non-linear time trend

Suppose that the longitudinal observations, $y_{i}$ follow a Gaussian distribution with mean $μ_{i}$ and precision parameter $τ$ , i.e., $y_{i} \sim N (μ_{i}, τ^{- 1})$ with $μ_{i} = η_{i}^{L} (s, t)$ with $t$ a time index and $s$ a space index. Also assume three competing causes and an exponential model for all the cause-specific hazard functions that is for observed time $t_{i}$ with event indicators $d_{i} = {d_{i 1}, d_{i 2}, d_{i 3}}$ , $h_{ij} (t) = λ_{ij}^{d_{ij}}$ with $λ_{ij} = exp (η_{i}^{(j)}), j = 1, 2, 3$ .

Now define the associated linear predictors as follows:

\begin{matrix} η^{L} (s, t) & = & β^{L} X^{L} (t) + v + f (t) + u (s) \\ η^{(1)} (s, t) & = & β^{(1)} X^{(1)} (t) + γ^{(1)} η^{L} (s, t) \\ η^{(2)} (s, t) & = & β^{(2)} X^{(2)} (t) + γ^{(2)} η^{L} (s, t) \\ η^{(3)} (s, t) & = & β^{(3)} X^{(3)} (t) + γ^{(3)} η^{L} (s, t), \end{matrix}

(3.4)

where $f (t)$ is an unknown non-linear function of time (i.e., a spline component) with hyperparameters $θ_{time}$ and $u (s)$ are a spatial field (e.g., Besag (Besag, 1974) or BYM (Besag et al., 1991) models for discrete regions or SPDE (e.g., Matern) (Gelfand et al., 2010; Lindgren et al., 2011; Krainski et al., 2018) model for continuous space) with hyperparameters $θ_{space}$ . Now each cause-specific hazard function is associated with the longitudinal process through its own scaling $γ^{(j)}$ .

Now define $θ = {θ_{time}, θ_{space}, τ}$ as the vector of hyperparameters. Define the latent field as follows:

X = {β^{L}, β^{(1)}, β^{(2)}, β^{(3)}, η^{L}, η^{(1)}, η^{(2)}, η^{(3)}, f (.), u (.)},

then this spatial competing risks joint model is an LGM following Section 3.

4 Application: SANAD trial

The SANAD trial (Marson et al., 2007) was a randomized trial study to investigate the effect of seizure control drugs in patient with epilepsy for whom carbamazepine (CBZ) was considered as the standard treatment protocol. Two treatment failure types, UAE and ISC were considered as competing events. The patients were randomly assigned to CBZ, lamotrigine (LTG), gabapentin, oxcarbazepine and topiramate. In the main analysis, LTG appeared to be superior to CBZ but since then various criticisms regarding this finding have risen mainly due to the difference in titration speed and therefore the dosage. In Williamson et al. (2008) this data was analysed using a competing risks random-effects joint model with a changepoint at time $t = 500$ . In Figure 2 the mean trajectories for the two treatment groups CBZ and LTG are presented. It is clear that the mean trajectories are first different between the two treatments and second the changepoint at $t = 500$ seems quite subjective since for CBZ the trajectory does not seem to have a constant linear trend after $t = 500$ . We therefore include a spline component in the form of a random walk order two model for the dosage trajectory over time. This spline model is a stochastic spline, that is, converges to the true trajectory as the number of knots increase. We thus used 50 knots over the range of observed times. Also, we did not enforce equal scaling for the random intercepts and slopes.

Figure 2:

Mean trajectories for CBZ and LTG treatment groups

4.1 Bayesian competing risks joint model

We propose the following competing risks joint model with non-linear longitudinal trajectories for the SANAD trial data at time $t$ :

\begin{matrix} η^{L} (t) & = & β_{0} + β_{1} CBZ + CBZ f_{CBZ} (t) + LTG f_{LTG} (t) \\ η^{(UAE)} (t) & = & β_{0} + β_{1}^{(1)} CBZ + γ_{1} η^{L} (t) \\ η^{(ISC)} (t) & = & β_{0} + β_{1}^{(2)} CBZ + γ_{2} η^{L} (t), \end{matrix}

(4.1)

where $f_{CBZ}$ is a random walk order two model (Lindgren and Rue, 2008) with hyperparameter $τ_{f_{CBZ}}$ on the interaction of the time and $CBZ$ , the CBZ treatment indicator, and similarly for LTG. We assume a Gaussian longitudinal distribution with noise precision $τ_{L}$ , and Weibull models for the cause-specific hazard functions with hyperparameters $α_{UAE}, α_{ISC}$ . The priors for the latent field are all Gaussian which are as follows:

\begin{matrix} β_{0} \sim N (0, τ_{β_{0}}^{- 1}), β_{1} \sim N (0, τ_{β_{1}}^{- 1}), β_{1}^{(1)} \sim N (0, τ_{β_{1}^{(1)}}^{- 1}), \\ β_{1}^{(2)} \sim N (0, τ_{β_{1}^{(2)}}^{- 1}), γ_{1} \sim N (0, τ_{γ_{1}}^{- 1}), γ_{2} \sim N (0, τ_{γ_{2}}^{- 1}), \end{matrix}

with penalizing complexity priors (Simpson et al., 2017) for most of the hyperparameters, such that

\begin{matrix} σ_{L} = τ_{L}^{- 0.5} \sim Exp (λ_{L}), σ_{β_{0}} = τ_{β_{0}}^{- 0.5} \sim Exp (λ_{β_{0}}), σ_{β_{1}} = τ_{β_{1}}^{- 0.5} \sim Exp (λ_{β_{1}}), \\ σ_{β_{1}^{(1)}} = τ_{β_{1}^{(1)}}^{- 0.5} \sim Exp (λ_{β_{1}^{(1)}}), σ_{β_{1}^{(2)}} = τ_{β_{1}^{(2)}}^{- 0.5} \sim Exp (λ_{β_{1}^{(2)}}), σ_{γ_{1}} = τ_{γ_{1}}^{- 0.5} \sim Exp (λ_{γ_{1}}), \\ σ_{γ_{2}} = τ_{γ_{2}}^{- 0.5} \sim Exp (λ_{γ_{2}}), σ_{f_{CBZ}} = τ_{f_{CBZ}}^{- 0.5} \sim Exp (λ_{CBZ}), σ_{f_{LTG}} = τ_{f_{LTG}}^{- 0.5} \sim Exp (λ_{CBZ}), \\ 10 ln (α_{UAE}) \sim N (0, τ_{UAE}^{- 1}), 10 ln (α_{ISC}) \sim N (0, τ_{ISC}^{- 1}) . \end{matrix}

(4.3)

The user-defined parameter $λ_{[.]}$ is defined through the probability $P (τ_{[.]}^{- 0.5} > 1)$ , as an example $P (τ_{[.]}^{- 0.5} > 1) = 0.01$ produces $λ_{[.]} = ln (0.01)$ . Note that the exponential prior on the standard deviation $σ_{[.]} = τ_{[.]}^{- 0.5}$ translates to a Gumbel type 2 prior on $τ_{[.]}$ . The motivation for employing penalizing complexity priors for the precision hyperparameters are founded in the fact that the usual priors for the variance components, that is, independent inverse-gamma priors as in Huang et al. (2018), overfits in the sense that it cannot contract to the simpler model in which the respective model component has trivial variance. This is especially important in the case of joint models since the effect of overfitting is exacerbated by the influence of the shared random effect on all the linear predictors.

4.2 Results

The estimated mean dosage trajectories for the two treatements, CBZ and LTG, on a log scale are presented in Figure 3. The effect of the joint model on the longitudinal trajectories is clear when compared to Figure 2, in which the dropout data is not used. The dropout process thus affects the longitudinal trajectories in a meaningful way. The estimated competing risks joint model is as follows:

\begin{matrix} {\hat{η}}^{L} (t) & = & 2.438 (2.319; 2.556) + 0.305 (0.145; 0.465) CBZ + CBZ f_{CBZ} (t) + LTG f_{LTG} (t) \\ {\hat{η}}^{(UAE)} (t) & = & 2.438 (2.319; 2.556) - 0.726 (- 1.104; - 0.371) CBZ - 1.1 (- 1.192; - 1.02) η^{L} (t) \\ {\hat{η}}^{(ISC)} (t) & = & 2.438 (2.319; 2.556) - 0.11 (- 0.406; 0.175) CBZ - 0.998 (- 1.087; - 0.928) η^{L} (t), \end{matrix}

Figure 3:

Mean estimated dosage trajectories for the CBZ and LTG treatment groups in bold (with $95 %$ credible intervals)

with $τ_{f_{CBZ}} = 14.879$ , $τ_{f_{LTG}} = 7.055$ , $α_{UAE} = 0.84 (0.05; 0.915)$ and $α_{ISC} = 1.21 (1.19; 1.33)$ .

The negative association between the cause-specific hazard functions of UAE and ISC, respectively, and the smoothed dosage level should be interpreted with the shape parameters, $α_{UAE}$ and $α_{ISC}$ of the Weibull models. For UAE, the estimated shape parameter is significantly smaller than $1$ which implies that as the dose increases over time as captured by ${\hat{η}}^{L} (t)$ , the linear predictor ${\hat{η}}^{(UAE)} (t)$ will decrease due to the negative association and even more so for treatment CBZ. The decreased ${\hat{η}}^{(UAE)} (t)$ together with $α_{UAE} = 0.84$ leads to the conclusion that the hazard of UAE is lower for LTG than for CBZ and increasing over time since higher doses correspond to larger times. The converse is true for the hazard of ISC since $α_{ISC}$ is significantly larger than $1$ . Due to the insignificant effect of treatment on the hazard of ISC, the two treatments under consideration provide similar seizure control.

Our findings are thus consistent with those of (Williamson et al., 2008) with the added advantage of the flexible estimated mean longitudinal trajectories for the two treatment dosages.

5 Simulation study: Discrete competing risks random-effects joint model

In this section we simulate longitudinal count data as well as time-event data with three competing risks (Causes $1, 2$ and $3$ ). We also included a non-linear time component in the longitudinal trajectory. The data is simulated from a competing risks joint random-effects model with a random intercept (similar to (3.3)), such that for $N$ individuals, each with $N_{i}, i = 1, . . ., N$ longitudinal counts and time to event $t_{i}$ with event indicator $d_{i} \in {0, 1, 2, 3}$ the model is as folllows:

\begin{matrix} y_{i} (t) & \sim & Poisson (exp (t^{1.2} + u_{i})) \\ h_{1} (t) & = & exp (0.3 u_{i} + 0.01 Age) \\ h_{2} (t) & = & exp (- 0.1 u_{i} + 0.015 Age) \\ h_{3} (t) & = & exp (0.2 u_{i} + 0.0003 Age), \end{matrix}

with age a discrete explanatory variable within the range $(15; 75)$ and $u_{i} \sim N (0, 1)$ is a random intercept which associates the longitudinal series with the cause-specific hazard functions, $h_{j} (t), j = 1, 2, 3$ . Note that we also include a non-linear trend in the longitudinal trajectory, $t^{1.2}$ . In the simulation setup, we induce a negative association between the hazard of event $2$ and the longitudinal series, whereas the hazard of events $1$ and $3$ are positively correlated with the longitudinal series. From the above, the hazard functions are constant over time and an exponential distribution is thus implied.

5.1 Results

In Figure 4 the longitudinal series and the mean estimated trajectory from a random walk order two model with precision $τ_{r}$ are presented. The results are presented in Table 1 and it is clear that this challenging model is recovered well.

Figure 4:

Simulated longitudinal data and estimated mean trajectory

Table 1:

Estimates from the discrete competing risks joint model using R-INLA

Parameter	True value	Estimate	$95 %$ credible interval
$σ_{u}$	$1$	$1.017$	$(0.925; 1.114)$
$β_{1}$	$0.01$	$0.012$	$(0.005; 0.018)$
$β_{2}$	$0.015$	$0.052$	$(0.009; 0.019)$
$β_{3}$	$0.0003$	$0.0001$	$(- 0.003; 0.005)$
$γ_{1}$	$0.3$	$0.285$	$(0.192; 0.369)$
$γ_{2}$	$- 0.1$	$- 0.108$	$(- 0.206; - 0.011)$
$γ_{3}$	$0.2$	$0.208$	$(0.062; 0.352)$
$τ_{r}$	$-$	$28.513$	$(11.709; 116.878)$

From the estimates of the covariate effects in the cause-specific hazard functions we can deduce various other quantities of interest such as the cause-specific survival-like function and CIF etc. In Figure 5 the survival-like functions of the three risks are illustrated. Certain functions in the R-INLA framework can be used to estimate the marginal distributions of these quantities.

Figure 5:

Cause-specific survival-like functions of the three risks with $95 %$ credible interval for risks 1,2 and 3 (from top to bottom)

6 Discussion

Here we presented a general method to apply a wide class of competing risks joint models using the R-INLA framework. The main stimulus for this research was the necessity for a proper non-linear competing risks joint model as motivated from the SANAD trial. Although we addressed this issue, we also developed a class of latent Gaussian competing risks joint models which enable the practitioner to use our method for various realistic situations which demand different types of competing risks joint models. We presented the details of two specific realistic situations, one being a competing risks model with three competing events and one longitudinal series of count data, the other a spatial competing risks joint model. The application of the non-linear competing risks joint model to the SANAD trial proved fruitful as this model properly captures the non-linear trajectory.

The possibility of a non-linear competing risks joint model based on a random walk order two model emerged through the R-INLA framework. The INLA methodology (Rue et al., 2009) has been proved to be a very efficient and popular computational platform for Bayesian models. The use of the R-INLA framework necessitated the formulation of the proposed model as a latent Gaussian model, which in turn facilitated the formulation of various other competing risks joint models that could arise in practical situations, such as multiple longitudinal series, discrete longitudinal distributions, spatial models and frailty models to name a few. Quantities of interest such as the cause-specific survival function and the all-cause survival function, the median survival time can be computed in R-INLA.

The possibilities of realistic competing risks joint models that can be inferred through the proposed methodology are vast and exciting. In particular, spatial/spatio-temporal random effects can be included in the cause-specific hazard functions and/or the longitudinal submodel to develop, for example, a spatial competing risks joint model. The works of (Hickey et al., 2018) on the SANAD trial, (Huang et al., 2010) on outlying values in the longitudinal submodel, (Rajeswaran et al., 2018) on multivariate longitudinal data and the model presented by (Elashoff et al., 2008) are some examples in the literature that can be applied with our approach through R-INLA. A thorough presentation of all the possibilities is not possible in this article but we believe that the general method presented herein catalyses future research, especially different instances of competing risks joint models motivated from a modeling perspective.

Footnotes

Declaration of Conflicting Interests

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Appendix

The code for this article is available at http://www.statmod.org/smij/archive.html.

References

Altshuler

(1970). Theory for the measurement of competing risks in animal experiments. Mathematical Biosciences , 6, 1–11.

Álvaro-Meca

Akerkar

Alvarez-Bartolome

Gil-Prieto

Rue

de Miguel

ÁG

(2013). Factors involved in health-related transitions after curative resection for pancreatic cancer. 10-years experience: a multi state model. Cancer Epidemiology , 37, 91–96.

Bakka

Rue

Fuglstad

Riebler

Bolin

Illian

Krainski

Simpson

Lindgren

(2018). Spatial modeling with R-INLA: A review. Wiley Interdisciplinary Reviews: Computational Statistics , 10, e1443. doi: 10.1002/wics.1443.

Besag

(1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society, Series B (Methodological), 36, 192–225.

Besag

York

Mollié

(1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics , 43, 1–20.

Crowder

(2001). Classical Competing Risks . Boca Raton, FL: Chapman and Hall/CRC.

Elasho

(2007). An approach to joint analysis of longitudinal measurements and competing risks failure time data. Statistics in Medicine , 26, 2813–2835.

Elasho

(2008). A joint model for longitudinal measurements and survival data in the presence of multiple failure types. Biometrics , 64, 762–771.

Fine

Gray

(1999). A proportional hazards model for the subdistribution of a competing risk. Journal of the American Statistical Association , 94, 496–509.

10.

Gelfand

Diggle

Guttorp

Fuentes

(2010). Handbook of Spatial Statistics. Boca Raton, FL: CRC Press.

11.

Gray

(2014). Cmprsk: Subdistribution analysis of competing risks.URL https://CRAN.R-project.org/package=cmprskhttps://CRAN.R-project.org/package=cmprsk (last accessed 30 March 2020).

12.

Hickey

Philipson

Jorgensen

Kolamunnage-Dona

(2018). A comparison of joint models for longitudinal and competing risks data, with application to an epilepsy drug randomized controlled trial. Journal of the Royal Statistical Society, Series A (Statistics in Society), 181, 1105–1123.

13.

Hoel

(1972) A representation of mortality data by competing risks. Biometrics , 28, 475–488.

14.

(2009). A bayesian approach to joint analysis of longitudinal measurements and competing risks failure time data. Statistics in Medicine , 28, 1601–1619.

15.

Huang

Elasho

(2010). A joint model of longitudinal and competing risks survival data with heterogeneous random effects and outlying longitudinal measurements. Statistics and Its Interface , 3, 185.

16.

Huang

Chen

Liang

Zangmeister

(2018). Joint model-based clustering of nonlinear longitudinal trajectories and associated time-to-event data analysis, linked by latent class membership: With application to AIDS clinical studies. Lifetime Data Analysis , 24, 699–718.

17.

Krainski

Gomez-Rubio

Bakka

Lenzi

Castro-Camilo

Simpson

Lindgren

Rue

(2018). Advanced Spatial Modeling with Stochastic Partial Differential Equations using R and INLA . Boca Raton, FL: Chapman and Hall/CRC.

18.

Lindgren

Rue

(2008). On the secondorder model for irregular locations. Scandinavian Journal of Statistics , 35, 691–700.

19.

Lindgren

Rue

Lindström

(2011). An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. Journal of the Royal Statistical Society, Series B (Statistical Methodology), 73, 423–498.

20.

Marson

Al-Kharusi

Alwaidh

Appleton

Baker

Chadwick

Cramp

Cockerell

Cooper

Doughty

Eaton

(2007). Carbamazepine versus gabapentin, lamotrigine, oxcarbazepine or topiramate for partial epilepsy: Results from arm a of the SANAD trial. The Lancet , 369, 1000–1015.

21.

Martino

Akerkar

Rue

(2011). Approximate Bayesian inference for survival models. Scandinavian Journal of Statistics , 38, 514–528.

22.

Moeschberger

David

(1971). Life tests under competing causes of failure and the theory of competing risks. Biometrics , 27, 909–933.

23.

Philipson

Sousa

Diggle

Williamson

Kolamunnage-Dona

Henderson

Hickey

Sudell

(2018). JoineR: Joint modelling of repeated measure- ments and time-to-event data . URL https://CRAN.R-project.org/package=joineRhttps://CRAN.R-project.org/package=joineR (last accessed 18 March 2020).

24.

Prentice

Kalbeisch

Peterson

Jr Flournoy

Farewell

Breslow

(1978). The analysis of failure times in the presence of competing risks. Biometrics , 34, 541–554.

25.

Rajeswaran

Blackstone

Barnard

(2018). Joint modeling of multivariate longitudinal data and competing risks using multiphase sub-models. Statistics in Biosciences , 10, 651–685.

26.

Rue

Martino

Chopin

(2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society, Series B (Statistical Methodology), 71, 319–392 URL www.r-inla.org.

27.

Rue

Martino

Chopin

(2011). INLA: Approximate Bayesian inference using integrated nested Laplace approximations . URL www.r-inla.org.www.r-inla.org.

28.

Schmidt

Gelfand

(2003). A Bayesian coregionalization approach for multivariate pollutant data. Journal of Geophysical Research: Atmospheres , 108, STS10.1-STS10.8.

29.

Sharabiani

Mahani

(2017). CFC. Cause-specic framework for competing-risk analysis . URL https://CRAN.R-project.org/package=CFChttps://CRAN.R-project.org/package=CFC (last accessed 30 March 2020).

30.

Simpson

Rue

Riebler

Martins

Sørbye

(2017). Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science , 32, 1–28.

31.

Van Niekerk

Bakka

Rue

(2019). Joint models as latent gaussianmodels-not reinventing the wheel. arXiv:1901.09365

32.

Williamson

Kolamunnage-Dona

Philipson

Marson

(2008). Joint modelling of longitudinal and competing risks data. Statistics in Medicine , 27, 6426–6438.

Competing risks joint models using R-INLA

Abstract

Keywords

1 Introduction

2 Joint models for competing risks and longitudinal data

2.1 Cause-specific competing risks joint model

Graphical illustration of the setup of the competing risks joint model

3.1 Latent Gaussian models and the integrated nested Laplace approximations method

3.3 Specific examples

3.3.1 Example 1. Poisson competing risks random-effects joint model

Figure 2:

Mean trajectories for CBZ and LTG treatment groups

Figure 3:

Mean estimated dosage trajectories for the CBZ and LTG treatment groups in bold (with 95 % credible intervals)

5.1 Results

Figure 4:

Simulated longitudinal data and estimated mean trajectory

Cause-specific survival-like functions of the three risks with 95 % credible interval for risks 1,2 and 3 (from top to bottom)

Footnotes

Declaration of Conflicting Interests

Funding

Appendix

References

Mean estimated dosage trajectories for the CBZ and LTG treatment groups in bold (with $95 %$ credible intervals)

Cause-specific survival-like functions of the three risks with $95 %$ credible interval for risks 1,2 and 3 (from top to bottom)