Estimating the association between blood pressure variability and cardiovascular disease: An application using the ARIC Study

Systolic blood pressure (SBP) is universally recognised as an important risk factor for cardiovascular disease (CVD) and is routinely included in risk scores for CVD risk prediction [3, 8]. The prognostic value of SBP is primarily based on the mean of measurements over multiple outpatient visits, whereas substantially less attention has been given to the variability of SBP across visits (i.e., visit-to-visit SBP variability). However, there is increasing evidence of associations between greater visit-to-visit SBP variability and CVD outcomes in community-based and clinical settings [11, 13, 16, 17]. Therefore, in addition to assessing mean clinic SBP levels over time, measurements of visit-to-visit SBP variability may improve the accuracy of CVD risk prediction, which is crucial for the optimization of patient care.

The standard deviation (SD), coefficient of variation (CV), average real variability (ARV), and SBP variability independent of the mean (VIM) across multiple visits have been widely used to quantify visit-to-visit SBP variability [11, 13, 16, 17, 25]. However, these measures are estimated on a person-by-person basis, and are therefore subject to considerable measurement error. This measurement error causes regression dilution bias in the estimated association between visit-to-visit SBP variability and CVD [9].

Mixed effects models have been proposed which allow the within-individual variability to differ between individuals. For example, Hedeker et al. introduced the mixed location scale model, where the within-individual variances may be assumed to depend on an additional random effect, as well as time-constant and time-varying covariates [6, 7]. More recently Goldstein et al. have explored the mixed location scale model using a Bayesian framework in the context of growth curve models[4].

We consider extensions of the mixed location scale model to estimate the association of SBP variability with CVD. Repeated SBP measurements are modelled using a mixed effects model with a random component on the within-individual standard deviation, thus allowing the borrowing of information across individuals. The model allows for each individual to have different within-individual variation, which describes the variability in their SBP measurements. We propose two methods for estimating the association of this SBP variability with CVD. The first is a two-stage method, where in stage one we obtain estimates of the usual level of SBP and the SBP variability for each individual from the mixed effects model, and in stage two these are included as covariates in a survival model for the time to the first CVD event. The second method is a joint model where the repeated SBP measurements and the time to first CVD event are modelled simultaneously using shared random effects [15].

There are advantages and disadvantages in the use of two-stage versus joint models. The heavy computational burden associated with simultaneous analysis of repeated measurements and time-to-event data makes the two-stage method a more feasible option for large data sets. The two-stage method could also be readily incorporated into a landmarking approach for dynamic risk prediction [23], for example in the dynamic model for CVD risk prediction recently proposed by Paige et al [12]. However, two-stage methods for the analysis of repeated measurements and time-to-event data have been shown to lead to bias in estimated covariate effects compared to joint models [21]. This bias is due to two factors: (1) regression dilution bias caused by residual measurement error in the estimates of random effects from the mixed effects model and (2) bias caused by informative truncation of the repeated measurements by the event of interest. In the two-stage approach this second source of bias can be removed by separating longitudinal follow-up from survival follow-up if data allows, as we describe in section 2.

The outline of the paper is as follows. In section 2 we describe the methods, including the naive method that has predominantly been used in the CVD literature, and our proposed two-stage and joint model methods. In section 3 we compare the methods using simulation studies and explore the accuracy of all methods in a variety of scenarios. We apply the two-stage and joint models to data from the Atherosclerosis Risk in Communities (ARIC) Study in section 4, and we present a final discussion in section 5.

The association of SBP with CVD has previously been investigated using mixed effects models to analyse repeated measurements of SBP [20, 1, 24]. In this section we first present the naive methods for estimating SBP variability that have been commonly used in the literature, and then explore how the mixed effects model can be extended to allow SBP variability to differ between individuals. We propose a two-stage method and a joint modelling method for estimating the association between SBP variability and CVD events. Although we focus here on the association between SBP variability and CVD, the methods and results presented would be relevant for any application where the association between the variability in a longitudinal outcome and a time-to-event is of interest. We consider estimating the association of the variability of a longitudinal outcome with a time-to-event, adjusted for the ‘usual level’ of the longitudinal outcome. This ‘usual level’ could be the average value, which may be assumed constant in time, or the current value or baseline value if time-dependence is to be taken into account.

