Livestock Research for Rural Development 31 (10) 2019  LRRD Misssion  Guide for preparation of papers  LRRD Newsletter  Citation of this paper 
The present investigation was carried out to analyze, compare and select the mathematical model of the lactation curve with a better adjustment for crossbred Criollo Limonero dairy cows. A total of 206 complete lactations (≥270 days) were evaluated, for a total of 12,845 daily milk production data from a herd of crossbreds cows Criollo Limonero x White Brahman and Criollo Limonero x Swiss Brown, under tropical dry forest conditions and grazing based feeding. Three empirical mathematical models were used to fit daily milk yield: Wood incomplete gamma function (WOOD), Logarithmic Quadratic (LOGQUAD) and multiple regression model of Ali and Shaeffer (ALI). The goodness of fit of the models were estimated through adjusted multiple coefficient of determination (R^{2}_{Adj} ), residual standard deviation (RSD) and the correlation between the real yield and the estimated milk yield (r). Average daily milk yield (5.41 kg/day), total lactation yield (1437 kg) and lactation length (262.94 days) were described; expected peak yield (6.50 kg/day), expected peak day (56th) and lactation persistence (93.23%) were estimated through mathematical expressions derived from the components of WOOD equation. ALI showed the highest values for R^{2}_{Adj} (0.974) and r (0.987), as well as the lower RSD (0.036), so it was selected as the best fit model to estimate the actual milk yield of Criollo Limonero crossbred cows.
Key words: daily yield, lactation persistence, peak day, peak yield
The cattle known as Tropical Criollo is a type of Bos taurus cattle whose genetic base comes from the biotypes and breeds introduced in America in the 15^{th} century (RosendoPonce and BecerrilPérez 2015). Many of these cattle types are an important part of grazing based milk production systems, as is the case of the Central American Criollos, Hartón del Valle de Colombia, Criollo Lageanos of Brazil, Carora and Criollo Limonero of Venezuela (Tewolde 1997).
High temperatures, solar radiation and relative humidity, as well as a marked seasonality of rain generally characterize the environmental conditions in which these animals have evolved. These situations determine a variable supply of low digestibility forages, with high concentrations of fiber and low protein and energy contents. These extreme conditions have determined a resistance to the tropical environment, that according to Tewolde (1997), confer to the Criollos outstanding characteristics with respect to other genotypes in terms of fertility, apparent resistance to parasites and environmental stress, ease management and carcass quality characteristics (RodasGonzález et al 2006; RodriguezVoigt et al 1997).
Production in dairy cows increases fast from calving to a peak production that is reached in a relatively short period of time and then begins to decrease gradually until drying. Animal genetics, management and environmental conditions mainly determine the natural dynamics of this process. The milk production throughout lactation can be expressed graphically in a lactation curve and can be described by the coefficients of a mathematical model; where the better adjustment to the actual shape, the greater will be the ability to predict milk production of the cattle breed studied (Batra 1986; Palacios et al 2016).
Lactation curve modeling is an important research tool to develop and validate standards in order to explain the main patterns of milk production associated with the stage of lactation, which is of great importance in the organization of management, feeding, selection and reproductive planning in genetic improvement programs (Silvestre et al 2006).
Contreras and Rincon (1978) report, maybe the first recorded data on Criollo Limonero cattle lactation; however, there is an apparent lack of available information on the shape and modeling of their lactation curve.
The present investigation was carried to analyze, compare and select the mathematical model of the lactation curve with a better adjust to productive behavior of crossbred Criollo Limonero dairy cows.
Milk production data analyzed are part of the records of Estación Experimental Local El Guayabo of the Instituto Nacional de Investigaciones Agrícolas (INIA) located in the south west of Venezuela. The herd was conformed of crossbred cows Criollo Limonero x White Brahman and Criollo Limonero x Swiss Brown. A total of 206 complete lactations (more than 270 days) were evaluated, for a total of 12,845 daily milk production data. The animals were hand milked twice a day (05:00 and 13:00), without calf and their performance was recorded.
The region is characterized as tropical dry forest with a precipitation of 1803 mm/year distributed seasonally in a bimodal way, average temperature of 27.9 °C and 91% of relative humidity (VergaraLópez et al 2006).
The feeding system was based on grazing of mixture pastures of Aleman grass (Echinochloa polystachya) and Tanner grass ( Brachiaria arrecta) managed with 1 day of occupation and 28 to 35 days of rest, depending on the season of the year. The average animal stocking rates was 2.0 AU/ha.
Three mathematical models were compared to describe the shape of the lactation curve of crossbred Criollo Limonero cows:
Wood incomplete gamma function (Wood 1967; WOOD), has been widely used to describe lactation curves in many species of ruminants (Quintero et al 2007). It is mathematically expressed with the following equation:
Where y_{t} represents milk production on day t, the parameter a is related to the initial production, b is the slope of the curve before the maximum production peak and c is the slope of the curve after the production peak (Barrios et al 1996).
Where y_{t} is the volume of production measured on day t of lactation, parameters a, b, and c are the parameters to be estimated from the curve, while d is the parameter that multiplies the natural logarithm of days in milk (t ) (Fernández et al 2004).
Ali and Schaeffer (1987) (ALI), described a multiple regression model to calculate the relative efficiencies of selection to change the shape of the lactation curve (Kettunen et al 2000). It is described with the expression:
Where yt is the milk yield on day t, a is a parameter associated with the peak yield, d and e are parameters associated with increasing slope, while b and c are associated with decreasing slope (Silvestre et al 2006).
The models described above were selected for this study, because they were frequently used in the studies of lactations in Criollo breeds and for achieve good adjustments in tropical regions (Fernández et al 2004; Barbosa et al 2007; Hernández and Ponce, 2008).
The goodness of fit of the models was evaluated according to the following criteria suggested by Cankaya et al (2011).
Adjusted multiple coefficient of determination (R^{2}_{adj}):
Residual standard deviation (RSD):
Where n is the number of observations, p is the number of parameters in the model and R^{2}=1−(RSS/TSS) ( RSS: residual sum of squares, TSS: total sum of squares). The coefficient of determination value (R^{2}) is a classic indicator measuring the proportion of total variation about the mean of y_{t} explained by the lactation curve model. However, this tends to increase as elements are added to the model, which can eventually overestimate the fact that the set of elements chosen in the model is adequate to explain the dependent variable. R^{2}_{adj} adjust and penalize the inclusion of coefficients in the model, which make more homogeneous comparisons between models with different parameters.
As an additional comparison parameter, the correlation between the real yield and the estimated milk yield (r) was determined as a way of quantifying the degree of association between the real and estimated values. The selection criteria used in this analysis was to select the model with the highest R^{2}_{adj} and lowest RSD.
Peak yield (PY), peak day on which the maximum yield is expected (PD), as well as the persistence of lactation (P) were determined by means of mathematical expressions derived from the components of the equation of Wood (Cankaya et al 2011; Quintero et al 2007).
The STATGRAPHICS (Ver. 18.1.09) Nonlinear Regression tool to run the models and determine the best fitting using the LevembergMarquardt adjust method was used. The parameters a, b, c, d, and e were estimated by nonlinear least squares through the same tool.
Table 1 shows the general information of the lactations evaluated in this study. The average daily milk production here obtained could be considered low compared to that reported for Criollo Limonero, but relatively higher than reported for other tropical criollos (RosendoPonce and BecerrilPérez 2015).
Table 1. Summary of the production of the lactations evaluated 

