Joint modelling of longitudinal processes and time-to-event outcomes in heart failure: systematic review and exemplar examining the relationship between serum digoxin levels and mortality

Abstract

Background

Joint modelling combines two or more statistical models to reduce bias and increase efficiency. As the use of joint modelling increases it is important to understand how and why it is being applied to heart failure research.

Methods

A systematic review of major medical databases of studies which used joint modelling within heart failure alongside an exemplar; joint modelling repeat measurements of serum digoxin with all-cause mortality using data from the Effect of Digoxin on Mortality and Morbidity in Patients with Heart Failure (DIG) trial.

Results

Overall, 28 studies were included that used joint models, 25 (89%) used data from cohort studies, the remaining 3 (11%) using data from clinical trials. 21 (75%) of the studies used biomarkers and the remaining studies used imaging parameters and functional parameters. The exemplar findings show that a per unit increase of square root serum digoxin is associated with the hazard of all-cause mortality increasing by 1.77 (1.34–2.33) times when adjusting for clinically relevant covariates.

Conclusion

Recently, there has been a rise in publications of joint modelling being applied to heart failure. Where appropriate, joint models should be preferred over traditional models allowing for the inclusion of repeated measures while accounting for the biological nature of biomarkers and measurement error.

Background

Heart failure is a condition where there are well documented inter-relationships between numerous physical, biochemical and imaging characteristics and outcomes. Many studies tend to examine these associations with outcomes using data from one point in time such as randomization in a trial or the start of a cohort study. This fails to account for changes in characteristics over time. Just as baseline values may be associated with outcomes, changes in variables are also associated with changes in outcomes e.g., falling levels of natriuretic peptides are associated with lower mortality. However, analysing the association between changes in variables and outcomes is often performed with traditional time to event models using change values and starting follow up for outcomes once change has occurred. More recently joint modelling, combining two or more statistical models to increase efficiency and reduce bias, has gained favour in the literature as a method of dealing with this issue. The most common type of joint modelling within medicine is the joint modelling of repeat measure longitudinal data (e.g., repeated measures of biomarkers over time) and time-to-event data (i.e., survival data) which are linked through an association structure via shared random effects [1,2,3,4]. This seeks to improve efficiency and reduce bias in respect of treatment effect, censoring and mortality when compared against traditional models. Joint models (JMs) of this type are formed of two sub models: a longitudinal model such as a linear mixed effect (LME) model (which allows the modelling of longitudinal changes in biomarkers or other characteristics like blood pressure) and a survival model such as a Cox Proportional Hazards (Cox PH) model to model the outcome e.g., mortality. The LME model allows for both fixed and random effects accounting for non-independence of repeated measures from the same patient, whilst also allowing for unevenly spaced measurement occasions, biological variances, and measurement error [5]. The survival model allows for covariates and typically includes an association parameter representing the association between the longitudinal and survival process [2, 3, 5,6,7]. JMs which use data from randomized controlled trials (RCTs) can also model the overall treatment effect as well as the treatment effect on both the longitudinal and survival models [3] i.e., the effect on the characteristic and the effect on the outcome. This is analogous to other JMs which are used in the cardiovascular literature, for example in recurrent events analyses where models that examine the effect of a treatment on recurrent hospitalisations while also estimating a treatment effect for a terminal event such as death [8].

Given the increasing use of JMs, the aim of this paper is to review the application of joint modelling in heart failure and to provide guidance on how to assess and interpret results of joint modelling. To achieve this, we conduct a systematic review to identify and critically review current applications of joint modelling within the heart failure population and then present a critical summary of how joint modelling can be applied to heart failure data sets with use of an illustrative example. We examine the association between changes in serum digoxin levels and mortality in the Effect of Digoxin on Mortality and Morbidity in Patients with Heart Failure (DIG) trial as prior studies have tried to examine the association between digoxin levels and outcomes, and suggested that higher levels at one month following randomization may be associated with higher mortality [9].

Methods

Systematic review: joint modelling applications within heart failure

Our systematic review was conducted following the Preferred Reporting Items for Systematic Reviews (PRISMA) framework [10] and the protocol is registered with PROSPERO, registration number: CRD42020210056. The aim of the review was to identify journal articles which employed joint modelling on an adult heart failure population to review how joint modelling was being applied to heart failure.

Searches

Our search strategy is provided in the figure S1, Additional file 1. Medline, Embase, Scopus and Google Scholar were searched, with the last search being conducted on 10^th December 2021.

