Journal of Dairy Science Vol. 85 No. 11 3107-3114
© 2002 by American Dairy Science Association ®
Test Day and Lactation Yield Predictions in Italian Simmental Cows by ARMA Methods
N. P. P. Macciotta*,
D. Vicario
,
G. Pulina* and
A. Cappio-Borlino*
* Dipartimento di Scienze Zootecniche, Università di Sassari, Via De Nicola 9, 07100 Sassari, Italia
Italian Association of Simmental Breeders, Via Nievo 19, 33100 Udine, Italia
Corresponding author:
Nicolò Pietro Paolo Macciotta; e-mail:
macciott{at}ssmain.uniss.it.
 |
ABSTRACT
|
|---|
Autoregressive Moving Average (ARMA) models, originally developed in the context of time series analysis, were used to predict Test Day (TD) yields of milk production traits in dairy cows. ARMA models are able to take into account both the average lactation curve of homogeneous groups of animals and the residual individual variability that may be explained in terms of probability models, such as Autoregressive (AR) and Moving Average (MA) processes. Milk, fat, and protein yields of 6000 Italian Simmental cows with 8 TD records per lactation were analyzed. Data were grouped according to parity (1st, 2nd, and 3rd calving) and fitted to a Box-Jenkins ARMA model in order to predict TD yields in five situations of incomplete lactations. Reasonable accuracies have been obtained for a limited horizon of prediction: average correlations among actual and predicted data were 0.85, 0.72, and 0.80 for milk, fat and protein yields when the first predicted TD was one step ahead (on average 42d) of the last actual record available. Cumulative 305-d yields were calculated using all actual (actual yields) or actual plus forecasted (estimated yields) daily yields. Accuracy of lactation predictions was remarkable even when only a few actual TD records were available, with values of 0.88 for milk and protein and 0.84 for fat for the correlations between actual and estimated yields when 6 out of 8 TD records were predicted. Accuracy rapidly increases with the number of actual TD available: correlations were about 0.96 for milk and protein and 0.93 for fat when 4 out of 8 TD records were predicted. In comparison with other prediction methods, ARMA models are very simple and can be easily implemented in data recording software, even at the farm level.
Key Words: ARMA models • milk production traits • test day
Abbreviation key: ACF = Autocorrelation function, AR = autoregressive, ARMA = Autoregressive Moving Average, ATD = Actual Test Day, AY = Actual yields, EDS = Estimation data set, MSEP = mean square error of prediction, PTD = Predicted Test Day, PY = Predicted yields, TD = Test Day, VDS = Validation data set
|
INTRODUCTION
|
|---|
The knowledge of a large number of test day (TD) records of milk production traits at regular intervals of time lactation represents an advantage for the genetic analysis by test day models, for control and management purposes, and for the estimation of cumulative 305-d yields, thus avoiding the use of extension factors. There are several situations in which few TD records are available, however, as with lactations in progress, or when some TD records are missing. Moreover, as TD recording is one of the most important expenses in genetic programs, especially for milk composition traits, there is a great deal of interest in reducing the number of TD recorded during lactation.
Several methods have been suggested to predict TD yields from a limited number of actual records available. Most of them are based on parametric lactation curves of homogeneous groups of cows combined with information on individual variation. For example, Perochon et al. (1996) developed models of the lactation curve which account for variations between animals in order to provide basis for forecasting individual yields. Jones (1997) proposed an empirical Bayesian method in which milk yields from a lactation in progress are combined with prior information gathered from herd mates. Schaeffer and Jamrozik (1996), in order to reduce the number of TD for fat and protein composition to be recorded, suggested a multiple-trait procedure that incorporated information about average lactation curves and covariances between TD yields for milk, fat and protein. Finally, Pool and Meuwissen (1999) developed phenotypic TD models with Legendre polynomials, able to interpolate and extrapolate missing records. These sophisticated approaches are able to give predictions with reasonable accuracy but at the expense of a great computational demand and mathematical effort.
An alternative empirical procedure for predicting TD records of milk production can be traits based on Time Series analysis. Some authors (Deluyker et al., 1990; Lark et al., 1999; Macciotta et al., 2000) showed that Autoregressive Integrated Moving Average (ARMA) models, developed in the Time Series context, can be easily adapted to describe sequences of milk yield records. Thus, an appropriate time series model could be used to forecast TD yields for each individual cow as soon as a few records become available for the current lactation.
The aim of this work is to check the suitability of ARMA methods to model and predict TD milk production data in cattle by using records of an official dairy recording scheme and simulating situations of incomplete lactations.
|
MATERIAL AND METHODS
|
|---|
Theoretical Basis of Time Series Analysis
A time series is a succession of real values of a variable Y (Y1, Y2,...,Yn) ordered by an index that refers to an interval of time. However, mathematical methods specifically developed for analyzing and modeling time series can be applied to any ordered succession of values regardless of the meaning of the ordering index (Hamilton, 1994). For example, if there are N lactations of a group of cows, each with n daily records, a succession of N x n TDs (Y1, Y2,...Yk, ...,YN x n) can be constructed where the index value (k) refers to a specific TD of a specific lactation (e.g., assuming that N = 100 and n = 8, there is a succession of 800 TDs in which at k = 94 corresponds the 6th TD record of the 12th lactation). The variability of this series can be decomposed into two main components. The first is a continuous component which occurs at intervals of the index (lags) multiples of n and that corresponds to the average lactation curve of the group of analyzed animals. The second is a stochastic component that expresses the variance of TD residuals around the mean lactation curve. The so-called ARMA seasonal model, originally suggested by Box and Jenkins (1970), can be adapted for describing a periodic succession as the one cited above; these models are able to account for both relationships among TD records within lactation (the nonseasonal part of the model) and among correspondent TD records of different lactations of the succession (the seasonal part of the model) by means of Autoregressive (AR) and Moving Average (MA) processes. In an AR process of order p, the value of Yk is modelled as a function of p previous values of Y, i.e., Yk-1, Yk-2,..., Yk-p. In a MA process of order q, Yk is regarded as the weighted sum of the q values of a random variable Z (Zk-1, Zk-2,..., Zk-q), where Zk measures the innovation of the series, i.e., the difference between the actual k-th value of Y and the corresponding value estimated by the model: Zk = Yk – 
k (Lark et al., 1999).
Data
Three data sets of 2000 cows each, grouped according to parity (first, second, and third calving) were extracted from the historical archive of the Italian Simmental cattle, a dual-purpose breed. Each cow had eight TD records for three milk production traits (milk, fat, and protein yield), with an average interval of 42 d between adjacent measurements. In each archive, lactations of different cows were randomly arranged in a vector in sequence (N = 1,...,2000), and within each cow TD measures were ordered according to their distance from calving (n = 1,...,8). Therefore the data set consisted of 16,000 TD records ordered in a series. Each data set was then split into: an estimation data set (EDS), consisting of the first 1000 lactations, used for estimating parameters; and a validation data set (VDS), consisting of the next 1000 lactations, which was used for making predictions by using parameter estimates with the EDS records.
Model Identification and Forecasting
The identification of a suitable ARMA model, i.e., the definition of the order of the AR and of MA parameters, both for the seasonal and nonseasonal parts, was performed by examining the pattern of the autocorrelation function (ACF) of the EDS data (Box and Jenkins, 1970; Diggle, 1990). According to the results of the diagnostic step, ARMA model parameters were then estimated by a maximum likelihood procedure (SAS, 1996) for all three traits within each parity. Goodness of fit was assessed by testing the statistical significance of parameters, by examining the ACF of residuals, and by calculating the percentage of the total variance explained by the model in analogy with the determination coefficient of multiple regression was determined (Piccolo, 1990; Chatfield, 1999).
Five plans of incomplete lactations were simulated in the VDS for each parity class and for each trait, with different TD records available: plan 1 on tests 1 and 2; plan 2 on tests 1, 2, and 3; plan 3 on tests 1, 2, 3, and 4; plan 4 on tests 1, 2, 3, 4, and 5; plan 5 on tests 1, 2, 3, 4, 5, and 6. ARMA models fitted to the corresponding EDS were then used to predict the missing TD. Consequently, six subsets of data were available for each cow of the EDS: all eight actual TDs, two actual and six predicted TDs, three actual and five predicted TDs, four actual and four predicted TDs, five actual and three predicted TDs, and six actual and two predicted TDs. Total milk, fat, and protein yields were calculated (CICPE, 1992) resulting in an actual yield (AY), calculated with all eight actual TDs, and five predicted yields (PY) obtained combining actual and predicted TDs. Goodness of TD record prediction was evaluated by calculating the mean square error of predictions (MSEP) and the correlations among actual (ATD) and predicted (PTD) TD records. Goodness of predictions for total yields was assessed by means of the following criteria: correlations between AY and PY; standard deviations of AY and PY; MSEP; correlation among AY and residuals; the prediction bias (mean AY-PY); and the ratio between the standard deviation of differences between AY and PY and the mean value of AY.
|
RESULTS AND DISCUSSION
|
|---|
Results of the diagnostic stage are reported for second parity cows and for milk yield only; similar results have also been obtained for first and third parities and for fat and protein yields. The ACF pattern of the 1000 cows of the EDS (Figure 1
) highlights the existence of a periodic component, as evidenced by values significantly different from zero at lag multiples of eight. Such a component can be identified with the amount of variability explained by the average lactation curve and modelled using a time series seasonal model with a season of eight lags. Furthermore, the rapid decrease of autocorrelation values for lags greater than one and the sinusoidal pattern of ACF are evidence of AR and/or MA processes of a low order (one or two). Actually, these processes are peculiar of longitudinal data because measurements close in time are more correlated than those far apart (Littel et al., 1998). In the specific case of lactation, they can be regarded as expression of the individual variation around the mean lactation curve due to environmental short-term effects such as dietary quality, intake, weather, minor injuries, and estrus (Ali and Schaeffer, 1987; Carvalheira et al., 1998; Wade et al., 1993) as well as genetic variations.
View larger version (12K):
[in this window]
[in a new window]
|
Figure 1. Autocorrelation function of TD milk yields for second parity cows. Dotted lines are plotted at zero plus and minus twice the standard errors for each coefficient.
|
|
On the basis of the aforementioned considerations, the most suitable model for the successions under study is the following seasonal ARMA model with a season of eight lags and AR and MA processes of order 1:
 | [1] |
