The geometry of the lactation curve based on Wood's equation: a two-step prediction

ABSTRACT

Lactation records from cows of the southwestern Paraná state, Brazil, form the dataset of this study. We applied the information-theoretic approach to evaluate the ability of the nonlinear Wood, Brody, Dijkstra, and Gamma functions to fit to these data by employing a two-step technique based on nonlinear mixed-effects models and generalized linear mixed-effects models. Wood's equation was fitted with the combination of a first-order autoregressive correlation structure and a variance function to account for heteroscedasticity. This version was the best choice to mimic lactation records. Some geometric attributes of Wood's model were deduced, mainly the ascending specific rate from parturition to peak milk yield and the descending specific rate as a measure of the lactation persistence of the milk yield at peak production. Breed and parity order of the cows were assumed as fixed effects to obtain a reliable model fitting process. Regardless of breed, first-order parity cows had greater persistency than their older counterparts, and the greater the ascending rate of milk yield from the parturition to the peak, the sharper the decrease in milk yield post-peak; therefore, the rates (absolute values) of ascending and descending phases correlated positively. Nonetheless, the actual estimated values of the descending phase rates are negative. Wood's equation was flexible enough to mimic either concave- and convex-shaped lactation profiles. The correlations between both peak milk yield and random estimates for β with total milk yield per lactation were positive. However, peak milk yield might not be the only variable used for ranking cows; the total milk yield integrates all information of the lactation profile through the estimated parameters of Wood's equation.

Keywords:
dairy cattle; generalized linear mixed-effects model; lactation models; nonlinear mixed-effects model; nonlinear parameter estimation

1. Introduction

Lactation functions provide an elegant example to study the applicability of mathematical models to explain nonlinear phenomena in animal science. Scientists have proposed several mathematical functions to mimic milk production records over time ( Brody et al., 1924 ; Wood, 1967 ; Grossman and Koops, 1988 ; Rook et al., 1993 ; Dijkstra et al., 1997 ; Pollott, 2000 ). Specialized dairy breeds typically exhibit convex-shaped lactation profiles with ascending, peak, and descending phases ( Brody et al., 1924 ). Nevertheless, there are literature reports that present the difficulty of some mathematical models to fit atypical lactation profiles, which include profiles with an insufficient number of records to characterize the entire lactation, concave-shaped profiles, as well as those profiles that present an irregular distribution of milk records over time. Furthermore, some studies have emphasized that genetics, diet, parity order, season, metabolic disorders, diseases, and other environmental aspects are responsible for the occurrence of atypical lactation profiles ( Wood, 1968 , 1970 , 1972 , 1976 , 1980 ; Grossman and Koops, 1988 ; Olori et al., 1999 ; Rekik and Ben Gara, 2004 ; Dijkstra et al., 2010 ; Hossein-Zadeh, 2016 ; Ahmed et al., 2019 ). These factors might explain why some models are unsuccessful in fitting some lactation records. Therefore, there is a need to investigate flexible models that are applicable to a wide range of situations to describe the nonlinear phenomenon of lactation. Among them, we can list some models that fitted lactation profiles adequately ( Ferreira et al., 2015 ), i.e., the models proposed by Brody et al. (1924) , Wood (1967) , and Dijkstra et al. (1997) , and we suggest the Gamma probability density function as a generalized form of Wood's model.

The development of tools from theories of nonlinear mixed-effects models and generalized linear mixed-effects models ( Pinheiro and Bates, 2000 ; Littell et al., 2006 ; Vonesh, 2012 ; Stroup, 2013 ) has provided instruments to describe several nonlinear phenomena in animal science ( Zanton and Heinrichs, 2009 ; Strathe et al., 2010 ; Albertini et al., 2012 ; Araujo et al., 2015 ; Vieira et al., 2018 , 2020 ; Rohem Júnior et al., 2020 ). The information-theoretic approach (I-T) provides the elements for evaluating the quality of fit of models by aiming at their predictive powers ( Buckland et al., 1997 ; Sober, 2002 ; Burnham and Anderson, 2004 ; Burnham et al., 2011a ). In this regard, we can use the geometric properties of mathematical functions to describe lactation records through a two-step prediction based on the statistical methods described. This approach can yield robust predictions about several biological mechanisms that explain the lactation course. Interesting examples rely on the quantitative representation of cell proliferation and death in the mammary gland, whose effects on milk secretion concur for the rise, peak, and fall observed in typical lactation records ( Knight and Peaker, 1984 ; Dijkstra et al., 1997 ; Capuco et al., 2001 ). The advantage of geometric functions and the two-step prediction based on the described statistical methods is the feasible estimation of parametric functions represented by geometric properties without increasing the number of parameters of a given model.

Therefore, our goal was to evaluate the quality of fit of Brody, Wood, Dijkstra, and Gamma models to describe lactation records from different cattle breeds and parity orders of the cows and, by deducing and using geometric functions, to predict essential attributes of the lactation function by employing a generalized two-step prediction technique.

2. Material and Methods

Lactation records from cows of the southwestern region of the Paraná state, Brazil, were used to develop the evaluation dataset, which was partially used in a previous study performed by our research group about curve fitting to lactation records of Holstein cows ( Ferreira et al., 2015 ). Here, we used a larger dataset and presented descriptive statistics by breed ( B ) and parity order ( PO ) ( Table 1 ). We used the entire dataset containing 6754 lactation profiles. We discarded no lactation records to challenge the ability of the equations studied to fit the data at hand.

Table 1
Breeds, parity orders, number of cows ( n _c ), number of time records ( n _t ), mean, standard deviation (SD), and range limits for the dataset

Breed	Parity order	n c	n t	Mean	SD	Minimum	Maximum
Girolando (1)	1	24	83	13.5	6.49	4.0	29.5
	2	4	25	17.4	6.35	6.0	27.6
Holstein (2)	1	3166	16903	22.8	8.37	0.5	70.0
	2	1502	9513	26.8	9.37	2.4	71.0
	3	600	3824	28.7	9.78	2.0	69.3
	4	212	1342	30.2	9.50	2.3	57.2
	5	80	526	31.4	9.90	5.5	61.8
	6	24	207	31.2	9.50	8.2	62.0
Jersey (3)	1	607	2720	15.2	6.21	2.0	37.5
	2	199	1151	17.9	6.79	1.0	40.5
	3	26	104	19.5	9.28	5.0	39.5
	4	3	12	15.1	3.21	12.7	24.5
Crossbreed (4)	1	187	742	14.9	7.44	2.5	53.1
	2	54	268	19.1	7.01	5.3	38.6
	3	10	39	18.8	6.65	8.0	31.9
	4	2	15	22.1	6.64	14.6	37.0
Brown Swiss (5)	1	35	209	19.0	7.68	3.8	38.5
	2	15	95	22.7	8.58	7.5	51.0
	3	4	32	22.9	10.56	7.3	46.0

2.1. Geometric attributes of Wood's equation

Wood (1967) reparameterized the gamma probability density function to yield a nonlinear equation that mimics records over the lactation course of dairy cows with a satisfactory fit ( Ferreira et al., 2015 ). The equation can take the form of the following expression:

Eq. (1) $${\mu }_t=\alpha t^{\beta }\exp \left(-\lambda t\right)$$

The parameter μ _t represents an expected value, that is, a population mean about random lactation records over time ( t , days or weeks) from one cow or, preferably, from a random sample of cows taken from a larger population. Parameters α , β , and λ are constants, and μ _t is a nonlinear function. Isolated, λ is a fractional rate of decline expressed as reciprocal of time units ( t ⁻¹), and a formal differential equation ( Dijkstra et al., 2005 ) reveals that the specific fractional rate is time-dependent, actually:

Eq. (2) $${\kappa }_t={\mu }_t^{-1}\cdot d{\mu }_t/dt=\lambda +\beta t^{-1}$$

Equation 2 represents a time trajectory instead of a single parameter that otherwise would offer a direct biological interpretation for both theoretical and practical purposes. In fact, λ is the asymptotic fractional rate for ${\kappa }_t\left(\underset{t\to \infty }{lim}{\kappa }_t=\lambda \right)$ . Therefore, Wood's model presents some attributes revealed by analytic geometry that may be of use for animal scientists ( Figure 1 ).

Figure 1
Geometric attributes of Wood's model.

We choose equation 1 because of its flexibility and feasibility to mimic lactation records as demonstrated by several authors ( Dijkstra et al., 2010 ; Ferreira et al., 2015 ; López et al., 2015 ; Ahmed et al., 2019 ). The geometric properties of equation 1 are nonlinear functions of its parameters. The time to peak milk production ( t _p ) and peak milk production ( $\left({\mu }_{t_p}\right)$ ) have been already described ( Wood, 1967 ):

Eq. (3) $$t_p=\beta {\lambda }^{-1}\text{and}$$

Eq. (4) $${\mu }_{t_p}=\alpha t_p^{\beta }\exp \left(-\beta \right)$$

These two attributes are time-independent. The inflection point ( t _i ) of equation 1 is one of two possible solutions (i.e., roots) obtained by equating the second derivative of equation 1 to zero. These two solutions are also time-independent:

Eq. (5) $$t_i=\left(\beta \pm {\beta }^{0.5}\right){\lambda }^{-1}$$

Only the root t _i > t _p is valid. All parameters must be subjected to the same constraints described by Wood (1967) . As a corollary of the parametric space for β , that is, β ≤ 1, lactation records with ascending and descending phases about peak milk yield are those with a constrained parametric space, i.e. 0 < β ≤ 1; whereas those records that show only a descending phase starting from the onset of lactation at parturition are those for which the parametric space is constrained to β ≤ 0 (concave profiles). We remind the reader that equation 1 resumes to the exponential model for β = 0. The cumulative or total milk production ( TMμ _t ) is the area under the curve used to describe the lactation records. Here, we provide an analytical solution that yields the definite integral as follows:

Eq. (6) $$TM{\mu }_t={\displaystyle \int {}_l^u\alpha t^{\beta }\exp }\left(-\lambda t\right)dt=-\alpha {\lambda }^{-\left(\beta +1\right)}\left(\Gamma \left[\beta +1,u\lambda \right]-\Gamma \left[\beta +1,l\lambda \right]\right)$$

Numerically, the upper limit ( u ) of the integral was set at 305 d, and the lower limit ( l ) can be a small value greater than zero (10⁻⁶) for numerical integration. The incomplete gamma function is $\Gamma [\beta +1,\lambda t]=\Gamma [\beta +1]–\gamma [\beta +1,\lambda t]$ . The gamma function $\Gamma \left[\beta +1\right]={\displaystyle {\int }_0^{\infty }{t}^{\beta }}$ exp(− t ) dt is a numerical result from a definite integral, and the lower incomplete gamma function is given by $\gamma \left[\beta +1,\lambda t\right]={\displaystyle {\int }_0^{\infty }{v}^{\beta }}$ exp(− λv ) dv . The average milk production was computed as ${\overline{\mu }}_{0,305}=TM{\mu }_t/305$ .

The geometry of the function depicted in Figure 1 can yield a fractional rate constant of decline with reference to peak milk yield. The dashed descending line tangent to the inflection point after the time for peak milk yield is the geometric place defined by equation 7a:

