A method to construct a points system to predict cardiovascular disease considering repeated measures of risk factors

Cox model with time-dependent variables

Let T be a non-negative random variable denoting the observed failure time, which is the minimum value of the true event time T^∗ and the censoring time C (non-informative right censoring). In other words, T = min(T^∗, C). In addition, we define δ as the event indicator, which takes the value 1 if T^∗ ≤ C and 0 otherwise. On the other hand, let W be the vector of baseline covariates and Y(t) the vector of time-dependent covariates, assuming a defined value for t ≥ 0. With these data, the Cox model with time-dependent variables takes the following form (risk function): ${\displaystyle h\left(t|w,y\left(t\right)\right)=h_0\left(t\right)exp\left\{{\gamma }^{T}w+{\alpha }^{T}y\left(t\right)\right\},}$ where h₀(t) is the baseline risk function, and γ and α are the vectors of the regression coefficients for the baseline and time-dependent covariates, respectively (Andersen & Gill, 1982).

The estimation of the model parameters is based on the partial likelihood function (Andersen & Gill, 1982). On the other hand, we have to corroborate whether the functional form of the covariates in the model is linear. This should be performed using graphical methods (Martingale residuals against the covariate of interest). Finally, we have to assess whether the model fits the data well, through the analysis of the Cox-Snell residuals (graphical test).

The classical Cox regression model (with no time-varying covariates), deletes α and y(t) from the above expression. Furthermore, the model has to verify the following condition (proportional hazard assumption): ${\displaystyle log\left(\frac{h\left(t|w\right)}{h_0\left(t\right)}\right)={\gamma }^{T}w.}$

Points system in the Framingham Heart Study

We summarize the steps of the method developed by the Framingham Heart Study to adapt a Cox regression model with p covariates to risk charts (Sullivan, Massaro & D’Agostino, 2004):

1. Estimate the parameters of the model: $\hat{\gamma }$.

2. Organize the risk factors into categories and determine reference values:

   1. Continuous risk factor (e.g., age): set up contiguous classes and determine reference values for each. Example for age: 18–30 [24], 30–39 [34.5], 40–49 [44.5], 50–59 [54.5], 60–69 [64.5] and ≥70 years [74.5]. In brackets is the reference value. The Framingham Heart Study researchers recommend mid-points as acceptable reference values, and for the first and last class the mean between the extreme value and 1st (first class) or 99th percentiles (last class).

   2. Binary risk factors (e.g., gender, 0 for female and 1 for male): the reference value is again either 0 or 1.

      Let W_ij denote the reference value for the category j and the risk factor i, where i = 1, …, p and j = 1, …, c_i (total number of categories for the risk factor i).

3. Determine the referent risk factor profile: the base category will have 0 points in the scoring system and it will be denoted as W_iREF, i = 1, …p.

4. Determine how far each category is from the base category in regression units: calculate ${\hat{\gamma }}_i\cdot \left(W_{ij}-W_{iREF}\right)$, i = 1, …, p and j = 1, …, c_i.

5. Set the fixed multiplier or constant B: the number of regression units equivalent to 1 point in the points system. The Framingham Heart Study generally uses the increase in risk associated with a 5-year increase in age.

6. Determine the number of points for each of the categories of each risk factor: the closest integer number to ${\hat{\gamma }}_i\cdot \left(W_{ij}-W_{iREF}\right)∕B$.

7. Determine risks associated with point totals: $1-{{\hat{S}}_0\left(t\right)}^{exp\left\{{\sum }_{i=1}^p\left({\hat{\gamma }}_i\cdot W_{iREF}\right)+B\cdot Points-{\sum }_{i=1}^p{\hat{\gamma }}_i\cdot \widehat{{\stackrel{̄}{w}}_i}\right\}}$, where ${\hat{S}}_0\left(t\right)$ is calculated through the Kaplan–Meier estimator.

Joint models for longitudinal and time-to-event data