Screening

Articles were screened by two reviewers and full text was accessed for relevant articles. To capture all available articles no date limit was set and only English language articles were included. Only full text journal articles where joint modelling was applied to an adult heart failure population were considered for inclusion. Data were extracted by two reviewers.

Exemplar: joint modelling of serum digoxin concentration and all-cause mortality

To demonstrate applications of JMs on heart failure data, the ‘The Effect of Digoxin on Mortality and Morbidity in Patients with Heart Failure’ (DIG) [11] trial was used. The dataset was obtained from the Biologic Specimen and Data Repository Information and Coordinating Center (BIOLINCC) under application #9257.

Statistical methods

Only data from patients on the treatment arm with at least one measurement of serum digoxin concentration (SDC) was used. SDC measurements were right skewed and therefore a square root (sqrt) transformation was applied. For this illustrate example, only patients with no missing covariates were included.

The JM Package was used to fit all joint models, this package allows the fitting of joint models of longitudinal and time-to-event data in R under a maximum likelihood approach [12].

Time must be modelled on the same scale for both models, and was modelled in the form of months (28 day calendar month) since randomisation; for SDC, time was taken as the specimen time. While the JM package allows for non-linear effects of time; for simplicity and ease of interpretation only linear terms were included.

The JM package requires an LME as fitted by the ‘LME’ function from the nlme package for the longitudinal sub-model [12]. For this example, both an unadjusted and adjusted LME model were fitted. With all models using sqrt SDC as the response variable. The unadjusted model included random intercepts as random effects. The adjusted model included the main effects of: estimated Glomerular Filtration Rate (eGFR), patient reported self-adherence, hours since last dose of the study drug and dose as fixed effects and included random intercepts and slopes for random effects. Full model equations for the LMEs and all other models are included in Table S1, Additional file 1.

Time-to-Event models for the JM package can be fit using either the ‘coxph’ or ‘survreg’ functions from the Survival package. For simplicity, Cox PH models fitted by the ‘coxph’ function were used. Like the LME both an unadjusted and adjusted model were fitted. The adjusted model containing the covariates of age, sex, ejection fraction, New York Heart Association class, history of hypertension, ischemic etiology of heart failure and body mass index. These covariates were selected on the basis of clinical relevance and prior knowledge of factors associated with outcomes in heart failure. The outcome examined was all-cause mortality.

Three JMs were constructed from both the unadjusted and adjusted LMEs and Cox PH Models as previously defined. Table 1 summarises the formulation of the JMs.

Table 1 Formulation of JMs Included in Exemplar

As an additional analysis the JMs were compared against traditional models. The traditional models were Cox PH models using first and last measurements of sqrt SDC as a covariate and an extended Cox PH model including sqrt SDC as a time varying covariate. All models were adjusted for the same clinical covariates as the adjusted time-to-event models from the JMs. The model fit of the JMs were compared against each other using the Akaike Information Criterion (AIC). Likewise, the model fit of the traditional models were compared against each other using AIC. All models were compared for performance using a discrimination index: c-index for the traditional models and a dynamic discrimination index for the JMs. The dynamic discrimination index was obtained using the function ‘dynCJM’ from the JM package. Based on the time dependant discrimination index proposed by Antolini et al. the dynamic discrimination index in this context provides a single statistic to summarise the discrimination power of the model over the follow-up time and is calculated from a weighted average of time-dependant AUCs which is comparable to the well-known c-index. Like the c-index it does not take into account censoring [13, 14]. The parameter estimates and standard errors from the model were also compared. One hundred bootstrap samples were used to internally validate the comparison of the discrimination index.

For descriptive purposes, categorical variables are represented as percentages, continuous variables are represented as median (IQR). JM association parameters are represented as hazard ratios (HRs) and 95% confidence intervals (CI). A time dependent association parameter is the hazard of all-cause mortality per one unit increase of sqrt SDC at any time point. A time dependent slope association parameter is the hazard of all-cause mortality per one standard deviation increase in the slope of sqrt SDC at a time point (known as the instantaneous or current slope). A p-value of less than 0.05 is considered statistically significant.

All statistical analysis was conducted using R Version 4.0 [15] and JM package version 1.4–8 [12].

Ethical approval was not required for this systematic review and exemplar. All methods were carried out in accordance with relevant guidelines and regulations.

Results

