Joint Models of Longitudinal and Time-to-Event Data with More Than One Event Time Outcome: A Review

1 Introduction

In clinical studies, measurements are often recorded about subjects at each follow-up visit; these response data give rise to longitudinal data. Subsequently, times to one or more clinically significant events are also recorded. The longitudinal data might be censored by one of these clinical events; for example, if the event was death or treatment failure. A growing field of research has emerged that seeks to jointly model these two outcomes —so-called joint modelling. When the outcome processes are correlated, joint modelling has been empirically demonstrated to lead to improved efficiency and reduced bias [1, 2, 3], improved prediction [4], and be applicable to outcome surrogacy [5]. The literature is extensive, with comprehensive reviews given by Hogan and Laird [6], Tsiatis and Davidian [7], Diggle et al. [8], Sousa [9], Proust-Lima et al. [10], and Gould et al. [11].

The classical joint model, from which most research has spawned, involves a single continuous longitudinal outcome and a single right-censored event time. Notwithstanding this simplicity, the joint modelling methodology has been recently extended to generalize both submodels. For the longitudinal submodel, developments include the incorporation of multiple outcomes [12], binary [13], count [14], and ordinal [15] outcomes, and extensions of the classical error and random effects distribution assumptions [16]. For the time-to-event submodel, extensions have involved the modelling of interval- [17] and left-censored [18] data, discrete event times [19], competing risks [20], parametric models [21], spline models [22], and subject- and institutional-level frailty effects [23]. Commensurate with this methodological research, there has been an increase in use of joint models in a wide –range of clinical settings [23, 24, 25, 26] and development of several mainstream statistical software packages [27, 28, 29, 30, 31, 32, 33, 34].

Due to current trends towards personalized medicine, models that utilise all available information more efficiently are of considerable value. In health research, patients may often experience multiple, recurrent, or a succession of clinical events, thus potentially admitting more than one event time. In this article, we comprehensively review the methodological literature for joint models involving multivariate event time data. Although the primary focus is on the ubiquitous shared random effects models, we also describe the growing framework of joint latent class models. Our review encapsulates multiple events, recurrent events (either in the presence of a terminal event, or not) and succession of events data. Although competing risks data can also be considered as multivariate time-to-event data, we do not review these models here as each subject still only admits a single event time. Furthermore, competing risks joint models have been extensively reviewed elsewhere in the joint model literature [35].

2 Longitudinal data submodels

Let Yik(tijk) denote the j-th observed value of the k-th longitudinal outcome for subject i, measured at time tijk, for i=1,…,N,k=1,…,K, and j=1,…,nik. In some cases, only a single longitudinal outcome (i. e. K=1) is considered, which greatly simplifies the model. We will consider both univariate (K=1) and multivariate (K>1) longitudinal data in this review, depending on the methodology presented in each article, but exclusively reserve the subscript k to denote multivariate cases.

In the framework of joint models involving more than one event time, the corresponding longitudinal measurements have predominantly been continuous. However, some models have considered binary and count data (Table 1). As noted earlier, some models have also considered multiple longitudinal outcomes. For a full review of joint models involving multivariate longitudinal outcomes, see Hickey et al. [12]. Król et al. [36] also considered left-censored longitudinal measurements, which is pertinent to biomarker measurements that involved minimum detection thresholds. There are a plethora of modelling approaches for multivariate longitudinal data [37]. In most cases, a generalized linear mixed model (GLMM) [38] is specified. Namely,

where hk(⋅) denotes a known one-to-one link function for the k-th outcome, E is the expectation operator, and μik(⋅) is the linear predictor:

where

and Xik(1)(tijk) and Zik(tijk) are vectors of (possibly time-varying) covariates for subject i associated with fixed and random effects respectively, which can vary by outcome, βk(1) is a vector of fixed effects parameters for the k-th outcome, and bik is a vector of subject-specific random effects for the k-th outcome. We denote the stacked vector of subject-specific random effects for all K outcomes by bi=(bi1⊤, bi2⊤,…, biK⊤)⊤. Some authors have considered including spline terms in Xik(1)(tijk) to capture complex functional forms between the outcome and measurement time [25, 39]. On the other hand, Dantan et al. [40] specified a segmented GLMM with a random change-point, which was intrinsically linked to the time-to-event submodel through one of the transition hazard functions. Random change-points were shown to be particularly useful for capturing changes in the longitudinal trajectory of the outcome following a clinical (pre-)diagnosis.

Generally, for continuous longitudinal outcomes, independent and identically distributed normal errors are assumed. However, extensions to robust skew-normal distributed errors have also been proposed [39]. Subject-specific random effects are generally modelled as being multivariate normally distributed, reducing to a normal distribution in the case of a random-intercepts only model. Different modelling approaches have also been considered. Notably, Huang et al. [41] adopted discrete independent probability distributions. Njagi et al. [14] considered over-dispersed data, and proposed conjugate Beta and Gamma random effects for binary and count outcomes respectively. Several authors who considered multivariate longitudinal outcomes have proposed capturing the cross-sectional association between repeated measures through a correlated errors structure rather than a correlated random effects structure, i. e. Yik(tij)= μik(tij)+εijk, with εij.∼NK(0,Σ) and bik∼Nvk(0,Ψk) [39, 42, 43, 44]. This allows for separate estimation of correlation between repeated measures and between different longitudinal outcomes.