Variable 
Mean 
Variation 
SEM 
Day milk yield (kg/day) 
5.41 
11.3 
0.04 
Total lactation yield (kg) 
1437 
41.3 
35.5 
Lactation length (days) 
262.9 
23.6 
3.70 
Perozo et al (1978), mention a mean of 1850 kg of milk per lactations of 280 days, while Contreras and Rincón (1978), report 1642 kg/milk per lactation, and Zambrano et al (2006) recorded 1964 kg/milk. Observed values for total lactation production in this investigation (1437 kg/milk), were lower than those reported. Probably the inclusion of the Bos indicus component in the crossing, seeking resistance to the environment, is affecting the milk production capacity of crossbred Criollo Limonero.
Another aspect to be consider in this discrepancy is the high variability associated with this type of measurements. Several authors report variabilities that can even exceed 50%, attributable to the effect of variations in diet due to environmental seasonality and to the joint evaluation of several lactations (Ribas and Pérez 1989; Fernández et al 2004; Hernández et al 2005).
The length of lactation observed in this study is related to what has been reported in other studies in Criollo Limonero crossbred cows (Zambrano et al 2006). The lactation duration is an important aspect to consider in this type of livestock that produce under unfavorable climatic conditions, because can compromise the persistence of the reproductive characteristics of cows and conditions of the calves at birth (RosendoPonce and BecerrilPérez 2015).
The projected value to peak yield was 6.50 kg/milk/day that is expected to be obtained at the 56th (55.7) day after calving and lactation has an estimated persistence of 93.2%. These estimates, in conjunction with the data reported in Table 1, support that Criollo Limonero crossbreds can express acceptable production levels under extreme conditions of feeding, managing and environment (> 5 kg/animal/day), with prolonged lactations whose production has high persistence after reaching the maximum yield (>90%).
The curve described by the average milk production of the lactations evaluated as a whole (Figure 1), presents the classic behavior described for this variable, with an accelerated initial rise after calving followed by a progressive decay until final drying. (Batra 1986). Fast rise to the PY, followed by a low rate of decline in this curve are associated with the estimated high P and correspond to cows with high productive potential, even more under the feeding conditions without supplementation, a situation that is common in most of tropical milk production systems.
Figure 1. Milk production curve in joint evaluated lactations in crossbred Criollo Limonero cows 
Table 2 and Figure 2 show the estimated parameters for the models evaluated and their graphic representation, respectively.
Table 2. Estimated parameters for the models evaluated 