We use simulation studies to explore the bias in parameter estimates due to imprecise estimation of SBP variability, and due to ignoring correlation between SBP and SBP variability in the mixed effects model. We consider two scenarios. In Scenario 1 each individual has the same number of repeated measurements, thus bias in the two-stage approaches arises solely from measurement error in the estimates of the usual level and variability of SBP (c.f. Figure 1b). In Scenario 2 the repeated measurement process is truncated by the event times, so that bias in the two-stage approaches arises from both measurement error and informative truncation of the repeated measurements (c.f. Figure 1a). We generated simulated data sets to explore the effect of the number of measurements per individual, the value of the association parameters and the extent of the correlation between the usual level and the variability.

We take a simple model for the repeated measurements with a random intercept only and random residual error standard deviation

where the residual error $\epsilon_{ij}\sim N(0,\sigma_{i}^{2})$ and we take a bivariate normal distribution for $b_{0i}$ and $\log\sigma_{i}$,

Here $b_{0i}$ is the usual level and $\sigma_{i}$ represents the variability in the longitudinal outcome. Based on results from the ARIC study, we take $\mu_{0}=120$, $\tau_{0}=15$, $\mu_{\sigma}=2$ and $\tau_{\sigma}=0.5$. Event times are drawn from a Weibull distribution with shape parameter $k$ and scale parameter $\lambda$. A loglinear effect of the covariates on $\lambda$ gives a proportional hazards model,

We introduce administrative censoring at 20 years. To maintain an event rate of 20% before censoring we take $\gamma_{0}=-10.26$ and $k=2$ with default association parameters $\alpha_{0}=0.02$ and $\alpha_{\sigma}=0.05$, similar to those found for the ARIC data set. For Scenario 1 we used all repeated measurements for all individuals. For Scenario 2 we assumed repeated measurements were taken equidistantly between baseline and 18 years follow-up and discarded all repeated measurements following event times. In each simulation set-up we generated 1000 data sets, each consisting of 1500 individuals.

We analysed each generated data set using the two-stage method from section 2.3.1, with (i) naive and (ii) mixed model estimates of SBP usual levels and variabilities, and the joint models from section 2.3.2. For the mixed effects and joint models we fitted models ignoring the correlation between SBP and SBP variability (mixed effects model LMM1 and joint model JM1 used equations (3) and (4)), and allowing for the correlation (mixed effects model LMM2 and joint model JM2 used equations (3) and (5)). The joint models were fitted using 15 time intervals for the baseline hazard. Model convergence was checked using the Gelman-Rubin statistic as modified by Brooks and Gelman [2]. We also fitted the survival models using the true values of $b_{0i}$ and $\sigma_{i}$ as covariates, representing the best achievable results from the time-to-event model when all covariates are known precisely. For all scenarios we calculated the mean and standard deviation of the estimated log hazard ratios and the root mean squared errors (RMSE) and coverage probabilities at the 95% level.

Table LABEL:table:n shows results for $n=4$, $n=7$ and $n=10$ measurements per individual for Scenario 1. With 4 measurements per individual using the naive method leads to a slight bias in the log hazard ratio (logHR) for the usual level, but substantial negative bias and associated loss of coverage in the logHR for variability, which is most likely caused by regression dilution due to imprecise measurement of longitudinal variability. Using the mixed effects model LMM1 leads to a positive bias in both the usual level logHR and the variability logHR, as we would expect from the arguments given in Section1 of the Supplementary Materials, because correlation between these quantities was not accounted for in their estimation. Using method LMM2, which accounts for the correlation, results in minimal bias, even with only 4 measurements per individual. For both two-stage models the coverage is a little below the nominal level, but the RMSE is larger for the variability logHR in LMM2 due to greater variation in effect estimates. Results for the joint models are similar to the two-stage models in terms of the bias, but coverage is closer to the nominal levels. Results for n=7 and n=10 are similar, with the naive model giving reduced but still substantial bias in the variability logHR with a greater number of measurements.