In some cases, a semiparametric paradigm has been adopted. Within the Bayesian framework, Tang et al. [44] and Tang and Tang [39] assumed a Dirichlet process prior for the random effects, removing the need to assume a fixed parametric form, which is therefore robust to potential misspecification. Li et al. [45] suggested a time-dependent vector of random effects, which are independently and identically distributed according to an unknown multivariate distribution. The longitudinal submodels are also specified as marginal proportional rates models; namely, as eq. (1) with hk(.) given by the exponential link function, and linear link functions are also suggested [46]; the time-dependent fixed effects are absorbed into an unspecified smooth baseline function.

Following Henderson et al. [47], additional autocorrelation structure can be incorporated into the model by augmenting eq. (3) to include a zero-mean stationary Gaussian process term. However, such models are met with a substantially increased computational burden, so it is not unexpected that very few methodological articles have considered this extension [48, 49]. Zhang et al. [49], as well as considering correlation for W1i(k)(t) in eq. (3), also allowed for correlation of errors within an outcome over time by letting εik=(εi1k,…,εink)⊤ have zero-mean multivariate normal distribution with a u-lag correlation function given by

A summary of the longitudinal data submodels used in joint models involving multivariate time-to-event data is given in Table 1.

3 Time-to-event data submodels

Let Tig∗ denote the g-th event time for the i-th subject (i=1,…,n). Also, let Ci be a censoring time for the subject such that we actually observe Tig=min(Tig∗,Ci). Typically, continuous event times are observed. Two exceptions were Huang et al. [41], who considered discrete event times, and Rouanet et al. [50] who allowed one of the semi-competing event times to be interval-censored. For each subject i, let the vector Xi(2)(t), which may be time-varying, denote the observed covariate data, and βg(2) denote the coefficient parameters associated with these covariates for the g-th event time. Similarly, for models involving a third submodel (e. g. a joint model of longitudinal data, recurrent events, and a terminal event), we will use the notation Xi(3)(t) and β(3), as appropriate. However, in practice, there will be an overlap between baseline measurements in Xi(1)(t),Xi(2)(t), and Xi(3)(t). Specification of the time-to-event model depends on the type of multivariate event time data and the association structure that gives rise to the joint model. These are described below and succinctly summarized in Table 2 and Table 3. We will denote the association parameters by γg, and any extra random effects terms by θi.

4 Model estimation

Several different estimation approaches have been utilized to fit the models described above (Table 3). Loosely, these methods can be separated as either likelihood maximisation or Bayesian model fitting.

Extending the original joint model developments of Wulfsohn and Tsiatis [67], the expectation-maximization (EM) algorithm has been used in some cases. In the case of Han et al. [51], the latent class membership, longitudinal data submodel random effects, and the time-to-event submodel frailty terms were treated as missing data. In the case of Kim et al. [58], only the random effects were treated as missing data, and recursive formulae used to reduce the number of model parameters required for estimation. Król et al. [36] used penalized maximum likelihood estimation using the Marquardt algorithm, with the penalization performed to obtain smooth estimates of the baseline hazard and intensity functions. Rouanet et al. [50] also utilized the Marquardt algorithm, with the number of latent classes selected according to the Bayesian Information Criterion. Dantan et al. [40] reported using a Newton-Raphson-like algorithm. Huang et al. [41] used automatic differentiation—a numerical technique for simultaneously evaluating a function and its derivatives—with a Newton-Raphson algorithm, which was purportedly faster than the EM algorithm. Njagi et al. [14] and Efendi et al. [54] used a partial marginalisation approach [68] whereby the conjugate random effects are analytically integrated out, and the normal random effects are numerically integrated using standard software. Efendi et al. [54] then exploited the ideas of Heagerty and Zeger [69] to establish marginal effects. Liu et al. [56] and Liu and Huang [57] reported using numerical likelihood maximisation via standard software. Standard errors of all these aforementioned model fits can be estimated from the inverse of the observed information matrix; however, Han et al. [51] reported using the bootstrap method.

Zhang et al. [49] proposed a two-stage estimation strategy. In the first stage, the covariance parameters were estimated from the repeated measures marginal likelihood function, with the mean function estimated by a weighted moving average. In the second stage, the expected likelihood function for the time-to-event data were maximized by an EM algorithm, with Gibbs sampling implemented for the high-dimension numerical integration, and a Newton-Raphson step used for the M-step. Shen et al. [48] developed a two-stage conditional estimating equations approach for model fitting, followed by a bootstrap approach for standard error estimating. As a precursory step, the authors reframed the time-to-event submodel from an intensity function to a conditional rate function. For models that accounted for informative observation times, generalized estimating equations in a semiparametric framework was the standard approach, which yielded consistent estimators [45, 46, 61]. In these cases, theoretical results have been derived on the asymptotic normality, which is subsequently used to make inference on the estimated parameters.

