# Bayesian estimation of two-part joint models for a longitudinal semicontinuous biomarker and a terminal event with INLA: Interests for cancer clinical trial evaluation

Denis Rustand, Janet van Niekerk, Havard Rue, Christophe Tournigand, Virginie Rondeau, Laurent Briollais

<figure><img src="image_1.png" /></figure>

## To cite this version:

Denis Rustand, Janet van Niekerk, Havard Rue, Christophe Tournigand, Virginie Rondeau, et al. Bayesian estimation of two-part joint models for a longitudinal semicontinuous biomarker and a terminal event with INLA: Interests for cancer clinical trial evaluation. Biometrical Journal, 2023, 65 (4), pp.2100322. 10.1002/bimj.202100322. hal-04036903

HAL Id: hal-04036903

https://hal.science/hal-04036903

Submitted on 20 Mar 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Bayesian estimation of two-part joint models for a longitudinal semicontinuous biomarker and a terminal event with INLA: Interests for cancer clinical trial evaluation

Denis Rustand<fnref n="1" />,<fnref n="2" /> ID

Janet van Niekerk<fnref n="2" /> ID

Håvard Rue<fnref n="2" /> |

Christophe Tournigand<fnref n="3" />

Virginie Rondeau<fnref n="1" />

Laurent Briollais<fnref n="4" />

## Correspondence

Denis Rustand, Biostatistic Team,
Bordeaux Population Health Center,
ISPED, Centre INSERM U1219, Bordeaux,
France.

Email: denis@rustand.fr

## Funding information

Government of Canada > Canadian Institutes of Health Research, Grant/Award Number: CIHR project grants; Government of Canada > Natural Sciences and Engineering Research Council of Canada, Grant/Award Number: NSERC Discovery grants; Institut National Du Cancer, Grant/Award Number: INCA EVALUATES OPE 2017-0680

<figure><img src="image_3.png" /></figure>

This article has earned an open data badge “Reproducible Research” for making publicly available the code necessary to reproduce the reported results. The results reported in this article could fully be reproduced.

## Abstract

Two-part joint models for a longitudinal semicontinuous biomarker and a terminal event have been recently introduced based on frequentist estimation. The biomarker distribution is decomposed into a probability of positive value and the expected value among positive values. Shared random effects can represent the association structure between the biomarker and the terminal event. The computational burden increases compared to standard joint models with a single regression model for the biomarker. In this context, the frequentist estimation implemented in the R package **frailtypack** can be challenging for complex models (i.e., a large number of parameters and dimension of the random effects). As an alternative, we propose a Bayesian estimation of two-part joint models based on the Integrated Nested Laplace Approximation (INLA) algorithm to alleviate the computational burden and fit more complex models. Our simulation studies confirm that **INLA** provides accurate approximation of posterior estimates and to reduced computation time and variability of estimates compared to **frailtypack** in the situations considered. We contrast the Bayesian and frequentist approaches in the analysis of two randomized cancer clinical trials (GERCOR and PRIME studies), where **INLA** has a reduced variability for the association between the biomarker and the risk of event. Moreover, the Bayesian approach was able to characterize subgroups of patients associated with different responses

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2023 The Authors. *Biometrical Journal* published by Wiley-VCH GmbH.

<footnote n="1">¹Biostatistic Team, Bordeaux Population Health Center, ISPED, Centre INSERM U1219, Bordeaux, France</footnote>

<footnote n="2">²Statistics Program, Computer, Electrical and Mathematical Sciences and Engineering Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Kingdom of Saudi Arabia</footnote>

<footnote n="3">³Hopital Henri Mondor, Creteil, France</footnote>

<footnote n="4">⁴Lunenfeld-Tanenbaum Research Institute, Mount Sinai Hospital, Dalla Lana School of Public Health (Biostatistics), University of Toronto, Toronto, Ontario, Canada</footnote>
to treatment in the PRIME study. Our study suggests that the Bayesian approach using the **INLA** algorithm enables to fit complex joint models that might be of interest in a wide range of clinical applications.

KEYWORDS

Bayesian estimation, computational efficiency, INLA, solid tumors cancer, two-part joint model

# 1 | INTRODUCTION

Estimation of joint models for longitudinal and time-to-event data was initially introduced using maximum likelihood estimation (Wulfsohn & Tsiatis, 1997; Henderson et al., 2000; Song et al., 2002; Chi & Ibrahim, 2006). It was further developed within the Bayesian framework in situations where maximum likelihood estimation with asymptotic assumptions faces nonidentifiability issues. It allows flexible and more complex association structures and can handle multiple longitudinal outcomes (Andrinopoulou & Rizopoulos, 2016). Bayesian joint models can be fitted with the R package **JMbayes** (Rizopoulos, 2016), which has been used in many biomedical researches (Lawrence Gould et al., 2015), among other packages (e.g., rstanarm; Muth et al., 2018). Bayesian estimation is usually based on Markov Chain Monte Carlo (MCMC) techniques (Hanson et al., 2011; R. Brown & G. Ibrahim, 2003; Xu & Zeger, 2001; Rizopoulos & Ghosh, 2011), which can have slow convergence properties. The Integrated Nested Laplace Approximation (INLA) algorithm has been recently introduced as an alternative to MCMC techniques for latent Gaussian models (LGM) (Rue et al., 2009; Martins et al., 2013). Many statistical models for spatial statistics, time series, etc., can be formulated as LGMs. A key feature of INLA is to provide approximations of the posterior marginals needed for Bayesian inference very efficiently and that still remain very accurate compared to MCMC methods (Rue et al., 2017). By formulating complex joint models as LGMs, INLA can be used to fit these models as developed recently (Van Niekerk et al., 2019; Van Niekerk, Bakka, Rue, & Schenk, 2021). For the two-part joint model (TPJM), INLA is yet to be used for inference.

TPJMs for a longitudinal semicontinuous biomarker and a terminal event have been recently introduced (Rustand et al., 2020) and applied to the joint analysis of survival times and repeated measurements of the sum of the longest diameter of target lesions (SLD), which is a biomarker representative of tumor burden in cancer clinical trials. The TPJM uses a conditional TPJM that decomposes the biomarker distribution into a binary outcome (zero vs. positive value) fitted with a logistic mixed effects model and a continuous outcome (positive values only) fitted with either a linear mixed effects model on the log-transformed outcome or a log-normal mixed effects model (Rustand et al., 2021). The “conditional” form of the two-part model gives the effect of covariates on the mean biomarker value conditional on a positive value in the continuous part. An alternative marginal model has recently been proposed to get the effect of covariates on the (unconditional) mean of the biomarker. A drawback of the marginal two-part model is that it may lead to arbitrary heterogeneity and provide less interpretable estimates on the conditional mean of the biomarker among positive values (Smith et al., 2014). In this article, we focus on the conditional TPJM, simply referred to as TPJM in what follows. The association with the survival model can be specified in terms of shared random effects, that is, random effects that are shared between the different components of the models including the binary and continuous parts of the model and the survival component. The TPJM is particularly interesting for cancer clinical trials evaluation because it can help in characterizing subgroups of patients who can benefit from a specific treatment (i.e., patients who had a complete tumor removal and no regrowth of the tumor would have a zero value of the biomarker and this is captured by the binary part of the model) and subgroups that do not respond or disbenefit from treatment, as it is sometimes observed when patients develop a resistance to chemotherapies (this will be captured by the continuous part of the model). These features cannot be modeled adequately by alternative strategies that assume a unique distribution for the zeros and the positive tumor size values (e.g., Tobit or Tweedie model; see Kurz, 2017). Besides, the effect of treatment on patients’ subgroups can be missed and may lead to biased conclusions as illustrated in Rustand et al. (2020).

An important limitation of the TPJM is the estimation procedure that requires a numerical approximation of the random effects distribution, which can lead to long computation times and convergence issues with high-dimensional parameter settings and complex association structures between the different components of the TPJMs. In this article, we propose an efficient Bayesian estimation procedure for the TPJM which relies on the INLA algorithm, as implemented in the R
package **INLA**. In practice, we used the R package **INLAjoint** that facilitates the fit of joint models with **INLA** and make it more user-friendly, as **INLA** was not developed specifically for this class of models. The Bayesian inference is compared to the frequentist estimation of the TPJM available in the R package **frailtypack** (Król et al., 2017). The remainder of the article is structured as follows: in Section 2, we describe the TPJM and introduce the frequentist and Bayesian estimations. In Section 3, we present a simulation study to assess the performance of these two estimation strategies in terms of bias, coverage probability, computation time, and convergence rate. An application to two randomized clinical trials each comparing two treatment strategies in patients with metastatic colorectal cancer is proposed in Section 4, and we conclude with a discussion in Section 5.

## 2 | ESTIMATION OF THE CONDITIONAL TWO-PART JOINT MODEL

### 2.1 | Model specification

Let $Y_{ij}$ denote the biomarker measurement for individual $i$ ($i = 1, \dots, n$) at visit $j$ ($j = 1, \dots, n_i$), $T_i$ denotes the survival time, and $\delta_i$ the censoring indicator for individual $i$. We use a logistic mixed effects model for the probability of a positive value of the biomarker and a linear mixed effects model for the conditional expected biomarker value. A proportional hazards survival model specifies the effect of covariates on survival time, adjusted for the individual heterogeneity captured in the biomarker model. The complete model is defined as follows:

