Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Ann Appl Stat. Author manuscript; available in PMC 2011 September 19.
Published in final edited form as:
Ann Appl Stat. 2010 September 1; 4(3): 1517–1532.
doi:  10.1214/10-AOAS339
PMCID: PMC3175771



In many medical studies, patients are followed longitudinally and interest is on assessing the relationship between longitudinal measurements and time to an event. Recently, various authors have proposed joint modeling approaches for longitudinal and time-to-event data for a single longitudinal variable. These joint modeling approaches become intractable with even a few longitudinal variables. In this paper we propose a regression calibration approach for jointly modeling multiple longitudinal measurements and discrete time-to-event data. Ideally, a two-stage modeling approach could be applied in which the multiple longitudinal measurements are modeled in the first stage and the longitudinal model is related to the time-to-event data in the second stage. Biased parameter estimation due to informative dropout makes this direct two-stage modeling approach problematic. We propose a regression calibration approach which appropriately accounts for informative dropout. We approximate the conditional distribution of the multiple longitudinal measurements given the event time by modeling all pairwise combinations of the longitudinal measurements using a bivariate linear mixed model which conditions on the event time. Complete data are then simulated based on estimates from these pairwise conditional models, and regression calibration is used to estimate the relationship between longitudinal data and time-to-event data using the complete data. We show that this approach performs well in estimating the relationship between multivariate longitudinal measurements and the time-to-event data and in estimating the parameters of the multiple longitudinal process subject to informative dropout. We illustrate this methodology with simulations and with an analysis of primary biliary cirrhosis (PBC) data.

Key words and phrases: Joint models, shared random parameter models, informative dropout, regression calibration

1. Introduction

Recently, many studies collect longitudinal data on a panel of biomarkers, and interest is on assessing the relationship between these biomarkers and time to an event. For example, Allen et al. (2007) examined the relationship between five longitudinally collected cytokines measured from serum plasma and survival. Interest focused on whether the values of these multiple cytokines are associated with survival. In another example, patients with primary biliary cirrhosis are followed longitudinally and interest is on examining whether multiple longitudinally biomarkers are prognostic for a poor clinical outcome. Important features in studies of this type are that there may be a relatively large number of biomarkers and that these biomarkers are subject to sizable measurement error due to laboratory error and biological variation.

Various authors have proposed joint modeling approaches for a single longitudinal measurement and time-to-event data [Tsiatis, DeGruttola and Wulfsohn (1995); Wulfsohn and Tsiatis (1997); Tsiatis and Davidian (2004); Henderson, Diggle and Dobson (2000), among others]. There is also limited work on joint models for a few longitudinal measurements and time-to-event data [Xu and Zeger (2001a, 2001b); Huang et al. (2001); Song, Davidian and Tsiatis (2002); Ibrahim, Chen and Sinha (2004); Brown, Ibrahim and DeGruttola (2005); Chi and Ibrahim (2006)]. However, these methods are difficult to implement when the number of longitudinal biomarkers is large since most of these approaches require integrating over the vector of all random effects to evaluate the joint likelihood of the multivariate longitudinal and time-to-event data. This paper proposes an approach for jointly modeling multivariate longitudinal and discrete time-to-event data which easily accommodates many longitudinal biomarkers.

Fieuws and Verbeke (2005) and Fieuws, Verbeke and Molenberghs (2007) have proposed an approach for modeling multivariate longitudinal data whereby all possible pairs of longitudinal data are separately modeled and are then combined in a final step. We use a similar approach along with a recent regression calibration approach for jointly modeling a single series of longitudinal measurements and time-to-event data [Albert and Shih (2009)] to implement the joint modeling approach proposed in this paper. Recently, Fieuws et al. (2008) have proposed a discriminant analysis based approach for using multivariate longitudinal profiles to predict renal graft failure. At the end of their discussion, they mention that a more elegant approach, which has not yet been developed, would involve a joint model for the many longitudinal profiles and time-to-event data. This paper presents such an approach.

We describe the approach in Section 2. We show the advantages of this approach using simulation in Section 3. We illustrate the methodology with an analysis of primary biliary cirrhosis data (PBC) in which we simultaneously examine the relationship between multiple longitudinal biomarker in Section 4. A discussion follows in Section 5.

2. Modeling approach

