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

Genet. Sel. Evol. 38 (2006) 343–358 343 c INRA, EDP Sciences, 2006 DOI: 10.1051/gse:2006008 Original article Estimation of non-linear growth models by linearization: a simulation study using a Gompertz function Kaarina V∗,IsmoS´ , Marja-Liisa S´ -A, Esa A. M¨  MTT Agrifood Research Finland, Biotechnology and Food Research, Biometrical Genetics, FIN-31600 Jokioinen, Finland (Received 6 July 2005; accepted 27 January 2006) Abstract – A method based on Taylor series expansion for estimation of location parameters and variance components of non-linear mixed effects models was considered. An attractive property of the method is the opportunity for an easily implemented algorithm. Estimation of non-linear mixed effects models can be done by common methods for linear mixed effects models, and thus existing programs can be used after small modifications. The applicability of this algorithm in animal breeding was studied with simulation using a Gompertz function growth model in pigs. Two growth data sets were analyzed: a full set containing observations from the entire growing period, and a truncated time trajectory set containing animals slaughtered prematurely, which is common in pig breeding. The results from the 50 simulation replicates with full data set indicate that the linearization approach was capable of estimating the original parameters satisfactorily. However, estimation of the parameters related to adult weight becomes unstable in the case of a truncated data set. Gompertz function /non-linear mixed effects /variance components /breeding values / likelihood approximation 1. INTRODUCTION Non-linear functions are particularly suited to model growth data, because predictions outside the data range can be made more reliably than by linear models, and the entire growth process can be described by few parameters. For example, growth data models commonly apply the Gompertz function, where the estimated parameters can have biological meaning. Non-linear models are, however, more complicated to solve than linear models, and several algorithms ∗Corresponding author: [email protected] Article published by EDP Sciences and available at http://www.edpsciences.org/gse or http://dx.doi.org/10.1051/gse:2006008 344 K. Vuori et al. have been proposed to estimate the parameters and variance components of non-linear mixed effects models [4]. In animal production research, the Bayesian framework has received much attention in growth curve analysis [2, 11]. This popularity is due to the use of Markov chain Monte Carlo methods that allow the solution of numerically complicated posterior density integration and calculation of confidence interval estimates. The cost of this procedure is, however, intensive calculations and the need to assure a sampling equilibrium [3]. Another possibility is to approximate the likelihood function using linearization [4, 10, 17, 18] or numerical integration [10]. Also, both of these alternatives are computationally difficult, yet, linearization may enable the inference of linear mixed effects models. All linearization methods in the literature are quite similar. Early methods use first-order Taylor series expansion of non-linear functions around expectation of the random effects, and are solved by either maximum likelihood (ML) or generalized least squares (GLS) estimation [4]. Lindstrom and Bates [9] suggested a more accurate method of making the expansion around current estimates of the random effects. Subsequent research has focused more on the second-order Taylor series expansion of integrals invoked by the Laplacian approximation [10, 17, 18]. Because of generality and familiar formulation, the most interesting choice of approximation is based on the second-order Taylor series expansion with respect to random effects that were presented by Wolfinger and Lin [18]. They gave two alternative approaches to select points of expansion: a zeroexpansion method using expected values, and an EBLUP-expansion method using the empirical best linear unbiased predictors of the random effects. Both approximations lead to algorithms that iteratively fit mixed linear models to the suitably transformed data using either ML or restricted maximum likelihood (REML). Therefore, they allow the use of commonly applied methods for linear mixed effects models, and the use of existing programs after small modifications. A similar algorithm was proposed by Breslow and Clayton [3] in the context of generalized linear mixed models. Because conditions to functionality of the approximation methods are difficult to identify, Wolfinger and Lin [18] recommended simulation studies for assessing the performance of the methods in diverse kinds of non-linear models and data sets. The aim of this work was to describe and examine the performance of the EBLUP-expansion method for the Gompertz function applied to the analysis of growth in the pig through simulation. The EBLUP-expansion is recommended especially for cases where the variance components are large, which may be the case for an adult weight parameter for pigs. Also, Lindstrom and Bates [9] Estimation of non-linear growth models 345 suggested that the expansion around the expected zero value may lead to poor estimates when substantial inter-individual variation exists. We chose to examine the method through the analysis of two data sets. The first analysis tested general performance of the EBLUP-expansion technique, and the second analysis tested performance of the method for incomplete data. Incomplete data are common in pig production, where the adult weight is unavailable due to an earlier slaughter age. 2. MATERIALS AND METHODS 2.1. Simulations The Gompertz function has been shown to fit pig growth data, such as live weight and protein retention, well [14–16]. We assumed that weights of an individual ifollowed the Gompertz model: yij =αexp(−βexp(−κtij)) +eij,j=1,...,ni where niis the number of observations for individual i,yij is the observed weight at age tij (in days), α,βand κare the parameters of the Gompertz function, and eij is the random residual. The parameters have biological meaning: αis the adult weight, κis the rate of exponential decay of the initial growth rate, and βis the logarithm of the ratio of birth weight to adult weight. Each of the parameters α,βand κcan be described by a linear mixed effects model. In this study, we will consider a sire model, although notation could be for an animal model. The full model for observation jof animal iis yij =(xαijbα+zs,αisα+zp,αipα) exp(−(xβijbβ+zs,βisβ+zp,βipβ) exp(−(xκijbκ+zs,κisκ+zp,κipκ)tij)) +eij,(1) where (bα,bβ,bκ)T=bis a d×1-vector of fixed effects, (sα,sβ,sκ)T=sis a l×1-vector of random additive genetic sire effects and (pα,pβ,pκ)T=pis a q×1-vector of random animal effects other than sire. Vectors x,zsand zpare from the design matrices of fixed, random sire and random animal effects X, Zsand Zp, respectively. It is assumed that           s p e           ∼N                    0 0 0           ,          G00 0P0 00R                     . 346 K. Vuori et al. Here, G=G0⊗A,whereAis a matrix of additive relationships between sires and G0isa3×3 genetic covariance matrix for the Gompertz parameters. Similarly, P=P0⊗Iq,whereP0is a 3 ×3 covariance matrix, i.e. the random animal effects pwere identically and independently distributed for the animals. Furthermore, the residuals were assumed to be independently distributed and homoscedastic, e∼N(0,Inσ2 e). The Gompertz function coefficients were generated to simulate pig growth. The model had one fixed effect with two levels, and the random effects were the genetic sire effect and the animal effect other than the sire. The first fixed effect level had values 210, 5 and 0.017 for the parameters α,βand κ, respectively. The other fixed effect level had values 220, 4.7 and 0.016 for α,βand κ, respectively. The random effects were assumed to be normally distributed with mean zero and block diagonal covariance matrices. Variance and covariance components in the matrices for the genetic sire effect and for the animal effect are shown in Tables I and II. The residual variance σ2 ewas one. These parameters approximated the variances calculated by the NLMIXED procedure of the SAS program that fitted the Gompertz model to growth performance data of Finnish pigs [12]. Simulation of the random sire effect required taking into account the pedigree. The pedigree had three generations of animals with 10 unrelated founder grandsires. Each of the 10 grandsires was mated with 20 unrelated dams that produced one son each, i.e., 20 half-sibs. The half-sibs were mated with unrelated dams to produce 24 progeny per sire. Only the last generation of animals had records. Thus, the data included 4800 tested animals. Two data sets were made: a complete set, and a truncated time trajectory set. The complete data contained 30 equally-distanced observations per animal between 50 and 253 days. The truncated time trajectory data contained slaughter weights up to 115 kg, which is similar to the common slaughter weight in pigs, and occurs at about 120 days of age. Consequently, the number of observations was reduced from 30 to about 11 per animal, i.e., almost two thirds of the data were discarded. 2.2. Method to estimate the values of the growth parameters The non-linear mixed model considered was y=f(X,b,Zs,s,Zp,p)+e, where yis an n×1-vector of observations, fis the Gompertz function, and e is an n×1-vector of random residuals. Vectors b,sand p, with matrices X, Zs Estimation of non-linear growth models 347 Table I. Relative bias, relative standard deviation (Rel. SD) and relative mean squared error (Rel. MSE) (as percent from the true value) for (co)variance components of genetic sire effects from the 50 replicates of full and truncated time trajectory data. Subscripts α,βand κdenote the three parameters in the Gompertz function. Full data Truncated data Parameter True 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. Relative bias, relative standard deviation (Rel. SD) and relative mean squared error (Rel. MSE) (as percent from the true value) for (co)variance components of the animal effect from the 50 replicates of full and truncated time trajectory data. Subscripts α, βand κdenote the three parameters in the Gompertz function. Full data Truncated data Parameter True 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. Vuori et al. and Zp, were defined as before. For the random effects, denote uT=(sTpT), Z=ZsZpand D=G0⊗A0 0P 0⊗Iq. Now the distribution assumptions were u e∼N0 0,D0 0R . Although Ris diagonal here, any form is allowed, so the general form Rwill be used hereinafter. Unknown elements of covariance matrices G0,P0and R are denoted by parameter vector θ. The maximized likelihood function was L(b,θ|y)=(2π)−n 2|R|−1 2(2π)−l+q 2|D|−1 2 exp −1 2(y−f(X,b,Z,u))TR−1(y−f(X,b,Z,u)) −1 2uTD−1udu.(2) Only in some cases is the closed form of (2) found, so the integral is often solved numerically. However, numerical methods for the non-linear functions are usually slow to converge and numerically unstable. Instead, the integral may be approximated by quadratic Taylor-series expansion of the exponent. The second-order expansion was made about the EBLUP before integration of the likelihood function (see Appendix). This gave approximation to the logarithm of the likelihood function (2): l∗(b,θ|y)=−1 2nln (2π)−1 2ln(|R||I+Z∗TR−1Z∗D|) −1 2(y−f(X,b,Z,˜ u))TR−1(y−f(X,b,Z,˜ u)) −1 2˜ uTD−1˜ u,(3) where Z∗=∂f∂uT|u=˜u and ˜ uis the empirical BLUP-estimate of random effects. For the Gompertz function and two random effects in the model, Z∗had elements ∂f ∂αi =exp(−βiexp(−κitj)) cαi ∂f ∂βi =αiexp(−βiexp(−κitj)) (−exp(−κitj)) cβi(4) ∂f ∂κi =αiexp(−βiexp(−κitj)) (−βiexp(−κitj)) (−tj)cκi where cis zsor zpdepending on the random effect differentiated (see (1)). Estimation of non-linear growth models 349 Pinheiro and Bates [10] used the approximation (3) in estimation of parameters by the Laplacian approximation. However, no straightforward generalization to the REML-estimation was presented. Wolfinger and Lin [18] developed formula (3) further. Denote V=Z∗DZ∗T+R|u=˜u . Then, l∗(b,θ|y)=−1 2nln(2π)−1 2ln |V| −1 2(y−f(X,b,Z,˜ u)+Z∗˜ u)TV−1(y−f(X,b,Z,˜ u)+Z∗˜ u), where 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 to a similar estimation function for variance component estimation presented by Lindstrom and Bates [9], although through different derivation. 2.2.1. Estimation of the fixed and random effects Assume that the variance component vector θis known. Maximum likelihood estimation for the parameters band uleads to solving equations: X∗TR−1(y−f(X,˜ b,Z,˜ u)) =0 Z∗TR−1(y−f(X,˜ b,Z,˜ u)) =D−1˜ u,(5) where X∗=∂f∂bTb=˜ band ˜ bis the estimate of fixed effects b.Elements of X∗are similar to Z∗, except that coefficient cin (4) is replaced by xdue to the differentiated fixed effect. However, in order to arrive to these simple equations, dependency of Von bthrough Z∗has to be ignored. On the basis of arguments made by Bates and Watts [1], Wolfinger and Lin [18] justified this by appealing to intrinsic non-linearity instead of non-linearity of the parameters. Denote Y=y−f(X,˜ b,Z,˜ u)+X∗˜ b+Z∗˜ u. Equations (5) can now be written 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 similar to the mixed model equations (MME) for the linear models. Thus, already established methods for solving linear models can be used to analyse the pseudo-data Ycreated from the original data ywith ˜ band ˜ uequal to their most recent estimates. 350 K. Vuori et al. 2.2.2. Estimation of the variance components After finding estimates of the location parameters, profile likelihood can be used to estimate the variance components by setting b=˜ b(θ). The logarithmic likelihood function of the parameter vector θcan be written with the pseudodata as l∗ ML(θ)=−1 2nln(2π)−1 2ln |V|−1 2(Y−X∗˜ b)TV−1(Y−X∗˜ b).(7) Differentiating equation (7) with respect to θgives −1 2tr V-1 ∂V ∂θj+1 2(Y−X∗˜ b)TV-1 ∂V ∂θj V-1(Y−X∗˜ b).(8) Maximum likelihood estimates of variance components are found by equating (8) to zero and solving for θ. Instead of the ML-estimates, REML-estimates are commonly used in practise. These estimates account for losses in degrees of freedom caused by the estimation of fixed effects b[5]. The logarithmic likelihood function 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) Differentiation with respect to θand equating to zero gives −1 2tr P∂V ∂θj+1 2(Y−X∗˜ b)TV-1 ∂V ∂θj V-1(Y−X∗˜ b)=0,(10) where P=V−1−V−1X∗(X∗TV−1X∗)−1X∗TV−1. Solutions in θare REMLestimates of variance components. 2.2.3. EBLUP-algorithm The approximate ML-solutions of location parameters and variance components can be obtained by iteratively solving the equations (6) and (8) until convergence. Correspondingly, the REML-solutions for the EBLUP-expansion are obtained by iteratively solving the equations (6) and (10). Hence, the algorithm fits the linear mixed effects model Y=X∗b+Z∗u+efor the pseudo-data Y and the working vectors X∗and Z∗,whereu∼N(0,G(θ)) and e∼N(0,R(θ)). Estimation of non-linear growth models 351 2.3. Implementation We chose to implement the REML-based EBLUP-algorithm, because programs to solve the linear mixed effects models are available and commonly used by animal breeders. MiX99 [13] was used to solve the mixed model equations (6), and DMU [6], modified for random regression by Kettunen et al. [7], was used to solve the REML estimates of covariance components (10). The capability of fitting fixed and random regression models is crucial for implementation, because the coefficients in X∗and Z∗can have any values. Implementation of the linearization procedure required the Gompertz function formulas to be included in MiX99. However, there was no need to make changes to the variance component estimation program. Starting values for both the location parameter effects and the variance components had to be assigned before first iteration. A natural choice was to initialize random effects with the expected value zero. However, initial values for fixed effects were derived with the model function of growth curve and available data. When the Gompertz model is used, only the asymptotic weight parameter has a natural initial value, which is the maximum value of the dependent variable. In the other cases, complex equations were derived in order to have a stable algorithm. Initial values for covariance matrices of the genetic sire effect and animal effect were diagonal matrices having values diag{100,10,1}. The initial value for the residual variance was 100. Additionally, variance components were reparametrized for computational reasons, because the variance component κwas close to zero. Convergence was improved by scaling the time before every round of the EBLUP-algorithm. Each time of measurement tij was multiplied by a scaling factor c,whichwas set equal to the most recent estimate of κ. Consequently, the variance component estimate for the scaled parameter κ∗was 1 c2Var(κ), and thus larger than the original parameter κwhen c<1. Convergence of the EBLUP-algorithm was assumed when the relative round to round change was less than 10−3. Furthermore, within every iteration of the EBLUP-algorithm, the location parameters were iterated until the relative difference between right-hand and left-hand sides of the MME was less than 1×10−7. Covariance component estimates were calculated by the Expectation Maximization (EM) -algorithm, and convergence was assumed when the round to round change was less than 5 ×10−7. The results are from 50 simulation replicates. Relative bias, relative standard deviation (Rel. SD) and relative mean squared error (Rel. MSE), as percentage from the true value, were calculated for the difference of two levels of fixed effect and for the variance component parameter estimates. The relative bias 358 K. Vuori et al. Here Z∗=∂f∂uT|u=˜u and ˜ uis the empirical BLUP-estimate of the random effects. In addition, the linear term in the expansion vanishes, because the first derivative of the function at ML-solutions is zero. Also (y−f(X,b,Z,˜ u))TR−1f(X,b,Z,˜ u) is assumed to be negligible, because the residual vector (y−f(X,b,Z,˜ u))TR−1has mean zero. Now, approximation for the likelihood function Lis L∗(b,θ|y)=(2π)−n 2|R|−1 2(2π)−l+q 2|D|−1 2 exp −1 2(y−f(X,b,Z,˜ u))TR−1(y−f(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−f(X,b,Z,˜ u))TR−1(y−f(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−f(X,b,Z,˜ u))TR−1(y−f(X,b,Z,˜ u)) −1 2˜ uTD−1˜ u and the logarithm of L∗(b,θ|y)is l∗(b,θ|y)=−1 2nln (2π)−1 2ln(|R||I+Z∗TR−1Z∗D|) −1 2(y−f(X,b,Z,˜ u))TR−1(y−f(X,b,Z,˜ u)) −1 2˜ uTD−1˜ u.