$$ 
\left\{ 
\begin{array}{l} 
\eta_{Bij} = \text{Logit}(\text{Prob}(Y_{ij} > 0)) = \mathbf{X}_{Bij}^{\top} \boldsymbol{\alpha} + \mathbf{Z}_{Bij}^{\top} \mathbf{a}_i \quad \text{(Binarypart),} \\ 
\eta_{Cij} = \text{E}[\log(Y_{ij})|Y_{ij} > 0] = \mathbf{X}_{Cij}^{\top} \boldsymbol{\beta} + \mathbf{Z}_{Cij}^{\top} \mathbf{b}_i \quad \text{(Continuouspart),} \\ 
\lambda_i(t) = \lambda_0(t) \exp(\eta_{Si}) = \lambda_0(t) \exp(\mathbf{X}_i^{\top} \boldsymbol{\gamma} + \mathbf{a}_i^{\top} \boldsymbol{\varphi}_a + \mathbf{b}_i^{\top} \boldsymbol{\varphi}_b) \quad \text{(Survivalpart),} 
\end{array} 
\right. 
$$

where $\mathbf{X}_{Bij}$, $\mathbf{X}_{Cij}$, and $\mathbf{X}_i$ are vectors of covariates associated with the fixed effects $\boldsymbol{\alpha}$, $\boldsymbol{\beta}$, and $\boldsymbol{\gamma}$, respectively. Similarly, $\mathbf{Z}_{Bij}$ and $\mathbf{Z}_{Cij}$ are vectors of covariates associated with the random effects $\mathbf{a}_i$ and $\mathbf{b}_i$ in the binary and continuous parts. The random effects are shared in the survival model, with association parameters $\boldsymbol{\varphi}_a$ and $\boldsymbol{\varphi}_b$, respectively. These two vectors of random effects follow a multivariate normal distribution $N(0, \mathbf{Q}_{ab}^{-1})$, with covariance matrix $\mathbf{Q}_{ab}^{-1}$, defined as

$$ \mathbf{Q}_{ab}^{-1} = \begin{bmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ab} & \Sigma_{bb} \end{bmatrix}. $$

They account for both the association between the three components of the model and the correlation between the repeated measurements in the longitudinal process (observations are independent conditional on the random effects). We use a log-transformation of the biomarker to account for the positivity constraint and right-skewness in the continuous part of the model. The joint distribution assumes that the vectors of random effects underlie both the longitudinal and survival process, the joint distribution of the observed outcomes for individual $i$ is defined by

$$ 
\begin{aligned} 
p(T_i, \delta_i, Y_i | \mathbf{a}_i, \mathbf{b}_i; \mathbf{\Theta}) &= p(T_i, \delta_i | \mathbf{a}_i, \mathbf{b}_i; \mathbf{\Theta}) \prod_{j=1}^{n_i} p(Y_{ij} | \mathbf{a}_i, \mathbf{b}_i; \mathbf{\Theta}) \\ 
&= p(T_i, \delta_i | \mathbf{a}_i, \mathbf{b}_i; \mathbf{\Theta}) \prod_{j=1}^{n_i} p(Y_{ij} | Y_{ij} > 0; \mathbf{a}_i, \mathbf{b}_i; \mathbf{\Theta}) p(Y_{ij} > 0; \mathbf{a}_i, \mathbf{b}_i; \mathbf{\Theta}) 
\end{aligned} 
$$

with $\mathbf{\Theta}$ the full parameter vector, including the parameters for the binary, continuous, and survival outcomes, the baseline hazard function, and the random effects covariance matrix, such that the full conditional distribution is given by

$$ p(T, \delta, Y | \mathbf{a}, \mathbf{b}; \mathbf{\Theta}) = \prod_{i=1}^{n} p(T_i, \delta_i, Y_i | \mathbf{a}_i, \mathbf{b}_i; \mathbf{\Theta}). $$
The likelihood contribution for the ith subject can be formulated as follows:

$$
\begin{align*}
L_i(\Theta | \mathbf{Y}_i, T_i, \delta_i) &= \int_{\mathbf{a}_i} \int_{\mathbf{b}_i} \prod_{j=1}^{n_i} \exp \left( \mathbf{X}_{Bij}^{\top} \boldsymbol{\alpha} + \mathbf{Z}_{Bij}^{\top} \mathbf{a}_i \right)^{U_{ij}} \left( 1 - \frac{\exp(\mathbf{X}_{Bij}^{\top} \boldsymbol{\alpha} + \mathbf{Z}_{Bij}^{\top} \mathbf{a}_i)}{1 + \exp(\mathbf{X}_{Bij}^{\top} \boldsymbol{\alpha} + \mathbf{Z}_{Bij}^{\top} \mathbf{a}_i)} \right) \\
&\quad \times \left\{ \frac{1}{\sqrt{2\pi\sigma_{\varepsilon}^2}} \exp \left( - \frac{(\log(Y_{ij}) - \mathbf{X}_{Cij}^{\top}\boldsymbol{\beta} - \mathbf{Z}_{Cij}^{\top}\mathbf{b}_i)^2}{2\sigma_{\varepsilon}^2} \right) \right\}^{U_{ij}} \\
&\quad \times \lambda_i(T_i | \mathbf{a}_i, \mathbf{b}_i)^{\delta_i} \exp \left( - \int_0^{T_i} \lambda_i(t | \mathbf{a}_i, \mathbf{b}_i) dt \right) p(\mathbf{a}_i, \mathbf{b}_i) d\mathbf{b}_i d\mathbf{a}_i,
\end{align*}
$$

where $U_{ij} = I[Y_{ij} > 0]$, $\delta_i = I[T_i$ is not censored] and $\lambda_i(t|\mathbf{a}_i, \mathbf{b}_i) = \lambda_0(t|\mathbf{a}_i, \mathbf{b}_i) \exp\{\mathbf{X}_i(t)^\top \mathbf{\gamma} + \mathbf{a}_i^\top \boldsymbol{\varphi}_a + \mathbf{b}_i^\top \boldsymbol{\varphi}_b\}$.

## 2.2 | Bayesian estimation of the TPJM

Define $\mathbf{D} = \{T_i, \delta_i, Y_{ij} : i = 1, \dots, n; j = 1, \dots, n_i\}$ the observation variables. The goal of the Bayesian inference is to estimate the posterior distribution $\pi(\Theta|\mathbf{D})$. The joint posterior distribution $\pi(\Theta|\mathbf{D})$ is given by Bayes theorem as

$$ \pi(\Theta|\mathbf{D}) = \frac{p(\mathbf{D}|\Theta)\pi(\Theta)}{\pi(\mathbf{D})} \propto p(\mathbf{D}|\Theta)\pi(\Theta), $$

where $p(\mathbf{D}|\Theta)$ is the likelihood and $\pi(\Theta)$ is the joint prior. The marginal likelihood $\pi(\mathbf{D}) = \int_{\Theta} p(\mathbf{D}|\Theta)\pi(\Theta)d\Theta$ acts as a normalizing constant. The posterior marginal distribution of each parameter is then obtained by integrating out the other parameters of the model. In many cases, the posterior distribution is not analytically tractable and sampling-based methods like MCMC can be used. Approximate methods like INLA provide exact approximations to the posterior at a lower cost than sampling-based methods. The INLA methodology is based on the assumption that the statistical model is a latent Gaussian model, which we show in the next section for the TPJM.

## 2.3 | Formulation of the TPJM as a latent Gaussian model (Gaussian priors)

We now decompose $\Theta$ into the Gaussian latent field $\mathbf{u} = (\eta_B, \eta_C, \eta_S, \mathbf{a}, \mathbf{b}, \boldsymbol{\alpha}, \boldsymbol{\beta}, \gamma, \lambda, \varphi)$ and the set of hyperparameters $\Theta = (\Theta_1, \Theta_2)$. The hyperparameters $\Theta_1$ pertain to the latent field precision structure, while $\Theta_2$ contains the likelihood hyperparameters. The hyperparameters can follow any distribution and do not need to be Gaussian. Note that the first $\sum_{i=1}^n n_i + \sum_{i=1}^n n_i + n$ elements of $\mathbf{u}$ are the linear predictors of the TPJM and the rest of the elements are the latent unobserved variables. The linear predictors are included in the latent field so that each observation depends on the latent Gaussian field $\mathbf{u}$ only through one single element, which greatly simplifies the computations needed in the INLA algorithm (see Rue et al., 2009, 2017). We assume $\mathbf{a}_i, \mathbf{b}_i | Q_{ab} \sim \mathcal{N}(0, Q_{ab}^{-1})$, where $Q_{ab}$ is the precision matrix of the shared random effects with corresponding hyperparameters $\tau_{aa}, \tau_{ab}, \tau_{bb}$. We also assume $\boldsymbol{\alpha} \sim \mathcal{N}(0, \tau_{\alpha} \mathbf{I})$, $\boldsymbol{\beta} \sim \mathcal{N}(0, \tau_{\beta} \mathbf{I})$, $\gamma \sim \mathcal{N}(0, \tau_{\gamma} \mathbf{I})$, and $\varphi \sim \mathcal{N}(0, \tau_{\varphi} \mathbf{I})$. Let $\lambda$ denote a vector of coefficients associated with a random walk order one or order two used to approximate the log-baseline hazard function $\log(\lambda_0(t))$ of the survival model. These models are stochastic spline models with precision parameter $\tau_{\lambda}$. Thus, the latent field $\mathbf{u}$ is multivariate Gaussian with zero mean and precision matrix $Q(\Theta_1)$, that is,

$$ \mathbf{u}|\Theta_1 \sim \mathcal{N}(0, Q^{-1}(\Theta_1)). $$

Note that $Q(\Theta_1)$ is a sparse matrix indexed by a low dimension of parameters $\Theta_1$. This implies that the latent field $\mathbf{u}$ is a Gaussian Markov random field. The distribution of the observation variables $\mathbf{D}$ is denoted by $p(\mathbf{D}|\mathbf{u}, \Theta)$, and they are
conditionally independent given the Gaussian random field **u** and hyperparameters **θ** that is,

$$D|u, \theta \sim \prod_{i=1}^{n} p(d_i|u_i, \theta).$$

Then the posterior of $(u, \theta)$ can be written as

$$\pi(u, \theta|D) \propto \pi(\theta)\pi(u|\theta) \prod_{i=1}^{n} p(d_i|u_i, \theta), \\ \propto \pi(\theta)|Q(\theta_1)|^{n/2} \exp \left[ \frac{1}{2} u^\top Q(\theta_1) u + \sum_{i=1}^{n} \log\{p(d_i|u_i, \theta)\} \right].$$

This construction then shows that the TPJM is in fact an LGM since the latent field is a Gaussian Markov field and each data contribution depends on only one element of the latent field.

The main aim of INLA is then to approximate the posterior marginals $\pi(u_i|D), i = 1, \dots, n, \pi(\theta|D)$ and $\pi(\theta_m|D), m = 1, \dots, dim(\theta)$, as presented in the next section.

## 2.4 Integrated nested Laplace approximation

The INLA methodology introduced by Rue and Held (2005) is a major contribution to achieving efficient Bayesian inference, especially for complex or large models. INLA uses a unique combination of Laplace approximations and conditional distributions to approximate the joint posterior density as well as the marginals of the latent field and hyperparameters. It is thus not a sampling-based method like MCMC and such.

For the sake of brevity, the INLA methodology can be presented in the following three steps:

1. Approximate the marginal posterior distribution of hyperparameters using the Laplace approximation,

$$\pi(\theta|D) = \frac{\pi(u, \theta|D)}{\pi(u|\theta, D)} \approx \frac{\pi(\theta)\pi(u|\theta)p(D|u, \theta)}{\tilde{\pi}_G(u|\theta, D)}|_{u=u^*(\theta)},$$

where $\tilde{\pi}_G(u|\theta, D)$ is the Gaussian approximation of $\pi(u|\theta, D)$ at the mode $u^*(\theta)$ of the latent field for a given configuration of $\theta$. The marginal posterior $\pi(\theta_m|D)$ can be approximated by integrating $\theta_{-m}$ out in the previous approximation (while a good approximation of $\pi(u_i|\theta, D)$ is required to approximate the posterior marginal $\pi(u_i|D)$).

2. Approximate the conditional posterior distributions of the latent field,

$$\pi(u_i|\theta, D) \propto \frac{\pi(u, \theta|D)}{\pi(u_{-i}|u_i, \theta, D)},$$

using a Gaussian approximation (option 1), a Laplace approximation in a similar way as mentioned in step 1 (option 2), or using a “simplified Laplace approximation” (Rue et al., 2009), which corrects the Gaussian approximation for location and skewness by expanding the numerator and denominator up to a third-order Taylor series expansion in the Laplace approximation (option 3).

3. Use numerical integration to approximate the marginal posterior distributions of the latent field,

$$\pi(u_i|D) \approx \sum_{h=1}^{H} \tilde{\pi}(u_i|\theta_h^*, D) \tilde{\pi}(\theta_h^*|D) \Delta_h,$$

from steps 1 and 2. The integration points $\{\theta_1^*, \dots, \theta_H^*\}$ are selected from a rotation using polar coordinates and based on the density at these points, and $\Delta_h$ are the corresponding weights. The approximation of the posterior marginal for each element of the latent field and each hyperparameter, using numerical integration, forms the “integrated” part of the INLA algorithm while the first two steps above correspond to the “nested Laplace” approximation steps of INLA.
## 2.5 | Priors for the hyperparameters, $\theta$

From the formulation of the TPJM as an LGM, the prior for the hyperparameters, $\pi(\theta)$, should be specified. This prior can assume any form while keeping the TPJM an LGM. Amidst the debate about priors, Simpson et al. (2017) proposed a framework to construct principled priors for hyperparameters, namely penalizing complexity (PC) priors. These priors are derived based on the distance from a complex model to a simpler (base) model, with a user-defined parameter that informs the strength of contraction towards the simpler model. This parameter defines whether the PC priors are vague, weakly informative, or strongly informative based on the departure from the base model measured by the Kullback-Leibler distance. It is based on the principle of parsimony, simplifying the interpretation of the results by ensuring that the priors do not overfit. For example, the PC prior for the precision of the random walk model (baseline hazard) is derived to contract towards infinity, which in turn results in zero variance. The base model thus is a constant level model (constant baseline hazard, point mass at zero), while the complex model stems from finite precision as a random spline with zero mean.

In our case, we have various precision hyperparameters, {$\tau_{aa}, \tau_{ab}, \tau_{bb}, \tau_{\alpha}, \tau_{\beta}, \tau_{\gamma}, \tau_{\phi}, \tau_{\lambda}$}. We assign weakly informative priors to the fixed effects such that $\tau_{\alpha} = \tau_{\beta} = \tau_{\gamma} = \tau_{\phi} = 10^{-3}$. We thus need to formulate priors for the elements of $Q_{ab}$ and $\tau_{\lambda}$. For all these hyperparameters (precision and correlation parameters), we assume the respective PC priors as given in Simpson et al. (2017).

As an illustration, we give the details for the precision of the random intercept in the binary part of the TPJM (i.e., a) assuming the model proposed in the next section, $\tau_a = \Sigma_{aa}^{-1}$. The PC prior is derived as

$$ \pi(\tau_a) = \frac{\rho}{2} \tau_a^{-3/2} \exp(-\rho \tau_a^{-1/2}), $$

with the user-defined scaling parameter $\rho = -\frac{\ln(v)}{w}$. This parameter is chosen based on the desired tail behavior (or strength of contraction towards the base model $\sigma_a = \tau_a^{-1/2} = 0$) in the sense that $v$ and $w$ are such that

$$ P[\sigma_a > w] = v, \quad w > 0, 0 < v < 1. $$

Larger values of $v$ and $w$ result in higher prior density away from the base model, whereas smaller values of $v$ place more density closer to the base model. The same principle is used for specifying $\pi(\tau_{\lambda})$.

# 3 | SIMULATION STUDY

## 3.1 | Settings

We designed simulation studies to compare the performances of **INLA** and **frailtypack** in terms of bias of the parameter estimates, coverage probabilities, computation time, and convergence rates. The main factor driving the performance is the model complexity defined by the number of parameters. In particular, the number of correlated random effects defines the dimension of the integration that needs to be numerically approximated. We propose four simulation scenarios where the parameter values of the simulation models and proportion of zeros are based on the results obtained from the real data analyses. The first scenario includes a random intercept in the binary and continuous parts of the TPJM that are correlated. The second simulation scenario includes an additional random effect for the individual deviation from the mean slope in the continuous part, thus three correlated random effects. These two first simulation scenarios include 200 individuals in each dataset, corresponding to a small sample size commonly seen in randomized clinical trials while the third simulation scenario includes 500 individuals in each dataset. The last simulation scenario includes a sample size of 200 individuals and natural cubic splines for the biomarker specified with the *ns()* function in R, with a knot at the sample median of observed times. This generates three bases for the continuous part (one intercept and two slopes separated by the knot) and one basis for the binary part (intercept). The bases are then assumed correlated random effects in the model. For each scenario, we generate 1000 datasets; we first sample the positive longitudinal biomarker repeated measurements from a Gaussian distribution and include the zero values sampled from a binomial distribution. The association between the probability of zero value and the positive values is given by the correlated random effects. Survival times for the terminal event are generated from an exponential baseline hazard function with a scale of 0.2, and administrative censoring is
assumed to occur at the end of the follow-up (4 years). The rate of zeros is 8% (SD = 1%) for the first three scenarios and 5% for the last scenario, which is in between what we observed in our two real datasets (12% of zeros in application 1 and 4% in application 2). A zero value observation corresponds to a patient who experienced a complete disappearance of his/her target lesions and thus is extremely informative about treatment effect. The model for data generation is given by

$$ 
\left\{ 
\begin{array}{l} 
\text{Logit}[\text{Prob}(Y_{ij} > 0)] = \alpha_0 + a_i + \alpha_1 \cdot \text{time}_j + \alpha_2 \cdot \text{trt}_i + \alpha_3 \cdot \text{time}_j \cdot \text{trt}_i, \\ 
\text{E}[\log(Y_{ij})|Y_{ij} > 0] = \beta_0 + b_{0i} + (\beta_1 + b_{1i}) \cdot \text{time}_j + \beta_2 \cdot \text{trt}_i + \beta_3 \cdot \text{time}_j \cdot \text{trt}_i + \varepsilon_{ij}, \\ 
\lambda_i(t|Y_{ij}) = \lambda_0(t) \exp(\gamma \cdot \text{trt}_i + \varphi_a \cdot a_i + \varphi_{b_0} \cdot b_{0i} + \varphi_{b_1} \cdot b_{1i}) 
\end{array} 
\right. 
$$

$$ \begin{bmatrix} a_i \\ b_{0i} \\ b_{1i} \end{bmatrix} \sim \text{MVN} \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ab} & \Sigma_{bb} \end{bmatrix} \right) \equiv \text{MVN} \left( \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_a^2 & \sigma_{ab0} & \sigma_{ab1} \\ \sigma_{ab0} & \sigma_{b0}^2 & \sigma_{b0b1} \\ \sigma_{ab1} & \sigma_{b0b1} & \sigma_{b1}^2 \end{bmatrix} \right). $$

