Estimation of non-linear growth models by linearization : a simulation study using a Gompertz function

Gene . Sel. E ol. 38 (2006) 343–358 343 c INRA, EDP Sciences, 2006 DOI: 10.1051/gse:2006008 O iginal a icle Es ima ion o non-linea g ow h models by linea iza ion: a simula ion s udy using a Gompe z unc ion Kaa ina V∗,IsmoS´ , Ma ja-Liisa S´ -A, Esa A. M¨  MTT Ag i ood Resea ch Finland, Bio echnology and Food Resea ch, Biome ical Gene ics, FIN-31600 Jokioinen, Finland (Recei ed 6 July 2005; accep ed 27 Janua y 2006) Abs ac – A me hod based on Taylo se ies expansion o es ima ion o loca ion pa ame e s and a iance componen s o non-linea mixed effec s models was conside ed. An a ac i e p ope y o he me hod is he oppo uni y o an easily implemen ed algo i hm. Es ima ion o non-linea mixed effec s models can be done by common me hods o linea mixed effec s mod- els, and hus exis ing p og ams can be used a e small modi ica ions. The applicabili y o his algo i hm in animal b eeding was s udied wi h simula ion using a Gompe z unc ion g ow h model in pigs. Two g ow h da a se s we e analyzed: a ull se con aining obse a ions om he en i e g owing pe iod, and a unca ed ime ajec o y se con aining animals slaugh e ed p e- ma u ely, which is common in pig b eeding. The esul s om he 50 simula ion eplica es wi h ull da a se indica e ha he linea iza ion app oach was capable o es ima ing he o iginal pa- ame e s sa is ac o ily. Howe e , es ima ion o he pa ame e s ela ed o adul weigh becomes uns able in he case o a unca ed da a se . Gompe z unc ion /non-linea mixed effec s / a iance componen s /b eeding alues / likelihood app oxima ion 1. INTRODUCTION Non-linea unc ions a e pa icula ly sui ed o model g ow h da a, because p edic ions ou side he da a ange can be made mo e eliably han by linea models, and he en i e g ow h p ocess can be desc ibed by ew pa ame e s. Fo example, g ow h da a models commonly apply he Gompe z unc ion, whe e he es ima ed pa ame e s can ha e biological meaning. Non-linea models a e, howe e , mo e complica ed o sol e han linea models, and se e al algo i hms ∗Co esponding au ho : kaa ina. uo i@m . i A icle published by EDP Sciences and a ailable a h p://www.edpsciences.o g/gse o h p://dx.doi.o g/10.1051/gse:2006008 344 K. Vuo i e al. ha e been p oposed o es ima e he pa ame e s and a iance componen s o non-linea mixed effec s models [4]. In animal p oduc ion esea ch, he Bayesian amewo k has ecei ed much a en ion in g ow h cu e analysis [2, 11]. This popula i y is due o he use o Ma ko chain Mon e Ca lo me hods ha allow he solu ion o nume ically complica ed pos e io densi y in eg a ion and calcula ion o con idence in e - al es ima es. The cos o his p ocedu e is, howe e , in ensi e calcula ions and he need o assu e a sampling equilib ium [3]. Ano he possibili y is o ap- p oxima e he likelihood unc ion using linea iza ion [4, 10, 17, 18] o nume i- cal in eg a ion [10]. Also, bo h o hese al e na i es a e compu a ionally diffi- cul , ye , linea iza ion may enable he in e ence o linea mixed effec s models. All linea iza ion me hods in he li e a u e a e qui e simila . Ea ly me hods use i s -o de Taylo se ies expansion o non-linea unc ions a ound expec a ion o he andom effec s, and a e sol ed by ei he maximum likelihood (ML) o gene alized leas squa es (GLS) es ima ion [4]. Linds om and Ba es [9] sug- ges ed a mo e accu a e me hod o making he expansion a ound cu en es- ima es o he andom effec s. Subsequen esea ch has ocused mo e on he second-o de Taylo se ies expansion o in eg als in oked by he Laplacian app oxima ion [10, 17, 18]. Because o gene ali y and amilia o mula ion, he mos in e es ing choice o app oxima ion is based on he second-o de Taylo se ies expansion wi h espec o andom effec s ha we e p esen ed by Wol inge and Lin [18]. They ga e wo al e na i e app oaches o selec poin s o expansion: a ze o- expansion me hod using expec ed alues, and an EBLUP-expansion me hod using he empi ical bes linea unbiased p edic o s o he andom effec s. Bo h app oxima ions lead o algo i hms ha i e a i ely i mixed linea models o he sui ably ans o med da a using ei he ML o es ic ed maximum likeli- hood (REML). The e o e, hey allow he use o commonly applied me hods o linea mixed effec s models, and he use o exis ing p og ams a e small modi ica ions. A simila algo i hm was p oposed by B eslow and Clay on [3] in he con ex o gene alized linea mixed models. Because condi ions o unc- ionali y o he app oxima ion me hods a e difficul o iden i y, Wol inge and Lin [18] ecommended simula ion s udies o assessing he pe o mance o he me hods in di e se kinds o non-linea models and da a se s. The aim o his wo k was o desc ibe and examine he pe o mance o he EBLUP-expansion me hod o he Gompe z unc ion applied o he analysis o g ow h in he pig h ough simula ion. The EBLUP-expansion is ecommended especially o cases whe e he a iance componen s a e la ge, which may be he case o an adul weigh pa ame e o pigs. Also, Linds om and Ba es [9] Es ima ion o non-linea g ow h models 345 sugges ed ha he expansion a ound he expec ed ze o alue may lead o poo es ima es when subs an ial in e -indi idual a ia ion exis s. We chose o exam- ine he me hod h ough he analysis o wo da a se s. The i s analysis es ed gene al pe o mance o he EBLUP-expansion echnique, and he second anal- ysis es ed pe o mance o he me hod o incomple e da a. Incomple e da a a e common in pig p oduc ion, whe e he adul weigh is una ailable due o an ea lie slaugh e age. 2. MATERIALS AND METHODS 2.1. Simula ions The Gompe z unc ion has been shown o i pig g ow h da a, such as li e weigh and p o ein e en ion, well [14–16]. We assumed ha weigh s o an indi idual i ollowed he Gompe z model: yij =αexp(−βexp(−κ ij)) +eij,j=1,...,ni whe e niis he numbe o obse a ions o indi idual i,yij is he obse ed weigh a age ij (in days), α,βand κa e he pa ame e s o he Gompe z unc- ion, and eij is he andom esidual. The pa ame e s ha e biological meaning: αis he adul weigh , κis he a e o exponen ial decay o he ini ial g ow h a e, and βis he loga i hm o he a io o bi h weigh o adul weigh . Each o he pa ame e s α,βand κcan be desc ibed by a linea mixed effec s model. In his s udy, we will conside a si e model, al hough no a ion could be o an animal model. The ull model o obse a ion jo animal iis yij =(xαijbα+zs,αisα+zp,αipα) exp(−(xβijbβ+zs,βisβ+zp,βipβ) exp(−(xκijbκ+zs,κisκ+zp,κipκ) ij)) +eij,(1) whe e (bα,bβ,bκ)T=bis a d×1- ec o o ixed effec s, (sα,sβ,sκ)T=sis a l×1- ec o o andom addi i e gene ic si e effec s and (pα,pβ,pκ)T=pis a q×1- ec o o andom animal effec s o he han si e. Vec o s x,zsand zpa e om he design ma ices o ixed, andom si e and andom animal effec s X, Zsand Zp, espec i ely. I is assumed ha           s p e           ∼N                    0 0 0           ,          G00 0P0 00R                     . 346 K. Vuo i e al. He e, G=G0⊗A,whe eAis a ma ix o addi i e ela ionships be ween si es and G0isa3×3 gene ic co a iance ma ix o he Gompe z pa ame e s. Simila ly, P=P0⊗Iq,whe eP0is a 3 ×3 co a iance ma ix, i.e. he andom animal effec s pwe e iden ically and independen ly dis ibu ed o he animals. Fu he mo e, he esiduals we e assumed o be independen ly dis ibu ed and homoscedas ic, e∼N(0,Inσ2 e). The Gompe z unc ion coefficien s we e gene a ed o simula e pig g ow h. The model had one ixed effec wi h wo le els, and he andom effec s we e he gene ic si e effec and he animal effec o he han he si e. The i s ixed effec le el had alues 210, 5 and 0.017 o he pa ame e s α,βand κ, espec- i ely. The o he ixed effec le el had alues 220, 4.7 and 0.016 o α,βand κ, espec i ely. The andom effec s we e assumed o be no mally dis ibu ed wi h mean ze o and block diagonal co a iance ma ices. Va iance and co a iance componen s in he ma ices o he gene ic si e effec and o he animal effec a e shown in Tables I and II. The esidual a iance σ2 ewas one. These pa ame- e s app oxima ed he a iances calcula ed by he NLMIXED p ocedu e o he SAS p og am ha i ed he Gompe z model o g ow h pe o mance da a o Finnish pigs [12]. Simula ion o he andom si e effec equi ed aking in o accoun he pedi- g ee. The pedig ee had h ee gene a ions o animals wi h 10 un ela ed ounde g andsi es. Each o he 10 g andsi es was ma ed wi h 20 un ela ed dams ha p oduced one son each, i.e., 20 hal -sibs. The hal -sibs we e ma ed wi h un- ela ed dams o p oduce 24 p ogeny pe si e. Only he las gene a ion o ani- mals had eco ds. Thus, he da a included 4800 es ed animals. Two da a se s we e made: a comple e se , and a unca ed ime ajec o y se . The comple e da a con ained 30 equally-dis anced obse a ions pe animal be ween 50 and 253 days. The unca ed ime ajec o y da a con ained slaugh e weigh s up o 115 kg, which is simila o he common slaugh e weigh in pigs, and occu s a abou 120 days o age. Consequen ly, he numbe o obse a ions was e- duced om 30 o abou 11 pe animal, i.e., almos wo hi ds o he da a we e disca ded. 2.2. Me hod o es ima e he alues o he g ow h pa ame e s The non-linea mixed model conside ed was y= (X,b,Zs,s,Zp,p)+e, whe e yis an n×1- ec o o obse a ions, is he Gompe z unc ion, and e is an n×1- ec o o andom esiduals. Vec o s b,sand p, wi h ma ices X, Zs Es ima ion o non-linea g ow h models 347 Table I. Rela i e bias, ela i e s anda d de ia ion (Rel. SD) and ela i e mean squa ed e o (Rel. MSE) (as pe cen om he ue alue) o (co) a iance componen s o gene ic si e effec s om he 50 eplica es o ull and unca ed ime ajec o y da a. Subsc ip s α,βand κdeno e he h ee pa ame e s in he Gompe z unc ion. Full da a T unca ed da a Pa ame e T ue Rel. Bias (%) Rel. SD (%) Rel. MSE (%) Rel. Bias (%) Rel. SD (%) Rel. MSE (%) σ2 α10.0 −1.9 15.9 25.0 4.7 41.3 169.5 σ2 β0.01 0.3 15.2 2.27e–02 7.1 13.6 2.31e–02 σ2 κ3.0e–07 3.0 15.6 7.41e–07 −17.1 33.4 4.16e–06 σαβ −0.06 6.0 55.1 1.8 −19.9 84.9 4.5 σακ −3.0e–04 −9.7 64.5 1.25e–02 71.4 196.2 0.1 σβκ 1.0e-05 21.7 65.7 4.70e–04 −8.9 75.0 5.59e–04 Table II. Rela i e bias, ela i e s anda d de ia ion (Rel. SD) and ela i e mean squa ed e o (Rel. MSE) (as pe cen om he ue alue) o (co) a iance componen s o he animal effec om he 50 eplica es o ull and unca ed ime ajec o y da a. Subsc ip s α, βand κdeno e he h ee pa ame e s in he Gompe z unc ion. Full da a T unca ed da a Pa ame e T ue Rel. Bias (%) Rel. SD (%) Rel. MSE (%) Rel. Bias (%) Rel. SD (%) Rel. MSE (%) σ2 α90.0 2.98e–02 2.3 4.5 18.8 10.6 416.5 σ2 β0.09 −0.2 1.8 2.88e–03 2.2 2.9 1.16e–02 σ2 κ3.7e–06 0.2 2.1 8.91e–07 −0.9 3.0 3.47e–07 σαβ −0.54 −1.9 9.5 0.5 −36.7 24.5 10.4 σακ −3.75e–03 −1.5 6.5 1.64e–03 −13.3 25.5 3.05e−02 σβκ 9.0e–05 −0.9 9.9 8.72e–05 4.5 16.9 2.71e–04 348 K. Vuo i e al. and Zp, we e de ined as be o e. Fo he andom effec s, deno e uT=(sTpT), Z=ZsZpand D=G0⊗A0 0P 0⊗Iq. Now he dis ibu ion assump ions we e u e∼N0 0,D0 0R . Al hough Ris diagonal he e, any o m is allowed, so he gene al o m Rwill be used he eina e . Unknown elemen s o co a iance ma ices G0,P0and R a e deno ed by pa ame e ec o θ. The maximized likelihood unc ion was L(b,θ|y)=(2π)−n 2|R|−1 2(2π)−l+q 2|D|−1 2 exp −1 2(y− (X,b,Z,u))TR−1(y− (X,b,Z,u)) −1 2uTD−1udu.(2) Only in some cases is he closed o m o (2) ound, so he in eg al is o en sol ed nume ically. Howe e , nume ical me hods o he non-linea unc ions a e usually slow o con e ge and nume ically uns able. Ins ead, he in eg al may be app oxima ed by quad a ic Taylo -se ies expansion o he exponen . The second-o de expansion was made abou he EBLUP be o e in eg a ion o he likelihood unc ion (see Appendix). This ga e app oxima ion o he loga- i hm o he likelihood unc ion (2): l∗(b,θ|y)=−1 2nln (2π)−1 2ln(|R||I+Z∗TR−1Z∗D|) −1 2(y− (X,b,Z,˜ u))TR−1(y− (X,b,Z,˜ u)) −1 2˜ uTD−1˜ u,(3) whe e Z∗=∂ ∂uT|u=˜u and ˜ uis he empi ical BLUP-es ima e o andom effec s. Fo he Gompe z unc ion and wo andom effec s in he model, Z∗had elemen s ∂ ∂αi =exp(−βiexp(−κi j)) cαi ∂ ∂βi =αiexp(−βiexp(−κi j)) (−exp(−κi j)) cβi(4) ∂ ∂κi =αiexp(−βiexp(−κi j)) (−βiexp(−κi j)) (− j)cκi whe e cis zso zpdepending on he andom effec diffe en ia ed (see (1)). Es ima ion o non-linea g ow h models 349 Pinhei o and Ba es [10] used he app oxima ion (3) in es ima ion o pa ame- e s by he Laplacian app oxima ion. Howe e , no s aigh o wa d gene aliza- ion o he REML-es ima ion was p esen ed. Wol inge and Lin [18] de eloped o mula (3) u he . Deno e V=Z∗DZ∗T+R|u=˜u . Then, l∗(b,θ|y)=−1 2nln(2π)−1 2ln |V| −1 2(y− (X,b,Z,˜ u)+Z∗˜ u)TV−1(y− (X,b,Z,˜ u)+Z∗˜ u), whe e V−1=R−1−R−1Z∗D(I+Z∗TR−1Z∗D)−1Z∗TR−1and |V|=|R||I+ Z∗TR−1Z∗D|[5]. This led o a simila es ima ion unc ion o a iance compo- nen es ima ion p esen ed by Linds om and Ba es [9], al hough h ough di - e en de i a ion. 2.2.1. Es ima ion o he ixed and andom effec s Assume ha he a iance componen ec o θis known. Maximum likeli- hood es ima ion o he pa ame e s band uleads o sol ing equa ions: X∗TR−1(y− (X,˜ b,Z,˜ u)) =0 Z∗TR−1(y− (X,˜ b,Z,˜ u)) =D−1˜ u,(5) whe e X∗=∂ ∂bTb=˜ band ˜ bis he es ima e o ixed effec s b.Elemen s o X∗a e simila o Z∗, excep ha coefficien cin (4) is eplaced by xdue o he diffe en ia ed ixed effec . Howe e , in o de o a i e o hese simple equa ions, dependency o Von b h ough Z∗has o be igno ed. On he basis o a gumen s made by Ba es and Wa s [1], Wol inge and Lin [18] jus i ied his by appealing o in insic non-linea i y ins ead o non-linea i y o he pa ame- e s. Deno e Y=y− (X,˜ b,Z,˜ u)+X∗˜ b+Z∗˜ u. Equa ions (5) can now be w i en as X∗TR−1X∗X∗TR−1Z∗ Z∗TR−1X∗Z∗TR−1Z∗+G−1˜ b ˜ u=X∗TR−1Y Z∗TR−1Y.(6) This is simila o he mixed model equa ions (MME) o he linea models. Thus, al eady es ablished me hods o sol ing linea models can be used o analyse he pseudo-da a Yc ea ed om he o iginal da a ywi h ˜ band ˜ uequal o hei mos ecen es ima es. 350 K. Vuo i e al. 2.2.2. Es ima ion o he a iance componen s A e inding es ima es o he loca ion pa ame e s, p o ile likelihood can be used o es ima e he a iance componen s by se ing b=˜ b(θ). The loga i hmic likelihood unc ion o he pa ame e ec o θcan be w i en wi h he pseudo- da a as l∗ ML(θ)=−1 2nln(2π)−1 2ln |V|−1 2(Y−X∗˜ b)TV−1(Y−X∗˜ b).(7) Diffe en ia ing equa ion (7) wi h espec o θgi es −1 2 V-1 ∂V ∂θj+1 2(Y−X∗˜ b)TV-1 ∂V ∂θj V-1(Y−X∗˜ b).(8) Maximum likelihood es ima es o a iance componen s a e ound by equa ing (8) o ze o and sol ing o θ. Ins ead o he ML-es ima es, REML-es ima es a e commonly used in p ac- ise. These es ima es accoun o losses in deg ees o eedom caused by he es ima ion o ixed effec s b[5]. The loga i hmic likelihood unc ion is now l∗ REML(θ)=−1 2nln(2π)−1 2ln |V|−1 2ln |X∗TV−1X∗| −1 2(Y−X∗˜ b)TV−1(Y−X∗˜ b).(9) Diffe en ia ion wi h espec o θand equa ing o ze o gi es −1 2 P∂V ∂θj+1 2(Y−X∗˜ b)TV-1 ∂V ∂θj V-1(Y−X∗˜ b)=0,(10) whe e P=V−1−V−1X∗(X∗TV−1X∗)−1X∗TV−1. Solu ions in θa e REML- es ima es o a iance componen s. 2.2.3. EBLUP-algo i hm The app oxima e ML-solu ions o loca ion pa ame e s and a iance compo- nen s can be ob ained by i e a i ely sol ing he equa ions (6) and (8) un il con- e gence. Co espondingly, he REML-solu ions o he EBLUP-expansion a e ob ained by i e a i ely sol ing he equa ions (6) and (10). Hence, he algo i hm i s he linea mixed effec s model Y=X∗b+Z∗u+e o he pseudo-da a Y and he wo king ec o s X∗and Z∗,whe eu∼N(0,G(θ)) and e∼N(0,R(θ)). Es ima ion o non-linea g ow h models 351 2.3. Implemen a ion We chose o implemen he REML-based EBLUP-algo i hm, because p o- g ams o sol e he linea mixed effec s models a e a ailable and commonly used by animal b eede s. MiX99 [13] was used o sol e he mixed model equa- ions (6), and DMU [6], modi ied o andom eg ession by Ke unen e al. [7], was used o sol e he REML es ima es o co a iance componen s (10). The ca- pabili y o i ing ixed and andom eg ession models is c ucial o implemen- a ion, because he coefficien s in X∗and Z∗can ha e any alues. Implemen- a ion o he linea iza ion p ocedu e equi ed he Gompe z unc ion o mulas o be included in MiX99. Howe e , he e was no need o make changes o he a iance componen es ima ion p og am. S a ing alues o bo h he loca ion pa ame e effec s and he a iance com- ponen s had o be assigned be o e i s i e a ion. A na u al choice was o ini ialize andom effec s wi h he expec ed alue ze o. Howe e , ini ial al- ues o ixed effec s we e de i ed wi h he model unc ion o g ow h cu e and a ailable da a. When he Gompe z model is used, only he asymp o ic weigh pa ame e has a na u al ini ial alue, which is he maximum alue o he dependen a iable. In he o he cases, complex equa ions we e de i ed in o de o ha e a s able algo i hm. Ini ial alues o co a iance ma ices o he gene ic si e effec and animal effec we e diagonal ma ices ha ing alues diag{100,10,1}. The ini ial alue o he esidual a iance was 100. Addi ionally, a iance componen s we e epa ame ized o compu a ional easons, because he a iance componen κwas close o ze o. Con e gence was imp o ed by scaling he ime be o e e e y ound o he EBLUP-algo i hm. Each ime o measu emen ij was mul iplied by a scaling ac o c,whichwas se equal o he mos ecen es ima e o κ. Consequen ly, he a iance compo- nen es ima e o he scaled pa ame e κ∗was 1 c2Va (κ), and hus la ge han he o iginal pa ame e κwhen c<1. Con e gence o he EBLUP-algo i hm was assumed when he ela i e ound o ound change was less han 10−3. Fu he mo e, wi hin e e y i e a ion o he EBLUP-algo i hm, he loca ion pa ame e s we e i e a ed un il he ela i e diffe ence be ween igh -hand and le -hand sides o he MME was less han 1×10−7. Co a iance componen es ima es we e calcula ed by he Expec a ion Maximiza ion (EM) -algo i hm, and con e gence was assumed when he ound o ound change was less han 5 ×10−7. The esul s a e om 50 simula ion eplica es. Rela i e bias, ela i e s anda d de ia ion (Rel. SD) and ela i e mean squa ed e o (Rel. MSE), as pe cen age om he ue alue, we e calcula ed o he diffe ence o wo le els o ixed effec and o he a iance componen pa ame e es ima es. The ela i e bias 358 K. Vuo i e al. He e Z∗=∂ ∂uT|u=˜u and ˜ uis he empi ical BLUP-es ima e o he an- dom effec s. In addi ion, he linea e m in he expansion anishes, be- cause he i s de i a i e o he unc ion a ML-solu ions is ze o. Also (y− (X,b,Z,˜ u))TR−1 (X,b,Z,˜ u) is assumed o be negligible, because he esidual ec o (y− (X,b,Z,˜ u))TR−1has mean ze o. Now, app oxima ion o he likelihood unc ion Lis L∗(b,θ|y)=(2π)−n 2|R|−1 2(2π)−l+q 2|D|−1 2 exp −1 2(y− (X,b,Z,˜ u))TR−1(y− (X,b,Z,˜ u)) −1 2˜ uTD−1˜ u−1 2(u−˜ u)TZ∗R−1Z∗+D−1(u−˜ u)du =(2π)−n 2|R|−1 2|D|−1 2Z∗R−1Z∗+D−1 −1 2 exp −1 2(y− (X,b,Z,˜ u))TR−1(y− (X,b,Z,˜ u)) −1 2˜ uTD−1˜ u (2π)−l+q 2Z∗R−1Z∗+D−1 1 2 ×exp −1 2(u−˜ u)TZ∗R−1Z∗+D−1(u−˜ u)du =exp −1 2nln(2π)−1 2ln |R|−1 2ln |D|−1 2ln Z∗R−1Z∗+D−1 −1 2(y− (X,b,Z,˜ u))TR−1(y− (X,b,Z,˜ u)) −1 2˜ uTD−1˜ u and he loga i hm o L∗(b,θ|y)is l∗(b,θ|y)=−1 2nln (2π)−1 2ln(|R||I+Z∗TR−1Z∗D|) −1 2(y− (X,b,Z,˜ u))TR−1(y− (X,b,Z,˜ u)) −1 2˜ uTD−1˜ u.