A Robust Statistical Model to Predict the Future Value of the Milk Production of Dairy Cows Using Herd Recording Data

The future value of an individual dairy cow depends greatly on its projected milk yield. In developed countries with developed dairy industry infrastructures, facilities exist to record individual cow production and reproduction outcomes consistently and accurately. Accurate prediction of the future value of a dairy cow requires further detailed knowledge of the costs associated with feed, management practices, production systems, and disease. Here, we present a method to predict the future value of the milk production of a dairy cow based on herd recording data only. The method consists of several steps to evaluate lifetime milk production and individual cow somatic cell counts and to finally predict the average production for each day that the cow is alive. Herd recording data from 610 Danish Holstein herds were used to train and test a model predicting milk production (including factors associated with milk yield, somatic cell count, and the survival of individual cows). All estimated parameters were either herd- or cow-specific. The model prediction deviated, on average, less than 0.5 kg from the future average milk production of dairy cows in multiple herds after adjusting for the effect of somatic cell count. We conclude that estimates of future average production can be used on a day-to-day basis to rank cows for culling, or can be implemented in simulation models of within-herd disease spread to make operational decisions, such as culling versus treatment. An advantage of the approach presented in this paper is that it requires no specific knowledge of disease status or any other information beyond herd recorded milk yields, somatic cell counts, and reproductive status.


INTRODUCTION
Herd managers cull cows to maintain optimal milk production and maintain herd profitability. For this reason, it is essential for herd managers to continually assess their cows to make timely and appropriate culling decisions.
Herd simulation models used to determine optimal disease-control strategies may benefit by including cow-and herd-specific information, as individual cows will differ in production performance and susceptibility to diseases, and hence their future milk production potential. To make appropriate disease control decisions, it is important to be able to accurately predict a cow's future production. This can be achieved by including prediction algorithms within simulation models.
There exist several methods for estimating the future value and hence the optimal time to cull a dairy cow. Many are based on retention pay-off (RPO) values estimated using marginal net returns (MNR) or expanded to dynamic programming (1)(2)(3)(4). RPO methods are appropriate for determining optimal culling times when complete information on incomes and expenses are included and the herd manager follows the recommendations of the model. This includes disease information and sundry costs. However, information about the health status of individual cows is often scarce, even in countries where individual cow-level health data are collected routinely, as in Denmark. Furthermore, herd managers are not solely motivated by economic incentives but also by habitual and social incentives (5), which may limit the enactment of optimal financial recommendations.
A disadvantage of dynamic programming is the dimensionality problem, in which the number of possible states of a cow grows exponentially with the number of state variables (i.e., production potential, parity, dry off week, and pregnancy) (4). Therefore, a concrete prediction of the future potential of a cow based on its characteristics (such as production and survivability) may provide a useful alternative for cow-and herd-specific decisions, and hence better meet the requirements and expectations of herd managers.
Here, we present a robust method for estimating the future average production (FAP) from any given point in the lifetime of an individual dairy cow. Using this approach, FAP is estimated every time new herd recording data becomes available allowing cows to be ranked according to their expected future profitability. In this way, individual cow culling decisions can be both timely and evidence based.

MATERIALS AND METHODS
Data on milk production, somatic cell counts, and demographic data were obtained from the Danish Cattle database (hosted by SEGES P/S, Aarhus N, Denmark) from which 610 herds with Holstein dairy cows were randomly selected among approximately 3,000 Danish herds participating in milk recording. Cows in each of the selected herds were herd tested either 6 or 11 times per year. Herd testing involved the measurement of milk yield (in kilograms), milk fat (percent) and milk protein (percent), and somatic cell count (expressed as the number of somatic cells per milliliter of milk). Most (91%) of the selected herds had herd testing carried out on 11 occasions throughout the year.
Determining FAP consists of several steps to determine: the average milk production, somatic cell count, survival, and correlation structures between these measures at the herd level. Here, we present how the future values of milk yield and somatic cell count were predicted, and how these values were combined with the predicted lifetime, in order to give the expected future average production (FAP).

