Nonlinear mixed effects models for comparing growth curves for Guzerá cattle Modelos

The objective of the present work was to evaluate the accuracy of the fitted Gompertz and von Bertalanffy models for male and female Guzerá cattle, respectively. Four production regions in Northeast Brazil were included in the models as a fixed effect, and the animals were included as a random effect. In addition, the coefficients of the growth models in the production regions were compared. The accuracy of the fit equations was assessed with the Akaike information criterion, Bayesian information criterion, mean absolute deviation, mean squared error, and coefficient of determination. Confidence intervals were used for comparing the production regions. The Guzerá males in the Gado–Algodão and Serra Geral da Bahia production regions were statistically equal in asymptotic weight, and the animals in the Itapetinga-Valadares and Mata-Agreste regions had equivalent maturity rates. The Guzerá females in the Itapetinga-Valadares and Serra Geral da Bahia regions had the same asymptotic weight. The maturity rates in Itapetinga-Valadares were equal to those estimated for Mata-Agreste and Serra Geral da Bahia. The inclusion of the fixed effect of the production region and the random effect of the animals in the models improved the fit quality and increased the possibility of generating growth curves for each region. Bahia possuem peso assintótico estatisticamente igual, enquanto que a taxa de maturidade é equivalente para animais das regiões Itapetinga-Valadares e Mata-Agreste. As fêmeas da raça Guzerá das regiões de Itapetinga-Valadares e Serra Geral da Bahia possuem o mesmo peso assintótico. A taxa de maturidade em Itapetinga-Valadares é igual a estimada para Mata-Agreste e Serra Geral da Bahia. Ao incluir o efeito fixo de região de produção e aleatório de animal nos modelos, houve melhora na qualidade de ajuste e a possibilidade de gerar curvas de crescimento para cada região.


INTRODUCTION
Since its introduction in the late 19th century, the Guzerá breed has been shown to have the ability to adapt to the various climatic conditions in Brazil, especially in the Northeast Region (Santos et al., 2012). Because of the region's rusticity and meat and milk production potential, the Guzerá is a popular breed of Zebu cattle and an important resource for the cattle breeders in the tropical regions (Peixoto et al., 2014). Yield and meat quality are directly dependent on herd growth. Thus, animal growth must be evaluated for the purpose of managing the economic profitability of livestock (Souza et al., 2011). The criteria for selecting the best strains for use in genetic improvement programs have been informed by previous research on beef cattle growth (Bonin et al., 2014). In most species, body growth can be assessed by sigmoid curves and fitted to non-linear models . These models are suitable for the descriptions of biological processes because they better fit the functional relationship between specimen age and weight (Gonçalves et al., 2011). Nonlinear mixed-effects models contain fixed-and random-effects factors (Glória, 2014). They are effective for longitudinal data and can be applied to unbalanced or incomplete data (Pereira, 2014). These models are also used for estimating the average behavior and variability between individuals in a population (Pinheiro & Bates, 2000). In northeastern Brazil, there are several ecosystems with different climate and vegetation types and animal densities (Arruda & Sugai, 1994). With such environmental diversity, utilize a single equation to model cattle growth in the Northeast Region might not be appropriate. Thus, the fit of the models and growth curves for each production region or group of regions can enable the identification of growth patterns that more accurately represent the reality of each region. In addition, it increases the efficiency of the strategies for the genetic improvement, nutrition, and management of the cattle herds . The objective of the present work was therefore to evaluate the accuracy of the fitted Gompertz and von Bertalanffy models for the Guzerá males and females, respectively. The fixed effect of four Northeastern Brazil production regions and the random effect of the animals were added to the models. The estimated coefficients of the models for the production regions were compared.