Table LABEL:table:ntrunc shows results for $n=4$, $n=7$ and $n=10$ measurements per individual for Scenario 2. The two-stage model LMM2 now also gives negative bias in estimating the variability logHR due to informative truncation of the repeated measurements by the CVD event. The joint model JM2, which accounts for this informative truncation, gives consistent results. Our results suggest that the bias incurred in the two-stage approach by informative truncation is more considerable here than the bias incurred by error in the estimates of SBP variability from the LMM models.

In Section 2 of the Supplementary Materials we present further results investigating performance of the models for different levels of association between the longitudinal trajectories and the CVD event (Supplementary Table 1) and for different levels of correlation between the usual level and the variability (Supplementary Table 2). All results are given for datasets with $n=4$ measurements. In brief, all methods performed well with no association between times-to-events and longitudinal variability. But when the associations between the time-to-event and the usual level and variability of the longitudinal trajectories was substantially stronger, all the methods struggled to give consistent results with only 4 repeated measurements. The joint models performed the best in this case with the least bias and highest coverage (Supplementary Table 1). When there is negative correlation between usual levels and variabilities the direction of bias is reversed. Bias in the usual level logHR using methods LMM1 and JM1 increases with increasing $\rho$. But for the variability logHR the pattern is less clear because of the interplay between regression dilution bias and the bias incurred by ignoring correlation in the mixed effects model (Supplementary Table 2).

In summary, the method which overall resulted in the least bias was JM2. For parameter values similar to those observed in the data example method LMM2 performed well. But for all methods biases were observed with strong association parameters with $n=4$ measurements per individual, so more measurements per individual would be required in this scenario. The interplay between biases due to multiple causes is complex, and depends on the correlation parameter $\rho$. The RMSE, however, is generally slightly higher for LMM2 than for LMM1 and for JM2 than for JM1, suggesting a trade-off between the bias caused by ignoring the correlation for models LMM1 and JM1 and the higher variance induced by the increased complexity of models LMM2 and JM2.

We illustrate our methods using data from the ARIC study [22]. Briefly, 15,792 mostly black and white adults aged 45-64 years were enrolled into the ARIC study between 1987 and 1989 via probability sampling from 4 U.S. communities: Washington County, Maryland; Forsyth County, North Carolina; Minneapolis, Minnesota, suburbs; and Jackson, Mississippi. Participants underwent five examinations during 25-years of follow-up (i.e., Visit 1, 2, 3, 4, and 5 examinations), with an annual contact by telephone. In the current analysis, we used SBP measurements from Visit 1 (1987-1989) through Visit 4 (1996-1998).

We analysed both the full ARIC dataset and a reduced dataset where longitudinal follow-up was truncated at the last systolic SBP measurement taken at Visit 4 (at $t_{sep}=11.9$ years from baseline) and the time origin for the survival follow-up was taken at the same time-point (Figure 1b). For the reduced dataset, to have enough measurements to estimate the standard deviation of repeated measurements using the naive method we restricted to individuals with three or four non-missing measurements recorded and also to individuals who did not experience an event before $t_{sep}$. In total, our full data set consisted of 13,161 individuals and our reduced data set consisted of 10,019 individuals. Supplementary Table 3 shows the baseline characteristics of the reduced and full datasets.

We analysed both the full and reduced data-sets using two-stage and joint models. For the two-stage approaches we used the repeated measurements of SBP to estimate the usual SBP level and the SBP variability using the naive method and the linear mixed effects models without and with accounting for the correlation between SBP and SBP variability (we again label these models LMM1 and LMM2 respectively). We also fitted joint models both without and with correlation between SBP and SBP variability (again JM1 and JM2 respectively). We applied the mixed effects models/sub-models in four ways

