# A Location-Scale Joint Model for Studying the Link Between the Time-Dependent Subject-Specific Variability of Blood Pressure and Competing Events

Léonie Courcoul, Christophe Tzourio, Mark Woodward, Antoine Barbieri,

Hélène Jacqmin-Gadda

► To cite this version:

Léonie Courcoul, Christophe Tzourio, Mark Woodward, Antoine Barbieri, Hélène Jacqmin-Gadda. A Location-Scale Joint Model for Studying the Link Between the Time-Dependent Subject-Specific Variability of Blood Pressure and Competing Events. Statistics in Medicine, 2025, 44 (20-22), pp.e70244. 10.1002/sim.70244. hal-05482173

HAL Id: hal-05482173

https://u-bordeaux.hal.science/hal-05482173v1

Submitted on 30 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

OPEN ACCESS

# A Location-Scale Joint Model for Studying the Link Between the Time-Dependent Subject-Specific Variability of Blood Pressure and Competing Events

Léonie Courcoul<fnref n="1" /> | Christophe Tzourio<fnref n="1" /> | Mark Woodward²,<fnref n="3" /> | Antoine Barbieri<fnref n="1" /> | Hélène Jacqmin-Gadda<fnref n="1" />

**Correspondence:** Léonie Courcoul (leonie.courcoul@u-bordeaux.fr; courcoul.leonie498@gmail.com)

**Received:** 22 May 2024 | **Revised:** 23 July 2025 | **Accepted:** 18 August 2025

**Funding:** This work was supported by the Agence Nationale de la Recherche (Grant No. ANR-21-CE36).

**Keywords:** blood pressure | cardio and cerebrovascular diseases | competing events | heterogeneous variance | joint model | location-scale model

## ABSTRACT

Given the high incidence of cardio and cerebrovascular diseases (CVD), and their association with morbidity and mortality, their prevention is a major public health issue. A high level of blood pressure is a well-known risk factor for these events, and an increasing number of studies suggest that blood pressure variability may also be an independent risk factor. However, these studies suffer from significant methodological weaknesses. In this work, we propose a new location-scale joint model for the repeated measures of a marker and competing events. This joint model combines a mixed model including a subject-specific and time-dependent residual variance modeled through random effects, and cause-specific proportional intensity models for the competing events. The risk of events may depend simultaneously on the current value of the variance, as well as, the current value and the current slope of the marker trajectory. The model is estimated by maximizing the likelihood function using the Marquardt-Levenberg algorithm. The estimation procedure is implemented in an R-package and is validated through a simulation study. This model is applied to study the association between blood pressure variability and the risk of CVD and death from other causes. Using data from a large clinical trial on the secondary prevention of stroke, we find that the current individual variability of blood pressure is associated with the risk of CVD and death. Moreover, the comparison with a model without heterogeneous variance shows the importance of taking into account this variability in the goodness-of-fit and for dynamic predictions.

## 1 | Introduction

Cardiovascular diseases, such as ischaemic heart disease, and cerebrovascular events are two leading causes of death. Moreover, these diseases often lead to acquired physical disability or to dementia. Medical care and disability management following this type of disease generate significant societal, human, and financial distress [1]. Given the frequency of cardio and cerebrovascular diseases (CVD) and their dramatic consequences at

the individual and societal level, the identification of modifiable risk factors is essential to implement prevention programs. Hypertension (high values of blood pressure) is a well-known risk factor for these diseases. More recently, some studies have suggested that the visit-to-visit variability of blood pressure could be associated with an increased risk of stroke and cardiovascular events independently of the level of blood pressure [2–4]. These studies have used the individual empirical standard deviation, or some other measure of variation (e.g., the coefficient of

<footnote n="1">¹Univ. Bordeaux, INSERM, Bordeaux Population Health, U1219, France | ²The George Institute for Global Health, Imperial College London, UK |</footnote>

<footnote n="3">³The George Institute for Global Health, University of New South Wales, Sydney, Australia</footnote>

<footnote>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.</footnote>

<footnote>© 2025 The Author(s). *Statistics in Medicine* published by John Wiley & Sons Ltd.</footnote>
variation) or extreme value (e.g., the maximum), of blood pressure as an explanatory variable in a Cox model for the event risk. However, they were exposed to methodological issues. A first strategy consists of calculating the empirical standard deviation of blood pressure on all available measurements [5]. This strategy induces conditioning on the future, likely leading to bias because measurements after the current time (and sometimes after the event time) are used to predict the event at the current time [6, 7]. A second strategy consists of computing the standard deviation of blood pressure on the measurements collected over an initial period of the study, keeping in the sample only the individuals who did not have the event before the end of this period in order to predict the risk beyond this period [7]. This could induce selection bias and certainly creates loss of power. To avoid these issues, the standard deviation of blood pressure can be considered as a time-dependent variable and calculated using only measurements before the event. Nevertheless, this approach neglects the measurement error of the standard deviation, which is a serious issue when the number of measurements differs between individuals, and requires imputation of the standard deviation at all event times. These limitations may introduce bias [8].

Recently, De Courson et al. [7] compared these different approaches and obtained contradictory results depending on the estimator used for the variability. Furthermore, they pointed out that blood pressure and its standard deviation are endogenous variables for which the Cox model is not suitable [9]. Thus, they also considered a joint model approach combining a mixed model to fit individual trajectories of standard deviation of blood pressure measurements and a proportional hazard model to assess the impact of the current variability on the event risk. In the joint model, the current value of standard deviation of blood pressure is included as a time-dependent explanatory variable in the time-to-event model. This allows the evaluation of the impact of the longitudinal data on the event risk without bias, contrary to the two stage estimation [10–12]. However, the interpretation of the association between the event risk and the variability computed as the standard deviation of measurements observed since the beginning is difficult. Indeed, this global variability encompasses the individual time trend of the blood pressure while the clinical question is the following: is irregularity among blood pressure measures a risk factor independently of current values and possibly time-trend of blood pressure?

Location-scale mixed models have been introduced to investigate the heterogeneity of intra-subject variability for longitudinal data by introducing random effects in the residual variance modeling [13]. For studying the association between the variability of a biomarker and a clinical event, Gao et al. [14] and Barrett et al. [15] have proposed a joint model combining a mixed model including a subject-specific random effect for the residual variance and a proportional hazard model for the event risk. However, the considered dependence structure is quite restrictive since, in their models, the event risk depends only on the random effects and not on time-dependent characteristics of the marker trajectory, such as the current value or the current slope. In addition, none of them consider time-dependent subject-specific variability of the maker and they do not handle competing events. However, it is essential to account for competing death from other causes because mortality and cardiovascular risk increase with age and may be both associated with blood pressure.

The objective of our work was, therefore, to propose a new location-scale joint model accounting for both time-dependent individual variability of a marker and competing events. To do this, we extended the model proposed by Gao et al. [14] and Barrett et al. [15] to include a time-dependent variability, competing events, a more flexible dependence structure between the event and the marker trajectory, and more flexible baseline risk functions. In contrast to the previous works we propose a frequentist estimation approach which is implemented in the R-package LSJM.

This paper is organized as follows. Section 2 describes the model and the estimation procedure using a robust algorithm for maximizing the likelihood. Section 3 presents a simulation study to assess the estimation procedure performance. In Section 4, the model is applied to the data from the Perindopril Protection Against Stroke Study (PROGRESS) clinical trial, a blood-pressure lowering trial for the secondary prevention of stroke [16]. Finally, Section 5 concludes this work with some elements of discussion.

## 2 | Method

Let us consider a sample of $N$ individuals. For each individual $i \in \{1, \dots, N\}$, we consider the $n_i$-vector of repeated measures $Y_i = (Y_{i1}, \dots, Y_{in_i})^\top$ with $Y_{ij}$ the value of the longitudinal outcome of individual $i$ at time $t_{ij}$ ($j = 1, \dots, n_i$). Assuming $K$ competing events, we denote $T_i = \min(C_i, T_{ik}^*, k = 1, \dots, K)$ the observed time with $T_{ik}^*$ the real time for the event $k$ ($k = 1, \dots, K$) and $C_i$ the censoring time for the $i$th individual. We then denote $\delta_i \in \{0, 1, \dots, K\}$ the individual event indicator such as $\delta_i = k$ if the competing event $k \in \{1, \dots, K\}$ occurs and $\delta_i = 0$ otherwise.

### 2.1 | Joint Model With Time-Dependent Individual Variability

We propose joint modeling for a longitudinal outcome and competing events using a shared random-effect approach. Joint models allow simultaneous analysis of longitudinal data and clinical events. They combine a mixed model for repeated measures of exposure and a time-to-event models. Functions of the random effects from the mixed model are included as explanatory variables in the time-to-event models to account for the association between the two outcomes. The longitudinal submodel is defined by a linear mixed-effect model with heterogeneous variance:

$$ \begin{cases} Y_{ij} = Y_i(t_{ij}) = \tilde{y}_i(t_{ij}) + \epsilon_{ij} = X_i(t_{ij})^\top \beta + Z_i(t_{ij})^\top b_i + \epsilon_i(t_{ij}), \\ \epsilon_{ij}(t_{ij}) \sim \mathcal{N}(0, \sigma_i^2(t_{ij})) \quad \text{with} \quad \log(\sigma_i(t_{ij})) = O_i(t_{ij})^\top \mu + M_i(t_{ij})^\top \tau_i \end{cases} \quad (1) $$

$\tilde{y}_i(t)$ is the value at time $t$ of a continuous time process that represents the true underlying value of the variable $Y$ at any time (without measurement error); $X_i(t_{ij})$ and $O_i(t_{ij})$ are two vectors of explanatory variables for subject $i$ at visit $j$, possibly including fixed variables, functions of time, interactions between fixed variables and functions of time and time-dependent variables assumed to be known at any time; $Z_i(t_{ij})$ and $M_i(t_{ij})$ are, respectively sub-vectors of $X_i(t_{ij})$ and $O_i(t_{ij})$ including a column of 1 and functions of time; $\beta$ and $\mu$ are vectors of fixed effects while $b_i$ and $\tau_i$ are the vectors of subject-specific random effects such as
$$ \begin{pmatrix} b_i \\ \tau_i \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \Sigma_b & \Sigma_{\epsilon b} \\ \Sigma_{\tau b} & \Sigma_{\tau} \end{pmatrix} \right) $$

The risk function for the event $k \in \{1, \dots, K\}$ is defined by

