An approach for jointly modeling multivariate longitudinal measurements and discrete time-to-event data
aa r X i v : . [ s t a t . A P ] N ov The Annals of Applied Statistics (cid:13)
Institute of Mathematical Statistics, 2010
AN APPROACH FOR JOINTLY MODELING MULTIVARIATELONGITUDINAL MEASUREMENTS AND DISCRETETIME-TO-EVENT DATA By Paul S. Albert and Joanna H. Shih
Eunice Kennedy Shriver National Institute of Child Health and HumanDevelopment and National Cancer Institute
In many medical studies, patients are followed longitudinally andinterest is on assessing the relationship between longitudinal measure-ments and time to an event. Recently, various authors have proposedjoint modeling approaches for longitudinal and time-to-event datafor a single longitudinal variable. These joint modeling approachesbecome intractable with even a few longitudinal variables. In this pa-per we propose a regression calibration approach for jointly modelingmultiple longitudinal measurements and discrete time-to-event data.Ideally, a two-stage modeling approach could be applied in which themultiple longitudinal measurements are modeled in the first stage andthe longitudinal model is related to the time-to-event data in the sec-ond stage. Biased parameter estimation due to informative dropoutmakes this direct two-stage modeling approach problematic. We pro-pose a regression calibration approach which appropriately accountsfor informative dropout. We approximate the conditional distribu-tion of the multiple longitudinal measurements given the event timeby modeling all pairwise combinations of the longitudinal measure-ments using a bivariate linear mixed model which conditions on theevent time. Complete data are then simulated based on estimatesfrom these pairwise conditional models, and regression calibration isused to estimate the relationship between longitudinal data and time-to-event data using the complete data. We show that this approachperforms well in estimating the relationship between multivariate lon-gitudinal measurements and the time-to-event data and in estimatingthe parameters of the multiple longitudinal process subject to infor-mative dropout. We illustrate this methodology with simulations andwith an analysis of primary biliary cirrhosis (PBC) data.Received August 2009; revised February 2010. Supported by the Intramural Research Program of the National Institutes of Health,
Eunice Kennedy Shriver
National Institute of Child Health and Human Development.
Key words and phrases.
Joint models, shared random parameter models, informativedropout, regression calibration.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Applied Statistics ,2010, Vol. 4, No. 3, 1517–1532. This reprint differs from the original in paginationand typographic detail. 1
P. S. ALBERT AND J. H. SHIH
1. Introduction.
Recently, many studies collect longitudinal data on apanel of biomarkers, and interest is on assessing the relationship betweenthese biomarkers and time to an event. For example, Allen et al. (2007)examined the relationship between five longitudinally collected cytokinesmeasured from serum plasma and survival. Interest focused on whether thevalues of these multiple cytokines are associated with survival. In anotherexample, patients with primary biliary cirrhosis are followed longitudinallyand interest is on examining whether multiple longitudinally biomarkers areprognostic for a poor clinical outcome. Important features in studies of thistype are that there may be a relatively large number of biomarkers and thatthese biomarkers are subject to sizable measurement error due to laboratoryerror and biological variation.Various authors have proposed joint modeling approaches for a singlelongitudinal measurement and time-to-event data [Tsiatis, DeGruttola andWulfsohn (1995); Wulfsohn and Tsiatis (1997); Tsiatis and Davidian (2004);Henderson, Diggle and Dobson (2000), among others]. There is also limitedwork on joint models for a few longitudinal measurements and time-to-eventdata [Xu and Zeger (2001a, 2001b); Huang et al. (2001); Song, Davidianand Tsiatis (2002); Ibrahim, Chen and Sinha (2004); Brown, Ibrahim andDeGruttola (2005); Chi and Ibrahim (2006)]. However, these methods aredifficult to implement when the number of longitudinal biomarkers is largesince most of these approaches require integrating over the vector of all ran-dom effects to evaluate the joint likelihood of the multivariate longitudinaland time-to-event data. This paper proposes an approach for jointly mod-eling multivariate longitudinal and discrete time-to-event data which easilyaccommodates many longitudinal biomarkers.Fieuws and Verbeke (2005) and Fieuws, Verbeke and Molenberghs (2007)have proposed an approach for modeling multivariate longitudinal datawhereby all possible pairs of longitudinal data are separately modeled andare then combined in a final step. We use a similar approach along with arecent regression calibration approach for jointly modeling a single series oflongitudinal 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 approachfor using multivariate longitudinal profiles to predict renal graft failure. Atthe end of their discussion, they mention that a more elegant approach,which has not yet been developed, would involve a joint model for the manylongitudinal profiles and time-to-event data. This paper presents such anapproach.We describe the approach in Section 2. We show the advantages of thisapproach using simulation in Section 3. We illustrate the methodology withan analysis of primary biliary cirrhosis data (PBC) in which we simultane-ously examine the relationship between multiple longitudinal biomarker inSection 4. A discussion follows in Section 5.
ULTIVARIATE LONGITUDINAL AND TIME-TO-EVENT DATA
2. Modeling approach.
Define T i to be a discrete event-time which cantake on discrete values t j , j = 1 , , . . . , J , and Y ij to be a binary indictor ofwhether the i th patient is dead at time t j . Then J i = P Jj =1 (1 − Y ij ) = J − Y i · ,where Y i · = P Jj =1 Y ij indicates the number of follow-up measurements beforethe event or the end of follow-up at time t J . Longitudinal measurements aremeasured at times t , t , . . . , t J i . Denote X i = ( X i , X i , . . . , X iJ i ) ′ , X i =( X i , X i , . . . , X iJ i ) ′ , . . . , X P i = ( X P i , X P i , . . . , X P iJ i ) ′ as the P biomark-ers measured repeatedly at j = 1 , , . . . , J i time points. Further, define X ∗ pi =( X ∗ pi , X ∗ pi , . . . , X ∗ piJ i ) ′ as the longitudinal measurements without measure-ment error for the p th biomarker and X ∗ i = ( X ∗ i , X ∗ i , . . . , X ∗ P i ). We considera joint model for multivariate longitudinal and discrete time-to-event datain which the discrete event time distribution is modeled as a linear functionof previous true values of the biomarkers without measurement error on theprobit scale. Specifically, P ( Y ij = 1 | Y i ( j − = 0; X ∗ i ) = Φ α j + P X p =1 α p X ∗ pi ( j − ! , (1)where i = 1 , , . . . , I , j = 2 , , . . . , J i , Y i is taken as 0, α j governs the base-line discrete event time distribution and α p measures the effect of the p thbiomarker ( p = 1 , , . . . , P ) at time t j − on survival at time t j . Specifically,(1) allows for examining the effect of multiple “true” biomarker values attime j − j − j th timepoint.The longitudinal data is modeled assuming that the fixed and random ef-fect trajectories are linear. Specifically, the multivariate longitudinal biomark-ers can be modeled as X pij = X ∗ pij + ε pij , (2)where X ∗ pij = β p + β p t j + γ pi + γ pi t j , (3)where β p and β p are the fixed effect intercept and slope for the p thbiomarker, and γ pi and γ pi are the random effect intercept and slope for the p th biomarker on the i th individual. Denote β = ( β , β , β , β , . . . , β P ,β P ) ′ and γ i = ( γ i , γ i , γ i , γ i , . . . , γ P i , γ P i ) ′ . We assume that γ i isnormally distributed with mean and variance Σ γ , where Σ γ is a 2 P by 2 P dimensional variance matrix, and ε pij are independent error termswhich are assumed to be normally distributed with mean 0 and variance σ pε ( p = 1 , , . . . , P ).Alternative to (1), where the probability of an event over an intervaldepends on the true biomarker values at the beginning of the interval, theevent-time process could be adapted to depend on the random effects of themultivariate longitudinal process [e.g., γ pi can replace X ∗ pi ( j − in (1)]. P. S. ALBERT AND J. H. SHIH
Difficulty in joint estimation.
Conceptually, model (1)–(2) can beestimated by maximizing the likelihood L = I Y i =1 Z γ i · · · Z ( P Y p =1 h ( X pi | γ pi , γ pi ) ) (4) × ( J i Y j =2 (1 − r ij ) ) ( r i ( J i +1) ) J i
Estimation.
We propose a two stage regression calibration approachfor estimation, which can be described as follows. In the first stage, multi-variate linear mixed models can be used to model the longitudinal data. Inthe second stage, the time-to-event model is estimated by replacing the ran-dom effects with corresponding empirical Bayes estimates. There are threeproblems with directly applying this approach. First, estimation in the firststage is complicated by the fact that simply fitting multivariate linear mixedmodels results in bias due to informative dropout; this is demonstrated byAlbert and Shih (2009) for the the case of P = 1. Second, as described inSection 2.1, parameter estimation for multivariate linear mixed models canbe computationally difficult when the number of longitudinal measurements( P ) is even moderately large. Third, calibration error in the empirical Bayesestimation needs to be accounted for in the time-to-event model. The pro-posed approach will deal with all three of these problems.The bias from informative dropout is a result of differential follow-upwhereby the longitudinal process is related to the length of follow-up. Thatis, in (2)–(3), patients with large values of X ∗ pij are more likely to havean early event when α p > p = 1 , , . . . , P . There would be no bias ifall J follow-up measurements were observed on all patients. As proposedby Albert and Shih (2009) for univariate longitudinal data, we can avoidthis bias by generating complete data from the conditional distribution of ULTIVARIATE LONGITUDINAL AND TIME-TO-EVENT DATA X i = ( X i , X i , . . . , X P i ) given T i , denoted as X i | T i . Since X i | T i under model(2)–(3) does not have a tractable form, we propose a simple approximationfor this conditional distribution. Under model (2)–(3), the distribution of X i | T i can be expressed as P ( X i | T i ) = Z h ( X i | γ i , T i ) g ( γ i | T i ) d γ i . (5)Since T i and the values of X i are conditional independent given γ i , h ( X i | γ i , T i ) = h ( X i | γ i ), where h ( X i | γ i ) = Q Pp =1 h ( X pi | γ pi , γ pi ). The distri-bution of X i | T i can be expressed as a multivariate linear mixed model if weapproximate g ( γ i | T i ) by a normal distribution. Under the assumption that g ( γ i | T i ) is normally distributed with mean µ T i = ( µ T i , µ T i , µ T i , µ T i ,. . . , µ P T i , µ P T i ) ′ and variance Σ ∗ γ T i , and by rearranging mean structure pa-rameters in the integrand of (5) so that the random effects have mean zero, X i | T i corresponds to the following multivariate linear mixed model: X pij | ( T i , γ ∗ ip T i , γ ∗ ip T i ) = β ∗ p T i + β ∗ p T i t j + γ ∗ ip T i + γ ∗ ip T i t j + ε ∗ pij , (6)where i = 1 , , . . . , I , j = 1 , , . . . , J i , and p = 1 , , . . . , P . The parameters β ∗ p T i and β ∗ p T i are intercept and slope parameters for the p th longitudi-nal measurement and for patients who have an event time at time T i orwho are censored at time t J . In addition, the associated random effects γ ∗ iT i = ( γ ∗ i T i , γ ∗ i T i , γ ∗ i T i , γ ∗ i T i , . . . , γ ∗ iP T i , γ ∗ iP T i ) ′ are multivariate normalwith mean and variance Σ ∗ γ T i , and the residuals ε ∗ pij are assumed to havean independent normal distribution with mean zero and variance σ ∗ εp . Thus,this conditional model involves estimating separate fixed effect interceptand slope parameters for each potential event-time and for subjects who arecensored at time t J . Likewise, separate random effects distributions are esti-mated for each of these discrete time points. For example, the intercept andslope fixed-effect parameters for the p th biomarker for those patients whohave an event at time T i = t is β ∗ p t and β ∗ p t , respectively. Further, theintercept and slope random effects for all P biomarkers on those patientswho have an event at time T i = t , γ ∗ it , is multivariate normal with mean and variance Σ ∗ γ t . A similar approximation has been proposed by Albertand Shih (2009) for univariate longitudinal data ( P = 1).Recall that by generating complete data from (6) we are able to avoid thebias due to informative dropout. However, when P is large, direct estimationof model (6) is difficult since the number of elements in Σ ∗ γ T i grows quadrat-ically with P . For example, the dimension of the variance matrix Σ ∗ γ T i is 2 P by 2 P for P longitudinal biomarkers. Fieuws and Verbeke (2005) proposedestimating the parameters of multivariate linear mixed models by formu-lating bivariate linear mixed models on all possible pairwise combinationsof longitudinal measurements. In the simplest approach, they proposed fit-ting bivariate linear mixed models on all (cid:0) P (cid:1) combinations of longitudinal P. S. ALBERT AND J. H. SHIH biomarkers and averaging “overlapping” or duplicate parameter estimates.Thus, we estimate the parameters in the fully specified model (6) by fitting (cid:0) P (cid:1) bivariate longitudinal models that only include pairs of longitudinalmarkers. Fieuws and Verbeke (2005) demonstrated with simulations thatthere is little efficiency loss using their approach relative to a full maximum-likelihood based approach. Fitting these bivariate models is computation-ally feasible since only four correlated random effects are contained in eachmodel. (That is, Σ ∗ γ T i is a four by four-dimensional matrix for each dis-crete event-time T i .) Duplicate estimates of fixed effects and random-effectvariances from all pairwise bivariate models are averaged to obtain final pa-rameter estimates of the fully specified model (6). For example, when P = 4there are ( P −
1) = 3 estimates of β ∗ p T i , β ∗ p T i , σ ∗ εp for the p th longitudinalbiomarker that need to be averaged.Model (6) is then used to construct complete longitudinal pseudo data setswhich in turn are used to estimate the mean of the posterior distribution ofan individual’s random effects given the data. Specifically, multiple completelongitudinal data sets can be constructed by simulating X pij values from theapproximation to the distribution of X i | T i given by (6) where the parametersare replaced by their estimated values. Since the simulated data sets havecomplete follow-up on each individual, the bias in estimating the posteriormean of γ i caused by informative dropout will be much reduced.The posterior mean of distribution γ i given the data can be estimated byfitting (2)–(3) to the generated complete longitudinal pseudo data. However,similar to fitting the conditional model (6), fitting model (2)–(3) is difficultdue to the high dimension of Σ γ . Thus, we again use the pairwise estimationapproach of Fieuws and Verbeke (2005), whereby we estimate the parametersof (2)–(3) by fitting all pairwise bivariate models and averaging duplicateparameter estimates to obtain final parameter estimates. For each generatedcomplete longitudinal pseudo data set, the estimate of the posterior mean,denoted as b γ i = ( b γ i , b γ i , . . . , b γ P i , b γ P i ) ′ , can be calculated as b γ i = Σ γ Z ′ i V − i ( X i − Z i b β ) , (7)where Z i is a P J × P design matrix corresponding to the fixed and randomeffects in (2)–(3), where Z i = diag( A ′ , A ′ , . . . , A ′ ) | {z } P times , A = (cid:16) t t ······ t J (cid:17) , and V i is the variance of X i . Estimates of X ∗ pij , denoted as b X ∗ pij , are obtainedby substituting ( b β p , b β p , b γ pi , b γ pi ) for ( β p , β p , γ pi , γ pi ) in (3).To account for the measurement error in using b γ i as compared with using γ i in (1), we note that P ( Y ij = 1 | Y i ( j − = 0; b X ∗ i ) = Φ (cid:18) α j + P Pp =1 α p b X ∗ pij q { P Pp =1 α p ( b X ∗ pij − X ∗ pij ) } (cid:19) , (8) ULTIVARIATE LONGITUDINAL AND TIME-TO-EVENT DATA where Var { P Pp =1 ω p ( b X ∗ pi ( j − − X ∗ pi ( j − ) } = R ′ ij Var( b γ i − γ i ) R ij , R ij = ( ω ,ω t j − , ω , ω t j − , . . . , ω p , ω p t j − ), Var( b γ i − γ i ) = Σ γ − Σ γ Z ′ i { V − i Z i − V − i × Z i QZ ′ i V − i } Z i Σ γ , and where Q = ( P Ii =1 Z ′ i V − i Z i ) − [Laird and Ware (1982);Verbeke and Molenberghs (2000)]. Expression (8) follows from the fact that E [Φ( a + V )] = Φ[( a + µ ) / √ τ ], where V ∼ N ( µ, τ ).In the second stage, α j ( j = 1 , , . . . , J ) and α p ( p = 1 , , . . . , P ) can beestimated by maximizing the likelihood L = I Y i =1 " J i Y j =2 { − P ( Y ij = 1 | Y i ( j − = 0; b X ∗ i ) } (9) × P ( Y i ( J i +1) = 1 | Y iJ i = 0; b X ∗ i ) J i 2. Fit the proposed estimation procedure.3. Iterate 500 times between steps 1 and 2. The bootstrap standard error isthe sample standard deviation of the 500 bootstrap estimates. The 95%confidence intervals were constructed using the percetile method (limitsare 2.5 and 97.5 percentiles of the bootstrap distribution).2.3. Incorporating covariate dependence. Covariates can be incorporatedin (3) by adding them directly into the multivariate linear mixed model(6). Specifically, if X ∗ pij = Z i η p + β p + β p t j + γ pi + γ pi t j , where Z i is avector of covariates with η p being parameters for the p th biomarker, then P ( X i | T i , Z i ) = R h ( X i | γ i , Z i ) g ( γ i | T i ) d γ i and X i | T i can be approximated bya multivariate linear mixed model with Z i η p being added to the right sideof (6). Estimation then proceeds as described in Section 2.2 Although moredifficult, covariates can also be incorporated into (1). If P ( Y ij = 1 | Y i ( j − = 0 , X ∗ i , Z i ) = Φ α j + Z i ζ + P X p =1 α p X ∗ pi ( j − ! , (10)then P ( X i | T i , Z i ) = R h ( X i | γ i , Z i ) g ( γ i | T i , Z i ) d γ i , and under the assumptionthat g ( γ i | T i , Z i ) is normally distributed with variance not depending on Z i (which we found to be the case in simulations not shown), then X i | T i , Z i canbe approximated by a multivariate linear mixed model with Z i ζ ∗ pT i added tothe right-hand side of (6). Extensive simulations showed that the conditionaldistribution of γ i | T i , Z i is nearly normally distributed over a wide range ofparameter values, with slight departures from normality found when therelationship between the longitudinal process and event-time processes isvery strong (i.e., α p ’s are very large in magnitude). The multivariate linearmodel approximation is flexible in that it allows the regression parametersto vary with T i . A more parsimonious model would be to constrain theparameters such that ζ ∗ pT i does not vary with T i . 3. Simulations. We conduct a simulation study to examine the statis-tical properties of the proposed approach. The approach is examined forthe situation where P = 3, I = 300, and σ pε = 0 . 75 for p = 1 , M = 10 with a model where the X ∗ pij ’s are assumed to beknown (true model) and with a model where the observed X pij ’s are used in(1) (observed model). Although the true values are never actually observedin practice, we examine the true model as a benchmark in comparing theother models. Table 1 shows that the proposed approach results in nearlyunbiased estimates of α , α and α , whereas the model which uses the ob-served observations (which are subject to measurement error) has severe biasfor estimating α and α (the two parameters which are nonzero). Although ULTIVARIATE LONGITUDINAL AND TIME-TO-EVENT DATA Table 1 Estimates of α j = α , α , α and α from model (1) with σ ε = σ ε = σ ε = 0 . and with P = 3 , J = 5 , and I = 300 . Randomeffects are generated under a diagonal covariance matrix. Further, weassume that t j = j and all individuals who are alive at t = 5 areadministratively censored at that time point. The means (standarddeviations) from 500 simulations are presented Parameters True values Truth Proposed Observed α − . − . − . − . . . . α . 40 0 . 408 0 . 405 0 . . . . α . 00 0 . 00 0 . . . . α . 40 0 . 405 0 . 400 0 . . . . β . . . β . . β . . . β − . . β . . β . . σ b . 25 0 . . σ b . 25 0 . . σ b . 25 0 . . σ b . 25 0 . . σ b . 25 0 . . σ b . 25 0 . . the variability of parameter estimates is larger for the proposed approachas compared with the observed approach, the root mean squared errors aresubstantially smaller for α and α for the proposed approach. For exam-ple, the root mean squared errors for α is 0.089 for the proposed approachand 0.184 for the observed model. Results when the “true” values for the P. S. ALBERT AND J. H. SHIH markers (markers without measurement errors) are assumed known are alsopresented in Table 1 (column labeled Truth). As expected, estimates areunbiased for this gold standard case. A comparison of the gold standardcase with the proposed approach shows efficiency loss. For example, the rel-ative efficiency for estimating α , α and α with the proposed approachversus the gold standard is 0.45 [(0 . / . ], 0.57 and 0.24, respectively.We conducted additional simulations with different parameter values. In allcases tried, the mean squared errors were substantially smaller using theproposed approach as compared with using the observed values (data notshown).Table 1 also presents the average estimated intercept and slope for each ofthe three longitudinal biomarkers. The results show that estimates of thesefixed effects are nearly unbiased for the proposed approach. 4. Example. We examine the effect of multiple longitudinal biomark-ers 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 inflam-matory destruction of the small bile ducts within the liver, which eventuallyleads to cirrhosis of the liver, followed by death. Various biomarkers suchas biliribin, prothrombin time and albumin were collected longitudinally,and interest is on examining whether these biomarkers relate to the natu-ral history of disease. Of major interest was whether these biomarkers areprognostic for transplantation-free survival (time to either transplantationor death). A total of 312 patients had a baseline measurement and werefollowed longitudinally at 6 months and at yearly intervals thereafter.For our application, we focused on the first 4 years of follow-up for anumber of reasons. First, individual changes in the biomarkers appearedto be close to linear over this time period. Figure 1 shows 4 examples ofcomplete follow-up in which the changes in log-transformed biliribin appearto be linear over the first 4 years of follow-up and not very linear over thewhole range of follow-up. For each panel, the solid line is a least-squaresregression line for the first four years of follow-up, while the dashed line is acorresponding line using all the follow-up data. The patterns over the wholefollow-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 timeinterval is that, even though the curves are nonlinear, they can adequately beapproximated as linear functions over a short time interval (i.e., a nonlinearfunction can be locally approximated by a linear function). Second, themethodology makes the assumption that the effect of the biomarkers onprognosis is constant over the follow-up period [i.e., α p parameters in (1) do ULTIVARIATE LONGITUDINAL AND TIME-TO-EVENT DATA Fig. 1. Plot of log-transformed biliribin values versus follow-up time for 4 patients. Eachplot shows an example of complete follow-up in which the changes over time appear to belinear over the first 4 years of follow-up and nonlinear over the entire length of follow-up.For each panel, the solid line is a least-squares regression line for the first four years offollow-up, while the dashed line is a corresponding line using all the follow-up data. Table 2 The effect of log-transformed biliribin, prothrombin time and albumin ontransplantation-free survival. The analysis is based on fitting model (1)with X i ( j − replacing X ∗ i ( j − . Time-to-event is modeled as a discretetime process with possible event times at 0.5, 1, 2, 3 and 4 years, where α reflects the baseline discrete-time event distribution for the intervals 0to 0.5 and 0.5 to 1 years. The subsequent yearly intervals are characterizedby α , α , α and α . 95% confidence intervals were estimated usingthe bootstrap with the percentile method (500 bootstrap samples) Parameters Estimate 95% confidence interval α − . − . 10 to − . α − . − . 51 to − . α − . − . 73 to − . α − . − . 91 to − . α − . − . 78 to − . . 58 0 . 45 to 0 . − . − . 76 to − . . 86 0 . 79 to 3 . P. S. ALBERT AND J. H. SHIH Table 3 The effect of log-transformed biliribin, prothrombin time and albumin( p = 1 , 2 and 3, respectively) on transplantation-free survival using theproposed approach with M = 10 . We fit model (1) with the truelongitudinal measurements following (2)–(3). Time-to-event is modeled asa discrete time process with possible event times at 0.5, 1, 2, 3 and 4years. Standard errors (SE) were estimated using the bootstrap. 95%confidence intervals were estimated using the bootstrap with the percentilemethod (500 bootstrap samples) Parameters Estimate 95% confidence interval α − . − . 57 to − . α − . − . 06 to 24 . α − . − . 67 to 27 . α − . − . 88 to 27 . α − . − . 05 to 27 . . 34 0 . 02 to 1 . − . − . 38 to − . . − . 78 to 22 . β . 50 0 . 22 to 0 . β . 26 0 . 09 to 0 . β . 26 0 . 56 to 1 . β − . − . 05 to − . β . 36 1 . 06 to 2 . β . 02 0 . . σ b . 99 0 . 45 to 1 . σ b . 27 0 . 05 to 0 . σ b . 03 0 . 01 to 0 . σ b . 01 0 . 005 to 0 . not vary over time]. This assumption is more reasonable over the shorter 4year interval rather than the entire follow-up period.Table 2 shows parameter estimates from fitting model (1) with the ob-served data as covariates instead of the true values. Although both standarderrors and 95% confidence intervals were estimated using the bootstrap,only the 95% confidence intervals are presented since the bootstrap esti-mates were not normally distributed for many of the parameter estimates.The results demonstrate a statistically significant positive effect of biliribinand prothrombin time and a negative effect of albumin on transplantation-free survival. However, it should be recognized that these parameter esti-mates may be distorted due to the measurement error in these longitudinalbiomarker measurements.Using the proposed approach, we initially fit model (1)–(3) which incor-porated a random intercept and slope term for each of the three biomarkers.However, the random effect for slope for prothrombin time and albumin wereestimated as nearly zero. Thus, we re-fit the model without a random slope ULTIVARIATE LONGITUDINAL AND TIME-TO-EVENT DATA effect for these two biomarkers. Table 3 shows parameter estimates from theproposed approach with 95% confidence intervals estimated using the boot-strap (as in the analysis with the observed biomarkers presented in Table 2,we do not present the parameter estimates of the standard errors). Exceptfor the effect of biliribin, estimates of the other two biomarkers are substan-tially larger in magnitude under the proposed approach than when ignoringmeasurement error and using the observed data (Table 2). This is consistentwith the common phenomenon that ignoring measurement error attenuatesparameter estimates. In terms of inference, the effect of prothrombin time onshort-term prognosis is no longer statistically significant with the proposedapproach, while the effects of biliribin and albumin on prognosis are sta-tistically significant with both approaches. In the PBC analysis, estimatesof σ pε were 0.31, 0.12 and 0.11 for log-transformed values of biliribin, al-bumin and prothrombin time, respectively. The smaller absolute values forparameter estimates of albumin and prothrombin time using the observedmarkers (Table 1) relative to estimates for these markers using the proposedapproach (Table 2) can be attributed to attenuation due to measurementerror, since, in these cases, the residual variances are substantially largerthan the between-subject variations.We also conducted analyses where we adjusted for treatment effect andage in the discrete-time survival model [ Z i is treatment group or age inmodel (10)]. As discussed in Section 2.3, we constrained the parameters ζ ∗ pT i so that they did not vary with T i (results were similar when we did notconstrain the parameters). When we adjusted for treatment group (withtreatment 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: − . 84 to 1.14). The estimates of other parameters were almost identical tothose presented in Table 3. When we adjusted for age in (10), we estimatedthe α coefficient corresponding to age as 0.018 (95% CI: 0.01 to 0.129), whereage was scaled in units of a year. Although age was statistically significant,the effects of log bilirubin, albumin and prothtime on transplantation-freesurvival were similar to those for the unadjusted model. Specifically, theregression coefficients ( α coefficients) corresponding to three markers are0.394 (95% CI: 0.07 to 2.21), − . 65 ( − . 51 to − . 52) and 5.76 ( − . 28 to30.97).Table 3 also shows the estimated fixed effect intercept and slope for thethree longitudinal biomarkers for the proposed approach. The estimates sug-gest that biliribin is increasing, while albumin and prothrombin time arenearly constant over time.The joint modeling approach is important in this application for a num-ber of reasons. First, survival models which use observed biomarkers canresult in attenuated estimates of risk. The proposed approach allowed us P. S. ALBERT AND J. H. SHIH to account for the measurement error in investigating the effect of multi-ple biomarker measurements on the short-term prognosis of PBC patientsin terms of transplantation-free survival. With the proposed approach, wefound that the “true” biliribin and albumin values at the beginning of aninterval had a sizable and statistically significant effect on the probability ofeither a transplantation or death in the subsequent interval. Second, the pro-posed approach allows us to appropriately make inference about changes inthe three “true” biomarkers over time. As stated before, the largest changeover time was in biliribin which sizably increased over time. When makingthese longitudinal inferences, not appropriately modeling the relationshipbetween the multiple biomarkers and survival may lead to bias due to infor-mative dropout [Wu and Carroll (1988)]. 5. Discussion. We proposed an approach for jointly modeling multivari-ate longitudinal and discrete time-to-event data. Unlike likelihood-based ap-proaches which require high-dimensional integration to evaluate the jointlikelihood, this approach only requires fitting bivariate random effects mod-els. This methodology uses recent methodology for fitting multivariate lon-gitudinal data with bivariate linear mixed models proposed by Fieuws etal. (2005, 2007). They discussed the simple averaging of duplicate parame-ters estimates as we did in implementing the proposed approach. They alsoproposed a pseudo-likelihood approach which involves maximizing the sumof likelihoods from bivariate models across all (cid:0) P (cid:1) combinations of pairwiselongitudinal markers. Although this later approach may provide some minorefficiency gain over simple averaging, it would be substantially more com-plicated to implement in our setting. Further, one of the advantages of thepseudo-likelihood approach is that it provides an analytic expression for theasymptotic variance of the parameter estimates. Unfortunately, this asymp-totic 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 byFieuws et al. (2008) for predicting renal graft failure based on multivariatelongitudinal profiles. Both approaches model the conditional distribution ofthe multivariate longitudinal profiles given the failure time. However, thereare major differences between the two approaches. Fieuws et al. model theconditional distribution of the longitudinal measurements given the failuretime and then use Bayes rule to estimate the probability of failure giventhe longitudinal profiles. In our approach, we use an approximation of theconditional distribution of X i | T i under the joint model of the multivariatelongitudinal and time-to-event data in order estimate the parameters of thisjoint model.We demonstrated the feasibility of the proposed approach with threebiomarkers ( P = 3). However, the approach can easily accommodate a large ULTIVARIATE LONGITUDINAL AND TIME-TO-EVENT DATA number of longitudinal profiles since it simply involves fitting (cid:0) P (cid:1) bivariatemodels. The relationship between the multivariate longitudinal and event-time data is governed by expression (1). However, other functional relation-ships are possible with this approach. For example, we could relate the twoprocesses by averages of “true” longitudinal biomarkers either across timeor across different biomarkers. Alternatively, the approach could be formu-lated so that the event-time process depends on the individual’s interceptand slope for each of the P longitudinal biomarkers.The multivariate longitudinal profiles are modeled as multivariate linearmixed models in (2) and (3), which was appropriate for the analysis of thePBC data. However, the methodology could be extended to allow for moreflexible nonlinear modeling of marker profiles. This would involve approx-imating the conditional probability X i | T i in (5) where h ( X i | γ i ) follows amultivariate nonlinear mixed model, rather than the linear mixed modeldiscussed in our paper. In the nonlinear case, we could approximate (5) by(6), where (6) would be a nonlinear mixed model with parameters indexedby T i rather than the linear mixed model presented. However, unless there isbiological rational for a particular nonlinear mixed model, it may be difficultto choose a reasonable model in most practical situations.In our formulation, we assumed that event times are only administrativelycensored after a fixed follow-up at the end of the study. For the situationin which patients are censored prematurely, dropout times can be imputedbased on a model fit using patients who had the potential to be followedover the entire study duration.In this article methodology was developed using a discrete-time survivalmodel with calibration error being incorporated by using a probit link func-tion. This approach led to an analytically tractable form for incorporatingcalibration 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 intervalsis low (the estimated probability of an event during each of the five timeintervals is 0.05, 0.03, 0.08, 0.07 and 0.07). For continuous-time survivalmodels such as the Cox model, incorporating calibration is more difficultsince there is no simple analytic solution. That said, we could use a Coxmodel if we do not account for the calibration error in replacing the randomeffect by their empirical Bayes estimators. Although not the case in the PBCstudy, in situations where the within-subject variation is small relative tothe between-subject sources of variation, the calibration error will be smalland there will be only a small amount of bias induced by not accounting forthe calibration error. Acknowledgments. We thank the Center for Information Technology,NIH, for providing access to the high-performance computational capabili-ties of the Biowulf cluster computer system. We thank the Editor, Associate P. S. ALBERT AND J. H. SHIH Editor and reviewer for their constructive comments which led to an im-proved manuscript. REFERENCES Abramowitz, M. and Stegun, I. (1974). Handbook of Mathematical Functions . Dover,New York. Albert, P. S. and Shih, J. H. (2009). On estimating the relationship between longi-tudinal measurements and time-to-event data using regression calibration. Biometrics DOI: 10.1111/j.1541-0420.2009.01324.x. Allen, C., Duffy., S., Teknos, T. Islam, M., Chen, Z., Albert, P. S., Wolf, G.Y. and Van Waes, C. (2007). A prospective study of serial measurements of NF- αβ related serum cytokines as biomarkers of response and survival in patients withadvanced oropharyngeal squamous cell carcinoma receiving chemoradiation therapy. Clinical Cancer Research Chi, Y. Y. and Ibrahim, J. G. (2006). Joint models for multivariate longitudinal andmultivariate survival data. Biometrics Brown, E. R., Ibrahim, J. G. and DeGruttola, V. (2005). A flexible B-spline modelfor multiple longitudinal biomarkers and survival. Biometrics Doran, H. C. and Lockwood, J. R. (2006). Fitting value-added models in R. Journalof Educational and Behavioral Statistics Efron, B. and Tibshirani, R. J. (1993). An Introduction to the Boostrap . Chapman andHall, New York. MR1270903 Fieuws, S. and Verbeke, G. (2005). Pairwise fitting of mixed models for the joint mod-elling of multivariate longitudinal profiles. Biometrics Fieuws, S., Verbeke, G. and Molenberghs, G. (2007). Random-effects models formultivariate repeated measures. Stat. Methods Med. Res. Fieuws, S., Verbeke, G., Maes, B. and Vanrenterghem, Y. (2008). Predicting renalgraft failure using multivariate longitudinal profiles. Biostatistics Henderson, R. Diggle, P. and Dobson, A. (2000). Joint modeling of measurementsand event time data. Biostatistics Huang, W. H., Zeger, S. L., Anthony, J. C. and Garrett, E. (2001). Latent variablemodel for joint anlaysis of multiple repeated measures and bivariate event times. Amer.Statist. Assoc. Ibrahim, J. G., Chen, M. and Sinha, D. (2004). Bayesian methods for jointly modelingof longitudinal and survival data with applications to cancer vaccine trials. Statist.Sinica Laird, N. M. and Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics Murtaugh, P. A., Dickson, E. R., Van Dam, G. M., Malincho, M., Grambsch, P.M., Langworthy, A. L. and Gips, C. H. (1994). Primary billary cirrhosis: Predictionof short-term survival based on repeated patient visits. Hepatology Song, X., Davidian, M. and Tsiatis, A. S. (2002). An estimator for the proportionalhazards model with multiple longitudinal covariates measured with error. Biostatistics Tsiatis, A. A. and Davidian, M. (2004). Joint modeling of longitudinal and time-to-eventdata: An overview. Statist. Sinica Tsiatis, A. A., DeGruttola, V. and Wulfsohn, M. S. (1995). Modeling the relation-ship of survival to longitudinal data measured with error. Applications to survival andCD4 counts in patients with AIDS. J. Amer. Statist. Assoc. Venables, W. N., Smith, D. M. and the R Development Core Team (2008). AnIntroduction to R. Version 2.8.1 (2008-12-22). Verbeke, G. and Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data .Springer, New York. MR1880596 Wei, G. C. G. and Tanner, M. A. (1990). A Monte-Carlo implementation of the E–Malgorithm and the poor man’s data augmentation algorithm. J. Amer. Statist. Assoc. Wu, M. C. and Carroll, R. J. (1988). Estimation and comparison of changes in thepresence of informative right censoring by modeling the censoring process. Biometrics Wulfsohn, M. S. and Tsiatis, A. A. (1997). A joint model for survival and longitudinaldata measured with error. Biometrics Xu, J. and Zeger, S. L. (2001a). Joint analysis of longitudinal data comprising repeatedmeasures and times to events. Appl. Statist. Xu, J. and Zeger S. L. (2001b). The evaluation of multiple surrogate endpoints. Bio-metrics Eunice Kennedy ShriverNational Institute of Child Healthand Human DevelopmentNational Institutes of HealthBethesda, Maryland 20892USAE-mail: [email protected]@mail.nih.gov