Livestock Research for Rural Development 25 (8) 2013 | Guide for preparation of papers | LRRD Newsletter | Citation of this paper |
The objective of this study was to analyze buffalo growth based on body weight, Longissimus dorsi muscle area (AOL), and fat deposition over the hip (FOH) using random regression mixed models of first (FORRM) and second order (SORRM), each with nine different variance structures. Ten measurements for each trait were taken on 26 animals during the first performance test (93 days test plus 23 d adaptation period) developed for buffaloes in Colombia. Computations were performed using the lme procedure of the R program. Preliminary analyses determined that an SORRM was appropriate for body weight (BW) and FOH and an FORRM was suitable for AOL. The maximum likelihood ratio (MLR), the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) were used to compare models.
The best models were an SORRM with homogeneous residual variances for BW, an FORRM with heterogeneous animal residual variances for AOL, and an SORRM with heterogeneous residual variances among farms times an exponential function of age for FOH. Heterogeneity of residual variances was likely due to environmental differences among farms, and to genetic differences among buffaloes not accounted for by FORRM and SORRM. Fixed intercepts with the best models for each trait were 227 ± 7.90 kg for BW, 34.8 ± 0.99 cm^{2} for AOL, and 4.19 ± 0.229 mm for FOH. Fixe linear regression coefficients were 1.29 ± 0.073 g / d for BW, 0.0584 ± 0.0042 cm^{2 }/ d for AOL, and 0.0035 ± 0.0032 mm / d for FOH. The fixed quadratic regression coefficient indicated that BW rate decreased after one year of age whereas FOH rate continued to increase until the end of the test. Random regression coefficients suggested that there was considerable variability among trait curves for individual buffaloes, particularly for FOH. The evaluated residual variance structures did not eliminate completely the heterogeneity of variances.
Keywords: body weight, meet quality, performance test, ultrasound
Experiments with repeated measures are common in livestock research (ZooBell et al 2003; Wang and Goonewardene 2004). In most cases, multiple observations per experimental unit are taken over time. The usual assumptions of independence and homogeneity of variances are usually not valid for data analysis in such experiments because the measurements made in the same animal are often correlated with each other and the variances between measurements may be different. Mixed model methodology allows for correct and efficient data analysis of experiments with repeated measures through modeling of the covariance structure, considering the correlations between repeated measures and the presence of heterogeneous variances (Littell et al 1998; Wang et al 2006).
Growth is one of the factors of greatest economic importance in beef production systems. Mixed models can be used to describe growth of farm animal species. Incorporating random effects in the models permits accounting for variability among growth curves from individual animals within a population (France et al 1996; El Halimi 2005). This is useful in breeding programs because it helps to choose the fastest growing animals without altering their adult weight (Malhado 2008). In addition, it allows to classify the productivity of a breed, a production system or a given population, and to measure genetic changes from one generation to another (Agudelo et al 2008).
The objective of this study was to analyze buffalo growth based on body weights, Longissimus dorsi muscle areas, and fat thicknesses over the hip using random regression mixed models with different residual variance structures using data from the first performance test developed in Colombia.
The performance test was conducted at the Experimental Center of the University of Antioquia, located in the municipality of Barbosa, Antioquia, Colombia. This area is classified as subtropical humid forest. It has an average altitude of 1300 m above sea level, an average annual temperature of 23 °C, and an average annual precipitation of 1800 mm per year.
A total of 26 buffaloes were used in the performance test. Animals came from four breeding farms where milking was not practiced. Farms were located in the municipalities of Buenavista and Ayapel (Córdoba province), Cimitarra (Santander) and Norcasia (Caldas). The performance test was conducted between July and October 2009 for a period of 116 days. There was an adjustment period of 23 d, and the evaluation period lasted 93 d. Animals entered the test at an average age of 285 ± 25.72 d with an average weight of 261 ± 29 kg, and they reached an average weight of 342 ± 39.08 kg at the end of the test.
Animals were confined in 16 m^{2} individual pens, with cement floor, without bedding, and with a covered area of 4 m^{2} where feeders and drinkers were located. The diet consisted of fresh cut Maralfalfa grass (Pennisetum sp.) plus two kilograms of a feed supplement per day. The supplement ingredients were corn (50%), corn protein (10%), extruded soybeans (15%), soybean meal (10%), extruded corn (10%), mineral salt containing 8% phosphorus (2% ), calcium carbonate (1%), and 2% of a vitamin and mineral premix. The pasture had 9% CP, 4.0 Mcal / kg GE, 37% NDF, and 75% ADF. The supplement had 17% CP, 2.84 Mcal / kg GE, and 75% TDN.
A total of 10 measurements per trait were taken over the 116 d of the performance test, 3 measurements during the adjustment period and 7 measurements during the test period. Animals were weighed (BW) after 12 h of fasting. Area of the Longissimus dorsi (AOL) and fat thickness over the hip (FOH) were measured by ultrasound, using an Akila-Pro equipment (Esaote Europe BV, Maastricht, The Netherlands) with a 3.5-MHz, 18cm transducer. Images were measured with the Eview program (Pie Medical, Maastricht, The Netherlands; Pie Medical Equipment B. V. 1996). Images were collected for AOL between the 12^{th} and the 13^{th} rib and for FOH at the junction of the biceps femoris and gluteus medius between the ischium and illium and parallel to the vertebral column (Jorge et al 2005; Urdapilleta et al 2005).
Data were first analyzed using additive random regression mixed models of first order (FORRM; Model 1) and of second order (SORRM; Model 2), both with homogeneous residual variance structures. These 2 base models were as follows:
Where: γ_{ij}: is a BW, AOL and FOH record measured in the i-th animal at the j-th age; β_{0}, β_{1 } y β_{2} fixed intercept, linear regression coefficient, and quadratic regression Coefficient averaged over all animals; χ_{j}: is the j-th age; ε_{ij}: is the random residual due to effects not considered in the model, where ε_{ij } ~ Ν (0, σ^{2}), and β_{0i}, β_{1i } y β_{2} random intercept and linear regression coefficient for the i-th animal, representing random deviations from the β_{0 } y β_{1} fixed regression coefficients, respectively, where:
Animals were assumed to be unrelated, thus covariances between random intercepts and random linear regression coefficients among different animals were assumed to be zero.
Linear mixed models of first, second and third order with one (β_{0}), two (β_{0}, β_{1}) and three (β_{0}, β_{1 } y β_{2}) random parameters were compared using Likelihood Ratio Tests in preliminary analyses. It was determined that an SORRM was appropriate for BW and FOH and an FORRM with two random parameters (β_{0}, β_{1}) was appropriate for AOL.
Subsequently, 9 heterogeneous residual variance structures (Table 1) were evaluated within models 1 and 2. Computations were carried out with the lme procedure of the nlme library of the R program (R Development Core Team 2008). The nlme library variance functions were used to model the 9 heterogeneous residual variance structures (Pinheiro and Bates 2000). The variance covariate used for all structures was the discrete variable age of animal. Structures 3 and 4 took into account the variance among individual animals and among farms, respectively. The R function used for these 2 residual structures was varIdent as suggested by Zuur et al (2009) for discrete variance covariates. For structure 5, the residual variance was modeled as times the absolute value of age of animal raised to the power of 2δ using R function varPower, where δ is unknown and must be estimated. Parameter δ has no restrictions, thus it can take any value (Pinheiro and Bates 2000; Zuur et al 2009). For structure 6, residual variances were calculated as multiplied by an exponential function of 2δ times the age of animal using R function varExp. Notice that this structure allows age to be zero, thus for cases where the variance covariate may be zero this structure is a good choice (Zuur et al 2009).
Residual variance structure 7 is a constant plus a variance covariate power function. It was computed with R function varConsPower. According to (Pinheiro and Bates 2000), structure 7 works better than the exponential when the variance covariate takes values close to zero. Structures 8, 9, 10 and 11 were combination of previous residual covariance structures. Structure 8 allows residual variance structure 6 to vary among individual animals, whereas structure 9 permits residual variance structure 7 to vary among individual animals. Similarly, structure 10 allows residual variance structure 6 to vary among farms, while structure 11 let residual variance structure 7 vary among farms. Variances were estimated using restricted maximum likelihood procedures (Harville 1977; Searle et al 1992) and computed by the lme procedure of the R program.
Because heteroscedasticity may not have a single origin, residual diagnostic plots were generated according to various factors (farm, animal, age, and measurement). Considering that farm, animal, age, and measurement are discrete variables, boxplot-type graphs were used.
Animals in the performance test came from farms that had different climatic conditions and used different production systems (direct grazing of different forage species). Thus, it was necessary for them to get adjusted to the conditions of the performance test (climate, confinement, fed cut grass instead of direct grazing). To evaluate the importance of the adjustment period, the 9 models were run first with all 10 measurements for the 3 traits, and then with only the 7 measurements of the evaluation period (i.e., without the first 3 measurements of the adaptation period).
Although residuals per animal decreased significantly, there were still heterogeneous variances, possibly a consequence of carry over environmental effects from their farms of origin not eliminated by the adjustment period, and perhaps genetic differences among animals not accounted for by the additive random regression models used here. Buffaloes entered the test with a maximum age difference of 73 d and an initial weight difference of 91 kg (215 to 306 kg). This maximum weight difference increased to 168 kg at the end of the test (244 to 412 kg). In addition, individual differences were observed in growth curves. For example, animals 14 and 21 grew slowly at the beginning and then the gain increased. On the other hand, animal 12 grew faster at the beginning and slowed down at the end, while animal 6, had a similar growth rate throughout the evaluation period. These differences in growth patterns may help explain the higher residual variances obtained from these individuals (Figure 1).
Estimates for fixed regression genetic effects were similar for all models. The estimate for , associated with initial weight, fluctuated between 208 ± 7.28 and 214 ± 6.76 kg when all measurements were considered, and between 225 ± 7.86 and 228 ± 7.90 kg without measurements from the adjustment period. The estimate for , associated with growth rate, ranged from 0.975 ± 0.048 to 1.09 ± 0.056 kg / d with all measurements considered, and between 1.289 ± 0.073 and 1.331 ± 0.066 kg / d without the measurements of the adjustment period. Lastly, the estimate for ranged from -0001 ± 0.0002 to -0.0016 ± 0.00021 kg / d for all measurements, and between -0.003 ± 0.0003 and -0.0033 ± 0.0002 kg / d without the measurements of the adjustment period. By eliminating the measurements of the adjustment period, became larger, indicating
A linear mixed model of second order was used by Ramírez et al (2011), to describe growth of 123 buffaloes castrated males from 12 months to slaughter. These authors estimated fixed effects lower than those found in this study, with values of 184 kg, 0.948 kg / d -0.00061 kg / d for β_{0}, β_{1 } y β_{2}, respectively. These results can be explained by differences in physiological status (castrated vs. non-castrated buffaloes) and the feeding system (grazing vs. confinement).
The average weight at the start of the test (261 kg ± 29) was similar to those reported by Lourenço et al (2010) in Eastern Amazon Embrapa (Belém, Pará, Brazil) for Murrah buffaloes in 2 weight gain performance tests (265 and 273 kg). This suggests that the genetic potential of the buffaloes used in the performance test here was comparable to that of buffaloes in the 2 tests in Brazil, a country with one of the largest buffalo populations with a well-established genetic improvement program. Malhado et al (2008), evaluated the weight at different ages of Mediterranean buffaloes using the database of buffalo breeding program in Brazil (PROMEBUL), obtaining an adjusted weight at 365 days of age similar to that obtained in this research (302 kg), confirming the genetic potential of buffaloes used in this test.
The daily gains obtained by Lourenço et al (2010) were 911 and 969 g in the first and second trials, during the 224 day trial (after the adjustment period), which were lower than that obtained here. The daily gains obtained during the adaptation period (70 days) by Lourenço et al (2010) were very low (257 and 285 g), indicating the importance of taking this period into account.
Crudeli et al (2007) reported 282 kg weaning weight at seven months for 170 Mediterranean buffaloes in Argentina, which is higher than the weight at the start of performance test here. However, these authors also reported a 335 kg average weight at 18 months of age. This is lower than that the average weight of 343 kg obtained in the performance test at 14 mo of age here. Lower weights were also reported by Jorge et al (2005) for Murrah buffaloes reared on pasture (251 kg at 12 months), by Nogueira et al (1989) for buffaloes reared at pasture in the northwest of São Paulo, Brazil (301 and 310 kg at 12 months for Mediterranean and Jafarabadi buffaloes, respectively) and Barbosa et al (1988) for Mediterranean buffaloes (299 kg al 12 months).
Random effects yielded similar values for the standard deviation of b_{0} and b_{1}; and for the covariance between b_{0} and b_{1} for all models, with a higher degree of similarity when measurements for the adjustment period were discarded (Table 2). However, the standard deviation of residuals varied widely across models when all measurements were considered (from 2.72 kg for model 3 to 2.45e^{+08} kg for model 9), whereas they differed minimally among models when measurements from the adjustment period were eliminated (Table 2).
Ramírez et al (2011), reported in a group of buffalos tested from 12 months of age to slaughter higher values for standard deviation of b_{0} and b_{1}_{ } (41.6 and 2.87, respectively), indicating less variability between animals in the performance test. By contrast, the standard deviation for residuals reported by Ramírez et al (2011) were higher, which can be explained by the greater control of environmental conditions on the performance test and the use of models that take into account heterogeneous residual variances.
As for BW, analyses for AOL were performed taking into account all measurements and eliminating those from the adjustment period. Similar impact of measurements from the adjustment period on residuals were observed for AOL, thus only results from analyses without measurements from the adjustment period are presented here.
Box-plot graphs showed an increase in standardized residual values when plotted against fitted values. Similarly, residual values increased with age of the animal. Box-plots also revealed heterogeneity of residual variability among farms, measurements, and animals. Thus, unlike BW, eliminating measurements from the adjustment period was insufficient to reduce residuals and improve the homogeneity of variances. When comparing models with heterogeneous residual variance structure, significant differences (Table 3) were found with the LRT between models with structures 3 (P<0.0288), 8 (P<0.0312) and 9 (P<0.0309) and model 1 (homogeneous variances). Model 3, which assumed a different residual variance for each animal showed the lowest AIC and BIC values (Table 3), thus it was considered to be the best model to analyze AOL. There was no clear evidence of heterogeneity with this model (Figure 2, Panel A) and there was higher uniformity by age (Figure 2, Panel B). However, although there was a reduction of the size of the residuals by farm, measurements, and animal (Figure 2, Panels C, D and E), heterogeneity of residual variances still remained.
As previously indicated, animals were exposed to different climatic, feeding, and management conditions in their farms of origin. Buffaloes in each of the 4 farms were also likely to be genetically different subpopulations. Thus, although animals entered at a similar age, initial AOL ranged from 25.1 to 45.2 cm^{2}, a difference that was maintained until the end of the test, with values between 35.6 and 53.9 cm^{2}. Similarly, there were differences among animals in tissue deposition, which may help explain the heterogeneous residual variances. For example, animals 3, 6 and 7 showed no muscle deposition at the start of the test, and then they exhibited a high rate of growth, whereas animals 2 and 8 had high growth rate at the beginning of the test, and subsequently their growth rate became exceedingly slow. The highest residual variances were observed for these 5 animals (Figure 2).
Comparable fixed regression genetic parameters were observed among models. The estimate for β_{0}, associated with the initial AOL, was between 34.5 ± 1.01 and 34.9 ± 0.99 cm^{2}, and the estimate for β_{1} associated with the growth rate, was between 0.0572 ± 0.004 and 0.0620 ± 0.0042 cm^{2}/d. Ramírez et al (2011), reported 20.3 cm^{2} at 12 months for buffalo males, which is lower than the AOL at the start of performance test here, but the growth rate (β_{1}) was higher (0.124 cm^{2 } / d) than that obtained in this work. The higher growth rate may have been due to a longer evaluation period (from 12 months to slaughter), allowing for more muscle deposition produced at an older age (Berg and Butterfield 1976).
Similar random regression genetic effects were also obtained across models. The standard deviation of b_{0} (σ b_{0}) and the standard deviation of b_{1} (σ b_{1}) related to individual animals ranged from 4.91 to 5.03 cm^{2}, and from 0.0159 to 0.0172 cm^{2 }/ d, respectively. On the other hand, the standard deviation for residual effects was substantially different across models, with values ranging from 0.001 to 4.80 cm^{2}. Ramírez et al (2011) reported a lower (σ b_{0}) (3.71), which indicates a greater variability between animals in the performance test and a higher standard deviation for residuals (4.92), which can be explained by the greater control of environmental conditions on the performance test and the use of models that take into account heterogeneous residual variances.
Inclusion of measurements taken during the adjustment period had similar effects on models as for BW and AOL. Thus, only FOH results obtained without measurements from the adjustment period are shown. Box-plots of residuals estimated for FOH with model 2 assuming homogeneity of variances suggested that there was a violation of the assumption of homogeneity of residual variances. Box-plot graphs showed an increase in standardized residual values when plotted against fitted values (Figure 3, Panel A). Similarly, residual values increased with age of the animal, coinciding with the greater residual observed in the last three measurements (Figure 3, Panels B, D). Heterogeneity of residual variance also existed among farms (Figure 3, Panel C). Farm 1 had the largest residual, perhaps because the body condition (closely associated with fat deposition) of buffaloes from this farm entering the performance test was highly variable. Residual variability among animals also differed widely (Figure 3, Panel E). Animal 23 had the highest variance and it also was the buffalo that deposited fat at the fastest rate.
Table 4 shows the estimates of regression parameters obtained by the models with different variance structures and comparison criteria. All models were significantly different from model 2 (homogeneous variances). The best model (lowest AIC and BIC values) was model 10. This model took into account the variance for each farm and an exponential variance for age. Figure 5 shows the graphs of residuals for FOH, estimated by model 10. Figures 4 and 5 (Panel A) show that heterogeneity of residual variances was substantially corrected. However, a clear heterogeneity of variance due to age (Panel B) was still observed, and this heterogeneity increased with age and measurement (Panel D). This increase in residual variance for FOH can be explained by an increase in fat deposition as the animal grows (Putrino et al, 2006). Heterogeneity of variance due to farm (Panel C) and animal (Panel E) was only partially corrected by model 10.
In agreement with BW and AOL, animals diverged widely in FOH at the beginning of the test, with values ranging from 2.53 to 6.40 mm. Although residuals per animal decreased significantly, there were still heterogeneous variances, possibly a consequence of carry over environmental effects from their farms of origin not eliminated by the adjustment period, and perhaps genetic differences among animals not accounted for by the additive random regression models used here. The growth of muscle and fat appears to vary widely between groups of animals (Shuey et al 1993; Crews and Kemp 2001; Silva et al 2002) while the growth of bone appears to be similar regardless of the type of cattle (Silva et al 2002).
Similar values for fixed and random regression coefficients were obtained for all models. Estimates for β_{0} ranged from 4.19 ± 0.23 to 4.48 ± 0.28 mm and for β_{2} fluctuate between 0.00005 ± 0.00002 and 0.00009 ± 0.00002 mm. Estimates for β_{1} were low and ranged from negative to positive values. Seven models yielded negative values indicating fat loss (models 2, 3, 4, 5, 7, 8, and 9), whereas 3 models (6, 10, and 11) had positive values indicating that animals had fat gains during the performance test. Estimates for β_{1} ranged from -0.0042 ± 0.0045 to 0.0035 ± 0.0033.
Values obtained for were higher than those reported by Ramírez et al (2011), who obtained a value of 0.65 mm for male buffalo at 12 months old. These authors reported a higher growth rate (0.042 mm / d) and negative (-0.0004 mm / d). These differences may be explained by the age of the animals tested, since the process of accumulation in animal tissue has certain rules of priority: first grow the viscera, then the bone and muscle tissue and finally the adipose tissue (Cavalho et al 2003). In this study, the animals entered an age of 285 days and ended with 401 days of age (13.4 months), while Ramírez et al (2011) assessed the animals from 12 months of age until slaughter.
The evaluated residual variance structures did not eliminate completely the heterogeneity of variances. Heterogeneity was likely due to environmental differences among farms where animals originated from, and to genetic differences among buffaloes not accounted for by the random regression additive genetic models used here. When BW was analyzed without the 3 weights taken during the adjustment period, residuals were substantially reduced and weight gain increased, showing the importance of allowing buffaloes sufficient time to adjust to the environment at the performance test station. There was no reduction of residuals for AOL and FOH when the adjustment period was eliminated. The best models to describe buffalo growth in this performance test were model 2 with homogeneous residual variances for BW, model 1 with different residual variances for each animal for AOL, and model 2 with heterogeneous residual variances among farms times an exponential function of age for FOH. Estimates of genetic fixed and random regression coefficients were similar across models, however residual variances varied considerably. Rate of BW change decreased after one year of age, while rates of muscle (AOL) and fat (FOH) deposition continued to increase until the end of the performance test. Growth curves for individual animals showed considerable variation, particularly for fat.
Project funded by the Colombian Ministry of Agriculture and Development, Universidad de Antioquia, Corporación Universitaria Lasallista, and Asociación Colombiana de Criadores de Búfalos. Agreement beetwen Universidad de Antioquia and Fundación Universitaria San Martín for training graduate studentsThe authors thank all the people who contributed to the project, especially undergraduate and MS students from the Universidad de Antioquia for their work at the Experimental Station, and CODI-Sostenibilidad 2013.
Agudelo G D A, Cerón-Muñoz M F y Restrepo L F B 2008 Modelación de las funciones de crecimiento aplicadas a la producción animal. Revista Colombiana de Ciencias Pecuarias, 21:39-58, from http://rccp.udea.edu.co/index.php/ojs/article/viewFile/330/327
Akaike H 1974 A new look at the statistical model identification. Transactions on. Automatic. Control 19 (6):716-723, from http://www.unt.edu/rss/class/Jon/MiscDocs/Akaike_1974.pdf
Barbosa C, Nogueira J R e Mattos J C A 1988 Desempenho de bubalinos da raça Mediterrâneo (leiteiros) para a produção de carne. Comunicações científicas da Faculdade de Medicina Veterinária e Zootecnia da Universidade de São Paulo. 12: 173-235.
Berg R T and Butterfield R M 1976 New concepts of cattle growth. Ed. Sidney Univ. Press. Australia, 240, from http://ecommons.library.cornell.edu/handle/1813/1008.
Carvalho P A, Bonnecarrère S L M, Pires C C, Viégas J, Velho J P e Paris W 2003 Composição corporal e exigências líquidas de proteína e energia para ganho de peso de bezerros machos de origem Leiteira do nascimento aos 110 dias de idade. Revista Brasileira de Zootecnia. 32 (6): 1484-1491, from http://www.scielo.br/pdf/rbz/v32n6/18438.pdf
Crews D H Jr and Kemp R A 2001 Genetic parameters for ultrasound and carcass measures of yield and quality among replacement and slaughter beef cattle. Journal of Animal Science 79:3008–3020, from http://www.journalofanimalscience.org/content/79/12/3008.full.pdf
Crudeli G, Pochon D, Olazarri M, Monzón N, Chaparro L, Flores S, Patiño E and Cedrés J 2007 Morphometric evaluation of male Mediterranean buffaloes in Northern Corrientes, Argentina. Italian Journal of Animal Science 6 (Suppl. 2): 1281-1283, from http://www.aspajournal.it/index.php/ijas/article/view/ijas.2007.s2.1281
El
Halimi R 2005 Nonlinear mixed-effects models and nonparametric inference. A Method
Based on Bootstrap for the Analysis of Non-normal Repeated Measures Data in
Biostatistical Practice. PhD Dissertation, University of Barcelona, Spain,
from
http://www.tdx.cat/bitstream/handle/10803/1556/TESIS_RACHID_EL_HALIMI.pdf?sequence=1
France J, Dijkstra J and Dhanoa M 1996 Growth functions and their application in animal sciences. Annales de Zootechnie. 45: 165-174, from http://animres.edpsciences.org/index.php?option=com_article&access=doi&doi=10.1051/animres:19960637&Itemid=129
Harville D A 1977 Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association. 72 (358):320-338, from http:// http://www.lce.esalq.usp.br/arquivos/aulas/2012/LCE5872/Harville1977.pdf
Jorge A M, Andrighetto C, Domingues D M, Golfetto C M e Freitas V A D 2005 Características Quantitativas da Carcaça de Bubalinos de Três Grupos Genéticos Terminados em Confinamento e Abatidos em Diferentes Estádios de Maturidade. Revista Brasileira de Zootecnia. 34 (6): 2376-2381, from http://www.scielo.br/pdf/rbz/v34n6s0/a25v3460.pdf
Littell R C, Henry P R and Ammerman C B 1998 Statistical analysis of repeated measures data using SAS procedures. Journal of Animal Science 76 (4): 1216-1231, from http://journalofanimalscience.org/content/76/4/1216.full.pdf+html
Lourenço J J de B, da Costa N M, Araújo C V, Dutra S, Rossetto G A, Nahúm B de S, Silveira de M J C and Brandão L de M 2010 Sistema silvipastoril na produção sustentável de búfalos para carne na pequena propriedade da amazônia oriental. Bioclimatologia Animal, Quinta-Feira 3 , from. http://www.bioclimatologia.ufc.br/lourenco.pdf
Malhado C H M, Ramos A A, Carneiro P L S, Souza J C, Wechsler F S, Eler J P, Azevêdo D M M R y Sereno J R B 2008 Modelos no lineales para describir el crecimiento de bufalinos de la raza Murrah. Archivos de Zootecnia. 57 (220): 497-503, from http://www.uco.es/organiza/servicios/publica/az/php/img/web/26_18_40_11ModelosMalhado.pdf
Nogueira J R, Barbosa C e Mattos J C A 1989 Peso ao nascer e desenvolvimento ponderal de bubalinos das raças Mediterrâneo e Jafarabadi. Boletim de Indústria Animal, Nova Odessa, SP. 46 (2): 193-198.
Pie Medical Equipment B V 1996 Eview – Echo Image Viewer. Version 1.0. Maastricht, The Netherlands.
Pinheiro J C and Bates D M 2000 Mixed-Effects Models in S and S-PLUS. Springer, NY. 523, from http://books.google.com.co/books/about/Mixed_Effects_Models_in_S_and_S_Plus.html?id=N3WeyHFbHLQC&redir_esc=y
Putrino S M, Leme P R, Silva S da L, Alleoni G F, Lanna D P D e Grossklaus C 2006 Exigências líquidas de proteína e energia para ganho de peso de novilhos Nelore alimentados com dietas contendo grão de milho úmido e gordura protegida. Revista Brasileira de Zootecnia. 35 (1): 301-308, from http://www.scielo.br/pdf/rbz/v35n1/28373.pdf
Ramírez E J, Mesa J A, Agudelo D A, Bolívar D M and Cerón-Muñoz M F 2011 Using mixed models to describe growth in buffaloes. Livestock Research for Rural Development. 23 (9), from http://www.lrrd.org/lrrd23/9/rami23189.htm
R Development Core Team 2008 R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, from http://spatial- analyst.net/book/node/330
Searle S R, Casella G and McCulloch C E 1992 Variance Components. John Wiley & Sons, New York.
Schwarz G 1978 Estimating the dimension of a model. The Annals of Statistics. 6 (2): 461-464, from http://www.andrew.cmu.edu/user/kk3n/simplicity/schwarzbic.pdf
Shuey S A, Birkelo C P and Marshall D M 1993 The relationship of the maintenance energy requirement to heifer production efficiency. Journal of Animal Science 71: 2253-2259, from http://www.journalofanimalscience.org/content/71/8/2253.full.pdf+html
Silva F F, Valadares Filho S C, Vinhas I L C, Veloso M C, Diniz V R F, Cecon P R, Rodrigues P P V e Kling de Moraes E B 2002 Composição Corporal e Requisitos Energéticos e Protéicos de Bovinos Nelore, Não-Castrados, Alimentados com Rações Contendo Diferentes Níveis de Concentrado e Proteína. Revista Brasileira de Zootecnia. 31(1): 503-513, from http://www.scielo.br/pdf/rbz/v31n1s0/10333.pdf
Urdapilleta T J, Piva L J F, Kroef T A e dos Santos M G 2005 Relação entre Medidas Ultra-Sônicas e a Espessura de Gordura Subcutânea e a Área de Olho de Lombo na Carcaça em Bovinos de Corte. Revista Brasileira de Zootecnia. 34 (6): 2074-2084, from http://www.scielo.br/pdf/rbz/v34n6/27262.pdf
Wang Z and Goonewardene L A 2004 The use of Mixed models in the analysis of animal experiments with repeated measures data. Canadian Journal of Animal Science 84:1-11, from http://pubs.aic.ca/doi/pdf/10.4141/A03-123
Wang Z, Nkrumah J D, Li C, Basarab J A, Goonewardene L A, Okine E K, Crews Jr, D H and Moore S S 2006 Test duration for growth, feed intake, and feed efficiency in beef cattle using the GrowSafe System. Journal of Animal Science 84: 2289-2298, from http://www.journalofanimalscience.org/content/84/9/2289.full.pdf+html
ZoBell D R, Goonewardene L A, Olson K C, Stonecipher C A and Weidmeier R D 2003 Effects of feeding wheat middlings on production, digestibility, ruminal fermentation and carcass characteristics in beef cattle. Canadian Journal of Animal Science 83 (3): 551-557, from http://pubs.aic.ca/doi/pdf/10.4141/A02-086
Zuur A F, Leno E N, Walker N J, Saveliev A A and Smith G M 2009 Mixed Effects Models and Extensions in Ecology with R. Springer, New York. 553 p, from http://www.highstat.com/book2.htmReceived 17 March 2013; Accepted 1 July 2013; Published 1 August 2013