Quantifying inter-group variability in lactation curve shape and magnitude with the MilkBot^® lactation model

Development of the MilkBot^® model

Rook, France & Dhanoa (1993) describe a general lactation model consisting of a growth process multiplied by a death process and a scalar. They stipulate that both growth and death processes should be monotonic, with the growth function rising from zero to approach one, and the death function decreasing from one to approach zero. These constraints have the important effect of making the scalar largely responsible for magnitude of the curve with the growth function dominating shape of the rising portion of a normal lactation curve, and the death function dominating the decline in late lactation, or decay portion of the curve. Note that the growth function returns the cumulative growth as a function of time, not the rate of growth. This may be inferred from the constraints. The multiplicative variant of the Pollot mechanistic model  (Pollott, 2000) follows this general model, but the frequently cited Wood (1967) lactation model, while otherwise similar to the general model, has no upper limit to the growth function, which results in the scaling of the lactation curve being shared between the scalar and the growth process.

In developing the MilkBot^® model, we accept the constraints of the Rook general model and begin by extending it backwards, defining a Growthrate Function, C(t) which models the rate of creation of lactational capacity. This means the Growth function will be the integral of the Growthrate function, and Growthrate will be the derivative of Growth. Equation (1) states the general model as defined, and extended, where Y(t) represents milk production on day t of a lactation with scalar a and G(t), D(t), and C(t) representing growth, death, and Growthrate functions. (1)${\displaystyle Y\left(t\right)=aG\left(t\right)D\left(t\right)=a{\int }_{-\infty }^{t}C\left(t\right)dtD\left(t\right).}$ To maintain conformity with standard practices, the model’s independent variable, t, represents DIM, with t = 1 at parturition, but in a further extension of the general model, and contrary to current standard practice, we note that parturition and the start of lactation are separate events, though related. We define the start of lactation as the midpoint of growth in production capacity, and an offset variable, c, which is the time between parturition and that growth midpoint (the start of lactation). Therefore G(t) = 1/2 when t = c by definition.

Though C(t) is not strictly a probability distribution, it is similarly constrained, with the full-range integral equal to 1. The most prominent probability distribution, the normal distribution, is an obvious candidate for C(t), which would require the rate of creation of udder capacity to rise in a bell curve to a peak near parturition, and then fall again. This has the advantage of being at least moderately credible in relation to observable growth of udder tissue, and easily described by two familiar parameters corresponding to the standard deviation of the Normal curve and the offset location parameter.

The integral of the Normal distribution, the Cumulative Normal function, can only be calculated as an improper integral, which adds considerable mathematical complication. Therefore, simply as a concession to ease of computation, MilkBot^® uses an approximation, the Laplace distribution, for C(t), because it is easily integrated, while yielding a curve of similar shape. The Laplace distribution takes different mathematical forms for left and right sides of the curve, but since the peak is expected to be at or near parturition, and we are unconcerned with the shape of growth before milking commences, the left side of the function can safely be dropped (Eq. (2)). This is equivalent mathematically to postulating that growth rate of udder capacity decreases exponentially from a peak that defines the start of lactation, and with enough high quality data it might be possible to discern whether an exponential or Gaussian growth rate function is closer to reality, but differences between the two are very small after the two weeks of lactation. The resulting growth function (Eq. (3)) ends up mathematically equivalent to the Mitscherlich function described by Rook, France & Dhanoa (1993), but with a different parameterization which emphasizes similarity with the initial Gaussian model. (2)${\displaystyle C\left(t\right)=\frac{1}{2b}exp\left(\frac{c-t}{b}\right)\text{where}t>c\text{(Laplace distribution)}}$ (3)${\displaystyle G\left(t\right)={\int }_{-\infty }^{t}C\left(t\right)dt=1-\frac{exp\left(\frac{c-t}{b}\right)}{2}.}$

