A dynamic Bayesian nonlinear mixed-effects model of HIV response incorporating medication adherence, drug resistance and covariates
Yangxin Huang, Hulin Wu, Jeanne Holden-Wiltse, Edward P. Acosta
aa r X i v : . [ s t a t . A P ] A p r The Annals of Applied Statistics (cid:13)
Institute of Mathematical Statistics, 2011
A DYNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODELOF HIV RESPONSE INCORPORATING MEDICATIONADHERENCE, DRUG RESISTANCE AND COVARIATES By Yangxin Huang, Hulin Wu, Jeanne Holden-Wiltse andEdward P. Acosta
University of South Florida, University of Rochester, University ofRochester and University of Alabama
HIV dynamic studies have contributed significantly to the under-standing of HIV pathogenesis and antiviral treatment strategies forAIDS patients. Establishing the relationship of virologic responseswith clinical factors and covariates during long-term antiretroviral(ARV) therapy is important to the development of effective treat-ments. Medication adherence is an important predictor of the effec-tiveness of ARV treatment, but an appropriate determinant of ad-herence rate based on medication event monitoring system (MEMS)data is critical to predict virologic outcomes. The primary objectiveof this paper is to investigate the effects of a number of summarydeterminants of MEMS adherence rates on virologic response mea-sured repeatedly over time in HIV-infected patients. We developeda mechanism-based differential equation model with considerationof drug adherence, interacted by virus susceptibility to drug andbaseline characteristics, to characterize the long-term virologic re-sponses after initiation of therapy. This model fully integrates viralload, MEMS adherence, drug resistance and baseline covariates intothe data analysis. In this study we employed the proposed model andassociated Bayesian nonlinear mixed-effects modeling approach to as-sess how to efficiently use the MEMS adherence data for prediction ofvirologic response, and to evaluate the predicting power of each sum-mary metric of the MEMS adherence rates. In particular, we intendto address the questions: (i) how to summarize the MEMS adherencedata for efficient prediction of virologic response after accounting forpotential confounding factors such as drug resistance and covariates,Received May 2009; revised September 2009. Supported in part by NIAID/NIH Grant AI080338 and MSP/NSA Grant H98230-09-1-0053 to Y. Huang, and NIH Grants AI50020, AI078498, AI078842 and AI087135 to H.Wu.
Key words and phrases.
Bayesian mixed-effects models, confounding factors, HIV dy-namics, longitudinal data, MEMS adherence assessment, time-varying drug efficacy, virusresistance.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Applied Statistics ,2011, Vol. 5, No. 1, 551–577. This reprint differs from the original in paginationand typographic detail. 1
HUANG ET AL.and (ii) how to evaluate treatment effect of baseline characteristicsinteracted with adherence and other clinical factors. The approachis applied to an AIDS clinical trial involving 31 patients who hadavailable data as required for the proposed model. Results demon-strate that the appropriate determinants of MEMS adherence ratesare important in order to more efficiently predict virologic response,and investigations of adherence to ARV treatment would benefit frommeasuring not only adherence rate but also its summary metric as-sessment. Our study also shows that the mechanism-based dynamicmodel is powerful and effective to establish a relationship of virologicresponses with medication adherence, virus resistance to drug andbaseline covariates.
1. Introduction.
The revolution in human immunodeficiency virus (HIV)treatment has brought diagnostic tests that can accurately measure levels ofHIV in blood. Resulting data show (plasma) viral load (HIV-1 RNA copies orRNA copies) to be an important predictor of the risk of progression to AIDS.The antiretroviral (ARV) agents, which include potent protease inhibitors(PIs) are, however, not a cure for HIV infection. While many patients benefitfrom ARV treatment, others do not benefit or only experience a temporarybenefit. There are several reasons why treatment fails, of which poor patientadherence to ARV therapy is a leading factor [Ickovics and Meisler (1997);Paterson et al. (2000)]. Thus, assessment of adherence within AIDS clinicaltrials is a critical component of the successful evaluation of therapy out-comes. Maintaining adherence may be particularly difficult when the drugregimen is complex or side-effects are severe, as is often the case for currentHIV therapy [Ickovics and Meisler (1997)].The measurement of adherence remains problematic; a standard definitionof adherence and completely reliable measures of adherence are lacking.Nevertheless, there has been substantial progress in both of these areas inthe past few years. First, it appears that higher levels of adherence areneeded for HIV disease than other diseases to achieve the desired therapeu-tic benefit [Paterson et al. (2000)]. Second, better appreciation of the valueand limitations of different adherence measurements has been addressed[Bova et al. (2005)]. In AIDS clinical trials adherence to medication regimenis currently measured by two methods: by use of questionnaires (patientself-reporting and/or face-to-face interview) and by use of electronic com-pliance monitoring (Medication Event Monitoring System [MEMS]) caps.The MEMS is often used as an objective adherence measure. It consists ofa computer chip in the cap of a medication bottle that records each timethe bottle is opened. The results can be downloaded, printed out and ana-lyzed. It demonstrates that medication-taking patterns are highly variableamong patients [Kastrissios et al. (1998)] and that they often give a moreprecise measure of adherence than self-report [Arnsten et al. (2001)]. How-ever, MEMS data are also subject to error and are not widely available in
YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL the clinical setting. Adherence assessment by self-report is usually evalu-ated by a patient’s ability to recall their medication dosing during a specifictime interval. Finally, it is important to note that the measurement of viralload levels is of special utility as an indirect measure of adherence in HIVtherapeutics. It has been argued that this is not an adherence measure be-cause other factors may influence viral load (drug resistance, etc.). However,there is a tight correlation between viral load and adherence [Paterson et al.(2000); Haubrich et al. (1999)].Viral dynamic models can be formulated through ordinary differentialequations (ODE), but there has been only limited development of statisticalmethodologies for assessing their agreement with observed data. Currentlythere also are substantial knowledge gaps between theoretical HIV dynam-ics and the role of many clinical factors. In developing long-term dynamicmodeling, this paper will address these problems by utilizing time-specificinformation, such as drug adherence and susceptibility factors, on the bio-logical mechanism of HIV dynamics to achieve more realistic and accuratecharacterization of the relationship between clinical/drug factors and viro-logic response. Several studies [Arnsten et al. (2001); Levine et al. (2006)]investigated the association between virologic responses and adherence as-sessed by both questionnaire and MEMS data. The results indicated thatthe MEMS cap adherence data may not be correlated better to virologic re-sponse compared to the questionnaire adherence data unless the MEMS capdata are summarized in an appropriate way. Further, Huang et al. (2008),Labb´e and Verotta (2006), Liu et al. (2007) and Vrijens et al. (2005) modeledthe relationship between virologic response and adherence rate using ques-tionnaire data and MEMS data averaged by each interval between studyvisits or weekly basis, but no significant differences were found in predictingvirologic response. Along with this line, this paper will investigate differentdeterminants of the adherence rate based on MEMS data from an AIDS clin-ical trial study [Hammer et al. (2002)] and compare their performance forpredicting a virological response. We employed the proposed mechanism-based dynamic model to assess how to efficiently use the adherence databased on MEMS to predict virological response. In particular, we intendto address the questions (i) how to summarize the MEMS adherence datafor efficient prediction of virological response after accounting for potentialconfounding factors such as drug resistance and baseline covariates, and (ii)how to evaluate treatment effect of baseline characteristics interacted withMEMS adherence and other clinical factors.The purpose of this paper is to describe a reparameterized ODE dynamicmodel (with identifiable parameters) which fully integrates viral load, medi-cation adherence, drug resistance and baseline covariates data from an AIDSclinical trial study into the analysis. Thus, our dynamic model will be ableto characterize sustained suppression or resurgence of the virus as arising HUANG ET AL. from intrinsic viral dynamics, and/or influenced by factors such as drug sus-ceptibility and adherence during the treatment period of the clinical trial.The Bayesian nonlinear mixed-effects (BNLME) modeling approach [David-ian and Giltinan (1995)] is employed to estimate dynamic parameters andidentify significant clinical factors and/or covariates on virologic responseto ARV treatment. The rest of this article is organized as follows. Section2 introduces reparameterized viral dynamic models with time-varying drugefficacy which incorporates the effects of drug adherence, drug resistanceand baseline covariates, and briefly describes the BNLME modeling ap-proach, implemented via Markov chain Monte Carlo (MCMC) procedures,followed by defining the deviance information criterion (DIC) for compari-son of models. In Section 3 we summarize the motivating data set from anAIDS clinical trial study including the data of plasma viral load, medicationadherence from MEMS cap, drug resistance and baseline covariates; the pro-posed methodology is applied to these data and the results are presented.The method is evaluated via a simulation study in Section 4. Finally, weconclude the article with some discussions in Section 5.
2. HIV dynamic mechanism-based ODE models and statistical approaches.
This section aims to introduce long-term viral dynamic models based ona system of ODE with time-varying coefficients but without closed-formsolutions, and to investigate associated methodologies to demonstrate theapplication of these models to an AIDS clinical trial study. Long-term vi-ral dynamic models can be used to describe the interaction between cellssusceptible to target cells ( T ), infected cells ( T ∗ ) and free virus ( V ) byconsidering time-varying drug efficacy [Huang and Wu (2006); Huang et al.(2006)]. These three compartments (variables) are described as follows.HIV virions ( V ) will infect target cells ( T ) and turn them into infectedcells ( T ∗ ) at an infection rate k . Due to the intervention of antiviral drugs,we assume that drugs reduce the infection rate in the infected cells ( T ∗ ) by1 − γ ( t ) [0 ≤ γ ( t ) ≤ δ after producingan average of N virions per cell during their lifetimes, and free virions areremoved from the system at rate c . In addition to the dynamics describingvirus infection, we have to specify the dynamics of the uninfected cell pop-ulation. The simplest assumption is that uninfected cells are produced ata constant rate λ at which new T cells are generated from sources withinthe body, such as the thymus and die at a rate d T . Thus, the HIV dynamicmodel, after initiation of antiviral therapy, can be written as ddt T ( t ) = λ − d T T ( t ) − [1 − γ ( t )] kT ( t ) V ( t ) ,ddt T ∗ ( t ) = [1 − γ ( t )] kT ( t ) V ( t ) − δT ∗ ( t ) , (1) YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL ddt V ( t ) = N δT ∗ ( t ) − cV ( t ) , where the time-varying parameter γ ( t ) (as defined below) quantifies thetime-varying drug efficacy. If the regimen is not 100% effective [i.e., 0 ≤ γ ( t ) < γ ( t ) = γ (an unknownconstant), the model (1) becomes the model developed by Perelson andNelson (1999). In particular, when γ ( t ) = 0 (the drug has no effect), themodel (1) reduces to the model in the publications [Bonhoeffer et al. (1997);Nowak et al. (1995, 1997, 2000); Stafford et al. (2000)]; while γ ( t ) = 1 (thedrug is 100% effective), the model (1) reverts to the model discussed byNowak and May (2000) and Perelson and Nelson (1999).However, it is challenging to estimate all the parameters in the model (1)and to conduct inference because the model (1) is not a priori identifiable(i.e., multiple sets of parameters obtain identical fits to the data), givenonly viral load measurements [Cobelli et al. (1979)]. To obtain a modelwith a priori identifiable parameters [Labb´e and Verttoa (2006)], this paperinvestigates mechanism-based reparameterized ODE models to quantify thelong-term viral dynamics with ARV treatment and the associated statisticalmethods for model fitting.2.1. Reparameterized model with time-varying drug efficacy.
Followingthe studies [Perelson and Nelson (1999); Nowak and May (2000); Labb´e andVerttoa (2006)], we reparameterize the model (1) using the rescaled variables e T ( t ) = ( d T /λ ) T ( t ), e T ∗ ( t ) = ( δ/λ ) T ∗ ( t ), e V ( t ) = ( k/d T ) V ( t ). These yield therescaled version as follows: ddt e T ( t ) = d T λ ddt T = d T { − e T ( t ) − [1 − γ ( t )] e T ( t ) e V ( t ) } ,ddt e T ∗ ( t ) = δλ ddt T ∗ ( t ) = δ { [1 − γ ( t )] e T ( t ) e V ( t ) − e T ∗ ( t ) } , (2) ddt e V ( t ) = kd T ddt V ( t ) = c { R e T ∗ ( t ) − e V ( t ) } , where R = kN λ/ ( cd T ) represents the basic reproductive ratio for the virus,defined as the number of newly infected cells that arise from any one infectedcell when almost all cells are uninfected [Nowak and May (2000)]. Note thatthe rescaled model (2) has fewer parameters than the ‘original’ model (1).The identifiability of the model (2) is guaranteed [Cobelli et al. (1979);Labb´e and Verttoa (2006)] and parameters of the model can be uniquelyidentified. If R <
1, then the virus will not spread, since every infected cellwill on average produce less than one infected cell. If, on the other hand,
R >
1, then every infected cell will on average produce more than one newly
HUANG ET AL. infected cell and the virus will proliferate. For the HIV virus to persist inthe host, infected cells must produce at least one secondary infection, and R must be greater than unity [Nowak and May (2000)].Assuming steady state before the beginning of drug therapy, initial con-ditions for the model can now be expressed as simple functions of the initialconditions for viral load ( e V ): e T = 1 / (1 + e V ) , e T ∗ = e V / (1 + e V ) , e V = e V (0)[Cobelli et al. (1979); Labb´e and Verttoa (2006)]. The assumption of initialsteady state is necessary to guarantee identifiable (none of the models re-ported or referenced here is identifiable if the initial states are unknown),and is often justified by the clinical trial protocol. For example, in ACTG398, individual patients were taken off the drug before the initiation of thenew therapy (washout period to eliminate the effect of previously adminis-tered drugs and to guarantee that all individuals started from steady-stateconditions). Finally, viral load [ V ( t )] in model (1) is related to an equationoutput of viral load amount [ e V ( t )] in model (2) as follows: V ( t ) = ρ e V ( t ),where ρ , which is equivalent to a volume of distribution of pharmacokinet-ics, is a viral load scaling (proportionality) factor (10,000 copies/ml) to beestimated from the data [Nowak and May (2000)]. The set of ODE (2) willbe used to construct the BNLME model.2.2. Time-varying drug efficacy model.
Within the population of HIVvirions in a human host, there is likely to be genetic diversity and corre-sponding diversity in susceptibility to the various ARV agents. In clinicalpractice, genotypic or phenotypic tests can be performed to determine thesensitivity of HIV-1 RNA to ARV agents before a treatment regimen isselected. Here we use the phenotypic marker, IC [Molla et al. (1996)], toquantify agent-specific drug susceptibility. Because experimental data track-ing development of resistance suggest that the resistant fraction of the viralpopulation grows exponentially, we propose a model of log-linear functionto approximate the within-host changes over time in IC as follows: IC ( t ) = log (cid:18) S + S f − S t f t (cid:19) for 0 < t < t f ,log( S f ) for t ≥ t f , where S and S f are respective values of IC ( t ) at baseline and time point t f at which the resistant mutations dominate. In our study, t f is the timeof virologic failure which is observed from clinical studies. If S f = S , nonew drug resistant mutation is developed during treatment. Although morecomplicated models for median inhibitory concentration have been proposedbased on the frequencies of resistant mutations and cross-resistance patterns[Wainberg et al. (1996); Bonhoeffer, Lipsitch and Levin (1997)], in clinicalstudies or clinical practice it is common to collect IC values only at base-line and failure time as designed in ACTG 5055 [Acosta et al. (2004)] and YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL ACTG 398 [Hammer et al. (2002); Pfister et al. (2003)]. Thus, given that IC is only measured at baseline and at the time of treatment failure, thisfunction may serve as a good approximation in terms of data availability.Poor adherence to a treatment regimen is one of the major causes oftreatment failure [Ickovics and Meisler (1997)]. The following model is usedto represent adherence for a time interval T k < t ≤ T k +1 : A ( t ) = (cid:26) , if all doses are taken in ( T k , T k +1 ] ,r k , if 100 r k % doses are taken in ( T k , T k +1 ] , where 0 ≤ r k <
1, with r k indicating the adherence rate computed for eachassessment interval ( T k , T k +1 ] between study visits based on the question-naire or MEMS data; T k denotes the k th adherence assessment time.In most viral dynamic studies, investigators assumed that either drug effi-cacy was constant over treatment time [Perelson and Nelson (1999); Wu andDing (1999)] or antiviral regimens had perfect effect in blocking viral repli-cation [Ho et al. (1995); Perelson et al. (1996)]. However, the drug efficacymay change as concentrations of ARV drugs and other factors (e.g., drugresistance) vary during treatment. A simple pharmacodynamic sigmoidal E max model for dose–effect relationship is [Gabrielsson and Weiner (2000)] E = E max CEC + C , (3)where E max is the maximal effect that can be achieved, C is the drug concen-tration, and EC is the drug concentration that induced an effect equivalentto 50% of the maximal effect. Many different variations of the E max modelhave been developed by pharmacologists to model pharmacodynamic effects. E max models include the sigmoid E max model, the ordinary E max model andcomposite E max models [Gabrielsson and Weiner (2000); Davidian and Gilti-nan (1995)]. The ordinary E max model describes agonistic and antagonistic(inhibitory) effects of a drug, the sigmoid E max model is more flexible forthe steepness or curvature of the response–concentration curve compared tothe ordinary E max model, and composite E max models are used for multi-ple drug effects. More detailed discussions on E max models can be found inthe book by Gabrielsson and Weiner (2000) and the paper by Huang et al.(2003). Here we employ the following modified E max model to represent thetime-varying drug efficacy for two ARV agents within a class, γ ( t ) = A ( t ) / IC ( t ) + A ( t ) / IC ( t ) φ + A ( t ) / IC ( t ) + A ( t ) / IC ( t ) , (4)where A k ( t ) and IC k ( t ) ( k = 1 ,
2) are the adherence profile of the drug asmeasured by MEMS data and the time-course of median inhibitory con-centrations for the two drugs, respectively; φ = exp( β + β w + β w ); w HUANG ET AL. and w are observed baseline viral load and CD4 cell count, respectively; β = ( β , β , β ) T are unknown covariate effect parameters to be estimatedfrom clinical data. If β = β = 0 (without considering effect of covariates), φ = exp( β ) can be used to quantify the conversion between in vitro and invivo IC which is the case discussed by Huang et al. (2003). If γ ( t ) = 1,the drug is 100% effective, whereas if γ ( t ) = 0, the drug has no effect. Notethat, if A k ( t ), IC k ( t ), w and w are measured or obtained from a clinicalstudy and β can be estimated from clinical data, then the time-varying drugefficacy γ ( t ) can be estimated for the whole period of ARV treatment.2.3. Bayesian modeling approaches.
A number of studies investigatedvarious statistical methods, including Bayesian approaches, to fit viral dy-namic models and to predict virological responses [Han et al. (2002); Huanget al. (2006); Perelson et al. (1996); Wu et al. (1998); Wu and Ding (1999)].The Bayesian approach to viral dynamic modeling is particularly appealingfrom a biological perspective, as it allows informative prior distributions tobe incorporated. From a statistical estimation point of view, a Bayesian ap-proach is preferable because of the difficulties which are often encounteredfrom a classical approach when models involve the large numbers of param-eters, and complex nonlinearity of the subject-specific models. A Bayesiannonlinear mixed-effects (BNLME) model allows us to incorporate prior infor-mation at the population level into the estimates of dynamic parameters forindividual subjects. We briefly summarize the main concepts in the Bayesianapproach to inference and the presentation is, of course, far from exhaus-tive [Davidian and Giltinan (1995); Gelfand and Smith (1990); Huang et al.(2006); Wakefield et al. (1994); Wakefield (1996)].In reference to the model (2), we denote the number of subjects by n and the number of measurements on the i th subject by m i . Let µ =(log c, log δ, log d T , log ρ, log R, β , β , β ) T , Θ = { θ i , i = 1 , . . . , n } , θ i = (log c i , log δ i , log d T i , log ρ i , log R i , log φ i ) T and Y = { y ij , i = 1 , . . . , n ; j = 1 , . . . , m i } .Let f ij ( θ i , t j ) = log ( V ( θ i , t j )), where V ( θ i , t j ) is proportional to the nu-merical solution of e V ( t ) in the differential equations (2) for the i th subjectat time t j . Let y ij ( t ) and e i ( t j ) denote the repeated measurements of viralload in log scale and a measurement error with mean zero, respectively.Note that log-transformation of dynamic parameters and viral load is usedto make sure that estimates of dynamic parameters are positive and to sta-bilize the variance and convergence, respectively. The BNLME model canbe written in the following three levels [Gelfand and Smith (1990); Davidianand Giltinan (1995); Huang and Wu (2006); Wakefield (1996)]. Level
1. Within-subject variation: y i = f i ( θ i ) + e i , e i | σ , θ i ∼ N ( , σ I m i ) , (5) YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL where y i = ( y i ( t ) , . . . , y im i ( t m i )) T , f i ( θ i ) = ( f i ( θ i , t ) , . . . , f im i ( θ i , t m i )) T , e i = ( e i ( t ) , . . . , e i ( t m i )) T . Level
2. Between-subject variation: θ i = W i µ + b i , [ b i | Σ ] ∼ N ( , Σ ) , (6)where b i are random effects with mean zero. It is noteworthy that, for β = ( β , β , β ) T , no log-transformation is required as they are not nec-essarily positive. W i = ( I , J i , J i ), where I is an identity matrix and J si = (0 , , , , , w si ) T ( s = 1 , i = 1 , , . . . , n ) are 6 × w i and w i being (standardized) individual baseline viral load (in log scale)and CD4 cell count, respectively. For β = ( β , β , β ) T , we are only interestedin estimating them at population level. Thus, the individual parameter φ i is related to them as follows, log φ i = β + β w i + β w i + b i , where b i isthe last element of b i ( i = 1 , , . . . , n ). Level
3. Hyperprior distributions: σ − ∼ Ga ( a, b ) , µ ∼ N ( η , Λ ) , Σ − ∼ Wi ( Ω , ν ) , (7)where the mutually independent Gamma ( Ga ), Normal ( N ) and Wishart( Wi ) prior distributions are chosen to facilitate computations [Davidian andGiltinan (1995)]. The hyper-parameters a, b, η , Λ , Ω and ν can be determinedfrom previous studies and literature.The Bayesian approach is developed in the presence of observations whosevalue is initially uncertain and described through a probability distribution,which depends on some parameters. In the applications we assume thatthe researcher has some knowledge about at least some of the parameterswhich often represent characteristics of interest describing the process. TheBayesian approach incorporates this information through prior distributioninto observed data to obtain its posterior distribution. While computationof the posterior distribution involves solving multidimensional integrals, theintroduction of Markov chain Monte Carlo (MCMC) methods such as theGibbs sampler and Metropolis–Hastings algorithm opened the way to analy-sis of complex models through decomposition and sampling from full condi-tional distributions; see Gamerman (1997) and Gilks et al. (1995) for generaltheory and implementation details. Some more specific discussion of theBayesian dynamic modeling approach, including the choice of the hyper-parameters, the iterative MCMC algorithm and the implementation of theMCMC procedures can be found in publications by Huang and Wu (2006)and Wakefield (1996). The Bayesian approach was developed and tailoredas required by the unique features of the proposed HIV dynamic models.The basic principles of these proposed methodologies were well establishedin the statistical literature [Gamerman (1997); Gilks et al. (1995); Wakefield(1996)], but the applications of these methods in this paper are nonetheless HUANG ET AL. innovative within the context of a system of nonlinear ODE of time-varyingcoefficient, but without a closed-form solution.The progress in Bayesian posterior computation due to MCMC procedureshas made it possible to fit increasingly complex statistical models [Huangand Wu (2006); Wakefield (1996)] and entailed the wish to determine thebest fitting model in a class of candidates. Thus, it has become more andmore important to develop efficient model selection criteria. A recent publi-cation by Spiegelhalter et al. (2002) suggested a generalization of the Akaikeinformation criterion (AIC) [Akaike (1973)] and related also to the Bayesianinformation criterion (BIC) [Schwarz (1978)] that is deviance informationcriterion (DIC). In this paper we demonstrate its usefulness to compareBNLME models for longitudinal HIV dynamics discussed previously. Forcompleteness, a brief summary of DIC follows. More detailed discussion ofDIC and its properties can be found in publications by Spiegelhalter et al.(2002) and Zhu and Carlin (2000).Assume that the distribution of the data, Y , depends on the parametervector Ψ . Most recently, Spiegelhalter et al. (2002) suggested examining theposterior distribution of the deviance statistics defined by D ( Ψ ) = − p ( Y | Ψ ) + 2 log g ( Y )for Bayesian model comparison, where p ( Y | Ψ ) is the likelihood function,that is, the conditional joint probability density function of the observeddata Y given the parameter vector Ψ , and g ( Y ) denotes a fully specifiedstandardizing term that is a function of the data alone (which thus has noimpact on model selection). Based on the posterior distribution of D ( Ψ ),DIC consists of two components as follows: DIC = ¯ D + p D = 2 ¯ D − D ( ¯ Ψ ) , (8)where ¯ D = E Ψ | Y [ D ( Ψ )] = E Ψ | Y [ − p ( Y | Ψ )] and p D = ¯ D − D ( ¯ Ψ ) is theeffective number of parameters, defined as the difference between the poste-rior mean of the deviance and the deviance evaluated at the posterior mean¯ Ψ of the parameters. As with other model selection criteria, we cautionthat DIC is not intended for identification of the ‘correct’ model, but rathermerely as a method of comparing a collection of alternative formulations.In our model with different baseline characteristics and/or the MEMS ad-herence summary metrics, DIC can be used to identify the most significantcovariate and MEMS adherence summary metrics in contribution to viro-logic response. Under the model (5), in the absence of any standardizingfunction g ( Y ), the deviance is D ( σ − , Θ ) = − p ( Y | σ − , Θ )(9) = n X i =1 σ − ( y i − f i ( θ i )) T ( y i − f i ( θ i )) − log σ − n X i =1 m i . YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL As discussed above, our MCMC approach to estimating DIC first draws { ( σ − , Θ ) ( g ) } Gg =1 values from the posterior distribution, and then calculatescorresponding { D ( g ) } Gg =1 values from (9), where G is the number of samplesof posterior distribution. Finally, we estimate DIC as 2 ¯ D − D (¯ σ − , ¯ Θ ), where¯ D = G P Gg =1 D ( g ) , ¯ σ − = G P Gg =1 ( σ − ) ( g ) , ¯ θ i = G P Gg =1 θ ( g ) i and ¯ Θ = { ¯ θ i , i =1 , . . . , n } .
3. Analysis of AIDS clinical data.
Motivating application and observed data.
The subject sample inour analysis was drawn from the AIDS Clinical Trials Group (ACTG) 398study, a randomized, double-blind, placebo-controlled, 4-Arm trial study ofamprenavir (APV) as part of several dual protease inhibitor (PI) regimens insubjects with HIV infection in whom initial PI therapy had failed. Subjectsin all arms received APV (PI), three reverse transcriptase inhibitors (RTI):efavirenz (EFV), abacavir (ABC) and adefovir dipivoxil (ADV) plus a secondPI or placebo: Arm A (saquinavir = SQV), Arm B (indinavir = IDV), Arm C(nelfinavir = NFV) and Arm D (placebo matched for one of these three PIs).Subjects are HIV-infected individuals with prior exposure to approved PIsand who have exhibited loss of virologic suppression as reflected by a plasmaHIV-1 RNA concentration of ≥ IC ) and baselinecovariates to link plasma drug concentration to the long-term changes inHIV-1 RNA observation after initiation of therapy. In the model we incor-porate the two clinical factors (drug adherence measured by MEMS dataand drug susceptibility) and baseline viral load and CD4 cell count into afunction of treatment efficacy (see Section 2.2).Because phenotype sensitivity testing was performed only on a subset ofrandomly selected subjects, the number of subjects available for our analysiswas greatly reduced. We chose to consider only the subjects within Arm Cfor our analysis because this arm afforded the greatest number of subjects( n = 31) with available phenotypic drug susceptibility data on the two PIs(APV and NFV) and had available MEMS adherence data, as required for HUANG ET AL. our model. A summary of measurements of data to be used in our analysisis briefly described below.
Plasma viral load : Plasma viral load was measured in copies/ml at de-signed study time by the ultrasensitive reverse transcriptase–polymerasechain reaction HIV-1 RNA assay (Roche Molecular Systems). Only mea-surements taken while on protocol-defined treatment were used in the anal-ysis. The exact day of viral load measurement (not predefined study week)was used to compute study day in our analysis. A log transformation wasused in the analysis of viral load data. The graph in Figure 1 (up-left panel)shows the viral load trajectories of those subjects. Note that some viral loadmeasurements at designed study time were not observed due to laboratoryand other problems (for example, viral load measurement was not observedat week 12 for the subject displayed in Figure 1). Medication adherence : Medication adherence was measured by two meth-ods: by the use of questionnaires and by the use of MEMS [Pfister et al.(2003)]. Subjects completed an adherence questionnaire at study weeks. Thequestionnaire was completed by the study participant and/or by a face-to-face interview with study personnel. For MEMS, an MEMS cap was used tomonitor APV and EFV compliance only. The subjects were asked to bringtheir medication bottles and caps to the clinic at each study visit, where capdata were downloaded to computer files and stored for later analysis. TheMEMS adherence rate for APV was determined as the sum of positive dosingevents divided by the sum of prescribed dosing events during the specifiedtime interval. In our analysis, we assumed that NFV had the same MEMS
Fig. 1.
The profiles of viral load measurements (in log scale) from the 31 patients(up-left panel) and one trajectory of viral load (solid curve) and associated adherence rates(stairsteps) over time from the thirteen summary measures of MEMS data with the APVdrug for one subject from ACTG398. YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL adherence rate as APV since both APV and NFV were prescribed with thesame dosing schedule (twice daily), a prescribed AM and PM dosing periodwas defined for each subject and, hence, the bottles were opened twice perday [Pfister et al. (2003)]. As discussed previously, this study focuses mainlyon investigating optimal strategy to summarize adherence rates determinedby MEMS data for efficient production of virologic responses. For the MEMSdata analysis, it was not possible to model daily adherence rates and insteadthe adherence rate was computed with the following scenarios to considereffects of both interval length and time frame (delay of timing) for MEMSassessment.To determine the best summary metric of the MEMS adherence rate, weevaluated different assessment interval lengths (averaging adherence dosingevents over 1, 2 or 3 week intervals) and different assessment time frames(fixing the assessment interval times to end either immediately or 1, 2 or3 weeks prior to the next measured viral load). Table 1 summarizes theMEMS assessment interval notation and definitions for the 13 scenarios. Asan example, M2.2 in Table 1 denotes an MEMS adherence interval lengthof 2 weeks fixed to end 2 weeks prior to the next viral load measurement;for instance, the MEMS adherence rate for a subject at study week 8 (day56) was calculated as the number of nominal dosing events divided by thenumber of prescribed dosing events over study days 29–42 and this value wasused to represent adherence from the previous study visit to the study visitat the day 56 for modeling. The case M serves as a reference and averages allthe available MEMS data between viral load measurements. As an example, Table 1
Summary of the MEMS interval definitions and other information
Adherence interval definitionCase MEMS adherence Time frame length Interval Example for week 8interval name (weeks prior to viral load length (day 56), adherencemeasurement) computed over days HUANG ET AL.
Fig. 2.
The baseline ( ◦ ) and failure time ( × ) IC for APV/NFV drugs (upper panel)with baseline IC mean, standard deviation (SD) and coefficient of variation (CV), re-spectively, and the baseline viral load in log scale and baseline CD4 cell count (lowerpanel) with mean, SD and CV, respectively, for the 31 subject-specific individuals from theACTG398 study. Note that for a subject, if a single measurement of IC is observed atbaseline only, there is no × sign appearing in the upper panel of the plot. the viral load (in log scale) and adherence rates over time from the thirteencases of MEMS data with APV drug for the one representative subject arepresented in Figure 1. Phenotypic virus susceptibility to drug : The phenotypic virus resistance todrug were retrospectively determined from baseline samples. Patients wereselected to have samples assayed based on receiving study treatment forat least 8 weeks and having available sample. Some patients had virologicfailure and phenotypic susceptibility testing done on samples at the timeof failure. Testing was done via the recombinant virus assay (PhenoSense,ViroLogic Inc., South San Francisco, CA). For analysis, we used the phe-notype marker, IC [Molla et al. (1996)], to quantify agent-specific drugresistance. We refer to this marker as the median inhibitory concentration.The baseline ( ◦ ) and failure time ( × ) IC ’s of 31 subject-specific individu-als for the APV and NFV drugs are displayed in Figure 2 (upper panel) andare used to construct IC ( t ). Note that some subjects have only baseline IC due to the fact that they maintained viral suppression or dropped outfrom the study. If no IC measurement is observed at failure time for asubject, IC ( t ) becomes a constant in this case. YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL Baseline characteristics : The baseline viral load in log scale (VL) andthe baseline CD4 cell count were chosen as covariates in the model for dataanalysis. The log-transformation of viral load is used to stabilize the varianceof measurement error and estimation algorithm. The baseline characteristicsof 31 subject-specific individuals with mean, standard deviation (SD) andcoefficient of variation (CV) are displayed in Figure 2 (lower panel). To avoidvery small (large) estimates which may be unstable, we standardized thesecovariate values. For baseline log (RNA), for instance, each log (RNA)value is subtracted by mean (4.71) and divided by standard deviation (0.70).3.2. Model fitting and parameter estimation results.
In this section weapply the BNLME modeling approach to fit the data described in Section3.1. Based on the discussion in Section 2, the prior distribution for µ wasassumed to be N ( η , Λ ) with Λ being a diagonal matrix. Following the ideaof Huang and Wu (2006) for prior construction, as an example we discussthe prior construction for log δ . The prior constructions for other parametersare similar and so are omitted here.Ho et al. (1995) reported viral dynamic data on 20 patients; the logarithmof the average death rate of infected cells (log δ ) is − δ with − .
84 and − .
33, respectively.Following these two studies, Nowak et al. (1995) estimated log δ = − . δ from these studies are − . − . − .
33 and − . δ from thesestudies approximately follow a symmetric normal distribution. Thus, wechose a normal distribution N ( − .
0, 100.0) as the prior for log δ (the largevariance indicated that we used a noninformative prior for log δ ). Similarly,the values of the hyper-parameters at population level are chosen as follows[Ho et al. (1995); Nowak et al. (1995); Nowak and May (2000); Perelsonet al. (1996, 1997); Perelson and Nelson (1999); Verotta (2005); Wei et al.(1995)]: a = 4 . , b = 9 . , ν = 10 . , Λ = diag(100 . , . , . , . , . , . , . , . , η = (1 . , − . , − . , . , . , . , . , . T , Ω = diag(2 . , . , . , . , . , . . We decide that one long chain is run for MCMC implication with con-siderations of the following two issues: (i) a number of initial “burn-in”simulations are discarded, since from an arbitrary starting point it wouldbe unlikely that the initial simulations came from the stationary distribu-tion targeted by the Markov chain; (ii) one may only save every k th ( k HUANG ET AL. being an integer) simulation sample to reduce the dependence among sam-ples used for parameter estimation. Because the antiviral response modelinginvolves numerically solving nonlinear differential equations, thus compu-tational burdens would be more pronounced with the Bayesian approachesvia MCMC procedure. Utilizing efficient computer algorithms are critical inthis regard. Therefore, we are going to adopt these strategies in our MCMCimplementation using FORTRAN code that calls a differential equation sub-routine solver (DIVPRK) in the IMSL library (1994), which uses the Runge–Kutta–Verner fifth-order method. The computer codes are available fromthe corresponding author upon request. An informal check of convergenceis conducted based on graphical techniques according to the suggestion ofGelfand and Smith (1990). Based on the results, we propose that, after aninitial number of 20,000 burn-in iterations, every 4th MCMC sample wasretained from the next 80,000 samples. Thus, we obtained 20,000 samplesof targeted posterior distributions of the unknown parameters.We fitted the model to the data from 31 subjects discussed in Section3.1 using the proposed BNLME modeling approach. We incorporate thetwo clinical factors, drug adherence assessed by MEMS cap data and drugsusceptibility (phenotype IC values), as well as baseline covariates into afunction of drug efficacy. For model fitting, adherence rates were determinedfrom MEMS data with 13 different scenarios. For model fitting and thepurpose of comparisons, we set up a control model as the one without usingany drug adherence, resistance and covariate information which correspondsto setting γ ( t ) = 2 / (exp( β ) + 2) with IC ( t ) = 1, A ( t ) = 1 and w = w = 0;for this case, our model reverts to that discussed by Nowak and May (2000)and Perelson and Nelson (1999). The other 13 models are specified based onthe combination of drug resistance ( IC ), baseline covariate data and 13different adherence summary metrics listed in Table 1.In order to assess how adherence rates, determined from 13 different sce-narios, interacted with drug susceptibility and covariates to contribute to vi-rologic response, we fitted the models to all 13 scenarios as well as the controlmodel and compared the fitting results. We found based on the DIC criterion(see Figure 3) that, overall, the model with adherence rate determined fromMEMS dosing events, taken time frame length of 2 weeks prior to a viralload measurement with over either a 2 week assessment interval (M2.2) or a3 week assessment interval (M2.3), provided the best fits to the observeddata, compared to the other 12 models for most subjects. The referencemodel with adherence rate averaged by all the available MEMS data be-tween viral load measurements gave a moderate fit to the observed data.We clearly see that all models fit the early viral load data well, but thecontrol model, lacking factors for subject-specific drug adherence and sus-ceptibility as well as baseline covariates, failed to fit viral load reboundsand fluctuations, and provided the worst fitting results for the majority of YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL Fig. 3.
Comparison of the DIC values for the models from 13 different determinants ofMEMS adherence, interacted by drug resistance and covariates, with the control model.The two horizontal lines represent the DIC values for the control model and the referencemodel with adherence rate determined by case M, respectively. subjects. For the purpose of illustration, the model fitting curves from thecontrol model (solid curves), the best fit model (M2.2: dotted curves) andthe reference model (M: dashed curves) are displayed in Figure 4 for thethree representative subjects.For the purpose of comparison, Figure 5 presented the population posteriormeans and the corresponding 95% equal-tail credible intervals (CI) of the eightparameters for the control model, the best fit model (M2.2) and the referencemodel (M). For the six dynamic parameters ( c, δ, d T , ρ, R, β ), it is shownthat the population estimates for the control model have higher clearance Fig. 4.
The estimate of viral load trajectory from the model fitting with the 3 different de-terminants of adherence: (i)
Control model (solid curves), (ii)
M2.2 model (dotted curves)and (iii) reference (M) model (dashed curves) for the three representative subjects. Theobserved values are indicated by circles. HUANG ET AL.
Fig. 5.
A summary of the estimated posterior means ( ◦ ) of population parameters andthe corresponding equal-tail credible intervals (CI) for the models from 3 differentdeterminants of adherence. rate of free virions ( c ), lower death rate of infected cells ( δ ), higher death rateof target T cells ( d T ), smaller viral load scaling factor ( ρ ), higher basic repro-ductive ratio for the virus ( R ) and larger φ than those for the best fit model(M2.2) and reference model (M), while the population estimates for the bestfit model and reference model are generally similar. These differences mayresult from the effects of drug adherence interacted with drug resistanceand covariates in the models. For the other two covariate effect parame-ters ( β , β ) which are relevant to treatment effect, we will discuss themseparately in Section 3.4. In terms of the individual parameter estimates, alarge between-subject variation in the estimates of all individual dynamicparameters was observed (data not shown here). Overall, the coefficient ofvariation ranges from 15.4% to 88.9% for all parameters.3.3. Effects of adherence rates determined by different MEMS summarymetrics.
Figure 3 in Section 3.2 displayed a comparison of the DIC val-ues for the models from 13 different determinants of MEMS adherence, in-teracted by drug resistance and covariates, with the control model. Theobserved patterns shown in Figure 3 provided information to answer thefollowing questions: (i) what MEMS assessment interval length is best and(ii) what MEMS assessment time frame (delay effect of timing) is best?
YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL We can see that when the time frame for MEMS assessment is fixed, mod-els with a 2 week MEMS assessment interval length generally outperformmodels with an assessment interval length of 1 or 3 weeks except for thetime frame length with 3 weeks prior to viral load measurement where themodel with a 3 week MEMS assessment interval length performs best.Regardless of the assessment interval length, models which assess compli-ance 2 weeks prior to viral load measurement generally outperform modelswhich assess compliance immediately before viral load measurement, 1 weekbefore or 3 weeks before viral load measurement. Overall, the model with aMEMS assessment interval length of 2 weeks measured from 4 to 2 weeksprior to viral load measurement (M2.2) was significantly a better predicatorof viral load over time than any other models, with the exception of theM2.3 model which shows no significant difference from the M2.2 model interms of DIC values.3.4.
Treatment effects of baseline characteristics interacted with clinicalfactors.
Figure 5 summarized the population posterior means and the cor-responding 95% equal-tail CI of the covariate effect parameters β and β for the best fit model (M2.2) and the reference model (M). It can be seenthat estimates of β (coefficient of baseline viral load) are negative, whileestimates of β (coefficient of baseline CD4 cell count) are positive. In fact,other models also provided the same scenarios for the estimates of these twoparameters (not shown here). As an example, we report results based on thebest fit model (M2.2). We can observe from Figure 5 that the estimates of β and β are ˆ β = − .
67 with 95% CI ( − − β = 0 .
719 with95% CI (0.371, 1.058). It indicates that, according to antiviral drug efficacymodel (4), the baseline viral load ( ˆ β = − .
67) has a significant positiveeffect on drug efficacy γ ( t ), while the baseline CD4 cell count ( ˆ β = 0 . γ ( t ) since the corresponding 95% credibleintervals do not contain zero for both parameters. These findings could sug-gest us with the following different ways. The lowest value of γ ( t ) [highest φ as displayed in Figure 6(b)] occurs in the subjects with the best prognosis(higher baseline CD4 cell count and lower baseline viral load). Alternatively,the highest value of γ ( t ) (lowest φ ) occurs in those with the worst prognosis(lower baseline CD4 cell count and higher baseline viral load). A possibleexplanation is that there is a floor effect of viral load (or ceiling/floor effectof CD4 cell count) that is not captured in the model. Further, given thatbaseline CD4 cell count and viral load are jointly used to make treatmentdecisions and are known to be negatively correlated as shown in Figure 6(a),the result based on the combination of baseline viral load and CD4 countin γ ( t ) indicates that the baseline CD4 cell count and viral load have theopposite effect on drug efficacy which might be intuitively understandable. HUANG ET AL.
Fig. 6. (a)
Correlation between (standardized) baseline log (RNA) and CD4 cell count.The correlation coefficient ( r ) and p-value are obtained from the Spearman rank correlationtest. The line is a robust (MM-estimator) linear regression fit. (b) Estimated φ values forthe 31 subject-specific individuals.
4. A simulation study.
In this paper we investigated the association be-tween virologic outcomes and medication adherence with confounding fac-tors based on the data from 31 subjects. As both one referee and the asso-ciate editor suggested that a simulation study may be useful to evaluate howour method performs, in this section we conduct a limited simulation studyhere due to intensive computations involved. The scenario we consider is asfollows.We simulate a clinical trial with 31 patients receiving antiviral treatment.For each patient, we assume that the designs of this experiment, in particu-lar, the sampling times for viral load, were the same as those in the ACTG398 study. The data for the phenotype marker (baseline and failure IC ’s),medication adherence and the baseline viral load/CD4 cell count were takenfrom the ACTG 398 study, where medication adherence was calculated bythe M2.2 summary measure. The “true” values of unknown parameters werethe same as those estimated from the data set of 31 subjects which were re-ported in Section 3. With generated individual true parameters based on theequation (6), known data [ IC ( t ) , A ( t ) , w and w ], we generated randomsamples for response (viral load) based on model (2). The values of hyper-parameters are chosen to be the same as those in Section 3. For each simu-lated data set, we fit the model using the Bayesian approach. The MCMCtechniques consisting of a series of Gibbs sampling and M–H algorithms werethe same as those in the real data analysis. We performed 50 replicationsand obtained the mean estimates (ME) of population parameters togetherwith the corresponding relative bias (RB), which is the difference betweenthe mean estimate and the true value of the parameter divided by absolutevalue of the true parameter, and the standard error (SE), defined as the YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL Table 2
The true values (TV) of parameters and mean estimates (ME) of population parameterswith 50 replications as well as the corresponding relative bias (RB), defined as × ( ME − TV ) / | TV | , and standard error (SE), defined as × √ MSE / | TV | Parameter TV ME RB ( % ) SE ( % ) log c .
767 0 .
771 0 .
522 9 . δ − . − . − .
685 7 . d T − . − . − .
367 5 . ρ .
433 0 . − .
39 18 . R .
040 1 .
100 5 .
769 3 . β − . − .
604 0 .
421 4 . β − . − .
665 0 .
746 10 . β .
719 0 . − .
921 13 . square root of mean-squared error divided by the absolute value of the trueparameter.In Table 2 we summarize the true values (TV) of parameters and the MEof population parameters with 50 replications as well as the correspond-ing RB and SE. The percentage is based on the absolute value of the trueparameter. It can be seen from Table 2 that the RB (%) for populationparameter estimates are very small, ranging from 0.522 to 10.39, and theSE (%) ranges from 3.089 to 18.03. The simulation results indicate that ourmethod with considering the M2.2 model performs reasonably well in termsof estimates of parameters except for the viral load proportionality factorlog ρ which has larger RB and SE. That is, our method produces a substan-tially biased estimate and may severely underestimate log ρ . This may beexplained by the fact that it is probably caused from inaccurate numericalsolutions to the system of ODE (2) which was used to construct the BNLMEmodel.
5. Concluding discussion.
In developing long-term dynamic modeling,this paper introduced a dynamic mechanism specified by a system of time-varying ODE to (i) establish a link between success of ARV therapy invirologic response and MEMS adherence confounded by drug resistance andbaseline covariates, (ii) fully integrate viral load, MEMS adherence, drug re-sistance and baseline covariates data into the statistical inference and analy-sis, and (iii) provide a powerful tool to evaluate the effects of MEMS adher-ence determined by a different summary metric on virologic response usingthe BNLME modeling approach. This approach cannot only combine priorinformation with current clinical data for estimating dynamic parameters,but also deal with complex dynamic systems. Thus, the results of estimateddynamic parameters based on this model should be more reliable and rea-sonable to interpret long-term HIV dynamics. Our models are simplified HUANG ET AL. with the main goals of retaining crucial features of HIV dynamics and, atthe same time, guaranteeing their applicability to typical clinical data, inparticular, long-term viral load measurements. The proposed model fittedthe clinical data reasonably well for most patients in our study, although thefitting for a few patients was not completely satisfactory because of unusualviral load fluctuation patterns for these subjects.We have explored the practical performance of DIC for the comparisonof developed models. DIC, a Bayesian version of the classical deviance formodel assessment, is particularly suited to compare Bayesian models whoseposterior distribution has been obtained using MCMC procedures and canbe used in complex hierarchical models where the number of unknowns oftenexceeds the number of observations and the number of free parameters is notwell defined. This is in contrast to AIC and BIC, where the number of freeparameters needs to be specified [Zhu and Carlin (2000)]. Overall, combinedwith more traditional residual analysis and posterior predictive model checksas discussed in this paper, DIC appears to offer a comprehensive frameworkfor comparison and evaluation within a complex model class.Several studies investigated the association between virologic responsesand adherence assessed by MEMS data only without considering other con-founding factors such as drug resistance using standard modeling methodsincluding Poisson regression [Knafl et al. (2004)], logistic regression [Vrijenset al. (2005)] and the linear mixed-effects model [Liu et al. (2007)]. In thisarticle we employed the proposed dynamic model and associated BNLMEmodeling approach to assessment of effects of adherence determinants basedon MEMS dosing events in predicting virologic response. In particular, weinvestigated (i) how to summarize the MEMS adherence data for efficientprediction of virological response after accounting for potential confound-ing factors such as drug resistance and baseline covariates, and (ii) how toevaluate treatment effect of baseline characteristics interacted with MEMSadherence and other clinical factors. Note that a further study in comparingthe performance of these different methods may be important and war-ranted, although some challenges are observed in terms of different modelstructures and data characteristics.The results indicate that the best summary metric for prediction of viro-logic response based on DIC criterion is the adherence rate determined byMEMS dosing events averaged over an assessment interval of 2 or 3 weeks,and 2 weeks prior to the next measured viral load observation (denoted byM2.2 or M2.3). We found that the best MEMS adherence predictor (M2.2)of the effectiveness of ARV medications on virologic response is consistentwith that reported in Huang et al. (2008) in which, however, the next bestMEMS adherence predictor (M1.2) is different from what is obtained in thispaper. This difference may be due to the various reasons as follows. In thestudy by Huang et al. (2008), (i) it directly applied the model (1) to fit data
YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL and, thus, some assumptions were made due to parameter unidentifiable is-sues; (ii) the analysis used the mean of the sum of squared deviations as acriterion to evaluate model fitting results; (iii) it assumed IC data wereextrapolated linearly to the whole treatment period instead of a log-linearextrapolation offered in this paper which is considered more reasonable bi-ologically; and (iv) it did not incorporate baseline covariates in the model.In addition, the superiority of the M2.2 model, associated with the MEMSadherence rate based on time frame length of 2 weeks prior to a viral loadmeasurement with over a 2 week assessment interval, may be explained bythe fact that it probably reflects how long it takes for resistance mutationsto first arise and then come to dominate the plasma population of a virus. Aspointed out by an anonymous referee, this finding may also be interpretedas follows. Low adherence two weeks prior to the viral load measurementmay not have had sufficient time for viral rebound to occur.In this paper we set up a connection between subject-specific baselinecharacteristics with interaction of clinical factors and drug efficacy. We alsofound that, according to antiviral drug efficacy model (4), the baseline viralload had a positive effect on drug efficacy, while the baseline CD4 cell counthad a negative effect on it. Our results may be explained by the fact thatfor those patients with higher baseline viral load, the drug efficacy needsto be higher than that for those with lower baseline viral load. Therefore,a strong treatment is recommended for those patients with higher baselineviral load. On the other hand, patients with higher CD4 cell count may needlower drug efficacy so that a more potent ARV drug regimen is not necessaryfor these patients to avoid side-effects of drug use. The results may suggestthe benefit of initiating ARV therapy with a lower baseline viral load and/ora higher baseline CD4 count. These results coincide with those investigatedby Notermans et al. (1998) and Wu et al. (2005) whose results were ob-tained using correlation analysis. Note that given the estimated parameters,the subject with both a high baseline CD4 cell count and a relatively highbaseline viral load [upper right quadrant of Figure 6(a)] has a very different φ than that with a similar baseline CD4 cell count, but a low baseline viralload [upper left quadrant of Figure 6(a)]. It is possible that the subject inthe upper right quadrant was more recently infected (hence the higher base-line CD4 cell count) or perhaps with a drug resistant virus and would notbe a candidate for a regimen with a “less potent drug efficacy.”Our findings need to be interpreted in light of the study limitations. First,in the ACTG 398 study, because phenotype sensitivity testing was performedonly on a subset of randomly selected subjects, we chose 31 patients whohave available data for analysis in this paper. Second, due to reasons such aslost caps and malfunction of caps, there were inaccurate MEMS data acrossthe treatment period which may not reflect actual adherence profile for in-dividual patients and, thus, the data quality could have some impact on the HUANG ET AL. results. Third, because of technical limitations, the undetectable values of vi-ral load were replaced with 25 copies/ml for analyses, which could introducesome bias due to a cluster of ties of data points. Finally, this paper combinednew technologies in mathematical modeling and statistical inference with ad-vances in HIV/AIDS dynamics and ARV therapy to quantify complex HIVdisease mechanisms. The complex nature of HIV/AIDS ARV therapy willnaturally pose some challenges including missing data and measurement er-ror in clinical factors and covariates. These complicated problems, which arebeyond the scope of this article, may be addressed, for example, using thejoint model method [Wu (2002)] and other techniques [Carroll et al. (1995)],and are warranted for further investigation. Nevertheless, these limitationswould not offset the major findings from this study.As the Associate Editor pointed out, we assumed that the distributionsof the random error and random effects are normal, which is a common as-sumption in the literature for statistical inference. However, due to the na-ture of AIDS clinical data, it is possible that the data may contain outliersand/or depart from normality and, thus, statistical inference and analysiswith normal assumption may lead to misleading results [Verbeke and Lesaf-fre (1996); Ghosh et al. (2007)]. Specially non-normal characteristics such asskewness with heavy right or left tail may appear often in virologic responses.Thus, a normality assumption may be too restrictive to provide an accuraterepresentation of the structure that is often present in repeated measuresand clustered data. Thus, it is of practical interest to investigate nonlinearmodels with a skew-normal distribution or t distribution for (within-subject)random error and random effects which are more robust to outliers and skew-ness than those with a normal distribution. In our recent study [Huang andDagne (2010)] we addressed a Bayesian approach to nonlinear mixed-effectsmodels in conjunction with the HIV dynamic model and relaxed the nor-mality assumption by considering both random error and random-effects tohave a multivariate skew-normal distribution. The proposed model providesflexibility in capturing a broad range of non-normal behavior and includesnormality as a special case. The results suggest that it is very important toassume models with a skew-normal distribution in order to achieve robustand reliable results, in particular, when the data exhibit skewness. We areactively applying this methodology into the data investigated in this paperand will report the results in a future study.In summary, the mechanism-based dynamic model is powerful and effi-cient to characterize relations between antiviral response and medicationadherence, drug susceptibility as well as baseline characteristics, althoughsome biological assumptions are required. It is important to find a way toincorporate subject-specific information with regard to drug susceptibility,medication adherence and baseline characteristics in predicting long-termvirologic response. Since each of these factors may only contribute a very YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL small portion to virologic response and they may be confounded throughcomplicated interactions, the appropriate modeling of the combination ef-fects of these factors is critical to efficiently utilize the information in viro-logic response predictions. The viral dynamic model and associated statis-tical approaches discussed here provide a good avenue to fulfill this goal. Inparticular, MEMS adherence rate summarized by an optimal way in termsof assessing both interval lengths and time frame lengths prior to viral loadmeasurement is an important factor that significantly determines the effec-tiveness of ARV treatment and needs to be taken into account in analy-sis of virologic responses. Our results demonstrate that MEMS adherencedata may not predict virologic response well unless the MEMS cap data aresummarized in an appropriate way as reported in Section 3.3. Additionally,although this paper concentrated on HIV dynamics, the basic concept oflongitudinal dynamic systems and the proposed methodologies in this paperare generally applicable to dynamic systems in other fields such as biology,medicine, engineering or PK/PD studies as long as they meet the relevanttechnical specification—a system of ODE. Acknowledgments.
The authors are extremely grateful to the Editor,an Associate Editor and one referee for their insightful comments and con-structive suggestions that led to an improvement of the article. We grate-fully acknowledge ACTG 398 study investigators for allowing us to use theclinical data from their study. The authors are indebted to Dr. Susan L.Rosenkranz from Frontier Science & Technology Research Foundation andSDAC of Harvard School of Public Health for her informative discussionsand data preparations. REFERENCES
Acosta, E. P., Wu, H., Walawander, A., Eron, J., Pettinelli, C., Yu, S., Neath,D., Ferguson, E., Saah, A. J., Kuritzkes, D. R., Gerber, J. G., for the AdultACTG 5055 Protocol Team (2004). Comparison of two indinavir/ritonavir regi-mens in treatment-experienced HIV-infected individuals.
Journal of Acquired ImmuneDeficiency Syndromes Akaike, H. (1973). Information theory and an extension of the maximum likelihoodprinciple. In (B. N. Petrov andF. Cs´aki, eds.) 267–281. Akad´emiai Kiad´o, Budapest. MR0483125
Arnsten, J. H., Demas, P. A., Farzadegan, H. et al. (2001). Antiretroviral therapyadherence and viral suppression in HIV-infected drug users: Comparison of self-reportand electronic monitoring.
Clin. Infect. Dis. Bonhoeffer, S., Lipsitch, M. and
Levin, B. R. (1997). Evaluating treatment protocolsto prevent antibiotic resistance.
Proc. Natl. Acad. Sci. USA Bonhoeffer, S., May, R. et al. (1997). Viral dynamics and drug therapy.
Proc. Natl.Acad. Sci. USA Bova, C. A., Fennie, K. P., Knafl, G. J. et al. (2005). Use of electronic monitoringdevices to measure antiretroviral adherence: Practical considerations.
AIDS Behav. HUANG ET AL.
Carroll, R. J., Ruppert, D. and
Stefanski, L. A. (1995).
Measurement Error inNonlinear Models . Chapman and Hall, London. MR1630517
Cobelli, C., Lepschy, A. and
Jacur, G. R. (1979). Identifiability of compartmentalsystems and related structural properties.
Math. Biosci. Davidian, M. and
Giltinan, D. M. (1995).
Nonlinear Models for Repeated MeasurementData . Chapman and Hall, London.
Gabrielsson, J. and
Weiner, D. (2000).
Pharmacokinetic and Pharmacodynamic DataAnalysis: Concepts and Applications . Apotekarsocieteten, Stockholm.
Gamerman, D. (1997).
Markov Chain Monte Carlo: Stochastic Simulation for BayesianInference . Chapman and Hall, London. MR2260716
Gelfand, A. E. and
Smith, A. F. M. (1990). Sampling-based approaches to calculatingmarginal densities.
J. Amer. Statist. Assoc. Ghosh, P., Branco, M. D. and
Chakraborty, H. (2007). Bivariate random effectmodel using skew normal distribution with application to HIV-RNA.
Statist. Med. Gilks, W. R., Richardson, S. and
Spiegelhalter, D. J., eds. (1995).
Markov ChainMonte Carlo in Practice . Chapman and Hall, London. MR1397966
Hammer, S. M., Vaida, F. et al. (2002). Dual vs single protease inhibitor therapyfollowing antiretroviral treatment failure: A random trial.
J. Amer. Med. Assoc.
Han, C., Chaloner, K. and
Perelson, A. S. (2002). Bayesian analysis of a populationHIV dynamic model. In
Case Studies in Bayesian Statistics (C. Gatsoiquiry, R. E.Kass, A. Carriquiry, A. Gelman, D. Higdon, D. K. Pauler and I. Verdinellinis, eds.) Haubrich, R. H., Little, S. J., Currier, J. S. et al. (1999). The value of patient-reported adherence to antiretroviral therapy in predicting virologic and immunologogicresponse.
AIDS Ho, D. D., Neumann, A. U., Perelson, A. S., Chen, W., Leonard, J. M. and
Markowitz, M. (1995). Rapid turnover of plasma virions and CD4 lymphocytes inHIV-1 infection.
Nature
Huang, Y. and
Dagne, G. (2010). Skew-normal Bayesian nonlinear mixed-effects modelswith application to AIDS studies.
Statist. Med. Huang, Y., Holden-Wiltse, J., Rosenkranz, S. L., Eron, J. J., Hammer, S. M.,Mellors, J. W. and
Wu, H. (2008). Modeling adherence of protease inhibitors forprediction of virologic response in HIV-1 infected patients: Comparison between ad-herence measurements by MEMS and questionnaire. Technical report, Department ofEpidemiology and Biostatistics, Univ. South Florida, Tampa, FL.
Huang, Y., Liu, D. and
Wu, H. (2006). Hierarchical Bayesian methods for estimation ofparameters in a longitudinal HIV dynamic system.
Biometrics Huang, Y., Rosenkranz, S. L. and
Wu, H. (2003). Modeling HIV dynamics and an-tiviral responses with consideration of time-varying drug exposures, sensitivities andadherence.
Math. Biosci.
Huang, Y. and
Wu, H. (2006). A Bayesian approach for estimating antiviral efficacy inHIV dynamic model.
J. Appl. Statist. Ickovics, J. R. and
Meisler, A. W. (1997). Adherence in AIDS clinical trial: A frame-work for clinical research and clinical care.
J. Clin. Epid. IMSL MATH/LIBRARY (1994).
FORTRAN Subroutines for Mathematical Applications ,Vol. 2. Visual Numerics, Houston.YNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL Kastrissios, H., Suarez, J. R., Katzenstein, D. et al. (1998). Characterizing pat-terns of drug-taking behavior with a multiple drug regimen in an AIDS clinical trial.
AIDS Knafl, G. J., Fennie, K. P., Bova, C. et al. (2004). Electronic monitoring device eventmodeling on an individual-subject basis using adaptive Poisson regression.
Statist. Med. Labb´e, L. and
Verotta, D. (2006). A nonlinear mixed effect dynamic model incorporat-ing prior exposure and adherence to treatment to describe long-term therapy outcomein HIV-patients.
J. Pharmacokinet. Pharmacodyn. Levine, A. J., Hinkin, C. H., Marion, S. et al. (2006). Adherence to antiretroviralmedications in HIV: Differences in data collected via self-report and electronic moni-toring.
Health Psychol. Liu, H., Miller, L. G., Golin, C. E. et al. (2007). Repeated measures analyses ofdose timing of antiretroviral medication and its relationship to HIV virologic outcomes.
Statist. Med. Molla, A., Korneyeva, M. et al. (1996). Ordered accumulation of mutations in HIVprotease confers resistance to ritonavir.
Nat. Medic. Notermans, D. W., Goudsmit, J., Danner, S. A. et al. (1998). Rate of HIV-1 declinefollowing antiretroviral therapy is related to viral load at baseline and drug regimen.
AIDS Nowak, M. A., Bonhoeffer, S., Clive, L., Balfe, P., Semple, M., Kaye, S.,Tenant-Flowers, M. and
Tedder, R. (1995). HIV results in the frame.
Nature
Nowak, M. A., Lloyd, A. et al. (1997). Viral dynamics of primary viremia and an-tiretroviral therapy in simian immunodeficiency virus infection.
J. Virol. Nowak, M. A. and
May, R. M. (2000).
Virus Dynamics: Mathematical Principles ofImmunology and Virology . Oxford Univ. Press, Oxford. MR2009143
Paterson, D. L., Swindells, S., Mohr, J., Brester, M., Vergis, E. N., Squler,C., Wagener, M. M. and
Singh, N. (2000). Adherence to protease inhibitor therapyand outcomes in patients with HIV infection.
J. Internal Med.
Perelson, A. S., Neumann, A. U., Markowitz, M., Leonard, J. M. and
Ho, D. D. (1996). HIV-1 dynamics in vivo: Virion clearance rate, infected cell life-span, and viralgeneration time.
Science
Perelson, A. S., Essunger, P., Cao, Y., Vesanen, M., Hurley, A., Saksela, K.,Markowitz, M. and
Ho, D. D. (1997). Decay characteristics of HIV-1-infected com-partments during combination therapy.
Nature
Perelson, A. S. and
Nelson, P. W. (1999). Mathematical analysis of HIV-1 dynamicsin vivo.
SIAM Rev. Pfister, M., Labbe, L., Hammer, S. M., Mellors, J., Bennett, K. K., Rosenkranz,S. L. and
Sheiner, L B. and AIDS
Clinical Trial Group Protocol 398 Inves-tigators (2003). Population Pharmacokinetics and Pharmacodynamics of Efavirenz,Nelfinavir, and Indinavir: Adult AIDS Clinical Trial Group Study 398.
AntimicrobialAgents and Chemoterapy Schwarz, G. (1978). Estimating the dimension of a model.
Ann. Statist. Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and van der Linde, A. (2002).Bayesian measures of model complexity and fit.
J. Roy. Statist. Soc. Ser. B Stafford, M. A. et al. (2000). Modeling plasma virus concentration during primaryHIV infection.
J. Theoret. Biol. HUANG ET AL.
Verbeke, G. and
Lesaffre, E. (1996). A linear mixed-effects model with heterogeneityin random-effects population.
J. Amer. Statist. Assoc. Verotta, D. (2005). Models and estimation methods for clinical HIV-1 data.
J. Comput.Appl. Math.
Vrijens, B., Goetghebeur, E., de Klerk, E., Rode, R., Mayer, S. and
Urquhart,J. (2005). Modeling the association between adherence and viral load in HIV-infectedpatients.
Statist. Med. Wainberg, M. A. et al. (1996). Effectiveness of 3TC in HIV clinical trials may be duein part to the M184V substition in 3TC-resistant HIV-1 reverse transcriptase.
AIDS (Suppl.), S3–S10. Wakefield, J. C., Smith, A. F. M., Racine-Poon, A. and
Gelfand, A. E. (1994).Bayesian analysis of linear and non-linear population models using the Gibbs sampler.
Appl. Statist. Wakefield, J. C. (1996). The Bayesian analysis to population pharmacokinetic models.
J. Amer. Statist. Assoc. Wei, X., Ghosh, S. K., Taylor, M. E., Johnson, V. A., Emini, E. A., Deutsch,P., Lifson, J. D., Bonhoeffer, S., Nowak, M. A., Hahn, B. H., Saag, M. S. and
Shaw, G. M. (1995). Viral dynamics in human immunodeficiency virus type 1infection.
Nature
Wu, H., Ding, A. A. and de Gruttola, V. (1998). Estimation of HIV dynamic param-eters.
Statist. Med. Wu, H. and
Ding, A. A. (1999). Population HIV-1 dynamics in vivo: Applicable modelsand inferential tools for virological data from AIDS clinical trials.
Biometrics Wu, H., Huang, Y. et al. (2005). Modeling antiretroviral response: Effects of drugpotency, pharmacokinetics, adherence and drug resistance.
Journal of Acquired ImmuneDeficiency Syndromes Wu, L. (2002). A joint model for nonlinear mixed-effects models with censoring andcovariates measured with error, with application to AIDS studies.
J. Amer. Statist.Assoc. Zhu, L. and
Carlin, B. P. (2000). Comparing hierarchical models for spatio-temporallymisaligned data using the deviance information criterion.
Stat. Med. Y. HuangDepartment of Epidemiologyand BiostatisticsCollege of Public Health, MDC 56University of South FloridaTampa, Florida 33612USAE-mail: [email protected]
H. WuJ. Holden-WiltseDepartment of Biostatisticsand Computational BiologySchool of Medicine and DentistryUniversity of RochesterRochester, New York 14642USAE-mail: [email protected] [email protected]