Define Ti to be a discrete event-time which can take on discrete values tj, j = 1, 2, …, J, and Yij to be a binary indictor of whether the ith patient is dead at time tj. Then Ji=j=1J(1Yij)=JYi·, where Yi·=j=1JYij indicates the number of follow-up measurements before the event or the end of follow-up at time tJ. Longitudinal measurements are measured at times t1, t2, …, tJi. Denote X1i = (X1i1, X1i2, …, X1iJi)′, X2i = (X2i1, X2i2, …, X2iJi)′, …, XPi = (XPi1, XPi2, …, XPiJi)′ as the P biomarkers measured repeatedly at j = 1, 2, …, Ji time points. Further, define Xpi*=(Xpi1*,Xpi2*,,XpiJi*) as the longitudinal measurements without measurement error for the pth biomarker and Xi*=(X1i*,X2i*,,XPi*). We consider a joint model for multivariate longitudinal and discrete time-to-event data in which the discrete event time distribution is modeled as a linear function of previous true values of the biomarkers without measurement error on the probit scale. Specifically,

P(Yij=1[mid ]Yi(j1)=0;Xi*)=Φ(α0j+p=1PαpXpi(j1)*),

where i = 1, 2, …, I, j = 2, 3, …, Ji, Yi1 is taken as 0, α0j governs the baseline discrete event time distribution and αp measures the effect of the pth biomarker (p = 1, 2, …, P) at time tj–1 on survival at time tj. Specifically, (1) allows for examining the effect of multiple “true” biomarker values at time j – 1 on the probability of an event between the (j – 1) th and jth time point.

The longitudinal data is modeled assuming that the fixed and random effect trajectories are linear. Specifically, the multivariate longitudinal biomarkers can be modeled as




where βp0 and βp1 are the fixed effect intercept and slope for the pth biomarker, and γpi0 and γpi1 are the random effect intercept and slope for the pth biomarker on the ith individual. Denote β = (β10, β11, β20, β21, …, βP0, βP1)′ and γi = (γ1i0, γ1i1, γ2i0, γ2i1, …, γPi0, γPi1)′. We assume that γi is normally distributed with mean 0 and variance Σγ, where Σγ is a 2P by 2P dimensional variance matrix, and εpij are independent error terms which are assumed to be normally distributed with mean 0 and variance σpε2(p=1,2,,P).

Alternative to (1), where the probability of an event over an interval depends on the true biomarker values at the beginning of the interval, the event-time process could be adapted to depend on the random effects of the multivariate longitudinal process [e.g., γpi1 can replace Xpi(j1)* in (1)].

2.1. Difficulty in joint estimation

Conceptually, model (1)–(2) can be estimated by maximizing the likelihood

L=[product]i=1Iγi(...){[product]p=1Ph(Xpi[mid ]γpi0,γpi1)}×{[product]j=2Ji(1rij)}(ri(Ji+1))Ji<Jf(γi)dγi,

where rij = P (Yij = 1|Yi(j − 1) = 0), h (Xpi|γpi0, γpi1) is the product of Ji univariate normal density functions each with mean Xpij* and variance σpε2, and f (γ) is a multivariate normal density with mean zero and variance Σγ. When P = 1, (4) can be maximized by numerical integration techniques such as a simple trapezoidal rule or Gaussian quadrature [Abramowitz and Stegun (1974)]. However, these methods are not feasible for even a few longitudinal biomarkers. Alternative Monte Carlo methods such as Monte Carlo EM [Wei and Tanner (1990)] are possible, but these methods do not perform well for even moderately high dimensional random effects (say, P > 2). In the next subsection we develop an alternative approach which is easy to implement with a large number of longitudinal biomarkers.

2.2. Estimation

We propose a two stage regression calibration approach for estimation, which can be described as follows. In the first stage, multivariate linear mixed models can be used to model the longitudinal data. In the second stage, the time-to-event model is estimated by replacing the random effects with corresponding empirical Bayes estimates. There are three problems with directly applying this approach. First, estimation in the first stage is complicated by the fact that simply fitting multivariate linear mixed models results in bias due to informative dropout; this is demonstrated by Albert and Shih (2009) for the the case of P = 1. Second, as described in Section 2.1, parameter estimation for multivariate linear mixed models can be computationally difficult when the number of longitudinal measurements (P) is even moderately large. Third, calibration error in the empirical Bayes estimation needs to be accounted for in the time-to-event model. The proposed approach will deal with all three of these problems.

The bias from informative dropout is a result of differential follow-up whereby the longitudinal process is related to the length of follow-up. That is, in (2)–(3), patients with large values of Xpij* are more likely to have an early event when αp > 0 for p = 1, 2, …, P. There would be no bias if all J follow-up measurements were observed on all patients. As proposed by Albert and Shih (2009) for univariate longitudinal data, we can avoid this bias by generating complete data from the conditional distribution of Xi = (X1i, X2i, …, XPi) given Ti, denoted as Xi|Ti. Since Xi|Ti under model (2)–(3) does not have a tractable form, we propose a simple approximation for this conditional distribution. Under model (2)–(3), the distribution of Xi|Ti can be expressed as