Models 
Parameters (SEM) 

a 
b 
c 
d 
e 

WOOD 
3.83 (0.19) 
0.18 (0.02) 
0.003 (0.00) 

LOGQUAD 
2.89 (0.31) 
0.028 (0.00) 
0.00003 (0.00) 
1.27 (0.12) 

ALI 
7.11 (2.41) 
3.58 (3.97) 
0.73 (1.13) 
0.47 (1.34) 
0.25 (0.19) 
Figure 2. Actual and empirical estimated lactation curves for crossbred Criollo Limonero cows 
The shape of curves shows that the selected models have high adjustments in the explanation of the variability under study, with a high predictive capacity (>R^{2}) (Table 3). These high values of determination have been widely reported in the literature by different authors (Cankaya et al 2011; Macciotta et al 2005; Vargas et al 2000).
Table 3. Goodness of fit parameters for the evaluated models 

Models 
R^{2} 
R^{2}_{adj} 
RSD 
r 
WOOD 
0.962 
0.961 
0.045 
0.981 
LOGQUAD 
0.972 
0.972 
0.038 
0.986 
ALI 
0.975 
0.974 
0.036 
0.987 
R^{2}: determination coefficient; R^{2} _{
adj}: Adjusted coefficient of 
As can be observed, all the models showed a high relation between the predicted and actual values, with highly significant correlation values (r) in the range of 0.981 and 0.987 (Table 3). Likewise, it is possible to appreciate the values of the adjustment parameters (R^{2}, R^{2}_{adj} and RSD), corroborate the power of models to predict the behavior of the lactation curve.
The analysis to residuals between actual and estimated values is shown in Table 4. The standardized bias and kurtosis values describe a fairly symmetric residual distribution. On the other hand, the data observed for the DurbinWatson statistic confirm the presumption of independence of this, and finally the AndersonDarling statistic was used to verify the normality of the residuals distribution, which allowed confirming this last assumption. Preference was given to the use of this statistic over other more classic ones such as the ShapiroWilk or KolmogorovSmirnovLilliefors tests to check the normality of the residuals because AndersonDarling has been described as a more solid test (Pedrosa et al 2015), especially in small sample sizes, which adjusts to the characteristics of this investigation.
Table 4. Residuals analysis for the evaluated models 

