Comparison of different lactation curve models in Holstein cows raised on a farm in the south-eastern Anatolia region

The objective of this study was to compare the goodness of fit of seven mathematical models (including the gamma function, the exponential model, the mixed log model, the inverse quadratic polynomial model and their various modifications) on daily milk yield records. The criteria used to compare models were mean R, root mean squared errors (RMSE) and difference between actual and predicted lactation milk yields. The effect of lactation number on curve parameters was significant for models with three parameters. Third lactation cows had the highest intercept post-calving, greatest incline between calving and peak milk yield and greatest decline between peak milk yield and end of lactation. Latest peak production occurred in first lactation for all models, while third lactation cows had the earliest day of peak production. The R values ranged between 0.590 and 0.650 for first lactation, between 0.703 and 0.773 for second lactation and between 0.686 and 0.824 for third lactation, depending on the model fitted. The root mean squared error values of different models varied between 1.748 kg and 2.556 kg for first parity cows, between 2.133 kg and 3.284 kg for second parity cows and between 2.342 kg and 7.898 kg for third parity cows. Lactation milk yield deviations of Ali and Schaeffer, Wilmink and Guo and Swalve Models were close to zero for all lactations. Ali and Schaeffer Model had the highest R for all lactations and also yielded smallest RMSE and actual and predicted lactation milk yield differences. Wilmink and Guo and Swalve Models gave better fit than other three parameter models.


Introduction
The graphical representation of daily milk yield against time after calving is called a lactation curve (SHERCHAND et al., 1995).Knowledge of the shape of the lactation curve in dairy cattle is useful in making management and breeding decisions (WILMINK, 1987;SWALVE, 1995;OLORI et al., 1999;SWALVE and GUO, 1999;BOHLSEN et al., 2003;AMIN, 2003).Feeding management programs might be arrange according to the shape of curve.While the rising portion of the curve suggests that a higher plane of nutrition should be supplied to cows, a declining portion of the curve suggests a lower plane of nutrition (SHERCHAND et al., 1995;TEKERLI et al., 2000).Furthermore, mathematical models estimated from the milk yield data may be used to predict future milk yields of an individual cow for the purpose of culling (SHERCHAND et al., 1995).Daily milk production in dairy cattle increases to a peak between 40 and 70 days post partum and declines afterward until the cow is dried-off (BAFFOUR-AWUAH et al., 1996;OLORI et al., 1999;SWALVE and GUO, 1999).Different mathematical models have been developed to describe the shape and type of lactation curve and to predict lactation milk yield (ALI and SCHAEFFER, 1987;OLORI et al., 1999;HORSTICK and DISTL, 2002;GRZESİAK et al., 2003).The incomplete gamma function first proposed by WOOD (1967) is the most popular model used to describe lactation curve.However, ROWLANDS et al. (1982) reported that Wood's curve don't give a consistently good fit to data from high yielding dairy cows around the time of peak milk yield.ALI and SCHAEFFER (1987) reported that a regression model of yields on day in lactation (linear and quadratic) and on log of 305 divided by day in lactation (linear and quadratic) were better than gamma function for predicting lactation milk yield.GROSSMAN et al. (1986) added sine and cosine terms to incomplete gamma function in order to account for seasonal variations.On the other hand, the use of more complicated nonlinear models to describe lactation curves has some disadvantages.Firstly, they require greater amount of data to estimate the lactation parameters accurately, and secondly, larger dataset requires additional mathematical techniques and may cause more computational problems.The objective of this study is to compare seven mathematical models using daily milk yield records by lactation number for a Holstein herd.