P(Xi[mid ]Ti)=h(Xi[mid ]γi,Ti)g(γi[mid ]Ti)dγi.

Since Ti and the values of Xi are conditional independent given γi, h (Xi|γi, Ti) = h (Xi|γi), where h(Xi[mid ]γi)=[product]p=1Ph(Xpi[mid ]γpi0,γpi1). The distribution of Xi|Ti can be expressed as a multivariate linear mixed model if we approximate g(γi|Ti) by a normal distribution. Under the assumption that g(γi|Ti) is normally distributed with mean μTi = μ01Ti, μ11Ti, μ02Ti, μ12Ti, …, μ0PTi, μ1PTi)′ and variance γTi*, and by rearranging mean structure parameters in the integrand of (5) so that the random effects have mean zero, Xi|Ti corresponds to the following multivariate linear mixed model:

Xpij[mid ](Ti,γip0Ti*,γip1Ti*)=βp0Ti*+βp1Ti*tj+γip0Ti*+γip1Ti*tj+εpij*,

where i = 1, 2, …, I, j = 1, 2, …, Ji, and p = 1, 2, …, P. The parameters βp0Ti* and βp1Ti* are intercept and slope parameters for the pth longitudinal measurement and for patients who have an event time at time Ti or who are censored at time tJ. In addition, the associated random effects γiTi*=(γi10Ti*,γi11Ti*,γi20Ti*,γi21Ti*,,γiP0Ti*,γiP1Ti*) are multivariate normal with mean 0 and variance γTi*, and the residuals εpij* are assumed to have an independent normal distribution with mean zero and variance σεp*2. Thus, this conditional model involves estimating separate fixed effect intercept and slope parameters for each potential event-time and for subjects who are censored at time tJ. Likewise, separate random effects distributions are estimated for each of these discrete time points. For example, the intercept and slope fixed-effect parameters for the pth biomarker for those patients who have an event at time Ti = t3 is βp0t3* and βp1t3*, respectively. Further, the intercept and slope random effects for all P biomarkers on those patients who have an event at time Ti = t3, γit3*, is multivariate normal with mean 0 and variance γt3*. A similar approximation has been proposed by Albert and Shih (2009) for univariate longitudinal data (P = 1).

Recall that by generating complete data from (6) we are able to avoid the bias due to informative dropout. However, when P is large, direct estimation of model (6) is difficult since the number of elements in γTi* grows quadratically with P. For example, the dimension of the variance matrix γTi* is 2P by 2P for P longitudinal biomarkers. Fieuws and Verbeke (2005) proposed estimating the parameters of multivariate linear mixed models by formulating bivariate linear mixed models on all possible pairwise combinations of longitudinal measurements. In the simplest approach, they proposed fitting bivariate linear mixed models on all (P2) combinations of longitudinal biomarkers and averaging “overlapping” or duplicate parameter estimates. Thus, we estimate the parameters in the fully specified model (6) by fitting (P2) bivariate longitudinal models that only include pairs of longitudinal markers. Fieuws and Verbeke (2005) demonstrated with simulations that there is little efficiency loss using their approach relative to a full maximum-likelihood based approach. Fitting these bivariate models is computationally feasible since only four correlated random effects are contained in each model. (That is, γTi* is a four by four-dimensional matrix for each discrete event-time Ti.) Duplicate estimates of fixed effects and random-effect variances from all pairwise bivariate models are averaged to obtain final parameter estimates of the fully specified model (6). For example, when P = 4 there are (P − 1) = 3 estimates of βp0Ti*,βp1Ti*,σεp*2 for the pth longitudinal biomarker that need to be averaged.

Model (6) is then used to construct complete longitudinal pseudo data sets which in turn are used to estimate the mean of the posterior distribution of an individual’s random effects given the data. Specifically, multiple complete longitudinal data sets can be constructed by simulating Xpij values from the approximation to the distribution of Xi |Ti given by (6) where the parameters are replaced by their estimated values. Since the simulated data sets have complete follow-up on each individual, the bias in estimating the posterior mean of γi caused by informative dropout will be much reduced.