Note that in the first simulation scenario, we consider only a random intercept in the continuous part and the covariance matrix of the random effects is only defined by its first two lines and columns while in the second and third scenarios we use the full matrix. For the last scenario, the evolution of the biomarker over time in the continuous part is captured by two spline bases, each one has an interaction with treatment and is associated with a random effect. The baseline hazard function in the survival part of the model is approximated by a random walk model with **INLA** (Martino et al., 2011) such that for *m* bins on the time axis,

$$ \lambda_k - \lambda_{k-1} \sim N(0, \tau_{\lambda}), $$

where the PC prior (see Section 2.5) is used for the prior of $\tau_{\lambda}$.

The random walk order one model is a stochastic smoothing spline that smooths based on first-order differences. The number of bins is not influential (as opposed to knots of other splines) since an increase in the number of bins only results in an estimate closer to the stochastic model. In the simulations and applications, we use the random walk order two model that provides a smoother spline since the smoothing is then done in the second order. (See Van Niekerk, Bakka, & Rue, 2021, for more details on the use of these random walk models as Bayesian smoothing splines.) This approximation is different from **frailtypack** that uses cubic M-splines with five knots. A penalization ensures that the baseline hazard is smooth (a smoothing parameter is chosen using an approximate cross-validation criterion from a separate Cox model). The Levenberg-Marquardt algorithm, a robust Newton-like algorithm maximizes the log-likelihood function with **frailtypack** (Marquardt, 1963). The convergence of the algorithm depends on three conditions: The difference between the log-likelihood, the estimated coefficients, and the gradient of the log-likelihood of two consecutive iterations must be under $10^{-3}$. These convergence criteria avoid spurious convergence, making this algorithm more reliable than classical alternatives (e.g., Expectation-maximization or Broyden-Fletcher-Goldfarb-Shanno; see Philipps et al., 2021). We use a Monte Carlo approximation for the approximation of the integrals over the random effects in the likelihood function, with 1000 integration points, which is a reasonable trade-off between the precision of the approximation and computation time (using 2000 integration points doubled the computation time with negligible improvements of the results). The simulation studies are performed with 80 CPUs, **frailtypack** uses message passing interface (MPI) for parallel computation while the conjunction of **INLA** with the **PARDISO** library provides a high-performance computing environment with parallel computing support using OpenMP (Schenk & Gärtner, 2004). In practice, the 80 CPUs are mainly useful to reduce the computation time with **frailtypack** because the computation time with **INLA** is very low regardless of the number of threads because of the small sample size and number of hyperparameters.

## 3.2 | Results

In the results, we are comparing a Bayesian and a frequentist method, which each has usually its own evaluation criteria. Frequentist bias is used to evaluate the results from **frailtypack**, while the plausibility of the result based on 95% credible intervals is used to evaluate the results from **INLA** (Hespanhol et al., 2019). However, we are interested in the Bayesian approximation of the maximum likelihood estimation (i.e., noninformative priors) and therefore we provide an interpretation in this context (Walker, 1969). These 95% credible intervals are obtained by truncating the tails of the posterior
<figure><img src="image_4.png" /><figcaption>FIGURE 1 Simulations with two correlated random effects (the black dot is the mean value with 95% confidence intervals and the true value is the red vertical bar).</figcaption></figure>

distributions of each parameter with **INLA**, while 95% confidence with **frailtypack** is obtained with the inverse Hessian of the model, assuming Gaussian distribution for each parameter.

### 3.2.1 | Scenario 1: Two correlated random effects

The results are displayed in Figure 1 and detailed in Table S1 of the Supporting Information. The fixed effect parameters from the binary and continuous parts are properly estimated with both algorithms, with similar precision and coverage probabilities (CP) close to the expected 95% level. The parameter for the treatment effect in the survival part ($\gamma_1 = 0.2$) is associated with a larger variability with **frailtypack** ($\hat{\gamma}_1 = 0.22$, SD = 0.36, CP = 96%) compared to **INLA** ($\hat{\gamma}_1 = 0.19$, SD = 0.28, CP = 95%). The true value of the standard deviation of the random intercept in the binary part ($\hat{\sigma}_a = 1$) is within the 95% credible interval with **INLA** ($\hat{\sigma}_a = 0.93$, SD = 0.19, CP = 92%), with a slightly lower posterior mean value compared to **frailtypack**'s estimate ($\hat{\sigma}_a = 0.96$, SD = 0.21, CP = 94%). The random intercept's standard deviation in the continuous part is found similar to both algorithms with a slightly lower coverage probability of the true value with **frailtypack** (87%), but the correlation between the random intercepts of the binary and continuous parts ($\hat{\rho}_{ab} = 0.5$) has a reduced variability estimate with **INLA** ($\hat{\rho}_{ab} = 0.47$, SD = 0.11, CP = 94%) compared to **frailtypack** ($\hat{\rho}_{ab} = 0.49$, SD = 0.16, CP = 93%). The main difference observed is the estimation of the parameters for the association of the random effects with the risk of event, which links the biomarker to the terminal event. The association involving the random intercept from the binary part ($\hat{\phi}_a = 1$) has much lower variability with **INLA** ($\hat{\phi}_a = 0.99$, SD = 0.13, CP = 99%) and is unbiased with good coverage with **frailtypack** ($\hat{\phi}_a = 1.06$, SD = 0.89, CP = 94%). The association involving the random intercept from the continuous part ($\hat{\phi}_b = 1$) is biased upwards with **frailtypack** with large variability ($\hat{\phi}_b = 1.40$, SD = 1.28, CP = 92%), while **INLA**'s posterior estimate recovers the true value ($\hat{\phi}_b = 1.07$, SD = 0.15, CP = 98%). This could be due to the small sample size problems that cause more convergence issues under the frequentist framework. Although **INLA** yields accurate posterior estimates with small variability for these parameters, the coverage probabilities are higher than the expected 95%. The computation times are much lower with **INLA** (6 seconds per model, SD = 1) compared to **frailtypack** (74 seconds per model, SD = 29). Finally, 9% of the 1000 models did not reach convergence with **frailtypack**.

### 3.2.2 | Scenario 2: Three correlated random effects

The results are displayed in Figure 2 and detailed in Table S2 of the Supporting Information. With an additional random effect parameter compared to scenario 1, the fixed effects parameters are still properly estimated with **INLA**. The coverage
<figure><img src="image_5.png" /><figcaption>FIGURE 2 Simulations with three correlated random effects (the black dot is the mean value with 95% confidence intervals and the true value is the red vertical bar).</figcaption></figure>