Systematic review: use of joint modelling in heart failure

Figure 1 shows the PRISMA flow diagram, with 28 studies meeting the criteria for inclusion. Table S2, Additional file 1 outlines the data sources of the 28 studies which met the inclusion criteria, the earliest included study being published in 2014 and between 4–7 studies being published per year from 2017 to the last search (10^th December 2021).

Fig. 1

Prisma 2020 flow chart

From the included studies, 25 (89%) used data from cohort studies and the remaining 3 (11%) studies used data from clinical trials. It is worth noting that 10 (36%) of the cohort studies used data from the Bio SHiFT study [16], likely because of the study design with its focus on repeated measurements of biomarkers.

From the studies, 22 (78%) exclusively included patients with heart failure, 6 of which specified patients with heart failure with reduced ejection fraction (HFrEF). The remaining 6 (22%) studies exclusively included patients with acutely decompensated heart failure. There were also studies which further selected patients on specific characteristics such as patients implanted with cardiac devices such as Cardioverter-Defibrillators (ICD) or Cardiac Resynchronization Therapy Devices (CRT-D), patients with advanced heart failure and patients who had undergone transcatheter mitral-valve repair.

Rationale

The most common rationale for using joint modelling was to assess the association of a biomarker with the hazard of an event. Other rationale included: joint modelling as a sensitivity analysis, reduction of bias due to censoring / mortality, comparison of prognostic models e.g., Weibull survival models and JMs, personalised prognostication using JMs, accounting for measurement error, different follow-up times and efficiency through combining data (i.e., smaller standard errors [5, 17]).

Longitudinal data

Table 2 summarises the longitudinal data used in the included studies; 21 (75%) studies used biomarkers with the most common biomarker being N-Terminal Pro-Brain Natriuretic Peptide (NT-proBNP). Some studies included multiple biomarkers in longitudinal sub-models, and some used multiple JMs of different biomarkers. The remaining studies used imaging parameters such as Left Ventricular Ejection Fraction (LVEF) and functional parameters such as health status, physical activity and depression for their longitudinal data. All but two studies specified their longitudinal sub-models as a linear mixed effects model.

Table 2 Summary of Longitudinal Data of Included Studies

Time-to-Event (Survival data)

Many studies included multiple events for their survival data through use of a composite outcome or multiple JMs. Table 3 shows that composite outcomes were the most common, but the events of composite outcomes varied by patient population as shown in Table S3, Additional file 1. The second most common event was all-cause mortality. Most models utilised Cox PH models for their survival sub-models with only two studies specifying a parametric Weibull model.

Table 3 Survival End Points of Included Studies

Missing data

Common joint modelling packages such as JM and JMBayes allow for both uneven spacing and missing longitudinal measurements. Both these packages require all covariates from both longitudinal and survival sub-models to be complete. This common limitation resulted in 13 (46%) of the included studies using imputation methods to complete missing data.

JM

All included studies modelled their JMs with R. The two most common packages being JM and JMBayes, with 7 (25%) studies using the JM package and 10 (36%) using JMBayes and another 3 (10%) studies specified both packages. Three studies used custom code, one study used the joineRML package and the remaining 3 (10%) did not specify the package used.

The R Packages used show both use of frequentist and Bayesian analysis. While the frequentist method use the maximum likelihood approach and is more comparable to more traditional models the Bayesian approach typically relies on Markov chain Monte Carlo (MCMC) sampling algorithms and may improve analysis by using related historical information and allowing for more flexible estimation [45].

Presentation of results

Generally, the results from the JMs included the hazard ratio of the associated longitudinal outcome of interest on the time-to-event outcome; this was typically either the association of the value of the longitudinal outcome or the slope of the outcome on the time-to-event model.

The longitudinal sub-model is often presented as a coefficient or a graph of the average change in the longitudinal outcome of interest over time, these graphs are commonly split into groups of subjects e.g., those who did or did not experience the time-to-event outcome of interest. An example of this is illustrated in Fig. 2 from the Vark et al. study [42].

Fig. 2

Average Estimated Biomarker Pattern, Combined with Individual Biomarker Measurements. During Follow-up in Patients with and without the Primary Endpoint from the van Vark et al. Study. Legend: Reprinted with permission from van Vark et al. © 2017 The American College of Cardiology Foundation [42]