The posterior mean of distribution γi given the data can be estimated by fitting (2)–(3) to the generated complete longitudinal pseudo data. However, similar to fitting the conditional model (6), fitting model (2)–(3) is difficult due to the high dimension of Σγ. Thus, we again use the pairwise estimation approach of Fieuws and Verbeke (2005), whereby we estimate the parameters of (2)–(3) by fitting all pairwise bivariate models and averaging duplicate parameter estimates to obtain final parameter estimates. For each generated complete longitudinal pseudo data set, the estimate of the posterior mean, denoted as [gamma with circumflex]i = ([gamma with circumflex]1i0, [gamma with circumflex]1i1, …, [gamma with circumflex]Pi0, [gamma with circumflex]Pi1)′, can be calculated as


where Zi is a PJ × 2P design matrix corresponding to the fixed and random effects in (2)–(3), where Zi=diag(A,A,,A)Ptimes,A=(11(...)1t1t2(...)tJ), and Vi is the variance of Xi. Estimates of Xpij*, denoted as X^pij*, are obtained by substituting ([beta]p0, [beta]p1, [gamma with circumflex]pi0, [gamma with circumflex]pi1) for (βp0, βp1, γpi0, γpi1) in (3).

To account for the measurement error in using [gamma with circumflex]i as compared with using γi in (1), we note that

P(Yij=1[mid ]Yi(j1)=0;X^i*)=Φ(α0j+p=1PαpX^pij*1+Var{p=1Pαp(X^pij*Xpij*)}),

where Var{p=1Pωp(X^pi(j1)*Xpi(j1)*)}=RijVar(γ^iγi)Rij, Rij = (ω1, ω1tj−1, ω2, ω2tj−1, …, ωp, ωptj−1), Var(γ^iγi)=γγZi{Vi1ZiVi1×ZiQZiVi1}Ziγ, and where Q=(i=1IZiVi1Zi)1 [Laird and Ware (1982); Verbeke and Molenberghs (2000)]. Expression (8) follows from the fact that E[Φ(a+V)]=Φ[(a+μ)/1+τ2], where V ~ N (μ, τ2).

In the second stage, α0j (j = 1, 2, …, J) and αp (p = 1, 2, …, P) can be estimated by maximizing the likelihood

L=[product]i=1I[[product]j=2Ji{1P(Yij=1[mid ]Yi(j1)=0;X^i*)}]×P(Yi(Ji+1)=1[mid ]YiJi=0;X^i*)Ji<J,

where P(Yij=1[mid ]Yi(j1)=0,X^i*) is given by (8). Thus, we propose the following algorithm for estimating α0j and αp (p = 1, 2, …, P):

  1. Estimate the parameters of model (6) by fitting (P2) bivariate models to each of the pairwise combinations of longitudinal measurements and averaging duplicate parameter estimates. The bivariate models can be fit in R [Venables, Smith and the R Development Core Team (2008)] using code presented in Doran and Lockwood (2006).
  2. Simulate complete longitudinal pseudo measurements (i.e., Xpij for p = 1, 2, …, P, i = 1, 2, …, I, j = 1, 2, …, J) from model (6) with model parameters estimated from step 1.
  3. Estimate the parameters in model (2)–(3) without regard to the event time distribution from complete longitudinal pseudo measurements (simulated in step 2) by fitting all possible (P2) bivariate longitudinal models and averaging duplicate model parameter estimates.
  4. Calculate [gamma with circumflex]i using (7) and X^pij* using (3) with γi replaced by [gamma with circumflex]i and β being replaced by [beta] estimated in step 3.
  5. Estimate α0j (j = 2, 3, …, J) and αp (p = 1, 2, …, P) using (8) and (9).
  6. Repeat steps 2 to 5 M times and average [alpha]0j and [alpha]p to get final estimates.

We choose M = 10 in the simulations and data analysis since this was shown to be sufficiently large for univariate longitudinal modeling discussed in Albert and Shih (2009). Asymptotic standard errors of [alpha]0j and [alpha]p cannot be used for inference since they fail to account for the missing data uncertainty in our procedure. Standard errors and 95% confidence intervals of parameter estimates using the bootstrap [Efron and Tibshirani (1993)] are as follows:

  1. Construct a bootstrap sample of size I, by resampling event-time and multivariate longitudinal data with replacement (( Tib,X1ib,X2ib,,Xpib)) from the (Ti, X1i, X2i, …, Xpi).
  2. Fit the proposed estimation procedure.
  3. Iterate 500 times between steps 1 and 2. The bootstrap standard error is the sample standard deviation of the 500 bootstrap estimates. The 95% confidence intervals were constructed using the percetile method (limits are 2.5 and 97.5 percentiles of the bootstrap distribution).

2.3. Incorporating covariate dependence

