Extension of the SAEM algorithm for nonlinear mixed models with two levels of random effects
This article focuses on parameter estimation of multi-levels nonlinear mixed effects models (MNLMEMs). These models are used to analyze data presenting multiple hierarchical levels of grouping (cluster data, clinical trials with several observation p…
Authors: Xavi`ere Panhard, Adeline Samson (MAP5)
Extension of the SAEM algorithm for nonlinear mixed mo dels with t w o lev els of random eets Xa vière P anhard 1 ∗ A deline Samson 2 1 INSERM U738, P aris, F rane; Univ ersit y P aris 7, UFR de Médeine, P aris, F rane 2 Lab oratoire MAP5, UMR CNRS 8145, Univ ersit y P aris 5, F rane * Corresp onding author: INSERM U738, Univ ersité P aris 7, UFR de Médeine, 16 rue Hu hard, 75018 P aris, F rane Abstrat This artile fo uses on parameter estimation of m ulti-lev els nonlinear mixed eets mo d- els (MNLMEMs). These mo dels are used to analyze data presen ting m ultiple hierar hial lev els of grouping (luster data, linial trials with sev eral observ ation p erio ds,...). The v ariabilit y of the individual parameters of the regression funtion is th us deomp osed as a b et w een-sub jet v ariabilit y and higher lev els of v ariabilit y (for example within-sub jet v ariabilit y). W e prop ose maxim um lik eliho o d estimates of parameters of those MNLMEMs with t w o lev els of random eets, using an extension of the SAEM-MCMC algorithm. The extended SAEM algorithm is split in to an expliit diret EM algorithm and a sto hasti EM part. Compared to the original algorithm, additional suien t statistis ha v e to b e appro ximated b y relying on the onditional distribution of the seond lev el of random ef- fets. This estimation metho d is ev aluated on pharmaokineti ross-o v er sim ulated trials, mimi king theoph yllin onen tration data. Results obtained on those datasets with either the SAEM algorithm or the F OCE algorithm (implemen ted in the nlme funtion of R soft- w are) are ompared: biases and RMSEs of almost all the SAEM estimates are smaller than the F OCE ones. Finally , w e apply the extended SAEM algorithm to analyze the pharma- okineti in teration of tenofo vir on atazana vir, a no v el protease inhibitor, from the ANRS 107-Puzzle 2 study . A signian t derease of the area under the urv e of atazana vir is found in patien ts reeiving b oth treatmen ts. 1 KEYW ORDS: Multilev el nonlinear mixed eets mo dels; SAEM algorithm; Multiple p erio ds; Cross-o v er trial; Bio equiv alene trials. 1 In tro dution The use of non-linear mixed eets mo dels (NLMEMs) inreases in sev eral elds su h as agron- om y , forestry , linial trials, p opulation pharmaokinetis (PK) and pharmao dynamis or viral dynamis to mo del longitudinal data. In some settings, data an presen t m ultiple hierar hial lev els of grouping, leading to m ultiple nested lev els of v ariabilit y . F or instane, w e ma y study patien ts that are group ed in medial servies that are themselv es group ed in to hospitals. In this artile, w e all m ultilev el non-linear mixed eets mo dels (MNLMEMs) the mo dels that desrib e su h data. MNLMEMs represen t a natural extension of mo dels with single v ariabilit y lev el, and they ha v e reen tly b een sub jet to a great deal of atten tion in statistial literature. In the eld of forestry , Hall and Clutter ( 1) analyze longitudinal measures of yield and gro wth that are mea- sured on ea h tree within a plot. In the eld of agronom y , Rek a y a et al. ( 2) onsider milk yield data where ea h o w is observ ed longitudinally during its rst three latations. Jarézi et al. (3) p erform geneti analyses of gro wth measuremen ts in b eef attle a kno wledging the fat that sev eral o ws ome from the same sire. Another example is p opulation PK, where onen tration measuremen ts ma y b e tak en with sev eral patien ts o v er sev eral distint time in terv als, that are often named p erio ds or o asions. That grouping pattern is used for instane in PK ross-o v er trials. In NLMEMs with only one lev el of v ariabilit y , often orresp onding to b et w een-sub jet v ari- abilit y , the analysis results in the estimation of the xed eets parameters and of the b et w een- sub jet v ariabilit y of the parameters, also alled in ter-sub jet v ariabilit y . When there is more than one lev el of grouping, the higher lev els of v ariabilit y an b e estimated. In the sp ei ase where the seond lev el of grouping orresp onds to m ultiple p erio ds of measuremen t, this v ariabil- it y is alled within-sub jet v ariabilit y (or in tra-sub jet v ariabilit y , or in ter-o asion v ariabilit y), and orresp onds to the v ariation of the individual parameters aross the dieren t study p erio ds or units. In the on text of pharmaokinetis, Karlsson et al. ( 4 ) demonstrate the imp ortane of mo deling this seond lev el of v ariabilit y in t w o-lev els NLMEMs. They sho w that negleting it resulted in biased estimates for the xed eets. 2 The parameter estimation of NLMEMs is not trivial b eause the lik eliho o d of NLMEMs annot b e expressed in a losed form due to the non-linearit y of the regression funtion in the random eets. Therefore, sev eral estimation metho ds ha v e b een prop osed. The First Order Conditional Estimates (F OCE) algorithm p erforms a rst order linearization of the regression funtion with resp et to the random eets (5; 6). The implemen tation of the F OCE algorithm in the NONMEM soft w are and in the nlme funtion of Splus and R enables the estimation of b oth b et w een- and within-sub jet v ariabilities. F rom our pratie, the main dra wba k of this metho d is ho w ev er that it do es not alw a ys on v erge when one estimates sim ultaneously the b et w een- and the within-sub jet v ariabilities on sev eral parameters. F urthermore, this linearization-based metho d annot b e onsidered as fully established in theory . F or instane, V onesh ( 7) and Ge et al. (8) giv e examples of sp ei designs resulting in inonsisten t estimates, su h as when the n um b er of observ ations p er sub jet do es not inrease faster than the n um b er of sub jets or when the v ariabilit y of random eets is to o large. Sev eral estimation metho ds ha v e b een prop osed as alternativ es to linearization algorithms. A ommon metho d to handle n umerial in tegration is the adaptativ e Gaussian quadrature (A GQ) metho d. An estimation algorithm of NLMEM parameters based on this lassial A GQ metho d has b een prop osed b y Pinheiro and Bates (9) and is implemen ted in the SAS pro edure NLMIXED (10). Ho w ev er, the A GQ metho d requires a suien tly large n um b er of quadrature p oin ts implying an often slo w on v ergene with v ery high omputational time. F urthermore, t w o-lev els NLMEM an b e implemen ted in the NLMIXED pro edure, but to our kno wledge, the on v ergene is diult to obtain in pratie (3 ). Impro v emen ts up on this metho d are th us needed. The seond alternativ e to linearization is the use of the Exp etation-Maximization (EM) algorithm (11 ) in order to estimate mo dels with missing or non-observ ed data su h as random eets. Beause of the nonlinearit y of the mo del, sto hasti v ersions of the EM algorithm ha v e b een prop osed. W ei et al. (12 ); W alk er (13 ) and W u ( 14 ) prop ose MCEM algorithms, with a Mon te-Carlo appro ximation of the E-step. Ho w ev er the MCEM algorithm ma y ha v e omputa- tional problems (i.e slo w or ev en non on v ergene). As an alternativ e to address omputational problems, a sto hasti appro ximation v ersion of EM (SAEM) has b een prop osed in ( 15 ; 16 ), whi h requires the sim ulation of only one realization of the missing data for ea h iteration, sub- stan tially reduing the omputation time. Kuhn and La vielle (16 ) prop ose to om bine the SAEM algorithm with a Mon te-Carlo Mark o v Chain (MCMC) pro edure adapted to the NLMEMs, and 3 pro v e that the resulting estimates are on v ergen t and onsisten t. T o date, none of the EM-based algorithms are diretly appliable to the ase of m ultilev el NLMEMs and ha v e to b e adapted. The ob jetiv e of this pap er is to extend the SAEM algo- rithm to MNLMEMs with t w o lev els of v ariabilit y: b oth E and M steps need to b e adapted to in tegrate higher lev els of random eets. W e also prop ose estimates of the lik eliho o d and of the Fisher information matrix. W e ev aluate this algorithm on a PK example, more preisely a t w o-p erio ds one-sequene ross-o v er trials sim ulated mimi king theoph yllin onen tration data (9). W e also apply the SAEM algorithm to the PK in teration of t w o HIV moleules (tenofo vir and atazana vir) from a PK substudy of the ANRS 107-Puzzle 2 trial. After desribing the mo del and notations (setion 2), setion 3 desrib es the SAEM algorithm. Setion 4 rep orts the sim ulation study and its results. In Setion 5, w e study the PK in teration of tenofo vir on atazana vir in HIV patien ts. Setion 6 onludes the artile with some disussion. 2 Mo dels Let us denote y ij k the observ ation in unit k ( k = 1 , . . . , K ) for sub jet i ( i = 1 , · · · , n ) and at time t ij k ( j = 1 , · · · , n ik ) . F or instane, the dieren t units an b e the dieren t p erio ds in the ase of PK trials, or the dieren t paren ts in the ase of geneti analyses. W e assume, as a kno wn fat, t w o nonlinear funtions f and g su h that the t w o-lev els non-linear mixed eets mo del linking observ ations to sampling times an b e written as: y ij k = f ( t ij k , φ ik ) + g ( t ij k , φ ik ) ε ij k , ε ij k ∼ N (0 , σ 2 ) , where φ ik is the p -v etor of the parameters of sub jet i for unit k and ε ij k is the measuremen t error. W e h yp othesize that the errors ε ij k giv en φ ik are m utually indep enden t. W e assume that the individual parameters φ ik are random v etors and that for ea h unit k , φ ik an b e brok en 4 in to: φ ik = µ + β k + b i + c ik , (1) b i ∼ N (0 , Ω) , c ik ∼ N (0 , Ψ) , where µ + β k is the mean v alue for unit k , b i is the random eet of size p of sub jet i , and c ik is the random eet of size p of sub jet i and unit k . T o ensure the iden tiabilit y of the parameters, w e assume that β 1 = 0 , ie µ is the mean of the rst unit and β k represen ts the dierene (or eet) of the k th unit in omparison to this rst unit. The random eets ( b i ) and ( c ik ) are assumed to b e m utually indep enden t. The total v ariane of the parameters is th us brok en in to a b et w een-sub jet v ariane Ω and a within-sub jet v ariane Ψ . Finally , the individual parameters pK -v etor φ i = ( φ i 1 , . . . , φ iK ) of sub jet i is distributed with a Gaussian distribution with mean v etor ( µ, µ + β 2 , . . . , µ + β K ) and a pK × pK o v ariane matrix Γ equal to Γ = Ω + Ψ Ω . . . Ω Ω Ω + Ψ . . . . . . . . . . . . . . . Ω Ω . . . Ω Ω + Ψ . Let θ = ( µ, β , Ω , Ψ , σ 2 ) , the v etor of all the parameters of the mo del where β denotes the v etor of unit eet β = ( β 1 , . . . , β K ) . The aim of this pap er is to prop ose an estimation of θ b y maximizing the lik eliho o d of the observ ations y = ( y ij k ) ij k . Let us denote ˜ b i := µ + b i . The lik eliho o d of y an b e written as: p ( y ; θ ) = Z p ( y , φ, ˜ b ; θ ) d ( φ, ˜ b ) where p ( y , φ, ˜ b ; θ ) is the lik eliho o d of the omplete data ( y , φ, ˜ b ) , with φ = ( φ ik ) i =1 ,...,n,k =1 , ...,K and ˜ b = ( ˜ b 1 , . . . , ˜ b n ) . Beause of the nonlinearit y of the regression funtion f with resp et to the random eets φ ik , the lik eliho o d has no losed form. Therefore, the maximization of the lik eliho o d in θ , θ ∈ Θ , is a omplex problem. W e prop ose to use a sto hasti v ersion of the EM algorithm, whi h is presen ted in detail in the next setion. 5 3 Estimation algorithm 3.1 The SAEM algorithm The EM algorithm in tro dued b y Dempster et al. ( 11 ) is a lassial approa h to estimate param- eters of mo dels with non-observ ed or inomplete data. In t w o-lev els NLMEMs, the non-observ ed data are the individual parameters ( φ, ˜ b ) and the omplete data of the mo del is ( y , φ, ˜ b ) . Let us denote L c ( y , φ, ˜ b ; θ ) = log p ( y , φ, ˜ b ; θ ) the log-lik eliho o d of the omplete data. The priniple of the iterativ e EM algorithm is to maximize the funtion Q ( θ | θ ′ ) = E ( L c ( y , φ, ˜ b ; θ ) | y ; θ ′ ) where the exp etation is the onditional exp etation under the p osterior distribution p ( φ, ˜ b | y ; θ ′ ) , the maximization of Q b eing often easier than the diret maximization of the observ ed data log- lik eliho o d. Ea h iteration of the EM algorithm is omputed through t w o steps: the Exp etation step (E-step) and the Maximization step (M-step). A t the ℓ th iteration of the algorithm, the E step is the ev aluation of Q ( θ | b θ ℓ ) , while the M step up dates b θ ℓ b y maximizing Q ( θ | b θ ℓ ) . Let us sho w that the funtion Q an b e redued in the ase of a MNLMEM. First, let us note that as p ( ˜ b | y , φ ; θ ) = p ( ˜ b | φ ; θ ) , b y appliation of the Ba y es theorem w e ha v e: p ( φ, ˜ b | y ; θ ) = p ( ˜ b | φ ; θ ) p ( φ | y ; θ ) = n Y i =1 p ( ˜ b i | φ i ; θ ) p ( φ i | y i ; θ ) . (2) Seond, through the linearit y of the individual parameters mo del in equation ( 1), the p oste- rior distribution p ( ˜ b i | φ i ; θ ) of the i th sub jet is expliit: p ( ˜ b i | φ i ; θ ) is a Gaussian distribution N ( m ( φ i , θ ) , V ( θ )) of mean and v ariane equal to: m ( φ i , θ ) = V ( θ ) Ψ − 1 K X k =1 ( φ ik − β k ) + Ω − 1 µ ! , (3) V ( θ ) = (Ω − 1 + K Ψ − 1 ) − 1 . Due to the fatorization giv en in equation (2), funtion Q an b e rewritten as: Q ( θ | θ ′ ) = Z Z L c ( y , φ, ˜ b ; θ ) p ( ˜ b | φ ; θ ′ ) d ˜ b p ( φ | y ; θ ′ ) dφ. Beause of the expliit p osterior distribution of random eets ˜ b giv en in equation (3), the omputation of this onditional exp etation an b e split in to t w o parts : the omputation of 6 the in tegral with resp et to the p osterior distribution of ˜ b whi h has an analytial form, and the omputation of the in tegral with resp et to the p osterior distribution of φ whi h has no analytial form. Therefore the EM algorithm is split in to an expliit diret EM algorithm for the omputation of the rst in tegral and the use of a sto hasti v ersion of the EM algorithm for the omputation of the seond in tegral. Let us detail the expliit omputation of the rst in tegral, denoted b y R ( y , φ, θ , θ ′ ) R ( y , φ, θ , θ ′ ) = Z L c ( y , φ, ˜ b ; θ ) p ( ˜ b | φ ; θ ′ ) d ˜ b. This in tegral has an analytial form. Indeed, the omplete log lik eliho o d L c ( y , φ, ˜ b ; θ ) is equal to L c ( y , φ, ˜ b ; θ ) = − 1 2 K X k =1 n X i =1 n ik X j =1 log(2 π σ 2 g 2 ( t ij k , φ ik )) − 1 2 K X k =1 n X i =1 n ik X j =1 ( y ij k − f ( t ij k , φ ik )) 2 σ 2 g 2 ( t ij k , φ ik ) − nK 2 log(2 π det Ψ) − 1 2 K X k =1 n X i =1 ( φ ik − ˜ b i − β k ) t Ψ − 1 ( φ ik − ˜ b i − β k ) − n 2 log(2 π det Ω) − 1 2 n X i =1 ( ˜ b i − µ ) t Ω − 1 ( ˜ b i − µ ) . As the p osterior distribution p ( ˜ b | φ ; θ ) is kno wn (equation 3 ), R ( y , φ, θ , θ ′ ) is equal to R ( y , φ, θ , θ ′ ) = − 1 2 X i,j,k log(2 π σ 2 g 2 ( t ij k , φ ik )) 1 2 X i,j,k ( y ij k − f ( t ij k , φ ik )) 2 σ 2 g 2 ( t ij k , φ ik ) (4) − nK 2 log(2 π det Ψ) − nK 2 Ψ − 1 / 2 V ( θ ′ )Ψ − 1 / 2 − 1 2 X i,k ( φ ik − m ( φ i , θ ′ ) − β k ) t Ψ − 1 ( φ ik − m ( φ i , θ ′ ) − β k ) − n 2 log(2 π det Ω) − n 2 Ω − 1 / 2 V ( θ ′ )Ω − 1 / 2 − 1 2 n X i =1 ( m ( φ i , θ ′ ) − µ ) t Ω − 1 ( m ( φ i , θ ′ ) − µ ) . Therefore Q is redued to the omputation of the seond in tegral under the p osterior distri- bution p ( φ | y ; θ ) as follo ws: Q ( θ | θ ′ ) = Z R ( y , φ, θ , θ ′ ) p ( φ | y ; θ ′ ) dφ. (5) Giv en the non-linearit y of funtion f with resp et to φ , the p osterior distribution p ( φ | y ; θ ′ ) has no losed form and the funtion Q dened b y (5) is in tratable. Th us w e prop ose to use the 7 sto hasti v ersion SAEM of the EM algorithm prop osed b y Dely on et al. (15 ) whi h ev aluates the in tegral Q b y a sto hasti appro ximation pro edure. Let us detail this SAEM algorithm in the ase of t w o-lev els NLMEMs. Let us note that the quan tit y R ( y , φ, θ , θ ′ ) b elongs to the regular urv ed exp onen tial family , i.e, it an b e written as R ( y , φ, θ , θ ′ ) = − Λ( θ ) + h S ( y , φ, θ ′ ) , Φ( θ ) i , (6) where h ., . i is the salar pro dut, Λ and Φ are t w o funtions t wie on tin uously dieren tiable on Θ and S ( y , φ, θ ′ ) is kno wn as the minimal suien t statistis of the omplete mo del. Those statistis are detailed later. In this ase, the Q funtion is redued to Q ( θ | θ ′ ) = − Λ( θ ) + h Z S ( y , φ, θ ′ ) p ( φ | y ; θ ′ ) dφ , Φ( θ ) i , In this ase, at the ℓ th iteration, the SAEM algorithm pro eeds as follo ws: • Sim ulation step: sim ulation of the missing data ( φ ( ℓ ) i ) i under the onditional distribution p ( φ | y ; b θ ℓ ) • Sto hasti appro ximation step: omputation of a sto hasti appro ximation s ℓ +1 of E h S ( y , φ, b θ ℓ ) | y ; b θ ℓ i = R S ( y , φ, b θ ℓ ) p ( φ | y ; b θ ℓ ) dφ , using ( γ ℓ ) ℓ ≥ 0 , a sequene of p ositiv e n um b ers dereasing to 0: s ℓ +1 = s ℓ + γ ℓ ( S ( y , φ ( ℓ ) , b θ ℓ ) − s ℓ ) . • Maximization step: up date of the estimate b θ ℓ +1 : b θ ℓ +1 = ar g max θ ∈ Θ ( − Λ( θ )) + h s ℓ +1 , Φ( θ )) i ) . The suien t statistis of the omplete mo del (4) ev aluated during the SA step of the SAEM 8 algorithm are as follo ws: s 1 ,i,ℓ +1 = s 1 ,i,ℓ + γ ℓ K X k =1 φ ( ℓ ) ik − s 1 ,i,ℓ ! , i = 1 , . . . , N , s 2 ,k,ℓ +1 = s 2 ,k,ℓ + γ ℓ n X i =1 φ ( ℓ ) ik − s 2 ,k,ℓ ! , k = 1 , . . . , K, s 3 ,ℓ +1 = s 3 ,ℓ + γ ℓ n X i =1 m ( φ ( ℓ ) i , b θ ℓ ) t m ( φ ( ℓ ) i , b θ ℓ ) − s 3 ,ℓ ! , s 4 ,ℓ +1 = s 4 ,ℓ + γ ℓ K X k =1 n X i =1 φ ( ℓ ) ik − m ( φ ( ℓ ) i , b θ ℓ ) t φ ( ℓ ) ik − m ( φ ( ℓ ) i , b θ ℓ ) − s 4 ,ℓ ! , s 5 ,ℓ +1 = s 5 ,ℓ + γ ℓ X i,j,k y ij k − f ( t ij k , φ ( ℓ ) ik ) g ( t ij k , φ ( ℓ ) ik ) ! 2 − s 5 ,ℓ , The expression of the M step is obtained b y deriv ation of equation (4). The parameter estimates are as follo ws: b µ ℓ +1 = V ( b θ ℓ ) b Ψ − 1 ℓ 1 n n X i =1 s 1 ,i,ℓ +1 − K X k =1 b β k,ℓ ! + V ( b θ ℓ ) b Ω − 1 ℓ b µ ℓ , b β k,ℓ +1 = s 2 ,k,ℓ +1 n − b µ ℓ +1 , for k = 2 , . . . , K , b Ω ℓ +1 = V ( b θ ℓ ) + s 3 ,ℓ +1 n − ( b µ ℓ +1 ) t b µ ℓ +1 , b Ψ ℓ +1 = V ( b θ ℓ ) + s 4 ,ℓ +1 nK − 1 K K X k =1 ( b β k,ℓ +1 ) t b β k,ℓ +1 , c σ 2 ℓ +1 = s 5 ,ℓ +1 P n i =1 P K k =1 n ik . Comparing with the lassi SAEM algorithm for single-lev el NLMEMs, the extension of SAEM to the t w o-lev els NLMEMs is nally split in to an expliit EM algorithm and a sto hasti EM part. F urthermore, it requires the omputation of t w o in termediate quan tities (the ondi- tional exp etations m ( φ i , θ ) and v ariane V ( θ ) of the b et w een-sub jet random eets parameters b i ) as w ell as t w o additional suien t statistis ( S 3 and S 4 ), funtions of m ( φ i , θ ) . The M-step diers from the one of the lassi SAEM for single-lev el NLMEMs, esp eially for the estimation of the v ariane matrix Ω and Ψ whi h uses the additional quan tit y V ( θ ) . As pro v ed b y (15 ; 17 ), the on v ergene of the SAEM algorithm is ensured under the follo wing assumption: 9 Assumption ( A1 ): 1. F untions Λ and Φ are t wie on tin uously dieren tiable on Θ . 2. The log-lik eliho o d log p ( y ; θ ) is d times dieren tiable on Θ , where d is the dimension of S ( y , φ, θ ′ ) . 3. F untion ¯ s dened as ¯ s ( θ, θ ′ ) = Z S ( y , φ, θ ′ ) p ( φ | y ; θ ) dφ is on tin uously dieren tiable on Θ with resp et to its rst v ariable. 4. F or all ℓ in N , γ ℓ ∈ [0 , 1] , P ∞ ℓ =1 γ ℓ = ∞ and P ∞ ℓ =1 γ 2 ℓ < ∞ . F or a on v enien t step sizes sequene γ ℓ , the assumption ( A1 ) is trivially he k ed in our mo del. A hoie of step sizes sequene γ ℓ is presen ted in Setion 4 . Ho w ev er, the sim ulation step of the SAEM algorithm, whi h p erforms the sim ulation of the non-observ ed v etor φ under the p osterior distribution p ( φ | y ; θ ) annot b e diretly p erformed b eause the p osterior distribution is only kno wn up to a onstan t. In this ase, Kuhn and La vielle (16 ) prop ose to om bine the SAEM algorithm with a Mark o v Chain Mon te Carlo (MCMC) pro edure for the sim ulation step. This v ersion of the SAEM-MCMC algorithm an b e used for the estimation of MNLMEMs. The MCMC pro edure used in this ase is detailed in setion 3.2. As pro v ed b y (17), the on v ergene of the SAEM-MCMC algorithm is ensured under assumption ( A1 ) and the follo wing additional assumption: Assumption ( A2 ): F or an y θ in Θ , w e assume that the onditional distribution p ( . | y ; θ ) is the unique limiting distribution of a transistion probabilit y Π θ , that has the follo wing prop erties: 1. F or an y ompat subset V of Θ , there exists a real onstan t L su h that for an y ( θ, θ ′ ) in V 2 sup { φ,φ ′ }∈E | Π θ ( φ ′ | φ ) − Π θ ′ ( φ ′ | φ ) | ≤ L k θ − θ ′ k . 2. The transition probabilit y Π θ supplies an uniformly ergo di hain whose in v arian t proba- bilit y is the onditional distribution p ( φ | y ; θ ) , i.e. ∃ K θ ∈ R + , ∃ ρ θ ∈ ]0 , 1 [ | ∀ ℓ ∈ N k Π ℓ θ ( ·| φ ) − p ( · , ·| y ; θ ) k T V ≤ C θ ρ ℓ θ , 10 where k · k T V is the total v ariation norm. F urthermore, C = sup θ ∈ Θ C θ < ∞ and ρ = sup θ ∈ Θ ρ θ < 1 . 3. F untion S is b ound on E . A t iteration ℓ , the S-step of the SAEM-MCMC algorithm onsists th us in sim ulating φ ( ℓ ) with the transistion probabilit y Π ˆ θ ℓ φ ( ℓ − 1) | dφ ( ℓ ) . The assumption ( A2 .2) is the most deliate to he k, and it dep ends on the hoie of the MCMC algorithm. This is detailed in the next subsetion after presen ting the MCMC pro edure. In pratie, the SAEM algorithm b eing a sto hasti algorithm, there exists no deterministi on v ergene riterion whi h ould b e used to stop the iterations of the algorithm as so on as the on v ergene is rea hed. Therefore, w e reommend to implemen t the SAEM algorithm with a suien tly large n um b er of iterations and to graphially he k the on v ergene b y plotting the v alues of the SAEM estimates obtained along iterations v ersus the iterations. Su h a gure is desrib ed in Setion 4. 3.2 MCMC algorithm for the sim ulation step Let us detail the sim ulation step of the SAEM-MCMC algorithm, whi h p erforms the sim ula- tion of the missing data φ through a Mark o v hain whi h has p ( φ | y ; θ ) as unique stationary distribution. F or sub jet i , b y Ba y es form ula, this onditional distribution is prop ortional to p ( φ i | y i ; θ ) ∝ K Y k =1 p ( y ik | φ ik ; θ ) p ( φ i ; θ ) . W e prop ose to use a Metrop olis-Hastings (M-H) algorithm to sim ulate this Mark o v hain. Let us reall the priniple of this algorithm. A t iteration r of the M-H algorithm, giv en the urren t v alue φ ( r ) i of the Mark o v Chain, the M-H algorithm pro eeds as follo ws: 1. Sim ulate φ c i with a prop osal distribution q ( · , φ ( r ) i ) 2. Compute the aeptane probabilit y α ( φ c i , φ ( r ) i ) = Q K k =1 p ( y ik | φ c ik ; θ ) p ( φ c i ; θ ) Q K k =1 p ( y ik | φ ( r ) ik ; θ ) p ( φ ( r ) i ; θ ) q ( φ c i , φ ( r ) i ) q ( φ ( r ) i , φ c i ) 11 3. Sim ulate u with a uniform distribution U [0 , 1] 4. Up date the Mark o v hain with φ ( r +1) i = φ c i if u ≤ α ( φ c i , φ ( r ) i ) φ ( r ) i else The on v ergene of the M-H algorithm strongly dep ends on the hoie of the prop osal distribu- tion q . The on v ergene is ensured for some prop osal distributions su h as indep enden t ( q ( · , φ ( r ) i ) indep enden t of φ ( r ) i ) or symmetrial ( q ( · , φ ( r ) i ) = q ( φ ( r ) i , · ) ) prop osals (18 ). These prop osals are detailed b elo w. Giv en the dimension of φ , w e also onsider a Metrop olis-Hastings-Within-Gibbs algorithm, om bining b oth Gibbs algorithm and M-H pro edure. The adv an tage of the Gibbs algorithm is to redue the m ulti-dimensional sim ulation problem to the suessiv e sim ulations of one-dimension v etors. Finally , at iteration ℓ of the SAEM algorithm, giv en the urren t estimate b θ ℓ , w e om bine the three follo wing prop osal transitions: 1. the prior distribution of φ i , that is the Gaussian distribution N ( b µ ℓ + b β ℓ , b Γ ℓ ) , orresp onding to an indep enden t M-H algorithm, 2. the m ultidimensional random w alks N ( φ ( ℓ − 1) i , ρ b Γ ℓ ) (symmetri prop osal), where ρ is a saling v alue hosen to ensure a satisfatory aeptation rate, namely around 30% as prop osed in (19 ), 3. a suession of K p unidimensional Gaussian random w alks (symmetri prop osal), i.e ea h omp onen t of φ i is suessiv ely up dated, leading to a Metrop olis-Hastings-Within-Gibbs algorithm, where b Γ ℓ is equal to b Γ ℓ = b Ω ℓ + b Ψ ℓ b Ω ℓ . . . b Ω ℓ b Ω ℓ b Ω ℓ + b Ψ ℓ . . . . . . . . . . . . . . . b Ω ℓ b Ω ℓ . . . b Ω ℓ b Ω ℓ + b Ψ ℓ . Giv en the prop osal distributions, as previously detailed, and using the theoretial on v ergene results prop osed in (18), this h ybrid Gibbs algorithm on v erges and generates an uniformly 12 ergo di hain with p ( φ | y ; θ ) as the stationary distribution. Consequen tly , b y applying the on- v ergene theorem of Kuhn and La vielle ( 16) and under assumptions (A1) and (A2), w e pro v e that the estimate sequene ( b θ ℓ ) ℓ ≥ 0 pro dued b y the extended SAEM algorithm on v erges to w ards a (lo al) maxim um of the lik eliho o d p ( y ; θ ) . In pratie, the on v ergene of the MCMC algorithm is diult to v erify . As in Ba y esian inferene, the only on v ergene riteria existing for MCMC pro edure are graphial riteria. W e ha v e to he k if the estimate sequene explores a suieny large domain of the Mark o v hain. A on v ergene gure is presen ted and ommen ted in setion 4 . 3.3 Estimation of the Fisher information matrix and the lik eliho o d T o p erform statistial tests su h as W ald test or lik eliho o d ratio test, w e prop ose estimators of the Fisher information matrix and the lik eliho o d, resp etiv ely . As the Fisher information matrix has no losed form in MNLMEMs, w e prop ose to appro ximate it b y the Fisher information matrix of the m ulti-lev el linear mixed mo del dedued from the MNLMEM after linearization of the funtion f around the onditional exp etation of the individual parameters ( E ( φ i | y ; ˆ θ ) , 1 ≤ i ≤ n ) . The omputation of this linearized Fisher information matrix is diret and do es not need an y appro ximation. The estimation of the lik eliho o d of the MNLEM is based on an Imp ortane Sampling pro- edure, as prop osed b y Samson et al. ( 20) for NLMEMs. The Imp ortane Sampling pro edure has b een in tro dued to appro ximate the in tegral of the lik eliho o d with a smaller v ariane than with other Mon te Carlo metho ds. In this ase, an estimate of the on tribution p ( y i ; θ ) of the individual i to the lik eliho o d is \ p ( y i ; θ ) = 1 T T X t =1 p ( y i , φ ( t ) i ; θ ) q i ( φ ( t ) i ) where ( φ ( t ) i ) t =1 ,...,T are sim ulated using the individual instrumen tal distribution q i . As an in- dividual instrumen tal distribution q i , w e prop ose a Gaussian appro ximation of the individual onditional p osterior distribution p ( φ i | y i ; θ ) . 13 4 Sim ulation study: a PK example 4.1 Sim ulation settings The ob jetiv e of this sim ulation study is to illustrate the main statistial prop erties of the extended SAEM algorithm (bias, ro ot mean square errors, group omparison tests) and to ompare them to the F OCE algorithm. W e do not use the A GQ algorithm, sine pro edure NLMIXED do es not sueed, in pratie, in estimating su h omplex v ariane mo dels. W e use the PK data of orally administered theoph yllin to dene the p opulation mo del for the sim ulation study . These data are lassial ones in p opulation pharmaokinetis, often used for soft w are ev aluation (21 ). W e assume that onen trations an b e desrib ed b y a one-ompartmen t mo del with rst order absorption and rst order elimination: f ( t, φ ) = D K a V K a − C l e − C l V t − e − K a t where D is the dose, V is the v olume of distribution, K a is the absorption rate onstan t and C l is the learane of the drug elimination. These parameters are p ositiv e and distributed aording to a log-normal distribution. Th us, φ has the follo wing omp onen ts: φ = (log V , log K a , log AU C ) , with AU C = D /C l . W e assume iden tial sampling times for all sub jets: for all i in 1 , . . . , n , k = 1 , 2 , t ij k = t j for j = 1 , . . . , J . A dditiv e Gaussian random eets are assumed for ea h parameter with a diagonal o v ariane matrix Ω and a a diagonal o v ariane matrix Ψ . Let ω 2 = ( ω 2 V , ω 2 K a , ω 2 AU C ) and ψ 2 = ( ψ 2 V , ψ 2 K a , ψ 2 AU C ) denote the v etor of the v arianes of the random eets. A om bined error mo del is assumed b y setting g ( t, φ ) = 1 + f ( t, φ ) . W e set the dose for all sub jets to the v alue of 4 mg. F or all the parameters, the v alues are those prop osed b y P anhard and Men tré ( 22 ): log V = − 0 . 73 , log K a = 0 . 39 and log AU C = 4 . 6 1 , ω 2 V = 0 . 01 , ω 2 K a = 0 . 04 , ω 2 AU C = 0 . 04 , ψ 2 V = 0 . 002 5 , ψ 2 K a = 0 . 01 , ψ 2 AU C = 0 . 01 and σ 2 = 0 . 01 . W e generate n = 24 and n = 40 total n um b ers of sub jets with J = 1 0 blo o d samples p er sub jet, tak en at 15 min utes, 30 min utes, 1, 2, 3.5, 5, 7, 9, 12 and 24 hours after dosing. The individual data of one sim ulated trial are displa y ed in Figure 1. [Figure 1 ab out here.℄ 14 4.2 Ev aluation of estimates Our aim is to ev aluate and ompare the estimates pro dued b y the extended SAEM algorithm with those pro dued b y the nlme funtion of the R soft w are. W e t the sim ulation mo del and ompute the relativ e bias and relativ e ro ot mean square error (RMSEs) for ea h omp onen t of θ from 1000 repliations of the t w o trials desrib ed b elo w ( n = 24 and n = 40 total n um b er of sub jets). The sim ulation mo del inludes a treatmen t eet on all omp onen ts of θ . W e test the n ull h yp othesis { β log AU C = 0 } using the W ald test. W e also t the mo del where the treatmen t eet on log AU C is not estimated, and test the same n ull h yp othesis using the Lik eliho o d Ratio T est (LR T). The SAEM algorithm is implemen ted with 500 iterations. During the rst 200 iterations, a onstan t step size γ ℓ = 1 is hosen, in order to let the Mark o v hain explore the parameters do- main. Then during the last 300 iterations, the sto hasti appro ximation s heme is implemen ted with a step size equal to γ ℓ = 1 ℓ − 200 at iteration ℓ . This hoie of step size sequene v eries on v ergene assumption ( A1 .1). The ev olution of ea h SAEM parameter estimates is plotted against iterations (logarithmi sale) on Figure 2 . During the rst iterations of the SAEM algorithm, the estimate sequenes explore randomly some neigh b orho o ds of the initial v alues, through the Mark o v hain sim ulation. In partiular, these b eha viors are learly visible for the xed eet parameters ( µ and β ). After 200 iterations, the estimates on v erge then rapidly to a neigh b orho o d of the maxim um lik eliho o d, due to the sto hasti appro ximation s heme. In this example, the iteration n um b er has b een hosen su h that the on v ergene is learly attained b efore the last iteration. [Figure 2 ab out here.℄ The relativ e bias and RMSEs obtained on the 1000 datasets with n = 24 and n = 40 sub jets are presen ted in T able 1. [T able 1 ab out here.℄ The bias and the RMSEs of the xed eets ( µ ) are small with the SAEM algorithm and are almost half of those obtained with nlme (esp eially the RMSEs). The bias and RMSEs of the unit eet ( β ) are small and on the same order with b oth metho ds. F or the b et w een-sub jet 15 v ariabilit y parameters ( Ω ), the bias are redued with SAEM, while the RMSEs are of the same order with b oth metho ds. F or the within-sub jet v ariabilit y parameters ( Ψ ), the bias and the RMSEs are satisfatory , and on the same order with b oth metho ds. The bias and RMSE for σ 2 are small and satisfatory for b oth metho ds. The t yp e I error of the W ald test and of the LR T are ev aluated on the same 1000 datasets. F or n = 24 , the t yp e I errors are 6.0% and 6.5% for SAEM and nlme, resp etiv ely , for the W ald test, and 4.6% and 5.6% for SAEM and nlme, resp etiv ely , for the LR T. F or n = 40 , the t yp e I error are 5.6% and 5.4% for SAEM and nlme, resp etiv ely , for the W ald test, and 5.8% and 5.2% for SAEM and nlme, resp etiv ely , for the LR T. 5 Appliation to the p opulation pharmaokinetis of atazana vir with tenofo vir 5.1 Study p opulation: ANRS 107 - Puzzle 2 study The Puzzle 2 - ANRS 107 trial w as a randomized op en-lab el, m ultiple-dose study supp orted b y the F ren h Agene Nationale de Re her he sur le Sida (ANRS) with HIV-infeted patien ts in treatmen t failure with their an tiretro viral therap y . P atien ts w ere randomized to reeiv e for the rst t w o w eeks either their un hanged treatmen t with PIs and n uleoside rev erse transrip- tase inhibitors (NR TIs) (group 1) or un hanged treatmen t with NR TIs in om bination with atazana vir (300 mg QD) plus ritona vir (100 mg QD) as a substitute for the failing PI therap y (group 2). F rom w eek 3 (da y 15) to w eek 26, patien ts from either group swit hed to atazana vir (300 mg QD) plus ritona vir (100 mg QD) plus tenofo vir DF at 300 mg QD and NR TIs seleted aording to the baseline rev erse transriptase genot yp e of the HIV isolated in ea h patien t. In this pap er, w e analyze onen tration data obtained from 10 patien ts from group 2 who w ere inluded and measured at 1, 2, 3, 5, 8, and 24 h after administering drug during ea h treatmen t p erio d. Those exat dosing in terv als w ere reorded. The ob jetiv e of the substudy w as to measure the pharmaokineti parameters of atazana vir (administered with ritona vir) either b efore (da y 14 [w eek 2℄) or after (da y 42 [w eek 6℄) initiation of tenofo vir DF in HIV- infeted patien ts in order to detet pharmaokineti in terations of tenofo vir on atazana vir. Data of this substudy w ere analyzed using a nonlinear mixed eet mo del b y P anhard et al ( 23 ) 16 and a signian t eet of the o-administration of tenofo vir on the pharmaokineti parameters of atazana vir w as found using the nlme funtion of the Splus soft w are. The aim of the presen t analysis is to ev aluate the eet of tenofo vir on the PK parameters of atazana vir using the SAEM algorithm and the W ald test desrib ed in setion 4. W e use the one-ompartmen t mo del with zero-order absorption prop osed b y P anhard et al. (23 ) to desrib e atazana vir onen trations: f ( t, φ ) = F D T a C l (1 − e − C l V t ) 1 t
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment