Selection of a Mathematical Model to Generate Lactation
Curves Using Daily Milk Yields of Holstein Cows'
L. SHERCHAND,2 R. W. McNEW,~ W. KELLOGG,
and 2. B. JOHNSON
Department of Animal Science
University of Arkansas
ABSTRACT (Key words: lactation curve, Holstein
cows, mathematical models)
Mathematical descriptions of early
stages of lactation were investigated us- Abbreviation key: EMS = error mean
ing daily milk yields of 117 first, 78 squares.
second, 57 third, and 36 fourth lactations
of 120 Holstein cows fitted by 10
models. The measure of fit was the error The graphical representation of milk yield
mean squares, which were replaced by against time is a lactation curve. The curve, a
ranks to perform an analysis of variance concise summary of the pattern of milk yield
with lactation number, model, and period determined by the biological efficiency of a
as factors and with cows as replicates. cow, is used for selection and feeding manage-
The interaction of model and lactation ment. Models estimated from the data may be
number was significant for the fit of the used to predict future milk yields of an in-
entire lactation. A significant interaction dividual or of a herd for the purpose of culling
of model and period was obtained for the or keeping breeding stock. The shape of the
fit of three 30-d intervals. For the entire curve indicates to dairy producers the needed
lactation, the best fit for all four lacta- changes in feeding management; for example,
tions occurred from the diphasic logistic a rising portion of the curve indicates that
function, y = dl(1 - tanh2(bl(nk - cl))) + cows should be given a higher plane of nutri-
d2(1 - tanh2b(n - ~2))). the first 30
For tion, and a declining portion of the curve indi-
d, a modified gamma function gave the cates a lower plane of nutrition.
Computer capability on dairy farms has
best fit for the first lactation, the inverse been improved in recent years, which has per-
polynomial function for the second lacta- mitted automatic recording of daily milk
tion, and the quadratic log function for yields. Curves of lactation data can be gener-
the third lactation. The diphasic logistic ated, and curve parameters can be compared
function gave the best fit for the remain- for management decisions. The need exists to
ing two periods and w s not significantly find mathematical models that best fit lactation
different from the best fitting models for curves.
the first 30-d period. Hence, this function The problem faced at present is the sys-
may be useful to describe the lactation tematic lack of fit of mathematical models,
curve of Holstein cows for dairy herds in especially the gamma function (24), near peak
which the daily milk yield of individual yield. In this study, 10 mathematical models
cows is constantly monitored with a were selected, and each was fit to daily milk
computer. records for an entire lactation. Our objective
was to compare these 10 estimated models for
their adequacy of fit to the lactation curves of
cows, particularly during early lactation.
Received August 8. 1994.
Accepted June 12. 1995.
lSalaries and research support provided by state and MATERIALS AND METHODS
federal funds appropriated to The Arkansas Agricultural
Experiment Station, University of Arkansas. Journal Arti- Mathematical Models
*PO Box 5353, Kathmandu, Nepal. In the models described, y denotes daily
3Reprint requests. milk yield; n denotes time from parturition;
1995 J Diy Sei 78:2507-2513
2508 SHERCHAND ET AL.
and a, b, c, and d denote model parameters.
The first mathematical model to describe the
lactation curve of a dairy cow was proposed by For Model , the parameter a is approxi-
Brody et al. (3) using the exponential decline mately the initial milk yield just after calving,
function: b is the inclining slope parameter up to peak
yield, and c is the declining slope parameter.
y = aeXn. Dl The model estimates a peak milk yield of a@/
which occurs (Wc) weeks after calving
The drawback of this model was that it did not when n is weeks. Although this model can
generate a rise in the curve because of exclu- give an acceptable fit to milk yield data from a
sion of an inclining function. One year later, given lactation (10). it overpredicted the actual
Brody et al. (4) modified [l] to data during early and late lactation and under-
predicted data during midlactation (5, 9).
= =-bn - =-cn, Dave (6) fitted the quadratic model
y = a + bn - cn2 [61
which increased to a peak yield at (c - b)-'ln(c/
b), where In is the natural logarithm. This to first lactation data of Indian water buffalo.
model underestimated milk yield during mid- Cobby and Le Du (5) proposed a model in
lactation and overestimated yield near the peak which the exponential decline in Model  is
and during late lactation (5). replaced by a linear decline:
Sikka (20) proposed a parabolic exponential
model that produced a truncated bell curve for y = a - bn - aeXn. 171
This model predicts peak yield at c-l ln(ac/b)
y = a exp(bn - cn2). 131 with an approximately linear decline thereafter.
Model  fitted to their weekly mean data
where exp is the exponential function. This gave a smaller residual mean square than
model gave good fit for yield during first Wood's model for over half of the 36 lacta-
lactation (21) but did not fit at all before the tions analyzed. Rowlands et al. (17) compared
peak was attained, because the function is Models , , and  and concluded that
symmetric around peak yield. Model  described the initial rise in milk
The inverse polynomial model developed yield up to wk 5 better than Model  but
by Nelder (15) has been used to model lacta- peaked slightly earlier. They also observed that
tion data: Models  and  slightly underestimated,
and Model  slightly overestimated, maxi-
y = n/(a + bn + cn2). mum milk yield. Model  provided the best
141 position of maximum yield.
Madalena et al. (12) used simple linear
This model gave good fit for lactations that regression to describe daily lactation data. This
started at low yield and peaked earlier than model exhibited only a declining slope:
usual (11). Batra (1) observed that the inverse
polynomial function provided a better fit than y = a - bn. @I
the gamma function based on comparison of
R2 when weekly milk data were used. This Molina and Boschini (13) fitted lactation
function was superior to gamma, parabolic, data of Holstein cows to a linear modal model
and exponential functions for abstracting mean that is made up of two straight lines of equal
lactation curves using weekly milk data of but opposite slopes meeting at a peak:
Hariana cows (2, 21), but a gamma function
gave a better fit for cows having a lactation y = a - bln - cl. [91
length of 44 wk (21).
The gamma function proposed by Wood Singh and Gopal (22) referred to a linear
(24) has been used widely for the lactation cum log model for the lactation curve with the
curve of dairy cows. Wood's equation is equation
Journal of Dahy Science Vol. 78. No. 11, 1995
OUR INDUSTRY TODAY 2509
y = a - bn +c ln(n). [lo1 When Models PI, PI, VI, [lo], [ W , U41,
and [U] were fit to daily milk yield for the
Model [lo] is undefined at n = 0 and predicts entire first lactation of Holstein cows, Model
the peak yield to occur at time c/b after calving  had the best fit (smallest EMS) but was
followed by a slow, continuous decline. not different (P> .05) from that of Models ,
Singh and Gopal (22) proposed the quad- [71. and [I41 (19).
ratic cum log model Grossman and Koops (9) first proposed a
multiphasic logistic function, which gives total
y = a + bn + cn2 + d ln(n), [ll] milk by adding each phase of the lactation:
which they found to be superior to Models ,
, and [lo] when fitted to biweekly lactation
data of Indian dairy buffalo. The peak yield is
-b f .\lb2 - 8cd
4c Yn = milk yield on day n,
p = number of lactation phases,
Dhanoa (8) reparameterized Wood's model: tanh = hyperbolic tangent function,
aibi = peak yield for phase i,
y = anmce-Cn [I21 ci = peak day for phase i, and
25-1 = duration in days to attain about
where the parameter m represents time re- 75% asymptotic yield during
quired to reach peak yield. The correlation phase i.
between estimates of m and c in Model  is
smaller than the correlation between b and c in Advantages of the diphasic logistic function
Wood's model. @ = 2) over the gamma function include
Papajcsik and Bodero (16) constructed the smaller and more random residuals, prediction
following six models by pairing suitable in- of t t l 305-d yield, and meaningful functions
creasing functions, such as xb, 1 - exp(-x), of easily interpretable parameters that have
ln(x), and arctan(x), and decreasing functions, biological importance. The triphasic logistic
such as exp(-x) and l/cosh(x): function @ = 3) was preferred to the diphasic
function on the basis of smaller and less au-
y = a nb/cosh(cn), ~ 3 1 tocorrelated residuals (7). The superiority of
the triphasic model was due to lack of fit of
y = a(1 - ean)/cosh(cn), 1141 the diphasic model during the first few weeks
of lactation. The reason for poor fit of the
y = a arctan(bn)/cosh(cn), 1151 diphasic model is that, because of the nature of
the hyperbolic tangent, the diphasic model re-
quires both phases to be symmetric. A sym-
metric phase may not rise quickly enough to
y = a lnfin)/mh(cn), and ~ 7 1 account for a steep rise in early lactation.
Instead of including a third phase, the fit of
y = a arctanfin) e*" [181 the diphasic model might be improved by in-
troducing the power transformation of time,
where arctan and cosh are the arctangent and which Weigel et al. (23) applied to daily milk
hyperbolic cosine functions, respectively. records to account for multiple peaks generated
Papajcsik and Bodero (16) fitted monthly lacta- by the administration of recombinant bST be-
tion data to Models 111 to [lo] and [131 to  ginning at 100 d postpartum. Those researchers
and compared them on the basis of error mean also replaced the product %bi by a new
squares (EMS). Models  and  gave the parameter di. The model including only a sin-
best representation of the lactation curve. gle phase is the monophasic function
Journal of Diy Science Vol. 78, No. 11, 1995
2510 SHERCHAND ET AL.
sumption of the analysis. Computation was
done with PROC GLWSAS, and Type III
Their diphasic model is mean squares were used for tests. Least
squares rank means were used to compare the
When the interaction of model and another
factor was significant (P < .OS), we tried to
isolate the nature of the interaction to facilitate
The exponent k in Models  and  is to model comparison. We tested for interaction of
be estimated from the data. model and subsets of levels of the factor. From
During the process of screening models by this result, we identified those combinations of
fitting some lactation data, Models [ll, , , levels of the factor for which the interaction
161, 181, [91, , , and  were elimi- with model was not significant (P > .OS). We
nated from the study because of poor fit com- then averaged means across these levels for
pared with fit of other models. Model , each model to obtain the means used for com-
which is a reparameterized Wood’s model, was parison of models.
dropped from the investigation to avoid dupli-
cation of the model. RESULTS AND DISCUSSION
Lactation Records Models Fit to the Entlre Lactation
The analyzed data consisted of daily milk When analysis of variance for the effects of
records of 117 first, 78 second, 57 third, and lactation number and models was performed,
36 fourth lactations of 120 Holstein cows. interaction between lactation number and
These cows had been bred and reared in the model was significant (P < .OS). Plots of lacta-
University of Arkansas herd and had their first tion data showed that lactation 1 was different
lactation between 1982 and 1992. Three cows in shape than the other three lactations, which
completed their first lactation far earlier than appeared to be similar. Therefore, we tested
225 d and, hence, were not included in the first the interaction of model and lactations 2, 3,
lactation data. Lactation length varied from and 4 and obtained nonsignificance (P > .OS).
225 to 575 d, but longer lactations were trun- Because of this result, we averaged the least
cated to 305 d. Total lactation yield ranged squares means of these three lactations for
from 5402 to 11,015 kg. Models [41, [51, [71, presenting model comparisons (Table 1).
P O I , U11, 1131, [141, [151, [191, and 1201 were Models were ranked by their least squares rank
fit to the data for each lactation of each cow by means, which were separated using multiple t
nonlinear regression [PROC NLIN using DUD tests.
iteration; SAS (18)]. Each estimated model was The diphasic logistic function with power
used to obtain the EMS for the entire lactation transformation of time (Model ) gave the
and for the first three 30-d periods. The EMS best fit to daily milk records of all lactations,
was chosen as the criterion for assessing model followed by the quadratic cum log function
adequacy. The sources in an analysis of vari- (Table 1). The inverse polynomial function
ance were cow, as a random factor, and the (Model ) and the modified gamma function
fixed factors lactation number, model, period (Model ) had the poorest fits for all lacta-
(when appropriate), and their interactions; the tions, which implied that the inverse hyper-
residual was used as the error term for all tests bolic cosine of Model [I31 did not fit the
involving model. From the initial analysis of declining phase of the lactation curve better
the three early periods of lactation, the three- than the exponential term of the gamma func-
factor interaction was significant (P < .Ol); tion.
therefore, we analyzed separately each lacta- The ranking of the remaining six models
tion that included model, period, and their differed between the first and the other three
interaction. For each data file analyzed, the lactations. The monophasic logistic function
EMS were replaced by their ranks in an at- (Model ), which ranked eighth for fust
tempt to approximate better the normality as- lactation, ranked third for the other lactations.
Journal of Dairy Science Vol. 78, No. 11. 1995
OUR INDUSTRY TODAY 2511
TABLE 1. Rank orders of models based on least squares monly used monophasic functions, such as
rank means for lactation 1 and for the combination of those models described by [41, [51, 171, [lo],
lactations 2. 3, and 4.
[ll], , , , and . This result is
Lactation similar to the findings of Grossman and Koops
Mean of (9) and de Boer et al. (7) in which diphasic or
Rank 1 2, 3. and 4 triphasic functions fit better than the mono-
Fit of Models During the Early
Part of Lactation
To compare models for lack of fit during
the early part of lactation, the first 90 d of
lactation were arbitrarily partitioned into three
periods: first 30 d, d 31 to 60, and d 61 to 90.
An analysis of variance of ranks of EMS for
period and model indicated strong evidence (P
< .01) of an interaction of period and model for
the first three lactations; however, for the first
two lactations, the interaction of models and
periods 2 and 3 was not significant (P > .05).
In contrast, Model  had the third position in For the fourth lactation, the interaction of
fitting the curve for first lactation but placed period and model was not significant (P> .05).
sixth for the other lactations. These six models For a nonsignificant interaction, the least
were not different (P > .05) for lactations 2, 3, squares means were averaged across the ap-
and 4. propriate periods to obtain the means for
Superiority of the diphasic logistic model model comparison (Table 2).
with transformation of time (Model ) over For the first period of the first lactation, a
the nine other models for all four lactations modification of the gamma function Nodel
indicated that the lactation curve may be better [U]j i n which the increasing function of the
fitted by a diphasic function than the com- gamma type, nb, was replaced by arctanfin)
TABLE 2. Rank orders' of models based on least squares rank means for each period of lactation
Lactation 1 Lactation 2 Lactation 3 Lactation 4
Rank PI* P233 P1 P23 P1 P2 P3 P1234
1 [is13 [201. [41. ' [ll]. [20l1 [20P [2018
2 [41ab [Illab [15Iab [l I p ab [lip [lip [I 1 p
3 [lll'b " [ll]ab [19Ik [I5lak [41ak ab ab
4 [20p [S]k ab k akd ak PIab [1 9 p
5 [10P " [14Iab [sik [41t=k [SIkd [15Iab [I4]"b
6 [7Ik [7Ikd ab cd [I4]Ck PIe [4Iab [41ab
7 [191c kd k [I4Icd [SI&' &f ab [slab
8 [141c [15lCd [SIcd [l5lcd [7Icf [ 1SI&' [13Ib [71b
9 [SIC Cd [lO]d [IOId (101' [lolef [101b [IOIk
10 [Wd [13Id [13Ie d [I319 [W' VIb [13IC
~b.c-d*e.r~KM~del~ common superscripts within a column differ (P e .OS).
'Error mean squares were ranked from smallest to largest.
*Period (P): P1 = First 30 d of lactation, P2 = d 31 to 60 of lactation, and P3 = d 61 to 90 of lactation.
3P23 = (€2 + m)n.
4P123 = 8 1 + P2 + P3)/3.
Journal of Diy Science Vol. 78, No. 11, 1995
2512 SHERCHAND ET AL.
and the decreasing function, e*n, was replaced than fourth, which means that the increasing
by l/cosh(cn)-exhibited the best fit but was function, nb, and decreasing function, e-, of
not significantly different (P > .05) from that the gamma type did not fit well for the first 90
for Models , [ll]. and . The slow in- d of lactation. The lack of fit by this function
crease of milk yield after calving, the later during the early part of lactation in our study is
occumng peak than for later lactations, and the in close agreement with the results of others (5,
slow decline after the peak are the characteris- 14, 17).
tics of the first lactation curve. These proper-
ties are well defined by the arctangent and the Relatlonship Between Ranks of Models
inverse hyperbolic cosine of Model , espe- for Entire and Early Lactation
cially for the first 30 d of lactation. The modi-
fied gamma function (Model ) ranked last Model , which gave the best fit to the
among 10 models for each period of four daily milk yields of all lactations, appeared to
lactations, indicating that the shape of the provide the best fit to the daily milk output for
curve was not adequately represented by the the 31- to 90-d period of all four lactations.
combination of nb and inverse hyperbolic co- This diphasic model, however, did not fit the
sine functions. daily milk records of the first 30 d of the first,
The inverse polynomial function (Model second, and third lactations better than the nine
) ranked first for fitting period 1 of the monophasic functions (Table 2).
second lactation but was not different statisti- Model , which ranked fifth for fit of the
cally (P > .05) from Models [ll], , , daily milk information of the first lactation,
, and . This model ranked fifth for the had the best fit for period 1 of all four lacta-
same period of the third lactation. The quad- tions. The eighth ranking of this model for the
ratic cum log function (Model [l l]), which had mean of periods 2 and 3 of first and second
the third rank for the first and second lacta- lactations suggests that the increasing shape
tions during period 1. gave the best fit for the near the peak and the declining stage of lacta-
third lactation during the first 30-d period. tion curve were not better represented by the
This function was not significantly different (P inverse hyperbolic cosine than the log linear
> .05) from that of Models [151, [ 191, and . and exponential decline functions.
The diphasic logistic function (Model ) Model  exhibited the poorest fit for the
gave the best fit during periods 2 and 3 for all entire first lactation but had the second best fit
four lactations but was not significantly differ- during period 1. This result shows that the
ent (P > .05) from Model [ l l ] for these two early fit of the lactation curve is not necessar-
periods. This result suggests that the three ily an indicator of good fit for the entire lacta-
parameter models may not fit the early part of tion. Model , which had the poorest fit to
the lactation curve as well. the daily milk yield of the average of first,
An inconsistency was observed in the order- second, and third lactations, remained in the
ing of models, except for Models [lo], [l 11, last rank for all periods of four lactations. This
, and  for periods 2 and 3 for all ranking means that the composition of the
lactations. For example, Model  was increasing function, nb, and decreasing func-
ranked ninth for first lactation, seventh for tion, l/cosh(cn), of the model did not ade-
second and third, and fifth for fourth lactation. quately fit the shape of the lactation curve.
The ordering of 10 models was also inconsis-
tent during the first period of lactations 1, 2, CONCLUSIONS
The gamma function (Model ) showed The diphasic logistic function with power
lack of fit during the early part of lactation, transformation of time exhibited the best fit to
especially during the first 30 d after calving for the daily milk records of the entire lactation
first and second lactation. For the first period, for all four lactations in which milk yields
the gamma function ranked ninth for first lac- ranged from 5402 to 11,015 kg. The same
tation and eighth for second lactation. During goodness of fit holds for the early part of the
the second and third periods of all four lacta- lactation except for the first 30 d after calving,
tions, the gamma function did not rank better which was fit better, although not significantly
Journal of Diy Science Vol. 78, No. 11, 1995
OUR INDUSTRY TODAY 2513
so, by other models. The logistic function may 10 Kellogg. D. W., N. S. Urquhart, and A. J. Ortega.
be useful to describe the lactation curve of 1977. Estimating Holstein lactation curves with a
gamma curve. J. Dairy Sci. 60:1308.
Holstein cows for dairy herds in which the 11 Kumar, R., and P. N. Bhat. 1979. Lactation curves in
daily milk yield of individual cows is cons- Indian buffaloes. Indian J. Dairy Sci. 32:156.
tantly monitored with a computer. The precise .
12MadaIenq F. E ,M. L. Martinez, and A. F. Freitas.
fitting of the lactation curve will help dairy 1979. Lactation curves of Holstein-Friesian and
Holstein-Friesian x Gir cows. Anim. Prod. 29:lOl.
producers use records to achieve efficient 13 Molina, J. R., and C. Boschini. 1979. Adjustment of
breeding and feeding management of dairy the dairy curve of a Holstein herd with a linear modal
herds. model. Agron. Costarric. 3:167.
14Morant. S. V., and A. Gnanasakthy. 1989. A new
approach to the mathematical formulation of lactation
REFERENCES curves. Anim. Prod. 49:151.
15Nelder, J. A. 1966. Inverse polynomials, a useful
1 Batra, T. R. 1986. Comparison of two mathematical group of multi-factor response functions. Biometrics
models in fitting lactation curves for pureline and 22128.
crossline dairy cows. J. Anim. Sci. 66:405. 16 Papajcsik, I. A., and J. Bodero. 1988. Modelling lacta-
2 Bhat, P. N., R. Kumar, and R.C. Garg. 1981. Note on tion curves of Friesian cows in a subtropical climate.
comparative efficiency of various lactation curve Anim. Prod. 47:201.
functions in Hariana cattle. Indian J. Anim. Sci. 51: 17 Rowlands, G. I., S. Lucey, and A. M. Russell. 1982. A
102. comparison of different models of the lactation curve
3Brody. S., A. C. Ragsdale, and C. W. Turner. 1923. in dauy cattle. Anim. Rod. 35135.
The rate of decline of milk secretion with the advance 18 SAS@ User’s Guide: Statistics, Version 5 Edition.
of the period of lactation. 1. Gen. Physiol. 5:441. 1985. SAS Inst., Inc., Cary, NC.
4Brody. S., C. W. Turner, and A. C. Ragsdale. 1924. 19Sherchand. L.,R. W. McNew, J. M. Rakes, D. W.
The relation between the initial rise and the subse- Kellogg, and 2. B. Johnson. 1992. Comparison of
quent decline of milk secretion following parturition. lactation curves fitted by seven mathematical models.
J. Gen. Physiol. 6541. 1. Dairy Sci. 75(Suppl. 1):303.(Abstr.)
5Cobby, J. M., and Y.L.P.Le Du. 1978. On fitting 20 S i L. C. 1950. A study of lactation as affected by
curves to lactation data. Anim. Prod. 26:127. heredity and environment. J. Dairy Res. 17:231.
6 Dave, B. K. 1971. First lactation curve of Indian water 21 Singh. B., and P. N. Bhat. 1978. Models of lactation
buffalo. J K Res. J. 593.
N W curves for Hariana cattle. Indian J. Anim. Sci. 48:791.
7 D Boer, J. A., J. I. Weller, T. A. Gipson, and M.
e 22Sigh. R. P.. and R. Gopal. 1982. Lactation curve
Grossman. 1989. Multiphasic analysis of milk and analysis of buffaloes maintained under village condi-
yield curves of Israeli Holsteins. J. Dairy Sci. 72: tions. Indian J. Anim. Sci. 52:1157.
2143. 23 Weigel, K. A., B. A. Craig, T. R. Bidwell. and D. M.
8 Dhanoa, M. S. 1981. A note on an alternative form of Bates. 1992. Comparison of alternative diphasic lacta-
the lactation model of Wood. Anim. Prod. 32:349. tion curve models under bovine somatotropin adminis-
9Gmsman. M., W. J. Koops. 1988. Multiphasic
and tration. J. Dairy Sci. 75580.
analysis of lactation curves in dairy cattle. J. Dairy 24 Wood, P.D.P. 1967. Algebraic model of the lactation
Sci. 71:1598. w e in cattle. Nature (Laxi.) 216:164.
Journal of Dairy Science Vol. 78, No. 11. 1995