An introduction to joint models-applications in nephrology.

In nephrology, a great deal of information is measured repeatedly in patients over time, often alongside data on events of clinical interest. In this introductory article we discuss how these two types of data can be simultaneously analysed using the joint model (JM) framework, illustrated by clinical examples from nephrology. As classical survival analysis and linear mixed models form the two main components of the JM framework, we will also briefly revisit these techniques.

INTRODUCTION

In nephrology, a great deal of information, ranging from lab results to blood pressure, is measured repeatedly in patients over time during observational studies, clinical trials or simply as patients undergo routine monitoring. Often, data on events of clinical interest, such as dialysis initiation or mortality are collected alongside this longitudinal data. The joint model (JM) framework was developed for application in two main scenarios in which both types of data are available. The first scenario, commonly found in nephrology research, involves investigating whether a longitudinal variable is related to an event of clinical interest. As we will reveal below, traditional models analysing time-to-event data and longitudinal data separately may be inadequate for answering such questions. In those cases, modelling both types of data simultaneously is desirable and allows for the characterization of their relationship. The second scenario in which JMs are often applied concerns research questions where the primary interest lies in establishing the trajectory of a longitudinal variable measured over time. In these longitudinal studies, missing values may be introduced when patients drop out of the study, potentially leading to biased estimates. In such cases, JMs can be applied to correct the longitudinal trajectory for dropout. Below we present examples describing the use of JMs in both scenarios.

JMs have been extensively studied in the technical statistical literature and are still a field of ongoing research. In the past, JMs were seldom utilized in clinical research due to their complexity and the computational power required to run this type of model. Over recent years, however, due to advancements in computational power and the development of user-friendly software, JMs have gained in popularity, leading to their application in various fields of clinical research [1-3]. Nonetheless, the technique has not yet been fully embraced by the nephrology community, and although extensive statistical theory exists, little has been published on the application and interpretation of JMs in the nephrology domain. For these reasons, the aim of this article is to provide an introduction to joint modelling, illustrated by clinical examples from nephrology. As classical survival analysis and linear mixed models form the two main components of the JM, we will also briefly revisit these techniques.

The JM

The JM generally consists of two components or ‘submodels’: one to model the survival outcome and the other to model the longitudinal outcome. The JM estimates the two submodels jointly. This allows for characterization of the relationship between both types of data and brings the uncertainty in estimating each submodel together, thus ensuring such uncertainty is properly accounted for in the analysis. Below we describe the most commonly used models for analysing longitudinal and survival data and how these models are combined in the JM. We also explain why using conventional survival models may be unsuitable for the analysis of longitudinal data and argue why JMs should be used instead.

The Cox proportional hazards model

The Cox proportional hazards model is the most popular regression technique for survival analysis. This type of model is used to estimate the effect of a given variable on the risk, or ‘hazard’, of encountering an event of clinical interest within a defined time frame. The term ‘survival’ in this sense does not necessarily entail mortality, but may form any outcome of interest, such as the time to dialysis initiation, hospitalization or a decrease in estimated glomerular filtration rate (eGFR) by 50%. Cox regression can be applied to study the effect of two main types of variables on survival, specifically, time-invariant variables and time-varying variables. As these names imply, the value of a time-invariant variable does not change over time (e.g. patient sex), whereas the value of time-varying variables may change over time (e.g. eGFR). The effect of a longitudinally measured variable on a survival outcome may therefore be investigated by including it as a time-varying variable in the survival model [4]. Although this approach allows for estimation of the effect of the longitudinal variable on the survival outcome, it is not the optimal choice for several reasons [5]. First, the time-varying Cox model assumes that the longitudinal measurements are measured without error, accepting only ‘exogenous’ variables. An exogenous, or external, variable implies that its value is not affected by other variables in the model and that its future value is not affected by the survival outcome (e.g. patient age, outdoor temperature, environmental factors or calendar date). Most biological variables studied in the field of nephrology, such as eGFR, blood pressure and biomarkers, are ‘endogenous’ or internal by nature. The value of an endogenous variable is generated by the patient, measured with error (due to biological variation), and is only observable when the subject is alive. Second, it is not typically possible in nephrology to measure the longitudinal variable at all event times (e.g. at death or initiation of dialysis), as patients are often measured at intervals. Consequently, we do not know the exact underlying value of the variable between measurements or at the time of the event. Several methods are available to interpolate the underlying value. A common solution is to carry the last known value forward, creating a step-like function of the trajectory of the longitudinal variable, which may not accurately reflect the true trajectory followed over time. Such a method may introduce bias, especially when visit dates are infrequent or irregular [6]. JMs were developed to overcome the above-mentioned limitations when studying the effect of a longitudinal variable on survival, and provided the JM is correctly specified, have been demonstrated to increase efficiency and reduce bias in comparison with the Cox model [7]. As a side note, although Cox regression is the most popular survival submodel, it deserves mention that parametric (and flexible parametric) models may also be used in the joint modelling framework.