Bayesian estimation of standard univariate joint models has seen increased attention over recent years [28, 30], especially as it is a natural tool for dynamic prediction and model averaging [4]. Moreover, there are multiple disadvantages to the ubiquitous frequentist estimation approach, including but not limited to, computational challenges—something one would expect to be particularly burdensome in a multivariate framework, the dependence on asymptotic approximations, and the complexity of model assessment and comparison. In joint models involving multivariate longitudinal data, Liu and Li [70] compared the performance of Bayesian approaches to maximum likelihood approaches under different strengths of association, and demonstrated superiority of the Bayesian methods with respect to bias, root-mean square error, and coverage. Of the joint models involving multivariate event time data that were estimated using Bayesian statistics [25, 39, 42, 43, 44], Markov chain Monte Carlo (MCMC) methods were employed in all cases with default non-informative prior distributions chosen for the parameters. As noted earlier, Tang et al. [44] and Tang and Tang [39] also assumed a Dirichlet process prior for the random effects, removing the need to assume a fixed parametric form, which is therefore robust to potential misspecification. Tang and Tang [39] explored the sensitivity of results to prior distribution inputs, showing that good prior knowledge led to marginally improved estimation. The Gibbs sampling algorithm was used in all cases, with non-standard conditionals sampled using adaptive rejection or Metropolis-Hasting algorithms. Chi and Ibrahim [42] specifically noted that hierarchical centring [71], as well as some parameter transformations were used to facilitate convergence of the MCMC algorithms. The posterior conditional distributions for each parameter were derived analytically by all authors, except Musoro et al. [25], who exploited the automation provided by the OpenBUGS software. In all cases, assessment of convergence was made using general diagnostic methods; for example, examination of trace plots, autocorrelation plots, and the Gelman-Rubin statistics [72].

5 Software

The ability to fit the models discussed is severely limited by the availability of software packages or modifiable code. Several authors have made code available either in an appendix, an online supplement, or via an online code repository system (Table 3). However, many authors do not report what software was used, or make said code available. Only one article released their code in the form of a software package, namely Król et al. [36], which fits a joint model for a single longitudinal outcome, a recurrent events process, and a single terminal event, and which is available through the trivPenal() function in the R package frailtypack [73].

6 Clinical applications

Development of novel methodology of joint models of longitudinal data and multivariate event times data have predominantly been motivated by real-world clinical datasets. Here, we summarize the applications that have led to the models discussed in this review.

7 Discussion

The case for use of joint models has been made already [1, 74, 75]. Namely, when the longitudinal and event time processes are correlated they reduce the bias obtained from simpler methods, including separate models (e. g. separate LMMs, survival models, recurrent event models, and multistate models), or even the two-stage approach. There has been a myriad of extensions in the joint modelling framework over the past few years, including extensions to multivariate longitudinal data [12] and competing risks data [35]. Relatively fewer developments have been made pertaining joint models involving more than a single event time, which includes multiple events, recurrent events, and a succession of events. Yet, as shown, there are wide-ranging clinical applications for these models. In particular, motivation has stemmed from disease areas representing cancer, infection, cardiovascular disease, neurological disease, mental health, and respiratory disease. Moreover, data were derived from both randomized controlled trials and cohort studies.

The review presented here contributes to this narrowly focused but important topic in joint models by bringing together in a single place and juxtaposing the models and distributional assumptions, outcome types, estimation and software implementations alongside clinical applications. This is a research area of growing interest and clinical importance, and the extensions developed are necessary to appropriately analyze this complex data. However, we found that availability of mainstream statistical software to fit these models is severely limited, and this will ultimately pose problems, since the complexity of the models means that ad hoc programming is required. This is not unexpected as joint models are computationally difficult to fit; a problem that is exacerbated by the extension to joint models involving more than a single event time. In fact, Musoro and colleagues [25] noted that their ambitious attempt to fit a model to 4 longitudinal outcomes and 9 recurrent event outcome types was precluded by computational time; development of approaches that reduce this computational burden are therefore of paramount importance.

The extension of joint models to more than a single event time offers not only improved inference, but also opportunity for dynamic prediction. This has received growing interest in the classical joint model framework [4], but less so in the extension of multivariate event time data. Król et al. [36] developed dynamic prediction and predictive assessment tools for their recurrent events joint model. Others have also discussed prediction in the context of joint models involving multivariate event time data [14, 50, 66]. Dynamic prediction is easily encompassed in a Bayesian joint model framework. Despite this, the use of Bayesian methods for model fitting has been rather limited in the methodological developments of joint models involving multivariate event time data. Moreover, there is also limited research on the role of prior distribution selection. Research to-date has been predominantly technical, and more attention is required on the interpretability of these models in clinical applications. Moreover, the complexity of these models requires further development on diagnostics that will facilitate model selection, including the choice of association structure.