1. 1.

   With a random intercept only, not adjusting for baseline CVD risk factors

2. 2.

   Including a slope term in the model as both a fixed and random effect. The covariance matrix between the random effects and the residual log-standard deviation was taken to be unstructured, allowing for correlations between the random slopes and the random intercept and between the random slopes and the SBP variability.

3. 3.

   Adjusting the random intercept mixed effects model for baseline CVD risk factors (age, sex, diabetes status, smoking status, baseline total cholesterol and baseline HDL cholesterol)

4. 4.

   Adjusting the random intercept and slope mixed effects model for baseline CVD risk factors.

Time-to-event models were adjusted for age, diabetes status, smoking status, baseline total cholesterol, baseline HDL cholesterol and sex.

For the Bayesian estimation we used diffuse uniform prior distributions $U[0,100]$ for standard deviations, uniform $U[-1,1]$ prior distributions for correlation parameters and diffuse normal prior distributions $N(0,100^{2})$ for all other parameters. Priors were specified for the bivariate and trivariate normal distributions by expressing them as two and three conditional univariate normal distributions respectively. We used a burn-in of 1000 MCMC updates for the mixed effects models and 2000 MCMC updates for the joint models. We calculated results from 1000 sampled updates, and checked convergence using the Gelman-Rubin statistic as modified by Brooks and Gelman [2]. For the joint models there was a high degree of auto-correlation for the usual SBP and SBP variability hazard ratios, so MCMC chains were thinned by 4 updates.

Figure LABEL:fig:hist shows histograms of the estimated SBP usual level and SBP variability using the naive method. Also plotted are fitted normal and log-normal distributions respectively. The plots suggest that the assumption of a normal distribution for the random intercepts and a log-normal distribution for the residual SDs is appropriate. The naive usual level estimates are plotted against the naive variability estimates in Figure LABEL:fig:cor. There is a positive correlation with those with higher SBP tending to have greater SBP variability. The correlation between the naive estimates is 0.42, but the true correlation without measurement error is likely to be higher.

Estimates of usual levels and variabilities from the mixed effects models are compared with the naive estimates in Figure LABEL:fig:naivevslmm. The plotted lines are the lines of agreement between the two methods. Note that there is more difference between the methods for SBP variability estimates than for usual level estimates, with lower SBP variabilities being underestimated by the naive method and higher variabilities being overestimated.

Results from all models applied to the reduced ARIC data set are shown in table LABEL:table:ARIC. As expected, the joint models took considerably longer to run than the two-stage models with runtimes of between 17 and 20 hours compared to between 4 and 10 minutes respectively. All methods find an increased risk of CVD events with higher SBP variability. The methods which include random slopes find no evidence for an association between SBP gradient and CVD risk after adjusting for the usual level of SBP and SBP variability.

The estimated association of SBP variability with CVD events is smaller for the naive model, suggesting regression dilution bias in this parameter estimate. Method LMM1, which does not account for the correlation between the SBP usual level and variability, gives a higher estimate of the association with the usual level of SBP than method LMM2, which does account for the correlation, as would be expected from the arguments given in Section 1 of the Supplementary Materials. The estimate of the association with SBP variability is generally lower for LMM2 than LMM1, suggesting a slight bias in the LMM1 estimates. The standard errors are, however, considerably larger for LMM2 than for LMM1. There is little difference between the results from the joint models and the results from the two-stage LMM models. Adjusting for other CVD risk factors in the mixed effects model makes little difference to the estimated association with SBP usual level and variability.

Results from all models applied to the full ARIC data set are shown in table LABEL:table:ARICfull. For the models without random slopes the logHRs for the SBP usual level and SBP variability are similar to those from the truncated dataset. Standard errors are slightly smaller reflecting the increased sample size. But for the models with random slopes the logHRs for the effect of slope on the hazard of CVD are very different for the full dataset, with more positive slopes now leading to a lower risk of a CVD event when adjusted for the SBP usual level and SBP variability. This is counterintuitive to what we would expect; as higher SBP is associated with greater risk of CVD we would expect those with increasing SBP over time to also have a higher risk.

