# Analysis of Longitudinal and Survival Data: Joint Modeling, Inference Methods, and Issues

Authors: Lang Wu; Wei Liu; Grace Y. Yi; Yangxin Huang
Journal: Journal of Probability and Statistics
 (2012)

---

## Abstract

In the past two decades, joint models of longitudinal and survival data have received
much attention in the literature. These models are often desirable in the following situations:
(i) survival models with measurement errors or missing data in time-dependent
covariates, (ii) longitudinal models with informative dropouts, and (iii) a survival process
and a longitudinal process are associated via latent variables. In these cases, separate
inferences based on the longitudinal model and the survival model may lead to biased
or inefficient results. In this paper, we provide a brief overview of joint models for
longitudinal and survival data and commonly used methods, including the likelihood
method and two-stage methods.


---

## Body

## 1. Introduction

Longitudinal data and survival data frequently arise together in practice. For example, in many medical studies, we often collect patients’ information (e.g., blood pressures) repeatedly over time and we are also interested in the time to recovery or recurrence of a disease. Longitudinal data and survival data are oftenassociated in some ways. The time to event may be associated with the longitudinal trajectories. Separate analyses of longitudinal data and survival data may lead to inefficient or biased results. Joint models of longitudinal and survival data, on the other hand, incorporate all information simultaneously and provide valid and efficient inferences.Figure1 shows a longitudinal dataset in which CD4 cell counts are measured repeatedly over time in an AIDS study. Here, the time to event could be time to viral rebound, time to dropout, or time to death, depending on the research objectives. Data analysis can mainly focus on either the longitudinal data or the survival data or both. When the analysis focuses on longitudinal data, we often need to address informative dropouts since dropouts are very common in longitudinal studies. When the analysis focuses on survival data, we often need to incorporate time-dependent covariates such as CD4 since the times to event may be associated with the covariate trajectories. Sometimes, the main interest may lie in the association between the longitudinal process and survival process. In any of these cases, joint models are required to feature correlated longitudinal and survival data.CD4 measurements over time. (a) All subjects. (b) Five randomly selected subjects.
(a)(b)Typically, joint models for longitudinal and survival data are required in the following situations:(i)
survival models with measurement errors in time-dependent covariates;(ii)
longitudinal models with informative dropouts;(iii)
longitudinal and survival processes are governed by a common latent process;(iv)
the use of external information for more efficient inference.Joint models of longitudinal and survival data have attracted increasing attention over the last two decades. Tsiatis and Davidian [1] provided a nice overview of early work on joint models, including De Gruttola and Tu [2], Wulfsohn and Tsiatis [3], Henderson et al. [4], and Wang and Taylor [5], among others. More recent work includes Ding and Wang [6], Nathoo and Dean [7], Ye et al. [8], Albert and Shih [9], Jacqmin-Gadda et al. [10], Rizopoulos et al. [11], Wu et al. [12], Huang et al. [13], and Pan and Yi [14], among others. A typical model setting is to assume a mixed-effects model for the longitudinal data and a Cox model or an accelerated failure time (AFT) model for the survival data, with the two models sharing some random effects or variables. The likelihood method is often used, implemented by EM algorithms. Another common approach is based on two-stage methods, which are computationally simpler. Henderson et al. [4] allow different random effects in the longitudinal and survival models, but assume that the random effects are correlated. Bayesian methods have also been proposed [13, 15, 16].Since the literature on joint models is quite extensive, it is difficult to review all references here. In this paper, we provide a brief review of the joint model literature. In Section2, we describe a standard formulation of joint models. In Section 3, we review a commonly used two-stage method. In Section 4, we describe the standard likelihood method. In Section 5, we discuss some extensions of standard joint models. A real data example and a simulation study are presented in Section 6 to illustrate and evaluate the methods. We conclude the paper in Section 7 with discussion.


## 2. A Standard Formulation of a Joint Model

In this section, we consider a standard formulation of a joint model. In the literature, a typical setup is a survival model with measurement errors in time-dependent covariates, in which a linear mixed-effects (LME) model is often used to model time-dependent covariates to address covariate measurement errors and a Cox proportional hazards (PH) model is used for modelling the survival data. We focus on this setup to illustrate the basic ideas.Consider a longitudinal study withn individuals in the sample. The objective is to model the time to an event of interest or survival time. Time-varying covariates are used in the survival model to partially explain the variation in the event times. Let si be the survival time for individual i, i=1,2,…,n. Some individuals may not experience any events by the end of the study so their event times may be right censored. We assume that the censoring is random or noninformative. For individual i, let ci be the censoring time, and let δi=I(si≤ci) be the censoring indicator such that δi=0 if the survival time for individual i is right censored and δi=1 otherwise. The observed survival data are {(ti,δi),i=1,2,…,n}, where ti=min(si,ci).In survival models, some time-dependent covariates may be measured with errors. For simplicity, we consider a single time-dependent covariate. Letzij be the observed covariate value for individual i at time uij, subject to measurement errors, i=1,2,…,n; j=1,2,…,mi. Let the corresponding unobserved true covariate value be zij*. Denote zi=(zi1,…,zimi)T, and zi*=(zi1*,…,zimi*)T. In many cases, the longitudinal covariate measurements are terminated at the event or censoring time ti. For example, this is the case when the events are dropouts. In this case, we have uimi≤ti, and the covariate values after the event time ti are all missing. Let xi be covariates without measurement errors.We consider the following Cox model for the survival data:(2.1)λi(t)=λ0(t)exp⁡(zi*(t)β1+xiTβ2),i=1,…,n,
where λi(t) is the hazard function, λ0(t) is the unspecified baseline hazard function, and β=(β1,β2T)T are unknown parameters. Survival model (2.1) assumes that the hazard function λi(t) depends on the unobserved true covariate values zi*(t), rather than the observed but mismeasured covariate value zi(t).In Cox model (2.1), the time-dependent covariate value zi(t) should be available at any event time t. In practice, however, covariates are usually measured intermittently at times (say) {uij,j=1,2,…,mi} for individual i, with the measurement times possibly varying across individuals, leading to possible “missing covariates.” Moreover, covariates may be measured with errors. Therefore, it is common to have both missing data and measurement for errors in time-dependent covariates, which must be addressed when conducting inference based on the Cox model (2.1). For simplicity, we assume that the missing covariate data are missing at random [17].To address either missing data or measurement error or both, a standard approach is to model the time-dependent covariates. A common choice is the following LME model:(2.2)zi=Uiα+Viai+ϵi≡zi*+ϵi,i=1,…,n,
where Ui and Vi are known design matrices, α is a vector of fixed-effects parameters, ai is a vector of random effects, ϵi=(ϵi1,…,ϵimi)T is a vector of measurement errors, and the unobserved true covariates are zi*=Uiα+Viai. Often, we assume that(2.3)ai∼N(0,A),ϵi∼N(0,σ2I),
and ai and ϵi are independent, where A is a covariance matrix, σ2 is a parameter, and I is an identity matrix. Here, we focus on the case that the observations of the covariate process is truncated by the event time; that is, no covariate data are available after the event occurs (such as death or dropout).Note that the survival model (2.1) and the longitudinal model (2.2) are linked through the shared random effectsai. In some applications, not necessarily in the measurement error context, the shared random effects may be viewed as a latent process that governs both the longitudinal process and the survival process. The shared random effects induce the dependence between the longitudinal and survival processes, and this dependence suggests the need of joint modelling.There are two commonly used approaches for inference of joint models:(i)
two-stage methods,(ii)
likelihood methods.In the following sections, we describe these two approaches in detail. Other approaches for joint modes have also been proposed, such as those based on estimating equations, but we omit them here for space consideration.