Milk Yield
The value of milk depends on the amount of milk produced per day as well as protein and fat content. Daily milk yield was quantified as the number of kilograms of energy-corrected milk [ECM (6)], defined as: ECM = milk (0.122 fat + 0.077 protein + 0.249) (1) where milk is milk in kilograms, protein is protein in percent, and fat is fat in percent, the unit of ECM is kilograms. We required one parameter to describe the milk yield of every cow and preferably a parameter that was stable within and between lactations. To achieve this, individual cow lactation yields were fitted as value relative to the herd average 305-day milk yield lactation curve based on Ref. (7), but extended so that it accounted for day of conception. Curves were fitted for the first, second, and subsequent lactation periods per herd, as a function of days in milk (DIM): where each cow, i, is defined by a milk yield level α M ij that is proportional to the average cow in lactation j in herd k. The parameters a M jk , b M jk , c M jk , and d M jk describe the average cow lactation curve per lactation per herd. The superscript M denotes parameters describing milk production, and is included to differentiate with parameters later introduced with the superscript C which relate to individual cow somatic cell count. I P is the identity operator which is 1 if the cow is pregnant and 0 otherwise. DIM P is the day of pregnancy in units of DIM. This extension that includes the pregnancy date is important when evaluating lactations that are longer than 305 days. This was implemented by first fitting the herd parameters, that returns the average milk yield of a cow in the given herd, then fitting α M parameters for the individual cows.
Given the large variability in milk recordings for the individual cow, we decided that at least six data points per lactation per cow were required to fit individual levels of milk yield (αs). Of the first lactation cows, 90% had at least six samples taken per lactation. This parametrization reduced the description of the individual cow to a single intuitively understandable parameter (α M ), which is the relative production level compared to the average cow in the same herd.

Somatic Cell Count
The somatic cell count (SCC) is an indicator of infection status in the dairy cattle, and Danish milk companies penalize the price paid for milk according to the level of SCC. We shall later include the level of SCC of the individual cow in the FAP. We use the total test day SCC (tSCC) given by tSCC = SCC·milk to model SCC. This method of modeling the somatic cell count as total somatic cells significantly reduced the variance of the residuals compared to modeling as somatic cells per milliliter (8). A one-parameter per cow model of the somatic cell count inspired by the Wilmink curve (9) was fitted as: where α C ij is the tSCC level of each cow i, relative to the average cow in lactation j in herd k, and similar to the milk yield the average cow is described by the herd-level parameters a C jk , b C jk , c C jk , and d C jk . The price of milk in Denmark is adjusted using the number of somatic cells per milliliter, therefore, we reverted back to SCC per milliliter whenever calculating the value of a cow.

Correlation Structures between Lactations
To be able to predict future production levels of individual cows, correlation structures of levels of milk yield and tSCC, α M , and α C , between lactations were determined by fitting a linear regression model. That is, the correlation, ξ M jlk , for α M between lactations j and l in herd k, was determined as the regression coefficients as fitted by lm() using R (10). This means that if the correlation is 0.5, a cow, i, observed to produce 10% more than the average cow in lactation 1 will be predicted to produce 5% more than the average cow in lactation 2,