Using the former notation, we have the random variables vector $\left\{T,W,Y\left(T\right)\right\}$, where $Y\left(T\right)$ is only a time-dependent variable (longitudinal outcome) which has its values defined intermittently for t. In other words, for a subject $\left(i=1,\dots ,n\right)$, y(t) is only defined for $t_{ij}\left(j=1,\dots ,n_i\right)$, $y_i\left(t_{ij}\right)$, where 0 ≤ t_i1 ≤ t_i2 ≤ … ≤ t_ini. Now, we will denote as $m\left(t\right)$ the true and unobserved value of the longitudinal outcome at time t ($m_i\left(t\right)$ for the subject i). To assess the effect of $m\left(t\right)$ on the event risk, a standard option is to adjust a Cox regression model with one time-dependent covariate: ${\displaystyle h\left(t|M\left(t\right),w\right)=h_0\left(t^{\ast }\right)exp\left\{{\gamma }^{T}w+\alpha m\left(t\right)\right\},}$ where $M\left(t\right)$ for a subject i is defined as $M_i\left(t\right)=\left\{m_i\left(u\right);0\le u<t\right\}$, which denotes the history of the true unobserved longitudinal process up to time t. The other parameters in the expression follow the structure of the Cox regression model with time-dependent variables (see former section). The baseline risk function can be unspecified or can be approximated with splines or step functions (Rizopoulos, 2012).