Materials and Methods
The dataset used in this study was obtained from a farm in south-eastern Anatolia region of Turkey from January 2000 to December 2003.The cows were housed in semi-open free-stall barns and were milked three times per day in a milking parlour.The management and feeding systems applied to cows during the study were the same.The data on daily milk yield was transferred to computer automatically by using Afimilk Meters.The farm has a capacity of 1,000 dairy cows and the movement of the cows, production status and the visits of the veterinarians are being controlled in a computerized herd management program.In the milking parlour a carousel milking system for 50 cows is being used.The climate in the south-eastern Anatolia region is continental.According to General Directorate of Meteorological Services of Turkey, environmental temperatures ranged between 16.0 and 46.8 in summer months, and ranged between -3.3 and 20.5 in winter months during the dates of the study.As the weather is very hot during summertime sprinklers and ventilators are being used in the farm during the hot hours of the day.In the region aridity is severe.There is nearly no rainfall in July and August and the maximum precipitation falls in winter.
The initial dataset included 1,299 lactations of 976 Holstein cows.The data of cows having a lactation less than 305 days, were removed from the dataset.On the other hand, when the lactation of a cow was longer than 305 days, only 305 day milk yield was used in the analyses.The data of cows having illnesses, which affect milk yield (mastitis, lameness, ketosis, left displaced abomasum) were also excluded from the analyses.Finally, the edited dataset used for analyses included information from 330 first, 106 second and 54 third lactations of 432 cows.In the current study, seven models describing lactation curve were studied.The first model was WOOD'S (1967) incomplete gamma function given below: y t = a 1 t b1 e -c1t (1) where: y t = the daily milk yield on day t, a 1…6 = a parameter associated with the initial milk yield just after calving; b 1…6 and c 1…6 = parameters associated with the increasing and decreasing slopes of the lactation curve, respectively.The second model was the exponential model of WILMINK (1987): (2) The third model was the mixed log model of GUO and SWALVE (1997): y t =a 3 + b 3 t 0.5 + c 3 ln (t) (3) The inverse quadratic polynomial model of NELDER (1966) as applied to dairy cattle by YADAV et al. (1977) was: y t -1 = a 4 + b 4 t -1 + c 4 t (4) The model of GOODALL (1986), which takes into account the influence of season to the lactation curve, was: y t = a 5 t b5 e (-c5t+d5D)  (5) where d 5 = coefficient of season, D = 0 for October to March; D = 1 for April to October.The Grossman Model, which was modified from the Wood's Model by including sine and cosine terms to account for seasonal variations (GROSSMAN et al., 1986), was: y t =a 6 t b6 e -c6t (1 + u sin (x)+ v cos (x)) (6) where x = the day of year of the daily yield calculated as radian, u and v = the coefficients of the day of year.The final model was a regression model of yields on day in lactation (linear and quadratic) and on log of 305 day yield divided by day in lactation (linear and quadratic) fitted by ALI and SCHAEFFER (1987): y t =a 7 + b 7 (t/305) + c 7 (t/305) 2 + d 7 ln (305/t) + f 7 ln 2 (305/t) (7) where t = denotes day in milking, y t = denotes daily milk yield on day t; a 7 , b 7 , c 7 , d 7 and f 7 = the regression coefficients where a 7 is associated with peak yield, d 7 and f 7 are associated with the increasing slope of curve and b 7 and c 7 are associated with the decreasing slope for this model.Computation of lactation parameters was done with STATISTICA using the Quasi-Newton method.A one-way ANOVA method in SPSS 10.0 program was applied to analyse the effect of parity on lactation parameters of each model (SPSS, 2004).The criteria used to compare models were mean R 2 (the square of correlation coefficient between actual milk yield and predicted milk yield according to the model), root mean squared errors (RMSE) and difference between actual and predicted lactation milk yields.Models resulting in smaller RMSE and higher R 2 were considered to be superior.

