Information on future milk yield performance (even if only approximate) provides premises for making rational decisions and stimulates actions aimed at the realization of favorable prediction (e.g. appropriate animal care). Daily or weekly milk yield prediction allows early identification of low-producing animals, which can subsequently be culled. Perturbations in milk production may indicate health [1,2] or reproductive [3] problems. Farmers can predict milk production level at successive lactation stages based on lactation curve shape [4]. Production costs, expected profits and the economic value of a cow or cows can be estimated based on current milk prices. They are also used for making decisions about subsequent artificial inseminations of such cows, their further utilization and herd replacement. Consequently, the policy of advance planning can be implemented.

Therefore, it seems understandable that the possibility of predicting lactation curve shape of cows is a highly beneficial endeavor, especially for primiparous cows, whose real production performance is unknown despite the previous breeding value estimation. The recording of milk yield in young cows from the previous production periods on a farm is the basis for generating predictions that may be compared with the real values achieved by the next generations of primiparae. Thus, it enables the verification of selection efficiency over successive years.

One of the most popular methods of milk yield prediction is a regression model, especially Wood’s model, but also the so-called artificial intelligence methods (including artificial neural networks) or time series prediction models. Wood’s model, as a mathematical equation describing the relationship between milk yield and lactation duration is quite frequently applied to milk yield prediction in cows [5–7] and often serves as a good reference model for comparison with other prediction models. This model includes parameters associated with milk yield at peak lactation, the stage of milk yield increase and its decrease after the lactation peak.

Another method of lactational milk yield prediction is the seasonal auto-regressive integrated moving average (SARIMA) model, which is based on the auto-regressive integrated moving average (ARIMA) model, whose theoretical basis was developed by Box and Jenkins in 1976 [8]. Assuming that lactation is characterized by certain homogeneity (i.e. excluding the local level and trend, any part of the lactational yield series is very similar in its shape to every other part), the difference among them is a stationary, mixed autoregressive and moving-average process. The SARIMA models are quite frequently used for different prediction tasks, including agricultural ones such as annual milk production prediction in India [9,10], lactation yield prediction in sheep based on test-day records [11], the prediction of weekly and lactational milk yield in Egyptian cows [12] or weekly milk yield in Spanish goats [13].

The next prediction method based on machine-learning, i.e. the recognition of some characteristic patterns of a given phenomenon through the learning process, is a nonlinear autoregressive exogenous (NARX) artificial neural network applied in various domains, including animal breeding, e.g. milk yield prediction in cows during lactation [14–16].

Consequently, the aim of the present study was to verify the usefulness of unconventional SARIMA and NARX methods for the prediction of lactational milk yield in primiparous cows kept on a medium-sized farm and to compare their predictive performance with the reference Wood’s model.

Since the present study involved only the analysis of production records routinely collected on the farms, the approval from the Local Ethics Committee on Animal Research was not necessary. The dataset comprised monthly test-day records from primiparous Polish Holstein-Friesian Black-and-White cows kept on one of the dairy farms located in the West Pomeranian Province. The A4 method was used for milk recording. In this method, the amount of milk is determined at least 11 times per year (the so-called test-days) from each cow during 24 hours and the total milk sample is collected, which is also used for milk composition analysis [17]. The data were collected between 2009 and 2017 from a total of 965 primiparous cows, which accounted for approximately 29% to 35% of the annual number of cows on the farm. Primiparae were assigned to age-season groups with two age categories (since the primiparous cows under local conditions have slightly different lactation curve shapes; those calving before 26 months of age have lower lactational yields than the animals beginning production at an older age [18]).

The following age groups were distinguished: younger primiparae, calving at the age of 20 to 26 months, older primiparae, calving at the age of 27 to 32 months. In addition, four calving seasons were defined: i) primiparae calving in January, February and March, ii) primiparae calving in April, May and June, iii) primiparae calving in July, August and September, iv) primiparae calving in October, November and December.

A total of eight age-season groups were created in this way (whose sample-size allowed for estimating model parameters and training). Lactation was divided into ten lactation stages of approximately 30 days each (except for the first lactation stage [from 5 to 30 days in milk] and the tenth lactation stage [from 270 to 305 days in milk]). Only the milk yields of more than 5 kg were included in the analysis. The milk yields recorded during the successive lactation stages were averaged over each age-season group. For instance, the average milk yields (kg/d) from 18 younger primiparae in their first lactation stage that were recorded in the first season of 2016 are presented in Table 1. These milk yields were averaged giving 28.49 kg of milk for the group of younger primiparae in the first calving season of 2016. The same method was applied to calculate the milk yields for subsequent lactation stages and obtain the values for the whole averaged lactation in the primiparae from a given age-season group.

The period of ten lactation stages was treated as one production cycle with a seasonal lag of 10, which did not always correspond to the production year of a primiparous cow. If the cow started production in the third or fourth season, her subsequent milk yields were recorded in the next year. However, to avoid confusion, the year-wise nomenclature was adopted by treating the years 2009 to 2010 as the year 2009, the years 2010 to 2011 as the year 2010 etc. (Supplementary Table S1). The numbers of younger primiparae in the analyzed years and seasons ranged from 10 to 18 and those of older primiparae ranged from 10 to 20 (Supplementary Table S1).