MATERIALS AND METHODS
The data, which were provided by the Brazilian Association of Zebu Breeders, refer to the observations of age and body weight for 3,862 female and 2,899 male Guzerá cattle (Table 1). Each animal was weighed 6 to 9 times over time of 345 days. The cattle were raised in the Gado-Algodão, Mata-Agreste, Serra Geral da Bahia, and Itapetinga-Valadares regions, which have the highest Guzerá cattle densities (Arruda & Sugai, 1994). The Northeast Region, which accounts for a large part of Brazil and has greater diversity of soil, climate, and vegetation.
The Brazilian Northeast has contrasts between production regions, e.g., southeastern Bahia which is noted for its high production and efficiency. The other regions are characterized by family farming systems and low productivity. Of the 44 Brazilian cattle production regions identified and described by Arruda and Sugai (1994), 11 are located in the Northeast. To classify the production regions, they used the components of the production systems and edaphoclimatic conditions of the micro-regions, as defined by the Brazilian Institute of Geography and Statistics. Gado-Algodão, the name of this production region, comes from the integrated cattle-cotton production system. This region has an average annual rainfall of approximately 750 mm. It is generally hot and semi-arid with a seven-month dry season from June to December. Despite the severe climatic restrictions for agriculture and livestock activities, this region has good fertile soil. The herd yields are low, especially because of the long sevenmonth annual drought. The Mata-Agreste production region has mixed-purpose herds. Some areas are used mostly for dairy cattle to supply milk and dairy products to the states of Rio Grande do Norte, Paraíba, Pernambuco, Alagoas, and Sergipe. The rainfall in this region is higher in the east. It ranges from 750 to 2,000 mm per year. The rainiest period is from April to July. The average annual temperature is approximately 24 °C, with a low of 12 °C and a high of 36 °C. The region is therefore hot, humid, and semi-arid with a dry period that can last three to five months. Except for the sugarcane areas in Pernambuco and Alagoas, this region has one of the highest bovine densities in Northeastern Brazil. The Serra Geral da Bahia, which accounts for most of the "Atlantic plateau" in Bahia, is the most mountainous area in Bahia. It has an annual average maximum temperature of 24 °C and annual rainfall of 600 to 1,000 mm. The dry season in the winter lasts five to six months in the south and three to five months in the north. Two areas are noteworthy for their annual rainfall levels: Vitória da Conquista with 600 mm and Rio de Contas with 1,000 mm. Because of the limitation of natural resources, notably soil and climate, cattle raising is limited to the pre and postweaning phases, associated with cultivation of manioc, beans, corn, cotton and castor beans crops. The Itapetinga-Valadares production region is considered to be one of the main cattle backgrounding areas in Brazil. It is in the Southeast Region (Arruda & Sugai, 1994); however, it covers a large area of the south and southeast of Bahia. It extends from the Doce River valley in Espírito Santo and northeast of Minas Gerais to the Contas River valley in the southeast of Bahia. The average rainfall in this region ranges from 1,250 to 1,500 mm per year. The greater concentrations of rainfall occur in November, December, and January in the south region. This gradually changes across the region such that in the north, most of the rainfall is in February, March, and April. These climatic variations, including the short dry season, which is shorter than three months, make this region a regular supplier of cattle with the optimal weight for slaughter. This cattle-raising region has the highest yields of the Northeast Brazil and, consequently, an advantage in meat production. The models used for the males and females were based partly on the work of Alves (2016). The Gompertz and von Bertalanffy models were selected. Both models were fitted while considering three fit classes: traditional non-linear model (M1 and M4), non-linear model with a fixed effect of region in the all parameters (M2 and M5) and non-linear mixed effects model with a fixed effect of production region on all the parameters and a random effect of animals in β1 and β3 (M3 and M6). The growth curves for the males were generated from the Gompertz model: . The von Berttalanffy model: , was used for plotting the growth curves for the females. In this model, β 1 is the asymptotic weight that is interpreted as the weight of an adult animal, β 2 is the constant of integration without biological interpretation, and β 3 is the maturity rate. The estimation method of the models was the maximum likelihood with the algorithm created by Lindstrom and Bates (1990) was used to approximate the integral. The models were fitted in package nlme of the R software, in which the default assumes the homoscedasticity of the errors (R Core Team, 2019). To assess the accuracy of the fitted models, the criteria used by Silva et al. (2011) and Reis et al. (2014)  The confidence intervals of the parameters of the models fitted for each production region were used for testing the hypothesis of the equality of the parameters (asymptotic weight or maturity rate) for the two production regions. Thus, as was the case with the criterion used by Pereira (2016), the overlapping confidence intervals indicated that the related models had parameters with common estimates, i.e., statistically equal parameters. The 95% confidence intervals were constructed through the use of the Gompertz and von Bertalanffy models for the males and females, respectively.