Results and Discussion
The mean parameter estimates for first, second and third lactations are presented in Table 1 for models with three parameters and table 2 for models with four or five parameters.The effect of lactation number on curve parameters was significant for models with three parameters (P<0.05 for b 1 parameter, and P<0.01 for other parameters).This effect was also significant (P<0.01) for parameters of the Goodall Model, except d 5 parameter, which is the coefficient of season.In the Grossman Model, the effect of lactation number on the milk yield incline between calving and peak milk production (b 6 ) and on the coefficient of calving season (v) were not significant, while it was significant for milk yield at the beginning of lactation (a 6 ), decline in milk production from peak to the end of lactation (c 6 ) and coefficient of year of the yield obtained as measured radian (u).In the Ali and Schaeffer Model, this effect was significant only on parameters of d 7 and f 7 .According to the results of Wood, Wilmink, Goodall and Grossman Models, third lactation cows had the highest intercept post-calving (a 1 , a 2 , a 5 and a 6 , respectively).Third lactation cows also had the greatest incline between calving and peak milk yield (b 1 , b 2 , b 5 and b 6 ) and greatest decline between peak milk yield and end of lactation (c 1 , c 2 , c 5 and c 6 ).Similar results were reported by HORAN et al. (2005) and WILMINK (1987).However, some authors (KAYGISIZ, 1999;ORMAN and ERTUGRUL, 1999;TEKERLI et al., 2000) reported non-significant effects of lactation number on lactation parameters using the Wood Model.
Predicted peak milk yield and day of peak milk yield by lactation number for different models are given in table 3.There were considerable differences between models for the estimates of peak milk yield and day of peak yield.Wood Model had lower peak milk yield compared to other six models for lactation 1, 2 and 3.The highest predicted peak milk yield for lactations 1, 2 and 3 was obtained by the Nelder-Yadav Model.Nelder-Yadav Model had an earlier peak compared to other models for lactations 2 and 3.For first lactation Wood Model had an earlier peak milk yield.
There was a tendency for peak milk yield to increase with lactation number for all models analysed in the current study.Similar trends for peak milk production were reported by earlier studies (SWALVE and GUO, 1999;TEKERLI et al., 2000;HORAN et al., 2005).Latest peak production occurred in first lactation for all models, while third lactation cows had the earliest day of peak production.The tendency for peak production time to become later as the lactation number increased found in the current study, was in agreement with reports of Keown et al. (1986) and Tekerli et al. (2000) for Holstein cows.Rao and Sundaresan (1979) reported that latest day of peak production in first lactation cows might be explained by the milk secretory tissue in primiparous animals taking longer to reach its peak activity than in pluriparous animals.In the study of Horan et al. (2005), however, day of peak milk yield was reported to be significantly shorter in second parity cows compared to both first and third parity cows.The coefficient of determination (R 2 ) and root mean square errors (RMSE) of models by lactation number are presented in table 4. The R 2 values ranged between 0.590 and 0.650 for first lactation, between 0.703 and 0.773 for second lactation and between 0.686 and 0.824 for third lactation, depending on the model fitted.Guo and Swalve (1997) reported the R 2 values obtained from three different farms between 0.549 and 0.719 by the Wood Model, between 0.567 and 0.736 by the Wilmink Model, between 0.569 and 0.738 by the Guo and Swalve Model, and between 0.558 and 0.703 by the Ali and Schaeffer Model.Orman and Ertuğrul (1999) reported R 2 estimates by the Wood Model to be between 0.73 and 0.79, depending on the parity.The R 2 estimates of Soysal and Gürcan (2000) from two different farms, were 0.44 and 0.89 for the Wood Model, 0.45 and 0.82 for the Goodal Model and 0.60 and 0.94 for the Grossman Model.The R 2 estimates reported by Olori et al. (1999), were 0.66 by the Wood Model, 0.65 by the Nelder-Yadav Model, 0.67 by the Guo and Swalve Model and 0.69 by the Wilmink Model.Estimates of R 2 obtained in the current study for different Models were generally in the range of these literature estimates.According to R 2 values, a better fit was found for second and third lactations than for first lactation.Wilmink (1987) also reported that R 2 values of first parity cows were lower than in later parities.On the other hand the same tendency was not found in the study of Orman and Ertuğrul (1999).
Based on the R 2 values, the Ali and Schaeffer Model gave the best fit to the current dataset for first, second and third lactations (0.650, 0.773 and 0.824, respectively).QUINN et al. (2004) also reported that the Ali and Schaeffer Model gave better fit than the Wood, Wilmink, and Guo and Swalve Models.Contrary to the current study, the R 2 estimates from the Wood, Nelder-Yadav, Guo and Swalve and Wilmink Models were reported to be similar by OLORI et al. (1999).In the study of SOYSAL and GÜRCAN (2000), who compared the Wood, Goodal and Grossman Models, the highest R 2 value was reported for the Grossman Model.The root mean square error values of different models varied between 1.748 kg and 2.556 kg for first parity cows, between 2.133 kg and 3.284 kg for second parity cows and between 2.342 kg and 7.898 kg for third parity cows.As lactation number increased, RMSE values also increased in all the models.Similar result was also reported by ORMAN and ERTUĞRUL (1999) andSWALVE (1997) andGRZESIAK et al. (2003) also reported that the differences between estimated and true lactation milk yield in the Ali and Schaeffer, Wilmink, and Guo and Swalve Models were similar, and the superiority of these models to the Wood Model were significant.The differences between lactation number groups by models were generally small.First lactation cows had the highest differences for all models.
According to actual and predicted lactation milk yield differences, although the Wilmink, Guo and Swalve, and Ali and Schaeffer Models could be used to describe lactation milk yield, the Ali and Schaeffer Model could be chosen as it had the smallest lactation milk yield difference and it also gave the best fit for R 2 and RMSE.The Ali and Schaeffer Model had the highest R 2 for all lactations and also yielded the smallest RMSE and actual and predicted lactation milk yield differences.Hence it should probably be considered to describe the lactation curve of cows.The Goodal and Grossmann Models did not provide any advantage for R 2 , RMSE and actual and predicted lactation milk yield differences compared with three parameter models.
Although R 2 values of Nelder-Yadav Model were sufficiently high compared with other models, this model yielded considerably higher RMSE especially for lactation 3. The Wood Model had higher RMSE and actual and predicted lactation milk yield differences compared to the Wilmink and Guo and Swalve Models.Therefore, the Wilmink and Guo and Swalve Models could be chosen as the better models amongst the three parameter models.