Commonly used packages in R provide capability to visually represent a trajectory of the longitudinal measure and the resulting changes in survival probability as shown in Fig. 3 taken from Zhang et al. where the trajectory of the longitudinal measures (NT-proBNP) is plotted on the left and the survival probability with 95% confidence intervals are plotted on the right. This is useful when looking at individual patient trajectories, in their example Zhang et al. show how the probability of survival changes in response to changes of the trajectory of NT-proBNP and a narrowing of the confidence intervals can be observed with the increase of measurements of NT-proBNP [44]. These such plots while useful were not common amongst the included studies.

Fig. 3

Dynamic survival probabilities with 95% CI based on various measurements of NT- ProBNP for a patient whose values fell. Legend: Reprinted with permission from Zhang et al. © 2018 Elsevier B.V. All rights reserved [44]

Given this is a relatively new approach to analysing longitudinal data simultaneously with survival data studies often compare results from JMs against more traditional models such as a Cox PH model with only a singular measurement of the variable of interest.

Joint modelling outcomes

Generally, the JMs of the included studies performed favourably in terms of improving prognostication and identifying associations with adverse events. The Bio-SHiFT study being the most common data source explored a variety of biomarkers and imaging parameters including cell adhesion circulating biomarkers, fibrinolysis factors, renal markers, echocardiographic parameters, Micro Ribonucleic Acid (MiRNA’s), cardiac remodelling biomarkers and macrophage and neutrophil related biomarkers highlighting a variety of biomarkers and imaging parameters that were associated with the adverse events [23,24,25,26, 32, 35,36,37,38,39]. Whilst many of the tested biomarkers and parameters produced positive results, Van den berg et al. suggested that repeated measures of imaging parameters such as LVEF do not add any more value than single parameters due to the lack of change in those parameters over the observation period [39]. While this suggestion may be true for the population of the Bio-SHiFT study, both Alvarez-Alvarez et al. and Van den Berge suggested that for other heart failure populations repeated measurements of echocardiographic parameters such as LVEF can be useful, with Alvarez-Alvarez et al. investigating these parameters in a chronic heart failure population after CRT and Van den Berge exploring these parameters in an acute heart failure population [19, 40].

One of the key biomarkers which was explored in many studies was NT-ProBNP with Zhang et al., Castelvecchio et al. and Van Boven exploring its association with adverse events in a chronic heart failure population [28, 36, 44]. While these studies demonstrated the association between NT-ProBNP and adverse events, the Zhang et al. study suggested that the most recent value of NT-ProBNP had a similar predictive value as the serial measurements, but similar to the Van den Berg et al. study this may simply be due to the lack of change in values of NT-ProBNP within the study population and may not be generalisable [44].

Other key biomarkers which appeared in multiple studies were High sensitivity Troponin T (HsTnT), C-Reactive Protein (CRP), Cancer Antigen 125 (CA125), creatinine, Suppression of Tumorigenicity 2 (ST2), Galectin-3 (GAL-3) and Growth Differentiation Factor 15 (GDF-15) [22, 27, 33,34,35, 37, 41, 42], indicating that repeat measures of these markers are of interest within heart failure populations. Additionally, other less frequent markers included lactase dehydrogenase trends (LDH) [30], Red blood cell Distribution Width (RDW) [17] and ambulatory markers such as Systolic Blood Pressure (SBP), heart rate and haemoglobin [18, 21, 27].

Along with biomarkers, functional parameters were also of interest, with Kelly et al. taking a novel approach investigating the association of physical activity as reported by implanted devices i.e., by ICD or CRT-D [31] and Arnold et al. using the joint modelling of health status in the form of Kansas City Cardiomyopathy Questionnaire Overall Summary Score (KCCQ-OS) score and all-cause mortality as a sensitivity analysis to illustrate how censoring attenuated health status with respect to treatment effect [20].

Along with the echocardiographic parameters mentioned above imaging parameters were used in a total of four studies, using many of the parameters obtained from imaging such as LVEF, left ventricular end-diastolic diameter (LVED), left ventricular end-systolic diameter (LVES) and Tricuspid Regurgitation [40, 43].

Association between serum digoxin concentration and mortality

Baseline characteristics

Table 4 shows the baseline characteristics of included patients (n = 2012), with a median age of 64 years, 22% of patients being women, median ejection fraction was 29% and 35% of patients died.

Table 4 Baseline Characteristics According to Patient Status (Dead / Censored) and Overall Sample (n = 2012)

JM – longitudinal data sqrt SDC over time