Survival
To determine the likely future yield of a cow, we need to determine the probability of further lactation periods from any given lactation period in any given herd. The probability of further lactations is framed as survival to subsequent lactations which was determined by observing the fraction of cows producing milk in a subsequent lactation as a function of their milk level, α M , and pregnancy status. An individual cow's survival probability was parameterized as a logistic function to make it more smooth and robust especially in the case of few data points where, e.g., spline fitting might give non-monotone results particularly in the situation where a farm have few cows with large DIM: where S jky takes a sigmoid form and the terms A jky , B jky , ν jky , and M jky are shape parameters (the subscripts jky are not included in equation (4) for ease of reading). The parameter DIM !P is the number of days in milk for non-pregnant cows. S represents the probability of a cow surviving a minimum of one further lactation, given that she is not pregnant at DIM !P . The subscripts indicate that a parameter is determined for all herds k, lactation periods, j, and is categorized by milk yield, y. The number of lactation periods fitted depends on the data available for each herd. The milk yield is divided into groups based on the relative yield α M . The final estimate of survival S jk used in the model was based only on cows with a milk level, α M , above 1, meaning cows producing average or above average compared to the herd level. Cows with a milk level below 1 were observed to have a reduced probability of surviving to further lactation periods, but this may be due to selection bias by the herd manager or indication of other problems such as disease. By comparison cows with a milking level above 1 were observed to have a stable survival independent of further increase in milk level.

Future Yield
Appropriate predictions of the future milk yield of a dairy cow can be made by considering past and present production performance.
Given that the milk yield and somatic cell count levels of a cow (α M and α C ) change over time, we need to identify the optimal balance between past and present. So each recorded measurement of milk or SCC was evaluated against the average cow in that lactation in the same herd giving an α M or α C value, which is represented in the following by x. We, then, used exponential smoothing to continually update the prediction ofx over time: where λ is the smoothing factor, x t is our observed quantity at time t, andx t is the smoothed average at time t. Thex is both our best estimate of the present milk or SCC level of the cow and, therefore, also our best prediction of the future. The level of smoothing as described by the parameter λ was estimated by taking the average over all future measurementsx t>t n for all timepoints during the training period, t n , and finding the value of λ that minimized the squared difference betweenx t n and x t>t n (see Figure 1).
As an initial value,x 0 , we may use information about the yield of the dam. Where this information is not available, we use the value 1 to reflect that our best estimate is that the animal will be the same as the average animal in the herd.

Production Value of Milk
The previous sections describe the parameters required to predict the value of the future milk production of a dairy cow. Now we begin the process of combining this information into a single estimate of FAP.
First, we need to combine the ECM and the SCC to estimate the influence of the SCC on the profitability of the cow. There are several pricing corrections with regards to the SCC in Denmark. Here, we implement the version from the largest Danish milk processor, Arla Foods. Herd managers receive a corrected milk price based on the bulk milk somatic cell count. Prices for milk vary over time, and therefore, the correction schemes are determined in percentage. For this reason, we can consider a 10% reduction in the milk price as being equivalent to a 10% reduction in the milk yield.
The simplest implementation of price correction would be to take the SCC per measurement and directly apply the correction to the ECM produced by the individual cow. However, this method may overestimate the economic consequences, as an animal with a relatively high SCC would affect the price paid to the herd manager less in a large herd than in a small herd, because the SCC is diluted in the bulk tank.
Our proposed method emulates the removal of a single cow from the herd and determines the effect this will have on the bulk tank milk SCC. If removing a cow that produces more cells than are present in the bulk tank, this will change the bulk tank SCC by an amount proportional to the milk this cow produces. So the difference in SCC in the bulk, ∆ c SCC B when removing cow c is given by: Frontiers in Veterinary Science | www.frontiersin.org February 2017 | Volume 4 | Article 13 where we sum over all i cows milking on the same day in the same herd. This method provides an estimate of the contribution of cow c to the bulk tank somatic cell count. Changes in the bulk tank milk from a specific cow are typically small and will not have economic consequences providing the typical structure, which is stepwise, changes in intervals of SCC (i.e., a bulk milk somatic cell count below 200,000 cells per milliliter gives a 2% increase in milk price, whereas a bulk milk somatic cell count between 400,000 and 500,000 cells per milliliter gives a 4% reduction in the milk price). Therefore, we interpolate and extrapolate the values in the pricing scheme so that all changes have an impact (see Figure 2).
The production value of milk (PVM c ) from an individual cow, c, can now be calculated as: where ECM c is the energy-corrected milk production of the individual cow, f C is the interpolated price correction function in Figure 2, and ∑ i ECM i is the sum of lactating cows on a given day, equivalent to the total milk in the bulk tank. All further calculations will be measured in PVM, which has the unit of kilogram milk per day.