## 3. Two-Stage Methods

In the joint modelling literature, various two-stage methods have been proposed. A simple (naive) two-stage method is as follows.Stage 1.
Fit a LME model to the longitudinal covariate data, and estimate the missing or mismeasured covariates based on the fitted model.Stage 2.
Fit the survival modelseparately, with the missing or unobserved true covariate values substituted by their estimates from the first stage as if they were observed values and then proceed with the usual survival analysis.Main advantages of the two-stage methods, including the modified two-stage methods as described below, are the simplicity and that they can be implemented with existing software. The limitation of those methods is that they may lead to biased inference for several reasons. First, in the estimation of the longitudinal covariate model parameters, the truncations resulted from the events are not incorporated. That is, the longitudinal covariate trajectories of subjects who experience an event may be different from those who do not experience that event, so estimation of the parameters associated with the longitudinal covariate model in the first stage, based only on observed covariate data, may be biased. Second, theuncertainty of the estimation in the first stage is not incorporated in the second stage of the survival model estimation. Thus, standard errors of the parameter estimates of the survival model may be underestimated. Third, all information in the longitudinal process and the survival process is not fully combined in each model fitting to produce the most efficient estimates.The bias in the estimation of the longitudinal model parameters caused by ignoring the informative truncations from the events may depend on thestrength of the association between the longitudinal process and the survival process. The bias resulted from ignoring the estimation uncertainty in Stage 1 may depend on the magnitude of measurement errors in covariates. To address these issues, various modified two-stage methods have been proposed, leading to better two-stage methods.Self and Pawitan [18] considered a two-stage method in which the least-square method was used to fit individual longitudinal covariate trajectories; the resulting estimates were used to impute the “true” covariate values in the survival models and inference was then based on the usual partial likelihood. Tsiatis et al. [19] considered an approximation to the hazard function and then using the approximation to construct the partial likelihood. They replaced the true covariate zi*(t) by an empirical Bayes estimate of the conditional expectation E(zi*(t)∣ziH(t),t≤si), where ziH(t)={zi(u),u≤t} is the covariate history up to time t. They obtained the empirical Bayes estimate from a standard fit of the LME model to the covariate data up to time t for all subjects still at risk at time t. Similar two-stage methods were also proposed in Bycott and Taylor [20] and Dafni and Tsiatis [21].More recently, other two-stage methods have been developed in the literature. In the sequel, we review some of these recent methods. Following Prentice [22], we rewrite the survival model (2.1) as(3.1)λi(t;zi(t),xi)=λ0(t)E[exp⁡(zi*(t)β1+xiTβ2)∣zi(t),xi,ti>t],
which involves an intractable conditional expectation. Following Dafni and Tsiatis [21] and Ye et al. [8], we approximate the above conditional expectation by(3.2)E[exp⁡(zi*(t)β1+xiTβ2)∣zi(t),xi,ti>t]≈exp⁡[E(zi*(t)β1+xiTβ2∣zi(t),xi,ti>t)].
A two-stage method may then proceed as follows. In the first step, we estimate the conditional expectation E(zi*(t)β1+xiTβ2∣zi(t),xi,ti>t) by fitting the covariate model (2.2) to the observed longitudinal data and survival data. In the second step, we then substitute the conditional expectation (3.2) in (3.1) by its estimate from the first step and then we proceed with standard inference for the Cox model. Ye et al. [8] proposed two approaches for the first step, called risk set regression calibration (RRC) method and ordinary regression calibration (ORC) method, respectively. The idea is to fit the LME covariate model (2.2) to either the observed covariate data in the risk set or all observed covariate data.Note that the bias resulted from the naive two-stage method is caused by the fact that the covariate trajectory is related to the length of followup. For example, subjects who drop out early or die early may have different trajectories than those who stay in the study. Thus, much of the bias may be removed if we can recapture these missing covariate measurements due to truncation by incorporating the event time information. Albert and Shih [9] proposed to recapture the missing measurements by generating data from the conditional distribution of the covariate given the event time:(3.3)f(zi∣si;θ)=∫f(zi∣ai;θ)f(ai∣si;θ)dai,
where the covariate zi and event time si are assumed to be conditionally independent given the random effects ai, and θ contains all unknown parameters. They approximate the conditional density f(zi∣si;θ) using a LME model, and then use standard software to simulate missing data from f(zi∣si;θ). Once the missing measurements are simulated, the covariate model is then fitted to the “complete data,” which are used in the second step. The procedure is iterated several times to incorporate the missing data uncertainty. Thus, the idea is similar to a multiple imputation method with nonignorable missing data. Such an approach may reduce the bias resulted from truncations.To incorporate the estimation uncertainty in the first step, we may consider aparametric bootstrap method as follows.Step 1.
Generate covariate values based on the assumed covariate model, with the unknown parameters substituted by their estimates.Step 2.
Generate survival times from the fitted survival model.Step 3.
For each generated bootstrap dataset from Steps1 and 2, fit the models using the two-stage method and obtain new parameter estimates.Repeating the procedureB times (say, B=500), we can obtain the standard errors for the fixed parameters from the sample covariance matrix across the B bootstrap datasets. This Bootstrap method may produce more reliable standard errors than the naive two-stage method if the assumed models are correct.Two-stage methods have bearing with theregression calibration method in measurement error literature. Many of these two-stage methods may not completely remove biases. Moreover, they rely on certain assumptions and approximations. The validity of these assumptions and the accuracy of these approximations need to be further investigated.


## 4. Likelihood Methods

The likelihood method is perhaps the most widely used approach in the joint model literature. It provides a unified approach for inference, and it produces valid and the most efficient inference if the assumed models are correct. The likelihood method is based on the likelihood for both longitudinal data and survival data. However, since the likelihood function can be complicated, a main challenge for the likelihood method is computation.


### 4.1. The Likelihood