For a Death Function, D(t), we choose the same exponential decay function used by Wood, Rook, and others. Dematawewa, Pearson & VanRaden (2007) examined models for extended lactations and found that models featuring an exponential decline of this form outperformed other candidates for lactations exceeding 305 days in length. (4)${\displaystyle D\left(t\right)=exp(-dt).}$ Substitution of Eqs. (3) and (4) in Eq. (1) yields the MilkBot^® Model. (5)${\displaystyle Y\left(t\right)=a\left(1-\frac{exp\left(\frac{c-t}{b}\right)}{2}\right)exp(-dt)\left({\text{MilkBot}}^{\circledR }\text{ model}\right).}$

Experimental data

Groups of 50 lactations each were drawn from 21 randomly selected dairy herds, and then compared to see whether herd effects influence lactation curve shape. Herds were selected from a large DHIA database of herds predominantly in the eastern half of the USA and containing monthly milk weights from more than six million lactations in over 17,000 herds. This is the same database used earlier for developing standard breed-parity aggregate curves (Ehrlich, 2011). Initial selection chose 1,056 herds with at least 1,000 recorded lactations between January 2005 and June 2008. From these, 21 herds were selected randomly. All data points after 305 DIM were dropped, then all lactations with fewer than 7 recorded monthly test days rejected. Lactations were divided in two parity groups, with all lactations after the first in a single group. Each group was then ordered by calving date and one lactation chosen randomly as a starting point. Lactations were then selected forward and backwards in time until 50 lactations had been collected. If a herd could not provide sufficient qualified lactations in either group, that herd was replaced by another randomly selected herd. This resulted in a set of 21 herds, each with 2 parity groups of 50 complete lactations having sequential calving dates. Lactations were all fitted using the DairySight fitting engine (Ehrlich, 2012), and parameter values recorded. Data for each herd-parity group was recorded as comma-delimited data in a single text file, then imported into Mathematica,¹ which was used to generate tables and figures. The number of herds chosen was set to 21 because using the cutoff p < .05 in statistical testing leads to a 1/20 chance of falsely rejecting the null hypothesis. Therefore each herd could be compared to 20 others, with the expectation that on the average there would be one false-positive at p < .05 for each parameter compared. Greater differences suggest greater-than-random variability between herds.

Interpretation of parameters

Parameter “a” is the scale parameter. It can be expressed as kilograms/day, pounds/day, or similarly. Scale can be seen as the theoretical maximum daily yield, which approaches actual peak production as ramp approaches zero. It must be a positive number. Changing the model to a different unit of measure for milk output only requires applying the appropriate conversion to the scale parameter, while all other parameter values remain unchanged. Scale can be seen as describing the magnitude of a lactation, while the other three parameters describe the curve’s shape.

Parameter “b” is the ramp parameter, controlling the width of the Growthrate function, and so the rate of rise in milk production in early lactation. Ramp values are time, normally days, and must be a positive number. A simple thumb rule is that when t = r a m p, growth in production capacity is about 82% complete.

Parameter “c” is the offset parameter describing the offset in time between parturition and the start of lactation. Offset values are time (days), and may be positive, negative, or zero. The effect of offset is slight, except in the first few days of the lactation curve, so it is generally not possible to detect variation in the offset parameter in fitted lactations unless there are daily milk weights covering the first weeks of lactation. With a monthly interval between test points, values for the offset parameter have little value, and the model may be simplified by setting offset equal to zero or a constant near zero. Fixing offset to zero is equivalent to the common assumption that lactation begins at parturition.

Parameter “d” is the decay parameter, controlling the rate of senescence of production capacity. Decay is inverse-time (days⁻¹). It should be constrained to positive values under normal circumstances. An equivalent alternative expression for a first-order decay constant like the decay parameter is as a half-life, which suggests a definition of lactation persistency, which, rather than being tied to an arbitrary stage of lactation, is an attribute of the lactation as a whole. (6)${\displaystyle Persistence=\frac{0.693}{decay}.}$ By this definition decay is the time it would take for production to drop by half, if we were to ignore the growth side of the model. Since the Growth function approaches one in late lactation, decay, by this definition, is close to the actual half-life for milk production in late lactation. Fitted decay values are likely to be more normally distributed than persistence, therefore decay should be preferred for most statistical calculations, but may be converted to persistence afterwards.