Future Average Production-FAP
The future average production is comprised of two parts: the value of milk from the reminder of the current lactation and the value of milk from future lactations.

Value of Milk from the Remainder of the Current Lactation
The value of milk from the remainder of the current lactation can be predicted for two categories of cows: those that are pregnant and those that are non-pregnant. The methods for calculating their values are similar, yet the end points of the lactation are determined in distinct ways. If a positive pregnancy test has been recorded, then the previous insemination attempt is regarded as the day of pregnancy, DIM P . For these cows, the last day of the lactation, DIM L , is set to: If the cow has not yet tested pregnant, the length of the lactation period will be set to: In equations (8) and (9), 282 days equals the average gestation period for Holstein cows, the dry period is 56 days, and the average time before achieving pregnancy is 61 days.
After calculating DIM L for each cow, the estimate for the milk and SCC levels (theαs) are calculated using equation (5). For a first lactation dairy cow with no prior test records, theαs are given by the correlation with theαs of the dam. Using these values and equations (2) and (3), the milk and SCC production of the remainder of the current lactation is calculated. These values are used to calculate PVM values using an equation similar to equation (7), but using the fitted values from the training data.
We note that while we use the predicted values of the individual cow for milk and SCC, we assume that the total amount of milk and the bulk SCC from the herd remains constant, as it would otherwise be necessary to also predict these future values, which would be complicated. The predicted future production of the current lactation, FP C , of a cow is therefore: For simplicity, we have removed the subscripts ijk, but FP C is dependent on the cow, the lactation, and the herd. The predicted future average production for the remaining current lactation is then: where T RC equals the total time remaining in the current parity. We add the 56 days because the average production over time should also include the dry period.

Value of Milk from Future Lactations
The value of future lactations is calculated in a similar way to that of the current lactation. However, here we use a standardized lactation length and the predicted value of the milk and SCC levels from equation (5) that were pre-calculated for the FAP C . We then use the fitted correlations of the milk and SCC levels between lactations, ξ. We can predict the PVM given the predicted ECM and SCC in a similar way to equation (10).
In equation (13) DIM L = 61 + 282 − 56 = 287, which is a simplification of equation (9), but this average lactating period could also be inferred from data for the individual herd. We keep only the subscripts jl as we must specify the lactation period being considered, from lactation j, to future lactation period, l.
where T S is the standard time in a parity.

FAP
Individual cow FAP estimates used for decision-making are calculated by weighing the future production with the different probabilities of surviving from the present lactation [equation (4)] to possible future lactations: (15) where γ jlk is the fraction of cows surviving from lactation period j to lactation period l in herd k.
Explanations of each of the acronyms used in these calculations and typical values of parameters are provided in Table 1. In Table 2 we present a step-by-step description of the process of calculating FAP.

Validation of the Methodology
A subset of the herd recording data was taken to include only milk records for cows born after January 1, 1990. The model was trained on data before January 1, 2010, and FAP estimates were produced for the test period January 1, 2010 to January 1, 2013. Herds were excluded if they did not have an average of at least 50 lactating cows per milk recording date during the training period. This reduced the number of test herds to 388 herds.
Summary statistics of the included herds are provided in Table 3. We compared the FAP estimates calculated using the method described in this paper with a non-parametric approach which involved taking the mean of all future PVM for each cow. The future PVM was weighted with (281 + 61 − 56)/(281 + 61) to account for cows being dried off. The mean of future PVM measurements was compared to FAP, the last measured PVM, the mean of the last three PVM, and the lifetime mean of PVM for each cow. Statistical comparisons were performed, unless otherwise stated, using the Wilcoxon rank sum test with continuity correction. All statistical analyses were carried out using R version 3.1.1 (10).

