Reliable data from milk recording systems are important for herd management and genetic improvement in dairy cattle (e.g. Liu et al. Reference Liu, Reents, Reinhardt and Kuwan2000). Because of high costs of milk recording in Slovenia the standard 4-week a.m.-p.m. testing scheme (A4) was replaced with the alternate 4-week a.m.-p.m. (AT4) testing scheme (Klopčič et al. Reference Klopčič, Malovrh, Gorjanc and Kovač2004; Sadar et al. Reference Sadar, Podgoršek, Perpar, Logar, Ivanovič, Jeretina and Prevolnik2005; ICAR, 2006). Since then an overall drop of protein percentage has been detected (Sadar et al. Reference Sadar, Logar, Perpar, Podgoršek, Žabjek and Ivanovič2008) and a necessity for the refinement of models for the estimation of daily records has arisen. A possible reason for the observed drop could be the relatively small data sample used by Klopčič et al. (Reference Klopčič, Malovrh, Gorjanc and Kovač2004). The routine collection of data by regular and supervised milk control provides a larger dataset for the revision of their method.
Several studies considered different methods for the estimation of daily yield (DY) based on partial yield (PY) from either a.m. (AM) or p.m. (PM) milking weights. If the reliability of milking interval (MI) is questionable or the difference between the PYs is negligible, then the simplest method is to use the doubling method. This is a special case of the regression method, where adjustment factors for PY are estimated with a statistical model using the regression of DY on PY and potentially also other effects, such as MI and others. Although optimal in the least square sense, the regression inherently leads to the overestimation of low yields and underestimation of high yields (e.g. Liu et al. Reference Liu, Reents, Reinhardt and Kuwan2000; Klopčič et al. Reference Klopčič, Malovrh, Gorjanc and Kovač2004) which is the prime reason why Galton (Reference Galton1886) used the term regression for this phenomenon. A widely used method for the estimation of DY from PY is the method proposed by DeLorenzo & Wiggans (Reference Delorenzo and Wiggans1986) who derived adjustment factors for several intervals between milkings using the ratio between DY and PY, i.e., DY:AM and DY:PM. Cassandro et al. (Reference Cassandro, Carnier, Gallo, Mantovani, Contiero, Bittante and Jansen1995) compared daily (DY:AM and DY:PM) and partial (AM:DY and PM:DY) yield ratios and showed that adjustment factors for partial yield ratio (PYR) had better properties than adjustment factors for daily yield ratio (DYR). The relationship between DYR is nonlinear: the correlation is close to −1 when the means of ratios are close to 2 and decreases when the means shift towards the infinity or 0. On the other hand, PYR represent the proportions of DY and they always sum to 1, which is also manifested by identical correlations (in absolute value) between them and other variables. This property provides a possibility to double the sample size for the estimation of adjustment factors, i.e., both a.m. and p.m. records can be used to estimate both a.m. and p.m. adjustment factors (DeLorenzo & Wiggans, Reference Delorenzo and Wiggans1986).
The estimate of DY from PY is greatly affected by MI, while the effects of stage of lactation (S), parity (P) and their interactions with MI are usually of minor importance (Everett & Wadell, Reference Everett and Wadell1970a,Reference Everett and Wadellb,Reference Everett and Wadellc; Cassandro et al. Reference Cassandro, Carnier, Gallo, Mantovani, Contiero, Bittante and Jansen1995). Therefore, single milk testing schemes must develop adjustment factors at least for MI.
The purpose of this study was to 1) study the sources of variation affecting the DY, AM:PM, PM:DY, DY:AM and DY:PM and 2) compare the application of AM, PM, AM:DY, PM:DY, DY:AM and DY:PM records for the estimation of DY in the AT4 testing scheme in Slovenia.
Material and Methods
Data
Test-day records of milk yield, fat and protein percentage were collected from the central cattle database GOVEDO, which is hosted at and maintained by the Agricultural Institute of Slovenia (Logar et al. Reference Logar, Podogoršek, Jeretina, Ivanovič, Perpar, Kuipers, Klopčič and Cled2005). Data from regular and supervision milk recordings between March 2004 and April 2008 were used.
For each supervision milk control (the A4 recording method), a corresponding regular milk control carried out by a different controller at the previous milking (the AT4 recording method) was matched. Therefore, approximately 36 h of milk production were covered with three milkings, one from regular (y1) and two from supervision (y2 and y3) milk control (Fig. 1). From these three milkings four records were created as suggested by DeLorenzo & Wiggans (Reference Delorenzo and Wiggans1986) : (y1, y2), (y2, y1), (y2, y3), and (y3, y2). Each record was used to estimate DY from PY or yield ratios (PYR or DYR). Where only supervised control was available, two records were obtained. For each test day MI was calculated. For AT4 milking controls the starting time of the preceding milking was reported by breeder.
Fig. 1. Scheme of regular and supervision milk control.
A total of 518 766 test-day records from 123 503 lactations of 90 719 cows from 26 548 milk test days on 5078 farms were available. The largest portion of records was available for Holstein cows (38%) followed by Simmental cows (29%), Brown Swiss cows (19%), crossbreeds with Simmental (11%) and cows of other breeds (3%). Records with days in milk less than 5 d, MI shorter than 9 h longer than 15 h, and milking three times per day were excluded. In addition, records which did not correspond to the logical controls (ICAR, 2006) were also deleted. After data editing, 493 028 test-day records remained in the data set.
Statistical analyses
Preliminary analyses showed that daily fat and protein content could be estimated with equal accuracy as fat and protein yield. Thus a decision was made to develop models only for yields. Fat and protein content can be subsequently computed using estimated milk, fat and protein DY.
Three models were fitted using different dependent variable (Table 1). The DY model was used to estimate DY directly from PY taking into account the effects of MI, S and P. With this model the estimate of DY is obtained (directly) as the expected value given the PY, MI, S and P. A variant of this model (without the effects of S and P) is currently used in Slovenian milk testing scheme (Klopčič et al. Reference Klopčič, Malovrh, Gorjanc and Kovač2004). In the PYR model, partial yield ratios (AM:DY or PM:DY) were included as a dependent variable, whereas the DYR model used daily yield ratios (DY:AM or DY:PM) as a dependent variable. With these two models, the first step for the estimation of DY is the calculation of expected value for PYR (PYR model) or DYR (DYR model) given the MI, S and P. In the second step, DY estimate is calculated as PY divided (PYR model) or multiplied (DYR model) with the expected value of PYR or DYR, respectively.
Table 1. Models for the estimation of milk, fat and protein daily yield
† 
: daily yield; 
: partial yield from milking i (a.m. or p.m.); bi,0: intercept; bi,1, bi,2: regression coefficients; Si,j: stage (month) of lactation j (j=1,…,10); Pi,k: parity k (k=1 for the first parity and k=2 for the second and subsequent parities); ti,jklm: milking interval; ei,jklm: residual ~ Normal(0, 
)
Comparison between the models was based on the accuracy and bias. Accuracy was defined as the correlation (
) between the true (yA4) and estimated (
) DY. Bias was defined as the average error of estimated DY (
). Additionally, bias was evaluated also at the lower quartile and the upper quartile of the true DY for milk, fat and protein. Values for the lower and upper quartile were 12·9 and 22·4 kg, 0·537 and 0·906 kg and 0·446 and 0·742 kg for milk, fat and protein DY, respectively. Statistical analysis and graphical presentation were performed with SAS (SAS, 2002) and R (R Development Core Team, 2009) program.
Results and Discussion
Data structure and descriptive statistics by parity are presented in Table 2. Two groups of parities were formed since the preliminary analyses showed small difference between the second and the third parity. The larger portion of the data came from the second and later parities (71·3%) which had higher means and larger variability for milk, fat and protein yield. Cows in the first parity had on average 1·6 kg, 0·06 kg and 0·05 kg lower yields for milk, fat and protein, respectively. This was also reflected in PY. However, there were almost no differences between parities in means and variances for yield ratios. The equality of variances for yield ratios is of particular importance, since there is no need to take into account the heterogeneous variances in the model as has been suggested by DeLorenzo & Wiggans (Reference Delorenzo and Wiggans1986) and Liu et al. (Reference Liu, Reents, Reinhardt and Kuwan2000). As expected, the sum of AM:DY and PM:DY was equal to 1.
Table 2. Descriptive statistics by parityFootnote †
† n, number of records; DY: daily yield; AM: a.m. yield; PM: p.m. yield; AM:DY: a.m. to daily yield ratio; PM:DY: p.m. to daily yield ratio; DY:AM: daily to a.m. yield ratio; DY:PM daily to p.m. yield ratio
Milk, fat and protein yield at a.m. milking were higher on average than at p.m. milking by 0·40 kg, 0·01 kg and 0·01 kg, respectively. The main cause of these differences and, consequently, the differences in yield ratios could be attributed to the difference in the length of MI. The nightly MI was 26·7 min longer on average than the daily MI. Klopčič et al. (Reference Klopčič, Malovrh, Gorjanc, Kovač, Osterc and Stekar2001), who did a previous study in Slovenia, reported that the nightly MI was 2 min longer than the daily MI. The reason for the relatively small difference between the daily and nightly MI compared with our study is probably due to a small sample of farms from one region in the study of Klopčič et al. (Reference Klopčič, Malovrh, Gorjanc, Kovač, Osterc and Stekar2001).
There were considerable differences between tested models with regard to the coefficient of determination (Table 3). The DY model explained a high portion of variation in milk, fat and protein DY (0·9 and above), whereas the coefficient of determination was considerably lower in the PYR model and the DYR model, especially for fat yield. Cassandro et al. (Reference Cassandro, Carnier, Gallo, Mantovani, Contiero, Bittante and Jansen1995) performed the analysis with the PYR and the DYR models and reported coefficients of determination approximately half the size of those in our study. This might be the result of data available for both a.m. and p.m. MI in our study, while Cassandro et al. (Reference Cassandro, Carnier, Gallo, Mantovani, Contiero, Bittante and Jansen1995) did not know the exact time of the previous a.m. milking – they calculated the MI between a.m. and p.m. milking as 24 h minus the time before a.m. milking.
Table 3. Analysis of variance for daily yield and yield ratios of milk, fat and protein
*** P<0·001; ** P<0·01; * P<0·05; AM: a.m. yield; PM: p.m. yield; S: stage (month) of lactation; P: parity; MI: milking interval; SP: parity and stage of lactation interaction; yAT4: a.m. or p.m. yield; E: residual
As expected, the most important source of variation for the estimation of yield ratios was MI, whereas PY had the greatest effect on DY followed by MI (Table 3). Variation caused by S, P and their interaction was negligible in comparison with other effects for all traits. Although S, P and their interaction were in most cases significant (P<0·05), only MI was included in equations for the calculation of expected (estimated) values. The same was done for the effect of breed (results not shown), since preliminary analyses showed negligible variation between breeds in addition to other effects in the models.
Regarding the coefficient of determination the DY model seems to have the best fit, but such a comparison is not valid since the models differed in the dependent variable. Therefore the comparison of models was assessed using other criteria (Table 4). The model with the highest accuracy (correlation between the true and the estimated DY); the sd of the estimated DY close to the sd of the true DY; the smallest rmse; and the smallest bias (average error) gives the best fit to the data (e.g. Liu et al. Reference Liu, Reents, Reinhardt and Kuwan2000).
Table 4. Accuracy (correlation between the true (yA4) and estimated () daily yield, 
), sd (
), rmse (
) and bias (average error, 
) of estimated daily yield by model and a.m. or p.m. milking
Accuracies ranged between 0·95 and 0·98 (Table 4), which is similar to the results of Liu et al. (Reference Liu, Reents, Reinhardt and Kuwan2000). There were no differences between the models in the accuracy for a.m. or p.m. milking for all traits. The similarity of sd of the estimated DY and the true DY was higher in the PYR and the DYR models than in the DY model, while the rmse was very similar between all the models. The overall bias was practically zero for the DY and PYR models, while the DYR model overestimated DY. However, the analysis of bias at the lower and upper quartile of DY showed that the DY model overestimated and underestimated DY of all traits. The difference at the lower quartile was 0·33 kg, 0·03 kg and 0·01 kg for milk, fat and protein DY, respectively, using PY from a.m. milking. The difference at the upper quartile was −0·37 kg, −0·04 kg and −0·01 kg for milk, fat and protein DY, respectively, again using PY from a.m. milking. Similar biases were observed also with the use of PY from p.m. milking. Overestimation and underestimation with regard to low and high DY, respectively is commonly observed for the regression method (e.g. Liu et al. Reference Liu, Reents, Reinhardt and Kuwan2000; Klopčič et al. Reference Klopčič, Malovrh, Gorjanc and Kovač2004). With this method the mentioned bias cannot be removed with the addition of any effect to the model or even estimating separate models for each combination of effects in the model as has been advocated by Liu et al. (Reference Liu, Reents, Reinhardt and Kuwan2000).
The amount of bias was further analysed for DY of milk for each class of S (Fig. 2) and P (Fig. 3) for a.m. and p.m. milking separately. The bias by the DY model was clearly changing from negative in the first stages of lactation through 0 in the middle of lactation to positive at the last stages of lactation (Fig. 2). This was observed for both a.m. and p.m. milking. With the PYR model the bias was slightly negative in the early stages of lactation for a.m. milking and the opposite (slightly positive) for p.m. milking. In the later stages of lactation the bias for the PYR model was close to 0 for both a.m. and p.m. milking. As observed in the analysis of overall bias (Table 4) the DYR model generally overestimated DY from both a.m. and p.m. milking. The change of bias over the lactation for the PYR model was similar to that of the DYR model. The analysis of bias by parity (Fig. 3) showed slightly smaller values than by the stage of lactation. Generally, the DY was more often overestimated in the first than in the later parities when using a.m. milking and vice versa for p.m. milking. The highest bias by parity was observed for the DYR model, while the bias was lower for the other two models.
Fig. 2. Bias (average error) of estimated daily milk yield by stage (month) of lactation and model from a) a.m. and b) p.m. milking.
Fig. 3. Bias (average error) of estimated daily milk yield by parity and model from a) a.m. and b) p.m. milking.
Analysis of bias according to the amount of the true milk DY (Fig. 4) confirmed that the DY model systematically overestimated low (positive bias) DY and underestimated high (negative bias) DY, which is in accordance with the findings of Liu et al. (Reference Liu, Reents, Reinhardt and Kuwan2000) and Klopčič et al. (Reference Klopčič, Malovrh, Gorjanc and Kovač2004). Bias reached up to 1 kg for cows yielding more than 40 kg of milk/d. The lowest bias over the whole interval of the true milk DY was obtained with the PYR model, while the bias constantly increased with the value of the true milk DY for the DYR model. The oscillation of bias above the 40 kg/d milk DY is probably due to a small number of records.
Fig. 4. Bias (average error) of estimated daily milk yield by daily milk yield and model from a) a.m. and b) p.m. milking.
The results show that the PYR model is the most appropriate for the estimation of DY from yield at a.m. or p.m. milking. The appropriateness of this model stems from the fact that the modelled variable is not the actual daily yield but the partial yield ratio – the proportion of partial to daily yield. This variable is biologically more related to the problem of the estimation of DY than the actual yields. A possible objection to this model could be the fact that this variable is continuous but bounded between 0 and 1, for which a model with beta distribution would be more appropriate (e.g. Smithson & Verkuilen, Reference Smithson and Verkuilen2006). However, the ratio of partial yield to daily yield will be most of the time around 0·5 unless there is an error in the data. Occasional deviations, due to variation in MI, environmental effects, cows in heat or similar effects, constitute a symmetric distribution that can be well approximated with the normal distribution, which is the implied distribution of the least squares method used in our study to estimate model parameters (Table 1).
In conclusion, the most important sources of variation for DY were PY and MI. The factor that described the highest variability of yield ratios was MI. Other factors (S, P and their interaction) accounted for a smaller amount of variability in DY or yield ratios. The problem of overestimating low DY and underestimating high DY milk records from just a.m. or p.m. milking arises if DY is estimated directly from partial yields (our DY model). When DY is calculated via division of PY with the estimated partial yield ratios (AM:DY and PM:DY – our PYR model) the problem of overestimating low DY and underestimating high DY vanishes. Daily yield ratios (DY:AM and DY:PM – our DYR model) do not have such properties. To our knowledge these findings were not recognized in previous studies. Based on our results, the PYR model is the most appropriate for the estimation of daily yield from single a.m. or p.m. milking. The advantage of this model is due to modelling the variation in partial yield ratios that are biologically more stable for predicting daily yield than the actual yield values.
The authors thank Boris Ivanovič, Janez Jeretina and Marija Sadar for assistance in processing the data for this study. The comments of three anonymous reviewers improved the initial version of manuscript.