All the observed data are{(ti,δi,zi,xi),i=1,2,…,n}. Let θ=(β,α,σ,A,λ0) denote the collection of all unknown parameters in the models, where λ0={λ0(ti),i=1,2,…,n}. We assume that the censoring of survival data and the assessment process of longitudinal measurements are noninformative. The (overall) likelihood for all the observed data is given by(4.1)L(θ)=∏i=1n∫f(ti,δi∣zi*,λ0,β)f(zi∣ai,α,σ2)f(ai∣A)dai,
where(4.2)f(ti,δi∣zi*,λ0,β)=[λ0(ti)exp⁡(zi*(ti)β1+xiTβ2)]δi×exp⁡[-∫0tiλ0(x)exp⁡(zi*(ti)β1+xiTβ2)dx],f(zi∣ai,α,σ2)=(2πσ2)-mi/2exp⁡[-(zi-zi*)T(zi-zi*)2σ2],f(ai∣A)=(2π|A|)-1/2exp⁡[-(aiTA-1ai)2].Parameter estimation can then be based on the observed-data likelihoodL(θ) via the maximum likelihood method. Note that the baseline hazard λ0(t) in the Cox model is unspecified. It can be estimated using the nonparametric maximum likelihood method by assuming that λ0(t) takes discrete mass at each failure time ti. Thus, the dimension of the parameter vector λ0 is equal to the number of unique failure times. This converts the semiparametric Cox model to a parametric model, but it introduces a major challenge since standard asymptotic theory for the maximum likelihood estimators (MLEs) may not apply due to the infinitely dimensional nature of λ0.MLEs of the model parameters can either be obtained by a direct maximization of the observed data log likelihood or by using an EM algorithm. Since the observed data log likelihood involves an intractable integral, a direct maximization is often based on numerical integration techniques such as the Gaussian Hermite quadrature or Monte Carlo methods. These methods, however, can be quite computationally intensive if the dimension of the unobservable random effectsai is not low. The EM algorithm is known for its stability and generality, so it is widely used for likelihood inference of joint models [1, 3, 11, 23]. Since the E-step of an EM algorithm still involves an intractable integral, Monte Carlo methods or Laplacian approximations are often used to approximate the conditional expectation in the E-step. In the M-step, the Newton-Raphson method is often used.Hsieh et al. [24] noted that standard errors for estimators of the parameters (α,β,σ) based on the Fisher information matrix may be problematic, because of the semiparametric nature of the joint model. They recommended a bootstrap method to obtain standard errors.


### 4.2. Computational Issues

A main challenge in the likelihood inference for joint models is the computational complexity, since numerical methods or Monte Carlo methods can be very computationally intensive when the dimension of the random effectsai is not small. Moreover, convergence of the EM algorithms can sometimes be an issue. Tsiatis and Davidian [1] and Tsiatis and Davidian [25] proposed an alternative approach, a so-called conditional score method, which makes no distributional assumption for the random effects ai. Their method treats ai as “nuisance parameters” and conditions on an appropriate “sufficient statistic” for ai. Conditioning on this sufficient statistic would remove the dependence of the conditional distribution on the random effects ai. This approach leads to a set of unbiased estimating equations, similar to the usual partial likelihood score equations or generalized estimating equation (GEE). The resulting estimates are, under certain regularity conditions, consistent and asymptotically normal, although they may not be the most efficient. Moreover, their method is relatively simple to implement. Song et al. [26] considered an alternative approach to relax the normality assumption of ai. They assume that the random effects follow distributions in a class of smooth density family, including the standard normal density as a special case. They use the likelihood method for inference via an EM algorithm, but the computation is considerably more intensive than the conditional score method. Song and Wang [27] also proposed a local corrected score estimator and a local conditional score estimator, for which no distributional assumptions are needed for the underlying true covariates.Approximate but computationally more efficient methods for joint models have also appeared in the literature, such as those based on Laplace approximations (e.g., [11, 12, 28]). When the dimension of the integration in the likelihood for joint models is high, the Laplace approximations offer considerably computational advantages over numerical or Monte Carlo methods. Note that the order of the Laplace approximation error is O(mi-1), which cannot be made arbitrarily accurate for a given dataset, where mi is the number of within-individual measurements for individual i. Therefore, the Laplace approximation works well if the number of within-individual measurements is large. Approximate methods based on Taylor series approximations are similar to the Laplace approximation; that is, their performance improves as mi increases.Rizopoulos et al. [11] proposed to use the full exponential Laplace approximation in the E-step of the EM algorithm. Compared to the standard (first-order) Laplace approximation, the full exponential Laplace approximate method has approximation error of order O(mi-2) and requires a much smaller number of within-individual longitudinal measurements to produce reliable results. Lee et al. [28] suggested second-order Laplace approximations. However, these Laplace approximation methods cannot control the magnitude of the approximation errors, unlike Gaussian quadrature or Monte Carlo integration techniques.


## 4.1. The Likelihood

All the observed data are{(ti,δi,zi,xi),i=1,2,…,n}. Let θ=(β,α,σ,A,λ0) denote the collection of all unknown parameters in the models, where λ0={λ0(ti),i=1,2,…,n}. We assume that the censoring of survival data and the assessment process of longitudinal measurements are noninformative. The (overall) likelihood for all the observed data is given by(4.1)L(θ)=∏i=1n∫f(ti,δi∣zi*,λ0,β)f(zi∣ai,α,σ2)f(ai∣A)dai,
where(4.2)f(ti,δi∣zi*,λ0,β)=[λ0(ti)exp⁡(zi*(ti)β1+xiTβ2)]δi×exp⁡[-∫0tiλ0(x)exp⁡(zi*(ti)β1+xiTβ2)dx],f(zi∣ai,α,σ2)=(2πσ2)-mi/2exp⁡[-(zi-zi*)T(zi-zi*)2σ2],f(ai∣A)=(2π|A|)-1/2exp⁡[-(aiTA-1ai)2].Parameter estimation can then be based on the observed-data likelihoodL(θ) via the maximum likelihood method. Note that the baseline hazard λ0(t) in the Cox model is unspecified. It can be estimated using the nonparametric maximum likelihood method by assuming that λ0(t) takes discrete mass at each failure time ti. Thus, the dimension of the parameter vector λ0 is equal to the number of unique failure times. This converts the semiparametric Cox model to a parametric model, but it introduces a major challenge since standard asymptotic theory for the maximum likelihood estimators (MLEs) may not apply due to the infinitely dimensional nature of λ0.MLEs of the model parameters can either be obtained by a direct maximization of the observed data log likelihood or by using an EM algorithm. Since the observed data log likelihood involves an intractable integral, a direct maximization is often based on numerical integration techniques such as the Gaussian Hermite quadrature or Monte Carlo methods. These methods, however, can be quite computationally intensive if the dimension of the unobservable random effectsai is not low. The EM algorithm is known for its stability and generality, so it is widely used for likelihood inference of joint models [1, 3, 11, 23]. Since the E-step of an EM algorithm still involves an intractable integral, Monte Carlo methods or Laplacian approximations are often used to approximate the conditional expectation in the E-step. In the M-step, the Newton-Raphson method is often used.Hsieh et al. [24] noted that standard errors for estimators of the parameters (α,β,σ) based on the Fisher information matrix may be problematic, because of the semiparametric nature of the joint model. They recommended a bootstrap method to obtain standard errors.