The coefficients from the longitudinal sub-model of JM 2 (Table S4, Additional file 1) for time -0.004 (-0.005—-0.002) suggests that the sqrt root SDC decreases by 0.004 per month after adjusting for covariates. Figure 4 shows a representation of the predicted, and therefore adjusted, average trajectories of SDC over time from JM 2, divided into patients who died during follow-up and those who did not, it suggests that patients who died had on average higher levels of SDC.

Fig. 4

Average Trajectories of SDC by Patient Status as Predicted by JM 2 with Observed Values and Trajectories of SDC. Legend: a Average Trajectories of SDC by Patient Status as Predicted by JM 2 with Observed Values and Trajectories of SDC on original axis scale, b Average Trajectories of SDC by Patient Status as Predicted by JM 2 with Observed Values and Trajectories of SDC on scaled axis scaled for readability. Average Trajectories were predicted using JM2 for patients whose status were either alive or dead based on mean and mode characteristics (covariates) of each stratum of patients (Alive or Dead)

JM – hazard ratios

Table 5 illustrates the results from the survival sub-models of the JMs in terms of a HR (95% CI) and p-value. All the JMs are time dependent relative risk models with a baseline risk function and as such the HRs can be interpreted similar to HRs from a proportional hazards model such as a Cox PH model. Focusing on the time dependant association parameter, the unadjusted model (JM 1) with a HR of 5.32 (3.07 – 9.22) suggesting a fivefold increase in the hazard of all-cause mortality per unit increase of sqrt SDC. This association is attenuated when adjusted for clinical covariates in JM 2 with a HR of 1.77 (1.34—2.33). The time dependent slope parameter of JM 3 is above the significance threshold (p-value 0.092) indicating insufficient evidence to establish an association between the slope of sqrt SDC and all-cause mortality. Neither the HR for the time dependant parameter or the HR for the time dependant slope parameter of JM are above the significance threshold (p-values of 0.427 and 0.13, respectively) suggesting insufficient evidence to establish an association with either value or slope when the model is adjusted for both. The interpretation of the time dependant parameter of this model would be the hazard of all-cause mortality per unit increase of SDC for patients having the same slope. The interpretation of the time dependant slope parameter of this model would be the hazard of all-cause mortality per one SD increase in slope for patients having the same level of sqrt SDC.

Table 5 Event Summary of JMs Represented Hazard Ratios

JM – individual patient trajectories

Figure 5 shows the individual patient trajectories of a patient randomly selected (from patients with at least four measurements of sqrt SDC) at four different time points. These plots contain the longitudinal measurements of sqrt SDC as fitted by the JM on the left and the survival probability on the right, the dashed line indicating the last point the patient was known to be alive and the start of the survival curve, this point changes with each added measurement and as a result the survival curves are not directly comparable. However, these plots demonstrate how the measurements of sqrt SDC effects the survival probability and how the confidence intervals change over time with more measurements.

Fig. 5

Individual Patient Trajectory of sqrt SDC and Survival Probability from Randomly Sampled Patient as Predicted by JM 2

Additional (Comparative) analysis

Results from the comparative analysis suggest that JM 2 had the overall best performance of the JMs with the lowest AIC, highest log likelihood and joint highest discrimination index of all models. The lower HR from the extended Cox PH model suggests that it underestimated the HR of the sqrt SDC parameter likely due to the nature of SDC as a time varying covariate; SDC being an endogenous biological covariate subject to measurement error, biological variances, being able to change between measurements and finally requiring the subject to be alive at measurement. Underestimation of the association parameter has been previously demonstrated in simulation studies [46].

Internal validation of the discrimination index using 100 bootstrap samples showed that JM2 with and mean discrimination index of 0.66 (range 0.6—0.72) outperformed the extended Cox PH model with a mean discrimination index of 0.65 (range 0.62 – 0.67) 71% of the time with respect to the discrimination index and 66% of the time when compared to Cox PH last measurement model with a mean discrimination index of 0.65 (range 0.62—0.68). The full additional analysis is available in Additional file 1.

Discussion

In 2016, a systematic review by Sudell et al. showed an increase in use of JMs of longitudinal and time-to-event data over time. However, only 3 identified studies used ‘heart related’ data; the most common applications were to cancer and HIV/AIDS studies [47]. Developing on their search strategy.

