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

ABSTRACT

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 earlier-growing animals. The inclusion of random effects in models provides greater accuracy in the estimate of the growth curve.

Keywords:
Gompertz; modelling; phenotypic variation; rooster

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.

Material and Methods

The study was conducted after approval by the institutional Animal Use Committee (case no. 404/17). A database was used with live-weight records (collected fortnightly) of 166 birds of the Graúna Dourada (53 females and 38 males), Nordestina (19 females and 17 males), and Teresina (16 females and 23 males; Figure 1) ecotypes, located in Teresina, PI - Brazil (5°03'57.2" S, 42°42'09.2" W), with a minimum number of five weights per bird.

Figure 1
Residuals standardized for the models that best described the growth of naturalized chicken ecotypes.

The chicks were weighed at birth and then housed in proper metal cages equipped with a feeder, a drinker, and a heating source for the first weeks of life. The animals were divided into four age ranges, in which the supplied feed would meet their respective requirements. The diet provided in the pre-starter phase, from the 1st to the 30th day of life, was composed of 630 g/kg corn, 320 g/kg soybean meal, and 50 g/kg vitamin-mineral mix (Núcleo Fit Aves Pré-Inicial^®, Poli Nutri, Brazil). In the starter phase, from the 30th to the 60th day of life, the feed composition was 660 g/kg corn, 290 g/kg soybean meal, and 50 g/kg vitamin-mineral mix (Núcleo Fit Aves Inicial^®, Poli Nutri, Brazil). This was followed by the grower phase, from the 60th day of life until the appearance of the first signs of reproduction, which occurred at around 180 days of life (feed composition: 700 g/kg corn, 250 g/kg soybean meal, and 50 g/kg vitamin-mineral mix [Núcleo Fit Aves Engorda^®, Poli Nutri, Brazil]). Lastly, in the finisher phase, which started on the 180th day of life, the birds received a diet composed of 630 g/kg corn, 245 g/kg soybean meal, 85 g/kg calcitic limestone, and 40 g/kg vitamin-mineral mix (Núcleo Poli Macro Ovo 4%^®, Poli Nutri, Brazil). The pre-starter, starter, grower, and finisher diets had approximate protein (g/kg) and energy (MJ/kg) contents of 195 and 11.9, 185 and 12.1, 170 and 12.3, and 160 and 11.3, respectively.

The linear models can be described as follows:

$$y_{ij}=f\left({\beta }_i,x_{ij}\right)+{\epsilon }_{ij},{\epsilon }_{ij}~N\left(0,{\sigma }_e^2\right),$$

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:

$${\beta }_i=A_i\beta +B_i{b}_i,b_i~N\left(0,{\sigma }^{2}D\right),$$

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 σ²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 “non-linear 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).

Table 1
Models without and with random effect proposed to explain the growth pattern of naturalized chickens

Model	Random effect	Formula
von Bertalanffy	Absent	$\text{y}=\text{A}{(1-\text{Bexp}(-\text{kt}))}^3$
	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-\text{Bexp}(-\text{kt}))}^3$
	B	$\mathrm{y}=\mathrm{A}{(1-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))}^3$
	k	$\mathrm{y}=\mathrm{A}{(1-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))}^3$
	A and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
	B and k	$\mathrm{y}=\mathrm{A}{(1-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
	A, B, and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t})]}^3$