## 4.2. Computational Issues

A main challenge in the likelihood inference for joint models is the computational complexity, since numerical methods or Monte Carlo methods can be very computationally intensive when the dimension of the random effectsai is not small. Moreover, convergence of the EM algorithms can sometimes be an issue. Tsiatis and Davidian [1] and Tsiatis and Davidian [25] proposed an alternative approach, a so-called conditional score method, which makes no distributional assumption for the random effects ai. Their method treats ai as “nuisance parameters” and conditions on an appropriate “sufficient statistic” for ai. Conditioning on this sufficient statistic would remove the dependence of the conditional distribution on the random effects ai. This approach leads to a set of unbiased estimating equations, similar to the usual partial likelihood score equations or generalized estimating equation (GEE). The resulting estimates are, under certain regularity conditions, consistent and asymptotically normal, although they may not be the most efficient. Moreover, their method is relatively simple to implement. Song et al. [26] considered an alternative approach to relax the normality assumption of ai. They assume that the random effects follow distributions in a class of smooth density family, including the standard normal density as a special case. They use the likelihood method for inference via an EM algorithm, but the computation is considerably more intensive than the conditional score method. Song and Wang [27] also proposed a local corrected score estimator and a local conditional score estimator, for which no distributional assumptions are needed for the underlying true covariates.Approximate but computationally more efficient methods for joint models have also appeared in the literature, such as those based on Laplace approximations (e.g., [11, 12, 28]). When the dimension of the integration in the likelihood for joint models is high, the Laplace approximations offer considerably computational advantages over numerical or Monte Carlo methods. Note that the order of the Laplace approximation error is O(mi-1), which cannot be made arbitrarily accurate for a given dataset, where mi is the number of within-individual measurements for individual i. Therefore, the Laplace approximation works well if the number of within-individual measurements is large. Approximate methods based on Taylor series approximations are similar to the Laplace approximation; that is, their performance improves as mi increases.Rizopoulos et al. [11] proposed to use the full exponential Laplace approximation in the E-step of the EM algorithm. Compared to the standard (first-order) Laplace approximation, the full exponential Laplace approximate method has approximation error of order O(mi-2) and requires a much smaller number of within-individual longitudinal measurements to produce reliable results. Lee et al. [28] suggested second-order Laplace approximations. However, these Laplace approximation methods cannot control the magnitude of the approximation errors, unlike Gaussian quadrature or Monte Carlo integration techniques.


## 5. Bayesian Methods

Bayesian joint models have also been studied by various authors, including Faucett and Thomas [29], Xu and Zeger [15], Wang and Taylor [5], Law et al. [30], Ibrahim et al. [16], and Huang et al. [13]. Joint models may contain many unknown parameters, which may lead to potential problems in inference. A main advantage of Bayesian methods is that they can borrow additional information from similar studies or from experts and incorporate this information in the current analysis, in the forms of prior distributions for the current model parameters. Thus, Bayesian methods can be very useful for inference of joint models.For Bayesian joint models, the model parameters are assumed to follow some prior distributions, and inference is then based on the posterior distribution given the observed data. Letθ denote the collection of unknown parameters in the joint model, and let f(θ) denote the prior distribution. Let 𝒟={(ti,δi,zi,xi),i=1,2,…,n} denote all observed data. The joint posterior distribution for all unknown parameters θ and random effects a={ai,i=1,…,n} is then given by(5.1)f(θ,a∣D)∝∏i=1n[f(ti,δi∣zi*,xi,θ)f(zi∣ai,θ)f(ai∣A)]f(θ∣θ0),
where θ0 are known hyperparameters. Bayesian inference is then based on Monte Carlo samples drawn from the posterior distribution f(θ,a∣𝒟) using an MCMC algorithm such as the Gibbs sampler. For example, the posterior means and variances of the parameters can be estimated based on these Monte Carlo samples, and Bayesian inference can then be based on these estimated posterior means and variances. This Monte Carlo sampling can be done using the publically available WinBUGS software [31], which is quite general, flexible, and easy to use.Like other Bayesian methods, it is desirable to check if the final results are sensitive to the choices of prior distributions. Sometimes, in the absence of prior information, noninformative priors or flat priors may be desirable.


## 6. Other Joint Models

In the previous sections, we have focused on joint models based on a Cox model for right-censored survival data and a LME model for longitudinal data. Other models for survival data and longitudinal data can also be considered in joint models. For example, for survival data, we may consider accelerated failure time (AFT) models and models for interval censored data and models for recurrent events. For longitudinal data, nonlinear, generalized linear mixed models or semiparametric/nonparametric mixed models can be utilized. Although the different survival models and longitudinal models can be employed, basic ideas and approaches for inference remain essentially the same. In the following, we briefly review some of these joint models.


### 6.1. Joint Models Based on an LME Model and an AFT Model

In joint modelling of longitudinal and survival data, we can use the AFT model to feature survival data. Here, we focus on an AFT model with measurement errors in time-dependent covariates. For longitudinal data, we again consider LME models for simplicity. The description below is based on Tseng et al. [23]. A semiparametric AFT model can be written in a form similar to the Cox model:(6.1)hi(t)=h0[∫0texp⁡{-zi*(u)β}du]exp⁡{-zi*(t)β},
where hi(t) is the hazard function of the ith individual at time t, h0(t) is the baseline hazard function, and zi*(t) is the unobserved true covariate value at time t. For the observed measurements zi(t), we again consider the LME model (2.2).Tseng et al. [23] proposed a likelihood method using an EM algorithm. The likelihood for all observed data is given by(6.2)L(θ)=∏i=1n∫f(ti,δi∣zi*,h0,β)f(zi∣ai,α,σ2)f(ai∣A)dai,
where f(zi∣ai,α,σ2) and f(ai∣A) are the same as those for (4.1) and(6.3)f(ti,δi∣zi*,h0,β)=[h0{ϕ(ti;θ,ai)}∂ϕ(ti;z̅i*,β)∂ti]δiexp⁡{-∫0ϕ(ti;z̅i*,β)h0(u)du},
where z̅i* denotes the covariate history and ϕ is a known function.Handling the AFT structure in the joint modelling setting is more difficult than for the Cox model, sincef(ti,δi∣zi*,h0,β) is more complicated and the baseline function h0{ϕ(ti;z̅i*,β)} involves unknown quantities (β,α,A,ai), while this is not the case in the Cox model. One cannot use the point mass function with masses assigned to all uncensored survival times ti for the baseline hazard function h0. In other words, in Cox models, the baseline hazard h0 can be represented by a collection of parameters which are point masses, but this approach is not feasible for the AFT model because of its dependence on covariates via function ϕ(ti;z̅i*,β). To circumvent this, Tseng et al. [23] assumed the baseline hazard function h0 to be a step function, taking constant values between two consecutive failure times.Tseng et al. [23] used an Monte Carlo EM algorithm to obtain the MLEs. The framework is similar to that in the previous section. They used a Monte Carlo method to approximate the conditional expectations in the E-step. The M-step involves more complicated computations due to the complicated baseline hazard h0. To obtain the standard errors of the MLEs, the usual asymptotic formula based on Fisher information matrix may be questionable, so they used a bootstrap method.