We identified 28 studies by systematic review applying joint modelling within an adult heart failure population, and with use of an illustrative example have shown how to fit and interpret a JM. We have also shown how a JM approach can be used to examine the association between a biochemical test and outcomes in patients with heart failure.

Open-source software packages available in R such as JM and JMBayes make joint modelling more accessible reflected by the 20 (71%) studies using these packages. While these packages limit the JMs by way of underlying methodology [12, 48], if this methodology is not suitable custom code may be written as illustrated by the Hurst et al. study [30]. Both the JM and JMBayes packages also contain limitations around missing data in covariates, while the packages allow for missing longitudinal data, they do not allow for missing covariates in the sub models used to build the JMs, this results in the need to either use a complete case with regards to the covariates as shown in our exemplar or use imputation techniques such as multiple imputation as highlighted in the included studies.

Due the clinical nature of the included studies we found that studies often lacked details on the formulation of the JMs, e.g., the baseline risk function. Whilst this information could usually be derived by considering the packages used to fit the JMs, this information may be useful for reproducibility. We also identified that there was a lot of ambiguity around the origin of figures; whether or not they came from JMs or the individual components of the JMs e.g., a linear mixed effects model, modelled independently of the JM. We therefore suggest the need for clarity and transparency of the presentation of results from JMs.

It is also important that the results are easily understandable to a general audience. For example, the HR of a time dependent association parameter is intuitive but the HR of a time dependant slope parameter less so. Clinicians will often consider trends of biomarkers in day-to-day decision making so understanding these association parameters are key to relating them to clinical practice.

One driving motivation of the use of JMs was utilisation of repeated measures to inform prognosis and the comparison against a single measure. Most studies investigated biomarkers such as NT-ProBNP, CA125 and renal markers. However, JMs were not only limited to biomarkers; with such studies as van den Berg et al. investigating echocardiographic parameters [39], Arnold et al. focusing on health status [20] and Kelly et al. exploring physical activity as reported by an implanted device [31]. The use of these data illustrates the robustness of joint modelling. Another key rationale was the use of joint modelling to reduce bias due to censoring and mortality. Bias of this nature often occurs because subjects who are sicker are more likely to experience the event of interest or withdraw from the study earlier than those who are healthier leading to fewer longitudinal measurements [2]. To overcome this joint modelling provides a framework that acknowledges the underlying relationship between the longitudinal and event process through the use of a joint distribution [5]. The Arnold et al. study illustrates this bias visually highlighting how censoring likely attenuated heath status with respect to treatment effect [20]. Further, we have highlighted the use of joint models to handle missing not at random data through use of a joint distribution [5].

Only three studies used data from RCTs [20, 22, 27], as previously mentioned joint modelling can be used to reduce bias with respect to treatment effect. Whilst this highlights a potential gap in the literature it should be noted that during screening, we identified numerous studies using joint modelling as a sensitivity analysis with results consistent to those from separate longitudinal and survival models but were excluded from review as not enough details about the models were included for full appraisal.

Compared to cancer studies, there was a lack of focus on quality-of-life data with only one study including quality of life in the form of a Kansas City Cardiomyopathy Questionnaire and SF-36 scores [20], and one which included depression by means of patient health questionnaire 9 scores [29]. Whereas joint modelling with quality of life is much more prevalent in cancer studies [3, 20]. This highlights another area which may be of interest to future studies using joint modelling in heart failure.

The Bio-SHiFT Cohort made up 36% of studies primary data, illustrating how a study can be developed to fully use the capabilities of joint modelling; with frequent blood sampling and measurements of endpoints the study leads the way for further larger studies of this nature [16].

From the studies using data from the Bio-SHiFT Cohort we identified 3 studies which only selected baseline and the last two measurements closest to the endpoint [23, 38, 49]. While justified to investigate the trajectories before and after an event it should be noted that this kind of analysis could lead to bias and should only be conducted with proper justification.

Many of the included studies demonstrated how repeated measures added value with respect to both prognostication and model fit. The outcomes of the JMs illustrate how joint modelling can improve on traditional models and highlights the use of joint modelling to assess associations of various biomarkers, imaging parameters and functional parameters, and adverse outcomes as well as provide dynamic predictions. However, there were studies which stated repeat measurements did not add prognostic value or improve model fit. For example, Van den Berg et al. stated that repeated measurements of echocardiographic parameters were associated with adverse events but did not add prognostic value due to the lack of change in measurements over time [39]. Whist this may be true for the Bio-SHiFT cohort, both van den Berge et al. and Alvarez-Alvarez et al. illustrated that given the right context these repeat measurements can still add value [19, 40]. This highlights an important caveat regarding JMs, in that the cost may not outweigh the benefit of the JMs; whilst biomarkers are routinely collected at little added cost other parameters may be costly to collect and an understanding of the temporal patterns of these parameters prior to joint modelling is advisable.