probabilities are low with **frailtypack** for the slope and treatment by slope parameters in the continuous part ($\beta_1 = -0.3$ and $\beta_3 = 0.3$), while the parameter estimates remain unbiased. The variability for these two parameters is lower with **INLA** ($\hat{\beta}_1 = -0.30$, SD = 0.06, CP = 94% and $\hat{\beta}_3 = 0.30$, SD = 0.08, CP = 95%) compared to **frailtypack** ($\hat{\beta}_1 = -0.25$, SD = 0.12, CP = 44% and $\hat{\beta}_3 = 0.29$, SD = 0.15, CP = 47%). As observed in the first scenario, the treatment effect's posterior estimate in the survival model has lower variability with **INLA** ($\hat{\gamma}_1 = 0.19$, SD = 0.30, CP = 94%); moreover, the coverage probability for this parameter is lower than expected with **frailtypack** ($\hat{\gamma}_1 = 0.25$, SD = 0.53, CP = 82%). For the random effects covariance structure estimation, the posterior mean from **INLA** is slightly lower than the true value of the random intercept's standard deviation in the binary part ($\hat{\sigma}_a = 0.93$, SD = 0.15, CP = 95%) while **frailtypack**'s value is slightly higher than the true value ($\hat{\sigma}_a = 1.12$, SD = 0.29, CP = 91%). Overall, **INLA** has lower variability for the standard deviation and correlation terms and **frailtypack** has poor coverage for these parameters (e.g., 22% for $\sigma_{b1}$ and 23% for $\rho_{b0b1}$). The association parameters ($\varphi_a = 1$, $\varphi_{b0} = 1$, $\varphi_{b1} = 1$) are recovered well and have much lower variability with **INLA** ($\hat{\varphi}_a = 0.98$, SD = 0.12, CP = 99%; $\hat{\varphi}_{b0} = 1.10$, SD = 0.14, CP = 98%; $\hat{\varphi}_{b1} = 1.07$, SD = 0.14, CP = 98%) compared to **frailtypack** ($\hat{\varphi}_a = 0.89$, SD = 1.86, CP = 89%; $\hat{\varphi}_{b0} = 1.09$, SD = 1.82, CP = 89%; $\hat{\varphi}_{b1} = 1.46$, SD = 1.81, CP = 89%), but still with conservative coverage probabilities. Computation times remain much lower with **INLA** (6 seconds per model, SD = 1) compared to **frailtypack** (181 seconds per model, SD = 57) for which the time increased substantially when adding the third random effect. Moreover, the convergence rate of the model is reduced with **frailtypack** for this scenario (81%), because the model complexity increased.

### 3.2.3 | Scenario 3: Similar to Scenario 2 with increased sample size of $n = 500$

The results are displayed in Figure 3 and detailed in Table S3 of the Supporting Information. With $n = 500$ instead of $n = 200$, the standard deviations of the mean parameters are reduced overall. However, the coverage probabilities for the slope and treatment by slope parameters in the continuous part ($\beta_1 = -0.3$ and $\beta_3 = 0.3$) remain low with **frailtypack** ($\hat{\beta}_1 = -0.30$, SD = 0.08, CP = 45% and $\hat{\beta}_3 = 0.29$, SD = 0.10, CP = 46%), while the parameter estimates are unbiased. **INLA** has a lower variability and better coverage probabilities for these parameters ($\hat{\beta}_1 = -0.30$, SD = 0.04, CP = 92% and $\hat{\beta}_3 = 0.30$, SD = 0.05, CP = 95%) compared to **frailtypack**. For the random intercept's standard deviation in the binary part ($\sigma_a = 1$), the posterior mean from **INLA** is slightly closer to the true value ($\hat{\sigma}_a = 0.96$, SD = 0.11) compared
<figure><img src="image_6.png" /><figcaption>FIGURE 3 Simulations with *n* = 500 (the black dot is the mean value with 95% confidence intervals and the true value is the red vertical bar).</figcaption></figure>

to Scenario 2 while **frailtypack**'s estimate remains unbiased with a larger variability ($\hat{\sigma}_a$ = 1.01, SD = 0.14). The other parameters in the random effects covariance structure are well recovered with lower variability for the standard deviation and correlation terms overall with **INLA**. The association parameters ($\varphi_a = 1$, $\varphi_{b0} = 1$, $\varphi_{b1} = 1$) are recovered similarly as in Scenario 2 with **INLA** ($\hat{\varphi}_a = 0.91$, SD = 0.17, CP = 96%; $\hat{\varphi}_{b0} = 1.12$, SD = 0.18, CP = 94%; $\hat{\varphi}_{b1} = 1.09$, SD = 0.15, CP = 97%), while their variability and bias are reduced with **frailtypack** ($\hat{\varphi}_a = 0.99$, SD = 1.02, CP = 92%; $\hat{\varphi}_{b0} = 1.06$, SD = 0.96, CP = 89%, $\hat{\varphi}_{b1} = 1.08$, SD = 1.00, CP = 89%) but remains much higher compared to **INLA**. Note that the high coverage for the association parameters with **INLA** in the first two scenarios was likely explained by the small sample size, leading to a higher importance of the noninformative priors and resulting in high coverage. In this scenario with an increased sample size, the coverage is getting closer to the nominal 95% level. We observe an increased difference in computation times between **INLA** (11 seconds per model, SD = 1) and **frailtypack** (340 seconds per model, SD = 89) compared to scenario 2. Finally, the convergence rate of the model has improved with **frailtypack** under this scenario (96%).

### 3.2.4 | Scenario 4: Similar to Scenario 2 with two natural cubic splines for the trend in the continuous part (four correlated random effects, 200 individuals)

The last scenario illustrates the behavior of **INLA** with a more flexible trend in the continuous part. It includes four correlated random effects and 200 individuals. The results are displayed in Figure 4 and detailed in Table S4 of the Supporting Information. The estimation with **frailtypack** had increased convergence issues compared to previous scenarios, 46% of the 1000 models did not reach convergence. Moreover, fixed effects related to spline functions ($\beta_1 = -1$, $\beta_2 = -1$) have poor coverage and higher standard deviation with **frailtypack** ($\hat{\beta}_1 = -1.00$, SD = 0.13, CP = 88%; $\hat{\beta}_2 = -1.02$, SD = 0.16, CP = 85%) compared to **INLA** ($\hat{\beta}_1 = -0.99$, SD = 0.09, CP = 95%; $\hat{\beta}_2 = -1.00$, SD = 0.09, CP = 92%). Finally, the random effects variance and covariance parameters have higher variability and poorer coverage with **frailtypack** compared to **INLA**, in particular, the random intercept from the binary part ($\sigma_a = 0.5$) has a mean value of the true value ($\hat{\sigma}_a = 0.78$, SD = 0.25, CP = 89%) with **frailtypack**, while **INLA** is recovering this parameter with more accuracy ($\hat{\sigma}_a = 0.55$, SD = 0.10, CP = 98%). Note that this scenario is associated with a lower zero rate on average (5%), which may explain the poor result of **frailtypack** for this parameter.
<figure><img src="image_7.png" /><figcaption>FIGURE 4 Simulations with splines (the black dot is the mean value with 95% confidence intervals and the true value is the red vertical bar).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 5 Quantiles at 2.5%, 50%, and 97.5% of the baseline survival curves estimated with *frailtypack* and *INLA* in the simulations (Section 3).</figcaption></figure>

3.3 | Conclusions

Our method comparison suggests that the frequentist approach, implemented in **frailtypack**, reaches some limitations when fitting the more complex TPJMs, compared to the Bayesian approach implemented in **INLA**, which provides accurate approximations of posterior quantities. Convergence might be an issue, and estimation of the association parameters is highly variable with **frailtypack**. A representation of the baseline survival curves estimated under both scenarios is displayed in Figure 5. The median of the estimated survival curves is unbiased with both estimation strategies, but **INLA** has a reduced variability compared to **frailtypack**. Note that **INLA** is a deterministic algorithm, the notion of “convergence rate” does not apply and is only provided in the results for comparison with the iterative algorithm implemented
<table><caption>TABLE 1 Description of the GERCOR and PRIME study datasets.</caption><thead><tr><th rowspan="3">Study treatment</th><th colspan="2">GERCOR</th><th colspan="2">PRIME</th></tr><tr><th>Arm A</th><th>Arm B</th><th>Arm A</th><th>Arm B</th></tr><tr><th>FOLFIRI/FOLFOX6</th><th>FOLFOX6/FOLFIRI</th><th>FOLFOX4</th><th>Panitumumab/FOLFOX4</th></tr></thead><tbody><tr><td>Number of patients enrolled</td><td>109</td><td>111</td><td>593</td><td>590</td></tr><tr><td>Number of patients for the analysis</td><td>101</td><td>104</td><td>223</td><td>219</td></tr><tr><td>number of repeated measurements of the SLD</td><td>748</td><td>727</td><td>1192</td><td>1081</td></tr><tr><td>Number of zero values (%)</td><td>118 (16.2%)</td><td>56 (7.5%)</td><td>47 (3.8%)</td><td>52 (4.6%)</td></tr><tr><td>Number of death (%)</td><td>83 (82.2%)</td><td>82 (78.8%)</td><td>164 (73.5%)</td><td>164 (74.9%)</td></tr><tr><td>Median OS (years)</td><td>1.8 (1.4–2.3)</td><td>1.8 (1.5–2.2)</td><td>1.7 (1.5–1.9)</td><td>1.4 (1.3–1.7)</td></tr><tr><td>KRAS exon 2 at codons 12 and 13</td><td></td><td></td><td></td><td></td></tr><tr><td>Nonmutated</td><td></td><td></td><td>132 (59.2%)</td><td>128 (58.4%)</td></tr><tr><td>Mutated</td><td></td><td></td><td>91 (40.8%)</td><td>91 (41.6%)</td></tr><tr><td>Not available</td><td>101 (100%)</td><td>104 (100%)</td><td></td><td></td></tr></tbody></table>

in **frailtypack**. R codes including data simulation and model estimation with **frailtypack** and **INLA** assuming a conditional TPJM are available at github.com/DenisRustand/TPJM_sim/blob/master/TPJM_INLA.R and in the Supporting Information for Scenario 2, other scenarios follow trivially.

## 4 | APPLICATION

We applied the Bayesian TPJM to two cancer clinical trials, the GERCOR, and the PRIME studies. A comparison with **frailtypack** is provided only for the GERCOR data since this approach did not converge on the PRIME study. We used the same parameterization for **INLA** and **frailtypack**, as detailed in the simulation studies but with an increased number of integration points for the Monte Carlo method in **frailtypack** (i.e., 5000 points). In the context of a Bayesian approximation of the maximum likelihood estimation, we provide indications of the p-value for both **frailtypack** and **INLA** to ease the interpretation and the comparison of the results. These p-values are computed using the Z-score, assuming normal distributions. Each package has a specific criterion for model selection that takes into account the goodness of fit and the complexity of the model (i.e., the number of parameters). **Frailtypack** uses the likelihood-based cross-validation, LCV (Commenges et al., 2007), which accounts for its penalized likelihood while **INLA** is based on the deviance information criterion, DIC (Spiegelhalter et al., 2002).

### 4.1 | GERCOR study

#### 4.1.1 | Description