RESULTS AND DISCUSSION
A comparison of the traditional nonlinear models (M1 and M4) and the models with a fixed effect of production region (M2 and M5) showed that there was a small improvement in the quality of the fit criteria for both sexes (Table 2). A comparison of the M1 and M4 models and the M2 and M5 models showed that the AIC and BIC for M2 and M5 decreased by approximately 0.5%. The MAD and MSE were reduced by approximately 4% and 6%, respectively, and R 2 increased by approximately 2%.  (Table 2). When the parameter estimates for M3 and M6, the models with the best fit, were taken into account, there were considerable differences in the estimates of asymptotic weight and maturity rates for both sexes in the production regions (Tables 3 and 4). For the males, the estimated asymptotic weight was 351.4 to 580.4 (Table 3). For the females, it was 296.6 to 436.0 (Table 4). The estimates obtained from the M1 and M4 models, which did not consider the effect of the production region, were 453.2 and 441.0 for the males and females, respectively. These results show that the inclusion of the fixed effect of region in the models allowed for a more suitable description of animal growth in the production regions. The highest and lowest estimated asymptotic weight for the males and females was found in the Mata-Agreste and Gado-Algodão regions, respectively. The regions with a high estimated asymptotic weight had a low maturity rate. This negative association has been previously described (Oliveira et al., 2000;Garnero et al., 2005;Malhado et al., 2009).   Oliveira et al. (2000) used the von Bertalanffy model to analyze the fit of the growth curves for the Guzerá females. They estimated a curve for each animal and presented an average curve with an asymptotic weight equal to 453.2, a value similar to that (441.0) estimated in the M4 model, and a maturity rate of 0.065, which is higher than that (0.0024) found for the M4 model. However, the Gompertz and von Bertalanffy models yielded lower estimates of asymptotic weight, 337.6 and 357.8, respectively, and maturity rate, 0.0050 and 0.0030, respectively, for the data from the Nelore breed in the north of Brazil (Lopes et al., 2011). In a study of cattle raised in the State of Pernambuco, Santoro et al. (2005) estimated the asymptotic weight of the females (352.2) and males (463.2) of the Guzerá breed by using the von Bertalanffy and Gompertz models, respectively, without considering the effects of the production regions. The State of Pernambuco extends over parts of the Gado-Algodão and Mata-Agreste production regions. This is indicative of the great diversity in the edaphoclimatic conditions and livestock production systems.
The results of the present study indicate that the Gado-Algodão and Mata-Agreste regions have the greatest asymptotic weight discrepancies (Tables  3 and 4). Therefore, the use of traditional models to fit a single growth curve for animals from an entire state, such as Pernambuco, without consideration of the characteristics of each production region is not appropriate.
With the consideration of the overlap of the 95% confidence intervals for comparing the growth curves for males in the production regions, the males in the Gado-Algodão and Serra Geral da Bahia regions were found to have the same asymptotic weight, as well as to have the lowest estimates. The highest maturity rate was found in the Gado-Algodão region. In the case of the Mata-Agreste region, which had the greater asymptotic weight, the estimated maturity rate was equivalent to that in the Itapetinga-Valadares region where the lowest maturity rate was found (Table 3). The comparison of the models that were fitted for the females showed that the Gado-Algodão region had the lowest asymptotic weight and the highest maturity rate (Table 4). In contrast, the Mata-Agreste region had the highest asymptotic weight and the lowest maturity rate. The maturity rate of the Mata-Agreste females was statistically equivalent to that of the Itapetinga-Valadares females, whose maturity rate was equal to that of the Serra Geral da Bahia females. The regional variations in the growth parameters are a reflection of the effects of the climatic and other non-genetic factors on the physiology of the cattle (Silva et al., 2001(Silva et al., , 2004. Therefore, the superior results in the Mata-Agreste region could have been caused by the favorable climatic conditions, which were evident in the comparisons of the environmental characteristics of Gado-Algodão. Silva et al. (2001) attributed the differences in the male and female growth rates to hormonal factors and to the more extensive breeding system used for males.

CONCLUSION
The inclusion of the production region as a fixed effect and the animals as a random effect in the model improved the fit quality. It allowed for fitting a growth curve and estimating the specific confidence intervals for each production region. The Guzerá males in the Gado-Algodão and Serra Geral da Bahia production regions had statistically equivalent asymptotic weight. The maturity rates in the Itapetinga-Valadares and Mata-Agreste regions were equivalent. The highest and lowest asymptotic weight was found in the Mata-Agreste and Gado-Algodão regions, respectively. The males in Gado-Algodão had the highest maturity rate, and the Guzerá males in Itapetinga-Valadares had the lowest maturity rate. The Guzerá females in the Itapetinga-Valadares and Serra Geral da Bahia regions had the same asymptotic weight. The maturity rate in Itapetinga-Valadares was equal to the rates estimated for Mata-Agreste and Serra Geral da Bahia. Unlike the growth pattern of the males, the lowest maturity rate for the females was found in Mata-Agreste.

ACKNOWLEDGMENTS
The co-authors thank the Brazilian Association of Zebu Breeders for making the data available.