Models 
Standardized 
Standardized 
Anderson 
Durbin 
WOOD 
0.78 
0.44 
0.33 
2.16 
QUADLOG 
0.41 
0.71 
0.22 
2.14 
ALI 
0.71 
0.15 
0.36 
2.06 
a. Values within 2 and 2 express no significant
deviation to normal distribution.

The extreme similarity of the estimated empirical curves can be confirmed with the graphical analysis of the residuals (Figure 3). However, it can be seen how the differences between the models tend to be concentrated around the first 60 days of lactation, where WOOD and QUADLOG shows an important dispersion in their estimation of the production values, making underestimates at the beginning of the period followed by overestimations. This behavior has already been observed by other researchers in incomplete gamma function (Grossman and Koops 1988). In the ALI equation, the inclusion of the estimate at 305 days and the d and e parameters in the multiple regression model improve the estimation since they have a better control in the early phase of the lactation curve (Ali and Schaeffer 1987), improving their prediction performance in the initial values of the curve.
Figure 3. Residuals graphic for actual and estimated milk production values 
The authors want to acknowledge invaluable support for this research from the Instituto Nacional de Investigaciones Agrícolas (INIA, Venezuela), especially to Estación Experimental Local El Guayabo, and all their staff. We also want to thank Misses Angelica R and Andrea L for their help and support in data management.
Ali T E and Schaeffer L R 1987 Accounting for covariances among test day milk yields in dairy cows. Canadian Journal of Animal Science. 67(3), 637644. https://doi.org/10.4141/cjas87067
Barbosa S B P, Pereira R G A, Santoro K R, Batista, A M V and Neto R 2007 Lactation curve of crossbred buffalo under two production systems in the Amazonian region of Brazil. Italian Journal of Animal Science. 6(sup2), 10751078. https://doi.org/10.4081/ijas.2007.s2.1075
Barrios A, Rincón E, Ventura M, Huerta N, Fondevilla M y Surra J 1996 Factores que afectan la curva de lactancia de vacas mestizas en regiones tropicales. Revista de la Facultad de Agronomía (LUZ). 13(6), 741749.
Batra T R 1986 Comparison of Two Mathematical Models in Fitting Lactation Curves for Pureline and Crossline Dairy Cows. Canadian Journal of Animal Science. 66(2), 405414. https://doi.org/10.4141/cjas86042
Cankaya S, Unalan A and Soydan E 2011 Selection of a mathematical model to describe the lactation curves of Jersey cattle. Archives Animal Breeding. 54(1), 2735. https://doi.org/10.5194/aab54272011
Contreras R y Rincón E 1978 Curvas de lactancia de vacas Criollo Limonero en un ambiente de trópico húmedo. Revista de la Facultad de Agronomía (LUZ). 5(2), 458467.
Fernández L, Menéndez A y Guerra C 2004 Estudio comparativo de diferentes funciones para el análisis de la curva de lactancia en el genotipo Siboney de Cuba. Revista Cubana de Ciencia Agrícola. 38(4), 351359.
Grossman M, and Koops W J 1988 Multiphasic Analysis of Lactation Curves in Dairy Cattle. Journal of Dairy Science. 71(6), 15981608. https://doi.org/10.3168/jds.S00220302(88)797234
Hernández A, Ponce de León R, Gutiérrez M, García S, García R, Mora M y Guzmán G 2005 Efectos ambientales en la producción lechera de la raza bovina Mambí de Cuba. Revista Cubana de Ciencia Agrícola. 39(4). Retrieved from: http://www.redalyc.org/resumen.oa?id=193017719001
Hernández R y Ponce P 2008 Caracterización de la curva de lactancia y componentes lácteos del genotipo siboney de Cuba en una granja ganadera de la provincia de La Habana. Revista Científica, FCVLUZ. 18(3), 291295.
Kettunen A, Mäntysaari E A and Pösö J 2000 Estimation of genetic parameters for daily milk yield of primiparous Ayrshire cows by random regression testday models. Livestock Production Science. 66(3), 251261. https://doi.org/10.1016/S03016226(00)001664
Macciotta N P P, Vicario D and CappioBorlino A 2005 Detection of Different Shapes of Lactation Curve for Milk Yield in Dairy Cattle by Empirical Mathematical Models. Journal of Dairy Science. 88(3), 11781191. https://doi.org/10.3168/jds.S00220302(05)727843
Palacios A, GonzálezPeña F, Guerra, D, Espinoza, J L, Ortega, R, Guillén, A y Ávila N 2016 Curvas de lactancia individuales en vacas Siboney de Cuba. Revista Mexicana de Ciencias Pecuarias. 7(1), 1528.
Pedrosa I, JuarrosBasterretxea J, RoblesFernández A, Basteiro J y GarcíaCueto E 2015 Pruebas de bondad de ajuste en distribuciones simétricas? qué estadístico utilizar? Universitas Psychologica. 14(1), 245–254.
Perozo N, Labbe S, Abreu O y Díaz E 1978 Producción de leche de ganado criollo venezolano. Agronomía Tropical. 28(3), 205220.
Quintero J C, Serna J I, Lugo N A H, Solano R N y Muñoz M F C 2007 Modelos matemáticos para curvas de lactancia en ganado lechero. Revista Colombiana de Ciencias Pecuarias. 20(2), 149156.
Ribas M y Pérez B 1989 Pesajes mensuales de leche y la producción a los 244 días. Efectos ambientales en primera lactancia. Revista Cubana de Ciencia Agrícola. 23(1), 15.
RodasGonzález A, VergaraLópez J, Arenas de Moreno L, HuertaLeidenz N y Pirela M 2006 Características al sacrificio, rasgos de la canal y rendimiento carnicero de novillos Criollo Limonero sometidos a suplementación durante la fase de ceba a pastoreo. Revista Científica. FCVLUZ. 16(4), 371380.
RodriguezVoigt A, Noguera E, Rodríguez H L, HuertaLeidenz N O, MorónFuenmayor O y RincónUrdaneta E 1997 Crossbreeding dualpurpose cattle for beef production in tropical regions. Meat Science. 47(3), 177185. https://doi.org/10.1016/S03091740(97)000442
RosendoPonce A y BecerrilPérez C M 2015 Avance en el conocimiento del bovino criollo lechero tropical de México. Ecosistemas y recursos agropecuarios. 2(5), 233243.
Silvestre A M, PetimBatista F and Colaço J 2006 The Accuracy of Seven Mathematical Functions in Modeling Dairy Cattle Lactation Curves Based on TestDay Records From Varying Sample Schemes. Journal of Dairy Science. 89(5), 18131821. https://doi.org/10.3168/jds.S00220302(06)722500
Singh R P and Gopal R 1982 Lactation curve analysis of buffaloes maintained under village conditions. Indian Journal of Animal Sciences. Retrieved from: http://agris.fao.org/agrissearch/search.do?recordID=US201302588751
Tewolde A 1997 Los Criollos bovinos y los sistemas de producción animal en los trópicos de América Latina. Archivos Latinoamericanos de Producción Animal. 5 (Supl. 1), 1319.
Vargas B, Koops W J, Herrero M and Van Arendonk J A M 2000 Modeling Extended Lactations of Dairy Cows. Journal of Dairy Science. 83(6), 13711380. https://doi.org/10.3168/jds.S00220302(00)750053
VergaraLópez J, RodríguezPetit A, Navarro C y Atencio A 2006 Efecto de la suplementación con leucaena (Leucaena leucocephala Lam. De wit) sobre la degradabilidad ruminal del pasto alemán ( Echinochloa polystachya H.B.K. Hitch). Revista Científica, FCVLUZ. 16(6), 642647.
Wood P 1967 Algebraic Model of the Lactation Curve in Cattle. Nature. 216(5111), 164. https://doi.org/10.1038/216164a0
Zambrano S, Contreras G, Pirela M, Cañas H, Olson T and LandaetaHernández A 2006 Milk yield and reproductive performance of crossbred Holstein x criollo limonero cows. Revista Científica. FCVLUZ, 16(2), 155164.
Received 21 July 2019; Accepted 6 September 2019; Published 2 October 2019