It is a randomized clinical trial investigating two treatment strategies that included a total of 220 patients with metastatic colorectal cancer. The reference strategy (arm A) corresponds to FOLFIRI (irinotecan) followed by FOLFOX6 (oxaliplatin), while arm B involves the reverse sequence. Patients were randomly assigned from December 1997 to September 1999, and the date chosen to assess overall survival was August 30, 2002. Complete data are available on 205 individuals for data analysis. Among them, 165 (80%) died during the follow-up. There are 1475 repeated measurements for the biomarker; 174 of which are zero values (12%). A summary of the dataset structure is given in Table 1. Our model uses death as the terminal event and the repeated SLD measurements (in centimeters) as the semicontinuous biomarker. Additional baseline covariates collected at the start of the study are also included, including performance status (0/1/2), lung metastatic site (Y/N), previous adjuvant radiotherapy (Y/N), previous surgery (no surgery/curative/palliative), and metastases (metachronous/synchronous). The first analysis of this dataset (Tournigand et al., 2004) did not find any significant difference between the two treatment strategies using classic survival analysis methods (i.e., log-rank tests). A trivariate joint model has been applied to this study for the simultaneous analysis of the longitudinal SLD, recurrent events (progres-
<figure><img src="image_9.png" /><figcaption>FIGURE 6 Mean biomarker value according to treatment received. The 95% credible intervals are obtained by resampling from the posterior parameter distributions.</figcaption></figure>

sion of lesions not included in the SLD or new lesions), and the terminal event (Król et al., 2018). A flexible mechanistic model using an ordinary differential equation was proposed to fit the biomarker dynamics. The results show a greater decline of the SLD for treatment arm A compared to treatment arm B. Moreover, the model finds a strong association between the biomarker model and the risk of terminal event. However, the interpretation of this treatment effect is difficult due to the nonlinear transformation applied to the outcome (Box-Cox) and the use of a nonlinear mechanistic model. Finally, a conditional TPJM was recently proposed (Rustand et al., 2020), which showed a significant treatment effect on the positive values of the biomarker (and no treatment effect on the probability of zero value). The model was able to show that when taking into account this treatment effect on the biomarker; the risk of terminal event is not significantly different between the two treatment arms. In the results, the mean parameters and their standard deviation are obtained by taking the maximum likelihood estimates and the inverse Hessian matrix with **frailtypack** while the posterior mean and standard deviation of the posterior distribution were used with **INLA**. For the evolution of the SLD over time in the continuous part, we compared a model with a linear trend (fixed + random slope) with a more flexible model with two natural cubic splines with a knot at the median of observed times which corresponds to 6 months of follow-up (each spline is associated with a fixed and a random effect). The goodness-of-fit criterion was in favor of the flexible model with both **INLA** (DIC = 2787 with a linear trend and DIC = 2260 with splines) and **frailtypack** (LCV = 0.45 with a linear trend and LCV = 0.40 with splines). The splines are described in Figure S1 of the Supporting Information and the results with a linear trend are given in Table S5 of the Supporting Information. The first spline corresponds to an initial increase during the first year followed by a decrease for the rest of the follow-up, while the second spline has an initial decrease during the first year and then increases for the rest of the follow-up.

4.1.2 | Results

As presented in Table 2, the fixed effect parameter estimates in the binary, continuous, and survival parts are quite similar between the frequentist and Bayesian approaches. The evolution of the SLD over time conditional on treatment is displayed in Figure 6, where random effects have been integrated out to have population average trajectories. Both treatment arms have a similar evolution of the log SLD over time. The uncertainty is slightly reduced with **INLA** at the end of the follow-up compared to **frailtypack**. Moreover, the mean value of the log SLD with **INLA** gets slightly lower over time compared to **frailtypack**, but no significant difference is observed marginally despite a significant association between the second spline and treatment with **INLA** ($\hat{\beta}_{11} = 1.49$, SD = 0.54), while **frailtypack** finds a lower effect size but still significant ($\hat{\beta}_{11} = 0.85$, SD = 0.39). This effect of treatment indicates that treatment arm B may be associated with a stronger reduction of the log SLD in early follow-up followed by a higher increase of the log SLD in the late follow-up compared to treatment arm A. This difference among positive values does not translate into a significant difference in the marginal evolution of the log SLD (i.e., including both zeros and positives), suggesting that models mixing zeros and positive (e.g., Tobit model) may miss this effect. Figure S2 in the Supporting Information displays observed vs. fitted longitudinal trajectory of the log SLD for five patients representative of the dataset, illustrating the good fit of the model. These fitted trajectories are only displayed with **INLA** since the linear predictors for each observation are directly available as part of the model's output, while **frailtypack** does not provide such information.
<table><caption>TABLE 2 Application of the Bayesian and frequentist TPJMs with two natural cubic splines to the GERCOR study with INLA and frailtypack.</caption><thead><tr><th>Approach</th><th></th><th>INLA<br>Est.<sup>a</sup> (SD<sup>b</sup>)</th><th>frailtypack<br>Est. (SD)</th></tr></thead><tbody><tr><td colspan="4"><strong>Binary part (SLD&gt;0 versus SLD = 0)</strong></td></tr><tr><td>Intercept</td><td>α<sub>0</sub></td><td>5.37*** (0.62)</td><td>5.72*** (0.71)</td></tr><tr><td>Time (year)</td><td>α<sub>1</sub></td><td>-2.01*** (0.37)</td><td>-2.02*** (0.40)</td></tr><tr><td>Treatment (B/A)</td><td>α<sub>2</sub></td><td>-0.98 (0.72)</td><td>-1.35 (0.69)</td></tr><tr><td>Performance status (1 vs. 0)</td><td>α<sub>3</sub></td><td>2.03*** (0.57)</td><td>2.17*** (0.57)</td></tr><tr><td>Performance status (2 vs. 0)</td><td>α<sub>4</sub></td><td>1.56 (1.12)</td><td>1.70 (1.14)</td></tr><tr><td>Previous� (Y/N)</td><td>α<sub>5</sub></td><td>0.71 (0.71)</td><td>0.59 (0.71)</td></tr><tr><td>Lung (Y/N)</td><td>α<sub>6</sub></td><td>1.81** (0.65)</td><td>1.54* (0.62)</td></tr><tr><td>Time: treatment (B/A)</td><td>α<sub>7</sub></td><td>0.31 (0.44)</td><td>0.36 (0.47)</td></tr><tr><td colspan="4"><strong>Continuous part (E[log(Y<sub>ij</sub>)|Y<sub>ij</sub> &gt; 0])</strong></td></tr><tr><td>Intercept</td><td>β<sub>0</sub></td><td>2.26*** (0.13)</td><td>2.15*** (0.09)</td></tr><tr><td>Slope 1</td><td>β<sub>1</sub></td><td>-0.26 (0.20)</td><td>0.06 (0.15)</td></tr><tr><td>Slope 2</td><td>β<sub>2</sub></td><td>1.48*** (0.37)</td><td>1.93*** (0.29)</td></tr><tr><td>Treatment (B/A)</td><td>β<sub>3</sub></td><td>-0.23** (0.09)</td><td>-0.23** (0.07)</td></tr><tr><td>Performance Status (1 vs. 0)</td><td>β<sub>4</sub></td><td>0.35*** (0.09)</td><td>0.41*** (0.07)</td></tr><tr><td>Performance status (2 vs. 0)</td><td>β<sub>5</sub></td><td>0.40** (0.14)</td><td>0.46*** (0.10)</td></tr><tr><td>Previous_surgery (curative)</td><td>β<sub>6</sub></td><td>-0.47** (0.16)</td><td>-0.26* (0.11)</td></tr><tr><td>Previous_surgery (palliative)</td><td>β<sub>7</sub></td><td>-0.03 (0.12)</td><td>0.10 (0.09)</td></tr><tr><td>Previous� (Y/N)</td><td>β<sub>8</sub></td><td>-0.22 (0.10)</td><td>-0.23 (0.08)</td></tr><tr><td>Metastases (metachronous vs. synchronous)</td><td>β<sub>9</sub></td><td>0.35* (0.14)</td><td>0.25* (0.10)</td></tr><tr><td>Slope 1:treatment (B/A)</td><td>β<sub>10</sub></td><td>0.91** (0.28)</td><td>0.53* (0.24)</td></tr><tr><td>Slope 2:treatment (B/A)</td><td>β<sub>11</sub></td><td>1.49*** (0.54)</td><td>0.85* (0.39)</td></tr><tr><td>Residual S.E.</td><td>σ<sub>ε</sub></td><td>0.25*** (0.01)</td><td>0.25*** (0.01)</td></tr><tr><td colspan="4"><strong>Death risk</strong></td></tr><tr><td>Treatment (B/A)</td><td>γ<sub>1</sub></td><td>0.25 (0.24)</td><td>0.07 (0.23)</td></tr><tr><td>Performance status (1 vs. 0)</td><td>γ<sub>2</sub></td><td>0.94*** (0.22)</td><td>1.05*** (0.24)</td></tr><tr><td>Performance status (2 vs. 0)</td><td>γ<sub>3</sub></td><td>1.66*** (0.36)</td><td>1.79*** (0.39)</td></tr><tr><td>revious_surgery (curative)</td><td>γ<sub>4</sub></td><td>-0.95* (0.44)</td><td>-0.70 (0.45)</td></tr><tr><td>Previous_surgery (palliative)</td><td>γ<sub>5</sub></td><td>-0.50 (0.31)</td><td>-0.39 (0.32)</td></tr><tr><td>Metastases (metachronous vs. synchronous)</td><td>γ<sub>6</sub></td><td>0.83* (0.36)</td><td>0.71 (0.37)</td></tr><tr><td colspan="4"><strong>Association</strong></td></tr><tr><td>Intercept (binary part)</td><td>φ<sub>a</sub></td><td>0.05 (0.08)</td><td>-0.07 (0.15)</td></tr><tr><td>Intercept (continuous part)</td><td>φ<sub>b0</sub></td><td>0.98*** (0.22)</td><td>1.26* (0.53)</td></tr><tr><td>Slope 1 (continuous part)</td><td>φ<sub>b1</sub></td><td>0.57*** (0.17)</td><td>0.76** (0.29)</td></tr><tr><td>Slope 2 (continuous part)</td><td>φ<sub>b2</sub></td><td>0.06 (0.07)</td><td>0.02 (0.08)</td></tr><tr><td colspan="4"><strong>Random effects' standard deviation</strong></td></tr><tr><td>Intercept (binary part)</td><td>σ<sub>a</sub></td><td>3.08*** (0.37)</td><td>3.12*** (0.39)</td></tr><tr><td>Intercept (continuous part)</td><td>σ<sub>b0</sub></td><td>0.59*** (0.04)</td><td>0.57*** (0.03)</td></tr><tr><td>Slope 1 (continuous part)</td><td>σ<sub>b1</sub></td><td>1.52*** (0.15)</td><td>1.73*** (0.12)</td></tr><tr><td>Slope 2 (continuous part)</td><td>σ<sub>b2</sub></td><td>2.50*** (0.29)</td><td>2.97*** (0.25)</td></tr><tr><td>ρ<sub>ab0</sub></td><td></td><td>0.51*** (0.08)</td><td>0.53*** (0.06)</td></tr><tr><td>ρ<sub>ab1</sub></td><td></td><td>0.59*** (0.12)</td><td>0.70*** (0.07)</td></tr><tr><td>ρ<sub>ab2</sub></td><td></td><td>0.32* (0.16)</td><td>0.55*** (0.09)</td></tr><tr><td>ρ<sub>b0b1</sub></td><td></td><td>-0.14 (0.10)</td><td>-0.03 (0.05)</td></tr><tr><td>ρ<sub>b0b2</sub></td><td></td><td>-0.09 (0.12)</td><td>0.08 (0.05)</td></tr><tr><td>ρ<sub>b1b2</sub></td><td></td><td>0.72 (0.08)**</td><td>0.78*** (0.03)</td></tr></tbody></table>