The linear mixed-effects model

When individuals are followed over time, certain variables may change as they are measured repeatedly in the same individual. In some situations, these time-varying variables are measured at different times for each individual. As the value of these repeated measurements are usually related with each other within a patient, this should be accounted for. Linear mixed models are often applied to take this correlation into account [8]. The term ‘mixed’ in linear mixed models refers to the inclusion of both ‘fixed’ and ‘random’ effects (variables) in the model. A ‘fixed effect’ is a covariate that is assumed to have a constant mean effect across all individuals in a population, whereas a ‘random effect’ covariate allows the effect to vary in each individual. These random effects are used to model the correlated measurements within each individual. A linear mixed model may be seen as an extension of simple linear regression, consisting of an intercept and slope to model the effect of a given variable on the outcome. The linear mixed model adds random effects to this equation, which allows each individual to have his or her own intercept and/or slope. For example, if we were to study the trajectory of eGFR over time, a simple linear regression model would provide a single intercept and a slope for time, describing the mean linear trajectory of eGFR in a given population, whereas a linear mixed model would provide the eGFR trajectory over time for each individual through the random effects, as well as the population average eGFR trajectory through the fixed effect. These models rely on the assumption that the relationship, in this example between eGFR and time, forms a straight line (i.e. the linearity assumption). If this is not the case, then various data transformation or the inclusion of splines (described below) must be explored to model the trajectory appropriately. Linear mixed models can be used to model a continuous outcome and the analogous generalized linear mixed model can be used for a dichotomous outcome [9].

The JM

The JM generally consists of the above submodels, bringing both the survival submodel (modelling the survival outcome) and the linear mixed-effects submodel (modelling the longitudinal outcome) together, allowing us to characterize the relationship between both outcomes. Using the most popular approach, this is achieved through the so-called shared random effects [10]. Simply put, the random effects (i.e. individual trajectories) from the linear mixed model are included (i.e. shared) in the survival model, thus capturing the relationship between the longitudinal outcome and the survival outcome on an individual level. It is important to note that this relationship works in both directions; the survival outcome can be used to inform the longitudinal trajectory on dropout due to the event and, vice versa, the longitudinal variable can be linked to survival, allowing one to associate the longitudinal variable with the event of interest. As we demonstrate in the examples below, joint modelling can therefore benefit the analyses of both longitudinal and survival outcomes, providing insights into both.

APPLYING JMS TO ASSOCIATE A LONGITUDINAL VARIABLE WITH THE RISK OF AN EVENT

In our first scenario, we describe how JMs can be applied when researchers are interested in associating a longitudinal variable with the risk of an event. As described above, simply including a longitudinal variable as a time-varying covariate in a Cox model is subject to limitations and may warrant the use of JMs. In this scenario, the outcome of the longitudinal submodel is incorporated as a covariate in the survival submodel, allowing for the estimation of its effect on the event of interest. Examples of potential research questions include determining the association between blood pressure trajectory and the risk of mortality or the association between eGFR measured over time and the risk of dialysis initiation. JMs may be used to answer aetiological research questions (i.e. to determine the association between an exposure and an outcome), but they are also valuable in the prediction setting, as we demonstrate in the example below. Lastly, JMs can be applied in the mediation analysis framework to estimate the indirect effect of a longitudinal mediator on a time-to-event outcome [11, 12].

JMs and prediction