### 6.2. Joint Models with Interval Censored Survival Data

In the previous sections, we have focused on right censored survival data and assume that either theexact survival times or censoring times are observed. In practice, however, we often cannot observe the exact survival nor censoring times, but we only know that events have occurred over certain time intervals. Such survival data are called interval censored. For simplicity, we assume that all individuals are assessed at the same times. Again, let Si be the time to an event (survival time) for individual i, with observed value si. Let ri=(ri1,…,rim)T be the vector of event indicators such that rij=1 if subject i has an event occurred from time tj-1 to time tj, and let rij=0 otherwise, i=1,2,…,n; j=1,2,…,m. We assume that ri1=0 for all i. Let pij=P(tj-1≤Si<tj), and let(6.4)πij=P(tj-1≤Si<tj∣Si≥tj-1)=1-P(Si≥tj∣Si≥tj-1).
Then, we have pij=(1-πi1)(1-πi2)⋯(1-πi,j-1)πij. The probability function for the event indicator vector ri can be written as(6.5)f(ri)=∏j=1mpijrij=∏j=1mπijrij(1-πij)1-rij,
which is the probability function for a Bernoulli distribution. We can introduce observed error-prone covariate value zi, with true value zi*, and assume(6.6)log⁡{-log⁡(1-πij)}=βTzi*+γj,
where β and γ=(γ1,…,γm)T are unknown parameters. Then, we can write the probability function of ri as f(ri∣zi*,β,γ). Alternatively, we can assume that ri depends on zi* only through the random effects ai and writing the probability function of ri as f(ri∣ai,β,γ).Letf(zi∣ai,α,σ) be the conditional probability density function, given the random effects ai, and f(ai∣A) be the marginal probability density function for ai with covariance matrix A. Let θ denote the collection of all parameters in all models. Then, the likelihood for all the observed data can be written as(6.7)Lo(θ)=∏i=1n[∫f(zi∣ai,α,σ)f(ri∣ai,β,γ)f(ai∣A)dai].
MLE of parameters θ can be obtained by maximizing the observed data likelihood Lo(θ). Because the observed-data likelihood Lo(θ) can be difficult to evaluate due to its involvement of an intractable and possibly high-dimensional integral, one may proceed with Monte Carlo EM algorithms or other computationally more efficient approximate methods.


### 6.3. GLMM and NLME Models for Longitudinal Data

We have focused on LME models for modelling the longitudinal data. Other models for longitudinal data can also be considered. For example, one may consider nonlinear mixed-effects (NLME) models for modelling the longitudinal data in joint models [12, 32]. NLME models are often mechanistic models in the sense that they are typically based on the underlying mechanisms which generate the data. On the other hand, LME models are typically empirical models; that is, they are usually used to approximately describe the observed data without considering possible data-generation mechanisms. Thus, NLME models may be scientifically more desirable if such models exist. Similarly, for nonnormal longitudinal data, generalized linear mixed models (GLMMs) can be considered, which are special nonlinear models but are essentially empirical models as well.When the longitudinal models are nonlinear, the general ideas of the two-stage methods and likelihood methods for joint models can still be applied. The complication is that computation becomes more demanding, because of the nonlinearity of the longitudinal models.


### 6.4. Joint Models with Missing Data

For longitudinal data, missing values are very common. When missing data are nonignorable in the sense that the missingness probability may be related to the missing values or the random effects, the missing data process is often needed to be incorporated in inferential procedures in order to obtain valid results. For likelihood methods, it is straightforward to incorporate missing data mechanisms in joint model inference. However, the computation becomes even more challenging. Wu et al. [12, 32] considered the missing data problems for joint models, using Monte Carlo EM algorithms and Laplace approximations.


## 6.1. Joint Models Based on an LME Model and an AFT Model

In joint modelling of longitudinal and survival data, we can use the AFT model to feature survival data. Here, we focus on an AFT model with measurement errors in time-dependent covariates. For longitudinal data, we again consider LME models for simplicity. The description below is based on Tseng et al. [23]. A semiparametric AFT model can be written in a form similar to the Cox model:(6.1)hi(t)=h0[∫0texp⁡{-zi*(u)β}du]exp⁡{-zi*(t)β},
where hi(t) is the hazard function of the ith individual at time t, h0(t) is the baseline hazard function, and zi*(t) is the unobserved true covariate value at time t. For the observed measurements zi(t), we again consider the LME model (2.2).Tseng et al. [23] proposed a likelihood method using an EM algorithm. The likelihood for all observed data is given by(6.2)L(θ)=∏i=1n∫f(ti,δi∣zi*,h0,β)f(zi∣ai,α,σ2)f(ai∣A)dai,
where f(zi∣ai,α,σ2) and f(ai∣A) are the same as those for (4.1) and(6.3)f(ti,δi∣zi*,h0,β)=[h0{ϕ(ti;θ,ai)}∂ϕ(ti;z̅i*,β)∂ti]δiexp⁡{-∫0ϕ(ti;z̅i*,β)h0(u)du},
where z̅i* denotes the covariate history and ϕ is a known function.Handling the AFT structure in the joint modelling setting is more difficult than for the Cox model, sincef(ti,δi∣zi*,h0,β) is more complicated and the baseline function h0{ϕ(ti;z̅i*,β)} involves unknown quantities (β,α,A,ai), while this is not the case in the Cox model. One cannot use the point mass function with masses assigned to all uncensored survival times ti for the baseline hazard function h0. In other words, in Cox models, the baseline hazard h0 can be represented by a collection of parameters which are point masses, but this approach is not feasible for the AFT model because of its dependence on covariates via function ϕ(ti;z̅i*,β). To circumvent this, Tseng et al. [23] assumed the baseline hazard function h0 to be a step function, taking constant values between two consecutive failure times.Tseng et al. [23] used an Monte Carlo EM algorithm to obtain the MLEs. The framework is similar to that in the previous section. They used a Monte Carlo method to approximate the conditional expectations in the E-step. The M-step involves more complicated computations due to the complicated baseline hazard h0. To obtain the standard errors of the MLEs, the usual asymptotic formula based on Fisher information matrix may be questionable, so they used a bootstrap method.