$$ \lambda_{ik}(t) = \lambda_{0k}(t) \exp(W_{ik}^{\top} \gamma_k + \alpha_{1k} \tilde{y}_i(t) + \alpha_{2k} \tilde{y}'_i(t) + \alpha_{\sigma k} \sigma_i(t)) \quad (2) $$

with $\lambda_{0k}(t)$ the baseline risk function, $W_{ik}$ a vector of baseline covariates associated with the regression coefficient $\gamma_k$. The hazard rates at time $t$ depend on characteristics of the marker trajectory at time $t$: the underlying value of the marker at the current time $\tilde{y}_i(t)$, the slope of the marker trajectory at the current time $\tilde{y}'_i(t) = d\tilde{y}_i(t)/dt$ and/or the individual error variance of the marker at the current time $\sigma_i(t)$. These associations are, respectively, measured by parameters $\alpha_{1k}$, $\alpha_{2k}$ and $\alpha_{\sigma k}$. This markovian assumption could be relaxed by introducing a time lag. Different parametric forms for the baseline risk function can be considered, such as exponential, Weibull, or, for more flexibility, a B-splines base with $Q$ knots defined by:

$$ \lambda_{0k}(t) = \exp\left(\sum_{q=1}^{Q+4} n_{qk} B_q(t, v_k)\right) $$

where $B_q(t, v_k)$ is the $q$th basis function of B-splines with the knot vector $v_k$ and $n_{qk}$ is the associated parameter to be estimated.

## 2.2 | Estimation Procedure

Let $\theta$ be the set of parameters to be estimated including parameters of the Cholesky decomposition of the covariance matrix of the random effects, $\beta$, $\mu$, $\alpha^\top = (\alpha_{11}, \alpha_{21}, \alpha_{\sigma 1}, \dots, \alpha_{1K}, \alpha_{2K}, \alpha_{\sigma K})$, $\gamma^\top = (\gamma_1^\top, \dots, \gamma_K^\top)$ and the parameters of the $K$ baseline risk functions. Considering the frequentist approach, the parameters are estimated by maximizing the likelihood function. Assuming $T_i$ independent of $C_i$ given $W_i$ and marker trajectory characteristics, the contribution of individual $i$ to the marginal likelihood is defined by:

$$ \begin{aligned} \mathcal{L}_i(\theta; Y_i, T_i, \delta_i) &= \int p(Y_i, T_i, \delta_i | b_i, \tau_i; \theta) f(b_i, \tau_i; \theta) db_i d\tau_i \\ &= \int f(Y_i | b_i, \tau_i; \theta) \exp\left(-\sum_{k=1}^K \Lambda_{ik}(T_i | b_i, \tau_i; \theta)\right) \\ & \qquad \prod_{k=1}^K \lambda_{ik}(T_i | b_i, \tau_i; \theta)^{\frac{1}{\delta_i}} f(b_i, \tau_i; \theta) db_i d\tau_i \end{aligned} $$

with $f(b_i, \tau_i; \theta)$ a multivariate Gaussian density and $f(Y_i | b_i, \tau_i; \theta) = \prod_{j=1}^{n_i} f(Y_{ij} | b_j, \tau_i; \theta)$ where $f(Y_{ij} | b_i, \tau_i; \theta)$ is a univariate Gaussian density; $\Lambda_{ik}(T_i | b_i, \tau_i; \theta)$ is the cumulative risk function for event $k$ ($k \in \{1, \dots, K\}$) given by:

$$ \Lambda_{ik}(t | b_i, \tau_i; \theta) = \int_0^t \lambda_{ik}(u | b_i, \tau_i; \theta) du \quad (3) $$

In cohort studies, data are frequently left-truncated. Left-truncation arises as soon as the time scale is not the time since inclusion and the subjects are enrolled only if they are free of the event at inclusion [17]. This is the case in most studies where the time scale is age. To deal with left truncation (also called

delayed entry), the individual contribution to the likelihood must be divided by the probability to be free of any event at entry time $T_{0i}$:

$$ \mathcal{L}_i^{\text{DE}}(\theta; Y_i, T_i, \delta_i) = \frac{\mathcal{L}_i(\theta; Y_i, T_i, \delta_i)}{\int \exp(-\sum_{k=1}^K \Lambda_{ik}(T_{0i} | b_i, \tau_i; \theta)) f(b_i, \tau_i; \theta) db_i d\tau_i} $$

Because the integral on the random effects does not have an analytical solution, the integral is computed by a Quasi Monte Carlo (QMC) approximation [18], using deterministic quasi-random sequences. The approximation of the integral is defined by

$$ \mathcal{L}_i(\theta; Y_i, T_i, \delta_i) \approx \frac{1}{S} \sum_{s=1}^{S} p(Y_i, T_i, \delta_i | b_i^s, \tau_i^s; \theta) $$

where $(b_i^1, \dots, b_i^S)$ and $(\tau_i^1, \dots, \tau_i^S)$ are draws of a S-sample in the sobol sequel for the distribution $f(b_i, \tau_i; \theta)$. Moreover, to approximate the cumulative risk function given in equation (3), we use the Gauss-Kronrod quadrature approximation with 15 points [19].

Parameter estimation is obtained by maximizing the log-likelihood function $\ell(\theta; Y_i, T_i, \delta_i) = \log(\prod_{i=1}^N \mathcal{L}_i(\theta; Y_i, T_i, \delta_i))$. The maximization is performed using the marqLevAlg R-package based on the Marquardt-Levenberg algorithm [20]. The latter is a robust variant of the Newton-Raphson algorithm [21, 22], which iteratively updates the parameters $\theta$ to be estimated until convergence with the following formula at iteration $l+1$:

$$ \theta^{(l+1)} = \theta^{(l)} - \psi_l \tilde{H}(\theta^{(l)})^{-1} \nabla(\ell(\theta^{(l)})) $$

where $\theta^{(l)}$ is the set of parameters at iteration $l$, $\nabla(\ell(\theta^{(l)}))$ the gradient of the log-likelihood at iteration $l$ and $\tilde{H}(\theta^{(l)})$ the inflated Hessian matrix where the diagonal terms of the Hessian matrix $H(\theta^{(l)})$ are replaced by

$$ \tilde{H}(\theta^{(l)})_{ii} = H(\theta^{(l)})_{ii} + \phi_l \left[ (1 - \rho_l) \left| H(\theta^{(l)})_{ii} \right| + \rho_l \mathrm{tr}(\theta^{(l)}) \right] $$

The scalars $\psi_l$, $\phi_l$, and $\rho_l$ are internally determined at each iteration $l$ to ensure that $\tilde{H}(\theta^{(l)})$ be definite-positive. $\tilde{H}(\theta^{(l)})$ approaches $H(\theta^{(l)})$ when $\theta$ approaches $\hat{\theta}$ and insure improvement of the likelihood at each iteration. Stringent convergence criteria are used, relying on parameter and function stability, and the relative distance to the maximum computed from the first and second derivatives of the log-likelihood which must not exceed a threshold $\epsilon_d$: $\frac{\nabla(\ell(\theta^{(l)}))(H(\theta^{(l)}))^{-1}\nabla(\ell(\theta^{(l)})}{\epsilon_d}$, with $m$ the number of parameters. This algorithm was previously compared to other algorithms (EM, BFGS and L-BFGS-B) and the results showed that this algorithm was the most reliable [20]. The variances of the estimates are estimated by the inverse of the Hessian matrix computed by finite differences. The variances of the estimated parameters from the covariance matrix of the random effects are computed using the Delta-Method [23].

In the estimation procedure, the choice of the number of QMC draws is an important factor in practice. Indeed, increasing the number of QMC draws ensures good precision of parameter variances for statistical inference but at the cost of a considerable increase in computation time. To limit computation time and
ensure precise estimates of parameters and their standard error, we propose a two-step strategy. We first applied the Marquardt algorithm with a small number S1 of QMC draws (e.g., S1 = 500) until convergence is achieved. The first step provides a good parameter estimation in a reasonable time, but does not guarantee the accuracy of their estimated standard errors. In the second step, to improve the accuracy of the computation of the standard error of the estimates, a few additional iterations are performed with a higher number S2 > S1 of QMC draws (e.g., S2 = 5000) until the Hessian matrix is invertible. Note that for model selection based on likelihood or information criteria, using the results obtained in the first step is sufficient, given the good estimation of model parameters. However, for statistical inference on the final model, we recommend increasing the number of QMCs (step 2).

## 2.3 | Individual Predictions

We implemented the computation of individual probability of having event k between time s and s + t given that the subject i has not experienced any event before time s, all his/her marker measures collected until time s ($Y_i(s) = \{Y_i(t_{ij})|t_{ij} < s\}$), and the set of estimated parameters. The prediction is defined for subject i by

$$ 
\begin{aligned}
\pi_i(s, t; \hat{\theta}) &= P(s < T_i < s + t, \delta_i = k | T_i > s, \mathcal{Y}_i(s), \hat{\theta}) \\
&= \frac{\int \left[ \int_s^{s+t} g(u) \lambda_{ik}(u|b_i, \tau_i, \hat{\theta}) du \right] f(\mathcal{Y}_i(s)|b_i, \tau_i, \hat{\theta}) f(b_i, \tau_i|\hat{\theta}) db_i d\tau_i}{\int g(s)(s|b_i, \tau_i, \hat{\theta}) f(\mathcal{Y}_i(s)|b_i, \tau_i, \hat{\theta}) f(b_i, \tau_i|\hat{\theta}) db_i d\tau_i} \\
&\text{with } g(u) = \exp \left( -\sum_{c=1}^{K} \Lambda_{lc}(u|b_i, \tau_i, \hat{\theta}) \right)
\end{aligned}
\quad (4) $$

As previously, the integral over the random effect is computed by QMC approximation and the integral over time with the Gauss-Kronrod quadrature.

The corresponding 95% confidence interval of predictions is obtained by the following Monte Carlo algorithm.

For L large enough and l = 1, ..., L (L = 1000 for instance):

* Generate $\tilde{\theta}^{(l)} \sim N(\hat{\theta}, V(\hat{\theta}))$ where $V(\hat{\theta})$ is given by the inverse of the Hessian matrix at $\hat{\theta}$

* Compute $\tilde{\pi}_i^{(l)}(s, t; \tilde{\theta}^{(l)})$ from equation (4)

* Compute the 95% confidence interval from the 2.5th and 97.5th percentiles of the L-sample of $\tilde{\pi}_i^{(l)}(s, t; \tilde{\theta}^{(l)})$

## 2.4 | Software

The R-package LSJM has been developed for the estimation of the model, the prediction of the subject-specific random effects, and the computation of the individual predicted probabilities of events. The package allows estimation of a model with an unconstrained time trend for the marker trajectory, constant or time-dependent subject-specific residual variance, one or two events with exponential, Weibull or B-splines baseline risk functions, and a flexible dependence structure between the events and the marker (possibly including the current value, the

current slope and the current subject-specific residual variance). The development version of LSJM is available on Github at the following link: https://github.com/LeonieCourcoul/LSJM. The package also allows the estimation of the linear mixed effect model with heterogeneous variance (or location-scale mixed model) without jointly modeling the time-to-events. This preliminary analysis is required to obtain initial values for the joint model and check that some random-effect variances are not null. Indeed a joint model assuming, for instance, a dependence on the random slope or on the random variance would be unidentifiable if the variances of these random effects were null.

# 3 | Simulations

In order to evaluate the performance of the estimation procedure, we performed a simulation study using a design similar to the application data. In the next sections, we describe the simulation design according to ADEMP framework [24].

## 3.1 | Aim

The aim of our simulation was to assess the good convergence properties of the proposed estimation procedure for the location-scale joint model. By denoting $\theta$, the whole set of parameters of the location-scale joint model, we wanted to obtain: (1) consistent and unbiased estimated parameters $\hat{\theta}$, (2) a consistent estimate of the empirical standard error of $\hat{\theta}$, $\hat{SE}(\hat{\theta})$, and (3) a nominal coverage rate for the 95% confidence interval of $\hat{\theta}$.

## 3.2 | Data-Generating Model

The design of our simulation was based on the application on PROGRESS data, using the following model to generate the data:

$$ 
\left\{ 
\begin{array}{l}
Y(t_{ij}) = 142 + b_{0i} + (3 + b_{1i}) \times t_{ij} + \epsilon_i(t_{ij}) \\
\epsilon_i(t_{ij}) \sim N(0, \sigma_i^2(t_{ij})) \quad \text{with} \quad \log(\sigma_i(t_{ij})) = 2.4 + \tau_{0i} + (0.05 + \tau_{1i}) \times t_{ij}
\end{array}
\right.
\quad (5) \\
\text{with } u_i = (b_{0i}, b_{1i}, \tau_{0i}, \tau_{1i})^\top \sim N(0, \Sigma).
$$

Visit times were generated using a uniform distribution around fixed time points (± 1 month) between 0 and 5 years.

Competing event times $T_{ik}^*$ ($k = 1, 2$) were generated using the Brent's univariate root-finding method [25] according to the following proportional hazards models:

$$ 
\left\{ 
\begin{array}{l}
\lambda_{i1}(t) = 1.1^2 t^{1.1^2-1} \exp(-7 + 0.02 \tilde{y}_i(t) + 0.01 \tilde{y}_i'(t) + 0.07 \sigma_i(t)) \\
\lambda_{i2}(t) = 1.3^2 t^{1.3^2-1} \exp(-4 - 0.01 \tilde{y}_i(t) - 0.14 \tilde{y}_i'(t) + 0.15 \sigma_i(t))
\end{array}
\right.
\quad (6) $$

Individuals were censored at $C_i$ the last visit observed for subject i in the dataset. Finally, the observed time was defined by $T_i = \min(T_{i1}^*, T_{i2}^*, C_i)$. Measures of the marker Y posterior to $T_i$ were removed from the datasets.

We now describe the four main simulation scenarios considered. Scenario A was performed with a maximum of seven times of
measurements, at 0 year, 0.5 year, and then one per year until 5 years. The random effects of the mean $b_i$ were supposed to be independent from the random effects of the variance, $\tau_i$ using the following matrix of covariance:

$$ 
\Sigma = 
\begin{bmatrix}
\mathrm{Var}(b_{0i}) & \mathrm{Cov}(b_{0i}, b_{1i}) & 0 & 0 \\
\mathrm{Cov}(b_{0i}, b_{1i}) & \mathrm{Var}(b_{1i}) & 0 & 0 \\
0 & 0 & \mathrm{Var}(\tau_{0i}) & \mathrm{Cov}(\tau_{0i}, \tau_{1i}) \\
0 & 0 & \mathrm{Cov}(\tau_{0i}, \tau_{1i}) & \mathrm{Var}(\tau_{1i})
\end{bmatrix}
= 
\begin{bmatrix}
207.36 & -17.4 & 0 & 0 \\
-17.4 & 9.28 & 0 & 0 \\
0 & 0 & 0.0001 & -0.0006 \\
0 & 0 & -0.0006 & 0.0157
\end{bmatrix}
$$

Scenario B was generated with a maximum of 13 times of measurements, at 0 year, every 3 months the first year and then twice per year until 5 years (as in the PROGRESS trial) and the random effects of the mean $b_i$ were also supposed to be independent from the random effects of the variance, $\tau_i$. For Scenarios C and D, we have supposed correlated random effects using the following matrix of correlation:

$$ 
\Sigma = 
\begin{bmatrix}
\mathrm{Var}(b_{0i}) & \mathrm{Cov}(b_{0i}, b_{1i}) & \mathrm{Cov}(b_{0i}, \tau_{0i}) & \mathrm{Cov}(b_{0i}, \tau_{1i}) \\
\mathrm{Cov}(b_{0i}, b_{1i}) & \mathrm{Var}(b_{1i}) & \mathrm{Cov}(b_{1i}, \tau_{0i}) & \mathrm{Cov}(b_{1i}, \tau_{1i}) \\
\mathrm{Cov}(b_{0i}, \tau_{0i}) & \mathrm{Cov}(b_{1i}, \tau_{0i}) & \mathrm{Var}(\tau_{0i}) & \mathrm{Cov}(\tau_{0i}, \tau_{1i}) \\
\mathrm{Cov}(b_{0i}, \tau_{1i}) & \mathrm{Cov}(b_{1i}, \tau_{1i}) & \mathrm{Cov}(\tau_{0i}, \tau_{1i}) & \mathrm{Var}(\tau_{1i})
\end{bmatrix}
= 
\begin{bmatrix}
210.25 & -15.95 & 2.9 & -0.145 \\
-15.95 & 9.05 & -0.304 & 0.067 \\
2.9 & -0.304 & 0.1309 & -0.0206 \\
-0.145 & 0.067 & -0.0206 & 0.0141
\end{bmatrix}
$$

Scenario C was performed with the same visit times as in Scenario A and Scenario D with the same visit times as in Scenario B. Each scenario was considered with $N = 500$ and $N = 1000$.

The simulations were repeated 300 times for each scenario, using $S1 = 500$ and $S2 = 5000$ draws for the QMC integration approximation.

## 3.3 | Estimated Models

To estimate the parameters, we fitted the following model for all scenarios:

$$ 
\left\{
\begin{array}{l}
Y(t_{ij}) = \beta_0 + b_{0i} + (\beta_1 + b_{1i}) \times t_{ij} + \epsilon_i(t_{ij}) \\
\epsilon_i(t_{ij}) \sim \mathcal{N}(0, \sigma_i^2(t_{ij})) \quad \text{with} \quad \log(\sigma_i(t_{ij})) = \mu_0 + \tau_{0i} + (\mu_1 + \tau_{1i}) \times t_{ij} \\
\lambda_{ik}(t) = \lambda_{0k}(t) \exp(\alpha_{1k} \bar{y}_i(t) + \alpha_{2k} \bar{y}_i'(t) + \alpha_{\sigma k} \sigma_i(t))
\end{array}
\right.
\quad (7)
$$

with $\lambda_{0k}(t) = \kappa_k t^{\kappa_k - 1} e^{\varepsilon_{0k}}$. We assumed no constraint for the covariance matrix $\Sigma$ for Scenario C and Scenario D, whereas we assumed independence between $b$ and $\tau$ for Scenarios A and B.

## 3.4 | Performance Measures

For each of the 300 replications, we estimated model parameters, as well as their standard errors and 95% confidence intervals. Using 300 replications, we computed four performance measures: the bias, the asymptotic standard error (ASE), the empirical standard error (ESE), and the empirical coverage (CR) of 95% confidence intervals. Empirical coverage was computed as the proportion of replications of our simulation study in which the estimated 95% confidence interval contained the true parameter value. With 300 replications, the Monte Carlo error for the coverage rate in our simulations was $\sqrt{0.95 \times 0.05/300} = 1.26\%$, allowing the estimated CR to vary within [92.5, 97.5].

## 3.5 | Results

Tables 1–3 report the mean estimates, the bias, the empirical, and mean asymptotic standard errors of the estimated parameters and the coverage rates of their 95% confidence intervals for scenarios A, B, and C on 500 individuals. Results for scenario D and for larger sample size are in the Supporting Information (Tables S1 to S5). The estimation procedure provided satisfactory results for the four scenarios of simulations. Indeed, the bias was minimal, the mean asymptotic standard errors and the empirical standard errors were close, and the coverage rates of the 95% confidence interval were close to the nominal value and most often within the interval [92.5, 97.5]. We note that the bias is minimal from the first step, using $S1 = 500$ QMC draws, but the second step, using $S2 = 5000$ QMC draws, helps to reduce the difference between the mean asymptotic and the empirical standard deviations and thus improve the coverage rates. We observed a slight undercoverage for some parameters of the random-effects covariance matrix. This small bias mainly affects covariance parameters that are expected to have negligible impact on individual risk prediction. To assess the impact of small biases in the variance parameter estimates on the individual predictions, we computed the individual predictions of the disease risk between years 3 and 5 among subjects event free at 3 years and given measurements collected up to 3 years (using Equation (4)). On a set of 100 datasets generated according to scenario B, Figure S1.A shows that the individual predictions calculated from the estimated parameters are globally equivalent to those obtained using the true parameter values used to generate the data. Moreover, the individual predictions using all estimated parameters are very close to the predictions obtained by replacing the estimates of variance parameters by their true values (see Figure S1.B). Figure S1 confirms the overall quality of the estimated individual prediction and that the small bias on variance parameters has little impact on the individual predictions. Furthermore, on a simulated dataset, we computed also the 95% confidence interval of the predicted probability, using the approach described in Section 2.3 with $L = 1000$. The true prediction probability systematically belongs to the confidence interval (the coverage rate was 100%).

This simulation study also illustrates the impact of the choice of the number of QMC draws $S1$ and $S2$ on the computation time: for scenario A, the medians of computation time are around 13 min (25 iterations in median) and 9 min (1 iteration in median), respectively, for steps 1 and 2, on 10 cores.
<table><caption>TABLE 1 | Simulation results for scenario A with 500 subjects (7 measures, $b_i$ and $\tau_i$ independent).<sup>a</sup></caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2"></th><th rowspan="2">True value</th><th colspan="5">Step 1 (S1 = 500 QMC)</th><th colspan="5">Step 2 (S2 = 5000 QMC)</th></tr><tr><th>Mean</th><th>Bias</th><th>ESE</th><th>ASE</th><th>CR (%)</th><th>Mean</th><th>Bias</th><th>ESE</th><th>ASE</th><th>CR (%)</th></tr></thead><tbody><tr><td colspan="13" style="text-align: left; font-weight: bold;">Longitudinal submodel</td></tr><tr><td>Intercept</td><td>β<sub>0</sub></td><td>142</td><td>142.1</td><td>0.1</td><td>0.736</td><td>0.717</td><td>94.3</td><td>142.1</td><td>0.1</td><td>0.730</td><td>0.728</td><td>95.0</td></tr><tr><td>Slope</td><td>β<sub>1</sub></td><td>3</td><td>2.943</td><td>-0.057</td><td>0.288</td><td>0.253</td><td>90.0</td><td>2.945</td><td>-0.055</td><td>0.280</td><td>0.271</td><td>92.6</td></tr><tr><td rowspan="2">Variability</td><td>μ<sub>0</sub></td><td>2.4</td><td>2.396</td><td>-0.004</td><td>0.027</td><td>0.026</td><td>94.7</td><td>2.395</td><td>-0.005</td><td>0.027</td><td>0.027</td><td>95.3</td></tr><tr><td>μ<sub>1</sub></td><td>0.05</td><td>0.050</td><td>0</td><td>0.017</td><td>0.015</td><td>93.0</td><td>0.051</td><td>0.001</td><td>0.016</td><td>0.016</td><td>94.3</td></tr><tr><td rowspan="3">Σ<sub>b</sub></td><td>σ<sup>2</sup><sub>b0</sub></td><td>207.36</td><td>205.0</td><td>-2.36</td><td>17.31</td><td>16.17</td><td>91.97</td><td>205.1</td><td>-2.26</td><td>17.19</td><td>16.60</td><td>91.97</td></tr><tr><td>σ<sub>b0b1</sub></td><td>-17.28</td><td>-16.06</td><td>1.22</td><td>4.132</td><td>3.663</td><td>87.96</td><td>-16.04</td><td>1.24</td><td>3.936</td><td>4.020</td><td>92.64</td></tr><tr><td>σ<sup>2</sup><sub>b1</sub></td><td>9.28</td><td>9.351</td><td>0.071</td><td>1.673</td><td>1.308</td><td>84.95</td><td>9.279</td><td>-0.001</td><td>1.536</td><td>1.465</td><td>91.64</td></tr><tr><td rowspan="3">Σ<sub>τ</sub></td><td>σ<sup>2</sup><sub>τ0</sub></td><td>0.0001</td><td>0.004</td><td>0.0039</td><td>0.006</td><td>0.006</td><td>98.66</td><td>0.004</td><td>0.0039</td><td>0.007</td><td>0.006</td><td>98.33</td></tr><tr><td>σ<sub>τ0τ1</sub></td><td>-0.0006</td><td>-0.002</td><td>0.0014</td><td>0.005</td><td>0.005</td><td>92.64</td><td>-0.003</td><td>0.0024</td><td>0.005</td><td>0.006</td><td>95.99</td></tr><tr><td>σ<sup>2</sup><sub>τ1</sub></td><td>0.0157</td><td>0.016</td><td>0.0003</td><td>0.005</td><td>0.004</td><td>91.97</td><td>0.017</td><td>0.0013</td><td>0.005</td><td>0.005</td><td>93.31</td></tr><tr><td colspan="13" style="text-align: left; font-weight: bold;">Survival submodel 1</td></tr><tr><td>Current variance</td><td>α<sub>σ1</sub></td><td>0.07</td><td>0.064</td><td>-0.006</td><td>0.041</td><td>0.039</td><td>94.7</td><td>0.064</td><td>-0.006</td><td>0.041</td><td>0.039</td><td>95.3</td></tr><tr><td>Current value</td><td>α<sub>11</sub></td><td>0.02</td><td>0.020</td><td>0</td><td>0.008</td><td>0.007</td><td>95.7</td><td>0.020</td><td>0</td><td>0.008</td><td>0.007</td><td>95.3</td></tr><tr><td>Current slope</td><td>α<sub>21</sub></td><td>0.01</td><td>0.008</td><td>-0.002</td><td>0.072</td><td>0.066</td><td>94.3</td><td>0.007</td><td>-0.003</td><td>0.071</td><td>0.067</td><td>96.0</td></tr><tr><td>Weibull</td><td>√<span style="text-decoration: overline;">κ<sub>1</sub></span></td><td>1.1</td><td>1.099</td><td>-0.001</td><td>0.056</td><td>0.059</td><td>97.7</td><td>1.098</td><td>-0.002</td><td>0.056</td><td>0.059</td><td>98.0</td></tr><tr><td></td><td>ξ<sub>01</sub></td><td>-7</td><td>-6.884</td><td>0.116</td><td>1.257</td><td>1.199</td><td>94.3</td><td>-6.885</td><td>0.115</td><td>1.257</td><td>1.204</td><td>94.7</td></tr><tr><td colspan="13" style="text-align: left; font-weight: bold;">Survival submodel 2</td></tr><tr><td>Current variance</td><td>α<sub>σ2</sub></td><td>0.15</td><td>0.169</td><td>0.019</td><td>0.091</td><td>0.046</td><td>92.3</td><td>0.168</td><td>0.018</td><td>0.091</td><td>0.054</td><td>96.3</td></tr><tr><td>Current value</td><td>α<sub>12</sub></td><td>-0.01</td><td>-0.012</td><td>-0.002</td><td>0.014</td><td>0.009</td><td>96.3</td><td>-0.012</td><td>-0.002</td><td>0.015</td><td>0.010</td><td>96.7</td></tr><tr><td>Current slope</td><td>α<sub>22</sub></td><td>-0.14</td><td>-0.171</td><td>0.031</td><td>0.173</td><td>0.086</td><td>92.6</td><td>-0.167</td><td>0.027</td><td>0.170</td><td>0.095</td><td>94.7</td></tr><tr><td>Weibull</td><td>√<span style="text-decoration: overline;">κ<sub>2</sub></span></td><td>1.3</td><td>1.310</td><td>0.010</td><td>0.102</td><td>0.079</td><td>97.0</td><td>1.311</td><td>0.011</td><td>0.105</td><td>0.084</td><td>98.0</td></tr><tr><td></td><td>ξ<sub>02</sub></td><td>-4</td><td>-4.075</td><td>-0.075</td><td>1.389</td><td>1.366</td><td>96.0</td><td>-4.075</td><td>-0.075</td><td>1.388</td><td>1.401</td><td>95.7</td></tr></tbody></table>

Abbreviations: ASE: asymptotic standard error; coverage rate: coverage rate of the 95% confidence interval; ESE: empirical standard error.

## 3.6 | Additional Simulation Scenarios

Two additional misspecified scenarios (E and F) were performed to assess the robustness of the estimates. In these scenarios, only one event was considered. In Scenario E, the marker was generated with a quadratic time trend but estimated with a linear trend. This scenario was performed to assess the effect of a misspecification of the blood pressure time trend on the estimation of the variance effect on the event risk. As expected, the estimates of the fixed effects in the mixed model are biased, but the estimates of the model for the residual variance and of the time-to-event model are robust (Table S6).

In Scenario F, the variance of the marker was generated as constant with time, so without fixed and random slopes, but the estimated model assumed a linear time trend for the variance with fixed and random slopes. This scenario was considered to assess if the estimation algorithm was robust enough to estimate the location-scale joint model with time-varying variance when the variance in the data is constant. The estimates remain unbiased with nominal coverage rates. As expected, the estimates of the fixed slope and the variance of the random slope for the residual

variance are close to zero (Table S7). Finally, this scenario was also estimated considering a constant variance over time, matching the underlying data generation model. The results showed that the estimates were unbiased, and the coverage rates were consistent with nominal value (Table S8).

# 4 | Application

## 4.1 | Progress Clinical Trial

We applied the proposed methodology on the data from the PROGRESS clinical trial [16] designed to evaluate a blood-pressure lowering treatment in secondary prevention. PROGRESS is a multicenter, double-blind randomized clinical trial including patients with a history of stroke or transient ischemic attack within 5 years before inclusion. Patients were recruited between May 1995 and November 1997. The follow-up comprised five visits in the first year, then two visits each year until the end of the study or the occurrence of a major CVD (stroke, myocardial infarction, and cerebral hemorrhage) or death. At each visit, blood pressure was measured twice and,

<footnote><sup>a</sup>Results for 299 replicates with complete convergence over 300.</footnote>
<table><caption>TABLE 2 | Simulation results for scenario B with 500 subjects (13 measures, $b_i$ and $\tau_i$ independent).<sup>a</sup></caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2"></th><th rowspan="2">True value</th><th colspan="5">Step 1 (S1 = 500 QMC)</th><th colspan="5">Step 2 (S2 = 5000 QMC)</th></tr><tr><th>Mean</th><th>Bias</th><th>ESE</th><th>ASE</th><th>CR (%)</th><th>Mean</th><th>Bias</th><th>ESE</th><th>ASE</th><th>CR (%)</th></tr></thead><tbody><tr><td colspan="13" style="text-align: left; font-weight: bold;">Longitudinal submodel</td></tr><tr><td>Intercept</td><td>β<sub>0</sub></td><td>142</td><td>142.0</td><td>0</td><td>0.779</td><td>0.655</td><td>90.7</td><td>142.0</td><td>0</td><td>0.767</td><td>0.721</td><td>92.7</td></tr><tr><td>Slope</td><td>β<sub>1</sub></td><td>3</td><td>2.996</td><td>-0.004</td><td>0.252</td><td>0.201</td><td>90.0</td><td>3.000</td><td>0</td><td>0.248</td><td>0.235</td><td>93.3</td></tr><tr><td rowspan="2">Variability</td><td>μ<sub>0</sub></td><td>2.4</td><td>2.402</td><td>0.002</td><td>0.019</td><td>0.018</td><td>92.0</td><td>2.401</td><td>0.001</td><td>0.019</td><td>0.018</td><td>92.7</td></tr><tr><td>μ<sub>1</sub></td><td>0.05</td><td>0.047</td><td>-0.003</td><td>0.012</td><td>0.011</td><td>92.0</td><td>0.049</td><td>-0.001</td><td>0.012</td><td>0.012</td><td>94.0</td></tr><tr><td rowspan="3">Σ<sub>b</sub></td><td>σ<sup>2</sup><sub>b<sub>0</sub></sub></td><td>207.36</td><td>208.1</td><td>0.74</td><td>17.79</td><td>14.18</td><td>87.00</td><td>208.0</td><td>0.64</td><td>17.29</td><td>15.92</td><td>92.33</td></tr><tr><td>σ<sup>2</sup><sub>b<sub>0</sub>b<sub>1</sub></sub></td><td>-17.28</td><td>-15.74</td><td>1.54</td><td>4.165</td><td>2.954</td><td>77.00</td><td>-15.85</td><td>1.43</td><td>4.022</td><td>3.614</td><td>85.00</td></tr><tr><td>σ<sup>2</sup><sub>b<sub>1</sub></sub></td><td>9.28</td><td>9.236</td><td>-0.044</td><td>1.476</td><td>0.948</td><td>79.00</td><td>9.256</td><td>0.024</td><td>1.322</td><td>1.246</td><td>89.67</td></tr><tr><td rowspan="3">Σ<sub>τ</sub></td><td>σ<sup>2</sup><sub>τ<sub>0</sub></sub></td><td>0.0001</td><td>0.002</td><td>0.0019</td><td>0.003</td><td>0.002</td><td>97.33</td><td>0.002</td><td>0.0019</td><td>0.003</td><td>0.003</td><td>97.67</td></tr><tr><td>σ<sup>2</sup><sub>τ<sub>0</sub>τ<sub>1</sub></sub></td><td>-0.0006</td><td>-0.001</td><td>-0.0004</td><td>0.003</td><td>0.003</td><td>91.33</td><td>-0.001</td><td>-0.0004</td><td>0.003</td><td>0.004</td><td>93.67</td></tr><tr><td>σ<sup>2</sup><sub>τ<sub>1</sub></sub></td><td>0.0157</td><td>0.015</td><td>-0.0007</td><td>0.003</td><td>0.003</td><td>88.67</td><td>0.016</td><td>0.0003</td><td>0.003</td><td>0.003</td><td>94.67</td></tr><tr><td colspan="13" style="text-align: left; font-weight: bold;">Survival submodel 1</td></tr><tr><td>Current variance</td><td>α<sub>σ1</sub></td><td>0.07</td><td>0.063</td><td>-0.007</td><td>0.030</td><td>0.028</td><td>92.7</td><td>0.065</td><td>-0.005</td><td>0.029</td><td>0.028</td><td>93.7</td></tr><tr><td>Current value</td><td>α<sub>11</sub></td><td>0.02</td><td>0.020</td><td>0</td><td>0.006</td><td>0.007</td><td>96.3</td><td>0.020</td><td>0</td><td>0.006</td><td>0.007</td><td>96.3</td></tr><tr><td>Current slope</td><td>α<sub>21</sub></td><td>0.01</td><td>0.007</td><td>-0.003</td><td>0.057</td><td>0.053</td><td>96.3</td><td>0.007</td><td>-0.003</td><td>0.055</td><td>0.055</td><td>96.0</td></tr><tr><td>Weibull</td><td>√<sub>κ<sub>1</sub></sub></td><td>1.1</td><td>1.106</td><td>0.006</td><td>0.055</td><td>0.054</td><td>94.0</td><td>1.105</td><td>0.005</td><td>0.055</td><td>0.055</td><td>95.0</td></tr><tr><td>Weibull</td><td>ζ<sub>01</sub></td><td>-7</td><td>-6.912</td><td>0.088</td><td>0.992</td><td>1.042</td><td>95.7</td><td>-6.914</td><td>0.086</td><td>0.992</td><td>1.050</td><td>95.7</td></tr><tr><td colspan="13" style="text-align: left; font-weight: bold;">Survival submodel 2</td></tr><tr><td>Current variance</td><td>α<sub>σ2</sub></td><td>0.15</td><td>0.158</td><td>0.008</td><td>0.034</td><td>0.030</td><td>92.7</td><td>0.159</td><td>0.009</td><td>0.031</td><td>0.032</td><td>97.3</td></tr><tr><td>Current value</td><td>α<sub>12</sub></td><td>-0.01</td><td>-0.010</td><td>0</td><td>0.008</td><td>0.008</td><td>94.7</td><td>-0.010</td><td>0</td><td>0.008</td><td>0.008</td><td>95.3</td></tr><tr><td>Current slope</td><td>α<sub>22</sub></td><td>-0.14</td><td>-0.146</td><td>-0.006</td><td>0.064</td><td>0.061</td><td>94.0</td><td>-0.145</td><td>-0.005</td><td>0.063</td><td>0.064</td><td>95.0</td></tr><tr><td>Weibull</td><td>√<sub>κ<sub>2</sub></sub></td><td>1.3</td><td>1.299</td><td>-0.001</td><td>0.067</td><td>0.067</td><td>95.7</td><td>1.299</td><td>-0.001</td><td>0.066</td><td>0.068</td><td>96.3</td></tr><tr><td>Weibull</td><td>ζ<sub>02</sub></td><td>-4</td><td>-4.104</td><td>-0.104</td><td>1.173</td><td>1.111</td><td>93.7</td><td>-4.107</td><td>-0.107</td><td>1.173</td><td>1.139</td><td>93.7</td></tr></tbody></table>

Abbreviations: ASE: asymptotic standard error; coverage rate: coverage rate of the 95% confidence interval; ESE: empirical standard error.

as recommended [26], we analyzed the mean of the two measurements at each time. Prior to randomization, eligible patients were subjected to a 4-week run-in phase to test their tolerance to the treatment. At randomization, patients assigned to the control group stopped the treatment. In order to avoid an effect of the change of therapy at randomization, we removed the blood pressure measure at randomization. The model was estimated only on Non Asian subjects due to differences between Asians and non Asians with regard to CVD risk and risk factors and because the treatment was not exactly the same. Finally, the current study was conducted over 3710 Non Asian patients, 1856 for the controlled group and 1854 for the treatment group, and included 672 CVD and 150 deaths without CVD. There are 2525 (68%) men and 1185 (32%) women. The average age at entry in the study is 67 years old (sd = 9.8) with a minimum at 26 and a maximum at 91 years old.

## 4.2 | Specification of the Model

This study aimed to evaluate the impact of the blood pressure variability on the risk of CVD with or without death and death

due to other causes. To do so, we estimated the proposed joint model (Model CVCS+V, for current value, current slope and variance) with heterogeneous time-dependent variance defined by equations (1) and (2) using the time since the first blood pressure measurement post-randomization. The trajectory of blood pressure was described over time by a linear mixed-effect model. The individual time trend of the marker was modeled by a quadratic trend and the residual variance by a linear trend. The baseline hazard functions of both events were defined by B-splines with three interior knots located at the quantiles of the observed event times. According to the AIC, the model with three knots was better than models with 1 or 5 knots for each baseline hazard function (respectively 298651.9, 298659.8, and 298661.2). The model allowed the risk of each event to depend on the time-dependent intra-subject variability, the individual current value and the current slope. The longitudinal submodel and the variance submodel were adjusted for treatment group and survival submodels were adjusted for treatment group, age at baseline, sex (male versus female), and the baseline measurement of blood pressure (BP0) measured before the 4-week run-in phase, to reflect the historical blood pressure level:

<footnote><sup>a</sup>Results for 300 replicates with complete convergence over 300.</footnote>
<table><caption>TABLE 3 | Simulation results for Scenario C with 500 subjects (7 measures, $b_i$ and $\tau_i$ correlated).<sup>a</sup></caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">True value</th><th colspan="5">Step 1 (S1 = 500 QMC)</th><th colspan="5">Step 2 (S2 = 5000 QMC)</th></tr><tr><th>Mean</th><th>Bias</th><th>ESE</th><th>ASE</th><th>CR (%)</th><th>Mean</th><th>Bias</th><th>ESE</th><th>ASE</th><th>CR (%)</th></tr></thead><tbody><tr><td colspan="13"><strong>Longitudinal submodel</strong></td></tr><tr><td>Intercept</td><td>β<sub>0</sub></td><td>142</td><td>141.9</td><td>-0.001</td><td>0.833</td><td>0.742</td><td>92.3</td><td>141.9</td><td>-0.001</td><td>0.820</td><td>0.756</td><td>93.3</td></tr><tr><td>Slope</td><td>β<sub>1</sub></td><td>3</td><td>3.019</td><td>0.019</td><td>0.320</td><td>0.282</td><td>90.6</td><td>3.023</td><td>0.023</td><td>0.314</td><td>0.290</td><td>92.0</td></tr><tr><td rowspan="2">Variability</td><td>μ<sub>0</sub></td><td>2.4</td><td>2.401</td><td>0.001</td><td>0.035</td><td>0.033</td><td>93.6</td><td>2.399</td><td>-0.001</td><td>0.035</td><td>0.033</td><td>94.3</td></tr><tr><td>μ<sub>1</sub></td><td>0.05</td><td>0.050</td><td>0</td><td>0.016</td><td>0.015</td><td>92.6</td><td>0.050</td><td>0</td><td>0.016</td><td>0.016</td><td>93.7</td></tr><tr><td rowspan="8">Σ<sub>bτ</sub></td><td>σ<sup>2</sup><sub>b0</sub></td><td>210.25</td><td>209.7</td><td>-0.55</td><td>20.69</td><td>16.90</td><td>88.23</td><td>209.4</td><td>-0.85</td><td>20.39</td><td>17.72</td><td>91.33</td></tr><tr><td>σ<sub>b0b1</sub></td><td>-15.95</td><td>-15.43</td><td>0.52</td><td>4.734</td><td>3.942</td><td>88.59</td><td>-15.48</td><td>0.47</td><td>4.398</td><td>4.209</td><td>92.33</td></tr><tr><td>σ<sub>b0τ<sub>0</sub></sub></td><td>2.9</td><td>2.796</td><td>-0.104</td><td>0.610</td><td>0.521</td><td>88.93</td><td>2.812</td><td>-0.088</td><td>0.592</td><td>0.533</td><td>89.67</td></tr><tr><td>σ<sub>b0τ<sub>1</sub></sub></td><td>-0.145</td><td>-0.129</td><td>0.016</td><td>0.251</td><td>0.210</td><td>89.54</td><td>-0.115</td><td>0.03</td><td>0.240</td><td>0.220</td><td>93.00</td></tr><tr><td>σ<sup>2</sup><sub>b1</sub></td><td>9.05</td><td>9.181</td><td>0.131</td><td>1.745</td><td>1.371</td><td>85.91</td><td>9.084</td><td>0.034</td><td>1.599</td><td>1.479</td><td>91.33</td></tr><tr><td>σ<sub>b1τ<sub>0</sub></sub></td><td>-0.304</td><td>-0.295</td><td>0.009</td><td>0.178</td><td>0.157</td><td>91.95</td><td>-0.291</td><td>0.013</td><td>0.169</td><td>0.162</td><td>94.67</td></tr><tr><td>σ<sub>b1τ<sub>1</sub></sub></td><td>0.067</td><td>0.069</td><td>0.002</td><td>0.073</td><td>0.061</td><td>88.93</td><td>0.065</td><td>-0.002</td><td>0.068</td><td>0.063</td><td>92.00</td></tr><tr><td>σ<sup>2</sup><sub>τ<sub>0</sub></sub></td><td>0.1309</td><td>0.123</td><td>-0.0079</td><td>0.029</td><td>0.026</td><td>87.53</td><td>0.128</td><td>-0.0029</td><td>0.029</td><td>0.028</td><td>91.67</td></tr><tr><td colspan="13"><strong>Survival submodel 1</strong></td></tr><tr><td>Current variance</td><td>α<sub>σ1</sub></td><td>0.07</td><td>0.065</td><td>-0.005</td><td>0.051</td><td>0.046</td><td>93.3</td><td>0.067</td><td>-0.003</td><td>0.050</td><td>0.047</td><td>94.7</td></tr><tr><td>Current value</td><td>α<sub>11</sub></td><td>0.02</td><td>0.021</td><td>0.001</td><td>0.012</td><td>0.011</td><td>93.3</td><td>0.021</td><td>0.001</td><td>0.012</td><td>0.011</td><td>93.3</td></tr><tr><td>Current slope</td><td>α<sub>21</sub></td><td>0.01</td><td>-0.004</td><td>-0.014</td><td>0.077</td><td>0.076</td><td>93.6</td><td>-0.051</td><td>-0.061</td><td>0.817</td><td>0.470</td><td>94.0</td></tr><tr><td>Weibull</td><td>√<span style="text-decoration: overline;">κ<sub>1</sub></span></td><td>1.1</td><td>1.094</td><td>-0.006</td><td>0.052</td><td>0.054</td><td>96.3</td><td>1.095</td><td>-0.005</td><td>0.052</td><td>0.055</td><td>96.7</td></tr><tr><td>Weibull</td><td>ζ<sub>01</sub></td><td>-7</td><td>-7.133</td><td>-0.133</td><td>1.437</td><td>1.296</td><td>94.0</td><td>-6.958</td><td>0.042</td><td>3.204</td><td>2.612</td><td>94.7</td></tr><tr><td colspan="13"><strong>Survival submodel 2</strong></td></tr><tr><td>Current variance</td><td>α<sub>σ2</sub></td><td>0.15</td><td>0.166</td><td>0.016</td><td>0.079</td><td>0.051</td><td>91.3</td><td>0.165</td><td>0.015</td><td>0.077</td><td>0.057</td><td>95.0</td></tr><tr><td>Current value</td><td>α<sub>12</sub></td><td>-0.01</td><td>-0.016</td><td>-0.006</td><td>0.018</td><td>0.013</td><td>92.0</td><td>-0.012</td><td>-0.002</td><td>0.017</td><td>0.014</td><td>93.7</td></tr><tr><td>Current slope</td><td>α<sub>22</sub></td><td>-0.14</td><td>-0.154</td><td>-0.014</td><td>0.113</td><td>0.088</td><td>95.6</td><td>-0.178</td><td>-0.038</td><td>0.486</td><td>0.253</td><td>95.7</td></tr><tr><td>Weibull</td><td>√<span style="text-decoration: overline;">κ<sub>2</sub></span></td><td>1.3</td><td>1.314</td><td>0.014</td><td>0.078</td><td>0.069</td><td>94.3</td><td>1.314</td><td>0.014</td><td>0.074</td><td>0.072</td><td>94.7</td></tr><tr><td>Weibull</td><td>ζ<sub>02</sub></td><td>-4</td><td>-4.084</td><td>-0.084</td><td>1.758</td><td>1.524</td><td>94.3</td><td>-3.983</td><td>0.017</td><td>2.469</td><td>2.109</td><td>95.0</td></tr></tbody></table>

Abbreviations: ASE: asymptotic standard error; coverage rate: coverage rate of the 95% confidence interval; ESE: empirical standard error.

<sup>a</sup>Results for 298 replicates with complete convergence over 300 for step 1 and for 300 replicates with complete convergence over 300 for step 2.

$$
\left\{
\begin{array}{l}
y(t_{ij}) = \beta_0 + b_{0i} + (\beta_1 + b_{1i})t_{ij} + (\beta_2 + b_{2i})t_{ij}^2 + \beta_3 trt_i \\
\qquad + \beta_4 trt_i t_{ij} + \beta_5 trt_i t_{ij}^2 + \epsilon_i(t_{ij}) \\
\epsilon_i(t_{ij}) \sim \mathcal{N}(0, \sigma_i^2(t_{ij})) \\
\text{with } \log(\sigma_i(t_{ij})) = \mu_0 + \tau_{0i} + (\mu_1 + \tau_{1i})t_{ij} + \mu_2 trt_i \\
\lambda_{ik}(t) = \lambda_{0k}(t) \exp(\gamma_{0k} trt_i + \gamma_{1k} male_i + \gamma_{2k} age_i + \gamma_{3k} BP O_i \\
\qquad + \alpha_{1k} \tilde{y}_i(t) + \alpha_{2k} \tilde{y}'_i(t) + \alpha_{\sigma k} \sigma_i(t))
\end{array}
\right.
$$

slope (CVCS model). It was also compared to a location-scale joint model with a structure similar to that proposed by [15] and [14], featuring time-fixed (i.e., constant) subject-specific residual variance and event risks dependent on the random intercept, the random slope and the subject-specific residual variance.

## 4.3 | Results

The estimation was performed with S1 = 500 and S2 = 10000 draws of QMC to ensure a greater accuracy.

This model was compared with two classical joint model without heterogeneous variance, that is, $\sigma_i^2(t_{ij}) = \sigma^2$ for all $i = 1, \dots, N$ and $j = 1, \dots, n_i$. The first one allowed the risk of each event to depend only on the individual current value (CV model) and the second one on both the individual current value and current

The AIC from the complete model (298651.9) was clearly better than the AIC from the two joint models with a constant residual variance and a dependence either on the current value only (CV model: 302062.4) or on both the current value and the current slope (CVCS model: 301662.9) and was also better than the joint model assuming a time-fixed subject-specific residual variance (298872.2). This highlights the importance of taking into account a time-dependent subject-specific variance.
Table 4 provides estimates from the complete joint model and Table 5 presents the covariance matrix of the random effects and their standard errors computed through the Delta-Method. Blood pressure was lower for individuals from the treatment group ($\beta_3 = -9.354, p < 0.001$). The interactions between randomization and time are not statistically significant; however, when compared to a model without these interaction terms, the full

<table><caption>**TABLE 4** | Parameter estimates of the joint model on the progress clinical trial data (CVCS+V model).</caption><thead><tr><th>Parameter</th><th>Estimate</th><th>Standard error</th><th>p</th></tr></thead><tbody><tr><td colspan="4"><em>Survival submodel for CVD</em></td></tr><tr><td>BP current variance</td><td>0.042</td><td>0.017</td><td>0.012</td></tr><tr><td>BP current value</td><td>-0.0145</td><td>0.005</td><td>0.002</td></tr><tr><td>BP current slope</td><td>-0.093</td><td>0.021</td><td>&lt; 0.001</td></tr><tr><td>BP0</td><td>0.006</td><td>0.003</td><td>0.016</td></tr><tr><td>Treatment group</td><td>-0.222</td><td>0.087</td><td>0.011</td></tr><tr><td>Male</td><td>0.293</td><td>0.089</td><td>0.001</td></tr><tr><td>Age</td><td>0.041</td><td>0.005</td><td>&lt; 0.001</td></tr><tr><td colspan="4"><em>Survival submodel for Death</em></td></tr><tr><td>BP current variance</td><td>0.102</td><td>0.035</td><td>0.004</td></tr><tr><td>BP current value</td><td>-0.021</td><td>0.011</td><td>0.055</td></tr><tr><td>BP current slope</td><td>-0.118</td><td>0.061</td><td>0.051</td></tr><tr><td>BP0</td><td>-0.002</td><td>0.005</td><td>0.672</td></tr><tr><td>Treatment group</td><td>-0.070</td><td>0.182</td><td>0.700</td></tr><tr><td>Male</td><td>0.521</td><td>0.194</td><td>0.007</td></tr><tr><td>Age</td><td>0.053</td><td>0.010</td><td>&lt; 0.001</td></tr><tr><td colspan="4"><em>Longitudinal submodel</em></td></tr><tr><td colspan="4"><em>Blood Pressure Mean</em></td></tr><tr><td>Intercept</td><td>143.3</td><td>0.368</td><td>&lt; 0.001</td></tr><tr><td>Time</td><td>-0.108</td><td>0.266</td><td>0.684</td></tr><tr><td>Time<sup>2</sup></td><td>-0.111</td><td>0.064</td><td>0.084</td></tr><tr><td>Treatment group</td><td>-9.354</td><td>0.518</td><td>&lt; 0.001</td></tr><tr><td>Time × treatment group</td><td>0.678</td><td>0.364</td><td>0.063</td></tr><tr><td>Time × treatment group<sup>2</sup></td><td>0.005</td><td>0.087</td><td>0.951</td></tr><tr><td colspan="4"><em>Blood Pressure Residual Variance</em></td></tr><tr><td>Intercept</td><td>2.321</td><td>0.012</td><td>&lt; 0.001</td></tr><tr><td>Time</td><td>0.007</td><td>0.004</td><td>0.088</td></tr><tr><td>Treatment group</td><td>-0.033</td><td>0.014</td><td>0.018</td></tr></tbody></table>

Abbreviation: BP: Blood pressure.

model provides a better fit according to the likelihood ratio test. The variance of the residual error was heterogeneous between the subjects ($\hat{Var}(\tau_{0i}) = 0.13, se = 7e-3$) and was lower for treated patients ($\hat{\mu}_2 = -0.033, p = 0.018$). The risk of CVD events increased with age ($HR = 1.04$ for one year, $IC = [1.03; 1.05]$) and was higher for men ($HR = 1.34, IC = [1.13; 1.60]$). It was also higher for elevated blood pressure values prior to entering the study ($HR = 1.03$ per 5 mmHg, $IC = [1.00; 1.06]$). Adjusting for age, sex, treatment group, and blood pressure before inclusion, the risk of CVD increased with the current blood pressure variance ($HR = 1.04, IC = [1.01; 1.08]$). More surprisingly, the risk decreased when the current value of blood pressure increased ($HR = 0.93, IC = [0.89; 0.98]$ for 5 mmHg) or when the current slope of blood pressure increased ($HR = 0.91, IC = [0.87; 0.95]$). These negative associations are not due to the adjustment for blood pressure variance, as the same results are observed in the CVCS model (see Table S7). The risk of death from other causes was higher for older individuals ($HR = 1.05, IC = [1.03; 1.07]$ for one year) and for men ($HR = 1.68, IC = [1.15; 2.46]$), but it was not associated with the treatment group or the first measurement of blood pressure. More importantly, the death risk significantly increased with the current variance of blood pressure ($HR = 1.11, IC = [1.03; 1.19]$) and decreased when the current value ($HR = 0.90, IC = [0.808; 1.002]$ for 5 mmHg) and the current slope of blood pressure increased ($HR = 0.89, IC = [0.789; 1.002]$). As for CVD, these unexpected results regarding the current value and the current slope are also observed without adjustment on the current blood pressure variance (CVCS model in Table S7). However, these results must be interpreted cautiously considering that the study population consists only of subjects who have survived a stroke and are included in a clinical trial, thus benefiting from close monitoring.

## 4.4 | Goodness-of-Fit Assessment

To assess the fit of the time-to-event submodels, we computed for each event, the predicted cumulative hazard function at each event time by plugging the empirical Bayes estimates of the random effects in the formula for the risk function. Then we compared the mean of this predicted cumulative hazard function with its Nelson-Aalen estimator for the whole sample (Figure S2) and stratified according to sex and randomization group (Figure S3). These figures show that the joint model adequately fitted both risks and that the proportional risk assumption was valid for each categorical variable.

To evaluate the fit of the longitudinal submodel, we computed the empirical Bayes estimates of the random effects and the predicted marker values for each individual at their respective visit times. Figure S4 compares the mean predicted marker values at

<table><caption>**TABLE 5** | Covariance matrix (and standard errors) of the random effects computed using the Delta-Method.</caption><tr><td>Σ =</td><td><table><tr><td>Var(b<sub>0i</sub>)</td><td></td></tr><tr><td>Cov(b<sub>0i</sub>, b<sub>1i</sub>)Var(b<sub>1i</sub>)</td><td></td></tr><tr><td>Cov(b<sub>0i</sub>, b<sub>2i</sub>)Cov(b<sub>1i</sub>, b<sub>2i</sub>)</td><td>Var(b<sub>2i</sub>)</td></tr><tr><td>Cov(b<sub>0i</sub>, τ<sub>0i</sub>)Cov(b<sub>1i</sub>, τ<sub>0i</sub>)</td><td>Cov(b<sub>2i</sub>, τ<sub>0i</sub>)Var(τ<sub>0i</sub>)</td></tr><tr><td>Cov(b<sub>0i</sub>, τ<sub>1i</sub>)Cov(b<sub>1i</sub>, τ<sub>1i</sub>)</td><td>Cov(b<sub>2i</sub>, τ<sub>1i</sub>)Cov(τ<sub>0i</sub>, τ<sub>1i</sub>)</td><td>Var(τ<sub>1i</sub>)</td></tr></table></td><td>=</td><td><table><tr><td>214.4<sub>(6.1)</sub></td><td></td></tr><tr><td>-25.3<sub>(3.2)</sub>41.5<sub>(3.1)</sub></td><td></td></tr><tr><td>0.95<sub>(0.71)</sub> - 7.4<sub>(0.69)</sub></td><td>1.7<sub>(0.2)</sub></td></tr><tr><td>2.2<sub>(0.2)</sub> - 0.14<sub>(0.10)</sub></td><td>-0.04<sub>(0.02)</sub>0.13<sub>(7e-3)</sub></td></tr><tr><td>-0.24<sub>(0.07)</sub>0.25<sub>(0.05)</sub></td><td>-0.02<sub>(0.01)</sub> - 0.02<sub>(2e-3)</sub></td><td>0.01<sub>(1e-3)</sub></td></tr></table></td></tr></table>
each visit with the mean of the observed measurements, demonstrating that the joint model accurately captures the trajectory of blood pressure.

To highlight the impact of adding a subject-specific and time-dependent residual variance in the mixed model, we computed the individual predictions of the marker over time for some selected subjects. The predicted value of blood pressure corresponds to the conditional expectation given the random effects, defined by $\hat{\mathbb{E}}(Y_i(t)|\tilde{b}_i, \tilde{\tau}_i)$ and the prediction interval around these predicted values is given by $\hat{\mathbb{E}}(Y_i(t)|\tilde{b}_i, \tilde{\tau}_i) \pm 1.96\hat{\mathbb{V}}(Y_i(t)|\tilde{b}_i, \tilde{\tau}_i)$. For each subject, the empirical Bayes estimates of the random effects, denoted by $(\tilde{b}_i, \tilde{\tau}_i) = \underset{f(b_i, \tau_i|Y_i, T_i, \delta_i)}{\operatorname{argmax}} f(b_i, \tau_i|Y_i, T_i, \delta_i)$, corresponds to the mode of their estimated conditional posterior given the data.

They are computed by maximising $f(Y_i, T_i, \delta_i)f(b_i, \tau_i)$ with the Marquardt–Levenberg algorithm.

For some selected subjects still at risk at three years, Figure 1 presents the predicted values and their confidence intervals from the models with and without subject-specific residual variance (CVCS+V and CVCS). It shows that assuming a time-dependent and subject-specific residual variance allows a better fit of the uncertainty around the individual prediction.

## 4.5 | Risk Predictions

We compared the predictive abilities of models with and without time-dependent individual variability using AUC using a fivefold

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** | Prediction over time of the individual blood pressure and its prediction interval at 95% for 12 subjects. Model CVCS+V assumed a time-dependent subject-specific variability and the model CVCS assumes a homogeneous and constant variability. The black triangles are the observed measurements.</figcaption></figure>
cross-validation. The individual predictions of having CVD (or death from other causes) between three and five years for subjects free of any event at three years were computed using equation (4). The AUC was computed using the `timeROC` package [27]. The results are slightly better for the model with heterogeneous variability. We obtained, respectively, 0.559 (0.064) and 0.537 (0.064) for the risk of CVD and 0.624 (0.080) and 0.625 (0.080) for the risk of death. Note that between three and five years, only 52 deaths were observed, compared to 417 CVD during the same period.

To illustrate the effect of taking into account the current value of individual variance, we also computed the predicted risk of the events between three and five years for different subjects from both models, with and without time-dependent individual variability. We used the subjects selected for Figure 1 and

present their predictions obtained via the cross-validation procedure. Figures 2 and 3 show that, for both the risk of CVD and the risk of death, the predicted probability is higher with the complete model when the individual experienced the event between three and five years than with the model without the heterogeneous variability. Conversely, the predicted risk is smaller with the complete model when the individual do not experience the corresponding event.

5 | Discussion

In this work, we have proposed a new joint model with a subject-specific time-dependent variance that extends the models proposed by Gao et al. [14] and Barrett et al. [15]. Indeed,

<figure><img src="image_3.png" /><figcaption>**FIGURE 2** | Prediction of the risk of CVD between three and five years (with its 95% confidence interval indicated by dashed lines), for six patients at risk at three years, for Model CVCS+V (purple) and Model CVCS (orange). The dashed lines represent the observed time.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**FIGURE 3** | Prediction of the risk of Death between three and five years (with its 95% confidence interval indicated by dashed lines), for six patients at risk at three years, for Model CVCS+V (purple) and Model CVCS (orange). The dashed lines represent the observed time.</figcaption></figure>

this new model allows time and covariate dependent individual variance and a flexible dependence structure between the competing events and the longitudinal marker. In particular, the risk of events may depend on both the current value and the current slope of the marker, in addition to the subject-specific time-dependent standard deviation of the residual error. This is an important asset of the model given that, in most health research contexts, it is more sensible to assume that the event risk depends on the time-dependent current value or slope of the marker instead of only time-independent random effects. Moreover, accounting for competing events may be important in many clinical applications. Simulation study allows us to demonstrate the good performance of the estimation procedure and to study the impact of the choice of the number of QMC

draws for numerical integration. The model converged without bias and with good coverage rates. Moreover, the estimates of the time-to-event sub-model are quite robust to a misspecification of the marker trajectory. In addition, we provided an R-package that allows frequentist estimation with a robust estimation algorithm which had shown very good behaviour in our simulations and in a previous work with different models [20].

The analysis of the PROGRESS trial has shown that a high variability of blood pressure is associated with a high risk of CVD and death from other causes. Moreover, the individual residual variability depends on treatment group. The counterintuitive findings regarding the effects of current blood pressure value and slope are not attributable to the correlation between
individual slope and individual variance. In fact, similar results were observed in the standard joint model, which does not account for subject-specific variance. Additionally, the studied population is highly specific-consisting of individuals with a history of stroke and under tight clinical control-limiting the generalizability of these findings to other groups. Replicating this analysis in representative samples of the general population would be highly valuable for determining whether blood pressure variability could be a target for primary prevention of cardiovascular disease, mortality, or stroke.

The proposed model is flexible enough to accommodate nonlinear time trends for the marker, including polynomial forms or smooth functions based on a spline basis. Additionally, it can handle time-varying, subject-specific residual variance but can also be applied under the assumption of constant subject-specific residual variance. We have not considered more flexible time trend for the variance in the application because, as in most data sets, the number of repeated measurements in the PROGRESS trial is not large enough for estimating such a flexible model. However, the use of more flexible time trend is possible in the package.

In this work, we have supposed that the visit times were not informative and that missing measurements before the event were missing at random. In the PROGRESS clinical study, these hypotheses are quite plausible since visits were planned following a prespecified protocol and the rate of missed visits before the event was low (less than 3%). For application to observational studies, it could be useful to extend this approach to consider an informative observation process. However, such a model would require three submodels: A mixed model for the evolution of the marker, a submodel for repeated events to describe the visit process and a model for the competing events of interest. This model would rely on non-verifiable parametric assumptions and its estimation process would be much more cumbersome.

The proposed approach addressed both right censoring and left truncation, the two most common observation schemes for time-to-event data. Considering interval censoring and semi-competing events could represent a valuable enhancement. This extension would be useful when the exact time of onset of the main event is unknown (dementia for instance) and the competing event may arise after the main event (death). However, this would require modeling the three transition intensities and the interval censoring would significantly complicate the computation of the likelihood.

Such joint models with dependence on the heterogeneous variance (that can be viewed as an extension of the location-scale mixed model [13]) are of great interest to investigate the association between the variability of markers or risk factors and the risk of health events in various fields of medical research, possibly allowing to improve the prediction ability for the event. For instance, hypotheses have emerged about the link between emotional instability and the risk of psychiatric events, or the variability of glycemia and the prognosis of diabetes. Thanks to wearable devices, recent medical research studies often include frequent repeated measures of exposures or biomarkers, allowing the investigation of hypotheses regarding the variability.

## Author Contributions

Léonie Courcoul contributed to the definition of the model and design of the simulation, implemented all the codes, performed the simulation study and the real data analysis, and wrote the first draft of the manuscript. Hélène Jacqmin-Gadda and Antoine Barbieri supervised the work (model definition and implementation, simulation design and real data analysis), contributed to the interpretation of results and revised the manuscript. Christophe Tzourio and Mark Woodward made the data available, contributed to the design and interpretation of the real data analysis and revised the manuscript.

## Acknowledgments

Computer time for this article was provided by the computing facilities MCIA of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.

This work was funded by the French National Research Agency (grant ANR-21-CE36 for the project "Joint Models for Epidemiology and Clinical research").

This PhD program was supported within the framework of the PIA3 (Investment for the Future). Project reference: 17-EURE-0019.

## Conflicts of Interest

The authors declare no conflicts of interest.

## Data Availability Statement

Data may be available from contacting Professor Woodward markw@georgeinstitute.org.au, subject to legal considerations.

## References

1. d G. Pouvourville, "Coût de la Prise en Charge Des Accidents Vasculaires cérébraux en France," *Archives of Cardiovascular Diseases Supplements* 8, no. 2 (2016): 161-168, https://doi.org/10.1016/S1878-6480(16)30330-5.

2. E. Pringle, C. Phillips, L. Thijs, et al., "Systolic Blood Pressure Variability as a Risk Factor for Stroke and Cardiovascular Mortality in the Elderly Hypertensive Population," *Journal of Hypertension* 21, no. 12 (2003): 2251-2257.

3. P. M. Rothwell, S. C. Howard, E. Dolan, et al., "Prognostic Significance of Visit-To-Visit Variability, Maximum Systolic Blood Pressure, and Episodic Hypertension," *Lancet* 375, no. 9718 (2010): 895-905, https://doi.org/10.1016/S0140-6736(10)60308-X.

4. D. Shimbo, J. D. Newman, A. K. Aragaki, et al., "Association Between Annual Visit-to-Visit Blood Pressure Variability and Stroke in Postmenopausal Women," *Hypertension* 60, no. 3 (2012): 625-630, https://doi.org/10.1161/HYPERTENSIONAHA.112.193094.

5. M. H. Mehlum, K. Liestøl, S. E. Kjeldsen, et al., "Blood Pressure Variability and Risk of Cardiovascular Events and Death in Patients With Hypertension and Different Baseline Risks," *European Heart Journal* 39, no. 24 (2018): 2243-2251, https://doi.org/10.1093/eurheartj/ehx760.

6. P. K. Andersen and N. Keiding, "Interpretability and Importance of Functionals in Competing Risks and Multistate Models," *Statistics in Medicine* 31, no. 11-12 (2012): 1074-1088, https://doi.org/10.1002/sim.4385.

7. d H. Courson, L. Ferrer, A. Barbieri, et al., "Impact of Model Choice When Studying the Relationship Between Blood Pressure Variability and Risk of Stroke Recurrence," *Hypertension* 78, no. 5 (2021): 1520-1526, https://doi.org/10.1161/HYPERTENSIONAHA.120.16807.
8. R. L. Prentice, "Covariate Measurement Errors and Parameter Estimation in a Failure Time Regression Model," *Biometrika* 69, no. 2 (1982): 331-342, https://doi.org/10.1093/biomet/69.2.331.

9. D. Commenges and H. Jacqmin-Gadda, *Dynamical Biostatistical Models* (CRC Press, 2015).

10. D. Rizopoulos, *Joint Models for Longitudinal and Time-To-Event Data: With Applications in R* (CRC Press, 2012).

11. A. A. Tsiatis and M. Davidian, "Joint Modeling of Longitudinal and Time-To-Event Data: An Overview," *Statistica Sinica* 14, no. 3 (2004): 809-834.

12. R. Henderson, P. Diggle, and A. Dobson, "Joint Modelling of Longitudinal Measurements and Event Time Data," *Biostatistics* 1, no. 4 (2000): 465-480, https://doi.org/10.1093/biostatistics/1.4.465.

13. D. Hedeker and R. Nordgren, "MIXREGLS: A Program for Mixed-Effects Location Scale Analysis," *Journal of Statistical Software* 52, no. 12 (1999): 1-38, https://doi.org/10.18637/jss.v052.i12.

14. F. Gao, J. P. Miller, C. Xiong, J. A. Beiser, and M. Gordon, "The Ocular Hypertension Treatment Study (OHTS) Group A Joint-Modeling Approach to Assess the Impact of Biomarker Variability on the Risk of Developing Clinical Outcome," *Statistical Methods & Applications* 20, no. 1 (2011): 83-100, https://doi.org/10.1007/s10260-010-0150-z.

15. J. K. Barrett, R. Huille, R. Parker, Y. Yano, and M. Griswold, "Estimating the Association Between Blood Pressure Variability and Cardiovascular Disease: An Application Using the ARIC Study," *Statistics in Medicine* 38, no. 10 (2019): 1855-1868, https://doi.org/10.1002/sim.8074.

16. S. Mac Mahon, S. Neal, C. Tzourio, et al., "Randomised Trial of a Perindopril-Based Blood-Pressure-Lowering Regimen Among 6105 Individuals With Previous Stroke or Transient Ischaemic Attack," *Lancet* 358, no. 9287 (2001): 1033-1041, https://doi.org/10.1016/S0140-6736(01)06178-5.

17. R. Betensky and M. M, "Recognizing the Problem of Delayed Entry in Time-To-Event Studies: Better Late Than Never for Clinical Neuroscientists," *Annals of Neurology* 78, no. 6 (2015): 839-844, https://doi.org/10.1002/ana.24538.

18. J. Pan and R. Thompson, "Quasi-Monte Carlo Estimation in Generalized Linear Mixed Models," *Computational Statistics & Data Analysis* 51, no. 12 (2007): 5765-5775, https://doi.org/10.1016/j.csda.2006.10.003.

19. P. Gonnet, "A Review of Error Estimation in Adaptive Quadrature," *ACM Computing Surveys* 44, no. 4 (2012): 22:1-22:36, https://doi.org/10.1145/2333112.2333117.

20. V. Philipps, B. P. Hejblum, M. Prague, D. Commenges, and C. Proust-Lima, "Robust and Efficient Optimization Using a Marquardt-Levenberg Algorithm With R Package marqLevAlg," *R Journal* 13 (2021): 273, https://doi.org/10.32614/RJ-2021-089.

21. K. Levenberg, "A Method for the Solution of Certain Non-Linear Problems in Least Squares," *Quarterly of Applied Mathematics* 2, no. 2 (1944): 164-168, https://doi.org/10.1090/qam/10666.

22. D. W. Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear Parameters," *Journal of the Society for Industrial and Applied Mathematics* 11, no. 2 (1963): 431-441, https://doi.org/10.1137/0111030.

23. K. Meyer and D. Houle, "Sampling Based Approximation of Confidence Intervals for Functions of Genetic Covariance Matrices," in *Proceedings of the Association for the Advancement of Animal Breeding and Genetics*, vol. 20 (2013), 523-526, http://www.aaabg.org/aaabghome/AAABG20-papers/meyer20523.pdf.

24. T. P. Morris, I. R. White, and M. J. Crowther, "Using Simulation Studies to Evaluate Statistical Methods," *Statistics in Medicine* 38, no. 11 (2019): 2074-2102, https://doi.org/10.1002/sim.8086.

25. R. P. Brent, *Algorithms for Minimization Without Derivatives* (Prentice-Hall, 1973).

26. T. G. Pickering, J. E. Hall, L. J. Appel, et al., "Recommendations for Blood Pressure Measurement in Humans and Experimental Animals: Part 1: Blood Pressure Measurement in Humans: A Statement for Professionals From the Subcommittee of Professional and Public Education of the American Heart Association Council on High Blood Pressure Research," *Circulation* 111, no. 5 (2005): 697-716, https://doi.org/10.1161/01.CIR.0000154900.76284.F6.

27. P. Blanche, C. Proust-Lima, L. Loubère, C. Berr, J. F. Dartigues, and H. Jacqmin-Gadda, "Quantifying and Comparing Dynamic Predictive Accuracy of Joint Models for Longitudinal Marker and Time-To-Event in Presence of Censoring and Competing Risks," *Biometrics* 71, no. 1 (2015): 102-113, https://doi.org/10.1111/biom.12232.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section. **Data S1.** Supporting Information.