In nephrology, prediction models are considered essential for facilitating decision-making and informing patients about their prognosis. For practical reasons, prediction models typically use measurements taken at baseline to predict the outcome of interest, thus utilizing only a small fraction of the information available. As described above, JMs allow for the inclusion of all available longitudinal measurements in the form of a linear mixed submodel. Compared with models using a single baseline measurement, the inclusion of this additional information may better reflect disease progression and improve prediction accuracy for two reasons. First, predictions using longitudinally updated values are calculated in closer proximity in time to the date of the event. Second, as the JM uses multiple measurements, estimates will be less prone to error compared with a single measurement, especially if the longitudinal variable is noisy [13].

An additional property of JMs is their ability to associate a longitudinal variable with the risk of event on an individual level [14], allowing for personalized prediction. JMs use the longitudinal information captured at each follow-up visit to predict the outcome for each individual [15]. This is achieved using random effects, which are unique for each individual. Furthermore, it may be of interest to update the predicted prognosis of an individual as new measurements become available. By adding this new information to previous measurements, JMs are capable of updating the risk of an event occurring in the future, thus allowing for ‘dynamic’ and individualized predictions of survival during follow-up [16].

When modelling a longitudinal variable over time, its trajectory may not necessarily be linear. For instance, it has been demonstrated that eGFR trajectories often do not follow a linear path [17, 18]. The so-called splines may be included in linear mixed models to allow for such non-linear trajectories and are therefore ideal for modelling volatile trajectories for variables that do not follow a linear course over time. Splines may be included in both the fixed and random effects of the mixed model and are readily available in most software packages. In the setting of personalized prediction, including splines in the random effects provides a great deal of flexibility when modelling individual trajectories and may result in the improvement of prediction accuracy. Besides splines (fractional), polynomials can also be used to deal with non-linearity in the longitudinal model [19].

In prediction studies using JMs, several methods, or ‘parameterizations’, exist to associate the longitudinal outcome with the survival outcome in the JM. The simplest method is to associate the expected longitudinal value at a given point in time with the event of interest. However, other methods are also available, which may improve prediction accuracy compared with using only the current value. For instance, the slope of the trajectory, reflecting speed and direction, may also be used to associate the longitudinal variable with the outcome. Other parameterizations include using the area under the trajectory as a summary measure for the cumulative exposure up to a given point in time or even a combination of the above. Depending on the relationship between the longitudinal variable and the event of interest, these parameterizations may add predictive value to the model [20].

Example 1: predicting mortality using troponin T

As we outline above, JMs provide an interesting tool in the development of prediction models. JMs enable individualized and dynamic prediction with updated survival probabilities at each new patient visit, allowing for dynamic changes in treatment decisions as patients progress. Moreover, the addition of splines to the random effects offers flexibility when modelling non-linear predictors over time. In the following example we will demonstrate how to use the JM to dynamically predict survival using troponin T (TnT) measured over time in a cohort of 174 (936 measurements) CKD Stages 4 and 5 patients from the European QUALity Study on treatment in advanced chronic kidney disease (EQUAL) study. As a caveat, in this example we greatly simplify the development of prediction models, as methods regarding the choice of predictors (external), validation and measures of discrimination fall out of scope.

The linear mixed submodel

We first use a linear mixed model to establish the trajectory of TnT over time for each individual as well as for the population as a whole. To fulfil the linearity assumption, we log-transform TnT before entering this variable as the outcome in our model. We then add the time variable as a fixed effect to model the population average trajectory of TnT. We also include a random intercept and random slope for time in order to provide each individual with his/her own TnT trajectory. Other covariates may also be added to aid the prediction of the TnT trajectory, but for the sake of simplicity, we will disregard this for now. As TnT often follows a non-linear evolution, we introduce splines to the fixed and random effects for time. Figure 1 illustrates the final mixed model describing the population average and individual trajectories for TnT.

FIGURE 1

This figure depicts each TnT measurement (black dots) and the population (black line) and patient trajectories (red lines) of TnT over time. Non-linear individual trajectories are visible as a result of the inclusion of splines in the random slope for time.

The Cox proportional hazards submodel

As our primary interest is the prediction of patient survival, mortality is modelled as the time-to-event outcome in the Cox model. Patients were censored at 5 years or when they were lost to follow-up. As we wish to achieve the highest possible predictive accuracy for mortality, we also add other established predictors for mortality in our population. In this example, we add the predictors patient age, sex, having had a previous myocardial infarction, pre-existing diabetes mellitus and baseline eGFR. Note that these predictors are time invariant; their values do not change over time. If we wish to add additional endogenous time-varying predictors, then one should do so by adding additional outcomes to our linear mixed submodel (i.e. multivariate linear mixed models).