Our exemplar shows how joint modelling can be applied to older studies in order to maximise information from data that was sometimes collected but unused in expensive clinical trials. We highlighted how they compare to traditional models and how they can compete and improve on these models while also providing new clinical insight. The HR of the extended Cox PH model when compared to JM2 and the last measurement Cox PH Model suggests that it underestimated the association parameter, as previously stated likely due to the nature of SDC as a covariate. Whilst the extended Cox PH allows for repeated measures it does not account for measurement error, biological variance or that SDC may vary between time points or after the last observed measurement; this underestimation has previously been demonstrated in in simulation studies [46]. Joint modelling while allowing for repeated measurements of SDC can handle the biological endogenous nature of SDC providing better inferences [5]. Our results suggest that higher values of SDC rather than the slope is associated with higher mortality in patients with heart failure. Our work extends the findings of prior studies that have tried to examine this association with landmark methods which do not perform as well as shown by our exemplar analysis. The implications of this finding are that for patients, their SDC level should be kept as low as possible while still maintaining adequate dosing and in patients with high SDC consideration may have to be given to reducing the dose. There is however an issue of reverse causality, sicker patients may be prescribed higher doses and consequently have higher SDC. However, higher SDC could still act as an indicator of risk and should alert clinicians to reassess the patient and consider other therapies for their heart failure.

Our exemplar only included patients who had no missing covariate values. While this is satisfactory for an exemplar, it can lead to loss of information and possible bias in research studies and the best practice may be to use multiple imputation [50]. However, it should be noted that multiple imputation requires pooling for valid inference, which may cause issues with computational complexity and the need for pooling for dynamic predictions.

Our internal validation using 100 bootstrap samples showed JM2 outperformed the extended Cox PH and Cox PH last measurement models with respect to the discrimination index most of the time (71% and 65% respectively) within the bootstrap samples providing validation to the prognostic performance of the JM. However, the range of the discrimination indices of JM2 is wider than both other models suggesting more variability of discrimination with the joint model within the bootstrap samples. While our exemplar used a dynamic discrimination index for prognostic comparison against traditional models, we found that there was little consistency in methods used to compare JMs against each other and traditional models, highlighting a need for consistency when evaluating JMs. We would suggest that any model specifications or parameters are clearly described to allow any comparisons to be made in future research.

JM 3 in our exemplar included a slope parameter corresponding to the rate of change in sqrt SDC at a time point, known as the instantaneous or current slope. As previously stated, this parameter can be difficult to interpret and as such the JMBayes2 package offers the use of other slope parameters such as delta change i.e., change in the last month / year prior to the time point. This parameter should be easier to interpret and is likely to be more prognostic than currently used slope parameters [51].

Both our review and exemplar highlight the various output and figures that can be produced from a JM and show how powerful joint modelling can be, with applications for prognostication, research of the association of repeat measurements of biomarkers and an endpoint, sensitivity analysis and more.

While joint modelling has a variety of uses, it may be most beneficial in the presence of informative censoring or dropout, when incorporating time varying exogenous covariates such as biomarkers into survival models, and for prognostic modelling where dynamic predictions are useful. However joint modelling can be computationally complex and take longer to fit than traditional models. It should also be noted that inferences may only be valid where the joint model has been correctly specified both with respect to the sub models and baseline hazard function. Joint models are also only valid when conducted on the same population, this is to say that both the longitudinal and time-to-event responses need to come from the same group of subjects. Joint models may not provide better prognostic inference where there is limited variability in the repeated measures such as shown by van den Berg et al. [39].

Our exemplar has some limitations such as the limited number of repeat measurements. This may have affected the power to estimate the slope association parameter and overall accuracy of the model.

Conclusions

In conclusion, this hybrid systematic review with exemplar highlights the rise in the use of JMs within heart failure, and our exemplar illustrates how JMs can be fitted to existing datasets adding value by utilising information from the repeated measures collected. This highlights why JMs are an increasingly popular alternative to traditional models such as Cox PH and Extended Cox PH.