Mathematical manipulation of Eq. (1) allows calculation of some useful results. By setting the derivative equal to zero we can calculate peak day, t_peak, and from that peak milk, y_peak. Note that like persistence, t_peak and y_peak are attributes of the lactation as a whole, and less sensitive to any single data point than metrics based on comparing individual test days. (7)${\displaystyle t_{peak}=c-b\hspace{0.17em}log\left[\right(2bd)/(1+bd\left)\right]}$ (8)${\displaystyle \begin{array}{ccc}\hfill y_{peak}& \hfill =\hfill & Y\left(t_{peak}\right)\hfill \\ \hfill & \hfill =\hfill & aexp(-d(c-b\hspace{0.17em}\mathrm{Log}\left[2\right]\left)\right[a,b,c,d\left]\right)\hfill \\ \hfill & \hfill \hfill & (1-1/2exp((c-(c-b\hspace{0.17em}\mathrm{Log}\left[2\right]\left)\right[a,b,c,d\left]\right)/b\left)\right).\hfill \end{array}}$

Cumulative production between two days can be calculated by integration, including M305. (9)${\displaystyle \text{M305}=(a-aexp(-305d\left)\right)/d+(abexp(c/b\left)\right(-1+exp(-305(1/b+d\left)\right)\left)\right)/(2+2bd).}$

Comparison of herd-parity groups

Herd-parity groups were compared to establish whether MilkBot^® parameters or derived metrics (RMSE, M305, peak day, peak milk) could be used to detect statistical differences in lactation curves between groups. It can be inferred that if such differences exist between herds, they probably are caused by differences in management or environment, and this suggests that statistical analysis of fitted parameter values or derived metrics may be a means of monitoring herd performance and detecting abnormalities. Existing literature firmly establishes that decay should be lower in first-parity groups (higher persistence in heifers), and that there are likely to be differences in measures of magnitude such as scale, peak milk, and M305 between herds and parities. It is much less certain whether management and environment influence persistency (decay), ramp, offset, RMSE, and time of peak milk, or whether these might be fixed biologically, and insensitive to management.

Table 1 gives herd mean M305, predominant breed, and links to files containing data on 50 individual lactations for each of the 42 groups. Table 2 summarizes data for all herd-parity groups for each parameter or derived metric. These are overall means of the group mean and group standard deviation, and mean divergence score. A “divergence score” was calculated for each group by tabulating how many of the 20 other groups of matching parity had significantly different mean parameter value by Z-test at P < .05, so that 1.0 would be the expected divergence score if there are not differences among groups, and higher scores suggest divergence is more common than would be expected from random variation.

Figures 1–8 show box-and-whiskers plots for M305, scale, ramp, offset, decay, RMSE, peak day, and peak milk. In each plot, herds are ordered by herd-average M305, with the lowest-production herd (Herd A) at the bottom. Parity groups are differentiated by color. For each group, a gray diamond gives the Z-test 95% confidence interval around the mean. Colored bars show quartiles above and below the mean, and whiskers show the full range for the group. This means that groups with diamonds that overlap vertically are not significantly different by Z-test with 95% confidence. If diamonds do not overlap, the sample means differ at p < .05. Groups as small as 50 lactations with monthly test data are clearly sufficient to identify significant variability among herds, with the exception of offset and possibly peak day. Scale, M305, and peak milk seem to follow similar patterns, suggesting the correlation that is to be expected, as all are primarily measures of lactation magnitude.

