Non-linear mixed models in the study of growth of naturalized chickens

This study was conducted to examine the inclusion of random effects in non-linear models, identify the most suitable models, and describe the growth of naturalized chickens. Live-weight records of 166 birds of the Graúna Dourada, Nordestina, and Teresina ecotypes were estimated. The asymptotic weight (A), integration constant, related to animal initial weight (B), and the maturing rate (k) parameters of the non-linear Gompertz, Logistic, and von Bertalanffy models were estimated and adjusted using the Gauss-Newton method. Residual variance decreased by more than 50% when random effects were added to the model. The best fits in the estimate of the growth curve of females were obtained by associating the random effects with the three parameters of the Gompertz and Logistic models. The association of random effects with two parameters (asymptotic weight and maturing rate) and with the three parameters of the Logistic model provided the best fits for the males. The Teresina ecotype has the highest adult weight in both sexes, despite its slower growth. The opposite is true for the Graúna Dourada ecotype, formed by lighter and earliergrowing animals. The inclusion of random effects in models provides greater accuracy in the estimate of the growth curve.


Introduction
Chickens were introduced in Brazil at the time of its discovery and colonization. The chicken groups that were not subjected to any breeding method and that adapted to the rearing conditions and to the environment in which they were managed were named "naturalized" chickens.
Several ecotypes of the species were extinguished after the introduction of genetically improved breeds and lines from other countries. Some groups are being subjected to genetic conservation in teaching and research institutions. However, little information exists on the production rates and growth of those animals.
Non-linear models allow for a comparison of the growth curve of different genetic groups, making it possible to evaluate differences in animal growth caused by sex, management, and rearing environment. They also provide essential information to guide the sustainable preservation of animals at risk of extinction, such as the estimate of nutritional requirements and growth (Hruby et al., 1994;Selvaggi et al., 2015).
In studies on growth curves, some basic principles should be taken into account for the efficient use of non-linear models. One of such principles is that the data must present homogeneity of variance and residuals must not be correlated.
Nevertheless, longitudinal data derived from growth studies may exhibit different variances throughout the animal's life. Moreover, repeated measures from the same individual provide correlated residuals, compromising the efficient use of these models, which are assumed to be fixed (Guedes et al., 2004;Mazucheli et al., 2011).
An alternative to address this problem is the incorporation of random effects associated with the individuals into the model, thus characterizing it as a non-linear mixed model. The use of this type of model makes it possible to adjust flexible covariance structures capable of handling imbalanced data (Lindstrom and Bates, 1990).
Non-linear mixed models have been applied in studies on the growth pattern of quail (Karaman et al., 2013) and birds (Sofaer et al., 2013), to estimate metabolizable energy utilization by chickens (Romero et al., 2009), estimate the nutritional requirements of laying hens (Strathe et al., 2011), among others. However, there are no reports of the use of methodologies to describe the growth pattern of naturalized chickens.
Therefore, this study proposes to evaluate the inclusion of random effects in non-linear mixed models to describe the growth curve of naturalized chickens.
The linear models can be described as follows: in which y ij is the j-th observation of individual I, f is a non-linear function of a parameter vector β i and a vector x ij of prediction of y ij , and ε ij is a normal-distribution error. The parameter vector can vary between the individuals and is incorporated into the model as shown below: in which β is a vector of parameters of fixed effects associated with the population; b i is a vector of random effects associated with individuals i; A i and B i are incidence matrices of the fixed and random effects, respectively; and σ 2 D is a matrix of covariances between the random effects (Lindstrom and Bates, 1990).
To adjust the growth curve, the A, B, and k parameters of the Gompertz (Laird, 1965), Logistic (Nelder, 1961), and von Bertalanffy (Bertalanffy, 1957) models were estimated using the Gauss-Newton method, described by Hartley (1961) for non-linear models.
In the non-linear mixed models, the random effects are associated with the parameters. For the Gompertz, Logistic, and von Bertalanffy models, it is possible to compare seven possibilities of addition of random effects. The present study tested a model without random effects, named "nonlinear model of fixed effects"; three models with addition of only one random effect; three models with addition of two random effects; and one model with addition of random effect on the three parameters (Table 1).
The models were compared by using the following criteria for the choice of the model that best fit the growth curve: the mean squared error (MSE), calculated by dividing the sum of residual squares by the number of observations; the average absolute deviation of residuals (AAD), described by Sarmento et al. (2006), calculated as the sum of the module or absolute values of the residuals divided by the number of observations; and the coefficient of determination (R²), obtained from the calculation of the square of the correlation between observed and estimated weights.
Akaike's Information Criterion (AIC) (Akaike, 1974) and the Bayesian Information Criterion (BIC) (Schwarz, 1978) were also used for the choice of the best-fitting model. The values were obtained as follows: AIC = -2logL(θ̂) + 2 (p) and BIC = -2logL(θ̂) + ln(N)p, in which p represents the total number of parameters estimated by the model, N is the total number of observations, and logL(θ̂) is the logarithm of restricted likelihood. Ibiapina Neto et al. 4 Age (t i ) and weight (y i ) at the inflection point were calculated using the equations t i = (lnB)/k and y i = A/2 for the Logistic model; t i = (logB)/k and y i = A/e for the Gompertz model; and t i = loge(3b)/k and y i = 8A/27 for the von Bertalanffy model.
R software was used to estimate the parameters and to adjust the non-linear models, applying the nls (Nonlinear Least Squares) function of the Stats package and the nlme (Nonlinear Mixed-Effects Models) package. The maximum likelihood method was employed to estimate the mixed-model parameters, using the algorithm created by Lindstrom and Bates (1990) for integer approximation.
The same software was used to compare and group the similar models via Tocher's optimization method. For this, the values from the evaluation of fitting criteria of the adjusted models were subjected to the Tocher function in the BioTools package of R software (R Core Team, 2017). Pearson's correlations between the parameters were obtained using the Stats package of R software.