Table 1
Mean values and standard errors of lactation parameters for three parameter models (Mittelwerte und Standardfehler von Laktationsnummern für die Drei-Parameter-Modelle) a, b, c = means with different superscript within a line was significantly different (* P<0.05, ** P<0.01), SE = standard error

Table 4
Squared correlations and root mean square errors of models by lactation number (R 2 und RMSE von Modellen nach Laktationsnummern) GUO and SWALVE (1997)00)he Ali and Schaeffer Model provided the smallest RMSE for all parities, whileNelder-Yadav Model had highest RMSE.OLORI et al. (1999)also reported smaller RMSE in the Ali and Schaeffer Model than theWood, Nelder-Yadav, Wilmink and  Guo and Swalve Models.QUINN et al. (2004)reported that the Ali and Schaeffer Model gave the best fit to their data compared to the Wilmink and Guo and Swalve Models, based on the mean square prediction error.SOYSAL and GÜRCAN (2000), who compared the goodness of fit of the Wood, Goodal and Grossmann Models, reported the highest mean square error by the Goodal Model, while the Grossmann Model had smallest.However, in the study ofGUO and SWALVE (1997), absolute error computed from difference of true and estimated yield per test day by the Guo and Swalve Model, were reported to be smaller than those estimates by the Wood, Wilmink, and Ali and Schaeffer Models.Actual lactation (305-day) milk yield and predicted lactation milk yield deviations (predicted lactation milk yield minus actual lactation milk yield) for models investigated in the current study, are presented in Table5.Lactation milk yield deviations of the Ali and Schaeffer, Wilmink, and Guo and Swalve Models were close to zero for all lactations.However, predicted minus actual lactation milk yield differences were between 316.77 and 461.63 kg in the Wood Model, -402.54 and -355.45kg in the 512.75 and 289.50kg in the Goodal Model and 329.34 and 238.07 kg in the Grossmann Model, depending on the parity of cows.GUO

Table 5
Actual lactation milk yield and differences between predicted and actual lactation milk yields according to different models (kg) (Reale Laktationsmilchleistung und Unterschiede zwischen realen und vorhergesagten Laktationsmilchleistungen nach den verschiedenen Modellen)