(Continues)
<table><caption>TABLE 2 (Continued)</caption><thead><tr><th>Approach</th><th>INLA<br>Est.<sup>a</sup> (SD<sup>b</sup>)</th><th>frailtypack<br>Est. (SD)</th></tr></thead><tbody><tr><td><strong>Computation time (Intel Xeon Gold 6248 2.50 GHz)</strong></td><td></td><td></td></tr><tr><td>Eight CPUs</td><td>13 s</td><td>11,230 s</td></tr><tr><td>Eighty CPUs</td><td>10 s</td><td>2355 s</td></tr></tbody></table>

Abbreviations: INLA, integrated nested Laplace approximation; SD, standard deviation; SLD, sum of the longest diameter of target lesions; TPJMs, two-part joint models.

aPosterior mean.

bPosterior standard deviation.

***p < 0.001.

**p < 0.01.

*p < 0.05.

The hazard ratio (HR) of treatment arm B versus treatment arm A that evaluates the change in the risk of death was higher with **INLA** (HR = 1.28, Confidence Interval, CI 0.80–2.06) compared to **frailtypack** (HR = 1.07, CI 0.68–1.69), but not significant in both cases. The main difference between **INLA** and **frailtypack** is in the estimation of the parameters for the association between the two-part model for the biomarker and the survival model. There is a positive and significant association between the random intercept ($\hat{\phi}_{b_0}$ = 0.98, SD = 0.22) and the first spline ($\hat{\phi}_{b_{s1}}$ = 0.57, SD = 0.17) from the continuous part and the risk of event with **INLA**. This association has a higher effect size and much larger variability with **frailtypack** ($\hat{\phi}_{b_0}$ = 1.26, SD = 0.52 and $\hat{\phi}_{b_1}$ = 0.76, SD = 0.29). This is in line with our simulation results (Scenario 1) where the association structure was estimated with better precision with **INLA**. The computation time for **frailtypack** increases quickly with the sample size and the model complexity (number of parameters and dimension of the random effects). The model was estimated in 11,230 s with **frailtypack** with eight CPUs and this reduces to 2355 s with 80 CPUs, while it is estimated in less than 15 s with **INLA**. The differences found in the association structure estimates are important when assessing the relationship between the biomarker dynamics and the risk of event. For instance, let us assume a clinician is interested in the top 15% patients who had the largest SLD at baseline compared to the average patient. Their random effect $b_{0i}$ should be higher than 1 standard deviation, that is from Table 2, $b_{0i} > 0.59$ with **INLA** (respectively, $b_{0i} > 0.57$ with **frailtypack**). Conditional on $b_{0i} > 0.59$ (respectively, $b_{0i} > 0.57$), the mean values of the random effects can be derived by sampling from a conditional multivariate normal distribution with a correlation matrix given in Table 2. These conditional means are 2.38, 0.90, -0.28, and -0.25 for a, $b_0$, $b_{s1}$, and $b_{s2}$, respectively (2.51, 0.87, -0.07, and 0.35 with **frailtypack**). Therefore, these top 15% individuals increase their chance to have the terminal event (i.e., to die) measured by a hazard ratio of HR = $exp(0.05 * 2.38 + 0.98 * 0.90 + 0.57 * (-0.28) + 0.06 * (-0.25)) = exp(0.83) = 2.29$, CI= 2.15–2.40, compared to a patient who has an average longitudinal SLD profile. **Frailtypack** gives a similar hazard ratio but with higher uncertainty HR = $exp((-0.07) * 2.51 + 1.26 * 0.87 + 0.76 * (-0.07) + 0.02 * 0.35) = exp(0.87) = 2.40$, CI=1.68–3.30. The confidence intervals were obtained by sampling parameters from the inverse Hessian matrix with **frailtypack** and the posterior distribution of the parameters with **INLA**. Figure 7 shows the baseline survival curves obtained with **frailtypack** and **INLA**, the survival decreases faster with **INLA** and has a reduced uncertainty compared to **frailtypack**, although no significant differences are observed between the two estimation strategies.

## 4.2 | PRIME study

### 4.2.1 | Description

The Panitumumab Randomized Trial in Combination with Chemotherapy for Metastatic Colorectal Cancer to Determine Efficacy (PRIME) study is a more challenging application for fitting the TPJM because it includes information about the KRAS mutation status (exon 2 codons 12/13), which has been shown to impact the clinical response to treatment in metastatic colorectal cancer patients (Van Cutsem et al., 2008; Normanno et al., 2009; Bokemeyer et al., 2008). It is, therefore, an important risk modifier, and clinicians are interested to assess treatment by mutation interaction in order to tailor treatment to patients’ genetic risk (Marabelle et al., 2020). This dataset is freely available on ProjectDataSphere.org (PDS UID: Colorec_Amgen_2006_309).
<figure><img src="image_10.png" /><figcaption>FIGURE 7 Baseline survival curves and their 95% confidence and credible intervals obtained from the application of the TPJM to the GERCOR study with frailtypack and INLA, respectively.</figcaption></figure>

The PRIME study is a randomized clinical trial that compares the efficacy and safety of panitumumab (anti-EGFR) in combination with FOLFOX4 (chemotherapy) with those of FOLFOX4 alone in the first-line treatment of patients, according to KRAS exon 2 status (wild type or mutant type). Between August 2006 and February 2008, 1183 patients were randomly assigned to receive treatment arm A (FOLFOX4 alone) or treatment arm B (panitumumab + FOLFOX4). The data for analysis include a subset of 442 patients (i.e., 741 excluded from the publicly available dataset). There are 2372 repeated measurements of the SLD, 99 of which are zero values (4%). The small rate of zero measurements in the SLD distribution leads to large variability in the binary part; however, zeros correspond to patients with a complete shrinkage of their target lesions, which is a very relevant information for clinicians about treatment effect. The number of individual repeated measurements for this biomarker varies between 1 and 24 with a median of 5. The death rate is 74%, corresponding to 328 deaths. Summary statistics of the dataset are given in Table 1. Additional baseline covariates collected at the start of the study are also included, including metastases to the liver at study entry (Y/N), the number of baseline metastases sites (1/2/3/4+), age (<60/60-69/≥70), and baseline Eastern Cooperative Oncology Group (ECOG) performance status (fully active/symptoms but ambulatory/in bed less that 50% of the time). We used a global backward selection procedure for each component of the model to select the covariates to include in the final joint model. The conclusions of the study are presented in Douillard et al. (2013) and show the importance of taking into account the mutation status when assessing treatment effect. Among patients without mutated KRAS, treatment arm B was associated with a slightly significantly reduced risk of death compared to treatment arm A. For patients with mutated KRAS, treatment arm B was associated with a nonsignificant increase in the risk of death compared to treatment arm A. Unlike the first application, we assume a linear trend on the log scale for the continuous part of the TPJM because the model with spline functions was not fitting well the data (high variability of the posterior parameter estimates that reflects the noninformative priors), suggesting that the data is not informative enough to fit all the parameters properly (small sample size, short follow-up due to high death rate, and many added parameters due to interactions between splines, treatment, and KRAS mutation status).

### 4.2.2 | Results

As presented in Table 3, in the binary part of the TPJM, the intercept is very large ($\hat{\alpha}_0 = 16.50$, SD = 3.49), corresponding to a high probability of positive value at baseline. This probability is increased for patients with mutated KRAS ($\hat{\alpha}_3 = 6.18$, SD = 4.54) and patients receiving treatment arm B ($\hat{\alpha}_2 = 3.15$, SD = 3.68) but with large standard deviations so that these effects are not significant. The slope parameter with time is negative and significant ($\hat{\alpha}_1 = -8.31$, SD = 1.66), meaning that patients without mutated KRAS and receiving treatment A have a higher odds of zero SLD value over time, that is, complete response to treatment. This odds decreases, but not significantly, among patients with either mutated KRAS ($\hat{\alpha}_6 = 0.64$, SD = 2.36) or receiving treatment B ($\hat{\alpha}_5 = 1.27$, SD = 1.81) and in patients with both mutated KRAS and receiving treatment arm B ($\hat{\alpha}_7 = 0.73$, SD = 4.42).

In the continuous part of the TPJM, patients with the wild-type KRAS status and in treatment arm A are associated with a decrease in the SLD value over time conditional on a positive SLD value ($\hat{\beta}_1 = -1.64$, SD = 0.10). This reduction of
<table><caption>TABLE 3 Application of the Bayesian TPJM with shared random effects to the PRIME study with the R package R-INLA.</caption><thead><tr><th>Approach</th><th></th><th>R-INLA<br>Est.<sup>a</sup> (SD<sup>b</sup>)</th></tr></thead><tbody><tr><td colspan="3"><strong>Binary part (SLD&gt;0 versus SLD=0)</strong></td></tr><tr><td>Intercept</td><td>α<sub>0</sub></td><td>16.50<sup>***</sup> (3.49)</td></tr><tr><td>Time (year)</td><td>α<sub>1</sub></td><td>-8.31<sup>***</sup> (1.66)</td></tr><tr><td>Treatment (B/A)</td><td>α<sub>2</sub></td><td>3.15 (3.68)</td></tr><tr><td>KRAS (MT/WT)</td><td>α<sub>3</sub></td><td>6.18 (4.54)</td></tr><tr><td>Treatment (B/A):kras (MT/WT)</td><td>α<sub>4</sub></td><td>-1.26 (6.73)</td></tr><tr><td>Time: treatment (B/A)</td><td>α<sub>5</sub></td><td>1.27 (1.81)</td></tr><tr><td>Time: KRAS (MT/WT)</td><td>α<sub>6</sub></td><td>0.64 (2.36)</td></tr><tr><td>Time: treatment (B/A):kras (MT/WT)</td><td>α<sub>7</sub></td><td>0.73 (4.42)</td></tr><tr><td colspan="3"><strong>Continuous part E[log(Y<sub>ij</sub>)|Y<sub>ij</sub> &gt; 0]</strong></td></tr><tr><td>Intercept</td><td>β<sub>0</sub></td><td>2.55<sup>***</sup> (0.17)</td></tr><tr><td>Time (years)</td><td>β<sub>1</sub></td><td>-1.64<sup>***</sup> (0.10)</td></tr><tr><td>Treatment (B/A)</td><td>β<sub>2</sub></td><td>-0.22<sup>*</sup> (0.10)</td></tr><tr><td>KRAS (MT/WT)</td><td>β<sub>3</sub></td><td>-0.24<sup>*</sup> (0.11)</td></tr><tr><td>Liver metastases (Y/N)</td><td>β<sub>4</sub></td><td>0.56<sup>***</sup> (0.13)</td></tr><tr><td>Eastern Cooperative Oncology Group performance status (symptoms but ambulatory vs. fully active)</td><td>β<sub>5</sub></td><td>0.18<sup>*</sup> (0.07)</td></tr><tr><td>Eastern Cooperative Oncology Group performance status (in bed less than 50% of the time vs. fully active)</td><td>β<sub>6</sub></td><td>0.49<sup>**</sup> (0.17)</td></tr><tr><td>Baseline metastases sites (2 vs. 1)</td><td>β<sub>7</sub></td><td>0.08 (0.11)</td></tr><tr><td>Baseline metastases sites (3 vs. 1)</td><td>β<sub>8</sub></td><td>0.23<sup>*</sup> (0.11)</td></tr><tr><td>Baseline metastases sites (4+ vs. 1)</td><td>β<sub>9</sub></td><td>0.18 (0.12)</td></tr><tr><td>Treatment (B/A): KRAS (MT/WT)</td><td>β<sub>10</sub></td><td>0.23 (0.15)</td></tr><tr><td>Time: treatment (B/A)</td><td>β<sub>11</sub></td><td>0.93<sup>***</sup> (0.14)</td></tr><tr><td>Time: KRAS (MT/WT)</td><td>β<sub>12</sub></td><td>1.18<sup>***</sup> (0.15)</td></tr><tr><td>Time: treatment (B/A):kras (MT/WT)</td><td>β<sub>13</sub></td><td>-1.12<sup>***</sup> (0.20)</td></tr><tr><td>Residual S.E.</td><td>σ<sub>ε</sub></td><td>0.27<sup>***</sup> (0.01)</td></tr><tr><td colspan="3"><strong>Death risk</strong></td></tr><tr><td>Treatment (B/A)</td><td>γ<sub>1</sub></td><td>0.08 (0.17)</td></tr><tr><td>KRAS (MT/WT)</td><td>γ<sub>2</sub></td><td>0.18 (0.18)</td></tr><tr><td>Treatment (B/A):kras (MT/WT)</td><td>γ<sub>3</sub></td><td>0.07 (0.24)</td></tr><tr><td>Age (60-69 vs. &lt;60)</td><td>γ<sub>4</sub></td><td>0.10 (0.13)</td></tr><tr><td>Age (70+ vs. &lt;60)</td><td>γ<sub>5</sub></td><td>0.25 (0.14)</td></tr><tr><td>Liver metastases (Y/N)</td><td>γ<sub>6</sub></td><td>0.05 (0.23)</td></tr><tr><td>Eastern Cooperative Oncology Group performance status (symptoms but ambulatory vs. fully active)</td><td>γ<sub>7</sub></td><td>0.31<sup>*</sup> (0.12)</td></tr><tr><td>Eastern Cooperative Oncology Group performance status (in bed less than 50% of the time vs. fully active)</td><td>γ<sub>8</sub></td><td>0.85<sup>**</sup> (0.26)</td></tr><tr><td>Baseline metastases sites (2 vs. 1)</td><td>γ<sub>9</sub></td><td>0.13 (0.20)</td></tr><tr><td>Baseline metastases sites (3 vs. 1)</td><td>γ<sub>10</sub></td><td>0.33 (0.20)</td></tr><tr><td>Baseline metastases sites (4+ vs. 1)</td><td>γ<sub>11</sub></td><td>0.44<sup>*</sup> (0.21)</td></tr><tr><td colspan="3"><strong>Association</strong></td></tr><tr><td>Intercept (binary part)</td><td>φ<sub>a</sub></td><td>0.00 (0.01)</td></tr><tr><td>Intercept (continuous part)</td><td>φ<sub>b0</sub></td><td>0.47<sup>***</sup> (0.10)</td></tr><tr><td>Slope (continuous part)</td><td>φ<sub>b1</sub></td><td>0.10 (0.14)</td></tr></tbody></table>