In the above expression, we have used m(t) as the true unobserved longitudinal process. However, in our sample we have y(t); therefore, we will estimate m(t) using y(t) through a linear mixed effects model to describe the subject-specific longitudinal evolutions: ${\displaystyle \left\{\begin{array}{c}{y}_i\left(t\right)=m_i\left(t\right)+{\epsilon }_i\left(t\right)\hfill \\ \begin{array}{c}{m}_i\left(t\right)=x_{\mathbf{i}}^T\left(t\right)\beta +z_{\mathbf{i}}^T\left(t\right)b_{\mathbf{i}}\hfill \\ b_{\mathbf{i}}\sim N\left(\mathbf{0},D\right)\hfill \\ {\epsilon }_i\left(t\right)\sim N\left(0,{\sigma }^2\right)\hfill \end{array}\hfill \end{array},\right.}$ where β and b_i denote the vectors of regression coefficients for the unknown fixed-effects parameters and the random effects respectively, x_i(t) and z_i(t) denote row vectors of the design matrices for the fixed and random effects respectively, and ε_i(t) is the error term with variance σ². Finally, b_i follows a normal distribution with mean 0 and covariance matrix D, and independent of ε_i(t) (Rizopoulos, 2012).

The estimation of the parameters of the joint models is based on a maximum likelihood approach that maximizes the log-likelihood function corresponding to the joint distribution of the time-to-event and longitudinal outcomes (Rizopoulos, 2012).

Regarding the assumptions of the model, we have to assess them for both submodels (longitudinal and survival) using the residual plots. For the longitudinal part, we will plot the subject-specific residuals versus the corresponding fitted values, the Q–Q plot of the subject-specific residuals, and the marginal residuals versus the fitted values. On the other hand, for the survival part, we will plot the subject-specific fitted values for the longitudinal outcome versus the martingale residuals, and finally we will determine graphically whether the Cox-Snell residuals is a censored sample from a unit exponential distribution (Rizopoulos, 2012). Regarding the last component (random effects part) of the joint model for which we have indicated an assumption, other authors have showed that linear mixed-effects models are relatively robust to misspecification of this distribution (Verbeke & Lesaffre, 1997).

Predictions of the longitudinal biomarkers using these joint models for longitudinal and time-to-event data

Let $\left\{t_i,{\delta }_i,w_i,y_i\left(t_{ij}\right),{0\le t}_{ij}\le t_i,j=1,\dots ,n_i\right\},i=1,\dots n$ be a random sample of the random variables vector {T, Δ, W, Y}, using the former notation. A joint model has been fitted using this sample. Now, we are interested in predicting the expected value of the longitudinal outcome at time u > t for a new subject i who has a history up to the time t of the observed longitudinal marker $Y_i\left(t\right)=\left\{y_i\left(s\right);0\le s<t\right\}$: ${\displaystyle {\omega }_i\left(u|t\right)=E_Y\left\{y_i\left(u\right)|t_i^{\ast }>t,Y_i\left(t\right),w_{\mathbf{i}};\theta \right\}}$ where θ denotes the parameters’ vector of the joint model (Rizopoulos, 2011).

Rizopoulos developed a Monte Carlo approach to perform this task, based on Bayesian formulation. He obtained the following simulation scheme (Rizopoulos, 2011):

Step 1: Draw ${\theta }^{\left(l\right)}\sim N\left(\hat{\theta },\widehat{var}\left(\hat{\theta }\right)\right)$.

Step 2 : Draw $b_{\mathbf{i}}^{\left(l\right)}\sim \left\{b_{\mathbf{i}}|t_i^{\ast }>t,Y_i\left(t\right),w_{\mathbf{i}}\mathbf{;}{\theta }^{\left(l\right)}\right\}$.

Step 3: Compute ${\omega }_i^{\left(l\right)}\left(u|t\right)={\mathbf{x}}_{\mathbf{i}}^T\left(u\right){\beta }^{\left(l\right)}+{\mathbf{z}}_{\mathbf{i}}^T\left(u\right){\mathbf{b}}_{\mathbf{i}}^{\left(l\right)}$.

This scheme should be repeated L times. The estimation of the parameter is the mean (or median) of the calculated values (${\omega }_i^{\left(l\right)}\left(u|t\right),l=1,\dots L$) and the confidence interval is formed by the percentiles (95%: 2.5% and 97.5% percentiles) (Rizopoulos, 2011).

We highlight that these predictions have a dynamic nature; that is, as time progresses additional information is recorded for the patient, so the predictions can be updated using this new information.

Construction

We wish to determine the probability of having CVD with effect from a baseline situation (t = 0) up to a fixed point in time ($\tilde{t}$), given a series of risk factors measured at baseline and during this follow-up. To do this requires the following steps:

1. Adjust a Cox regression model with time-dependent variables. As we are unable to estimate a joint model with multiple longitudinal parameters (Rizopoulos, 2012), we use the classic extended Cox model (with no shared structure), which requires knowing the values of all the longitudinal parameters at any value of t. As this is not known because the parameters are recorded intermittently, we take the last value in time as a reference.

2. Use the procedure of the Framingham study to adapt the coefficients of the model obtained to a points system and determine the probabilities of CVD for each score up to the moment $\tilde{t}$. We then use these probabilities to construct risk groups that are easy for the clinician to understand (for example, in multiples of 5%) (Sullivan, Massaro & D’Agostino, 2004).

3. Adjust a joint model for longitudinal and time-to-event data for each longitudinal parameter recorded during the follow-up. This will also include all the baseline variables. These models are constructed to make predictions about the longitudinal parameters in new patients (statistical validation by simulation and potential utilization).

Statistical validation by simulation

Once the points system has been constructed, we wish to see whether the model determines the onset of CVD accurately in a different set of subjects (validation sample). In this validation sample we know the longitudinal markers up to the point t = 0 (record of cardiovascular risk factors in the clinical history which were measured before the baseline situation (t < 0)) and the value of the variables at baseline. With this information we determine the probability each subject has of experiencing an event, and we then compare this with what actually occurred; i.e., determine whether the model is valid. To determine this validity we follow these steps:

1. Determine L simulations of the longitudinal parameters at the time point $\tilde{t}$ using the models mentioned in step (3) of construction, from the history (t < 0) and the baseline variables (t = 0) (Rizopoulos, 2011). We will use these simulated values to construct a distribution of the points for each sample subject. Thus, each subject will have L values for the points variable (evaluating the points system using the simulated values and the baseline variables is sufficient), and for each lth simulation each patient will have a points score. In other words, each simulation will have a distribution of the points variable.

2. For each lth simulation adjust a classic Cox model (without time-dependent variables), using just the score obtained as the only explanatory variable. Determine the Harrell’s concordance statistic for each of these L models. These values will give us a distribution of values for this statistic, with which we calculate the mean (or the median) and the 2.5% and 97.5% percentiles (Rizopoulos, 2011). This way we construct a confidence interval for this statistic, which will indicate the discriminating capacity of the points system to determine which patients will develop CVD.

3. Calculate the median of the points distribution for each patient in the validation sample. Note that we do not use the mean as it could contain decimals and this has no sense when applying the scoring system. Using these medians, classify each patient in a risk group and compare the rate of events predicted by the points system in each group to the actual observed rate. The test used for this process will be Pearson χ² test.

The concordance statistic used has been reported to have various limitations (Lloyd-Jones, 2010). For example, it does not compare whether the estimated and observed risks are similar in the subjects. Accordingly, we have added the analysis of the differences between the expected events and the observed events, which minimises this particular problem. In addition, it is very sensitive to large hazard ratio values (≥9). Nonetheless, we have to consider that as all the variables are quantitative (not categorized), the hazard ratio values do not surpass this threshold. Accordingly, the joint analysis of the concordance index of Harrell and the differences between the expected and the observed events enables us to validate statistically by simulation of the proposed model.

Explanation of potential utilization

Once the points system has been validated statistically the clinician can then apply the system to determine the cardiovascular risk in a new patient, and take any necessary measures to reduce this risk. The healthcare professional will already have historical information about the longitudinal parameters (t < 0) and information about the baseline situation (t = 0) of the new patient. The steps to be followed by the clinician are:

1. Determine the value of each longitudinal parameter at the time $\tilde{t}$. To do this we apply the models obtained in step (3) of construction to the history and the baseline situation of the new patient, in order to determine L simulations for each longitudinal parameter, similar to what was done in the validation process. For each lth simulation we determine the score corresponding to the profile of cardiovascular risk factors obtained (simulated and baseline information values). This will give us a points distribution for the new patient.

2. Determine the median and the 2.5% and 97.5% percentiles of the points vector constructed above. The median will be the estimation of the score for the new patient and the percentiles will define the confidence interval (Rizopoulos, 2011). As each score has an associated risk, the healthcare professional will be able to know the probability of CVD at time $\tilde{t}$, together with its confidence interval. Finally, the clinician will know the values of the biological parameters at $\tilde{t}$ of the median of the points system. This way the clinician will be able to see which of these parameters has a score above normal; i.e., see the possible areas of intervention to reduce the cardiovascular risk.

3. The clinician now knows the cardiovascular risk and which parameters have a score above normal, so he or she can then design the best intervention for that patient. This presents a problem, as we need to know the value of each biological parameter at time $\tilde{t}$; i.e., the clinician knows an approximation based on simulations constructed from the patient history but does not know how the interventions will affect the cardiovascular risk.

From the previous step the clinician knows the parameters on which to act and the history of these parameters as well as the baseline situation. From these measurements the clinician can establish a realistic objective for the next patient visit at time $\widetilde{\stackrel{̃}{t}}\left(0<\widetilde{\stackrel{̃}{t}}<\tilde{t}\right)$. The clinician now inserts the desired value of the biological parameter at $\widetilde{\stackrel{̃}{t}}$ and determines its value at time $\tilde{t}$; i.e., determine L simulations for each cardiovascular risk factor using the previous models (step 3 of construction), adding a new value to the history ($\widetilde{\stackrel{̃}{t}}$).

These calculations will give the benefit of the intervention (estimation (mean or median) of the biological parameter at $\tilde{t}$) and the clinician will be able to see from the points system how the patient’s risk will be reduced.

Simulation on a data set

With the sole purpose of explaining how to use the method proposed here, we have simulated a data set upon which to apply each of the steps described above. Note that we are in fact going to simulate two data sets, one to construct the model and the other to validate it statistically via simulation. So that both sets are biologically plausible we have used estimations obtained in the Puras-GEVA cardiovascular study, which has been published in Medicine (Artigao-Ródenas et al., 2015).

Our data sets will include the following biological parameters: age (years), systolic blood pressure (SBP) (mmHg), HbA1c (%), atherogenic index, gender (male or female) and smoking (yes or no). Of these, the SBP, HbA1c and the atherogenic index will be present at baseline (t = 0) and in the follow-up for the construction sample (t > 0) or recorded in the clinical history for the statistical validation sample via simulation (t < 0). The choice to include these variables was based on the current cardiovascular risk scales (Conroy et al., 2003; National Heart, Lung, and Blood Institute, 2015), except for HbA1c, which is used instead of a diagnosis of diabetes mellitus in order to include another time-dependent parameter in the final model, in addition to which this way enables us to value the control of the diabetes mellitus (HbA1c < 6.5%) when preventing CVD.

For the main variable (time-to-CVD) we shall suppose that our cohort is used to predict CVD with a follow-up of 2 years. Note that the traditional cardiovascular risk scales use a time of 10 years (Conroy et al., 2003; National Heart, Lung, and Blood Institute, 2015). We have used this lower value because we are going to make predictions for the longitudinal parameters with effect from the baseline situation (t = 0) up to the prediction time and if we take a prediction value of 10 years the predictions for the longitudinal parameters will vary greatly and not allow us to make precise predictions about which patient will develop CVD, which would negate the usefulness of the method proposed here. Nevertheless, the fact that the predictions for the longitudinal parameters have a dynamic character (see Predictions of the longitudinal biomarkers using these joint models for longitudinal and time-to-event data) enables us to determine the risk at 2 years with greater precision whenever the patient attends the office of the healthcare professional. Note that the method proposed here has been developed for a theoretical time period $\tilde{t}$ but it can be applied for any time period. Nonetheless, generally speaking the longitudinal parameters would vary more over longer time periods, though this clearly depends on the nature of the data, both at the individual level and the population level (Rizopoulos, 2012).

The work used for our simulated data set developed and validated a predictive model of CVD (angina of any kind, myocardial infarction, stroke, peripheral arterial disease of the lower limbs, or death from CVD), to enable calculation of risk in the short, medium and long term (the risk associated with each score was calculated every 2 years up to a maximum of 14) in the general population (Artigao-Ródenas et al., 2015). Table 4 of this scoring system shows the importance of this question. For example, a patient with a score of 9 points has a probability of CVD at 2 years of 0.67%, whereas at 10 years this rises to 5.16% (Artigao-Ródenas et al., 2015: Table 4). If we regularly calculate the 2-year risk of CVD for our patient and the score remains the same then no new therapeutic action will be taken (risk < 1%), whereas if we only calculate the risk once every 10 years we will take aggressive therapeutic measures when the patient first attends the office, as the score will correspond to a cut point defined as high in the SCORE project (5% → one in 20 patients) (Conroy et al., 2003). We see, then, that a regular short-term prediction could lead to a change in the therapeutic decisions regarding prevention of CVD, provided of course that the possibility exists of calculating the risk regularly. As the risk table given in the Puras-GEVA study includes predictions for 4, 6, 8, 10, 12 and 14 years, we selected the lowest cut point because if we had to make predictions for a longer time the dispersion could have increased (Rizopoulos, 2012). This is why we chose this cut point of 2 years for the simulation.

The longitudinal follow-up measurements (construction sample) assumed that the patient attends the physician’s office once every 3 months for measurements of SBP, HbA1c and the atherogenic index. This is done until the end of the follow-up for each patient. The statistical validation sample using simulation supposes that there is a certain probability of having records in the clinical history of all the longitudinal parameters every 3 months for 5 years retrospectively (t < 0). The probability is different for each of the visits and will depend on each patient. In other words, we will have intermittent measurements of all these parameters from t = − 5 years to t = 0.

The Supplemental Information (S1) details all the mathematical formulae used to construct our data sets, always based on the Puras-GEVA study (Artigao-Ródenas et al., 2015). The simulation was done using R 2.13.2 and IBMS SPSS Statistics 19.

One could think that by managing a shorter time period of just 2 years there would be no variability in the cardiovascular risk factors. However, in S1 we can see that the models used show a temporal variability in the risk factors. If there were no variability in the factors, the models would contain the constant with a very small random error. In other words using this prediction time makes sense.

We decided to use a simulated data set as we did not have available any data set with real data. This way of explaining a new method has already been used by others working with joint models, as the only objective of the simulated data set is to explain how to apply the new method (Faucett & Thomas, 1996; Henderson, Diggle & Dobson, 2000; Wang & Taylor, 2001; Brown, Ibrahim & DeGruttola, 2005; Zeng & Cai, 2005; Vonesh, Greene & Schluchter, 2006; Rizopoulos & Ghosh, 2011).

Construction of the model

The parameters of the Cox model with time-dependent variables are shown in Table 1, and its adaptation to the points system with a prediction time of 2 years is reflected in Fig. 1. Table 2 shows the joint models for the longitudinal parameters. To avoid computational cost simple models were used: (1) linear equation for the survival part with all the predictors included (age, gender, smoking, and longitudinal marker) and (2) mixed linear model with polynomial degree 1 at (1, t), in both the fixed and the random parts. The baseline risk function was defined piecewise.

Statistical validation by simulation

Explanation of potential utilization

A new patient arrives at our office with the following characteristics: male, 83 years old, non-smoker, and taking pharmacological medication (one antihypertensive drug and one oral antidiabetic agent) and non-pharmacological measures (diet and exercise). His history of cardiovascular risk factors is available (Table 3).

Application of the new model gives a histogram of the cardiovascular risk score obtained for this patient (Fig. 3). This chart shows a high cardiovascular risk, as most of the simulations have around 16 points. The estimation of the score was 16 (95% CI [15–17]). The median score corresponded to a SBP of 160 mmHg, HbA1c of 5.0% and an atherogenic index of 6.76. Bearing in mind that the model contains factors upon which it is not possible to act (gender and age) that give the patient a minimum of 13 points, we should consider strategies to help the patient not to score in the other categories on the scale (Fig. 1).

The clinician can now see that if the patient complies with a series of interventions (pharmacological (add two antihypertensive drug → −20 mmHg; prescribe a statin →−40% atherogenic index) and non-pharmacological (reduce salt in the diet →−5 mmHg)), his longitudinal parameters after 3 months would be: SBP 120 mmHg (145 – 2 × 10 – 5 = 120 mmHg), atherogenic index 3.10 (5.17 – 40% = 3.10), and HbA1c 4.9% (same value because no intervention was done). Applying the model using the new information gives the cardiovascular risk at 2 years (Fig. 4). The estimation of the score is 15 (95% CI [14–15]) and the values that provide a median score are: SBP 124 mmHg, atherogenic index 4.85, and HbA1c 5.0%. Thus, the risk is reduced, as now the patient has 15 points (Fig. 1 and S3).