Assessing the association between trends in a biomarker and risk of event with an application in pediatric HIV/AIDS
aa r X i v : . [ s t a t . A P ] O c t The Annals of Applied Statistics (cid:13)
Institute of Mathematical Statistics, 2009
ASSESSING THE ASSOCIATION BETWEEN TRENDS INA BIOMARKER AND RISK OF EVENT WITHAN APPLICATION IN PEDIATRIC HIV/AIDS
By Elizabeth R. Brown University of Washington
We present a new joint longitudinal and survival model aimedat estimating the association between the risk of an event and thechange in and history of a biomarker that is repeatedly measured overtime. We use cubic B-splines models for the longitudinal componentthat lend themselves to straight-forward formulations of the slopeand integral of the trajectory of the biomarker. The model is appliedto data collected in a long term follow-up study of HIV infected in-fants in Uganda. Estimation is carried out using MCMC methods.We also explore using the deviance information criteria, the condi-tional predictive ordinate and ROC curves for model selection andevaluation.
1. Introduction.
In longitudinal studies it is common to monitor one ormore biomarkers repeatedly over time while following participants until theoccurrence of an event. Researchers are often interested in examining boththe repeated measures and the time to event to gain an understanding ofthe underlying disease process. Additionally, the risk of an event may notdepend solely on the level of the biomarker but also on the rate at which thatbiomarker is changing or its past average level. For example, two patientsmay present with the same biomarker value, but one patient’s biomarkertrajectory may be increasing while the other’s is remaining constant. Theirprognosis may appear to be be the same if only the current value is accountedfor, but in fact the patient with the increasing biomarker may be at higherrisk. In this paper we present a model to estimate the association betweenthe risk of an event and the current value, as well as the rate of change orhistory of a longitudinally sampled biomarker. We illustrate the approachfrom a sub-study of HIVNET 012 [Guay et al. (1999); Jackson et al. (2003)]
Received June 2008; revised April 2009. Supported by a Contract from the National Institutes of Health NIAID U01 AI46702.
Key words and phrases.
HIV/AIDS, disease progression, mother-to-child transmission,joint longitudinal and survival models, biomarker change.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Applied Statistics ,2009, Vol. 3, No. 3, 1163–1182. This reprint differs from the original in paginationand typographic detail. 1
E. R. BROWN of HIV disease progression in Ugandan children who acquired HIV vertically,either in utero, during delivery or via breastmilk. In this sub-study childrenwho tested positive before 18 months of age were followed until five yearsof age with blood samples drawn every 6 months. From these samples, thelab determined viral load, total lymphocyte count (TLC), CD4 percent andother related biomarkers. One aim was to determine the predictive valueof these measures for time until death. In this manuscript we explore theassociation between change in and history of these markers and the risk ofdeath. Overall, 128 children were followed after their first positive HIV testwith lab measurements taken at regular intervals. Of these 128 infants, 70died during follow-up.In estimating the association between trends in a biomarker or set ofbiomarkers and the risk of an event, we face two important and distinctchallenges. The first is selecting the correct model when the biomarker iscollected in discrete time with error. The second is to determine how toobtain different summaries of the biomarker(s) and include them in thetime to event model.The first issue and its resolution through joint modeling of the time toevent and longitudinal marker has been studied extensively as summarizedby Tsiatis and Davidian (2004) when dealing with the current level of abiomarker and the risk of an event. In summary, survival analyses withtime-dependent covariates can be biased if we simply include the raw mea-surements in the survival analysis [Prentice (1982)]. Joint longitudinal andsurvival models resolve this issue by modeling the biomarker process overtime and including subject-specific parameters from the longitudinal modelas covariates in the survival model. These same issues exist when trying toinclude other summaries of the biomarker process in the survival model, andmay, in fact, be exacerbated.Wulfsohn and Tsiatis (1997) and Faucett and Thomas (1996) introducedlikelihood-based methods for analyzing a longitudinal marker and its asso-ciation with the time to event simultaneously. This class of models is notrestricted to linking the time-varying value of the longitudinal marker tothe time-varying hazard through a regression on the current value of longi-tudinal model. They may instead include other information from the longi-tudinal trajectories summarized by the longitudinal model. Several authorshave proposed models that group the trajectories into latent classes, thenlink the hazard and the longitudinal model by the latent class [Lin et al.(2002); Proust-Lima, Letenneur and Jacqmin-Gadda (2007); Han, Slate andPena (2007); Proust-Lima et al. (2009)]. These models can help characterizelongitudinal trajectories that may indicate a higher risk for event. In thesemodels the latent class for the trajectory of the biomarker is defined basedon the entire follow-up for the biomarker, which may have the drawback of
RENDS IN A BIOMARKER AND RISK OF EVENT using information about the biomarker from the future to estimate risk inthe present.Other descriptors from the longitudinal model have also been used tolink the trajectory of the biomarker to risk of event. Yu, Taylor and San-dler (2008) developed individual risk prediction models for prostate cancerrecurrence based on PSA trajectories. PSA was modeled using a nonlin-ear exponential decay and exponential growth model. The current value aswell as slope (first derivative) of the PSA trajectory were included as time-varying covariates in the hazard model. Ye, Lin and Taylor (2008) presenteda two-stage regression calibration approach for modeling longitudinal andtime-to-event data. Their longitudinal model included smoothing splines atthe population level with individual deviations from the smoothing splineallowed through a mean 0 integrated Wiener process and subject-specificslopes and intercepts. Both the current value and slope of the subject-specifictrajectories were included as covariates in the hazard model.In this manuscript we extend the work of Brown, Ibrahim and DeGruttola(2005) who used cubic B-splines to model the impact of multiple biomarkerson time to event to include the slope and integral of the cubic B-spline mod-els as time-varying covariates in the hazard model. They showed that cubicB-splines provided flexibility in the biomarker model that simple parametricmodels could not. Because cubic B-spline trajectory models and their slopesand integrals are linear functions of the parameters that do not increaseexponentially with time, they avoid computational instability. Additionally,because the basis functions of the cubic B-splines weight more heavily onlocal (in time) information for estimating the value of the trajectory, esti-mates of the slope early in follow-up do not rely heavily on values of thebiomarker observed late in follow-up.The paper proceeds as follows. In the next section we review cubic B-splinemodels for longitudinal data and outline the model associating change andcumulative exposure with time to event. Next, we discuss estimation andmodel selection procedures. We then show an example from HIVNET 012.We conclude with a discussion.
2. The joint longitudinal and survival model.
In this section we describea joint longitudinal and survival model to estimate the association betweenthe rate of change in a biomarker or cumulative history of a biomarker andthe risk of an event. We begin with a description of the notation and areview of the longitudinal cubic B-spline model, including expressions forthe first derivatives and integrals. We then introduce the model linking thebiomarker and its rate of change or cumulative history to the risk of event.2.1.
The longitudinal model.
Let y ijl denote the i th, i = 1 , . . . , N, sub-ject’s j th, j = 1 , . . . , m i , observation of the l th, l = 1 , . . . , L, biomarker at E. R. BROWN time t ij < T , where T denotes the end of follow-up. We define an observa-tion at time t ij to be a function of the true underlying trajectory ψ ( t ij ) pluserror, y ijl = ψ l ( t ij ) + e ijl , where the errors are independent and normally distributed such that ( e ij , . . . ,e ijL ) ′ ∼ N l (0 , Σ), where N l ( a, b ) is the l -dimensional multivariate normal dis-tribution with mean vector a and covariance matrix b . Brown, Ibrahim andDeGruttola (2005) modeled the true, but unobserved, trajectory using cubicB-splines such that ψ l ( t ) = q X k =1 β ilk B k ( t ) , (2.1)where { B k ( · ) } is a ¯ q -dimensional basis for spline functions on [0 , T ] witha fixed knot sequence, u = ( u , . . . , u q +4 ) and β il = ( β il , . . . , β ilq ) ′ is a vec-tor of subject-specific parameters of length q that determine the shape ofthe i th subject’s trajectory. We assume a hierarchical model where β il ∼ N l ( b l + X ′ i α l , V l ). In this model the effect of the covariates is modeled atthe population level where α l is a vector of parameters of length p linkingthe vector of baseline covariates X i to the longitudinal outcome and b is thevector of length q of the mean of the coefficients for the k th basis functionwhen X i , . . . , X ip = 0. Brown, Ibrahim and DeGruttola (2005) showed themerits (better mixing in the Gibbs sampler and more flexible effect of thecovariate on the longitudinal outcome) of modeling the effect of covariateson the longitudinal model in this level of the hierarchy. As a further exten-sion, we allow for a covariance structure for the spline coefficients where V l is the q × q covariance matrix of the coefficients of the basis functions.An individual’s contribution to the likelihood of the longitudinal markercan be expressed as p ( Y i | Σ , β i ) ∝ | Σ | m i exp ( − m i X j =1 ( Y ij − ψ ( t ij )) ′ Σ − ( Y ij − ψ ( t ij )) ) , (2.2)where Y ij = ( y ij , . . . , y ijL ) ′ and ψ ( t ij ) = ( ψ ( t ij ) , . . . , ψ L ( t ij )) ′ .2.2. The joint model.
We next review the form of the joint model whenwe are only interested in modeling the relationship between the level of themarker at time t and the hazard at time t . We then show how the splinemodel of the trajectory can be used to describe the rate of change in andhistory of the biomarker. Next, we incorporate these functionals into thehazard, thus extending the joint model to estimate the impact of the rateof change and history on the risk of an event. RENDS IN A BIOMARKER AND RISK OF EVENT The usual joint longitudinal and survival model assumes proportionalhazards, and the effect of the trajectories of the biomarkers on the hazardis modeled as λ ( t ) = λ ( t ) exp( γ ′ ψ ( t ) + Z ′ i ζ ) , (2.3)where λ ( t ) is the baseline hazard at time t, γ is a parameter vector of length L linking the trajectory vector at time t to the hazard at time t , and ζ is avector of parameters of length p z linking the vector of baseline covariates Z i to the hazard. Taking a likelihood approach requires some specification ofthe baseline hazard. Here, we specify a piecewise constant hazard allowingfor approximately 8–10 events in each interval, where λ ( t ) = λ j , w j ≤ t < w j +1 , j = 1 , . . . , J, where the w j ’s are the jump points with w = 0 and w J +1 = ∞ .Then the cumulative hazard, Z s i λ ( u ) e γ ′ ψ ( u )+ z ′ i ζ du, can be rewritten as e z ′ i ζ J X j =1 H ij ( β, γ, λ ) , where H ij ( β, γ, λ ) = I { s i ≥ u j − } λ j Z min( u j ,s i ) u j − e γ ′ ψ ( u )+ z ′ i ζ du (2.4)and I { s i ≥ u j − } is an indicator function which equals 1 if the event timeoccurs in or later than the j th interval and 0 otherwise. The integral in (2.4)does not have an analytical solution for the trajectory defined by cubic B-splines. Instead, for computational ease and speed, we use Gaussian quadra-ture to approximate it.Extending the model to include the relationship between another functionor functions of the time-varying covariates and the hazard requires addinganother term to the model. If we are interested in the relationship betweenrate of change of the biomarker and the hazard, we include the first deriva-tive of ψ ( t ). If we are interested in the relationship between the cumulativehistory of the biomarker and the hazard, we include the integral of ψ ( t ).For clarity, we drop the subscripts for definition of the first derivative andintegral. The first derivative of the trajectory function as shown by de Boor[(2001), page 116] can be expressed as ψ ′ ( t ) = B (2) ∗ ( t ) ′ β, (2.5) E. R. BROWN where B (2) ∗ ( t ) is a vector of length q with the j th element equal to B (2) j ( t ) u j +3 − u j − B (2) j +1 ( t ) u j +4 − u j +1 , and B (2) ( t ) is the quadratic B-spline based on the same sequenceof knots as the original B-spline in (2.1) with the first and last knot removed.Equation (2.5) is written as a linear function of the elements of β i ; however,it is a quadratic B-spline and still retains many of the desirable propertiesof B-splines mentioned earlier.The general idea for calculating the integral of the cubic B-spline was laidout by de Boor [(2001), page 128]. We derived the integral of the trajectoryup to time t to be Z t ψ ( v ) dv = ψ ( − ( t ) = B ∗ (4) ( t ) ′ β, (2.6)where B ∗ (4) ( t ) is a vector of length q with j th element equal to P q +1 k = j +1 B (4) k ( t ) × ( u j +4 − u j ), and { B (4) k ( · ) } is the vector of q + 1-dimensional basis of a quar-tic B-spline based on the knot vector u augmented by two arbitrary knots, u < u and u q +5 > u q +4 .To link the rate of change and history of the trajectories to the hazard,we use the following regression model: λ ( t ) = λ ( t ) exp { γ ′ ψ ( t ) + γ ′ s ψ ′ ( t ) + γ ′ h ψ ( − ( t ) + Z ′ i ζ } , (2.7)where γ s is the L -length vector of parameters linking the L -length vectorof the slopes of the trajectories at time t , ψ ′ ( t ) to the hazard at time t and γ h is the L -length vector of parameters linking the L -length vector of thecumulative histories of the trajectories up to time t , ψ − ( t ) to the hazardat time t .As in the case where only the trajectory value at time t is linked to thehazard at time t , the cumulative hazard can be written as e z ′ i ζ J X j =1 H ij ( β, γ, γ s , γ h , λ ) , where H ij ( β, γ, γ s , γ h , λ )= I { s i ≥ u j − } λ j (2.8) × Z min( u j ,s i ) u j − exp { γ ′ ψ ( u ) + γ ′ s ψ ′ ( u ) + γ ′ h ψ ( − ( u ) } du and γ s = ( γ s , . . . , γ sL ) ′ and γ h = ( γ h , . . . , γ hL ) ′ . Because equation (2.8) doesnot have a closed form solution, we evaluate it numerically using Gaussianquadrature. RENDS IN A BIOMARKER AND RISK OF EVENT The distribution for an individual’s time to event, s i , given the trajectoryfunction and choice of hazard regression, is given by f ( s i , ν i | Y i ) = λ ( s i ) ν i exp { ν i ( γ ′ ψ ( s i ) + γ ′ s ψ ′ ( s i ) + γ ′ h ψ ( − ( s i ) + z ′ i ζ ) } (2.9) × exp ( − e z ′ i ζ J X j =1 H ij ( β, γ, γ s , γ h , λ ) ) , where ν i is the censoring indicator for subject i .
3. Estimation.
We can now express the i th subject’s contribution to thejoint likelihood function as f ( Y i , s i , ν i ) = f ( s i , ν i | Y i ) × f ( Y i ) ,f ( Y i , s i , ν i ) ∝ λ ( s i ) ν i exp { ν i ( γ ′ ψ ( s i ) + γ ′ s ψ ′ ( s i ) + γ ′ h ψ ( − ( s i ) + z ′ i ζ ) }× exp ( − e z ′ i ζ J X j =1 H ij ( β, γ, γ s , γ h , λ ) ) × | Σ | m i exp ( − m i X j =1 ( Y ij − ψ ( t ij )) ′ Σ − ( Y ij − ψ ( t ij )) ) . We then specify the prior distributions for the parameters in the likelihoodas follows. We assume ( γ, γ s , γ h ) ′ ∼ N L ( G , G ), Σ ∼ Wishart ν Σ ( S − ) and λ j ∼ gamma( d j , d j ). We may also specify priors on the hyperparameters of β , α l ∼ N p ( C l , C l ), b l ∼ N p ( A l , A l ) and V − l ∼ Wishart ν v l ( S − v l ). Here,Wishart ν ( S − ) denotes the Wishart distribution with ν degrees of freedomand scale matrix S − and gamma( a, b ) denotes the gamma distribution withshape parameter a and scale parameter b . The prior distributions were cho-sen to be as general as possible while still being proper and conjugate to thelikelihood when possible.We use the Gibbs sampler [Gelfand and Smith (1990)] to sample fromthe posterior distribution of the parameters. Because the full condition-als of γ , γ s , γ h and ζ are log-concave, we use adaptive rejection sampling(ARS) [Gilks and Wild (1992)] to sample from them. We use the slice sam-pler [Neal (2003)] to sample the random effects, β ij , j = 1 , . . . , q, i = 1 , . . . , N .The full conditionals of λ , Σ, V and b are sampled from directly. The esti-mation procedure is implemented in R [R Development Core Team (2006)]and C with code available from the author.
4. Model comparison.
We examine two statistics for model comparison,the deviance information criterion (DIC) [Spiegelhalter et al. (2002)] and theconditional predictive ordinate (CPO) [Gelfand, Dey and Chang (1992)].
E. R. BROWN
The DIC is a measure of the deviance penalized for the number of pa-rameters in the model which may be difficult to ascertain in hierarchicalmodels and is therefore estimated. The DIC is the sum of the deviance es-timated using the posterior estimates of the parameters, D ( ¯Θ), and twicethe effective number of parameters, p D . The effective number of parametersis estimated by p D = D (Θ) − D ( ¯Θ), where D (Θ) is the posterior mean ofthe deviance (the average of the deviances calculated using the estimatedparameters at each step of the MCMC sampler). For the model presentedin this paper, the DIC can be expressed as DIC = 2 1 G G X g =1 N X i =1 log { f ( s i , ν i , Y i | Θ ( g ) ) } − N X i =1 log { f ( s i , ν i , Y i | ¯Θ) } , where Θ ( g ) denotes the parameter samples at the g th, g = 1 , . . . , G , iterationof the Gibbs sampler and ¯Θ represents the means of the posterior samples.A smaller DIC indicates a better fit when comparing models.For the i th observation, the CPO statistic is defined as CPO i = f ( s i , ν i , Y i | D ( − i ) ) = Z f ( s i , ν i , Y i | Θ , D i ) π (Θ | D ( − i ) ) d Θ , (4.1)where Θ denotes the model parameters, D i denotes the i th patient’s covari-ate data and D ( − i ) denotes the covariate data for all the patients exceptthe i th patient. We cannot obtain a closed form solution for (4.1); however,Chen, Shao and Ibrahim [(2000), Chapter 10] show that a Monte Carloapproximation of (4.1) is \ CPO i = G G X g =1 f ( s i , ν i , Y i | Θ ( g ) ) ! − , where Θ ( g ) denotes the parameter samples at the g th, g = 1 , . . . , G , itera-tion of the Gibbs sampler. A large CPO value indicates a better fit. Wecan compare different models using the sums of the logs of the CPOs of theindividual observations, also known as the log pseudo-marginal likelihood(LPML). Models with greater LPML = P log( \ CPO i ) values will indicate abetter fit. Both the DIC and LPML are designed to measure a model’s pre-dictive ability, although the DIC is based on a penalized deviance approachand the LPML is based on a cross-validated approach.
5. Application.
HIVNET 012, conducted in Uganda, was a double-blindedcontrolled randomized trial of single dose nevirapine for the mother and new-born infant versus AZT administered only to the mother to prevent motherto child transmission (MTCT) of HIV. In spite of the success of nevaripinein reducing the risk of transmission, many infants still experienced MTCT
RENDS IN A BIOMARKER AND RISK OF EVENT of HIV. To better understand HIV infection in young children, 128 of thoseinfants were enrolled in a long term follow-up study and followed until age 5or death. CD4 cell percent and HIV viral load are known indicators of dis-ease progression. TLC is also being studied for its potential use in resource-poor settings where routine CD4 and viral load monitoring may be costprohibitive. In this section we examine the association between longitudinalmeasures of CD4 cell percent, viral load and TLC and time until death usingseparate models (one biomarker per model) in the HIVNET 012 long termfollow-up study using the methods proposed in the previous sections.Seventy infants died during follow-up. Jump points for the baseline hazardfunction were selected based on quantiles of the event times so that approxi-mately 8 events occurred between the jump points. Infants had between oneand thirteen longitudinal measures with a median of four. Overall, therewere a total of 594 measurements of CD4 percent, 603 measurements of vi-ral load and 763 measurements of TLC. 13% of the CD4 percent measures,7% of the viral load measures and 16% of the TLC measures are taken attime 0. In this analysis we placed the knots for the splines based on thequantiles of the data; therefore, this clumping of measurement times limitsthe number of knots we could select before we start placing multiple knotsat 0, making the slope undefined. Therefore, for CD4 percent, we fit modelswith q = 5 , . . . ,
10. For TLC, we fit models with q = 5 , . . . ,
9. For compari-son to potential models with more basis functions, we also fit models withequally spaced knots with q = 11 for CD4 percent and q = 10 for TLC. Forviral load, we can fit models with q = 5 , . . . ,
19. Additionally, we includedage at first positive HIV test as a covariate in the hazard model; therefore,the interpretation of ζ is the change in the log hazard associated with a 1month increase in the age of the infant at the time of the first positive HIVtest. We implemented the Gaussian quadrature procedure using 10 points.We also ran models with higher q using 50 points and found no differencein the estimates.Table 1 shows the estimated values of the parameters of interest along withtheir 95% credible intervals obtained from fitting the data to three versionsof the hazard model shown in equation (2.7), one with γ s = 0 and γ h = 0,one with γ h = 0 and one with γ s = 0, for a range of q . The Gibbs samplerwas run twice with different starting values and seeds for the random num-ber generator for 100,000 iterations for each model with a burn-in of 10,000.Samples from every 10th iteration were saved to reduce possible autocorre-lation. The resulting sample was used to compute parameter estimates andcredible intervals as well as DIC and LPML. Based on DIC, the minimumnumber of basis functions for a cubic B-spline, q = 5, is not adequate forany of the three longitudinal outcomes. Although not shown, models withequally spaced knots never outperformed the models with knots based onquantiles according to DIC or LPML. There were no meaningful differences E . R . B R O W N Table 1
Results from the models for the three markers of HIV disease progression γ γ s γ h ( × ζ DIC LPML VL q = 5 current 1 .
34 (0 . , . − .
08 ( − . , − .
01) 1644 − .
57 (1 . , .
18) 5 .
23 (1 . , . − .
07 ( − . , .
01) 1607 − .
72 ( − . , .
43) 5 .
56 (1 . , . − .
09 ( − . , − .
02) 1632 − q = 6 current 1 .
45 (0 . , . − .
08 ( − . , − .
01) 1635 − .
88 (1 . , .
53) 2 .
79 (1 . , . − .
07 ( − . , .
01) 1590 − .
95 (0 . , .
73) 4 .
26 ( − . , . − .
09 ( − . , − .
01) 1633 − q = 7 current 1 .
58 (1 . , . − .
08 ( − . , − .
01) 1649 − .
03 (1 . , .
70) 1 .
49 (0 . , . − .
07 ( − . , .
02) 1647 − .
20 (0 . , .
99) 3 .
28 ( − . , . − .
09 ( − . , − .
01) 1653 − q = 8 current 1 .
64 (1 . , . − .
09 ( − . , − .
01) 1646 − .
93 (1 . , .
62) 0 .
88 (0 . , . − .
09 ( − . ,
0) 1616 − .
22 (0 . , .
04) 3 .
84 ( − . , . − .
09 ( − . , − .
01) 1655 − q = 9 current 1 .
62 (1 . , . − .
08 ( − . , − .
01) 1662 − .
92 (1 . , .
66) 0 .
98 ( − . , . − .
08 ( − . , .
00) 1633 − .
26 (0 . , .
07) 3 .
27 ( − . , . − .
10 ( − . , − .
03) 1668 − q = 5 current − .
83 ( − . , − . − .
12 ( − . , − .
05) 1767 − − .
87 ( − . , − . − .
81 ( − . , . − .
13 ( − . , − .
05) 1764 − − . − . , − . − .
05 ( − . , . − .
13 ( − . , − .
05) 1747 − q = 6 current − .
81 ( − . , − . − .
12 ( − . , − .
04) 1738 − − .
91 ( − . , − . − .
47 ( − . , . − .
12 ( − . , − .
05) 1728 − − .
57 ( − . , − . − .
76 ( − . , . − .
13 ( − . , − .
05) 1722 − q = 7 current − .
96 ( − . , − . − .
12 ( − . , − .
05) 1674 − − .
98 ( − . , − . − .
38 ( − . , . − .
12 ( − . , − .
05) 1672 − − .
77 ( − . , − . − . − . , . − .
13 ( − . , − .
05) 1666 − R E N D S I NA B I O M A R K E R AN D R I S KO F E V E N T Table 1 (Continued) γ γ s γ h ( × ζ DIC LPML q = 8 current − .
01 ( − . , − . − .
13 ( − . , − .
05) 1690 − − .
03 ( − . , − . − .
26 ( − . , . − .
13 ( − . , − .
05) 1686 − − .
83 ( − . , − . − .
14 ( − . , . − .
13 ( − . , − .
05) 1684 − q = 9 current − .
94 ( − . , − . − .
12 ( − . , − .
05) 1678 − − .
94 ( − . , − .
59) 0 .
07 ( − . , . − .
12 ( − . , − .
05) 1678 − − .
74 ( − . , − . − .
95 ( − . , . − .
13 ( − . , − .
05) 1666 − q = 5 current − .
20 ( − . , . − .
11 ( − . , − .
04) 3653 − − .
31 ( − . , − . − .
77 ( − . , − . − .
11 ( − . , − .
04) 3617 − − .
22 ( − . , .
06) 0 .
12 ( − . , . − .
10 ( − . , − .
04) 3649 − q = 6 current − .
45 ( − . , − . − .
11 ( − . , − .
05) 3485 − − .
47 ( − . , − . − .
34 ( − . , . − .
11 ( − . , − .
04) 3474 − − .
70 ( − . , − .
38) 1 .
46 (0 . , . − .
11 ( − . , − .
04) 3458 − q = 7 current − .
39 ( − . , − . − .
11 ( − . , − .
04) 3521 − − .
40 ( − . , − . − .
11 ( − . , . − .
11 ( − . , − .
04) 3516 − − .
54 ( − . , − .
24) 0 .
96 ( − . , . − .
11 ( − . , − .
04) 3504 − q = 8 current − .
51 ( − . , − . − .
12 ( − . , − .
04) 3465 − − .
35 ( − . , − .
10) 0 .
02 ( − . , . − .
10 ( − . , − .
04) 3451 − − .
66 ( − . , − .
40) 1 .
24 (0 . , . − .
11 ( − . , − .
04) 3451 − q = 9 current − .
45 ( − . , − . − .
10 ( − . , − .
04) 3367 − − .
44 ( − . , − .
21) 0 .
01 ( − . , . − .
10 ( − . , − .
04) 3368 − − .
54 ( − . , − .
27) 0 .
77 ( − . , . − .
11 ( − . , − .
04) 3367 − E. R. BROWN in the estimates between the results from the two samplers. Convergence wasassessed using diagnostic tools provided in the CODA package [Plummer etal. (2006)].For viral load, DIC selected q = 6 for all three models. The LPML in-creased as q increased except for the model with γ and γ s where it selected q = 7; however, the value was not very different than that for q = 6. Here, wefocus on the results for q = 6. The point estimate of γ indicates an increasedrisk of death with increasing viral load such that a 10-fold increase in viralload is associated with a 16.3-fold increase in the hazard [95% credible inter-val (CI): (4.4, 74.5)]. The rate of change in viral load is also associated withrisk of death. At a known level of viral load, a unit increase per month in therate of change of log viral load is associated with a 2.8-fold increase in thehazard [95% CI: (1.15, 7.4)]. The estimate of the strength of this associationdecreases as q increases. However, increasing q also introduces more fluc-tuations in the estimate of the slope over time that may not be supportedby the data. Additionally, the history of viral load is also associated withrisk of death, with a one unit increase in the mean value of the trajectoryof log viral load times a unit change in the length (in months) of follow-upbeing associated with 1.04-fold increase in the hazard. In plainer terms, iftwo infants have been followed for 7 months, with a difference in mean val-ues equal to 1, the hazard ratio would be exp(7 ∗ . / .
35. If after 12months the difference in mean levels was still equal to 1, the hazard ratiowould be 1.68. Figure 1 shows the trajectory fits for viral load with q = 6.The model fits a variety of shapes suggested by the data. For example, theinitial steep rise in infants 2, 10 and 12 is fit well, as is the initial decreasein infants 6 and 9.For CD4 percent, the DIC selected q = 7 for all three models. The LPMLselected q = 7 for the current value and current value plus slope models and q = 8 for the current value plus history model (although the value was closeto that for q = 7). Here, we focus on the results for q = 7. A 10 unit decreasein CD4 percent is associated with a 2.6-fold increase in the hazard [95%CI: (1.86, 3.82)]. There is no suggestion of a further association betweenrisk of death and the change in CD4 percent or its history. Figure 2 showsthe fit of the spline model to the observed CD4 percent data for 12 infants.The majority of the infants in Figure 2 have an initial sharp decline in CD4percent which is fit well by the model without forcing decreases where thedata do not suggest it (infant 6).For TLC, the DIC and LPML selected q = 9. The density of measurementsnear time zero forced quantile-based knots very close together if not equalshortly after infection leading to undefined slopes for models with q > q = 9. A possible association between TLC andrisk of death was suggested in this model with a 1000 unit increase in TLC RENDS IN A BIOMARKER AND RISK OF EVENT Fig. 1.
Longitudinal profiles of viral load for 12 infants. Circles mark the observed values.The solid black line indicates the fit from the model. The red line represents the fittedcumulative value. The green line represents the fitted slope. A vertical grey line indicatestime of death if it was observed. E. R. BROWN
Fig. 2.
Longitudinal profiles of CD4 percent for 12 infants. Circles mark the observedvalues. The solid black line represents the fit from the model. The red line representsthe fitted cumulative value. The green line represents the fitted slope. A vertical grey lineindicates time of death if it was observed.
RENDS IN A BIOMARKER AND RISK OF EVENT being associated with a 34% [95% CI: (20, 49)] decrease in the risk of death.Figure 3 shows the fit from the longitudinal models for TLC.The models all estimate a similar association between age at first positivetest and risk of death, with a one month increase in age at first positive testassociated with an approximate 10% reduction in risk given the trajectoryof the longitudinal marker.Figure 4 shows the posterior estimate of the cumulative hazards for thethree markers when only the current value of the marker is included in themodel. TLC provide better fits to the Kaplan–Meier estimate in the firstyear, while viral load provides a better fit between 2 and 3 years.We propose using ROC curves for censored data [Heagerty, Lumley andPepe (2000)] as an additional comparison of joint longitudinal and survivalmodels. We plotted ROC curves for predicting death within 1 year basedon linear predictor, γψ ( t ) + γ s ψ ′ ( t ) + γ h ψ ( − ( t )) + Z ′ i ζ , at t = 6, 12 and18 months after infection (Figure 5). Taking t = 6 months as an example,we treat the baseline time for survival as 6 months, and calculate the ROCcurve for death within 12 months (18 months after infection) with the linearpredictor calculated at 6 months as the biomarker. Infants who died or werecensored before 6 months are not included. The same procedure was usedfor 12 and 18 months. These results do not suggest any large improvementin prediction when either the rate of change in or history of the biomarkerare included in the models, except possibly for including slope in the viralload model. This agrees with the results from the model and the DIC andLPML. Overall, TLC does not compare favorably to either viral load or CD4percent in predicting death within one year, and viral load appears to bethe best predictor.
6. Discussion.
Understanding the relationship between trends in abiomarker and risk of an event can yield important insight into the mecha-nisms of disease progression. Clearly, this is recognized scientifically, as the
Journal of the American Medical Association recently made an exception toits own policies on publishing data over five years old to report an analysisof the association between CD4 slopes estimated from linear mixed effectsmodels and time to development of AIDS or death in antiretroviral naiveHIV-infected adults [Mellors et al. (2007)]. Here, we present a model thatgoes several steps further to addressing that question by jointly modelingthe time to event and the slope and looking at local changes as opposedto long term trends. Additionally, we do not restrict ourselves to a linearmodel that assumes constant rate of change. This model is motivated andillustrated by a long-term follow-up study of HIV-infected children in Africa.However, it can be used in many settings where understanding the relation-ship between trends in a biomarker and time to an event is of interest. Other E. R. BROWN
Fig. 3.
Longitudinal profiles of TLC for 16 infants with q = 6 . Circles mark the observedvalues. The solid black line indicates the fitted trajectory from the model. The red linerepresents the fitted cumulative value. The green line represents the fitted slope. A verticalgrey line indicates time of death if it was observed. RENDS IN A BIOMARKER AND RISK OF EVENT . . . . . . Months since first positive HIV test C u m u l a t i v e ha z a r d VLCD4TLC
Fig. 4.
Fitted cumulative hazard curves for the selected models for TLC, viral load andCD4 percent. The solid and dashed stepped lines represent the Kaplan–Meier fit with 95%confidence intervals, with the small vertical lines representing censored observation times. examples may include cognitive function and death or subclinical coronarydisease and clinical coronary events.This model gives interpretable parameters, although these values are dif-ficult to translate into practice. For example, although we estimate that a10-fold (1 log unit) increase per month in viral load is associated with a2.8-fold increase in the hazard, we cannot easily measure the instantaneouschange in viral load in practice. To do so, we would have to measure it fre-quently, which is not feasible in practice, especially in resource-poor settings.However, this approach may indicate if there is more information available incollected information than just the current value and what that informationmight be (change versus average history, for example).We have proposed a model that provides additional insight into the bi-ological processes of disease progression by linking the hazard of the eventto the rate of change or cumulative history of a biomarker. This model ex-pands the possibilities of examining disease progression within the class ofjoint longitudinal and survival models.APPENDIX: SAMPLING FROM THE POSTERIORWe use Gibbs sampling to sample from the joint posterior distributionof the parameters: β , α , γ , γ s , γ h , λ , b o and V . The joint posterior doesnot have a closed form; however, given that the conditional posteriors eitherhave a closed form or are log-concave, implementation of the Gibbs sampleris straight-forward. Let D denote the data and rest denote the remaining E. R. BROWN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S en s i t i v i t y T L C V i r a l l oad CD pe r c en t Fig. 5.
ROC curves for predicting death within one year based on the hazard functioncalculated at 6, 12 and 18 months after initial positive HIV test for current value (solidline), current value and slope (dashed line) and current value and history (dotted line)models. The areas under the ROC curves are shown with each plot. parameters. Then at each iteration of the Gibbs sampler, we proceed asfollows:1. Let β il = ( β il , . . . , β ilq ) ′ , i = 1 , . . . , N, l = 1 , . . . , L . Then use ARS to sam-ple [ β il | rest , D ] from p ( β i | rest , D ) ∝ exp ( − " m i X j =1 ( Y ij − ψ β ( t ij )) ′ Σ − ( Y ij − ψ β ( t ij ))+ L X l =1 ( β il − b l − X ′ i α l ) ′ V − l ( β il − b l − X ′ i α l ) RENDS IN A BIOMARKER AND RISK OF EVENT × exp ( ν i ( γ ′ ψ ( s i ) + γ ′ s ψ ′ ( s i ) + γ ′ h ψ ( − ( s i ) + z ′ i ζ ) − e z ′ i ζ J X j =1 H ij ( β, γ, γ s , γ h , λ ) ) .
2. Sample[ V − l | rest , D ] ∼ Wishart S − v l + N X i =1 ( β il − b l − x ′ i α )( β il − b l − x ′ i α ) ′ ! − ,N + ν v l ! .
3. Sample[Σ − | rest , D ] ∼ Wishart S − + N X i =1 m i X j =1 ( Y ij − ψ ( t ij ))( Y ij − ψ ( t ij )) ′ ! − , N X i =1 m i + ν Σ ! .
4. Let b l = ( b l , . . . , b lq ) ′ , i = 1 , . . . , N, l = 1 , . . . , L . Then sample[ b l | rest , D ] ∼ N ( µ b l , Σ b l ) , where µ b l = Σ b l V − l N X i =1 ( β il − x ′ i α ) + A − l A l ! and Σ b l = ( N V − l + A − l ) − .
5. Use ARS to sample [ γ, γ s , γ h | rest , D ] from p ( γ | rest , D ) ∝ exp ( N X i =1 ν i ( γ ′ ψ ( s i ) + γ ′ s ψ ′ ( s i ) + γ ′ h ψ ( − ( s i ) + z ′ i ζ ) − e z ′ i ζ J X j =1 H ij ( β, γ, γ s , γ h , λ ) !) × exp (cid:26) −
12 (( γ, γ s , γ h ) ′ − g ) ′ g − (( γ, γ s , γ h ) ′ − g ) (cid:27) . E. R. BROWN
6. Sample [ λ j | rest , D ] = gamma( d j + n j , P Ni =1 e z ′ i ζ H ij ( β, γ, γ s , γ h ,
1) + d j ) , where n j is the number of events in the j th interval and H ij ( β, γ,
1) is H ij ( β, γ, γ s , γ h , λ ) evaluated with λ j = 1.7. Sample [ α | rest , D ] ∼ N p ( µ α , Σ α ) , where µ α = Σ α N X i =1 X i q V − l ( β il − b l ) + C − C ! and Σ α = N X i =1 X i q V − l q X ′ i + C − ! − , where q is a vector of ones of length q . Acknowledgments.
The author thanks an Associate Editor and the ref-erees for their comments that have greatly improved the presentation of thearticle. REFERENCES
Brown, E. R., Ibrahim, J. G. and
DeGruttola, V. (2005). A flexible B-spline modelfor multiple longitudinal biomarkers and survival.
Biometrics Chen, M.-H., Shao, Q.-M. and
Ibrahim, J. G. (2000).
Monte Carlo Methods in BayesianComputation . Springer-Verlag, New York. MR1742311 de Boor, C. (2001).
A Practical Guide to Splines (revised ed.). Springer-Verlag, NewYork. MR1900298
Faucett, C. L. and
Thomas, D. C. (1996). Simultaneously modelling censored survivaldata and repeatedly measured covariates: A Gibbs sampling approach.
Statist. Med. Gelfand, A. E., Dey, D. K. and
Chang, H. (1992). Model determination using predic-tive distributions with implementation via sampling-based methods (with discussion).In
Bayesian Statistics 4 (J. M. Bernardo, J. O. Berger, A. P. Dawid and A. F. M.Smith, eds.). Oxford University Press, Oxford. MR1380275
Gelfand, A. E. and
Smith, A. F. M. (1990). Sampling-based approaches to calculatingmarginal densities.
J. Amer. Statist. Assoc. Gilks, W. R. and
Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling.
Appl. Statist. Guay, L., Musoke, P., Fleming, T., Bagenda, D., Allen, M., Nakabiito, C.,Sherman, J., Bakaki, P., Ducar, C., Deseyve, M., Emel, L., Mirochnick, M.,Fowler, M., Mofenson, L., Miotti, P., Dransfield, K., Bray, D., Mmiro, F. and
Jackson, J. (1999). Intrapartum and neonatal single-dose nevirapine comparedwith zidovudine for prevention of mother-to-child transmission of HIV-1 in Kampala,Uganda: HIVNET 012 randomised trial.
Lancet
Han, J., Slate, E. H. and
Pena, E. A. (2007). Parametric latent class joint model for alongitudinal biomarker and recurrent events.
Statist. Med. Heagerty, P. J., Lumley, T. and
Pepe, M. S. (2000). Time-dependent ROC curvesfor censored survival data and a diagnostic marker.
Biometrics Jackson, J., Musoke, P., Fleming, T., Guay, L., Bagenda, D., Allen, M., Nakabi-ito, C., Sherman, J., Bakaki, P., Owor, M., Ducar, C., Deseyve, M., Mwatha,A., Emel, L., Duefield, C., Mirochnick, M., Fowler, M., Mofenson, L., Miotti,P., Gigliotti, M., Bray, D. and
Mmiro, F. (2003). Intrapartum and neonatal single-dose nevirapine compared with zidovudine for prevention of mother-to-child transmis-sion of HIV-1 in Kampala, Uganda: 18-month follow-up of the HIVNET 012 randomisedtrial.
Lancet
Lin, H., Turnbull, B. W., McCulloch, C. E. and
Slate, E. H. (2002). Latent classmodels for joint analysis of longitudinal biomarker and event process data: Applicationto longitudinal prostate-specific antigen readings and prostate cancer.
J. Amer. Statist.Assoc. . MR1947272 Mellors, J. W., Margolick, J. B., Phair, J. P., Rinaldo, C. R., Detels, R.,Jacobson, L. P. and
Munoz, A. (2007). Prognostic value of HIV-1 RNA, CD4 cellcount, and CD4 cell count slope for progression to AIDS and death in untreated HIV-1infection.
J. Amer. Med. Assoc.
Neal, R. M. (2003). Slice sampling.
Ann. Statist. Plummer, M., Best, N., Cowles, K. and
Vines, K. (2006). CODA: Conver-gence diagnosis and output analysis for MCMC.
R News http://CRAN.R-project.org/doc/Rnews/ . Proust-Lima, C., Joly, P., Dartigues, J. and
Jacqmin-Gadda, H. (2009). Joint mod-elling of multivariate longitudinal outcomes and a time-to-event: A nonlinear latent classapproach.
Comput. Statist. Data Anal. Prentice , R. L. (1982). Covariate measurement errors and parameter estimation in afailure time regression model.
Biometrika Proust-Lima, C., Letenneur, L. and
Jacqmin-Gadda, H. (2007). A nonlinear latentclass model for joint analysis of multivariate longitudinal data and a binary outcome.
Statist. Med. R Development Core Team (2006).
R: A Language and Environment for StatisticalComputing , R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. Available at . Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and van der Linde, A. (2002).Bayesian measures of complexity and fit.
J. Roy. Statist. Soc. Ser. B Tsiatis, A. and
Davidian, M. (2004). An overview of joint modeling of longitudinal andtime-to-event data.
Statist. Sinica Wulfsohn, M. S. and
Tsiatis, A. A. (1997). A joint model for survival and longitudinaldata measured with error.
Biometrics Ye, W., Lin, X. and
Taylor, J. M. G. (2008). Semiparametric modeling of longitudinalmeasurements and time-to-event data—A two-stage regression calibration approach.
Biometrics Yu, M., Taylor, J. M. G. and
Sandler, H. M. (2008). Individual prediction in prostatecancer studies using a joint longitudinal survival cure model.
J. Amer. Statist. Assoc.