## 6.2. Joint Models with Interval Censored Survival Data

In the previous sections, we have focused on right censored survival data and assume that either theexact survival times or censoring times are observed. In practice, however, we often cannot observe the exact survival nor censoring times, but we only know that events have occurred over certain time intervals. Such survival data are called interval censored. For simplicity, we assume that all individuals are assessed at the same times. Again, let Si be the time to an event (survival time) for individual i, with observed value si. Let ri=(ri1,…,rim)T be the vector of event indicators such that rij=1 if subject i has an event occurred from time tj-1 to time tj, and let rij=0 otherwise, i=1,2,…,n; j=1,2,…,m. We assume that ri1=0 for all i. Let pij=P(tj-1≤Si<tj), and let(6.4)πij=P(tj-1≤Si<tj∣Si≥tj-1)=1-P(Si≥tj∣Si≥tj-1).
Then, we have pij=(1-πi1)(1-πi2)⋯(1-πi,j-1)πij. The probability function for the event indicator vector ri can be written as(6.5)f(ri)=∏j=1mpijrij=∏j=1mπijrij(1-πij)1-rij,
which is the probability function for a Bernoulli distribution. We can introduce observed error-prone covariate value zi, with true value zi*, and assume(6.6)log⁡{-log⁡(1-πij)}=βTzi*+γj,
where β and γ=(γ1,…,γm)T are unknown parameters. Then, we can write the probability function of ri as f(ri∣zi*,β,γ). Alternatively, we can assume that ri depends on zi* only through the random effects ai and writing the probability function of ri as f(ri∣ai,β,γ).Letf(zi∣ai,α,σ) be the conditional probability density function, given the random effects ai, and f(ai∣A) be the marginal probability density function for ai with covariance matrix A. Let θ denote the collection of all parameters in all models. Then, the likelihood for all the observed data can be written as(6.7)Lo(θ)=∏i=1n[∫f(zi∣ai,α,σ)f(ri∣ai,β,γ)f(ai∣A)dai].
MLE of parameters θ can be obtained by maximizing the observed data likelihood Lo(θ). Because the observed-data likelihood Lo(θ) can be difficult to evaluate due to its involvement of an intractable and possibly high-dimensional integral, one may proceed with Monte Carlo EM algorithms or other computationally more efficient approximate methods.


## 6.3. GLMM and NLME Models for Longitudinal Data

We have focused on LME models for modelling the longitudinal data. Other models for longitudinal data can also be considered. For example, one may consider nonlinear mixed-effects (NLME) models for modelling the longitudinal data in joint models [12, 32]. NLME models are often mechanistic models in the sense that they are typically based on the underlying mechanisms which generate the data. On the other hand, LME models are typically empirical models; that is, they are usually used to approximately describe the observed data without considering possible data-generation mechanisms. Thus, NLME models may be scientifically more desirable if such models exist. Similarly, for nonnormal longitudinal data, generalized linear mixed models (GLMMs) can be considered, which are special nonlinear models but are essentially empirical models as well.When the longitudinal models are nonlinear, the general ideas of the two-stage methods and likelihood methods for joint models can still be applied. The complication is that computation becomes more demanding, because of the nonlinearity of the longitudinal models.


## 6.4. Joint Models with Missing Data

For longitudinal data, missing values are very common. When missing data are nonignorable in the sense that the missingness probability may be related to the missing values or the random effects, the missing data process is often needed to be incorporated in inferential procedures in order to obtain valid results. For likelihood methods, it is straightforward to incorporate missing data mechanisms in joint model inference. However, the computation becomes even more challenging. Wu et al. [12, 32] considered the missing data problems for joint models, using Monte Carlo EM algorithms and Laplace approximations.


## 7. Example and Simulation

### 7.1. An Illustrating Example

As an illustration, we consider an AIDS dataset which includes 46 HIV infected patients receiving an anti-HIV treatment. Viral load (i.e., plasma HIV RNA) and CD4 cell count were repeatedly measured during the treatment. The number of viral load measurements for each individual varies from 4 to 10. It is known that CD4 is measured with substantial errors. About 11% viral load measurements are below the detection limit after the initial period of viral decay. We call the viral load below the detection limit as “viral suppression.” We wish to check if the initial CD4 trajectories are predictive for the time to viral suppression.Letsi be the time to viral suppression, that is, the time from the start of the treatment to the first scheduled time when the viral load drops below the detection limit. The viral suppression times for patients whose viral loads never dropped below detection limit may be regarded as right censored, with the study end time as the censoring time. We employ two models, the Cox model (2.1) or the AFT model (6.1), to feature the time to viral suppression, where zi*(t) is the unobserved true CD4 cell count at time t for individual i. We do not consider additional covariates here.To address the measurement error in the time-dependent covariate CD4 cell count, we use the LME model to model the CD4 trajectories:(7.1)CD4ij=(α0+ai1)+(α1+ai2uij)+ϵij,
where the parameters, random effects and random errors, and the assumed distributions are the same as those described in the previous sections. The fixed parameters (α0,α1) are population intercept and slope of the CD4 process and (ai1,ai2) are individual variations from the population averages. To avoid very small (large) estimates, which may be unstable, we standardize the CD4 cell counts and rescale the original time t (in days) so that the new time scale is between 0 and 1. We estimate the model parameters using the joint model method and the two-stage method with/out bootstrap standard error correction. The number B=500 of the bootstrap samples is taken. For the joint model method, we consider the Cox model (2.1) and the AFT model (6.1) with h0(·) being the Weibull baseline risk function. On the other hand, only the Cox model (2.1) is employed for the two-stage method for comparison. These analyses may be implemented by using the functions coxph(), lme(), and jointModel() in R software.Table1 presents the resulting estimates of main parameters of interest and their standard errors. We can see from Table 1 that under either the Cox or the AFT survival models, both the two-stage and the joint model methods produce similar estimates for the covariate (CD4) longitudinal model. However, the two-stage method may underestimate the standard errors of the parameter estimates since it does not incorporate the survival information in the estimation procedure. Note that the parameters in the Cox model and the AFT model have different interpretations, due to different model formulations, so they are not directly comparable.Table 1
Analyses of the AIDS data under different models.
ModelMethodβ0β1α0α1σCox modelTwo-stageEstimate—0.315−0.2331.3450.605SE—0.2080.1130.154—P value—0.1290.040<0.001—BSE—0.237———Cox modelJoint modelEstimate—0.648−0.2011.3420.603SE—0.2340.1350.162—P value—0.0060.137<0.001—AFT modelJoint modelEstimate0.168−0.487−0.2371.3410.604SE0.2600.2890.1250.156—P value0.5170.0910.059<0.001—SE: standard error; BSE: bootstrap standard error. TheP values are from the Wald tests for testing if the corresponding parameters are zero or not.The parameterβ1 in the survival models measures the effect of the true time-dependent covariate CD4 values on event times, so its estimate and the associated P value can be used to check if the true CD4 values are predictive for the times to viral suppression. Since the covariate CD4 is measured with errors, addressing the measurement error is the main focus of joint modelling in this application. Thus, the estimation of β1 is of primary interest. For the joint model method, under either the Cox or the AFT models, there is some evidence that covariate CD4 is associated with the time to viral suppression, after measurement error has been addressed. It is seen that evidence of significance of the covariate effect is the strongest under the Cox model. On the other hand, the two-stage method may severely underestimate the covariate CD4 effect (the small value of β̂1). Moreover, the naive two-stage method underestimates the standard error of β̂1, due to failing to incorporate the estimating uncertainty from the first step. This underestimation of standard error is somewhat corrected by the bootstrap method.