(Continues)
<table><caption>TABLE 3 (Continued)</caption><thead><tr><th>Approach</th><th></th><th>R-INLA<br>Est.<sup>a</sup> (SD<sup>b</sup>)</th></tr></thead><tbody><tr><td colspan="3"><strong>Random effects' standard deviation</strong></td></tr><tr><td>Intercept (binary part)</td><td>σ<sub>a</sub></td><td>8.61</td></tr><tr><td>Intercept (continuous part)</td><td>σ<sub>b<sub>0</sub></sub></td><td>0.73</td></tr><tr><td rowspan="3">Slope (continuous part)</td><td>σ<sub>b<sub>1</sub></sub></td><td>0.73</td></tr><tr><td>ρ<sub>ab<sub>0</sub></sub></td><td>0.00</td></tr><tr><td>ρ<sub>ab<sub>1</sub></sub></td><td>0.76</td></tr><tr><td rowspan="3"><strong>Computation time (Intel Xeon E5-4627 v4 2.60 GHz)</strong></td><td>ρ<sub>b<sub>0</sub>b<sub>1</sub></sub></td><td>-0.23</td></tr><tr><td>Eight CPUs</td><td>46 s</td></tr><tr><td>Eighty CPUs</td><td>39 s</td></tr></tbody></table>

Abbreviations: INLA, integrated nested Laplace approximation; MT, mutant type; SD, standard deviation; SLD, sum of the longest diameter of target lesions; TPJMs, two-part joint models; WT, wild type.

aPosterior mean.

bPosterior standard deviation.

***p < 0.001.

**p < 0.01.

*p < 0.05.

SLD over time is attenuated in patients with mutated KRAS ($\hat{\beta}_{12}$ = 1.18, SD = 0.15) or receiving treatment B ($\hat{\beta}_{11}$ = 0.93, SD = 0.14). Patients with the KRAS mutation and who received treatment B have a similar SLD trend over time as patients with the KRAS mutation who received treatment A or patients who received treatment B but with the wild-type KRAS status because of the negative interaction term between time, treatment, and KRAS status ($\hat{\beta}_{13}$ = -1.12, SD = 0.20).

In the survival part, the model shows no significant difference between treatment arms for the risk of death ($\hat{\gamma}_1$ = 0.08, SD = 0.17). Besides, patients with mutated KRAS have a similar risk of death compared to patients with the wild type ($\hat{\gamma}_2$ = 0.18, SD = 0.18), so do patients with mutated KRAS receiving treatment B ($\hat{\gamma}_3$ = 0.07, SD = 0.24). The random effect from the binary part and the random slope from the continuous part is not associated with the risk of death ($\hat{\phi}_a$ = 0.00, SD = 0.01 and $\hat{\phi}_{b1}$ = 0.10, SD = 0.14), but the random intercept from the continuous part ($\hat{\phi}_{b0}$ = 0.47, SD = 0.10) has a positive and highly significant association with the risk of event. This means that conditional on a positive value, the individual deviation from the mean baseline value of the SLD is predictive of the risk of event. Similarly to the GERCOR study, we can compare the top 15% patients with the smallest SLD at baseline to the average patient, their risk of death is reduced by 39% (HR = 0.61, CI = 0.51-0.75).

In conclusion, we did not find a direct effect of treatment B versus A on the risk of death while the initial study (Douil-lard et al., 2013) finds a slightly significant improvement in overall survival for patients with wild-type KRAS status (HR = 0.78, CI = 0.62-0.99), likely because of the reduced sample size available for our analysis (publicly available dataset only includes 37% of the original set of patients). Interestingly, the analysis of the continuous part of the TPJM suggests that the reduction of the SLD over time conditional on a positive value is attenuated with treatment B compared to treatment A for patients with wild-type KRAS status. A graphical representation of the mean biomarker evolution over time according to KRAS mutation status and treatment received is depicted in Figure 8. It confirms the suggested significant difference between treatment arms for patients with wild-type KRAS status and shows no treatment effect for patients with mutant KRAS.

5 | DISCUSSION

In this article, we developed a Bayesian estimation approach based on the **INLA** algorithm for TPJMs for a longitudinal semicontinuous biomarker and a terminal event. We also provided a comparison with a frequentist alternative approach implemented into the **frailtypack** package, using small sample sizes as seen in cancer clinical trial evaluation. The frequentist estimation faced some limitations when fitting complex joint models with a high number of random effects and much increased computation time compared to **INLA**. The Bayesian estimation proposed in the R package **INLA** has
<figure><img src="image_11.png" /><figcaption>FIGURE 8 Mean biomarker value according to treatment received for patients with wild-type KRAS status (left) and mutant KRAS status (right). The 95% credible intervals are obtained by resampling from the posterior parameter distributions.</figcaption></figure>

been recently introduced to fit complex joint models (Van Niekerk, Bakka, Rue, & Schenk, 2021) but, to our knowledge, has never been proposed for TPJMs. Accounting for the semicontinuous nature of the biomarker, that is, the SLD, and being able to fit joint models with more complex association structures between the biomarker and the terminal event, can be quite relevant in clinical applications by providing critical insights into the direct and indirect effect of a treatment on the event of interest. This was illustrated in our simulations and applications to two randomized cancer clinical trials.

In our simulation studies, the estimation with **INLA** was found superior to **frailtypack** in terms of computation time and precision of the fixed effects estimation. The point estimates from **frailtypack** yielded closer results to the true values of the random effects' standard deviations, the residual error term, and the baseline hazard function than the posterior mean from **INLA**, even though **INLA** recovered all parameters well based on the estimated credible intervals.

Our first application to the GERCOR randomized clinical trial investigating two treatment lines to treat metastatic colorectal cancer shows some differences between the two estimation approaches. In line with our simulations, the variability of the association parameters estimates between the biomarker and the survival outcome is reduced with **INLA** compared to **frailtypack**. Moreover, the computation time is reduced by a factor of more than 200 with **INLA** compared to **frailtypack** for this application to GERCOR data. The second application to the PRIME study illustrates the fact that treatment response might depend on genetic alterations or tumor biomarker status (DNA/RNA/protein features). There is now a great interest in identifying subgroups of patients with specific patterns of responses; however, most methods provide an average effect of covariates. Instead, our model can distinguish complete responders (i.e., SLD = 0) from partial responders (i.e., SLD > 0). This leads also to an increase in model complexity as additional covariates and random effects are included in each submodel of the TPJM. The frequentist approach proposed in **frailtypack** might have convergence issues in that situation. Interestingly, the analysis of the continuous part of the TPJM suggested that the subgroup of patients with the KRAS mutation receiving treatment B had a similar decrease of the SLD over time compared to the KRAS mutation group receiving treatment A or patients who received treatment B with wild-type KRAS status. Therefore, the lack of response to the addition of anti-EGFR to FOLFOX4 chemotherapy was not fully explained by the KRAS mutation status. This could motivate further investigations of the interaction between KRAS mutation and anti-EGFR therapies to treat advanced colorectal cancer patients, in particular by including information on other somatic tumor mutations (e.g., BRAF or NRAS mutations).

Our work has several limitations. Our applications focused on clinical trials of very advanced cancers, which often have high death rates, short follow-up, and small proportions of complete responses (i.e., SLD = 0). In situations where we have a higher proportion of complete responders, the relative performances of **INLA** versus **frailtypack** could be different. The conclusions might be different for different settings (i.e., with a higher zero rate and reduced censoring). For instance, a meta-analysis evaluating the responses among non-Hodgkin's lymphoma patients estimated complete response rates (i.e., SLD = 0) ranging from 1.2% to 84% in the different pooled clinical trials (Mangal et al., 2018). We also notice that the two models estimated with **INLA** and **frailtypack** are not completely comparable because of the difference in the
approximation of the baseline hazard function. Besides the shared random effects, other association structures have also been proposed such as the current value association, that is, it uses the current level of the biomarker and is available in **frailtypack**. For the TPJMs, the current value of the biomarker is defined as $E[Y_{ij}] = \text{Prob}(Y_{ij} > 0)E[Y_{ij} | Y_{ij} > 0]$, which is nonlinear. It cannot be directly defined as part of the latent Gaussian model, and more work is warranted to include this development in **INLA**. So at the stage of development, **frailtypack** still provides more flexibility when specifying an association structure between the biomarker and the survival outcome. Besides, dynamic predictions of the event of interest are not yet implemented in **INLA** and thus require postcomputations but are a major component of **frailtypack**, available for a wide range of joint models. It would also be interesting to consider a Bayesian development for the marginal TPJM we recently proposed (Rustand et al., 2021). Finally, the definition of the hyperparameter prior distributions is an important aspect of Bayesian estimation. In this work, the PC priors provided a general setting for the priors since they provide a natural avenue to incorporate knowledge from the practitioner about the expected size of the parameter and they are constructed to be proper and avoid overfitting. Alternative prior choices for the hyperparameters can be used in **INLA** if the practitioner possesses motivation for it from an expert or prior knowledge.