where: 
1 and ß1 are the AR and MA nonseasonal parameters, respectively; 2 and ß2 are the AR and MA seasonal parameters, respectively; and Z is the innovation of the series, with E(Z) = 0. Given this structure, in the prediction step model [1] takes account of: i) the relationship among adjacent TDs and between TsD that are eight lags apart, i.e., among the same TD of different lactations; ii) the prediction error made in the estimation of the immediately previous TD and of the TD which is eight lags apart.
ARMA parameters for the three traits within each parity class are reported in Table 1. Goodness of fit is evidenced by the statistical significance of the parameters and further confirmed by the absence of autocorrelations among residuals. Seasonal parameters are largest, confirming the great importance of the average lactation curve of the group of animals; however, also values of AR nonseasonal parameters, i.e., the ones that account for the linkage among TD measurements within lactation, underlines the relevance of this component in the temporal pattern of milk production. The proportion of the original variance explained by the ARMA model was highest for milk yield (0.62, 0.71, and 0.72 for first, second, and third parity, respectively), lowest for fat (0.37, 0.49, and 0.45), and intermediate for protein (0.53, 0.63, and 0.65). In general, results confirm the higher suitability of the ARMA model for milk and protein yield in comparison with fat yield.
The mean square error of prediction for TD yields tends to decrease as the number of ATD increases, for all the three traits and all the three parities (Figures 2a, b, and c). Values of MSEP for milk production are of the same order or slightly lower than those obtained in previous studies with other methodologies (Perochon et al., 1996; Jones, 1997; Pool and Meuwissen, 1999). These results highlight the great suitability of time series models to detect short-term variations in the over time evolution of milk production. On the other hand, it must be noticed that the prediction error is not homogeneously distributed across all lactations, but is mainly localized in cases of very atypical patterns, such as those characterized by a sudden increase of production in the declining phase or by an exceptional persistency. A few examples will help to explain how the model works. Figure 3a shows a selected case of very good fit of the individual lactation curve of milk yield for a second parity cow: the continuous line is the actual lactation curve, whereas the dotted line is the predicted curve with the first two TD records available and last six to be predicted. A selected case of bad fit is reported in Figure 3b; the actual curve shows an atypical pattern, already observed in dairy cattle (Shanks et al., 1981), and that can be quite common, especially in herds characterized by high production levels. Such a great sensitivity may be regarded as a weak point of the model: however, it must be remembered that data sets were randomly built, considering only parity as a classification factor. It might therefore have happened that lactations included in the EDS were obtained in very different conditions (genetic merit of animals, management levels of the herd, season of calving and of production, geographical region of location of the herd) from those of the corresponding VDS. In view of a technical implementation, predictions can be improved by estimating ARMA parameters with cows as homogeneous as possible with those to be predicted. Furthermore, the high sensitivity of the model to irregular lactation patterns could be usefully exploited in control systems, as, for example, alert mechanisms in automated milk recording systems in order to detect animals affected by diseases.
View larger version (16K):
[in this window]
[in a new window]
|
Figure 2a. Mean square error of predictions of test day milk yields in the three parities (1st = black bars; 2nd = grey bars, and 3rd = white bars for the five plans of incomplete lactations (Plan 1: tests available 1, 2; plan 2: 1, 2, 3; plan 3: 1, 2, 3, 4; plan 4: 1, 2, 3, 4, 5; plan 5: 1, 2, 3, 4, 5, 6).
Figure 2b. Mean square error of predictions of test day fat yields in the three parities (1st = black bars; 2nd = grey bars, and 3rd = white bars for the five plans of incomplete lactations (Plan 1: tests available 1, 2; plan 2: 1, 2, 3; plan 3: 1, 2, 3, 4; plan 4: 1, 2, 3, 4, 5; plan 5: 1, 2, 3, 4, 5, 6).
Figure 2c. Mean square error of predictions of test day protein yields in the three parities (1st = black bars; 2 = grey bars, and 3rd = white bars for the five plans of incomplete lactations (Plan 1: tests available 1, 2; plan 2: 1, 2, 3; plan 3: 1, 2, 3, 4; plan 4: 1, 2, 3, 4, 5; plan 5: 1, 2, 3, 4, 5, 6), and errors for each coefficient.
|
|
View larger version (11K):
[in this window]
[in a new window]
|
Figure 3a. Actual (— = all eight actual test day records) and estimated (–––– = 2 actual and 6 predicted) individual lactation curve for milk yield of a second parity cow.
Figure 3b. Actual (—) and estimated (––––) individual lactation curve for milk yield of a second parity cow.
|
|
Correlations among actual and predicted TD yields are on average 0.85, 0.72, and 0.80 for the first PTD in each plan for milk, fat, and protein yields, respectively (Table 2). Values tend to decrease, even if an average value of 0.70 can be found for milk yield for PTD that are three lags apart from the last ATD.
Results of lactation predictions are reported in Tables 3(a, b, and c). Pearson correlations among actual and predicted milk yields (Table 3a) are remarkable in the last three plans, and are of the same order of magnitude as those reported in previous studies on dairy cattle (Schaeffer and Jamrozik, 1996); however, the correlation of 0.88 obtained when only two TDs out of eight are known must not be undervalued. Standard deviations of actual and estimated yields are quite similar, at least in the last two or three plans.
The average bias is negligible, and it underlines a slight tendency of the model to overestimate lactation yields, in the case of first calving, and to underestimate them for the other two parity classes. However, the existence of a bias is evidenced by values of correlations among actual yields and residuals that indicate a tendency of residuals to increase as AY increases. The situation in which the first four TD records are available represent a lower bound of good prediction, as can be inferred from the values of the mean square error of prediction. Finally, the last index, obtained as 100 times the ratio of the standard deviation of residuals to AY, was below eight, except for the first plans, which is considered a maximum threshold for a good prediction (Wilmink, 1987).
Results concerning fat and protein yields are similar, although they are a little worse for the former, thus confirming the lower predictability of this trait, probably due to its greater sensitivity to environmental variation.
|
CONCLUSIONS
|
|---|
ARMA models gave interesting results in terms of accuracy of prediction both for TD and lactation yields, comparable to those obtained with more sophisticated methodologies. TD forecasting was efficient for a limited horizon (i.e., when PTD were at not more than two or three lags far apart from the last ATD), whereas lactation yield extrapolation became interesting when at least three TD records were available. ARMA performances could be further improved by a more detailed choice of animals to be used for parameter estimation. With this aim, an implementation on a herd level with completed lactations as the data set of reference to estimate parameters seems to be an advisable option.
In this work, ARMA methodology has been used in a univariate approach, but time series models can also be employed in a multivariate system: therefore, the predictions of fat TD yield could be improved by considering also the correlation structure with milk and protein yield.
These methods were very easy to implement; we used commercial software on a PC with a Pentium I 133 MHz microprocessor, 64 Mb of RAM, and the estimation step took about 2 min, whereas the prediction step took 3 to 4 min.
|
ACKNOWLEDGEMENTS
|
|---|
Research was funded by the Italian Ministry of University and Research (grant 40%). The authors wish to acknowledge C. Zanolla for her technical support.
Received for publication January 2, 2002.
Accepted for publication May 8, 2002.
|
REFERENCES
|
|---|
Ali, T. E., and L. R. Schaeffer. 1987. Accounting for covariances among test day milk yields in dairy cows. Can. J. Anim. Sci. 67:637–644.
Box, G. E. P., and G. M. Jenkins. 1970. Time series analysis: Forecasting and control. Holden Day, San Francisco, CA.
Carvalheira, J. G. V., R. W. Blake, E. J. Pollak, R. L. Quaas, and C. V. Duran-Castro. 1998. Application of an autoregressive process to estimate genetic parameters and breeding values for daily milk yield in a tropical herd of lucerna cattle and in United States Holstein herds. J. Dairy Sci. 81:2738–2751.[Abstract]
Chatfield, C. 1999. The analysis of time series. An introduction. Chapman 
Hall, London, UK.
CICPE Comitee International par les Controles des performances en Elevage. 1992. Reglement international pour le controle laitiere ovin. Rome, Italy.
Deluyker, H. A., R. H. Shumway, W. E. Wecker, A. S. Azari, and L. D. Weaver. 1990. Modelling daily milk yield in holstein cows using time series analysis. J. Dairy Sci. 73:539–548.[Abstract]
Diggle, P. J. 1990. Time series, a biostatistical introduction. Oxford statistical science series, 5. Oxford University Press.
Hamilton, J. D. 1994. Time series analysis. Princeton University Press, Princeton NY.
Jones, T. 1997. Empirical Bayes prediction of 305-day milk production. J. Dairy Sci. 80:1060–1075.[Abstract]
Lark, R. M., B. L. Nielsen, and T. T. Mottram. 1999. A time series model of daily milk yields and its possible use for detection of a disease (ketosis). Anim. Sci. 69:573–582.
Littel, R. C., P. R. Henry, and C. B. Ammermann. 1998. Statistical analysis of repeated measures data using SAS procedures. J. Anim. Sci. 76:1216–1231.[Abstract/Free Full Text]
Perochon, L., J. B. Coulon, and F. Lescourret. 1996. Modelling lactation curves of dairy cows with emphasis on individual varability. Anim. Sci. 63:189–220.
Piccolo, D. 1990. Introduzione all’analisi delle serie storiche (Introduction to time series analysis). La Nuova Italia Scientifica eds., Roma.
Pool, M. H., and T. H. E. Meuwissen. 1999. Prediction of daily milk yields from a limited number of test days using test day models. J. Dairy Sci. 82:1555–1564.[Abstract]
SAS/ETS Software. User’s Guide. Version 6, 2nd edition. 1993. SAS Institute Inc., Cary, NC.
Schaeffer, L. R., and J. Jamrozik. 1996. Multiple-trait prediction of lactation yields for dairy cows. J. Dairy Sci. 79:2044–2055.[Abstract]
Shanks, R., P. J. Berger, and A. E. Freeman. 1981. Genetic aspects of lactation curves. J. Dairy Sci. 64:1852–1860.
Wade, K. M., R. L. Quaas, and L. D. Van Vleck. 1993. Estimation of the parameters involved in a first order autoregressive process for contemporary groups. J. Dairy Sci. 76:3033–3040.[Abstract/Free Full Text]
Wilmink, J. B. M. 1987. Comparison of different methods of predicting 305-day milk yield using means calculated from within-herd lactation curves. Liv. Prod. Sci. 17:1–17.
This article has been cited by other articles:
|
|
|
|
 
J. Vasconcelos, A. Martins, M. F. Petim-Batista, J. Colaco, R. W. Blake, and J. Carvalheira
Prediction of Daily and Lactation Yields of Milk, Fat, and Protein Using an Autoregressive Repeatability Test Day Model
J Dairy Sci,
August 1, 2004;
87(8):
2591 - 2598.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