### 7.2. A Simulation Study

In this section, we conduct a simulation study to compare the joint model method and the two-stage method with/out bootstrap standard error correction. We generate 500 datasets from the time-dependent covariate CD4 process (7.1) in the example of the previous section and the Cox model (2.1) with constant baseline hazard function λ0(t)≡1 with emphasis on the effect of the time-dependent covariate. The measurement time points used in the simulation are the same as those in the example of the previous section. The true values of model parameters, given in Table 2, are similar to those in the example, and the variance-covariance matrix A of the random effect ai is set to be diagonal. Again, we take B=500 bootstrap samples for the bootstrap standard error in the two-stage method.Table 2
Comparison of the two-stage Method and the joint likelihood method via a simulation study.
MethodParameterβ1α0α1σA11A22True value0.6−0.21.30.60.50.3Two-stageEst0.183−0.1831.3030.5980.5010.353ESE0.2160.1110.1640.0250.1140.209ASE0.2010.1110.164———BSE0.250—————Bias−0.4170.0170.003−0.002−0.0010.053MSE0.2210.0130.0270.00070.0130.046CR42.895.694.4———Joint modelEst0.6004−0.1751.2960.5980.4920.321ESE0.2560.1030.1610.0200.0920.156ASE0.2490.0990.163———Bias0.00040.025−0.004−0.002−0.0080.021MSE0.0660.0110.0260.00040.0080.025CR95.695.895.2———In Table2, we report the simulation results of averages parameter estimates (Est), empirical standard errors (ESE), averages of asymptotic standard errors (ASE), average of bootstrap standard errors (BSE), empirical biases (Bias), mean square errors (MSE), and coverage rates (CR) for the 95% confidence intervals. We can see from Table 2 that the two-stage and the joint model methods produce similar parameter estimates (close to true parameters) except the one for the covariate effect β1. In particular, the estimate β̂1 based on the joint model method is very close to its true value, while the estimate based on the two-stage method is about one-third of its true value, which indicates that the two-stage method may underestimate the time-dependent covariate effect severely. The joint model method provides smaller mean square errors and more reasonable coverage rates for the 95% confidence intervals than the two-stage method. Moreover, the two-stage method may underestimate the standard deviation of β̂1 and the bootstrap correction on this standard error seems plausible.From the above results, we see that the joint likelihood method produces less biased estimates and more reliable standard errors than the two-stage method. These results have important implications. For example, if one uses Wald-type tests for model selection, the likelihood method would give more reliable results. However, two-stage methods are generally simpler and computationally quicker to output estimates than likelihood methods. We can also compare the two methods with Bayesian methods. Note that, however, Bayesian methods are equivalent to the likelihood method when noninformative priors are used. We expect that Bayesian methods have similar performance to likelihood methods.


## 7.1. An Illustrating Example

As an illustration, we consider an AIDS dataset which includes 46 HIV infected patients receiving an anti-HIV treatment. Viral load (i.e., plasma HIV RNA) and CD4 cell count were repeatedly measured during the treatment. The number of viral load measurements for each individual varies from 4 to 10. It is known that CD4 is measured with substantial errors. About 11% viral load measurements are below the detection limit after the initial period of viral decay. We call the viral load below the detection limit as “viral suppression.” We wish to check if the initial CD4 trajectories are predictive for the time to viral suppression.Letsi be the time to viral suppression, that is, the time from the start of the treatment to the first scheduled time when the viral load drops below the detection limit. The viral suppression times for patients whose viral loads never dropped below detection limit may be regarded as right censored, with the study end time as the censoring time. We employ two models, the Cox model (2.1) or the AFT model (6.1), to feature the time to viral suppression, where zi*(t) is the unobserved true CD4 cell count at time t for individual i. We do not consider additional covariates here.To address the measurement error in the time-dependent covariate CD4 cell count, we use the LME model to model the CD4 trajectories:(7.1)CD4ij=(α0+ai1)+(α1+ai2uij)+ϵij,
where the parameters, random effects and random errors, and the assumed distributions are the same as those described in the previous sections. The fixed parameters (α0,α1) are population intercept and slope of the CD4 process and (ai1,ai2) are individual variations from the population averages. To avoid very small (large) estimates, which may be unstable, we standardize the CD4 cell counts and rescale the original time t (in days) so that the new time scale is between 0 and 1. We estimate the model parameters using the joint model method and the two-stage method with/out bootstrap standard error correction. The number B=500 of the bootstrap samples is taken. For the joint model method, we consider the Cox model (2.1) and the AFT model (6.1) with h0(·) being the Weibull baseline risk function. On the other hand, only the Cox model (2.1) is employed for the two-stage method for comparison. These analyses may be implemented by using the functions coxph(), lme(), and jointModel() in R software.Table1 presents the resulting estimates of main parameters of interest and their standard errors. We can see from Table 1 that under either the Cox or the AFT survival models, both the two-stage and the joint model methods produce similar estimates for the covariate (CD4) longitudinal model. However, the two-stage method may underestimate the standard errors of the parameter estimates since it does not incorporate the survival information in the estimation procedure. Note that the parameters in the Cox model and the AFT model have different interpretations, due to different model formulations, so they are not directly comparable.Table 1
Analyses of the AIDS data under different models.
ModelMethodβ0β1α0α1σCox modelTwo-stageEstimate—0.315−0.2331.3450.605SE—0.2080.1130.154—P value—0.1290.040<0.001—BSE—0.237———Cox modelJoint modelEstimate—0.648−0.2011.3420.603SE—0.2340.1350.162—P value—0.0060.137<0.001—AFT modelJoint modelEstimate0.168−0.487−0.2371.3410.604SE0.2600.2890.1250.156—P value0.5170.0910.059<0.001—SE: standard error; BSE: bootstrap standard error. TheP values are from the Wald tests for testing if the corresponding parameters are zero or not.The parameterβ1 in the survival models measures the effect of the true time-dependent covariate CD4 values on event times, so its estimate and the associated P value can be used to check if the true CD4 values are predictive for the times to viral suppression. Since the covariate CD4 is measured with errors, addressing the measurement error is the main focus of joint modelling in this application. Thus, the estimation of β1 is of primary interest. For the joint model method, under either the Cox or the AFT models, there is some evidence that covariate CD4 is associated with the time to viral suppression, after measurement error has been addressed. It is seen that evidence of significance of the covariate effect is the strongest under the Cox model. On the other hand, the two-stage method may severely underestimate the covariate CD4 effect (the small value of β̂1). Moreover, the naive two-stage method underestimates the standard error of β̂1, due to failing to incorporate the estimating uncertainty from the first step. This underestimation of standard error is somewhat corrected by the bootstrap method.