Gompertz	Absent	$\mathrm{y}=\mathrm{Aexp}(-\mathrm{Bexp}(-\text{kt}))$
	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-\mathrm{Bexp}(-\text{kt}))$
	B	$\mathrm{y}=\mathrm{Aexp}(-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	k	$\mathrm{y}=\mathrm{Aexp}(-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	A and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	B and k	$\mathrm{y}=\mathrm{Aexp}(-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	A, B, and k	$\text{y}=(\text{A}+\text{u1})\text{exp}(-(\text{B}+\text{u}2)\text{exp}(-(\text{k}+\text{u}3)\text{t}))$
Logistic	Absent	$\mathrm{y}=\mathrm{A}/(1+\mathrm{Bexp}(-\text{kt}))$
	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+\mathrm{Bexp}(-\text{kt}))$
	B	$\mathrm{y}=\mathrm{A}/(1+(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	k	$\mathrm{y}=\mathrm{A}/(1+\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	A and k	$\text{y}=(\text{A}+\text{u}1)/(1+\mathrm{Bexp}(-(\text{k}+\text{u}3)\text{t}))$
	B and k	$\mathrm{y}=\mathrm{A}/(1+(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	A, B, and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{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.

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(\widehat{\theta })+2(p)$ and $BIC=-2logL(\widehat{\theta })+\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(\widehat{\theta })$ is the logarithm of restricted likelihood.

Age (t_i) and weight (y_i) at the inflection point were calculated using the equations ${\text{t}}_{\text{i}}\text{=}(\text{lnB})/\text{k}$ and ${\text{y}}_{\text{i}}=\text{A}/2$ for the Logistic model; ${\text{t}}_{\text{i}}\text{=}(\text{logB})/\text{k}$ and ${\text{y}}_{\text{i}}=\text{A}/\text{e}$ for the Gompertz model; and ${\text{t}}_{\text{i}}=\text{loge}(3\text{b})/\text{k}$ and ${\text{y}}_{\text{i}}=8\text{A}/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).

Table 2
Estimates of the parameters and inflection point (age and weight) of non-linear models of fixed effects and non-linear models of mixed effects

Random effect in:	Gompertz model
	A		B		k		Age		Weight
	M	F	M	F	M	F	M	F	M	F
Absent	1,438a	1,261b	6.52a	5.70a	0.0214a	0.0197b	38	38	529	464
A	1,482a	1,269b	6.07a	5.25b	0.0201a	0.0186b	39	39	545	467
B	1,437a	1,261b	6.51a	5.96b	0.0214a	0.0197b	38	39	528	464
k	1,809a	1,532b	4.54a	4.29b	0.0147a	0.0141b	45	45	665	563
A and B	1,482a	1,269b	6.07a	5.25b	0.0201a	0.0186b	39	39	545	467
A and k	1,546a	1,532a	5.58a	4.29b	0.0186a	0.0141b	40	45	568	563
B and k	1,864a	1,261b	4.46b	5.69a	0.0141b	0.0197a	46	38	685	464
A, B, and k	1,808a	1,385b	4.54b	4.85a	0.0147b	0.0167a	45	41	665	509
Logistic model
Absent	1,310a	1,129b	40.48a	35.47b	0.0367a	0.0354b	101	101	655	565
A	1,335a	1,126b	36.19a	30.59b	0.0349a	0.0335b	103	102	668	563
B	1,310a	1,128b	40.47a	35.46b	0.0367a	0.0354b	101	101	655	564
k	1,465a	1,226b	27.42a	25.00b	0.0299a	0.0294a	111	109	733	613
A and B	1,336a	1,226b	36.18a	30.60b	0.0349b	0.0355a	103	96	668	613
A and k	1,362a	1,226b	39.78a	25.01b	0.0356a	0.0294b	103	109	681	613
B and k	1,465a	1,226b	27.43a	25.02b	0.0299a	0.0294a	111	110	733	613
A, B, and k	1,364a	1,152b	38.75a	29.25b	0.0352a	0.0322b	104	105	682	576
von Bertalanffy model
Absent	1,570a	1,365b	1.16a	1.09b	0.0155a	0.0147b	98	97	465	344
A	1,569a	1,365b	1.16a	1.09b	0.0155a	0.0147b	98	97	465	344
B	1,569a	1,365b	1.16a	1.09b	0.0155a	0.0147b	98	97	465	344
k	2,430a	2,071b	0.84a	0.81b	0.0084a	0.0081b	130	130	720	249
A and B	1,569a	1,413b	1.16a	1.01b	0.0155a	0.0133b	98	99	465	344
A and k	2,680a	2,071b	1.08a	0.81b	0.0139a	0.0081b	101	130	794	320
B and k	2,430a	1,365b	0.84b	1.09a	0.0084b	0.0147a	130	97	720	249
A, B, and k	1,670a	1,356b	1.08a	1.05b	0.0139b	0.0143a	101	96	495	320

• M - males; F - females.
• Means followed by different letters in the sexes differ from each other according to Duncan's test at the 5% probability level.

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.

Table 3
Evaluation criteria of the adjustment of non-linear models of fixed effects and non-linear models of mixed effects for males and females

Random effect in:	Gompertz model
	Females					Males
	MSE	AAD	R2	AIC	BIC	MSE	AAD	R2	AIC	BIC
Absent	12653	83.42	0.956	7354	7371	21028	106.9	0.945	6749	6766
A	6220	58.41	0.975	7170	7192	8408	69.57	0.983	6514	6536
B	12653	83.42	0.956	7356	7378	21028	106.9	0.945	6751	6772
k	6214	55.71	0.956	7181	7203	8774	67.23	0.959	6540	6561
A and B	6220	58.41	0.975	7174	7204	8409	69.57	0.982	6518	6548
A and k	6215	55.71	0.956	7185	7216	8200	67.79	0.981	6510	6540
B and k	12653	83.42	0.956	7360	7390	7635	63.19	0.954	6541	6571
A, B, and k	3339	42.88	0.988	7069	7113	8771	67.22	0.959	6550	6592
Logistic model
Absent	12105	78.16	0.923	7327	7345	20226	100.3	0.910	6729	6746
A	6096	55.41	0.941	7152	7174	8012	64.42	0.946	6490	6512
B	12105	78.16	0.923	7329	7351	20226	100.3	0.910	6731	6752
k	6138	57.13	0.930	7167	7189	8448	68.44	0.931	6519	6540
A and B	6096	55.41	0.941	7156	7187	8011	64.42	0.946	6494	6524
A and k	6138	57.12	0.930	7171	7202	4940	49.72	0.969	6428	6458
B and k	6137	57.12	0.930	7171	7202	8447	68.43	0.931	6523	6553
A, B, and k	2972	40.09	0.947	7024	7067	4916	49.55	0.967	6433	6476
von Bertalanffy model
Absent	13060	86.34	0.963	7372	7390	21623	109.7	0.944	6764	6781
A	13060	86.34	0.963	7374	7396	21623	109.7	0.944	6766	6787
B	13060	86.34	0.963	7374	7396	21623	109.7	0.944	6766	6787
k	6746	59.59	0.957	7228	7250	9662	71.76	0.959	6587	6608
A and B	6380	60.70	0.987	7190	7221	21623	109.7	0.944	6770	6800
A and k	6747	59.59	0.957	7232	7263	8589	70.60	0.991	6532	6562
B and k	13060	86.34	0.963	7378	7409	9663	71.76	0.959	6591	6621
A, B, and k	6486	61.34	0.988	7193	7215	8587	70.60	0.991	6538	6580

• MSE - mean squared error; AAD - average absolute deviation of residuals; AIC - Akaike Information Criterion; BIC - Bayesian Information Criterion.

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 described in 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).

Table 4
Groups established by Tocher's optimization method of different non-linear models to describe the growth of female chickens

Group	Model	Random effect	Formula
1	Gompertz	Absent	$\mathrm{y}=\mathrm{Aexp}(-\mathrm{Bexp}(-\text{kt}))$
	Gompertz	B	$\mathrm{y}=\mathrm{Aexp}(-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	Gompertz	B and k	$\mathrm{y}=\mathrm{Aexp}(-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Logistic	Absent	$\mathrm{y}=\mathrm{A}/(1+\text{Bexp}(-\text{kt}))$
	Logistic	B	$\mathrm{y}=\mathrm{A}/(1+(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	Bertalanffy	Absent	$\mathrm{y}=\mathrm{A}{(1-\mathrm{Bexp}(-\text{kt}))}^3$
	Bertalanffy	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-\mathrm{Bexp}(-\text{kt}))}^3$
	Bertalanffy	B	$\mathrm{y}=\mathrm{A}{(1-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))}^3$
	Bertalanffy	B and k	$\mathrm{y}=\mathrm{A}{(1-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
2	Gompertz	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-\mathrm{Bexp}(-\text{kt}))$
	Gompertz	k	$\mathrm{y}=\mathrm{Aexp}(-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Gompertz	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	Gompertz	A and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Logistic	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+\mathrm{Bexp}(-\text{kt}))$
	Logistic	k	$\mathrm{y}=\mathrm{A}/(1+\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Logistic	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	Logistic	A and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Logistic	B and k	$\mathrm{y}=\mathrm{A}/(1+(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Bertalanffy	k	$\mathrm{y}=\mathrm{A}{(1-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
	Bertalanffy	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))}^3$
	Bertalanffy	A and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
	Bertalanffy	A, B, and k	$\text{y}=(\text{A}+\text{u}1){(1-(\text{B}+\text{u}2)\exp (-(\text{k}+\text{u}3)\text{t}))}^3$
3	Gompertz	A, B, and k	$\text{y}=(\text{A}+\text{u}1)\exp (-(\text{B}+\text{u}2)\exp (-(\text{k}+\text{u}3)\text{t}))$
	Logistic	A, B, and k	$\text{y}=(\text{A}+\text{u}1)/(1+(\text{B}+\text{u}2)\exp (-(\text{k}+\text{u}3)\text{t}))$

Table 5
Groups established by Tocher's optimization method for different non-linear models to describe the growth of males

Group	Model	Random effect in:	Formula
1	Gompertz	Absent	$\mathrm{y}=\mathrm{Aexp}(-\mathrm{Bexp}(-\text{kt}))$
	Gompertz	B	$\text{y=Aexp(-(B+u2)exp(-kt))}$
	Bertalanffy	Absent	$\mathrm{y}=\mathrm{A}{(1-\mathrm{Bexp}(-\text{kt}))}^3$
	Bertalanffy	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-\mathrm{Bexp}(-\text{kt}))}^3$
	Bertalanffy	B	$\mathrm{y}=\mathrm{A}{(1-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))}^3$
	Bertalanffy	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))}^3$
2	Logistic	Absent	$\mathrm{y}=\mathrm{A}/(1+\mathrm{Bexp}(-\text{kt}))$
	Logistic	B	$\mathrm{y}=\mathrm{A}/(1+(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
3	Logistic	k	$\mathrm{y}=\mathrm{A}/(1+\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Logistic	B and k	$\mathrm{y}=\mathrm{A}/(1+(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
4	Logistic	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+\mathrm{Bexp}(-\text{kt}))$
	Logistic	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
5	Gompertz	A	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-\mathrm{Bexp}(-\text{kt}))$
	Gompertz	A and B	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-(\mathrm{B}+\mathrm{u}2)\exp (-\text{kt}))$
	Gompertz	A and k	$\text{y}=(\text{A}+\text{u1})\text{exp}(-\text{Bexp}(-(\text{k}+\text{u}3)\text{t}))$
6	Bertalanffy	k	$\mathrm{y}=\mathrm{A}{(1-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
	Bertalanffy	B and k	$\mathrm{y}=\mathrm{A}{(1-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
7	Bertalanffy	A and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1){(1-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))}^3$
	Bertalanffy	A, B, and k	$\text{y}=(\text{A}+\text{u}1){(1-(\text{B}+\text{u}2)\text{exp}(-(\text{k}+\text{u}3)\text{t}))}^3$
8	Logistic	A and k	$\text{y}=(\text{A}+\text{u}1)/(\text{1}+\text{Bexp}(-(\text{k}+\text{u}3)\text{t}))$
	Logistic	A, B, and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)/(1+(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
9	Gompertz	k	$\mathrm{y}=\mathrm{Aexp}(-\mathrm{Bexp}(-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
	Gompertz	A, B, and k	$\mathrm{y}=(\mathrm{A}+\mathrm{u}1)\exp (-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$
10	Gompertz	B and k	$\mathrm{y}=\mathrm{Aexp}(-(\mathrm{B}+\mathrm{u}2)\exp (-(\mathrm{k}+\mathrm{u}3)\mathrm{t}))$

Table 6
Estimates of parameters and evaluation criteria of the models that best described the growth of naturalized chicken ecotypes

Graúna Dourada ecotype
	Female (y11)		Female (y2)		Male (y3)		Male (y4)
	Estimate	t-value	Estimate	t-value	Estimate	t-value	Estimate	t-value
A	1,301.44c	30.502	1,117.50b	40.562	1,230.70c	36.622	1,241.03c	42.562
B	4.59b	26.912	35.31a	12.792	55.38a	9.962	52.86a	10.962
k	0.0145b	21.042	0.0332a	36.812	0.0394a	28.182	0.0382a	28.632
MSE	5,004.96		1,797.25		3,058.32		2,923.01
AAD	50.38		30.53		39.72		39.04
R2	0.9719		0.9692		0.9936		0.9904
AIC	4,317.98		4,158.59		3,196.08		3,192.50
BIC	4,357.00		4,197.62		3,221.27		3,228.48
Nordestina ecotype
A	1,467.22b	18.852	1,159.05b	25.052	1,530.53b	26.182	1,517.73b	21.472
B	5.09a	12.482	24.31b	7.362	20.45c	6.342	24.85b	7.092
k	0.0203a	14.592	0.0332a	19.642	0.0278b	13.662	0.0306b	14.652
MSE	3,564.84		3,590.21		9,587.60		6,405.60
AAD	47.25		47.74		78.38		65.12
R2	0.9929		0.9310		0.9177		0.9438
AIC	1,417.45		1,420.36		1,438.13		1,431.72
BIC	1,445.08		1,447.98		1,457.35		1,459.17
Teresina ecotype
A	2,136.56a	7.302	1,353.36a	20.562	1,636.01a	28.252	1,635.92a	28.252
B	3.75c	18.852	21.30b	8.082	27.35b	7.382	27.36b	7.382
k	0.0107c	7.652	0.0281b	19.072	0.0293b	16.872	0.0293b	16.872
MSE	8,265.63		5,572.70		8,199.86		8,199.68
AAD	66.53		53.79		66.65		66.65
R2	0.9105		0.9155		0.9308		0.9308
AIC	1,416.79		1,407.34		1,758.97		1,764.97
BIC	1,444.24		1,434.79		1,779.67		1,794.53

• MSE - mean squared error; AAD - average absolute deviation of residuals; AIC - Akaike Information Criterion; BIC - Bayesian Information Criterion.
• Means followed by common letters between the estimates of A, B, and k parameters in the columns do not differ significantly from each other according to the Scott-Knott test (P<0.05).
• ¹
  Gompertz model: y1=(A+u1)exp(−(B+u2)exp(−(k+u3)t)); Logistic model: y2=(A+u1)/(1+(B+u2)exp(−(k+u3)t)); Logistic model: y3=(A+u1)/(1+Bexp(−(k+u3)t)); Logistic model: y4=(A+u1)/(1+(B+u2)exp(−(k+u3)t)).
• ²
  P-value <0.0001.

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.

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 above-mentioned ecotypes did not follow the behavior shown by the group in which they were clustered.

Figure 2
Observed and estimated weight (in grams) of naturalized female chickens, 2018.

Figure 3
Observed and estimated weight (in grams) of naturalized male chickens, 2018.

Figure 4
Observed and estimated weight (in grams) of naturalized chickens, 2018.

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.

Table 7
Pearson correlation coefficients between the A, B, and k parameters, above the diagonal, and respective P-values, below the diagonal

Parameter	Gompertz model - Males
	A	B	k
A	1	−0.9864	−0.9944
B	6.15 × 10−6	1	0.9982
k	4.46 × 10−7	1.4 × 10−8	1
	Logistic model - Males
A	1	−0.9262	−0.9717
B	0.0009	1	0.9890
k	5.55 × 10−5	3.28 × 10−6	1
	von Bertalanffy model - Males
A	1	−0.7439	−0.7340
B	0.0343	1	0.9997
k	0.0381	9.37 × 10−11	1
	Gompertz model - Females
A	1	−0.9348	−0.9861
B	0.0006	1	0.9768
k	6.63 × 10−6	3.06 × 10−5	1
	Logistic model - Females
A	1	−0.8071	−0.6426
B	0.0154	1	0.9348
k	0.0857	0.0007	1
	von Bertalanffy model - Females
A	1	−0.9817	−0.9933
B	1.51 × 10−5	1	0.9968
k	7.49 × 10−7	8.24 × 10−8	1

The Gompertz model estimated age at the inflection point between 38 and 46 days of life (Table 2). The Logistic and von Bertalanffy models estimated age at the inflection point between 96 and 130 days.

The widest variation of weight estimates at the inflection point were obtained by the von Bertalanffy model, with values ranging between 249 and 794 g.

Discussion

The significant difference in asymptotic weight between the sexes indicates that the growth of males and females should be evaluated separately. Sexual dimorphism is a marked characteristic in the species. This phenomenon can also be observed in the results presented by Topal and Bolukbasi (2008) and in the study led by Rizzi et al. (2013), who evaluated the growth of Italian chickens.

The best model fits were obtained by associating the random effects with the A and k parameters for data of males. The inclusion of random effects on these parameters generated similar values of the evaluation criteria when they were associated with the A, B, and k parameters—for the Logistic model, in both cases.

By arranging the AIC values of the females in descending order, it is observed that the models allocated in the first group, resulting from cluster analysis using Tocher's optimization method, showed the 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², 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.

The R² values showed dissimilarity to most criteria in the choice of the best model and were thus disregarded. They did not demonstrate efficiency as a criterion for model selection. A similar trend can be observed in the results published by Sarmento et al. (2006) and Teixeira et al. (2012).

In this way, the model that best described the growth of females of the Graúna Dourada and Teresina ecotypes was the Logistic model with random effects associated with the three parameters. The Gompertz model with random effects associated with the three parameters was the one that best described the growth of Nordestina females, based on the evaluation criteria. The Gompertz model with random effect on the A and k parameters is the most suitable for the study of growth of laying hens (Galeano-Vasco et al., 2014).

The growth of Grauná Dourada and Nordestina males was best described by the Logistic model with random effects associated with the three parameters. Although the Logistic model with random effect on the A and k parameters exhibited the lowest values of the AIC and BIC criteria for the Teresina males, the AAD and MSE values were similar to those of the model with three random effects: 66.65 and 8200, respectively.

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 non-linear 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.

Acknowledgments

The authors thank the Núcleo de Conservação de Galinhas Naturalizadas do Meio-Norte do Brasil (NUGAN-MN), financed by Banco do Nordeste, for the provided data.

References

• Aggrey, S. E. 2009. Logistic nonlinear mixed effects model for estimating growth parameters. Poultry Science 88:276-280. https://doi.org/10.3382/ps.2008-00317
  » https://doi.org/10.3382/ps.2008-00317
• Akaike, H. 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control 19:716-723. https://doi.org/10.1109/TAC.1974.1100705
  » https://doi.org/10.1109/TAC.1974.1100705
• Bertalanffy, L. von. 1957. Quantitative laws in metabolism and growth. The Quarterly Review of Biology 32:217-231. https://doi.org/10.1086/401873
  » https://doi.org/10.1086/401873
• Galeano-Vasco, L. F.; Cerón-Muñoz, M. F. and Narváez-Solarte, W. 2014. Ability of non-linear mixed models to predict growth in laying hens. Revista Brasileira de Zootecnia 43:573-578. https://doi.org/10.1590/S1516-35982014001100003
  » https://doi.org/10.1590/S1516-35982014001100003
• Guedes, M. H. P.; Muniz, J. A.; Perez, J. R. O.; Silva, F. F.; Aquino, L. H. and Santos, C. L. 2004. Estudo das curvas de crescimento de cordeiros das raças Santa Inês e Bergamácia considerando heterogeneidade de variâncias. Ciência e Agrotecnologia 28:381-388. https://doi.org/10.1590/S1413-70542004000200019
  » https://doi.org/10.1590/S1413-70542004000200019
• Hartley, H. O. 1961. The modified Gauss-Newton method for the fitting of non-linear regression functions by least squares. Technometrics 3:269-280. https://doi.org/10.1080/00401706.1961.10489945
  » https://doi.org/10.1080/00401706.1961.10489945
• Hruby, M.; Hamre, M. L. and Coon, C. N. 1994. Growth modelling as a tool for predicting amino acid requirements of broilers. Journal of Applied Poultry Research 3:403-415. https://doi.org/10.1093/japr/3.4.403
  » https://doi.org/10.1093/japr/3.4.403
• Karaman, E.; Narinc, D.; Firat, M. Z. and Aksoy, T. 2013. Nonlinear mixed effects modeling of growth in Japanese quail. Poultry Science 92:1942-1948. https://doi.org/10.3382/ps.2012-02896
  » https://doi.org/10.3382/ps.2012-02896
• Laird, A. K. 1965. Dynamics of relative growth. Growth 29:249-263
• Lindstrom, M. J. and Bates, D. M. 1990. Nonlinear mixed effects models for repeated measures data. Biometrics 46:673-687. https://doi.org/10.2307/2532087
  » https://doi.org/10.2307/2532087
• Mazucheli, J.; Souza, R. M. and Philippsen, A. S. 2011. Modelo de crescimento de Gompertz na presença de erros normais heterocedasticos: um estudo de caso. Revista Brasileira de Biometria 29:91-101
• Nelder, J. A. 1961. The fitting of a generalization of the logistic curve. Biometrics 17:89-110. https://doi.org/10.2307/2527498
  » https://doi.org/10.2307/2527498
• R Core Team. 2017. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria
• Rizzi, C.; Contiero, B. and Cassandro, M. 2013. Growth patterns of Italian local chicken populations. Poultry Science 92:2226-2235. https://doi.org/10.3382/ps.2012-02825
  » https://doi.org/10.3382/ps.2012-02825
• Romero, L. F.; Zuidhof, M. J.; Renema, R. A.; Robinson, F. E. and Naeima, A. 2009. Nonlinear mixed models to study metabolizable energy utilization in broiler breeder hens. Poultry Science 88:1310-1320. https://doi.org/10.3382/ps.2008-00102
  » https://doi.org/10.3382/ps.2008-00102
• Sarmento, J. L. R.; Rezazzi, A. J.; Souza, W. H.; Torres, R. A.; Breda, F. C. and Menezes, G. R. O. 2006. Estudo da curva de crescimento de ovinos Santa Inês. Revista Brasileira de Zootecnia 35:435-442. https://doi.org/10.1590/S1516-35982006000200014
  » https://doi.org/10.1590/S1516-35982006000200014
• Schwarz, G. 1978. Estimating the dimensional of a model. The Annals of Statistics 6:461-464
• Selvaggi, M.; Laudadio, V.; Dario, C. and Tufarelli, V. 2015. Modelling growth curves in a nondescript italian chicken breed: an opportunity to improve genetic and feeding strategies. The Journal of Poultry Science 52:288-294. https://doi.org/10.2141/jpsa.0150048
  » https://doi.org/10.2141/jpsa.0150048
• Sofaer, H. R.; Chapman, P. L.; Sillett, S. T. and Ghalambor, C. K. 2013. Advantages of nonlinear mixed models for fitting avian growth curves. Journal of Avian Biology 44:469-478. https://doi.org/10.1111/j.1600-048X.2013.05719.x
  » https://doi.org/10.1111/j.1600-048X.2013.05719.x
• Strathe, A. B.; Lemme, A.; Htoo, J. K. and Kebreab, E. 2011. Estimating digestible methionine requirements for laying hens using multivariate nonlinear mixed effect models. Poultry Science 90:1496-1507. https://doi.org/10.3382/ps.2011-01345
  » https://doi.org/10.3382/ps.2011-01345
• Topal, M. and Bolukbasi, S. C. 2008. Comparison of nonlinear growth curve models in broiler chickens. Journal of Applied Animal Research 34:149-152. https://doi.org/10.1080/09712119.2008.9706960
  » https://doi.org/10.1080/09712119.2008.9706960
• Teixeira, M. C.; Villarroe, A. B.; Pereira, E. S.; Oliveira, S. M. P.; Albuquerque, I. A. and Mizubuti, I. Y. 2012. Curva de crescimento de cordeiros oriundos de três sistemas de produção na Região Nordeste do Brasil. Semina: Ciências Agrárias 33:2011-2018. https://doi.org/10.5433/1679-0359.2012v33n5p2011
  » https://doi.org/10.5433/1679-0359.2012v33n5p2011

Publication Dates

• Publication in this collection
  15 May 2020
• Date of issue
  2020

History

• Received
  26 Sept 2019
• Accepted
  12 Mar 2020