Eq. (7a) $${\mu }_t={\mu }_{t^{\text{'}}}+\left(t-t\text{'}\right)\cdot {\left(d{\mu }_t/dt\right)}_{t_i}$$

in which ( t ', ${\mu }_{t\prime }$ ) is the point where equation 7a crosses the tangent line to the point ( t _p , ${\mu }_{t_p}$ ), and ${\left(d{\mu }_t/dt\right)}_{t_i}$ is the value of the first derivative of equation (1) at t _i ( Figure 1a ). The tangent to ( t _p , ${\mu }_{t_p}$ ) is parallel to the abscissa axis, i.e., ${\mu }_{t\prime }={\mu }_{t_p}$ , and the abscissa coordinate t ' can be isolated in equation 7a and solved for t = t _i . Finally, if we divide both sides of equation 7a by ${\mu }_{t_p}$ and manipulate signs for κ _d because ${\left(d{\mu }_t/dt\right)}_{t_i}<0$ , for not changing the trend of equation 7a , we have the persistency trend over time ( P _t ) as follows:

Eq. (7b) $$P_t={\mu }_t{\mu }_{t_p}{}^{-1}={\mu }_{t\text{'}}{\mu }_{t_p}{}^{-1}+{\left(d{\mu }_t/dt\right)}_{t_i}{\mu }_{t_p}{}^{-1}\left(t-t\text{'}\right)=1-{\kappa }_d\left(t\text{'}-t\right)$$

Therefore, the specific rate constant of decline after peak milk yield arbitrarily becomes the positive value, as follows:

Eq. (8) $${\kappa }_d=-{\mu }_{t_p}{}^{-1}{\left(d{\mu }_t/dt\right)}_{t_i}$$

The new parameter κ _d can be d⁻¹ or w⁻¹ and dependent on t _i ( equation 5 ) as an intrinsic attribute of equation 1 . We can relate the fractional rate described by equation 8 to the phenomenon of stromal apoptosis in the mammary gland ( Stefanon et al., 2002 ; Dijkstra et al., 1997 , 2010 ). The relative rate of decline midway between peak milk yield and the arbitrary end point of lactation ( t _f = 305) is an additional rate computed from the parameters of Wood's equation, namely $r_d=2\beta \lambda {\left(\beta -\lambda t_f\right)}^{-1}-\lambda$ ( Dijkstra et al., 2010 ).

Another critical factor of the lactation might be the rate of DNA proliferation in the mammary gland, which would be responsible for the rise in milk production after parturition ( Dijkstra et al., 2010 ). We can geometrically generate a point estimator for an ascending rate of milk production by calculating the average milk production during the ascending phase until peak milk yield, as follows:

Eq. (9) $${\overline{\mu }}_{{10}^{-6},t_p}=\left(TM{\mu }_t-{\displaystyle \int {}_{t_p}^{305}}\alpha t^{\beta }\exp \left(-\lambda t\right)dt\right)/\left(t_p-{10}^{-6}\right),$$

for ${\displaystyle {\int }_{t_p}^{305}\alpha t^{\beta }}\exp \left(-\lambda t\right)dt=-\alpha {\lambda }^{-\left(\beta +1\right)}\left(\Gamma \left[\beta +1,305\lambda \right]-\Gamma \left[\beta +1,\lambda t_p\right]\right)$ . By computing the tangent of the δ angle ( Figure 1b ), as equation 10 , we have:

Eq. (10) $$\tan \delta =\left({\mu }_{t_p}-{\mu }_{t={10}^{-6}}\right)/\left(t_p-{10}^{-6}\right),$$

and the following ratio yields a constant, fractional ascendant rate of milk production given by equation 11 .

Eq. (11) $${\kappa }_{\alpha }=\tan {\delta /\overline{\mu }}_{{10}^{-6,t_p}}$$

This rate can be d⁻¹ or w⁻¹. The rate κ _α can be associated with the process of cell proliferation until peak milk yield ( Dijkstra et al., 1997 , 2010 ). The constraint 0 < β ≤ 1 must hold for a valid κ _α ; nonetheless, even though estimates for the fixed parameter β are negative, there are possible occurrences of positive values from the random cow effect over parameter β that may result in typical ascending-peak-descending (convex-shaped) lactation profiles.

2.2. Brody's lactation equation

Brody and colleagues proposed a model to mimic the standard course of milk flow of dairy cows, that is, the rise of milk flow after parturition and its steady decline after peak milk yield ( Brody et al., 1924 ). Their model describes a biphasic phenomenon, in which the rise follows a monomolecular change and, simultaneously, an exponential decay follows the peak milk yield. Mathematically, the final equation has the form:

Eq. (12) $${\mu }_t={\theta }_1\exp \left(-{\kappa }_{1}t\right)-{\theta }_2\exp \left(-{\kappa }_{2}t\right)$$

The mean μ _t has the same meaning as previously described, but θ ₁ and θ ₂ are scale parameters, κ ₁ (d⁻¹ or w⁻¹) is the characteristic constant of the decline of milk production after peak milk yield (“after the second month” in those authors' words), and κ ₂ (d⁻¹ or w⁻¹) is the characteristic diminishing constant of the rising course of lactation. It is possible that equation 12 mimics concave-shaped lactation profiles, but this will depend on final parametric estimates: if θ ₁ > 0 and θ ₂ ≤ 0, the equation resumes to an exponential decay. All derivations applicable to Wood's equation also apply to the Brody's equation.

2.3. Dijkstra's lactation equation

There is an elegant mechanistic equation based on the pool size of DNA in the cell population of the mammary gland during lactation. This model roots on cell proliferation and death and mimics the typical convex-shaped lactation profile as follows:

Eq. (13) $${\mu }_t={\mu }_0\exp \left({\theta }_T{\kappa }_3{}^{-1}\left(1-\exp \left(-{\kappa }_{3}t\right)\right)-{\kappa }_{4}t\right)$$

Parameter μ ₀ (kg · d⁻¹ or kg · w⁻¹) is the theoretical initial milk production at parturition, θ _T (d⁻¹ or w⁻¹) represents the specific rate of cell proliferation at parturition, κ ₃ (d⁻¹ or w⁻¹) is a decay parameter, and κ ₄ (d⁻¹ or w⁻¹) is the specific rate of cell death. Dijkstra et al. (1997) showed several attributes of equation 13 as general nonlinear functions of its parameters.

2.4. The Gamma density function used as a lactation equation

The gamma probability density function (PDF) models time-to-event data ( Mood et al., 1974 ; Stroup, 2013 ). Because Wood (1967) reparametrized the gamma density function as previously mentioned, we believe that it is worthy evaluating this model in its original form because the properties of the gamma PDF may be of use to describe lactation records. We can compute the mean and variance based on this lactation function. The model can have the following form:

Eq. (14) $${\mu }_t={\alpha }_0{\left(\lambda t\right)}^{r-1}\lambda \exp \left(-\lambda t\right)/\Gamma \left[r\right]$$

The mean μ _t represents the lactation course, and α and r are scale and shape parameters, respectively. Parameter λ (d⁻¹ or w⁻¹) has the same meaning as previously described. The constraints α ₀ > 0, r > 0, and λ > 0 must hold. The equivalence between equations 1 and 14 is trivial, because of β = r − 1, α = α ₀ λ ^r / Γ [ r ], and λ is the parameter in common. One can use the same rationale for Wood's model to obtain the nonlinear parametric functions of interest.

2.5. Model fitting: first-step prediction

We used the nlme function of R to fit equations 1 , 12 , 13 , and 14 to the lactation records. The gnls and nlme functions belong to the nlme package ( Pinheiro et al., 2017 ). The stochastic version of the models was set as follows:

Eq. (15) $$y_{t_{ijk}}=g\left({\Theta }_m,t_{ijk}\right)+e_{t_{ijk}}$$

The record $y_{t_{ijk}}$ corresponds to the milk production rate (kg · d⁻¹) for the k -th cow of the i -th breed, during the j -th parity order, and recorded at time t . The basic nonlinear function that fits $y_{t_{ijk}}$ is $g\left({\Theta }_m,t_{ijk}\right)={\mu }_{t_{ijk}}$ , ∀ μ _t related to equations 1 , 12 , 13 , and 14 . Parameter Θ _m represents the vector of parameters of the m -th model. However, prior to the fit of the model equations, we had to group lactation records by the random cow and cow × parity order intersecting effects by using the groupedData function of nlme. The resulting formula for the grouped data was y ~ Time|cow/cpo, i.e., y is daily milk yield, time is t , cow and cpo are the random cow and cow × parity order (cpo) intersection. Breed ( B ), parity order ( PO ), and their interaction ( BPO ) were the fixed effects associated with each parameter of μ _t for each model evaluated.

We modeled the variance associated with the error term ( $e_{t_{ijk}}$ ) according to four conditions:

Eq. (16) $${\sigma }_t^2={\sigma }^2,$$

Eq. (17) $${\sigma }_t^2={\sigma }^2{\left|g\left({\Theta }_m,t_{ijk}\right)\right|}^{2\psi },$$

Eq. (18) $${\sigma }_t^2={\sigma }^2{\left(\omega +{\left|g\left({\Theta }_m,t_{ijk}\right)\right|}^{\psi }\right)}^2,\text{and}$$

Eq. (19) $$cov\left(e_{t_{ijk}},e_{t{\text{'}}_{ijk}}|u_m^{\left(p\right)}\right)={\sigma }^2{\varphi }^{\left|t_{ijk}-t{\text{'}}_{ijk}\right|}$$

Parameters represent the residual variance ( σ ²), the power ( ψ ) that scales the absolute mean | g ( Θ _m , t_ijk )|, and the assumed positive correlation parameter (0 < ϕ < 1). The model fitted with the variance function in equation 17 challenges the fitted one based on the traditional homoscedastic assumption described by equation 16 . Parameter ω is the intercept of the VarConstPower function ( equation 18 ), which allowed modeling the variances when | g ( Θ _m , t_ijk )| = 0 or equal to a baseline or background value, and when the variance scales to the mean for | g ( Θ _m , t_ijk )| > 0 or greater than a baseline value ( Pinheiro and Bates, 2000 ; Rohem Júnior et al., 2020 ; Vieira et al., 2020 ). The fit of the model with the first-order autoregressive correlation structure given by equation 19 challenges the traditional assumption of independent time records on the same subject ( Vonesh, 2012 ); this was accomplished by fitting the corCAR1 function of nlme package ( Pinheiro and Bates, 2000 ). We also checked for possible diagonal and symmetric covariance structures for random effects by using the pdDiag (diagonal, D) and pdSymm (unstructured, U) covariance classes available in the R software ( Pinheiro and Bates, 2000 ).

The models were initially fitted with the gnls function with no random effects introduced. Thence, by using the nlme function, we gradually introduced random effects (cow and cow × parity order) in the different fixed parameters contained in g ( Θ _m , t_ijk ) creating several model versions. An example is the introduction of the random effects over parameter λ in Wood's equation, in sequence over parameters λ and β , and so on for all equations studied, which formed the different model versions. In addition, we assessed the quality of fit of the model solutions by using the I-T approach ( Buckland et al., 1997 ; Burnham and Anderson, 2004 , 2014 ; Burnham et al., 2011a , b ; Sober, 2002 ). This methodology is based on the corrected Akaike Information Criterion as AICc_m ( Akaike, 1974 ; Sugiura, 1978 ; Hurvich and Tsai, 1989 ; Cavanaugh, 1997 ), and the derived measures Δ _m ( m -th Akaike differences), p_m ( m -th model probabilities), and ER_m ( m -th evidence ratios) computed for each m −th feasible model version fitted. We discarded solutions that yielded heavy-tailed 0.95 confidence intervals (0.95CI) for variance components of the pdDiag or pdSymm structures, as well as for parameters of the correlation structure, the power-of-the-mean function, or both ( Pinheiro and Bates, 2000 ).

2.6. Generalized linear mixed-effects models: second step prediction

Only one solution satisfied all criteria for model selection ( Table 2 ). Therefore, we obtained the random effects as outputs of the nlme function and computed the nonlinear geometric functions of the parameters. In sequence, we fitted the values for each cow and cpo for all BPO , which generated new random variables to be analyzed by using the GLIMMIX procedure (SAS University Edition, SAS Systems Inc., Cary, NC, USA). We checked the patterns of the Pearson residuals to evaluate the quality of fit of the model and the assumed probability density function. The model definition was as follows:

Table 2
Estimate of the negative logarithm of the likelihood function ( $-{\widehat{\ell }}_m$ ) and derived information criteria ¹ for the fitted Brody, Dijkstra, Gamma, and Wood nonlinear mixed-effects models

Model 2	Fixed	Random	Variance 3	$-{\widehat{\ell }}_m$	AICc m	Δ m	p m	ER m	n m
Brody	PO	No effect	Homog	114556.3	229166.7	9301.0	→ 0	→ ∞	27
Brody	PO	No effect	VP	113542.8	227141.6	7275.9	→ 0	→ ∞	28
Djikstra	BPO	No effect	Homog	116077.2	232240.5	12374.8	→ 0	→ ∞	43
Djikstra	BPO	No effect	VP	114708.3	229504.6	9638.9	→ 0	→ ∞	44
Gamma	BPO	No effect	Homog	129032.2	258180.6	38314.9	→ 0	→ ∞	58
Gamma	PO	α	Homog	114547.9	229137.8	9272.1	→ 0	→ ∞	21
Gamma	PO	α	VP	113854.7	227753.5	7887.8	→ 0	→ ∞	22
Gamma	PO	D ( α , λ )	Homog	113267.2	226580.4	6714.7	→ 0	→ ∞	23
Gamma	PO	D ( α , N )	Homog	113164.7	226375.4	6509.7	→ 0	→ ∞	23
Gamma	PO	D ( α , λ )	VP	112545.2	225138.5	5272.8	→ 0	→ ∞	24
Wood	B	λ	Homog	121633.9	243303.9	23438.2	→ 0	→ ∞	18
Wood	PO	λ	Homog	121281.2	242604.5	22738.8	→ 0	→ ∞	21
Wood	B	λ	VP	120137.8	240319.6	20453.9	→ 0	→ ∞	22
Wood	B	λ	VCP	120041.7	240123.4	20257.7	→ 0	→ ∞	20
Wood	B	β	Homog	115326.5	230689	10823.3	→ 0	→ ∞	18
Wood	B	D ( α , β )	Homog	114854.3	229748.6	9882.9	→ 0	→ ∞	20
Wood	PO	β	Homog	114642.8	229327.6	9461.9	→ 0	→ ∞	21
Wood	PO	No effect	VP	114539.7	229123.3	9257.6	→ 0	→ ∞	22
Wood	PO	D ( α , β )	Homog	114493.4	229032.7	9167.0	→ 0	→ ∞	23
Wood	B	β	VCP	114292.6	228625.2	8759.5	→ 0	→ ∞	20
Wood	BPO	β	Homog	114092.6	228305.4	8439.7	→ 0	→ ∞	60
Wood	PO	β	VCP	113775.9	227597.8	7732.1	→ 0	→ ∞	21
Wood	PO	D ( α , β , λ )	Homog	113265	226580.1	6714.4	→ 0	→ ∞	25
Wood	PO	D (α, λ )	VCP	113251.5	226553.1	6687.4	→ 0	→ ∞	25
Wood	BPO	β	VP	113193.8	226509.8	6644.1	→ 0	→ ∞	61
Wood	B	D ( β , λ )	Homog	113224.4	226488.9	6623.2	→ 0	→ ∞	20
Wood	B	D ( α , β , λ )	Homog	113185.9	226415.8	6550.1	→ 0	→ ∞	22
Wood	BPO	D ( α , λ )	Homog	112928	225980.3	6114.6	→ 0	→ ∞	62
Wood	PO	D ( β , λ )	Homog	112863.1	225772.3	5906.6	→ 0	→ ∞	23
Wood	BPO	D ( α , β , λ )	Homog	112686.5	225501.3	5635.6	→ 0	→ ∞	64
Wood	PO	No effect	VP	112680.4	225404.8	5539.1	→ 0	→ ∞	22
Wood	BPO	D ( β , λ )	Homog	112376.2	224876.7	5011.0	→ 0	→ ∞	62
Wood	B	D ( β , λ )	VP	112193	224428.3	4562.6	→ 0	→ ∞	21
Wood	PO	D ( β , λ )	VCP	111904.9	223859.7	3994.0	→ 0	→ ∞	25
Wood	PO	D ( β , λ )	VP	111904.9	223857.8	3992.1	→ 0	→ ∞	24
Wood	B	U ( β , λ)	Homog	111776	223596.1	3730.4	→ 0	→ ∞	22
Wood	BPO	D ( β , λ )	VP	111396.7	222919.6	3053.9	→ 0	→ ∞	63
Wood	BPO	β	CORR	110917.1	221956.4	2090.7	→ 0	→ ∞	61
Wood	B	No effect	VP	110787	221620	1754.3	→ 0	→ ∞	23
Wood	BPO	β	VCP, CORR	109883.2	219892.7	26.6	1.67 · 10−6	5.97 · 105	63
Wood	BPO	β	VP, CORR	109871	219866.1	0.0	≅ 1	≅ 1	62

• ¹
  AICcm - corrected Akaike information criterion; Δ m - Akaike difference; p m - model probability; ER m - evidence ratio; and n m - number of parameters of the m -th model combination.
• ²
  Model combination described by the fixed part of the model, i.e., Brody, Dijkstra, Gamma or Wood, followed by the parameter(s) to which random effects were ascribed with a symmetrical ( U ) or diagonal ( D ) covariance structure applied to two or more parameters, parity order ( PO ), breed ( B ), and breed × parity order ( BPO ) effects. The homogeneous variance (Homog) and residual independence correspond to the traditional assumptions for fitting a given model.
• ³
  Variance functions: VarPower (VP), VarConstPower (VCP), and correlation (CORR).

Eq. (20) $${\eta }_{ijk}=\eta +B_i+P_j+B{P}_{ij}+c_k$$

The linear predictor is η_ijk and the Greek uppercase letters represent the i -th breed ( B_i ), the j -th parity order ( P_j ), and their interaction ( BP_ij ). We assumed that variables could be y_ijk | c_k ~ Gamma( μ_ijk , Φ ) or y_ijk | c_k ~ Normal( μ_ijk , σ ²), and the random effect of cow is identically and independently distributed as c_k ~ Normal(0, ${\sigma }_c^2$ ). The Greek capital Φ is a scale parameter, and the link functions used were as follows ( Stroup, 2013 ): the identity link or η_ijk = μ_ijk was used for the Normal distribution, whereas the inverse and log links, namely ${\eta }_{ijk}={\mu }_{ijk}{}^{-1}$ or η_ijk = log ( μ_ijk ), were used for the Gamma PDF. Because the Normal distribution and the Gamma distribution produced equal residual patterns, the alternative Gamma PDF was chosen as the best solution. Therefore, the Gamma was used as a generalizing distribution for variables in the domain (0, ∞).

3. Results

We challenged the traditional assumption of independence and homoscedasticity among errors. Initially, we fitted the nonlinear models corresponding to the general mean with fixed effects of breed ( B ), parity order ( PO ), and their interaction ( BPO ), without the correlation ( equation 18 ) among repeated measures, and with a simple variance function ( equation 16 ). The gradual introduction of fixed effects resulted in an improvement of the model quality of fit ( Table 2 ). The addition of a random effect was a necessary improvement, as well as the introduction of a first-order autoregressive correlation ( equation 19 or corCAR1 function, parameter ϕ ) and the power-of-the-mean variance function ( equation 17 , parameter ψ ). Nonetheless, the introduction of a single random effect to parameter β of equation 1 was the best solution after an exhaustive and careful choice among several possible and feasible model versions, given the data ( Tables 2 and 3 ). Therefore, the best model in the set of all models studied was formed by combining equation 1 fitted with 62 parameters related to the BPO interactions as the fixed effects, the variance function described by equation 17 , the corCAR1 function, and a random effect associated with the shape parameter β . We also presented the standard errors of the variance components of the random effects and their respective 0.95CI related to parameter β , namely, the standard deviations related to the cows ( σ_cow,β ), the intersecting cow × parity order ( σ_cpo,β ), and the residual error ( σ ). There is one parameter estimated from equation 17 and another parameter from equation 19 , namely ψ and ϕ , respectively ( Table 3 ).

Table 3
Fixed parameter estimates and respective standard errors (SE) of the chosen model according to the information-theoretic approach

Effect 1	αtβ exp(− λt ) 2	SEα	SEβ	SEλ
BPO 11	15.5 t 0.023 exp(−0.0021 t )	3.53	0.0714	9.70 · 10−4
BPO 12	16.2 t 0.094 exp(−0.0038 t )	6.79	0.1304	1.67 · 10−3
BPO 21	17.4 t 0.120 exp(−0.0022 t )	0.31	0.0051	5.00 · 10−5
BPO 22	22.5 t 0.111 exp(−0.0028 t )	0.41	0.0055	6.20 · 10−5
BPO 23	24.2 t 0.104 exp(−0.0028 t )	0.66	0.0083	9.40 · 10−5
BPO 24	26.7 t 0.080 exp(−0.0027 t )	1.23	0.0140	1.58 · 10−4
BPO 25	26.0 t 0.090 exp(−0.0029 t )	1.66	0.0194	2.32 · 10−4
BPO 26	26.8 t 0.074 exp(−0.0024 t )	2.58	0.0283	3.43 · 10−4
BPO 31	16.1 t 0.037 exp(−0.0017 t )	0.86	0.0153	1.55 · 10−4
BPO 32	16.2 t 0.079 exp(−0.0024 t )	0.96	0.0185	2.41 · 10−4
BPO 33	14.4 t 0.149 exp(−0.0038 t )	1.96	0.0459	7.83 · 10−4
BPO 34	16.7 t 0.046 exp(−0.0021 t )	9.74	0.2171	4.74 · 10−3
BPO 41	19.0 t −0.041 exp(−0.0010 t )	1.84	0.0289	3.18 · 10−4
BPO 42	18.2 t 0.053 exp(−0.0024 t )	1.84	0.0323	4.25 · 10−4
BPO 43	23.3 t 0.001 exp(−0.0025 t )	5.80	0.0801	1.18 · 10−3
BPO 44	13.2 t 0.212 exp(−0.0050 t )	7.24	0.1702	2.05 · 10−3
BPO 51	16.5 t 0.096 exp(−0.0025 t )	2.53	0.0444	5.03 · 10−4
BPO 52	15.2 t 0.178 exp(−0.0041 t )	2.99	0.0609	7.75 · 10−4
BPO 53	18.0 t 0.154 exp(−0.0043 t )	6.37	0.1135	1.58 · 10−3

• ¹
  Random effect standard errors (0.95 confidence intervals within parentheses) related to parameter β are as follows: cow effect, σ cowβ = 0.0561 (0.0542, 0.0580), cow × parity order effect, σ cpoβ = 0.0131 (0.0098, 0.0176), and residual standard error, σ = 0.94 (0.91, 0.97).
• ²
  Estimates of the first-order autoregressive correlation coefficient and the power-of-the-mean variance function (0.95 confidence intervals within parentheses) are ϕ = 0.67 (0.66, 0.68) and ψ = 0.57 (0.56, 0.58), respectively.

Equation 1 presented convergence in all cases. In contrast, the other models sometimes failed at convergence or yielded non-positive definite Hessians for variance-covariance matrices. Other problems observed were unreliable (heavy-tailed) confidence intervals for variance-covariance parameters, and we discarded those solutions as recommended by Pinheiro and Bates (2000) . Equations 12 and 13 yielded poor performances if compared with the best solution for equation 1 . Only model versions with fixed effects were feasible for those models.

The nonlinear geometric functions of the fixed parameters ( α , β , and λ ) of the best model ( equation 1 ) and respective standard errors were estimated for each BPO ( Table 3 ). The estimate obtained for parameter β of BPO ₄₁ was negative. The sign of the estimate of parameter β determines the shape of the lactation function, so that the curve was concave for β ≤ 0 ( Figure 2 , panel e) or presented a typical convex profile for 0 < β ≤ 1 ( Figure 2 , panels a, b, c, d, and f). A group of 385 cows presented a final negative random value for β ; therefore, 6369 cows presented typical convex-shaped lactation profiles. The curves fitted with fixed parameter β ≤ 0 did not present peak milk yields as observed for convex-shaped profiles, and we found those results for the first and third parity orders of the crossbreds, namely BPO ₄₁ and BPO ₄₃. Possibly, the low number of records and cows of the third and fourth parity orders limit the strength of the inferences taken.

Figure 2
Examples of expected lactation profiles for some breed ( i ) × parity order ( j ) interactions ( BPO_ij ).

Other attributes computed from the two-step prediction, such as ${\mu }_{t_p}$ , κ_α , κ_d , t_p , etc., are presented with standard errors of the means and respective 0.95CI for all BPO combinations ( Tables 4 and 5 ). Sixth parity order cows of the Holstein breed produced more milk per lactation than any other parity order for Holsteins or breeds, as one can check by comparing the confidence intervals reported. Therefore, Holsteins produced more milk than any other breed in the Southwest Paraná state. Fourth-, fifth-, and sixth-order Holstein cows presented the same peak milk yield, as demonstrated by overlapping confidence intervals ( Table 4 ). The combination of the two-step statistical tools allowed the identification of these attributes for ranking animals within breeds and groups.

Table 4
Attributes ¹ of the lactation records obtained from the two-step estimation based on Wood's equation: average milk production ( ${\overline{\mu }}_{0,305}$ , kg · d⁻¹), time to peak milk yield ( t _p , w), peak milk yield ( ${\mu }_{t_p}$ , kg · d⁻¹), cumulative or total milk production ( TMμ_t , kg), and their respective standard errors (SE) and lower (L) and upper (U) 0.99 confidence limits

BPO 2	${\overline{\mu }}_{0,305}$	$\pm {\text{SE}}_{\overline{\mu }}$	${\text{L}}_{\overline{\mu }}$	${\text{U}}_{\overline{\mu }}$	t p	$\pm {\text{SE}}_{t_p}$	${\text{L}}_{t_p}$	${\text{U}}_{t_p}$	${\mu }_{t_p}$	$\pm {\text{SE}}_{{\mu }_{t_p}}$	${\text{L}}_{{\mu }_{t_p}}$	${\text{U}}_{{\mu }_{t_p}}$	TMμ t	±SETM	LTM	UTM
BPO 11	12.1	0.33	11.3	13.0	2.2	0.35	1.5	3.3	17.2	0.49	16.0	18.5	3869	289.1	3124	4614
BPO 12	14.2	0.48	13.1	15.6	2.6	0.51	1.6	4.3	21.3	0.87	19.2	23.8	4526	298.7	3756	5295
BPO 21	21.7	0.08	21.4	21.9	7.2	0.08	7.0	7.4	25.0	0.07	24.9	25.2	6880	23.6	6820	6941
BPO 22	24.3	0.11	24.0	24.6	5.2	0.06	5.0	5.3	29.8	0.10	29.6	30.1	7862	23.8	7800	7923
BPO 23	25.2	0.12	24.9	25.5	4.8	0.06	4.6	4.9	31.2	0.11	30.9	31.4	8197	24.3	8134	8259
BPO 24	25.2	0.12	24.9	25.5	3.9	0.07	3.7	4.0	31.5	0.13	31.2	31.9	8218	25.6	8152	8284
BPO 25	25.0	0.13	24.7	25.3	4.1	0.10	3.8	4.4	31.6	0.15	31.2	32.0	8149	28.9	8075	8223
BPO 26	25.5	0.16	25.1	26.0	4.3	0.19	3.9	4.9	31.4	0.22	30.8	32.0	8410	40.5	8305	8514
BPO 31	14.6	0.09	14.4	14.9	3.2	0.08	2.9	3.4	18.2	0.09	17.9	18.4	4643	55.8	4499	4787
BPO 32	16.1	0.11	15.8	16.4	4.3	0.13	4.0	4.7	20.4	0.13	20.1	20.7	5138	56.7	4992	5284
BPO 33	16.2	0.13	15.9	16.5	5.7	0.30	5.0	6.5	21.7	0.20	21.2	22.2	5196	62.9	5034	5358
BPO 34	14.7	0.22	14.2	15.4	2.3	0.31	1.6	3.2	19.0	0.48	17.8	20.3	4783	98.5	4529	5037
BPO 41	13.3	0.14	13.0	13.7	1.7	0.13	1.3	2.0	18.7	0.28	18.0	19.4	4240	102.2	3977	4504
BPO 42	15.8	0.20	15.3	16.3	3.6	0.28	2.9	4.4	22.9	0.41	21.9	24.0	5155	104.1	4887	5423
BPO 43	15.8	0.22	15.2	16.4	1.8	0.22	1.3	2.5	24.7	0.64	23.1	26.4	5140	112.1	4851	5429
BPO 44	16.5	0.31	15.7	17.3	8.8	1.66	5.4	14.3	27.0	0.93	24.8	29.6	5506	142.6	5139	5874
BPO 51	17.5	0.55	16.2	19.1	5.2	0.55	4.0	6.8	21.8	0.52	20.5	23.2	5610	233.0	5010	6211
BPO 52	18.7	0.63	17.2	20.4	6.1	0.70	4.6	8.3	25.0	0.71	23.3	27.0	6073	235.4	5466	6679
BPO 53	19.1	0.67	17.5	21.0	5.1	0.72	3.6	7.4	26.4	0.85	24.4	28.8	6286	243.0	5660	6912

• ¹
  Additional parameters estimated as nonlinear functions of Wood's equation parameters.
• ²
  BPOij means i -th breed versus j -th parity order interaction, ∀ i , j: Girolando ( i = 1), Holstein ( i = 2), Jersey ( i = 3), Crossbreed ( i = 4), Brown Swiss ( i = 5); and parity orders are j = 1st, 2nd, …, 6th.

Table 5
Attributes ¹ of the lactation records obtained from the two-step estimation based on Wood's equation: specific ascending rate of milk production until peak milk yield ( κ _α , w⁻¹), time at the inflection point ( t _i , w), specific rate of milk production decline post peak milk yield ( κ _d , w⁻¹), their respective standard errors (SE), and lower (L) and upper (U) 0.99 confidence limits

BPO 2	κ α	$\pm {\text{SE}}_{{\kappa }_{\alpha }}$	${\text{L}}_{{}_{{\kappa }_{{}_{\alpha }}}}$	${\text{U}}_{{}_{{\kappa }_{{}_{\alpha }}}}$	t i	$\pm {\text{SE}}_{t_i}$	${\text{L}}_{t_i}$	${\text{U}}_{t_i}$	κ d	$\pm {\text{SE}}_{{\kappa }_d}$	${\text{L}}_{{\kappa }_d}$	${\text{U}}_{{\kappa }_d}$
BPO 11	0.177	0.0131	0.149	0.219	14.1	1.32	11.1	17.9	1.16	2.790	1.10	1.24
BPO 12	0.200	0.0180	0.163	0.261	12.9	1.48	9.6	17.4	1.73	6.820	1.57	1.92
BPO 21	0.117	0.0004	0.116	0.118	28.8	0.18	28.4	29.3	0.96	0.130	0.96	0.96
BPO 22	0.152	0.0007	0.150	0.154	21.5	0.14	21.1	21.8	1.24	0.220	1.24	1.25
BPO 23	0.157	0.0008	0.155	0.159	20.5	0.16	20.1	20.9	1.26	0.250	1.25	1.26
BPO 24	0.172	0.0010	0.169	0.174	18.1	0.18	17.7	18.6	1.27	0.300	1.26	1.28
BPO 25	0.176	0.0014	0.172	0.179	18.3	0.26	17.7	19.0	1.34	0.440	1.32	1.35
BPO 26	0.156	0.0020	0.151	0.161	20.3	0.51	19.0	21.6	1.12	0.580	1.11	1.14
BPO 31	0.133	0.0013	0.130	0.136	19.5	0.31	18.7	20.3	0.86	0.270	0.86	0.87
BPO 32	0.149	0.0017	0.145	0.154	20.7	0.37	19.7	21.6	1.11	0.480	1.10	1.12
BPO 33	0.167	0.0027	0.160	0.174	20.7	0.61	19.2	22.4	1.54	1.230	1.51	1.57
BPO 34	0.157	0.0047	0.146	0.170	15.7	1.23	12.8	19.2	1.08	1.550	1.04	1.12
BPO 41	0.098	0.0021	0.093	0.104	16.5	0.78	14.6	18.7	0.57	0.360	0.56	0.58
BPO 42	0.148	0.0045	0.137	0.161	18.8	0.87	16.7	21.1	1.12	1.310	1.09	1.16
BPO 43	0.199	0.0098	0.176	0.228	12.9	0.92	10.7	15.5	1.26	2.240	1.21	1.32
BPO 44	0.146	0.0061	0.132	0.164	21.5	2.11	16.7	27.7	1.70	4.400	1.59	1.82
BPO 51	0.146	0.0060	0.132	0.164	22.5	1.42	19.2	26.5	1.13	1.770	1.08	1.18
BPO 52	0.159	0.0073	0.142	0.181	21.0	1.42	17.7	25.0	1.64	3.820	1.54	1.74
BPO 53	0.181	0.0100	0.158	0.211	17.4	1.42	14.2	21.5	1.76	4.710	1.64	1.89

• ¹
  Additional parameters estimated as nonlinear functions of Wood's equation parameters. To be presented, actual estimates of κ d , Lκd , and Uκd were multiplied by 10², whereas actual SEκd estimates were multiplied by 10⁴.
• ²
  BPOij means i -th breed versus j -th parity order interaction, ∀ i , j: Girolando ( i = 1), Holstein ( i = 2), Jersey ( i = 3), Crossbreed ( i = 4), Brown Swiss ( i = 5); and parity order are j = 1st, 2nd, …, 6th.

Older Holstein cows reached t_p earlier than their first- and second-order counterparts did ( Table 4 ). Brown Swiss cows followed the same Holstein patterns, but stayed in second place as milk producers. Jersey and crossbred cows did not follow the same patterns presented by Holsteins and Brown Swiss cows. Second- and third-order Jersey cows produced more milk than the ones of the first parity order. Girolando cows were low in rank, but the dataset contained only first and second parity order cows ( Tables 1 and 4 ).

Holsteins presented an increasing pattern for κ_α from the first to the fifth parity order, but the rate reduced at the sixth parity order ( Table 5 ). The same happened for Jerseys and crossbreds from first to the fourth lactation, whereas Brown Swiss cows showed an increase from first to the third lactation. The second-order Girolando cows presented an ascending rate of milk production faster than their first-order counterparts did. Despite the β ≤ 0 estimates for BPO ₄₁ and BPO ₄₃, the convex profiles generated from positive random effects over the fixed β allowed the estimation of κ_α ( Table 5 ).

Holsteins and Brown Swiss cows of the first parity order generally reached the time at the inflection point later than their older counterparts did. The exception was the sixth parity order Holsteins ( Table 5 ). However, this pattern was not followed by the other breeds and the crossbred cows. In fact, Jersey cows increased t_i as they got second and third parity orders, but fourth-order Jersey cows reached the inflection point earlier.

The negative fractional descending rate was higher for first parity order cows, irrespective of breed or crossbreed ( Table 5 ). Older cows generally presented lower descending rates; therefore, they did not sustain the peak milk yield as first parity order cows did. There was a negative association between κ_α and the negative ( ${\kappa }_d\left(\widehat{\rho }=-0.896\right)$ , P<0.001). Given the values of the SE reported, we can depict that κ_α and κ_d estimates for the fixed effects were precise. The relative variation for the ascending rate ( $100{\text{SE}}_{\kappa \alpha }/{\kappa }_{\alpha }$ ) ranged from 0.33 to 8.97%, whereas the relative variation for κ_d ranged from 0.14 to 3.95%. On the other hand, the estimates for λ presented a greater relative variation for the fixed effects, i.e., from 2.18% to 221%. The introduction of a random effect on parameter λ was ineffective to reduce the AICc_m substantially if compared with the random effect ascribed to β ( Table 2 ). Therefore, only positive values occurred for λ as population estimates ( Table 3 ). The ascending and descending rates are nonlinear functions of the model parameters, and the integrated information resulted in precise estimates, which favored the direct comparisons; isolated parameters were less precise ( Tables 3 and 5 ).

For comparison purposes, we regressed the estimated least-squares means for r_d (w⁻¹) over the estimated least-squares means for κ_d (w⁻¹, Eq. 8 ) in their original scales, i.e., both rates are characteristic of the descending phase after peak milk yield and are both negative. Because of ${\widehat{r}}_d$ and ${\widehat{k}}_d$ domains (estimates < 0), we assumed r_d normally distributed, and we estimated the simple linear regression by least squares (PROC REG, SAS University Edition): ${\widehat{r}}_d=3.5\cdot {10}^{-3}+1.5\cdot {\widehat{k}}_d$ . Therefore, for the proper estimation of ${\widehat{r}}_d$ one must use $-{\widehat{k}}_d$ in the preceding regression ( Table 5 ). The SE of the intercept was 4.11 · 10⁻⁴, the SE of the slope was 3.19 · 10⁻², the standard error of the regression was 4.19 · 10⁻⁴, and the R ² = 0.9925. The estimated values for r_d and κ_d were precise and close to each other. Regardless of breed, cows of first parity order presented smaller arbitrarily positive κ_d estimates in the first half of the down ranked least-squares means contained in Table 5 .

The relationships among parameters α , β , and λ with, for instance, κ_α , κ_d , t_p , ${\mu }_{t_p}$ , and TM μ_t are mathematical. Therefore, the functional relationships established among those parameters explain, by and large, the correlation estimates among those same parameters. Nonetheless, the Pearson correlation for $TM{\mu }_t\times {\mu }_{t_p}$ in our study amounted to 0.982 (P<0.001) regardless of breed and parity order. This was done for those cows ( n_c = 6369) with an overall random estimate for β ∈ (0, 1). Nonetheless, the random β estimates for all 6754 cows presented a positive correlation with TMμ_t , namely 0.873 (P<0.001). Therefore, generally, even for cows with atypical concave profiles, the greater the random estimates of β , the greater the predictions for TMμ_t . In our data, some examples of the Pearson correlation coefficients were: −0.519 (P<0.001) for TMμ_t × κ_α, 0.035 (P = 0.005) for TMμ_t × κ_α, −0.828 (P<0.001) for κ_α × κ_d , and −0.896 (P<0.001) for κ_α × t_p .

4. Discussion

The technique of the nonlinear mixed-effects models allows the ranking of dairy cows within each breed and each parity order (ranks not shown). Therefore, one can replace the unstructured G matrix and ascribe a phenotypic value by associating to it genetic merits by the two-step estimation technique with a genomic relationship matrix in the linear mixed model ( Silva et al., 2017 ; Soares et al., 2017 ). An extension to the generalized linear mixed-effects model is possible ( Littell et al., 2006 ; Vonesh, 2012 ; Stroup, 2013 ). Nonetheless, we remind the reader that we fitted the random effects of cow ( σ_cow,β ) and cow × parity order ( σ_cpo,β ) successfully only to the shape parameter β , and this result may explain the fitting flexibility of Wood's equation over the other equations studied here. In addition, the combined powers of equations 17 and 19 accounted for variations attributable to either scale and correlations among repeated measures ( Tables 2 and 3 ). If not discounted, those variations can be understood as pseudorandom variations that otherwise would inflate the residual variation ( Vieira et al., 2018 ). The nlme function of R accounts for all these issues and its main advantage is to integrate all information in a single variance-covariance matrix. The overall benefits are an improved precision of parametric estimates and information loss minimization during the quantitative interpretation of nonlinear phenomena with mathematical models ( Pinheiro and Bates, 2000 ).

4.1. Nonlinear models and their quality of fit

The Akaike criterion ( Akaike, 1974 ) corrected for small samples ( Sugiura, 1978 ; Hurvich and Tsai, 1989 ; Cavanaugh, 1997 ), namely AICc_m , has been reported as one of the most important measures to evaluate the predictive power of mathematical models in a multiple hypotheses framework ( Buckland et al., 1997 ; Sober, 2002 ; Burnham and Anderson, 2004 ). The AICc_m and its derived measures constitute the I-T approach, which allows the comparison of the quality of fit of different models and, whenever necessary, establishes parameters for model averaging. If one takes the smallest AICc_m value to choose a model over the others, does an incomplete procedure within the I-T framework ( Burnham and Anderson, 2004 ; Vieira et al., 2018 ; 2020 ). The smallest AICc_m value indicates the lower information loss due to the fit of a given model among the feasible models evaluated. However, if one assumes that a given model is the best solution and ignores the others, depending on the p_m for the chosen model, its uncertainty can be large. This is why we assumed that a low uncertainty for a given model is 1 − p_m < 0.10 ( Burnham and Anderson, 2004 ; Vieira et al., 2018 ).

The frequentist statistical tools available nowadays allow the association of random factors to the fixed parameters of nonlinear models and the fit of the resultant models by maximum or restricted maximum likelihood. In addition, one can challenge the traditional assumptions of homoscedasticity and independence of the errors with these tools based on the theories of nonlinear and generalized linear mixed-effects models ( Pinheiro and Bates, 2000 ; Vonesh, 2012 ; Stroup, 2013 ). In our study, the use of the available statistical tools (nlme from R and GLIMMIX of SAS) provided considerable improvements in model fitting and prediction. Based on the I-T approach, we observed that Wood's equation was superior in terms of quality of fit and the best choice to represent the diversity in shape and form of the lactation profiles, given the data. Nonetheless, many studies use the R ² to compare models, but this measure is not sensitive to differences among nonlinear models, which make comparisons useless because of the difficulty to demonstrate differences. Therefore, models that actually differ based on I-T may be considered equal if one bases on R ² or, alternatively, its adjusted version for model selection ( Spiess and Neumeyer, 2010 ; Hossein-Zadeh, 2016 ).

4.2. Lactation curve models

In this study, we challenged the ability of some models to mimic actual lactation records over time. Although the abundant dataset (37810 time-records) might have favored both accuracy and precision of the estimates, considerable variation may have arisen due to different ethnic types, parity orders, low number of cows for some breed-parity order interactions, and incomplete or irregular lactation records to characterize an entire lactation (0-305 days). These sources of variation might have constrained some of the results and increased the difficulty of obtaining valid, convergent, and robust solutions for some models. Nonetheless, the dataset offered the opportunity for Wood's model to exhibit its flexible nature to fit both convex and concave lactation profiles. Concave profiles are rather frequent than one might expect; these concave profiles are characteristic of animals that did not present peak milk yield, and equation 1 mimicked those profiles when other models ( equations 12 - 14 ) failed. Within the framework of nonlinear mixed-effects models, Wood's model even was able to fit lactation profiles with incomplete records, which are needed to characterize an entire lactation trend ( Table 3 and Figure 2 , panel e).

Several authors faced difficulties with fitting models to atypical lactation profiles, that is, those profiles that deviate from the standard lactation convex-shaped trend ( Rekik and Ben Gara, 2004 ; Macciotta et al., 2005 ; López et al., 2015 ). Therefore, the greater flexibility of equation 1 provides the necessary accuracy and precision to the selection of cows based on ranking their performances as parametric attributes of the lactation profiles such as peak milk yield, mean milk production, ascending and descending rates, and more importantly, the total milk production per lactation.

There was a positive agreement between random β and total milk production estimates. Therefore, the use of the peak milk yield as a parameter for selection may result in the exclusion of cows that may produce high quantities of milk, even though the shape parameter β ≤ 0 for a given breed. Our findings favored the use of all lactation profiles from all cows. They avoided an equivocated exclusion of animals before the fit of the model, i.e., all animals were evaluated regardless of the occurrence of a peak milk yield or not ( Rekik and Ben Gara, 2004 ). The poor performance of the other models ( equations 12 − 14 ) may be explained by the no prior arbitrary exclusion of concave lactation records, as well as incomplete and irregular lactation profiles from our dataset. The advantage of the nonlinear mixed-effects models relies exactly upon the joint analysis of all lactation records employing a matrix of random effects, which encompasses the variation from the random sample of cows about the fixed parameters.

There are some reports that present a poor performance for equation 1 if compared with equations 12 and 13 . However, they used the nonlinear least-squares method of estimation to fit the models, did not consider random factors over the fixed parameters, nor did challenged the independence of residuals and homogeneity of variances ( Brody et al., 1924 ; Dijkstra et al., 1997 ; Hossein-Zadeh, 2016 ). Generally, convex-shaped lactation records are selected because they resemble the typical, standard lactation profile ( Morant and Gnanasakthy, 1989 ; Dijkstra et al., 1997 ; Macciotta et al., 2005 ). Nonetheless, the final shape and form of the nonlinear model depend on its parametric estimates ( Wood, 1967 ; Grossman and Koops, 1988 ; Morant and Gnanasakthy, 1989 ; Dijkstra et al., 1997 ; Macciotta et al., 2005 ). The sign and magnitude of the β estimate determine if the curve is flatter, presents a sharp peak, or decline monotonically at the onset of lactation at parturition ( Congleton Jr. and Everett, 1980 ; Morant and Gnanasakthy, 1989 ). In addition, it is important to note that the inflection point in the descending phase is part of a typical convex-shaped lactation profile ( Druet et al., 2003 ).

Congleton Jr. and Everett (1980) and Macciotta et al. (2006) believed that the lack of points at the ascending phase do not allow the proper characterization of the typical convex-shaped lactation profiles with a characteristic peak milk yield. Nonetheless, our dataset contains concave-shaped lactation records with several time points, as the better-characterized convex-shaped profiles do. This is why the accommodation of the random factor to parameter β resulted in the superiority of equation 1 to represent the data, because final random estimates can be either 0 < β ≤ 1 and β ≤ 0. Because of the factors mentioned previously, we believe that those variations occurred when cows presented the peak milk yield on the day of parturition or nearly after parturition. These types of lactation profiles occurred for groups BPO ₄₁ and BPO ₄₃, which might be associated with the poor dairy temperament of Zebu breeds and its crossbreeds with dairy breeds ( Bangar and Verma, 2017 ). Sometimes, even cows that exhibited a typical convex-shaped lactation profile in a previous lactation fail to present a peak milk yield in a further lactation, because of environmental factors such as metabolic disorders, mastitis, other diseases, and nutritional or other management errors ( Wood, 1968 , 1970 , 1972 , 1976 , 1980 ; Macciotta et al., 2005 ; Hossein-Zadeh, 2016 ; Ahmed et al., 2019 ). Nevertheless, the cows that presented typical convex-shaped lactation profiles were those that produced more milk in the entire course of lactation than their counterparts did.

4.3. Geometric attributes of Wood's model

The ascending phase of equation 1 is represented by κ_α and can be associated with the proliferation of secretory cells of the epithelial parenchyma in the mammary gland of cows ( Capuco et al., 2001 ), dairy does ( Knight and Peaker, 1984 ), and mice and rats ( Dijkstra et al., 1997 ). The DNA content of the mammary gland is richer in the ascending phase, which results in greater RNA transcription and the consequent translation into enzymes responsible for the synthesis of milk components such as galactosyl transferase, fatty acid synthase, and acetyl-CoA carboxylase ( Knight and Peaker, 1984 ; Boutinaud et al., 2004 ). The continuous milk production increase may be associated with the net number of secretory cells that result from the processes of cell mitosis or differentiation ( Capuco et al., 2001 ; Dijkstra et al., 1997 , 2010 ). Homeorhetic mechanisms orchestrate changes in the mammary gland and other metabolic organs of the animals to provide the machinery and nutrients for milk synthesis by the secretory tissue. The net results are the changes observed for milk yields in the course of the lactation, including changes in other variables that dairy animals go through during the dry and transition periods ( Bauman and Currie, 1980 ; Bauman et al., 1999 ; Bauman, 2000 ).

The smaller estimates for κ_α of primiparous cows within each breed observed in our study may be explained by the incomplete development of their mammary glands. First parity order cows are younger and still grow after the first calving; therefore, there is a competition for nutrients between the processes of growth, gestation, and the subsequent lactation, which delays the growth of the mammary gland and its full secretory potential observed in a mature cow ( Bauman and Currie, 1980 ; Bauman et al., 1999 ; Bauman, 2000 ). Indeed, the mammary glands of cows ( Capuco et al., 2001 ) and dairy does ( Knight and Peaker, 1984 ) continue to grow until and after the first parturition. The volume of secretory tissue in the mammary gland of second parity order dairy does is greater than that of the first parity order does, but the greater differentiation in the quantity of secretory parenchyma occurs during pregnancy ( Fowler et al., 1990 ; Knight and Wilde, 1993 ). Those biological facts may explain the greater κ_α estimates for multiparous cows and, consequently, their greater peak milk yields estimated in our study.

The peak milk yield can be used as a parametric reference for phenotypic selection of cows because of the good positive agreement between total milk production and peak milk yield in dairy cattle ( Grossman and Koops, 1988 ; Tekerli et al., 2000 ; Hossein-Zadeh, 2014 ; López et al., 2015 ) and buffaloes ( Hossein-Zadeh, 2016 ).

The descending phase occurs soon after the peak milk yield ( Figure 1 ) ( Wood, 1967 , 1976 ; Knight and Peaker, 1984 ; Dijkstra et al., 1997 , 2010 ; López et al., 2015 ). We used the inflection point at the descending phase to obtain an estimator for κ_d . This specific constant rate can be an index of persistency: the greater the absolute values, the lower the lactation persistence, i.e., the sharper the decline after peak milk yield. Nonetheless, those cows that presented lactation profiles with greater κ_α estimates also presented greater arbitrarily positive κ_d estimates (smaller persistency). We observed a negative Pearson correlation between these two specific rates, which clearly demonstrated the relationship. Wood (1969) observed that as cows grow older, they start with greater milk yield records at parturition, but their rate of decline in milk production reduces as lactation advances post-peak milk yield (persistency reduces = larger absolute κ_d ). Our results confirm that primiparous cows present flatter lactational profiles, which means a greater persistency if compared with their older counterparts ( Grossman and Koops, 1988 ; Bauman et al., 1999 ; Cobuci and Costa, 2012 ; Hossein-Zadeh, 2014 , 2016 ; López et al., 2015 ). Therefore, the geometrization of equation 1 yielded parameter estimates that agreed with the biological trends reported in the literature.

The milk production decline post-peak yield may be associated to the reduction in the numbers of secretory cells in the mammary glands of cows ( Capuco et al., 1997 ; Cobuci and Costa, 2012 ), dairy does ( Knight and Peaker, 1984 ; Quarrie et al., 1994 ; Dijkstra et al., 1997 ), and small rodents ( Walker et al., 1989 ; Strange et al., 1992 ; Dijkstra et al., 1997 ). The number of epithelial secretory cells recorded at the peak milk yield can be reduced by 40% until the end of lactation in cows ( Capuco et al., 2001 ) and by 30 ( Knight and Peaker, 1984 ) to 40% in dairy does ( Fowler et al., 1990 ). Capuco et al. (2001) reported that the number of secretory cells results from the balance between cell proliferation and death, i.e., the quantity of secretory tissue reduces as the cell death rate overcomes the cell proliferation rate. This resultant cell loss occurs by some sort of a programmed mechanism for cell death also called cell apoptosis ( Walker et al., 1989 ; Strange et al., 1992 ; Quarrie et al., 1994 ; Stefanon et al., 2002 ). The first report about apoptosis was for no lactating mice and rats during the period of mammary gland involution ( Walker et al., 1989 ; Strange et al., 1992 ). Nonetheless, Knight and Peaker (1984) detected the reduction in the number of cells of the secretory tissue after the eighth week of lactation in Saanen does. However, the association of this reduction in the number of secretory cells to the phenomenon of lactational apoptosis in ruminants was introduced by Quarrie et al. (1994) . Therefore, from the estimation of κ_α and κ_d , we can infer quantitatively about the anatomic and physiological changes during pregnancy and mainly after parturition. We can also speculate, based on our correlated κ_α and κ_d estimates, that the faster the cell proliferation occurs, the sooner the process of cell apoptosis prevails in the secretory tissue of the mammary gland. In addition, we can use those parametric estimates to predict the impact of long-term environmental constraints over milk production trends ( Wood, 1972 ), e.g., the attenuation or worsening of the effects of the oxidative stress over the secretory capacity of the mammary gland, and the consequent result in milk production ( Bernabucci et al., 2002 , 2005 ). Other solutions may be of use for short-term effects that dramatically alter the lactation trends ( Ahmed et al., 2019 ).

We can make a direct algebraic comparison, for instance, between equations 1 and 13 . If we take the root for t_i > t_p ( equation 21 ) among other possible root solutions to d ² μ_t / dt ² = 0 for μ_t from equation 13 ^1,2 , the time at the inflection point would be as follows:

Eq. (21) $$t_i={\kappa }_3{}^{-1}\log \left(0.5{\kappa }_4{}^{-2}{\theta }_T\left({\kappa }_3+2{\kappa }_4\pm \sqrt{{\kappa }_3\left({\kappa }_3+4{\kappa }_4\right)}\right)\right)$$

Thence, we can apply the same rationale for κ_d as described by equations 7b and 8 , with the properties described by Dijkstra et al. (1997) . In sequence, we can also obtain an equivalent specific ratio as described by equation 11 , because

Eq. (22) $${\kappa }_{\alpha }=\left({\mu }_{t_p}-{\mu }_0\right){\left({\displaystyle \int {}_0^{t_p}}{\mu }_0\exp \left({\theta }_T{\kappa }_3{}^{-1}\left(1-\exp \left(-{\kappa }_{3}t\right)\right)-{\kappa }_{4}t\right)dt\right)}^{-1}$$

In equation 22 , $t_p={\kappa }_3^{-1}\log \left({\kappa }_4^{-1}{\theta }_T\right)$ , ${\mu }_{t_p}={\mu }_0{\left({\kappa }_4/{\theta }_T\right)}^{{\kappa }_4/{\kappa }_3}\exp \left(k_3^{-1}\left({\theta }_T-{\kappa }_4\right)\right)$ , and the definite integral can be solved numerically. Unfortunately, the model described by equation 13 did not fit our dataset with the same flexibility as Wood's model did, and the comparison was only on a theoretical basis.

5. Conclusions

Based on the I-T approach, we demonstrated the superior quality of fit of Wood's equation combined with a power-of-the-mean function to account for heterogeneous variances, with a first-order autoregressive correlation structure, and by introducing the random cow and cow × parity order effects over the shape parameter ( β ). This choice provides the best-suited model for mimicking lactation records quantitatively from a large dataset of cows from different breeds and with several parity orders, given the data, particularly if the dataset contains convex- and concave-shaped lactation profiles. It also provides the possibility of extracting nonlinear, geometric functions based on Wood's equation parameters via a two-step prediction procedure to yield both accurate and precise estimates for ranking cows by predicting random effects over nonlinear functions of parameters. In this regard, the two-step prediction of the variables based on the geometric functions we deduced from Wood's equation ( κ_α , κ_d , ${\overline{\mu }}_{0,305}$ , t_i and TMμ_t ), together with previous ones already known ( t_p and ${\mu }_{t_p}$ ), may be of use as selection measures about cow performance and for quantitative associations or predictions for decision making at the farm level.

• 1
  $d{\mu }_t/dt={\mu }_0\left({\theta }_{\text{T}}\exp \left(-{\kappa }_{3}t\right)-{\kappa }_4\right)\exp \left({\kappa }_3^{-1}{\theta }_{\text{T}}\left(1-\exp \left(-{\kappa }_{3}t\right)\right)-k_{4}t\right)$ , for μ_t from equation 13 .
• 2
  $d^2{\mu }_t/d{t}^2={\mu }_0{\left({\theta }_{\text{T}}\exp \left(-{\kappa }_{3}t\right)-{\kappa }_4\right)}^2\exp \left({\kappa }_3^{-1}{\theta }_{\text{T}}\left(1-\exp \left(-{\kappa }_{3}t\right)\right)-k_{4}t\right)-{\mu }_0{\theta }_{\text{T}}{\kappa }_3\exp \left({\kappa }_3^{-1}{\theta }_{\text{T}}\left(1-\exp \left(-{\kappa }_{3}t\right)\right)-\left({\kappa }_3+{\kappa }_4\right)t\right)$ , for μ_t from equation 13 .

Acknowledgments

We wish to thank Prof. Fabyano Fonseca e Silva (Departamento de Zootecnia of Universidade Federal de Viçosa) for his critical reading of the manuscript. R.A.M. Vieira thanks the grants provided by the Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq (303974/2019-2) and Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro - FAPERJ (E-26/202.544/2019). J.G. Oliveira thanks the fellowships provided by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code 001, and CNPq (141916/2018-4). Finally, we dedicate, at heart, this paper to the memory of Prof. Douglas Sampaio Henrique (UTFPR): disciple, mentor, and above all, a true friend.

References

• Ahmed, B. A.; Laurence, P.; Pierre, G. and Olivier, M. 2019. Lactation curve model with explicit representation of perturbations as a phenotyping tool for dairy livestock precision farming. bioRxiv 661249. https://doi.org/10.1101/661249
  » https://doi.org/10.1101/661249
• Akaike, H. 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control 19:716-723. https://doi.org/10.1109/TAC.1974.1100705
  » https://doi.org/10.1109/TAC.1974.1100705
• Albertini, T. Z.; Medeiros, S. R.; Torres Júnior, R. A. A.; Zocchi, S. S.; Oltjen, J. W.; Strathe, A. B. and Lanna, D. P. D. 2012. A methodological approach to estimate the lactation curve and net energy and protein requirements of beef cows using nonlinear mixed-effects modeling. Journal of Animal Science 90:3867-3878. https://doi.org/10.2527/jas.2010-3540
  » https://doi.org/10.2527/jas.2010-3540
• Araujo, R. P.; Vieira, R. A. M.; Rocha, N. S.; Abreu, M. L. C.; Glória, L. S.; Rohem Júnior, N. M. and Fernandes, A. M. 2015. Long-term growth of body, body parts and composition of gain of dairy goat wethers. Journal of Agricultural Science 153:1321-1340. https://doi.org/10.1017/S0021859615000532
  » https://doi.org/10.1017/S0021859615000532
• Bangar, Y. C. and Verma, M. R. 2017. Non-linear modelling to describe lactation curve in Gir crossbred cows. Journal of Animal Science and Technology 59:3. https://doi.org/10.1186/s40781-017-0128-6
  » https://doi.org/10.1186/s40781-017-0128-6
• Bauman, D. E. 2000. Regulation of nutrient partitioning during lactation: homeostasis and homeorhesis revisited. p.311-328. In: Ruminant physiology. Cronjé, P. B., ed. CABI, Wallingford.
• Bauman, D. E. and Currie, W. B. 1980. Partitioning of nutrients during pregnancy and lactation: a review of mechanisms involving homeostasis and homeorhesis. Journal of Dairy Science 63:1514-1529. https://doi.org/10.3168/jds.S0022-0302(80)83111-0
  » https://doi.org/10.3168/jds.S0022-0302(80)83111-0
• Bauman, D. E.; Everett, R. W.; Weiland, W. H. and Collier, R. J. 1999. Production responses to bovine somatotropin in northeast dairy herds. Journal of Dairy Science 82:2564-2573. https://doi.org/10.3168/jds.S0022-0302(99)75511-6
  » https://doi.org/10.3168/jds.S0022-0302(99)75511-6
• Bernabucci, U.; Ronchi, B.; Lacetera, N. and Nardone, A. 2002. Markers of oxidative status in plasma and erythrocytes of transition dairy cows during hot season. Journal of Dairy Science 85:2173-2179. https://doi.org/10.3168/jds.S0022-0302(02)74296-3
  » https://doi.org/10.3168/jds.S0022-0302(02)74296-3
• Bernabucci, U.; Ronchi, B.; Lacetera, N. and Nardone, A. 2005. Influence of body condition score on relationships between metabolic status and oxidative stress in periparturient dairy cows. Journal of Dairy Science 88:2017-2026. https://doi.org/10.3168/jds.S0022-0302(05)72878-2
  » https://doi.org/10.3168/jds.S0022-0302(05)72878-2
• Boutinaud, M.; Guinard-Flament, J. and Jammes, H. 2004. The number and activity of mammary epithelial cells, determining factors for milk production. Reproduction Nutrition Development 44:499-508. https://doi.org/10.1051/rnd:2004054
  » https://doi.org/10.1051/rnd:2004054
• Brody, S.; Turner, C. W. and Ragsdale, A. C. 1924. The relation between the initial rise and the subsequent decline of milk secretion following parturition. The Journal of General Physiology 6:541-545. https://doi.org/10.1085/jgp.6.5.541
  » https://doi.org/10.1085/jgp.6.5.541
• Buckland, S. T.; Burnham, K. P. and Augustin, N. H. 1997. Model selection: an integral part of inference. Biometrics 53:603-618. https://doi.org/10.2307/2533961
  » https://doi.org/10.2307/2533961
• Burnham, K. P. and Anderson, D. R. 2004. Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods & Research 33:261-304. https://doi.org/10.1177/0049124104268644
  » https://doi.org/10.1177/0049124104268644
• Burnham, K. P. and Anderson, D. R. 2014. P values are only an index to evidence: 20th- vs. 21st-century statistical science. Ecology 95:627-630. https://doi.org/10.1890/13-1066.1
  » https://doi.org/10.1890/13-1066.1
• Burnham, K. P.; Anderson, D. R. and Huyvaert, K. P. 2011a. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behavioral Ecology and Sociobiology 65:23-35. https://doi.org/10.1007/s00265-010-1029-6
  » https://doi.org/10.1007/s00265-010-1029-6
• Burnham, K. P.; Anderson, D. R. and Huyvaert, K. P. 2011b. Erratum to: AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behavioral Ecology and Sociobiology 65:415. https://doi.org/10.1007/s00265-010-1084-z
  » https://doi.org/10.1007/s00265-010-1084-z
• Capuco, A. V.; Akers, R. M. and Smith, J. J. 1997. Mammary growth in Holstein cows during the dry period: quantification of nucleic acids and histology. Journal Dairy of Science 80:477-487. https://doi.org/10.3168/jds.S0022-0302(97)75960-5
  » https://doi.org/10.3168/jds.S0022-0302(97)75960-5
• Capuco, A. V.; Wood, D. L.; Baldwin, R.; McLeod, K. and Paape, M. J. 2001. Mammary cell number, proliferation, and apoptosis during a bovine lactation: relation to milk production and effect of bST. Journal of Dairy Science 84:2177-2187. https://doi.org/10.3168/jds.S0022-0302(01)74664-4
  » https://doi.org/10.3168/jds.S0022-0302(01)74664-4
• Cavanaugh, J. E. 1997. Unifying the derivations for the Akaike and corrected Akaike information criteria. Statistics & Probability Letters 33:201-208. https://doi.org/10.1016/S0167-7152(96)00128-9
  » https://doi.org/10.1016/S0167-7152(96)00128-9
• Cobuci, J. A. and Costa, C. N. 2012. Persistency of lactation using random regression models and different fixed regression modeling approaches. Revista Brasileira de Zootecnia 41:1996-2004. https://doi.org/10.1590/S1516-35982012000900005
  » https://doi.org/10.1590/S1516-35982012000900005
• Congleton Jr., W. R. and Everett, R. W. 1980. Error and bias of the incomplete gamma function to describe lactation curves. Journal of Dairy Science 63:101-108. https://doi.org/10.3168/jds.S0022-0302(80)82894-3
  » https://doi.org/10.3168/jds.S0022-0302(80)82894-3
• Dijkstra, J.; Forbes, J. M. and France, J. 2005. Quantitative aspects of ruminant digestion and metabolism. 2nd ed. CABI Pub., Cambridge.
• Dijkstra, J.; France, J.; Dhanoa, M. S.; Maas, J. A.; Hanigan, M. D.; Rook, A. J. and Beever, D. E. 1997. A model to describe growth patterns of the mammary gland during pregnancy and lactation. Journal of Dairy Science 80:2340-2354. https://doi.org/10.3168/jds.S0022-0302(97)76185-X
  » https://doi.org/10.3168/jds.S0022-0302(97)76185-X
• Dijkstra, J.; Lopez, S.; Bannink, A.; Dhanoa, M. S.; Kebreab, E.; Odongo, N. E.; Nasri, M. H. F.; Behera, U. K.; Hernandez-Ferrer, D. and France, J. 2010. Evaluation of a mechanistic lactation model using cow, goat and sheep data. Journal of Agricultural Science 148:249-262. https://doi.org/10.1017/S0021859609990578
  » https://doi.org/10.1017/S0021859609990578
• Druet, T.; Jaffrézic, F.; Boichard, D. and Ducrocq, V. 2003. Modeling lactation curves and estimation of genetic parameters for first lactation test-day records of French Holstein cows. Journal of Dairy Science 86:2480-2490. https://doi.org/10.3168/jds.S0022-0302(03)73842-9
  » https://doi.org/10.3168/jds.S0022-0302(03)73842-9
• Ferreira, A. G. T.; Henrique, D. S.; Vieira, R. A. M.; Maeda, E. M. and Valotto, A. A. 2015. Fitting mathematical models to lactation curves from holstein cows in the southwestern region of the state of Parana, Brazil. Anais da Academia Brasileira de Ciências 87:503-517. https://doi.org/10.1590/0001-3765201520130514
  » https://doi.org/10.1590/0001-3765201520130514
• Fowler, P. A.; Knight, C. H.; Cameron, G. G. and Foster, M. A. 1990. In-vivo studies of mammary development in the goat using magnetic resonance imaging (MRI). Journal of Reproduction and Fertility 89:367-375. https://doi.org/10.1530/jrf.0.0890367
  » https://doi.org/10.1530/jrf.0.0890367
• Grossman, M. and Koops, W. J. 1988. Multiphasic analysis of lactation curves in dairy cattle. Journal of Dairy Science 71:1598-1608. https://doi.org/10.3168/jds.S0022-0302(88)79723-4
  » https://doi.org/10.3168/jds.S0022-0302(88)79723-4
• Hossein-Zadeh, N. G. 2014. Comparison of non-linear models to describe the lactation curves of milk yield and composition in Iranian Holsteins. Journal of Agricultural Science 152:309-324. https://doi.org/10.1017/S0021859613000415
  » https://doi.org/10.1017/S0021859613000415
• Hossein-Zadeh, N. G. 2016. Comparison of non-linear models to describe the lactation curves for milk yield and composition in buffaloes ( Bubalus bubalis ). Animal 10:248-261. https://doi.org/10.1017/S1751731115001846
  » https://doi.org/10.1017/S1751731115001846
• Hurvich, C. M. and Tsai, C.-L. 1989. Regression and time series model selection in small samples. Biometrika 76:297-307. https://doi.org/10.1093/biomet/76.2.297
  » https://doi.org/10.1093/biomet/76.2.297
• Knight, C. H. and Peaker, M. 1984. Mammary development and regression during lactation in goats in relation to milk secretion. Quarterly Journal of Experimental Physiology 69:331-338. https://doi.org/10.1113/expphysiol.1984.sp002809
  » https://doi.org/10.1113/expphysiol.1984.sp002809
• Knight, C. H. and Wilde, C. J. 1993. Mammary cell changes during pregnancy and lactation. Livestock Production Science 35:3-19. https://doi.org/10.1016/0301-6226(93)90178-K
  » https://doi.org/10.1016/0301-6226(93)90178-K
• Littell, R. C.; Milliken, G. A.; Stroup, W. W.; Wolfinger, R. D. and Schabenberger, O. 2006. SAS^®for mixed models. 2nd ed. SAS Institute Inc., Cary, USA.
• López, S.; France, J.; Odongo, N. E.; McBride, R. A.; Kebreab, E.; AlZahal, O.; McBride, B. W. and Dijkstra, J. 2015. On the analysis of Canadian Holstein dairy cow lactation curves using standard growth functions. Journal of Dairy Science 98:2701-2712. https://doi.org/10.3168/jds.2014-8132
  » https://doi.org/10.3168/jds.2014-8132
• Macciotta, N. P. P.; Dimauro, C.; Catillo, G.; Coletta, A. and Cappio-Borlino, A. 2006. Factors affecting individual lactation curve shape in i talian river buffaloes. Livestock Science 104:33-37. https://doi.org/10.1016/j.livsci.2006.03.001
  » https://doi.org/10.1016/j.livsci.2006.03.001
• Macciotta, N. P. P.; Vicario, D. and Cappio-Borlino, A. 2005. Detection of different shapes of lactation curve for milk yield in dairy cattle by empirical mathematical models. Journal of Dairy Science 88:1178-1191. https://doi.org/10.3168/jds.S0022-0302(05)72784-3
  » https://doi.org/10.3168/jds.S0022-0302(05)72784-3
• Mood, A. M.; Graybill, F. A. and Boes, D. C. 1974. Introduction to the theory of statistics. 3rd ed. McGraw-Hill Kogakusha, LTD., Tokyo.
• Morant, S. V. and Gnanasakthy, A. 1989. A new approach to the mathematical formulation of lactation curves. Animal Production 49:151-162. https://doi.org/10.1017/S000335610003227X
  » https://doi.org/10.1017/S000335610003227X
• Olori, V. E.; Brotherstone, S.; Hill, W. G. and McGuirk, B. J. 1999. Fit of standard models of the lactation curve to weekly records of milk production of cows in a single herd. Livestock Production Science 58:55-63. https://doi.org/10.1016/S0301-6226(98)00194-8
  » https://doi.org/10.1016/S0301-6226(98)00194-8
• Pinheiro, J. C. and Bates, D. M. 2000. Mixed-effects models in S and S-PLUS. Springer-Verlag Inc., New York.
• Pinheiro, J. C.; Bates, D. M.; DebRoy, S. and Sarkar, D. 2017. nlme: Linear and nonlinear mixed effects models. R package version 3.1-130. Available at: <https://CRAN.R-project.org/package=nlme>. Accessed on: Dec. 20, 2017.
  » https://CRAN.R-project.org/package=nlme
• Pollott, G. E. 2000. A biological approach to lactation curve analysis for milk yield. Journal of Dairy Science 83:2448-2458. https://doi.org/10.3168/jds.S0022-0302(00)75136-8
  » https://doi.org/10.3168/jds.S0022-0302(00)75136-8
• Quarrie, L. H.; Addey, C. V. P. and Wilde, C. J. 1994. Local regulation of mammary apoptosis in the lactating goat. Biochemical Society Transactions 22:178S. https://doi.org/10.1042/bst022178s
  » https://doi.org/10.1042/bst022178s
• Rekik, B. and Ben Gara, A. 2004. Factors affecting the occurrence of atypical lactations for Holstein−Friesian cows. Livestock Production Science 87:245-250. https://doi.org/10.1016/j.livprodsci.2003.09.023
  » https://doi.org/10.1016/j.livprodsci.2003.09.023
• Rohem Júnior, N. M.; Silva, M. C.; Abreu, M. L. C.; Oliveira, J. G.; Glória, L. S.; Tedeschi, L. O. and Vieira, R. A. M. 2020. The transit of external markers throughout the ruminant digestive tract: 1. The fitting quality of models to marker profiles in feces using an information-theoretic approach. Animal Feed Science and Technology 261:114407. https://doi.org/10.1016/j.anifeedsci.2020.114407
  » https://doi.org/10.1016/j.anifeedsci.2020.114407
• Rook, A. J.; France, J. and Dhanoa, M. S. 1993. On the mathematical description of lactation curves. The Journal of Agricultural Science 121:97-102. https://doi.org/10.1017/S002185960007684X
  » https://doi.org/10.1017/S002185960007684X
• Silva, F. F.; Zambrano, M. F. B.; Varona, L.; Glória, L. S.; Lopes, P. S.; Silva, M. V. G. B.; Arbex, W.; Lázaro, S. F.; Resende, M. D. V. and Guimarães, S. E. F. 2017. Genome association study through nonlinear mixed models revealed new candidate genes for pig growth curves. Scientia Agricola 74:1-7. https://doi.org/10.1590/1678-992x-2016-0023
  » https://doi.org/10.1590/1678-992x-2016-0023
• Soares, A. C. C.; Guimarães, S. E. F.; Kelly, M. J.; Fortes, M. R. S.; Silva, F. F.; Verardo, L. L.; Mota, R. R. and Moore, S. 2017. Multiple-trait genomewide mapping and gene network analysis for scrotal circumference growth curves in Brahman cattle. Journal of Animal Science 95:3331-3345. https://doi.org/10.2527/jas.2017.1409
  » https://doi.org/10.2527/jas.2017.1409
• Sober, E. 2002. Instrumentalism, parsimony, and the Akaike framework. Philosophy of Science 69:S112-S123. https://doi.org/10.1086/341839
  » https://doi.org/10.1086/341839
• Spiess, A.-N. and Neumeyer, N. 2010. An evaluation of R²as an inadequate measure for nonlinear models in pharmacological and biochemical research: a Monte Carlo approach. BMC Pharmacology 10:6. https://doi.org/10.1186/1471-2210-10-6
  » https://doi.org/10.1186/1471-2210-10-6
• Stefanon, B.; Colitti, M.; Gabai, G.; Knight, C. H. and Wilde, C. J. 2002. Mammary apoptosis and lactation persistency in dairy animals. Journal of Dairy Research 69:37-52. https://doi.org/10.1017/S0022029901005246
  » https://doi.org/10.1017/S0022029901005246
• Strange, R.; Li, F.; Saurer, S.; Burkhardt, A. and Friis, R. R. 1992. Apoptotic cell death and tissue remodelling during mouse mammary gland involution. Development 115:49-58. https://dev.biologists.org/content/115/1/49
  » https://dev.biologists.org/content/115/1/49
• Strathe, A. B.; Danfaer, A.; Sørensen, H. and Kebreab, E. 2010. A multilevel nonlinear mixed-effects approach to model growth in pigs. Journal of Animal Science 88:638-649. https://doi.org/10.2527/jas.2009-1822
  » https://doi.org/10.2527/jas.2009-1822
• Stroup, W. W. 2013. Generalized linear mixed models: Modern concepts, methods and applications. Taylor & Francis Group, LLC, Boca Raton.
• Sugiura, N. 1978. Further analysis of the data by Akaike's information criterion and the finite corrections. Communications in Statistics - Theory and Methods 7:13-26. https://doi.org/10.1080/03610927808827599
  » https://doi.org/10.1080/03610927808827599
• Tekerli, M.; Akinci, Z.; Dogan, I. and Akcan, A. 2000. Factors affecting the shape of lactation curves of holstein cows from the Balikesir province of Turkey. Journal of Dairy Science 83:1381-1386. https://doi.org/10.3168/jds.S0022-0302(00)75006-5
  » https://doi.org/10.3168/jds.S0022-0302(00)75006-5
• Vieira, R. A. M.; Rohem Júnior, N. M.; Abreu, M. L. C.; Silva, M. C.; Oliveira, J. G.; Tedeschi, L. O. and Glória, L. S. 2020. The transit of external markers throughout the ruminant digestive tract: 2. The estimation of fiber digestibility, ruminoreticular fill, and related biases. Animal Feed Science and Technology 261:114420. https://doi.org/10.1016/j.anifeedsci.2020.114420
  » https://doi.org/10.1016/j.anifeedsci.2020.114420
• Vieira, R. A. M.; Rohem Júnior, N. M.; Gomes, R. S.; Oliveira, T. S.; Bendia, L. C. R.; Azevedo, F. H. V.; Barbosa, D. L.; Glória, L. S. and Rodrigues, M. T. 2018. The ontogenetic allometry of body morphology and chemical composition in dairy goat wethers. Animal 12:538-553. https://doi.org/10.1017/S1751731117001884
  » https://doi.org/10.1017/S1751731117001884
• Vonesh, E. F. 2012. Generalized linear and nonlinear models for correlated data: theory and applications using SAS^® SAS Institute Inc., Cary.
• Walker, N. I.; Bennett, R. E. and Kerr, J. F. 1989. Cell death by apoptosis during involution of the lactating breast in mice and rats. The American Journal of Anatomy 185:19-32. https://doi.org/10.1002/aja.1001850104
  » https://doi.org/10.1002/aja.1001850104
• Wood, P. D. P. 1967. Algebraic model of the lactation curve in cattle. Nature 216:164-165. https://doi.org/10.1038/216164a0
  » https://doi.org/10.1038/216164a0
• Wood, P. D. P. 1968. Factors affecting persistency of lactation in cattle. Nature 218:894. https://doi.org/10.1038/218894a0
  » https://doi.org/10.1038/218894a0
• Wood, P. D. P. 1969. Factors affecting the shape of the lactation curve in cattle. Animal Production 11:307-316. https://doi.org/10.1017/S0003356100026945
  » https://doi.org/10.1017/S0003356100026945
• Wood, P. D. P. 1970. The relationship between the month of calving and milk production. Animal Production 12:253-259. https://doi.org/10.1017/S0003356100038824
  » https://doi.org/10.1017/S0003356100038824
• Wood, P. D. P. 1972. A note on seasonal fluctuations in milk production. Animal Production 15:89-92. https://doi.org/10.1017/S0003356100011260
  » https://doi.org/10.1017/S0003356100011260
• Wood, P. D. P. 1976. Algebraic models of the lactation curves for milk, fat and protein production, with estimates of seasonal variation. Animal Production 22:35-40. https://doi.org/10.1017/S000335610003539X
  » https://doi.org/10.1017/S000335610003539X
• Wood, P. D. P. 1980. Breed variations in the shape of the lactation curve of cattle and their implications for efficiency. Animal Production 31:133-141. https://doi.org/10.1017/S0003356100024363
  » https://doi.org/10.1017/S0003356100024363
• Zanton, G. I. and Heinrichs, A. J. 2009. Evaluation of modeling procedure for fitting in situ feed degradation profiles. Journal of Animal Science 87:2080-2088. https://doi.org/10.2527/jas.2008-1180
  » https://doi.org/10.2527/jas.2008-1180

Publication Dates

• Publication in this collection
  16 Oct 2020
• Date of issue
  2020

History

• Received
  03 Feb 2020
• Accepted
  15 June 2020