Covariates can be incorporated in (3) by adding them directly into the multivariate linear mixed model (6). Specifically, if Xpij*=Ziηp+βp0+βp1tj+γpi0+γpi1tj, where Zi is a vector of covariates with ηp being parameters for the pth biomarker, then P (Xi|Ti, Zi) = ∫ h(Xi|γi, Zi)g(γi|Ti) dγi and Xi|Ti can be approximated by a multivariate linear mixed model with Zi ηp being added to the right side of (6). Estimation then proceeds as described in Section 2.2 Although more difficult, covariates can also be incorporated into (1). If

P(Yij=1[mid ]Yi(j1)=0,Xi*,Zi)=Φ(α0j+Ziζ+p=1PαpXpi(j1)*),

then P (Xi|Ti, Zi ) = ∫ h(Xi|γi, Zi )g(γi|Ti, Zi ) dγi, and under the assumption that g(γi|Ti, Zi) is normally distributed with variance not depending on Zi (which we found to be the case in simulations not shown), then Xi|Ti, Zi can be approximated by a multivariate linear mixed model with ZiζpTi* added to the right-hand side of (6). Extensive simulations showed that the conditional distribution of γi|Ti, Zi is nearly normally distributed over a wide range of parameter values, with slight departures from normality found when the relationship between the longitudinal process and event-time processes is very strong (i.e., αp’s are very large in magnitude). The multivariate linear model approximation is flexible in that it allows the regression parameters to vary with Ti. A more parsimonious model would be to constrain the parameters such that ζpTi* does not vary with Ti.

3. Simulations

We conduct a simulation study to examine the statistical properties of the proposed approach. The approach is examined for the situation where P = 3, I = 300, and σpε2=0.75 for p = 1, 2 and 3. The remaining parameters are presented in Table 1. We compare the proposed approach with M = 10 with a model where the Xpij* ’s are assumed to be known (true model) and with a model where the observed Xpij’s are used in (1) (observed model). Although the true values are never actually observed in practice, we examine the true model as a benchmark in comparing the other models. Table 1 shows that the proposed approach results in nearly unbiased estimates of α1, α2 and α3, whereas the model which uses the observed observations (which are subject to measurement error) has severe bias for estimating α1 and α3 (the two parameters which are nonzero). Although the variability of parameter estimates is larger for the are presented in proposed approach as compared with the observed approach, the root mean squared errors are substantially smaller for α1 and α3 for the proposed approach. For example, the root mean squared errors for α1 is 0.089 for the proposed approach and 0.184 for the observed model. Results when the “true” values for the markers meters Table 1. (markers without measurement errors) are assumed known are also presented in Table 1 (column labeled Truth). As expected, estimates are unbiased for this gold standard case. A comparison of the gold standard case with the proposed approach shows efficiency loss. For example, the relative efficiency for estimating α1, α2 and α3 with the proposed approach versus the gold standard is 0.45 [(0.06/0.089)2], 0.57 and 0.24, respectively. We conducted additional simulations with different parameter values. In all cases tried, the mean squared errors were substantially smaller using the proposed approach as compared with using the observed values (data not shown).

Table 1
Estimates of α0j = α0, α1, α2 and α3 from model (1) with σ1e = σ2e = σ3e = 0.75 and with P = 3, J = 5, and I = 300. Random effects are generated under a diagonal covariance matrix. Further, ...

Table 1 also presents the average estimated intercept and slope for each of the three longitudinal biomarkers. The results show that estimates of these fixed effects are nearly unbiased for the proposed approach.

4. Example

We examine the effect of multiple longitudinal biomarkers on the short-term prognosis for patients with primary biliary cirrhosis (PBC) using the PBC study conducted at the Mayo Clinic from 1974 to 1984 [Murtaugh et al. (1994)]. PBC is a chronic disease characterized by inflammatory destruction of the small bile ducts within the liver, which eventually leads to cirrhosis of the liver, followed by death. Various biomarkers such as biliribin, prothrombin time and albumin were collected longitudinally, and interest is on examining whether these biomarkers relate to the natural history of disease. Of major interest was whether these biomarkers are prognostic for transplantation-free survival (time to either transplantation or death). A total of 312 patients had a baseline measurement and were followed longitudinally at 6 months and at yearly intervals thereafter.