We have proposed two-stage and joint modelling methods to estimate the association between visit-to-visit SBP variability and CVD, both of which use extensions of standard mixed effects models to allow different individuals to have different SBP variabilities. Both methods reduce regression dilution bias in the estimated hazard ratio compared to naive methods that have previously been used where SBP variability is estimated on an individual-by individual basis. In addition, the joint modelling methods allow for informative truncation of the repeated measurement process by the event of interest. In practice the method of choice will depend on the likely extent of the biases incurred by the two-stage approach vs the computational tractability of the joint modelling approach.

An important assumption of the two-stage method is that truncation of the observation process is non-informative, otherwise the number of observations may depend on the underlying risk level, which may lead to bias. Informative truncation can be avoided in the two-stage approach if distinct periods of follow-up for longitudinal and survival data can be defined, but this may only be possible at the cost of discarding information from individuals with events during longitudinal follow-up and repeated measurements taken during survival follow-up. In the clinical literature analyses have been conducted both using survival follow-up recorded subsequently to the repeated observations of SBP and using survival follow-up recorded concurrently with SBP measurements. As an example of the former, Rothwell et al explored the use of 7 compared with 10 SBP measurements, but commented that when 10 SBP measurements were used survival follow-up was shorter[16].

Another important assumption of our method is that the measurement schedule is non-informative, which will not be satisfied in some studies. Some clinical studies investigating the association between SBP variability and CVD have used data which has not been collected for research purposes, such as data from electronic health records. For example, Gosmanova et al used data from around 3 million US veterans receiving healthcare from the US Veterans Health Administration[5] and Hippisley-Cox et al used SBP variability as a risk predictor in the most recent QRISK tool for predicting CVD based on around 8 million individuals attending general practices in the UK[8]. In such cases additional bias may be incurred by informative observation of the repeated measurements, where the act of measurement is more (or less) likely in individuals at greater risk of disease. In these cases it may be necessary to extend the joint model to include an additional sub-model for the recurrent event process of an SBP measurement being taken. However, for cohort studies with a pre-specified measurement schedule, such as the ARIC study analyse in this paper, such additional model complexity is unlikely to be warranted. We therefore leave this investigation for future research.

Results from our simulation studies suggest that the best model giving the least bias and coverage closest to nominal levels is the joint model which allows for correlation between the SBP usual level and the SBP variability. The naive method gives considerable bias in the association between SBP variability and CVD, as expected. Results previously reported in the CVD literature using the naive method could be subject to considerable bias in both the estimated association and its standard error. Our simulation results also indicate there may be a complex interplay between different sources of bias. While the use of mixed effects models and joint models can reduce regression dilution bias, model misspecification, such as not accounting for correlation between these variables, may introduce additional sources of bias.

Our analyses of the ARIC data found evidence of a positive association between higher SBP variability and the risk of CVD events. For this dataset results between the two-stage LMM approach and the joint modelling approach were similar, suggesting that the biases incurred by use of the two-stage approach were minimal in this example. However, use of the joint model comes at a cost of considerably longer computational time. This might be reduced by running multiple MCMC chains in parallel. If computational time is an issue then a possible modelling strategy would be to use the two-stage method for model selection and the joint model for inference. If there was also clinical interest in identifying explanatory variables associated with the longitudinal variability then the mixed effects models could be extended to incorporate a linear predictor in the mean of the $\log\sigma_{i}$ distribution, $\mu_{\sigma,i}=\zeta^{T}X^{\sigma}_{i}$. Similar models have previously been considered by Hedeker et al[6, 7] and Goldstein et al[4].