RESULTS
Line plots showing our estimates of FAP as a function of days in milk for an average cow of parities 1, 2, and 3 in a randomly selected herd are shown in Figure 3. In Figure 4 (left), FAP as a function of days in milk is shown for a parity 1 cow with herd average yield,α M , herd average yield less 10%, and herd average yield plus 10%. Similarly, in Figure 4 (right) FAP as a function of days in milk is shown for a parity 1 cow with herd average total somatic cell count level,α C , herd average total somatic cell count level less 10%, and herd total somatic cell count level plus 10%. First parity cows typically have a lower yield than higher parities; therefore, the FAP is lower in the beginning of the lactation. Later in each lactation, non-pregnant first parity cows have a higher survival probability, and therefore a higher FAP. The overall predictive value of FAP was estimated by comparing the predicted FAP of each dairy cow for each test day in the test dataset, compared to the known future value as described in the section 2.8. The results were averaged within herds, as shown in Figure 5.
In the majority of herds, FAP predicted future PVM within ±2 kg milk per day of the observed values. The SDs of the predictions were also estimated for each herd (Figure 5), which shows that the variance of FAP is small.  (α M andα C ). Solid lines are pregnant cows, dotted lines are non-pregnant. An increased milk production leads to a higher value of FAP (left), while increase in SCC leads to a decrease in FAP value (right). Colors represent change inα M andα C compared to the average cow (black: no change, red: increased by 10%, blue: decreased by 10%). When comparing the prediction of FAP with the latest value of PVM, the FAP predictions were significantly better than future PVMs compared to the latest PVM (P = 0.013). The FAP was on average 50% closer to the future values compared to the latest PVM and more than 70% closer to the future values compared to the average of the latest three and lifetime mean (significantly different with P < 10 −4 ). The SDs of the difference between FAP and mean of future PVM were also significantly smaller than that for the difference between the latest PVM and the mean of latest three PVM (P < 10 −8 ). The difference between FAP and the lifetime mean showed similar SD.

