doi:10.2478/v10014-010-0017-x COBISS: 1.01 Agris category code: U10 SIMPLE REPARAMETERIZATION TO IMPROVE CONVERGENCE IN LINEAR MIXED MODELS Gregor GORJANC Tina FLISAR 2, Jose Carlos MARTÍNEZ-ÁVILA 3, Luis Alberto GARCÍA-CORTÉS 3 Received October 08, 2010; accepted December 01, 2010. Delo je prispelo 08. oktobra 2010, sprejeto 01. decembra 2010. Simple reparameterization to improve convergence in linear mixed models Slow convergence and mixing are one of the main problems of Markov chain Monte Carlo (McMC) algorithms applied to mixed models in animal breeding. Poor convergence is to a large extent caused by high posterior correlation between variance components and solutions for the levels of associated effects. A simple reparameterization of the conventional model for variance component estimation is presented which improves McMC sampling and provides the same posterior distributions as the conventional model. Reparameterization is based on the rescaling of hierarchical (random) effects in a model, which alleviates posterior correlation. The developed model is compared against the conventional model using several simulated data sets. Results show that presented reparameterization has better behaviour of associated sampling methods and is several times more efficient for the low values of heritability. Key words: statistics / mixed model / Bayesian analysis / McMC / reparameterization / convergence 1 INTRODUCTION Mixed models are abundantly used in the field of animal breeding and genetics with the aim to infer genetic values of animals given some phenotypic and pedigree information (Henderson, 1984). In it simplest form the mixed model can be written as: y = Xb + Za + e, (1) where y is a vector of phenotypes, b is a vector of effects Enostavna reparametrizacija za izboljšanje konvergence linearnih mešanih modelov Počasna konvergenca je eden največjih problemov uporabe metode Monte Carlo z Markovimi verigami (McMC) za mešane modele na področju genetike in selekcije domačih živali. Slaba konvergenca je v veliki meri posledica visoke posteriorne korelacije med komponentami variance in rešitvami za ravni pripadajočih vplivov. Predstavljamo enostavno reparametriza-cijo običajnega modela, ki izboljša lastnosti metode McMC in daje enake posteriorne porazdelitve parametrov modela kot standardni pristop. Reparametrizacija temelji na standardizaciji hierarhičnih (naključnih) vplivov v modelu, kar posledično spremeni posteriorne korelacije med parametri. Oba pristopa smo primerjali na večjem setu simuliranih podatkov. Rezultati kažejo, da reparametrizacija vodi do bolj učinkovitih metod McMC vzorčenja in je nekajkrat bolj učinkovita za analizo lastnosti z nizko heritabiliteto. Ključne besede: statistika / mešani model / bayesovska analiza / McMC / reparametrizacija / konvergenca like sex, breed, age, etc., a is a vector of individual additive genetic effects and e residual, p(e | a2)~ N(o,Ia2), while X and Z are design matrices linking effects to phenotypic records. Pedigree information is included in the model hierarchically with prior distribution of individual additive genetic values, p(a |A,a2)~ N(o, Aa2). Henderson (1972) developed the so called mixed model equations (2) to efficiently obtain joint solutions for b and a, where G=A a2 and R = Ia2 : 1 Univ. ofLjubljana, Biotechnical Fac., Dept. of Animal Science, Groblje 3, SI-1230 Domžale, Slovenia,, Ph.D., e-mail: gregor.gorjanc@bf.uni-lj.si 2 The same address as 1 3 Departamento de Mejora Genética, Instituto Nacional de Investigación Agraria, Carretera de La Coruña, km 7, 28040 Madrid, Spain, Ph.D. G. GORJANC et al. X T R-1X X T R-1Z b I I ATR-1y ZT R-1X ZT R-1Z + G-1 a I I ZT K'y Use of mixed model equations assumes known variance components ca and ce . Standard procedure is to estimate these variance components using restricted maximum likelihood method (REML; Patterson and Thompson, 1971) and to use these estimates in mixed model equations (2) ignoring the error of estimation in variance components. Another approach to statistical inference, Bayesian approach, treats inference of all model parameters jointly. Although conceptually very appealing, Bayesian approach leads to formulas that are computationally intractable. This can be avoided by sampling methods such as Markov chain Monte Carlo (McMC; e.g., Gelman, et al., 2004). Wang et al. (1993) showed how McMC methods can be used with linear mixed models in animal breeding applications. In the case of linear mixed models all McMC computations follow from the posterior distribution (3): p(b, a, c My |R| 'exp(-i(y-Xb - Za f R '(y-Xb-Za). X (3) G 2 exp(- i a f G-1 a) where prior distributions for and both variance components ca2 and2 ae were assumed uniform (e.g., Gelman et al., 2004). Given that cra and a are a priori correlated due to the prior definition of a, the a posteriori correlation between them is expected to be high. This leads to high autocorrelation between consecutive samples, making McMC method inefficient. Autocorrelations can be really problematic with low or near zero values for some variance components (e.g. additive genetic variance). This is caused by the shrinkage of a towards zero and in a next round of sampling variance component will again be close to zero, which can make the sampler stuck for quite some time at the values near zero (Gelman et al., 2004). Chib and Carlin (1999) proposed block sampling of some parameters in (2) to improve convergence. Autocorrelation has also been alleviated by the use of centered models (Gelfand et al., 1995), parameter expanded models (Liu and Wu, 1999; Gelman et al., 2003; Gelman, 2004) and data augmentation based models (Meng and van Dyk, 1997; van Dyk and Meng, 2001). These methods have been applied both to accelerate the Expectation-Maximization (EM) algorithm and to alleviate the autocorrelation of McMC algorithms. In this work a reparameterization will be employed where additive genetic values will be a priori uncorrelated with c a . This approach will be compared against the conventional model of Wang et al. (1993). 2 METHOD Let us consider a simple animal model y = Xb + Za + e with the following distributional assumptions: p(y | b, a, o2)~ N(Xb + Za, Io'„ ) p(a\A,o2 ) ~ N (o, Ao2 ) p(e|o2 )~ N(o, Io2e ) (4) For this particular case and assuming uniform priors for b and both variance components, p(b) x const., pCa )x const., and p(cA )J J Zt c ■ u N c , t (8) IX T .-1 u A " u SIMPLE REPARAMETERIZATION TO IMPROVE CONVERGENCE IN LINEAR MIXED MODELS ?(«• 1 u-i,b'CT2'CT2'y)~ N1 I, (9) where both S{ and c,j are closely related with the conventional mixed model (2) but modified as: C = ( XT X X1 Zaa ZT Xaa ZT + A~We f \Ty S = X y (10) y"ay The full conditional distribution of ae can be sampled from scaled inverted chi-square distribution with n - 2 degrees of freedom as in the conventional model: 4r 2 |b, u,a\, y)~ (y - Xb - Zuoa ) (y - Xb - Zuaa )zn-2 . (11) 2 After some algebra the full conditional of aa is P A 1 b, u, A, y exp _ 2 u'U (y-X u' Z' Zu b) ■a„ + a n (12) 2 u Z Zu from which a truncated normal distribution can be recognized when presented in terms of aa with mean uTz (y-xb) —Zi- variance u' ztzu and truncation point at 0: L ib, u, a 2 jy ) ~ TN{ — b), —TT, , 0 (13) , , tj 1 _ zr Zu —T Z' Z— 2 When the full conditional distribution of aa does not involve the neighbourhood ofzero, it is a scaled non-central xX distribution with 1 degree of freedom, with a scale parameter u T Z u and noncentrality parameter 2 _ uTZT (y-Xb)(y-Xb)T Zu 2 _ _ T _T__ 2uT ZT Zúa, A a 1 b y ) ~ U zjzux (u)- (14) For cases where the posterior distribution of a is close to zero, the Metropolis-Hastings algorithm with positive proposal can be implemented, where the natural logarithm of the conditional density derived from (12) is: -2z