The reduction in the computation times with **INLA** was beyond our expectations. It improves drastically the applicability of the Bayesian estimation for complex models such as the TPJMs and other families of joint models, such as a bivariate joint model for recurrent events and a terminal event or a trivariate joint model for a longitudinal biomarker, recurrent events, and a terminal event, which are currently available in **frailtypack**. Finally, **INLA** can accommodate multiple longitudinal outcomes while **frailtypack** is currently limited to a single longitudinal outcome.

## ACKNOWLEDGMENTS

This publication is based on research using information obtained from www.projectdatasphere.org, which is maintained by Project Data Sphere, LLC. Neither Project Data Sphere, LLC nor the owner(s) of any information from the website has contributed to, approved, or are in any way responsible for the contents of this publication. The authors acknowledge the insightful and constructive comments made by the associate editor and two reviewers. These comments have greatly helped to sharpen the original submission.

## CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

## DATA AVAILABILITY STATEMENT

Data from the GERCOR study are not publicly available due to privacy or ethical restrictions. Data from the PRIME study are freely available on the Project Data Sphere platform, ProjectDataSphere.org (PDS UID: Colorec-Amgen_2006_309).

## OPEN RESEARCH BADGES

<figure><img src="image_12.png" /></figure>

This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available in the **Supporting Information** section.

This article has earned an open data badge “**Reproducible Research**” for making publicly available the code necessary to reproduce the reported results. The results reported in this article could fully be reproduced.

## ORCID

Denis Rustand https://orcid.org/0000-0001-9708-5220

Janet van Niekerk https://orcid.org/0000-0002-4334-2057

## REFERENCES

Andrinopoulou, E.-R., & Rizopoulos, D. (2016). Bayesian shrinkage approach for a joint model of longitudinal and survival outcomes assuming different association structures. *Statistics in Medicine*, *35*(26), 4813-4823.

Bokemeyer, C., Bondarenko, I., Hartmann, J., De Braud, F., Volovat, C., Nippgen, J., Stroh, C., Celik, I., & Koralewski, P. (2008). KRAS status and efficacy of first-line treatment of patients with metastatic colorectal cancer (mCRC) with FOLFOX with or without Cetuximab: The OPUS experience. *Journal of Clinical Oncology*, *26*(suppl 15), 4000-4000.

Chi, Y.-Y., & Ibrahim, J. G. (2006). Joint models for multivariate longitudinal and multivariate survival data. *Biometrics*, *62*(2), 432-445.

Commenges, D., Joly, P., Gegout-Petit, A., & Liquet, B. (2007). Choice between semi-parametric estimators of Markov and non-Markov multi-state models from coarsened observations. *Scandinavian Journal of Statistics*, *34*(1), 33-52.
Douillard, J.-Y., Oliner, K. S., Siena, S., Tabernero, J., Burkes, R., Barugel, M., Humblet, Y., Bodoky, G., Cunningham, D., Jassem, J., Rivera, F., Kocákova, I., Ruff, P., Blasińska-Morawiec, M., Šmakal, M., Canon, J. L., Rother, M., Williams, R., Rong, A., ... Patterson, S. D. (2013). Panitumumab-Folfox4 treatment and RAS mutations in colorectal cancer. *New England Journal of Medicine*, **369**(11), 1023-1034.

Hanson, T. E., Branscum, A. J., & Johnson, W. O. (2011). Predictive comparison of joint longitudinal-survival modeling: A case study illustrating competing approaches. *Lifetime Data Analysis*, **17**(1), 3–28.

Henderson, R., Diggle, P., & Dobson, A. (2000). Joint modelling of longitudinal measurements and event time data. *Biostatistics*, **1**(4), 465–480.

Hespanhol, L., Vallio, C. S., Costa, L. M., & Saragiotto, B. T. (2019). Understanding and interpreting confidence and credible intervals around effect estimates. *Brazilian Journal of Physical Therapy*, **23**(4), 290–301.

Król, A., Mauguen, A., Mazroui, Y., Laurent, A., Michiels, S., & Rondeau, V. (2017). Tutorial in joint modeling and prediction: A statistical software for correlated longitudinal outcomes, recurrent events and a terminal event. *Journal of Statistical Software*, **81**(3), 1–52.

Król, A., Tournigand, C., Michiels, S., & Rondeau, V. (2018). Multivariate joint frailty model for the analysis of nonlinear tumor kinetics and dynamic predictions of death. *Statistics in Medicine*, **37**(13), 2148–2161.

Kurz, C. F. (2017). Tweedie distributions for fitting semicontinuous health care utilization cost data. *BMC Medical Research Methodology*, **17**(1), 1–8.

Lawrence Gould, A., Boye, M. E., Crowther, M. J., Ibrahim, J. G., Quartey, G., Micallef, S., & Bois, F. Y. (2015). Joint modeling of survival and longitudinal non-survival data: Current methods and issues. Report of the DIA Bayesian Joint Modeling Working Group. *Statistics in Medicine*, **34**(14), 2181–2195.

Mangal, N., Salem, A. H., Li, M., Menon, R., & Freise, K. J. (2018). Relationship between response rates and median progression-free survival in non-Hodgkin's lymphoma: A meta-analysis of published clinical trials. *Hematological Oncology*, **36**(1), 37–43.

Marabelle, A., Fakih, M., Lopez, J., Shah, M., Shapira-Frommer, R., Nakagawa, K., Chung, H. C., Kindler, H. L., Lopez-Martin, J. A., Miller, W. H., Italiano, A., Kao, S., Piha-Paul, S. A., Delord, J.-P., McWilliams, R. R., Fabrizio, D. A., Aurora-Garg, D., Xu, L., Jin, F., ... Bang, Y.-J. (2020). Association of tumormutational burden with outcomes in patients with advanced solid tumours treated with Pembrolizumab: Prospective biomarker analysis of the multicohort, open-label, phase 2 keynote-158 study. *The Lancet Oncology*, **21**(10), 1353–1365.

Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. *Journal of the Society for Industrial and Applied Mathematics*, **II**(2), 431–441.

Martino, S., Akerkar, R., & Rue, H. (2011). Approximate Bayesian inference for survival models. *Scandinavian Journal of Statistics*, **38**(3), 514–528.

Martins, T. G., Simpson, D., Lindgren, F., & Rue, H. (2013). Bayesian computing with INLA: New features. *Computational Statistics & Data Analysis*, **67**, 68–83.

Muth, C., Oravecz, Z., & Gabry, J. (2018). User-friendly Bayesian regression modeling: A tutorial with Rstanarm and Shinystan. *Quantitative Methods for Psychology*, **14**(2), 99–119.

Normanno, N., Tejpar, S., Morgillo, F., De Luca, A., Van Cutsem, E., & Ciardiello, F. (2009). Implications for KRAS status and EGFR-targeted therapies in metastatic CRC. *Nature Reviews Clinical Oncology*, **6**(9), 519.

Philipps, V., Hejblum, B. P., Prague, M., Commenges, D., & Proust-Lima, C. (2021). Robust and efficient optimization using a Marquardt-Levenberg algorithm with R package marqLevAlg. *The R Journal*, **13**(2), 365–379.

Brown, E. R., & Ibrahim, J. G. (2003). A Bayesian semiparametric joint hierarchical model for longitudinal and survival data. *Biometrics*, **59**(2), 221–228.

Rizopoulos, D. (2016). The R package jmbayes for fitting joint models for longitudinal and time-to-event data using MCMC. *Journal of Statistical Software*, **72**, i07.

Rizopoulos, D., & Ghosh, P. (2011). A Bayesian semiparametric multivariate joint model for multiple longitudinal outcomes and a time-to-event. *Statistics in Medicine*, **30**(12), 1366–1380.

Rue, H., & Held, L. (2005). *Gaussian Markov random fields: Theory and applications*. CRC Press.

Rue, H., Martino, S., & Chopin, N. (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**(2), 319–392.

Rue, H., Riebler, A., Sørbye, S. H., Illian, J. B., Simpson, D. P., & Lindgren, F. K. (2017). Bayesian computing with INLA: A review. *Annual Review of Statistics and Its Application*, **4**(1), 395–421.

Rustand, D., Briollais, L., & Rondeau, V. (2023). A marginalized two-part joint model for a longitudinal biomarker and a terminal event with application to advanced head and neck cancers. [Manuscript in preparation].

Rustand, D., Briollais, L., Tournigand, C., & Rondeau, V. (2020). Two-part joint model for a longitudinal semicontinuous marker and a terminal event with application to metastatic colorectal cancer data. *Biostatistics*, **23**(1), 50–68.

Schenk, O., & Gärtner, K. (2004). Solving unsymmetric sparse systems of linear equations with PARDISO. *Future Generation Computer Systems*, **20**(3), 475–487.

Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sørbye, S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. *Statistical Science*, **32**(1), 1–28.

Smith, V. A., Preisser, J. S., Neelon, B., & Maciejewski, M. L. (2014). A marginalized two-part model for semicontinuous data. *Statistics in Medicine*, **33**(28), 4891–4903.

Song, X., Davidian, M., & Tsiatis, A. A. (2002). A semiparametric likelihood approach to joint modeling of longitudinal and time-to-event data. *Biometrics*, **58**(4), 742–753.

Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. (2002). Bayesian measures of model complexity and fit. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, **64**(4), 583–639.
Tournigand, C., André, T., Achille, E., Lledo, G., Flesh, M., Mery-Mignard, D., Quinaux, E., Couteau, C., Buyse, M., Ganem, G., Landi, B., Colin, P., Louvet, C., & de Gramont, A. (2004). Folfiri followed by Folfox6 or the reverse sequence in advanced colorectal cancer: A randomized GERCOR study. *Journal of Clinical Oncology*, **22**(2), 229-237.

Van Cutsem, E., Lang, I., D'haens, G., Moiseyenko, V., Zaluski, J., Folprecht, G., Tejpar, S., Kisker, O., Stroh, C., & Rougier, P. (2008). KRAS status and efficacy in the first-line treatment of patients with metastatic colorectal cancer (mCRC) treated with Folfiri with or without Cetuximab: The crystal experience. *Journal of Clinical Oncology*, **26**(suppl 15), 2-2.

Van Niekerk, J., Bakka, H., & Rue, H. (2019). *Joint models as latent Gaussian models-not reinventing the wheel*. arXiv. https://doi.org/10.48550/arXiv.1901.09365

Van Niekerk, J., Bakka, H., & Rue, H. (2021). Stable non-linear generalized Bayesian joint models for survival-longitudinal data. *Sankhya A*, **1-27**.

Van Niekerk, J., Bakka, H., Rue, H., & Schenk, O. (2021). New frontiers in Bayesian modeling using the INLA package in R. *Journal of Statistical Software*, **100**, 1-28.

Walker, A. M. (1969). On the asymptotic behaviour of posterior distributions. *Journal of the Royal Statistical Society: Series B (Methodological)*, **31**(1), 80-88.

Wulfsohn, M. S., & Tsiatis, A. A. (1997). A joint model for survival and longitudinal data measured with error. *Biometrics*, **53**, 330-339.

Xu, J., & Zeger, S. L. (2001). Joint analysis of longitudinal data comprising repeated measures and times to events. *Journal of the Royal Statistical Society: Series C (Applied Statistics)*, **50**(3), 375-387.

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Rustand, D., van Niekerk, J., Rue, H., Tournigand, C., Rondeau, V., & Briollais, L. (2023). Bayesian estimation of two-part joint models for a longitudinal semicontinuous biomarker and a terminal event with INLA: Interests for cancer clinical trial evaluation. *Biometrical Journal*, 2100322.

https://doi.org/10.1002/bimj.202100322