Results
The Gompertz model with random effects associated with the A and k parameters estimated similar asymptotic weights between the two sexes. However, in the other tested models, a significant difference was observed between males and females for this parameter (Table 2).
A and B y = (A + u1)/(1 + (B + u2)exp(-kt)) A and k y = (A + u1)/(1 + Bexp(-(k + u3)t)) B and k y = A/(1 + (B + u2)exp(-(k + u3)t)) A, B, and k y = (A + u1)/(1 + (B + u2)exp(-(k + u3)t)) y is body weight at age t; A is the asymptotic weight when t tends to plus infinite and is interpreted as weight at adult age; B is an integration constant, related to the initial weights of the animal; k is established by the initial values of y; t is interpreted as the time (age) when weight was measured; and u1, u2, and u3 are the random effects.
Non-linear mixed models in the study of growth of naturalized chickens Ibiapina Neto et al.

5
Growth rate also did not differ significantly between the sexes in the Logistic model when the random effects were associated with the k parameter and with the B and k parameters. The Gompertz and von Bertalanffy models, in turn, with random effects associated with the B and k parameters and A, B, and k parameters, showed higher growth rates in the females. The same result was obtained when the random effects were associated with the A and B parameters of the Logistic model. Males showed the highest growth rates in the other model variations.
The lowest MSE, AAD, AIC, and BIC values of the models that described the growth of females were 2972, 40.09, 7024, and 7067, respectively, for the Logistic model, and 3339, 42.88, 7069, and 7113 for the Gompertz model, both with the random effects associated with the three parameters of the curve (Table 3). For males, the lowest values of the evaluation criteria were obtained with the Logistic model with random effects associated with the three model parameters and with the A and k parameters.
When cluster analysis was performed applying Tocher's optimization method, the models were divided into three similar groups for females and into 10 groups for males (Tables 4 and 5), using the results  Table 3. Thus, the models of the group that showed the lowest values of the evaluation criteria were used to describe the growth of the birds (Table 6).
The Logistic model with random effects associated with the three model parameters was the most suitable to describe the growth of Graúna Dourada and Teresina females, as it showed the lowest values for the MSE, AAD, AIC, and BIC selection criteria. The Gompertz model, in turn, also with random effects associated with the three parameters, was the model that showed the lowest values of the above-mentioned parameters (3564.84, 47.25, 1417.45, and 1445.08, respectively), using the data of Nordestina females.
The Logistic model with random effects associated with the three model parameters also showed the lowest values for the MSE, AAD, and AIC selection criteria in describing the growth of Graúna Dourada and Nordestina males. The same model, now with random effects associated with the A and k parameters, was the one that best described the growth of Teresina males, as it showed the lowest AIC (1758.97) and BIC (1779.67) values. The other selection criteria presented values similar to those of the model with random effects associated with the three parameters. Non-linear mixed models in the study of growth of naturalized chickens Ibiapina Neto et al. Teresina chickens showed the highest asymptotic weight in both sexes, while the Graúna Dourada ecotype exhibited the lowest asymptotic weights for the males (Table 6). Adult weight in the females of the Graúna Dourada and Nordestina ecotypes did not differ significantly. The Graúna Dourada ecotype showed the highest initial weights and the highest growth rates.
The scales of the residuals scatterplot ( Figure 1) revealed that the Gompertz model had a wider range of error distribution in describing the growth of Graúna Dourada females. A similar result was seen for the Nordestina males using the Logistic model with random effect associated with the A and k parameters.
Graúna Dourada females reached the asymptotic weight at an age close to 150 days, whereas the Nordestina and Teresina females attained adult weight at approximately 200 days of age (Figures 2  and 4). The same was observed for the males (Figure 3). However, some animals of the abovementioned ecotypes did not follow the behavior shown by the group in which they were clustered.
Negative correlation coefficients were obtained between the A and B and A and k parameters in all evaluated models, for both sexes ( Table 7). The correlation coefficients obtained between the B and k parameters were positive in the evaluated models and for both sexes.  highest AIC values, which were higher than 7327. The second group was formed by the models with AIC values ranging from 7232 and 7152. The third group was formed by the Logistic and Gompertz models with random effect associated with the three parameters, which showed the lowest AIC values (7023.57 and 7069.48, respectively).
A similar result was observed for the clustering considering the data of males (Table 5). Group 1 contained the models with the highest AIC values, followed by groups 2, 6, 9, 7, 3, 5, 4. Lastly, group 8 included the models with the lowest AIC values.
It should be stressed that cluster analysis was processed using the values of the MSE, AAD, R 2 , AIC, and BIC model selection criteria, which were all considered in the formation of the groups. This is clearly observed as the AIC value of the model allocated in group 10 was between the values of the two models of group 9.
Despite showing similar values for this criterion, the model of group 10 was allocated in another group because it showed different MSE and AAD values when compared with the models of group 9.  There was a similarity in the dispersion of residuals in the scatterplots generated with the different models. The errors observed at the beginning of the curve had a wider range for the Logistic model with the addition of three random effects, for both sexes. They were also lower at the beginning of the curve and tended to increase with age.
Residual variance estimated by MSE decreased with the addition of random effects to the model, for both sexes. Similar results showing a reduction in residual variance after the addition of random effects to the model were found by Aggrey (2009) and Karaman et al. (2013). It can thus be affirmed that the inclusion of random effects in the model may generate more-reliable estimates compared with the nonlinear models of fixed effects, although the introduction of the random effect on parameter B did not lead to improvements in the estimate of the parameters.
By introducing the random effect, pointing out the difference between the individuals, individual curves are generated, as different individuals also grow differently.

Conclusions
The Teresina ecotype has the highest asymptotic weights for both sexes and is the group of slower growth when compared with the Graúna Dourada ecotype. The latter, in turn, is formed by lighter and earlier-growing animals. The inclusion of random effects in the Logistic and Gompertz models provides greater accuracy in the estimate of the growth curve.