In the study period, the animals were housed in a free-stall barn under similar production conditions over the studied years and fed a standard diet, consisting of concentrate, green forage or silage. They also had access to outside runs. The average milk yield of the primiparae during the analyzed period amounted to approximately 9,000 kg over the standard lactation (Supplementary Table S2); however younger cows had lower yield. The records from the years 2009 to 2015 were used for the identification and parameter estimation of Wood’s and SARIMA models and for the development of the NARX artificial neural networks (training and validation data sets). The data from 2016 were utilized for the comparison of the real values with the predictions generated by the SARIMA and NARX models and those produced by Wood’s reference model (testing data set).

The application of the SARIMA models or the NARX artificial neural networks may be an interesting alternative to other methods of average milk yield estimation at successive lactation stages, such as Wood’s model [11,30]. In the light of our study and the available literature [14–16,31], the usefulness of different statistical methods for solving problems related to milk yield prediction could be confirmed. The fitting of the SARIMA models required the determination of the p, d, q, and P, D, Q parameters based on the autocorrelation and PACFs, parameter estimation and the goodness-of-fit verification. The correlogram analysis itself is not an ideal tool for the accurate model parameter selection, since the correlograms do not always unequivocally indicate the adequate choice of the p, q, and P, Q parameters. In the selection of the SARIMA model for specific prediction, AIC and AIC_c can be used. It is worth noting that, in most cases, the models with the low AIC values were also characterized by the low MAE, RMSE, and MAPE values. As stated by Box et al [21], the fitting of the processes with the increasing values of the p and q or P and Q parameters until the moment, when the residuals of such a fitted process do not show autocorrelation, can be limited to the value of p+q/(P+Q)<3. The investigated SARIMA models did not differ considerably in the number of parameters and the AIC and AIC_c values were mainly associated with the above-mentioned errors. It seems that the ultimate selection of the best predicting SARIMA model can be based on RMSE, MAE, and MAPE. The AIC measure, including the number of model parameters, is not more accurate, since the models usually did not differ between each other by more than one or two parameters. In the case of artificial neural networks, the number of parameters was incomparably higher than that for the SARIMA models. Therefore, the comparison between SARIMA and NARX using AIC and AIC_c was not performed, similarly as for Wood’s model.

It can be stated that the predictions presented in our study reflected the course of real milk yields during the distinguished lactation stages quite well. Adopting the 30-day lactation stages, some deviations expressed in kilograms ranged from 0.37 kg to 210 kg. The average deviations for the whole lactation ranged between 22.5 kg and 802 kg of milk. The performed analysis showed that the average predictions at successive lactation stages with the lowest deviations from the real milk yields in younger cows were generated by the SARIMA models, especially during the first lactation periods. At later stages, Wood’s and NARX models were more accurate with an increasing prediction horizon. In the case of older primiparae, this relationship became apparent in season 4, but it was not so evident for the remaining seasons. Wood’s models were more accurate, especially in a quite problematic season 3, similarly as the NARX models. This observation was further confirmed by the low values of the Theil coefficients for five SARIMA models and three Wood’s models. The NARX networks had slightly higher Theil coefficients. In the case of SARIMA (four models), NARX (five models) and Wood’s models (five models), the value of the Theil coefficient was mainly affected by the error (I²_E) resulting from the insufficient convergence between the directions of prediction changes and the changes in the predicted variable. This indicates an insufficient prediction of turning points, which may lead to the prediction overestimation in some lactation periods (seasons 2 and 4 for younger cows and season 4 for older cows) and underestimation in others (season 1 and 3 for older cows) (Figures 2, 3). This finding corresponds to some extent with the results obtained by Murphy et al [15], who investigated lactational milk yield prediction at different time horizons and observed greater MAE and RMSE values for the NARX networks in the first prediction horizon (within 30 days in milk) but these errors decreased in the subsequent prediction periods. On the other hand, the static artificial neural networks were characterized by a general increase in the prediction accuracy with decreasing horizon (a reduction in RMSE from 12.0% to 10.7% with the prediction horizon changing from 305 to 10 days), although in the same prediction horizon range, the underestimation of lactation peak yield increased from 4.8% to 7.4%. However, when analyzing these predictions at 10-day intervals, the errors for NARX and the rest of the models were characterized by larger fluctuations – the neural networks had lower errors in some periods, and the remaining models in others. The predictions generated by the SARIMA, NARX, and Wood’s models in our study were more flexible, which was confirmed by the lowest contribution of the I²_B term to the overall Theil coefficient. Moreover, a very low contribution of the I²_O term (especially for the SARIMA models) indicated a low prediction bias. However, this coefficient predominated among the components of the Theil coefficient for three NARX and three Wood’s models, confirming a higher prediction bias for them.

A different trend in the error values was found in the prediction of monthly total dissolved solids in the Rio Grande [32], where the authors observed slightly lower RMSE and MAE values for the neural networks compared with the ARIMA model at the one-month prediction horizon, whereas this situation changed in the second and third month, when smaller errors were made by the ARIMA model. A similar relationship was found by Zhang et al [31], who recorded an increase in RMSE for the NARX model with lengthening prediction horizon. However, it should be remembered that the results from these studies are not directly comparable with ours and only some general trends can be observed.