DISCUSSION
Even though replacement models have been available since the 1960s (11) and the theory behind optimal replacement using RPOs was fully formed in the 1980s (1, 2), few herd managers use these models. One reason for this lack of uptake may be that optimized culling using RPO models require herd managers to exercise rational behavior. However, the behavior of individual herd managers is known not to be always rational, but is instead frequently governed by a range of other factors including prior beliefs, values, attitudes, knowledge, and social norms (12). The FAP model is data driven and adaptive so that a change in herd manager behavior is reflected in the FAP values. Changing culling strategy will change the survival curve, changing feed will likely change lactation curves, both of these changes would affect the FAP value. Furthermore, we propose that the FAP is a more readily understood metric for herd managers (since it is expressed in monetary units) compared with RPOs, which are presented to herd managers using metrics that are conceptually more difficult to understand.
The FAP is designed as a tool to support herd managers in making culling decisions on the basis of production data only. We have shown that FAP can accurately predict the future production of milk and SCC combined into a single value, based on herd test information only.
Most herd managers will attempt to inseminate first parity cows for a longer time than cows in later lactations. This behavior leads to first parity cows having a higher survival in the high DIM when open compared with later parity cows. This results in a higher FAP value for first parity cows, as shown in Figure 3. The number of additional times a herd manager will inseminate first parity cows compared with cows in later parities can be herd-specific, and therefore, the survival curves vary among herds.
The presented FAP calculation methodology can be improved and extended. The bias shown in Figure 5 is primarily driven by not knowing the exact values of factors that may vary between herds, such as the length of the dry period, as well as the average number of days from calving to conception. We expect that if such information were available or inferred from data, our FAP estimates may be improved. Thus, further information can be added where it is specifically known. Other factors that may be included in order to improve the model could be seasonal variation both on the day of milking and the day of calving. Factors related to diseases that influence the milk yield may also be included. However, it would then be necessary for herd managers to consistently check their stock for disease and to consistently record the presence of disease.
The genetic value of a cow and its potential future calves could also possibly be included in the FAP. There are several genetic indexes in Scandinavia, an example of which is the Nordic Total Merit index (13). Furthermore, there are many different genetic indexes worldwide [e.g., Ref. (14)], and such information will often enter into the decision process to cull an animal by the herd manager. It is possible that genetic indexes with predictions based on observed values (such as the FAP) are most valuable to a herd manager, but this is dependent on individual management practices.
If large amounts of data are available (e.g., daily measurements of milk production), it may be desirable to track the production using more sophisticated methods such as a Kalman filter, that predicts optimally based on the resent measurement and the observed noise in previous measurements (15). However, the herd recording data used in this study (including milk production and SCC) is collected for individual cows on an almost monthly basis and is, therefore, relatively sparse. For this reason, we opted for the method of exponential smoothing. This method is more robust, less computationally demanding, does not require initialization, and is also less sensitive to non-uniform variation of data over the lactation period.
The FAP is a framework to estimate future production and is as such not dependent on the specific parametrization of the lactation or SCC curves. Therefore, it is possible to use this framework with any lactation or SCC curve desired, though the number of parameters per herd and per cow should be considered when choosing formulations of these curves.
Surprisingly, the latest PVM (which is the most recent observation of the cow) was the second best predictor of the future value after FAP. The latest PVM was better at predicting the future than the average of the last three measurements or the lifetime mean. However, given that milk yield typically decreases toward the end of a lactation, we suspect that the precision of the latest PVM may be most precise in the beginning and more uncertain toward the end of a lactation.
Theoretically, the optimal model for culling is the RPO, which requires full economic parameters (1,4). However, there are large differences between herds for both income given by production values and expenses dependent on many parameters, including herd type, feeding strategies, and milking systems. We also note that many RPO methods include single diseases or production only, most likely due to the dimensionality problem and hence the high computational power demanded by these methods. However, not including all diseases such as mastitis, lameness, and udder malfunction may overestimate the maximum average value of a cow in the herd, leading to non-optimal culling decisions. Diseases are implicitly implemented in FAP, given that the milk production is continually updated via the exponential smoothing algorithm. However, we would expect correlation structures and survival to be different for diseased animals, and this information could be explicitly implemented if such information were available. Therefore, further improvements may be made by investigating how correlation structures for cows with diseases such as mastitis or paratuberculosis differ from healthy cows. Once such information is available, FAP can, in principle, include several diseases, with only little additional computational power required. This would make it possible to implement FAP within-herdspecific models in order to predict the future potential of individual cows for operational decisions, such as treatment and culling, while also taking into account the general strategic management within the herd. This will allow cow-and herd-specific decisions to be made, and is expected to simultaneously improve animal health and herd economy following implementation.
Another way of utilizing FAP predictions would be to compare marginal net production with the average healthy cow, not the optimal cow. This would allow a herd manager to compare the expected lifetime milk yield of given cow with that of a replacement.

CONCLUSION
We have shown how to predict the average future value corrected production FAP, which may be a useful option for making informed culling decisions for dairy cows, given sparse production and disease data and incomplete economic information. The FAP can be used on a day-to-day basis to directly rank cows for culling, or can be implemented in simulation models of disease spread to make inferences on operational decisions, such as culling and treatment, while examining the impact of management strategies of the herd. This is expected to improve animal health and dairy herd profitability.

AUTHOR CONTRIBUTIONS
Developed models, performed data analysis, wrote manuscript: KG. Assisted in developing of models, analysis of data, and drafting of manuscript: LC. Assisted in analysis of data and drafting of manuscript: CK, SN, TH, and NT.