The JM combines both the linear mixed model and the Cox model, allowing for estimation of the effect of TnT on mortality. Although TnT measured over time was the outcome in the linear mixed model, it now enters the JM as a predictor for mortality alongside the other predictors from the Cox submodel. In this example we will limit the model by linking the current value of TnT to mortality, but as mentioned previously, other parameterizations, such as the slope or the history of the TnT trajectory, may improve prediction compared with the TnT value alone. The standardized coefficients of the JM are presented in Table 1, suggesting a strong association between the level of TnT and mortality. Keeping in mind that we log-transformed TnT and that the coefficients in the survival submodel represent the log-hazard, we can calculate that a doubling in TnT levels results in a 2.6-fold increased risk of death (21.36 = 2.6).

Table 1.

The coefficients from the JM are presented below. TnT, which was the first outcome in the linear mixed submodel, now enters the survival submodel as a covariate (given in bold type), allowing for the estimation of the effect of TnT on mortality

JM	Standardized coefficients (95% CI)	P-value
Linear mixed submodel
Intercept	3.66 (3.55–3.77)	<0.0001
Time (per month)	0.02 (0.01–0.02)	<0.0001
Cox regression submodel
Patient age (per SD)	0.33 (0.06–0.6)	0.02
Sex (male versus female)	0.87 (0.30–1.43)	0.003
Diabetes mellitus (yes versus no)	0,23 (−0,32–0,76)	0.42
Myocardial infarction (yes versus no)	−0.92 (−1.57 to −0.22)	0.02
Baseline eGFR (per SD)	0.23 (−0.08–0.57)	0.15
Log TnT (per SD)	1.36 (1.03–1.71)	<0.001

Keeping in mind that we log-transformed TnT and that the coefficients in the survival submodel represent the log hazard, we can calculate that a doubling in TnT levels results in an 2.6-fold increased risk of death (21.36 = 2.6). In comparison, men have a 2.4-fold increased risk of death in this model (e0.87).

Dynamic prediction

In Figure 2, we illustrate how JMs can be applied to dynamically predict the survival probability in an individual using TnT measured over time. The sequence of figures shows how the JM updates the TnT trajectory as new TnT measurements become available, as well as the patient’s survival probability. In these figures, TnT levels tend to rise over time, resulting in a decline in the 12-month survival probability as the patient progresses.

FIGURE 2

Dynamic prediction of mortality using TnT in a single patient. On the left-hand side of the plot, each asterisk represents a measurements of TnT and the line represents the TnT trajectory modelled over time. On the right-hand side of the plot, the JM updates the survival probability as new TnT measurements become available. Here we present the predicted survival (and 95% CIs) 12 months after the baseline measurement and 12 months after the third and fifth TnT measurement. The survival probability in this patient declines visibly as TnT levels increase over time.

APPLYING JMS TO DEAL WITH INFORMATIVE CENSORING

JMs are often applied in research questions where the primary interest lies in establishing the ‘true’ trajectory of the longitudinal variable. In longitudinal studies where patients are followed over time and repeated measurements are available, missing values may be introduced when patients drop out of the study. In studies where these values are missing (completely) at random, the linear mixed model may be used to obtain unbiased estimates [21]. However, when the rate of dropout is related to the level of the longitudinal variable, dropout is deemed ‘informative’ and is also known as informative censoring [22-24]. This type of missing data is also referred to as ‘not missing at random’ or non-ignorable missingness. For instance, if one wishes to model the eGFR trajectory in a given population, then missing eGFR values will be introduced over time as patients drop out of the cohort due to death. As the level of renal function is related to this event, dropout is more likely to be in patients with a lower eGFR value and dropout is therefore deemed informative. Failure to take this into account may lead to biased estimates, especially when comparing groups with different rates of dropout. In such cases, JMs can address this issue by using the survival component to inform the longitudinal process of missingness. In other words, JMs are able to correct a longitudinal trajectory for informative dropout. We describe such an approach in the example below.

Example 2: eGFR decline in CKD patients with concomitant chronic heart failure