Ramp is higher for heifers, indicating a slower rise in production after calving. It can be speculated that inter-herd variability in ramp might be influenced by transition management. Some herds (C, I, L, N, R, T) show little difference in ramp between parity groups, and since transition heifers are often managed differently from cows, it is easy to develop hypotheses on possible causes. For example, the differences could relate to whether fresh heifers are group-housed with older animals since competition with older animals might limit feed consumption in heifers, and slow growth in production capacity. That the same herds tend to show a wide range in ramp values among mature cows adds another angle for speculation and research.

There is little difference in offset within or between herds, which is not surprising since only monthly milk weights were available and biologically significant differences in offset would have little effect past the first few days of lactation. Even so, in mature cows there are statistically significant differences in mean offset of as much as 0.3 days (8 h). It is very unlikely that this reflects biologically or economically significant differences in offset. Rather, with only monthly data points the fitting engine has difficulty differentiating offset from ramp, so makes the appropriate choice of attributing most but not all variability to ramp. This leads to correlation between offset and ramp as an artifact of the fitting process. One simple solution is to ignore offset when monthly test data is used. The untested possibility remains that true differences in offset would be found with daily milk weights.

Heifers have markedly lower decay (greater persistency) than mature cows, but there is also considerable variation between herds. For example, persistency of mature cows in Herd K is slightly better than heifers in Herd C. The difference between cows and heifers also varies widely among herds. Factors that might influence persistency include use of rBST (recombinant bovine somatotropin), nutrition, mastitis control, and many other variables in management, environment, and genetics. These could have considerable economic importance through the effect on overall production.

Distribution of RMSE is interesting because greater consistency either in the testing protocol generating production data or in actual production should lower RMSE. Mastitis, lameness, poor feed management, and many other problems would likely increase RMSE so that it may be useful as a general, if imperfect, measure of management quality and cow welfare. These are extremely difficult to measure by objective criteria, so even imperfect methodology may be useful. For example, Herd T has high production but also relatively high RMSE, which might reflect management problems, or just lower quality data collection. There are a few outliers, such as individual lactations in herds N and T with very high RMSE. For example herd N includes a single lactation with RMSE of 24.8 kg and milk weights recorded as (47.7, 83.1 84.4, 32.2, 8.2, 55.4, 77.2, 64.5, 34, 24.1) kg at approximately monthly intervals. This is obviously a highly abnormal lactation, and probably a cow with health problems. Some portion of RMSE should be attributed to bias inherent in the model design and fitting process, since no model or fitting process is perfect. This systemic error should be randomly distributed among herds, so the inter-herd variability in RMSE suggests at least that systemic error is not overwhelming. The remaining portion is sometimes called random error, reflecting the high variability that is expected even between consecutive milkings. This probably is a misnomer, since most of the variability probably has a cause, whether it is inconsistency in milking times, variability in feeds, meter error, weather, disease, or any of thousands of other possible causes. Clearly some of these are controllable, and some are not, but RMSE allows them to be measured, at least collectively.

Peak milk and peak day will correlate with the parameters from which they are calculated. It is possible that they may turn out to correlate better with certain variables of interest than individual parameters, but that remains to be seen. Peak milk shows a pattern similar to scale, but with more difference between cows and heifers because of the lower decay which is characteristic of heifers. Peak day is not greatly influenced by production level, and shows relatively little variation among herds, suggesting that it is more fixed biologically, and mostly insensitive to management and environment.

Peak milk is sometimes used with the thumb rule that each unit of peak milk translates to 250 units in M305, or some similar formula. This presupposes that lactation curve shape is the same in all herds, and is not influenced by whatever management changes might improve peak milk. We can get a feel for the practical effect of herd variability in curve shape by looking at its effect on this thumb rule. A rather complicated formula² can be derived to calculate how much change there would be in M305 if scale were increased enough to increase peak milk by one unit. This depends on values for the other three parameters, especially decay. For Herd C, for example, a unit of peak translates to 230 units of M305 for heifers, and 188 for cows while in Herd Q these are 267 and 221 units respectively.