INTRODUCTION
Lactation curves are valuable tools for dairy producers for management decision-making and selection. Knowledge of the lactation curve allows prediction of total milk yield from a single or several test days early in lactation. With such a knowledge, a dairy producer can make management decisions early based on individual production (Gipson & Grossman Reference Gipson and Grossman1990). There are various mathematical equations describing lactation curves in dairy cows, from the more empirical equations that relate input to output statistically with little consideration of the biology of lactation (e.g. Wood Reference Wood1967; Rook et al. Reference Rook, France and Dhanoa1993), to the more mechanistic ones which describe the lactation curve based on the biology of lactation (e.g. Dijkstra et al. Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997; Pollott Reference Pollott2000). Although the latter may produce parameters that have a more biological interpretation, they can be too complex for routine use outside research, and some of the mechanisms represented are as yet not fully understood. Growth functions which have been used to describe growth in various animal (Lopez et al. Reference Lopez, France, Gerrits, Dhanoa, Humphries and Dijkstra2000) and microbial (Lopez et al. Reference Lopez, Prieto, Dijkstra, Dhanoa and France2004) species can also be used as simple empirical equations to describe cumulative yield of milk over lactation by analogy with growth of body tissues and microbial populations. Differentials of these functions would give an expression for daily milk yield.
The objective of the current paper, therefore, was to compare the performance of standard growth functions (logistic, Gompertz, Schumacher and Morgan) in their differential form with the Wood equation (the most widely used lactation equation) and the Dijkstra equation (a simple mechanistic model) by fitting these equations to monthly milk records for an entire lactation from a commercial herd of Iranian Holstein cows.
MATERIALS AND METHODS
Data
For the study, data on 215 completed lactations by Iranian Holstein cows were obtained from a commercial herd belonging to Kenebist farm, Mashad, Iran. Parity is a factor that has been cited as affecting milk production and characteristics of the lactation curve (e.g. Pérochon et al. Reference Pérochon, Coulon and Lescourret1996). Therefore, the cows studied were of different parities (70 first, 70 second and 75 third). These animals were chosen as being representative of the herd and had undergone a minimum of 11 months of lactation and were free of any health disorders over this period. For uniformity of duration and to remove the effect of drying-off, only the first 11 months of lactation were considered. The effect of season was controlled by including only winter calving cows. Once-monthly observation of milk production was used for each lactation. The ingredients and composition of diets fed to the cows are presented in Table 1. Summary statistics of the lactational data by parity are shown in Table 2.
Table 1. Diet ingredients (g/kg DM) and composition
a Bergafat-F100 (palm oil), component (DM basis): 20 g glycerol/kg and 750 g free fatty acids/kg (myristic acid <35, palmitic acid 750–900, stearic acid 50–100, oleic acid 50–100 and linoleic acid <10 g/kg total fatty acids).
b Composition of trace minerals (g/kg DM) and vitamin mix: Ca, 196; P, 96; Mg, 19; Fe, 3; Na, 71; Cu, 0·3; Mn, 2; Zn, 3; Co, 0·1; I, 0·1; Se, 0·01; and Vit A, 500 000 IU/kg DM; Vit D, 100 000 IU/kg DM; Vit E, 100 IU/kg DM.
c Estimated from NRC (2001).
d N-free basis.
e Non-fibre carbohydrate=1000−(crude protein+neutral detergent fibre+ether extract+ash).
Table 2. Summary of the observed lactation yield data (kg) by parity
y avg, average yield; y 0, initial yield (kg/day); y m, peak yield (kg/day); Y, total yield; Y 305, 305-day yield.
Equations
The equations used to describe the lactation curves are presented in Table 3. The incomplete gamma function proposed by Wood (Reference Wood1967), which has been used widely to study lactation curves, was selected as the simplest equation. Although this equation can give an acceptable fit to milk yield data from a given lactation, it tends to over-predict during early and late lactation and under-predict data during mid-lactation (Grossman & Koops Reference Grossman and Koops1988; Sherchand et al. Reference Sherchand, McNew, Kellogg and Johnson1995). The logistic, Gompertz, Schumacher and Morgan equations are growth functions (Thornley & France Reference Thornley and France2007). A typical lactation curve rises to a peak before falling away, which is the same trajectory mapped by the slope of a sigmoidal growth function. Therefore, growth functions written in their differential form and expressed as a function of time have potential application as lactation equations and are simple to use. The equation of Dijkstra et al. (Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997) is the solution to a simple mechanistic model based on a set of differential equations representing cell proliferation and cell death in the mammary gland and is a four-parameter algebraic equation. Although the latter is based on growth biology of the mammary gland and is more attractive as it is more realistic physiologically, the model with its four parameters is not necessarily simple enough for routine statistical use.
Table 3. Equations used to describe the lactation curve
* y is the milk yield (kg/day), t is the time from parturition (day), and a, b, c and d (all >0) are parameters that define the scale and shape of the lactation curve.
For all equations, expressions for initial yield, time to peak, maximum yield, total yield over a finite lactation and relative rate of decline (persistency) at the midway point of the declining phase were derived analytically (see Thornley & France, Reference Thornley and France2007, p. 620 et seq., by way of example) and are presented in Table 4. These features of lactation were calculated for the different equations by using the related formulae, after parameter estimation. To determine whether an equation predicts initial, maximum and total yields accurately or if there is under/overprediction, the difference between predicted and observed (calculated from actual records) values were t-tested. The observed total yield was the yield that accumulated during the 11 months of lactation and was calculated by multiplying the milk yield observed with the number of days until next observation. The observed maximum yield was the maximum value of the 11 observations of each animal.
Table 4. Expressions for the features of the lactation curve for each equation
* Definition of terms: y 0, initial yield (kg/day); t m, time to peak yield (days); y m, peak yield (kg/day); t f, length of lactation (days); t h=(t m+t f)/2; r(t h), relative rate of decline at the point halfway between peak yield and end of lactation; Y, total yield (kg); γ(q, x), incomplete gamma function=
.
Statistical analyses
Each equation was fitted to monthly records of the 330-day lactation of each cow (total of 215 curves) using the PROC NLIN procedure in SAS (SAS Institute 1999) and the parameters were estimated. The models were tested for goodness of fit (quality of prediction) using a runs test, residual sum of squares (RSS), and calculation of mean square prediction error (MSPE), root of MSPE (RMSPE) expressed as a proportion of the observed mean (Theil Reference Theil1966), adjusted multiple coefficient of determination (R adj2) and Akaike's information criterion (AIC). In a runs test, a run is a sequence of residuals with the same sign (positive or negative), and for this test, the average residual of replicate observations was used for each day in milk. Clustering of residuals with the same sign and serial correlation may be indicative of inappropriate fitting of the model to experimental data or choice of model and a small number of runs of sign is obtained when the residuals are not randomly distributed, so residuals of the same sign tend to cluster together on some parts of the curve (Motulsky & Ransnas Reference Motulsky and Ransnas1987). Mean square prediction errors were calculated as the sum of squared differences between observed and predicted values divided by the number of experimental observations. RMSPE was also calculated (square root of MSPE divided by the observations mean) so that MSPE could be expressed in the same units as the observed and predicted values. The MSPE was broken down into mean bias or error in central tendency (ECT); slope bias or error due to regression (ER); and random error or error due to disturbance (ED). Error in central tendency indicates how the average of predicted values deviates from the average of observed values. Error due to regression (regression bias) measures deviation of the least squares regression coefficient from 1, the value it would have been if the predictions were completely accurate. Error due to disturbance represents the variation in observed values unexplained after the mean and the regression biases have been removed. R adj2 was calculated using the formula:
where R 2 is the multiple coefficient of determination [1−(RSS/TSS)], TSS is total sum of squares, n is the number of observations (data points) and p is the number of parameters in the equation. Note that R 2 was adjusted for the number of parameters in the equation to make a fair comparison between models. For simplicity, R adj2 will be reported as R 2. AIC was calculated as using the equation (Burnham & Anderson Reference Burnham and Anderson2002):
AIC is a good statistic for comparison of models of different complexity because it adjusts the RSS for number of parameters in the model. A smaller numerical value of AIC indicates a better fit when comparing models. All the above statistics were analysed for variance using the PROC GLM procedure of SAS (SAS Institute 1999) using the following model:
where Z ijk represents the different statistics, μ is the overall mean, M i is the fixed effect of lactation equation, P j is the fixed effect of parity, MP ij is the fixed interaction effect of lactation equation i with parity j, and e ijk is the random effect (i=1–6; j=1–3; k=1–70 for parities 1 and 2 but 1–75 for parity 3). The least square means (LSM) for each effect were compared using the Duncan test.
RESULTS AND DISCUSSION
In the present work, comparison of the equations was carried out according to three criteria: model behaviour when fitting the curves using nonlinear regression, statistical performance, and comparison of biologically meaningful parameters estimated by each equation.
Model behaviour
The relationship between predicted milk yield and days in milk (DIM) for first, second and third parity cows was investigated using six models as shown in Fig. 1. For comparison with predicted values, observed milk yield was plotted against DIM (Fig. 2) and the residuals were plotted against DIM to identify difficulties in prediction by a particular model (Fig. 3). When fitting curves by nonlinear regression using the six equations and utilizing the PROC NLIN procedure of SAS (SAS Institute 1999), non-convergence to a solution was observed with all equations, but its occurrence differed. When using growth functions, it was low (0·066, 0·08, 0·013 and 0·03 proportion of total curves for logistic, Gompertz, Schumacher and Morgan, respectively), but for the Wood and especially Dijkstra equations it was higher (0·096 and 0·305 of total curves, respectively). The range of possible starting values for parameters was selected by inspecting the final estimates for each parameter obtained with curves that were fitted without problems, after checking the uniqueness of these solutions. This method of selecting initial parameter estimates for the nonlinear modelling process was judged to be appropriate to achieve a reliable solution for most data sets. Following this approach, further sensitivity to starting values and convergence after a large number of iterations, which are considered symptomatic of an ill-conditioned (parameter values tending to 0 or to biologically unacceptable values) or inappropriate model, were observed more often for Wood and Dijkstra (0·09 and 0·11, respectively) than for the other equations. Thus, the results of model behaviour analysis showed superior ability of differentials of growth functions for fitting the lactation data used in the present study. However, it should be noted that there is not a single simple method to evaluate similarities and differences between nonlinear equations and to deal with the question of which equation should be used; therefore, selection of a model to explain a particular set of data should not be based entirely on model behaviour. Statistical measures are also important (Motulsky & Ransnas Reference Motulsky and Ransnas1987) and are presented in the following section.
Fig. 1. Plots of predicted milk yield (kg/day) against DIM for cows of first, second and third parity using the Wood (△), logistic (■), Gompertz (□), Schumacher (•), Morgan (×) and Dijkstra (▲) equations.
Fig. 2. Plots of observed milk yield (kg/day) against DIM for cows of first, second and third parity.
Fig. 3. Plots of residuals against DIM for cows of first, second and third parity using Wood (△), logistic (■), Gompertz (□), Schumacher (•), Morgan (×) and Dijkstra (▲) equations (residuals are calculated as predicted minus observed values).
Statistical evaluation
There are different statistical tests for ranking and evaluating models. Sometimes results from these different tests seem contradictory, so an overall assessment is needed in this situation (Fathi Nasri et al. Reference Fathi Nasri, Danesh Mesgaran, France, Cant and Kebreab2006). The number of runs of sign of residuals, RSS, MSPE, breaking down of MSPE into ECT, ER and ED, RMSPE (as an indicator of model accuracy), R 2 and AIC, as the most widely used statistical criteria for comparing models, were calculated and the results are shown in Tables 5 and 6. Number of runs of sign of residuals from fitting the different equations (Table 5) tended to be small for the Wood, Schumacher and Morgan equations with a high percentage of curves with five or fewer runs, indicating systematic under- or over-estimation. With all the other equations, the number of runs tended to be higher as a result of random distribution of the residuals over time. The number of runs of sign tended to increase with all equations as parity increased, which means the equations were able to fit lactation records of multiparous cows better than primiparous cows. Additionally, the other statistics (apart from RSS) showed there was significant difference between parities, so the equations fitted data from multiparous cows better than primiparous cows. Differences between the characteristics of the lactation curve of primiparous and multiparous cows are likely to be responsible for the significant difference between goodness of fit of the equations for the different parities. Primiparous cows are more persistent, and their lactation curves are flatter than those of multiparous cows, sometimes with little discernible rise to a peak (Wood Reference Wood1969). Scott et al. (Reference Scott, Yandell, Zepeda, Shaver and Smith1996) pointed out that lactation curves of multiparous cows are similar, except that the estimate of daily milk production is multiplied by a slightly higher factor for cows in their fourth lactation and beyond. There was no significant difference between the goodness of fit for cows of parities 2 and 3 in the present study.
Table 5. Proportion of curves for each number of runs of sign of the residualsFootnote * (in the range ⩽4 to ⩾8) observed when fitting each lactation equation (smaller number of runs indicates systematic under- or over-fitting)
* Residuals are calculated as predicted minus observed values.
Table 6. Some features and statistics for comparing goodness of fit of the different lactations (standard errors are in parentheses)
* Equations 1–6: Wood, logistic, Gompertz, Schumacher, Morgan and Dijkstra, respectively.
† Degrees of freedom.
‡ y 0, initial yield (kg/day); t m, time to peak yield (days); y m, peak yield (kg/day); r(t h) and r(t h)305, relative rate of decline at the point halfway between peak yield and end of lactation and between peak yield and 305 days of lactation, respectively; Y and Y 305, total yield and 305-day yield (kg), respectively.
§ RMSPE, square root of MSPE expressed as a percentage of the observed mean; R 2, adjusted multiple coefficient of determination; RSS, residual sum of squares; MSPE, mean square prediction error; ECT, error due to central tendency; ER, error due to deviation from regression line; ED, random error; AIC, Akaike information criteria.
After fitting each equation to raw milk yield data of cows at different parities (Fig. 1), LSM of the statistics was calculated (Table 6). None of these statistics was able to discriminate between the goodness of fit of different equations, so despite the results of model behaviour analysis, all equations were of the same rank or quality in the fitting of the lactation data. These results are in contrast to those of Dijkstra et al. (Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997) and Val-Arreola et al. (Reference Val-Arreola, Kebreab, Dijkstra and France2004), who found that the Dijkstra equation fitted better than more empirical equations such as the Wood equation. One of the main differences between this and previous studies is the frequency of measurement. Both the Dijkstra et al. (Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997) and Val-Arreola et al. (Reference Val-Arreola, Kebreab, Dijkstra and France2004) studies had higher frequency and number of cows in the analysis, which might have helped in finding significant differences between the equations they evaluated. The results, however, are in accordance with those of Pollot (Reference Pollott2000), who concluded (based on many dairy recording schemes which use a maximum of ten monthly test-day records per lactation) that the equations with more parameters such as Dijkstra et al. (Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997) are over-parameterized and will not yield better fits than Wood. This is because in these situations it is common to find no recordings, or at best one recording, taken before peak yield, and this makes estimating the cell proliferation phase of the lactation and the peak yield less accurate. Both Dijkstra et al. (Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997) and Val-Arreola et al. (Reference Val-Arreola, Kebreab, Dijkstra and France2004) used many more observations in their studies than just once per month, so it appears that frequency of measurements determines the preferred model.
Statistical evaluation showed a significantly lower initial yield (P<0·001) for the first parity than for second and third parity cows, as expected. Also, a significantly lower peak yield (P<0·001), a later time to peak (P<0·001) and a significantly greater persistency (P<0·001) were obtained for cows in their first parity than for those in later parities, which was in agreement with Stanton et al. (Reference Stanton, Jones, Everett and Kachman1992) and Dijkstra et al. (Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997). Total lactational milk yield was also significantly less (P<0·001) in primiparous than in multiparous cows, as expected.
There were no significant equation×parity effects (Table 6), which shows there was no tendency for an equation to fit a particular parity better than another. This result is in contrast to Dijkstra et al. (Reference Dijkstra, France, Dhanoa, Maas, Hanigan, Rook and Beever1997) who pointed out that, with smooth patterns of lactation data such as primiparous records, the Wood equation will yield a better fit than the Dijkstra equation and the probable reason was that in this case the number of cells and the enzymatic activity per cell, as required parameters of their equation, cannot both be defined uniquely from lactation data owing to the timing and frequency of pre-peak observations being insufficient.
Estimates of equation parameters
The features of lactation calculated from different equations by using the formulae in Table 4 and the parameter estimates are presented in Tables 6 and 7. Parameters a, b, c and d, which determine the scale and shape of each equation, were not compared, because their definition is different in each equation. Regarding the features of the different equations, all equations showed systematic deviation from actual milk yield in accordance with Gipson & Grossman (Reference Gipson and Grossman1989) and Vargas et al. (Reference Vargas, Koops, Herrero and Van Arendonk2000), especially at the beginning of lactation and peak yield. Initial milk production was forced to zero by using the Wood and Morgan equations, which is not strictly acceptable biologically and is a limitation for application of these equations, and was over-predicted by the logistic, Gompertz and Schumacher. However, the Dijkstra equation could predict it accurately. The Wood equation was one of the equations that gave lowest milk production values in early lactation, though it was expected to over-predict at this stage of lactation. All equations except Dijkstra under-predicted peak milk production (observed value taken as the highest average monthly yield and to occur at the mid-month time point), which is in accordance with the results of Scott et al. (Reference Scott, Yandell, Zepeda, Shaver and Smith1996) and Sherchand et al. (Reference Sherchand, McNew, Kellogg and Johnson1995) for the Wood equation. The Dijkstra equation provided the best prediction of maximum yield, which was similar to Val-Arreola et al. (Reference Val-Arreola, Kebreab, Dijkstra and France2004) who also found the Dijkstra equation to fit better to lactation data from Mexican dairy cows. Time to peak yield was over-predicted by the logistic, Gompertz and Schumacher equations and was under-predicted by Dijkstra, but Wood and Morgan could predict it accurately. Despite differences between equations with respect to initial milk yield, time to peak and peak yield, estimates of total milk yield were not significantly different between equations. So, these results showed that each lactation equation (including the Wood and Dijkstra) had some disadvantages in predicting lactation curve features but all equations were able to predict total milk yield accurately, and under the conditions of the present study the tested growth functions showed potential as candidates for the lactation equation.
Table 7. Parameter estimates for the different lactation equations (standard errors are in parentheses)
* Equations 1–6: Wood, logistic, Gompertz, Schumacher, Morgan and Dijkstra, respectively.
† a, b, c and d (all >0) are parameters that define the scale and shape of the lactation curve.
CONCLUSION
A mathematical function allows the lactation curve to be expressed in terms of a set of parameters that has to be estimated. Various functions have been used to study lactation in dairy cattle, with each function having advantages and disadvantages. The growth functions selected were the logistic, Gompertz, Schumacher and Morgan equations and their differentials were compared with the Dijkstra equation and a simpler equation (Wood equation). Based on criteria to measure goodness of fit, the results of this study showed that selected empirical growth functions were able to fit monthly lactation records satisfactorily. Although the Dijkstra equation had more convergence problems in fitting monthly data, the numerical tests were better for that model.