In this example we will use JMs to establish the eGFR trajectory in CKD patients with concomitant chronic heart failure (CHF). As CKD patients with CHF may have a high mortality risk, estimation of the eGFR trajectory may be biased by informative censoring caused by dropout due to death. We will therefore apply JMs to inform the eGFR trajectory on missingness caused by mortality, producing mortality-corrected eGFR trajectories for our population. In this example, we included 290 patients with CHF, with a total of 1237 measurements, derived from the EQUAL study, an observational European cohort of CKD 4 and 5 patients not on dialysis. To illustrate the magnitude of the effect that dropout may have on the trajectory, we will compare the eGFR trajectory obtained from the ‘naïve’ linear mixed model with that obtained from the JM.

We first model the trajectory of eGFR over time in our population using a linear mixed model. An important assumption of the linear mixed model is that the error residuals are normally distributed. To fulfil this assumption, we perform a log-transformation on eGFR before it enters our model as the outcome. Next, we model the population mean trajectory of eGFR over time by adding the time variable as a fixed effect (covariate). We also include a random intercept and a random slope for time to provide each individual with his/her own intercept and slope for their eGFR trajectory. The final linear mixed model is presented in Figure 3 and depicts the population average eGFR trajectory (reflecting the fixed effect for time) and the individual eGFR trajectories (reflecting the random intercepts and slopes for time).

FIGURE 3

In the linear mixed model, the population mean eGFR trajectory (black) reflects the fixed effect for time, whereas the individual eGFR trajectories (red) reflect the random intercept and random slope for time.

As we wish to adjust for dropout caused by mortality, we model the risk of death in the Cox proportional hazards submodel as the time-to-event outcome. Patients were censored when they were lost to follow-up or if they were still alive at 5 years. As we were solely interested in using the distribution of mortality in our population to inform the eGFR trajectory, no other covariates were included in this model. In such cases, the Cox model simplifies to a Kaplan–Meier curve.

The JM combines both the linear mixed submodel and the Cox submodel as described above into a single model. By combining the two, we allow the survival submodel to inform the eGFR trajectory on missingness caused by mortality. In other words, JM corrects the eGFR trajectory for informative dropout due to death. In Figure 4, we compare the ‘naïve’ trajectory from the linear mixed model and the death-corrected eGFR trajectory from the JM, demonstrating how the JM affects the estimation of the eGFR trajectory by correcting for dropout. The JM shows a slower decline in eGFR after correction for mortality compared with the unadjusted naïve mixed model slope. The JM provides a coefficient (−0.49) for the effect of the current value of eGFR at any time point for the risk of death in the survival submodel. Although in this example we chose to associate eGFR with survival using the current value of eGFR at any given time point, other parameterizations, such as eGFR slope, could also have been used. Keeping in mind that we log-transformed eGFR and that the coefficient represents the log-hazard, we can calculate that a doubling in eGFR levels results in a 29% decreased risk of death (hazard ratio 2−0.49 = 0.71). As the JM reflects the eGFR trajectory in the hypothetical situation that none of the patients died, this resulted in a slower mean decline after correction for dropout due to death.

FIGURE 4

This figure shows the population mean trajectory of eGFR over time. The trajectory is less steep after correcting for mortality in the JM (red line) compared with that from the naïve linear mixed model (blue line), reflecting the higher mean eGFR in patients that died earlier in the follow-up.

AVAILABLE SOFTWARE FOR JMS

Multiple software solutions are available for the implementation of JMs. R-packages include joineR [25], joineRML [26] and JM, which takes a frequentist approach to JMs [20], whereas JMBayes can handle multiple longitudinal variables [27, 28]. In Stata (StatCorp, College Station, TX, USA), merlin [29] and the stjm command can be used [30], and for SAS users (SAS Institute, Cary, NC, USA). the JMFit macro is freely available [31].

CONCLUSION

In this introductory article, we describe the conditions in which conventional analysis techniques are inadequate for characterizing the relationship between survival data and longitudinal data, necessitating the use of JMs. In the prediction setting, we demonstrate the flexibility offered by JMs in the modelling of longitudinal predictors and illustrate how researchers can implement JMs for the purposes of dynamic and individualized prognostics. Additionally, we show that JMs are invaluable in studies that focus on establishing the trajectory of a longitudinal outcome where informative censoring poses a risk.

CONFLICT OF INTEREST STATEMENT

None declared.