For our application, we focused on the first 4 years of follow-up for a number of reasons. First, individual changes in the biomarkers appeared to be close to linear over this time period. Figure 1 shows 4 examples of complete follow-up in which the changes in log-transformed biliribin appear to be linear over the first 4 years of follow-up and not very linear over the whole range of follow-up. For each panel, the solid line is a least-squares regression line for the first four years of follow-up, while the dashed line is a corresponding line using all the follow-up data. The patterns over the whole follow-up period are nonlinear curves and are not systematic over subjects, and therefore not easily characterized by a simple nonlinear mixed model. The reason why the linear assumption is reasonable over the shorter time interval is that, even though the curves are nonlinear, they can adequately be approximated as linear functions over a short time interval (i.e., a nonlinear function can be locally approximated by a linear function). Second, the methodology makes the assumption that the effect of the biomarkers on prognosis is constant over the follow-up period [i.e., αp parameters in (1) do not vary over time]. This assumption is more reasonable over the shorter 4 year interval rather than the entire follow-up period.

FIG. 1
Plot of log-transformed biliribin values versus follow-up time for 4 patients. Each plot shows an example of complete follow-up in which the changes over time appear to be linear over the first 4 years of follow-up and nonlinear over the entire length ...

Table 2 shows parameter estimates from fitting model (1) with the observed data as covariates instead of the true values. Although both standard errors and 95% confidence intervals were estimated using the bootstrap, only the 95% confidence intervals are presented since the bootstrap estimates were not normally distributed for many of the parameter estimates. The results demonstrate a statistically significant positive effect of biliribin and prothrombin time and a negative effect of albumin on transplantation-free survival. However, it should be recognized that these parameter estimates may be distorted due to the measurement error in these longitudinal biomarker measurements.

Table 2
The effect of log-transformed biliribin, prothrombin time and albumin on transplantation-free survival. The analysis is based on fitting model (1) with Xi(j − 1) replacing Xi(j1)*. Time-to-event is modeled as a discrete time ...

Using the proposed approach, we initially fit model (1)–(3) which incorporated a random intercept and slope term for each of the three biomarkers. However, the random effect for slope for prothrombin time and albumin were estimated as nearly zero. Thus, we re-fit the model without a random slope effect for these two biomarkers. Table 3 shows parameter estimates from the proposed approach with 95% confidence intervals estimated using the bootstrap (as in the analysis with the observed biomarkers presented in Table 2, we do not present the parameter estimates of the standard errors). Except for the effect of biliribin, estimates of the other two biomarkers are substantially larger in magnitude under the proposed approach than when ignoring measurement error and using the observed data (Table 2). This is consistent with the common phenomenon that ignoring measurement error attenuates parameter estimates. In terms of inference, the effect of prothrombin time on short-term prognosis is no longer statistically significant with the proposed approach, while the effects of biliribin and albumin on prognosis are statistically significant with both approaches. In the PBC analysis, estimates of σpε2 were 0.31, 0.12 and 0.11 for log-transformed values of biliribin, albumin and prothrombin time, respectively. The smaller absolute values for parameter estimates of albumin and prothrombin time using the observed markers (Table 1) relative to estimates for these markers using the proposed approach (Table 2) can be attributed to attenuation due to measurement error, since, in these cases, the residual variances are substantially larger than the between-subject variations.