We thank Angela Wood for helpful discussions. JKB was funded by the Medical Research Council grant numbers G0902100 and MR/K014811/1 and unit programme number MC_UU_00002/5. RP was funded by Medical Research Council grant number MR/N027485/1. This work was supported by the UK Medical Research Council (G0800270), British Heart Foundation (SP/09/002), UK National Institute for Health Research Cambridge Biomedical Research Centre, European Research Council (268834UK). The Atherosclerosis Risk in Communities study has been funded in whole or in part with Federal funds from the National Heart, Lung, and Blood Institute, National Institutes of Health, Department of Health and Human Services, under Contract nos. (HHSN268201700001I, HHSN268201700003I, HHSN268201700005I, HHSN268201700004I, HHSN2682017000021). The authors thank the staff and participants of the ARIC study for their important contributions.

Consider a linear regression of outcome $y$ on correlated covariates $x_{1}$ and $x_{2}$,

where $\varepsilon\sim N(0,\sigma_{y}^{2})$ is the residual error in $y$. Now let

where $\lambda$ and $x_{2}$ are uncorrelated, $\textrm{Var}(x_{1})=\sigma_{1}^{2}$, $\textrm{Var}(x_{2})=\sigma_{2}^{2}$ and $\rho$ is the correlation between $x_{1}$ and $x_{2}$. Re-writing the linear predictor,

As $\lambda$ and $x_{2}$ are uncorrelated we find

If $x_{1}$ and $x_{2}$ are estimated ignoring the correlation between them, effectively setting $\rho=0$, we would therefore expect to observe a positive bias in $\hat{\beta}_{2}$ when $\rho$ and $\hat{\beta}_{1}$ are both positive, because the second term in equation (7) is ignored. By symmetry, we would also expect to observe a positive bias in $\hat{\beta}_{1}$ for positive $\rho$ and positive $\hat{\beta_{2}}$.

In this section we present the results of further simulation studies investigating performance of the models for different levels of association between the longitudinal trajectories and the CVD event and for different levels of correlation between the usual level and the variability for Scenario 1.

Table 1 shows results for different values of the true logHR’s $\alpha_{0}$ and $\alpha_{\sigma}$. All results are given for datasets with $n=4$ measurements per individual. When $\alpha_{\sigma}=0$, i.e. there is no association between the longitudinal variability and the time-to-event, the naive method has negative bias in the usual level logHR due to regression dilution, but there is no bias in the variability logHR. By contrast, for methods LMM1 and JM1 there is slight positive bias in the variability logHR because the correlation between the usual level and the variability has been ignored. In this case, therefore, using the mixed effects model may lead to increased bias compared to the naive method because the model is misspecified. Results for larger values of the logHRs show substantial bias in both logHRs for the naive and LMM methods, with the naive method performing the worst. The joint models give little bias in the logHR for the usual level, but some negative bias for the variability logHR. Coverage probabilities are very much lower than the nominal 95% level, with those of the joint models being substantially closer than all other methods. These results suggest that for substantially larger effect sizes, or equivalently for substantially less variation in the usual levels or variabilities in the population, all methods require a greater number of measurements per individual to sufficiently reduce estimation bias.

Results for various values of the correlation parameter $\rho$ are shown in Table 2. Again, each individual has $n=4$ measurements. When there is negative correlation we find negative bias in both association parameters for the naive, LMM1 and JM1 methods, as we would expect from the arguments given in Section 1. Bias in the usual level logHR using the methods LMM1 and JM1, which ignore the correlation, increases with increasing $\rho$. But for the variability logHR the pattern is less clear because of the interplay between regression dilution bias and the bias incurred by ignoring correlation in the mixed effects model. For models LMM2 and JM2, which account for the correlation, we found minimal bias in estimated effects and coverage close to nominal values in all cases.

Table 3 shows baseline characteristics of individuals in the reduced and full ARIC datasets at the study baseline. The age and cholesterol levels are similar between the two datasets, but the reduced dataset has a lower proportion of males, smokers and diabetics than the full dataset.