## 7.2. A Simulation Study

In this section, we conduct a simulation study to compare the joint model method and the two-stage method with/out bootstrap standard error correction. We generate 500 datasets from the time-dependent covariate CD4 process (7.1) in the example of the previous section and the Cox model (2.1) with constant baseline hazard function λ0(t)≡1 with emphasis on the effect of the time-dependent covariate. The measurement time points used in the simulation are the same as those in the example of the previous section. The true values of model parameters, given in Table 2, are similar to those in the example, and the variance-covariance matrix A of the random effect ai is set to be diagonal. Again, we take B=500 bootstrap samples for the bootstrap standard error in the two-stage method.Table 2
Comparison of the two-stage Method and the joint likelihood method via a simulation study.
MethodParameterβ1α0α1σA11A22True value0.6−0.21.30.60.50.3Two-stageEst0.183−0.1831.3030.5980.5010.353ESE0.2160.1110.1640.0250.1140.209ASE0.2010.1110.164———BSE0.250—————Bias−0.4170.0170.003−0.002−0.0010.053MSE0.2210.0130.0270.00070.0130.046CR42.895.694.4———Joint modelEst0.6004−0.1751.2960.5980.4920.321ESE0.2560.1030.1610.0200.0920.156ASE0.2490.0990.163———Bias0.00040.025−0.004−0.002−0.0080.021MSE0.0660.0110.0260.00040.0080.025CR95.695.895.2———In Table2, we report the simulation results of averages parameter estimates (Est), empirical standard errors (ESE), averages of asymptotic standard errors (ASE), average of bootstrap standard errors (BSE), empirical biases (Bias), mean square errors (MSE), and coverage rates (CR) for the 95% confidence intervals. We can see from Table 2 that the two-stage and the joint model methods produce similar parameter estimates (close to true parameters) except the one for the covariate effect β1. In particular, the estimate β̂1 based on the joint model method is very close to its true value, while the estimate based on the two-stage method is about one-third of its true value, which indicates that the two-stage method may underestimate the time-dependent covariate effect severely. The joint model method provides smaller mean square errors and more reasonable coverage rates for the 95% confidence intervals than the two-stage method. Moreover, the two-stage method may underestimate the standard deviation of β̂1 and the bootstrap correction on this standard error seems plausible.From the above results, we see that the joint likelihood method produces less biased estimates and more reliable standard errors than the two-stage method. These results have important implications. For example, if one uses Wald-type tests for model selection, the likelihood method would give more reliable results. However, two-stage methods are generally simpler and computationally quicker to output estimates than likelihood methods. We can also compare the two methods with Bayesian methods. Note that, however, Bayesian methods are equivalent to the likelihood method when noninformative priors are used. We expect that Bayesian methods have similar performance to likelihood methods.


## 8. Discussion

We have provided a brief review of common joint models and methods for inference. In practice, when we need to consider a longitudinal process and an event process and suspect that the two processes may be associated, such as survival models with time-dependent covariates or longitudinal models with informative dropouts, it is important to use joint model methods for inference in order to avoid biased results. The literature on model selection for joint models is quite limited. In practice, the best longitudinal model can be selected based on the observed longitudinal data, and the best survival model can be selected based on the survival data, using standard model selection procedures for these models. Then, we specify reasonable link between the two models, such as shared random effects. To choose methods for inference, the joint likelihood method generally produces most reliable resultsif the assumed models and distributions are correct. On the other hand, the two-stage methods may be computationally simpler, and many existing models and methods for longitudinal data and survival data can be easily adapted. However, two-stage methods may not completely eliminate the biases in parameter estimates in some cases.When the longitudinal covariate process terminates at event times, that is, when the longitudinal values are unavailable at and after the event times such as deaths or dropouts, the covariates are sometimes calledinternal time-dependent covariates. Sometimes, however, longitudinal covariate information is available at and after the event times. For example, CD4 measurements may be still available after patients have been diagnosed with AIDS. Such covariates are sometimes called external covariates. In joint models, it is important to distinguish internal and external time-dependent covariates. In particular, for internal time-dependent covariates, joint models are more desirable since separate analysis in this case may lead to more severe bias.Survival models with measurement errors in time-dependent covariates have received much attention in the joint models literature. Another common situation is longitudinal models with informative dropouts, in which survival models can be used to model the dropout process. Both situations focus on characterizing the association between the longitudinal and survival processes. Some authors have also considered joint models in which the focus is on more efficient inference of the survival model, using longitudinal data as auxiliary information [15, 33, 34] or assume that the longitudinal process and the survival process are governed by a common latent process [4]. Nathoo and Dean [7] considered an interesting joint model in which an NLME model is used to model tree growth, with spatial correlation incorporated.Joint models can also be extended tomultivariate cases, in which more than one longitudinal processes and more than one event processes can be modelled simultaneously. Extensions are often conceptually straightforward, but computation and implementation can be more tedious than univariate cases. See Henderson et al. [4], Xu and Zeger [35], and Song et al. [26].Zeng and Cai [36] derived some asymptotic results for maximum likelihood estimators in joint analysis of longitudinal and survival data. They showed the consistency of the maximum likelihood estimators, derived their asymptotic distributions, and showed that the maximum likelihood estimators in joint analysis are semiparametrically efficient.Although there has been extensive research in joint models in the last two decades and the importance of joint models has been increasingly recognized, joint models are still not widely used in practice. A main reason is perhaps lack of software. Recently, Dimitris Rizopoulos has developed an R package calledJM that can be used to fit joint models with normal longitudinal responses and event times under a maximum likelihood approach. Various options for the survival model and optimization/integration algorithms are provided, such as Cox models and AFT models for survival data and the Gauss-Hermite integration methods and Laplace approximations.



---

*Source: 640153-2011-12-07.xml*