Table 3
The effect of log-transformed biliribin, prothrombin time and albumin (p = 1, 2 and 3, respectively) on transplantation-free survival using the proposed approach with M = 10. We fit model (1) with the true longitudinal measurements following (2)–( ...

We also conducted analyses where we adjusted for treatment effect and age in the discrete-time survival model [Zi is treatment group or age in model (10)]. As discussed in Section 2.3, we constrained the parameters ζpTi* so that they did not vary with Ti (results were similar when we did not constrain the parameters). When we adjusted for treatment group (with treatment group coded as 1 for D-penicillamine and 0 for placebo) in (10), we estimated the α coefficient corresponding to treatment as 0.051 (95% CI: −0.84 to 1.14). The estimates of other parameters were almost identical to those presented in Table 3. When we adjusted for age in (10), we estimated the α coefficient corresponding to age as 0.018 (95% CI: 0.01 to 0.129), where age was scaled in units of a year. Although age was statistically significant, the effects of log bilirubin, albumin and prothtime on transplantation-free survival were similar to those for the unadjusted model. Specifically, the regression coefficients (α coefficients) corresponding to three markers are 0.394 (95% CI: 0.07 to 2.21), −10.65 (−83.51 to −5.52) and 5.76 (−4.28 to 30.97).

Table 3 also shows the estimated fixed effect intercept and slope for the three longitudinal biomarkers for the proposed approach. The estimates suggest that biliribin is increasing, while albumin and prothrombin time are nearly constant over time.

The joint modeling approach is important in this application for a number of reasons. First, survival models which use observed biomarkers can result in attenuated estimates of risk. The proposed approach allowed us to account for the measurement error in investigating the effect of multiple biomarker measurements on the short-term prognosis of PBC patients in terms of transplantation-free survival. With the proposed approach, we found that the “true” biliribin and albumin values at the beginning of an interval had a sizable and statistically significant effect on the probability of either a transplantation or death in the subsequent interval. Second, the proposed approach allows us to appropriately make inference about changes in the three “true” biomarkers over time. As stated before, the largest change over time was in biliribin which sizably increased over time. When making these longitudinal inferences, not appropriately modeling the relationship between the multiple biomarkers and survival may lead to bias due to informative dropout [Wu and Carroll (1988)].

5. Discussion

We proposed an approach for jointly modeling multivariate longitudinal and discrete time-to-event data. Unlike likelihood-based approaches which require high-dimensional integration to evaluate the joint likelihood, this approach only requires fitting bivariate random effects models. This methodology uses recent methodology for fitting multivariate longitudinal data with bivariate linear mixed models proposed by Fieuws et al. (2005, 2007). They discussed the simple averaging of duplicate parameters estimates as we did in implementing the proposed approach. They also proposed a pseudo-likelihood approach which involves maximizing the sum of likelihoods from bivariate models across all (P2) combinations of pairwise longitudinal markers. Although this later approach may provide some minor efficiency gain over simple averaging, it would be substantially more complicated to implement in our setting. Further, one of the advantages of the pseudo-likelihood approach is that it provides an analytic expression for the asymptotic variance of the parameter estimates. Unfortunately, this asymptotic variance is not generalizable to the joint model with time-to-event data, making the pseudo-likelihood approach less attractive in our setting.

There are similarities between our approach and the recent approach by Fieuws et al. (2008) for predicting renal graft failure based on multivariate longitudinal profiles. Both approaches model the conditional distribution of the multivariate longitudinal profiles given the failure time. However, there are major differences between the two approaches. Fieuws et al. model the conditional distribution of the longitudinal measurements given the failure time and then use Bayes rule to estimate the probability of failure given the longitudinal profiles. In our approach, we use an approximation of the conditional distribution of Xi|Ti under the joint model of the multivariate longitudinal and time-to-event data in order estimate the parameters of this joint model.

We demonstrated the feasibility of the proposed approach with three biomarkers (P = 3). However, the approach can easily accommodate a large number of longitudinal profiles since it simply involves fitting (P2) bivariate models. The relationship between the multivariate longitudinal and event-time data is governed by expression (1). However, other functional relationships are possible with this approach. For example, we could relate the two processes by averages of “true” longitudinal biomarkers either across time or across different biomarkers. Alternatively, the approach could be formulated so that the event-time process depends on the individual’s intercept and slope for each of the P longitudinal biomarkers.

The multivariate longitudinal profiles are modeled as multivariate linear mixed models in (2) and (3), which was appropriate for the analysis of the PBC data. However, the methodology could be extended to allow for more flexible nonlinear modeling of marker profiles. This would involve approximating the conditional probability Xi|Ti in (5) where h(Xi|γi) follows a multivariate nonlinear mixed model, rather than the linear mixed model discussed in our paper. In the nonlinear case, we could approximate (5) by (6), where (6) would be a nonlinear mixed model with parameters indexed by Ti rather than the linear mixed model presented. However, unless there is biological rational for a particular nonlinear mixed model, it may be difficult to choose a reasonable model in most practical situations.

In our formulation, we assumed that event times are only administratively censored after a fixed follow-up at the end of the study. For the situation in which patients are censored prematurely, dropout times can be imputed based on a model fit using patients who had the potential to be followed over the entire study duration.

In this article methodology was developed using a discrete-time survival model with calibration error being incorporated by using a probit link function. This approach led to an analytically tractable form for incorporating calibration error (8). For the PBC study, little is lost by using a discrete-time model since the probability of an event in each of the five intervals is low (the estimated probability of an event during each of the five time intervals is 0.05, 0.03, 0.08, 0.07 and 0.07). For continuous-time survival models such as the Cox model, incorporating calibration is more difficult since there is no simple analytic solution. That said, we could use a Cox model if we do not account for the calibration error in replacing the random effect by their empirical Bayes estimators. Although not the case in the PBC study, in situations where the within-subject variation is small relative to the between-subject sources of variation, the calibration error will be small and there will be only a small amount of bias induced by not accounting for the calibration error.


We thank the Center for Information Technology, NIH, for providing access to the high-performance computational capabilities of the Biowulf cluster computer system. We thank the Editor, Associate Editor and reviewer for their constructive comments which led to an improved manuscript.


1Supported by the Intramural Research Program of the National Institutes of Health, Eunice Kennedy Shriver National Institute of Child Health and Human Development.


  • Abramowitz M, Stegun I. Handbook of Mathematical Functions. Dover; New York: 1974.
  • Albert PS, Shih JH. On estimating the relationship between longitudinal measurements and time-to-event data using regression calibration. Biometrics. 2009 doi: 10.1111/j.1541-0420.2009.01324.x. [PubMed] [Cross Ref]
  • Allen C, Duffy S, Teknos T, Islam M, Chen Z, Albert PS, Wolf GY, Van Waes C. A prospective study of serial measurements of NF-αβ related serum cytokines as biomarkers of response and survival in patients with advanced oropharyngeal squamous cell carcinoma receiving chemoradiation therapy. Clinical Cancer Research. 2007;13:3182–3190. [PubMed]
  • Chi YY, Ibrahim JG. Joint models for multivariate longitudinal and multivariate survival data. Biometrics. 2006;62:432–445. [PubMed]
  • Brown ER, Ibrahim JG, DeGruttola V. A flexible B-spline model for multiple longitudinal biomarkers and survival. Biometrics. 2005;61:64–73. [PubMed]
  • Doran HC, Lockwood JR. Fitting value-added models in R. Journal of Educational and Behavioral Statistics. 2006;31:205–230.
  • Efron B, Tibshirani RJ. An Introduction to the Boostrap. Chapman and Hall; New York: 1993.
  • Fieuws S, Verbeke G. Pairwise fitting of mixed models for the joint modelling of multivariate longitudinal profiles. Biometrics. 2005;62:424–431. [PubMed]
  • Fieuws S, Verbeke G, Molenberghs G. Random-effects models for multivariate repeated measures. Stat Methods Med Res. 2007;16:387–397. [PubMed]
  • Fieuws S, Verbeke G, Maes B, Vanrenterghem Y. Predicting renal graft failure using multivariate longitudinal profiles. Biostatistics. 2008;9:419–431. [PubMed]
  • Henderson R, Diggle P, Dobson A. Joint modeling of measurements and event time data. Biostatistics. 2000;1:465–480. [PubMed]
  • Huang WH, Zeger SL, Anthony JC, Garrett E. Latent variable model for joint anlaysis of multiple repeated measures and bivariate event times. Amer Statist Assoc. 2001;96:906–914.
  • Ibrahim JG, Chen M, Sinha D. Bayesian methods for jointly modeling of longitudinal and survival data with applications to cancer vaccine trials. Statist Sinica. 2004;14:863–883.
  • Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38:963–974. [PubMed]
  • Murtaugh PA, Dickson ER, Van Dam GM, Malincho M, Grambsch PM, Langworthy AL, Gips CH. Primary billary cirrhosis: Prediction of short-term survival based on repeated patient visits. Hepatology. 1994;20:126–134. [PubMed]
  • Song X, Davidian M, Tsiatis AS. An estimator for the proportional hazards model with multiple longitudinal covariates measured with error. Biostatistics. 2002;3:511–524. [PubMed]
  • Tsiatis AA, Davidian M. Joint modeling of longitudinal and time-to-event data: An overview. Statist Sinica. 2004;14:809–834.
  • Tsiatis AA, DeGruttola V, Wulfsohn MS. Modeling the relationship of survival to longitudinal data measured with error. Applications to survival and CD4 counts in patients with AIDS. J Amer Statist Assoc. 1995;90:27–37.
  • Venables WN, Smith DM. the R Development Core Team. An Introduction to R. Version 2.8.1 (2008-12-22) 2008.
  • Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. Springer; New York: 2000.
  • Wei GCG, Tanner MA. A Monte-Carlo implementation of the E–M algorithm and the poor man’s data augmentation algorithm. J Amer Statist Assoc. 1990;85:699–704.
  • Wu MC, Carroll RJ. Estimation and comparison of changes in the presence of informative right censoring by modeling the censoring process. Biometrics. 1988;45:939–955. [PubMed]
  • Wulfsohn MS, Tsiatis AA. A joint model for survival and longitudinal data measured with error. Biometrics. 1997;53:330–339. [PubMed]
  • Xu J, Zeger SL. Joint analysis of longitudinal data comprising repeated measures and times to events. Appl Statist. 2001a;50:375–387.
  • Xu J, Zeger SL. The evaluation of multiple surrogate endpoints. Biometrics. 2001b;57:81–87. [PubMed]