The correlation coefficients between the predicted and real milk yields in our study were slightly higher than those obtained by Abudu et al [32], who recorded the values ranging from 0.78 to 0.89; however, these correlations were almost always greater for the neural networks compared with ARIMA, similarly as in our study. In the work by Zhang et al [31] on the use of an automatic system of model configuration and optimization for milk production prediction, the correlation coefficients exceeded 0.84 for all the nine developed models (multinomial, adaptive multinomial, Legendre’s multinomial, cubic splines, log-squared, multiple linear regression, static artificial neural networks, surface fitting and NARX), which confirmed their good overall quality. At a longer prediction horizon (amounting to 365 days), the correlation coefficients for these models ranged from 0.85 to 0.98 (r = 0.98 for the best surface fitting model and r = 0.85 for the worst multiple linear regression). The accuracy of the NARX model increased with a decreasing prediction horizon (the correlation rose from 0.96 to 0.98 with the horizon changing from 365 to 10 days). Also, in the study by Murphy et al [15] on the use of the NARX and static artificial neural networks, among others, for milk production prediction, the correlation coefficients generally slightly increased with a decreasing prediction horizon (from 305 to 10 days) and they were slightly higher for NARX compared with the static artificial neural networks each time (0.97 vs 0.94 to 305 days, 0.97 vs 0.94 to 50 days, 0.98 vs 0.95 to 30 days, and 0.98 vs 0.95 to 10 days for the NARX model and the static artificial neural networks, respectively). However, such a trend was not observed for the last investigated model, i.e. multiple linear regression, for which the values of the correlation coefficient were almost invariant.

We obtained exceptionally worse predictions (compared with other seasons) in the third season for older cows, which is also depicted in Figure 3. This most probably resulted from the large variations in the real milk yields of cows calving during this season. It is difficult to explicitly determine the source of these fluctuations but, in general, a sudden breakdown caused by e.g. diseases, the decreased nutritive value of feeds, staff changes and other unpredictable factors, may (although not necessarily) result in milk yield changes and, consequently, negatively affect prediction accuracy. It can also be noticed that the lactation peak was slightly earlier in season 1 (especially for younger cows) with a more rapid decrease in milk yield, which may indicate inappropriate feeding. In the remaining seasons, a second (although not so evident) lactation peak could be observed, which in turn may reflect additional feeding (in this case, supplementation of the diet with green forage). Considerable milk yield variation occurred in season 4, which again may indicate inappropriate feeding. The predictions do not reflect perturbations visible in the course of real milk yield (especially for SARIMA), which allows us to suppose that they resulted from some unpredictable factors occurring on the farm.

The predictions produced by the SARIMA models, rather poorly reflect seasonal fluctuations in the analyzed series and the predictions are increasingly closer to the trend line with a lengthening prediction horizon. It is probably associated with the fact that seasonality is only modelled by the stationary component of the model (no seasonal differencing). The SARIMA models poorly detect non-linearity in time series data, whereas artificial neural networks are more accurate in this regard, similarly as Wood’s model. The use of artificial neural networks for time series prediction, including additional factors such as calving season, makes it possible to quickly develop a prediction model. The fitting and selection of the best SARIMA model seems to be more time-consuming and absorbing. It requires the verification of many elements (correlograms, parameter significance, the normal distribution of residuals etc.). There are no such requirements for the neural networks, and the use of the automatic network construction module greatly facilitates the development process. When fitting the SARIMA models, it was necessary (for higher accuracy) to divide the dataset into two age groups and four calving seasons, which gave a total of eight models for selection. If the dataset had been larger, it would have been possible to distinguish e.g. 13 age groups (20 to 32 months) instead of just two, which would have resulted in a more accurate prediction, but also the number of the SARIMA models to select from would have increased to 52. The selection itself caused the necessity of verifying many combinations of the p, q, and P, Q parameters. The advantage of the NARX networks in this regard is evident. Finally, the development of Wood’s models was easier compared with other models (NARX and SARIMA) presented in our study.

In conclusion, the applied SARIMA and NARX models predicted the average milk yields at successive stages of the standard lactation of primiparous cows quite well, similarly as the frequently used Wood’s model. It is difficult to indicate the most accurate model in our study, since the predictions did not differ significantly from the real values and among themselves. It is only possible to indicate a slight superiority of individual models. The SARIMA models slightly better predicted the values in the initial prediction periods, whereas NARX and Wood’s models produced more accurate predictions in the later periods, but the indication of the more accurate model is doubtful since the predictions did not differ significantly from the observed values or among each other. The application of the SARIMA model was more time-consuming than that of NARX and Wood’s model.

The application of the prediction methods presented in our study may, to some extent, allow farmers to estimate the milk yield of cows that begin milk production for the first time, providing the basis for the rational management of farm resources, a better definition of tasks and a clearer presentation of the reasons behind each decision made. However, the implementation of such models in user-friendly farm management software would be required first.