Investigating the effect of temporal, geographic, and management factors on US Holstein lactation curve parameters

We fit the Wood’s lactation model to an extensive database of test-day milk production records of US Holstein cows to obtain lactation-specific parameter estimates and investigated the effects of temporal, spatial, and management factors on lactation curve parameters and 305-d milk yield. Our approach included 2 steps as follows: (1) individual animal-parity parameter estimation with nonlinear least-squares optimization of the Wood’s lactation curve parameters, and (2) mixed-effects model analysis of 8,595,413 sets of parameter estimates from individual lactation curves. Further, we conducted an analysis that included all parities and a separate analysis for first lactation heifers. Results showed that parity had the most significant effect on the scale (parameter a ), the rate of decay (parameter c ), and the 305-d milk yield. The month of calving had the largest effect on the rate of increase (parameter b ) for models fit with data from all lactations. The calving month had the most significant effect on all lactation curve parameters for first lactation models. However, age at first calving, year, and milking frequency accounted for a higher proportion of the variance than month for first lactation 305-d milk yield. All parameter estimates and 305-d milk yield increased as parity increased; parameter a and 305-d milk yield rose, and parameters b and c decreased as year and milking frequency increased. Calving month estimates parameters a , b , c , and 305-d milk yield were the lowest values for September, May, June, and July, respectively. The results also indicated the random effects of herd and cow improved model fit. Lactation curve parameter estimates from the mixed-model analysis of individual lactation curve fits describe well US Holstein lactation curves according to temporal, spatial, and management factors.


INTRODUCTION
The lactation curve represents milk production as a reproducible pattern that can be expressed mathematically as functions of time.The functions allow for milk yield predictions at the farm or individual level.The use of lactation curve models to predict target milk yield is also a practical method to describe cow and herd milk production without previous milk production data (Rotz et al., 2013;Calsamiglia et al., 2018) and with only a small number of parameters (Macciotta et al., 2011).
The shape of the lactation curve, which can be interpreted through parameters that define the representative equation, provides insight for implementing effective production strategies on farms.In general, lactation curves have an ascending phase leading up to the peak and a descending phase following the peak.The physiology behind the lactation curve shape is related to cell division and cellular differentiation (Knight, 2000;Capuco et al., 2001) and gene expression (Atashi et al., 2020) throughout a lactation.Parametric models, such as the Wood's model (Wood, 1967), have been widely applied to characterize the behavior of milk production (Bouallegue and M'Hamdi, 2020).
Numerous factors may affect the values of the parameters that define a lactation curve and ultimately have a cumulative effect on milk productivity throughout the lactation.For example, the geographic region where the cow is located (Cole et al., 2011), the season of the year when the lactation starts (Torshizi, 2016), milking frequency (Stelwagen, 2001), and age at first calving (Atashi et al., 2020) have all been reported to alter the shape and scale of lactation curves.Thus, updating and obtaining lactation curve parameters that reflect factors known to influence milk production patterns remains a critical task for practical application of these models.

MATERIALS AND METHODS
No animals were used in this study, and ethical approval for the use of animals was thus deemed unnecessary.

Data
The data set used in this study was provided by the US Council on Dairy Cattle Breeding (https: / / www .uscdcb.com/ ) and contained individual test-day records for 10.15 million lactations from US Holstein cows located in 45 states from January 2006 to December 2016.Each lactation was composed of at least 10 test-day records (average lactation length of 342 d) and the calving date associated with the beginning of the lactation.

Lactation-Specific Curve Fitting
Lactation curves are mathematical equations that describe milk yield as a function of DIM (i.e., days postpartum).Wood's model (Wood, 1967), the most commonly used parametric model for describing lactation curves (Dijkstra et al., 2010;Bouallegue and M'hamdi, 2020), is a gamma function with 3 parameters (a, b, and c).It describes milk production as y = at b e −ct + ε, where y is the observed daily milk yield (kg), t is the number of days postpartum (or DIM), parameter a is a scale factor for the production level, b indicates the growth rate in milk yield until peak production, c describes the rate of decline after the peak (i.e., a greater parameter c is associated with lower persistence), and ε is the random residual (Wood, 1967(Wood, , 1972)).We used a nonlinear least-square method to fit each set of individual cow lactation test-day records to Wood's function using the nlme v3.1-152 package (Pinheiro et al., 2021) in R v4.0.3 (https: / / www .r-project .org/ ) and RStudio v1.3.1093 (www .rstudio .com).We collected the lactation-specific estimates for a, b, c, and the residual sum of squares (a measure of goodness of fit) from the nonlinear fitting procedure.Due to data errors and outliers, some of the resulting parameter estimates did not follow the typical shape of lactation curves.Thus, the following criteria were used to classify estimated lactation curve fits as atypical: (1) b or c parameters were negative (10.80%), which would result in an inverted lactation curve (Wood, 1976;Macciotta et al., 2011); (2) the peak of production occurred after 300 DIM (1.01%); (3) the residual sum of squares was greater than 1,000, which was found to produce unreliable fitting results (1.02%); (4) estimates of either a, b, or c parameters were more than 3 standard deviations away from their respective means (2.48%); and (5) estimated 305-d milk yield was greater than 5 standard deviations from its respective mean (0.01%).If one parameter set met any of the above criteria, we excluded the atypical lactation curve and all parameters associated with it, resulting in 8,599,526 lactations and their corresponding parameters retained for further analysis.The cumulative 305-d milk yield, commonly used to compare lactation performance, was calculated Figure 1 shows the distributions of lactation-specific model parameters and 305-d milk yield values for the subset of lactations retained for further analyses.

Linear Mixed-Effects Models of Lactation Curve Parameters and 305-d Milk Yield
Parameter estimates obtained from fitted individual lactation curves were used to investigate the effects of spatial (farm region), temporal (calving year and calving month), and management (milking frequency, age at calving for first lactation, and parity) factors on the lactation shape, as well as on the 305-d milk yield.We categorized the 45 states available in our data set into 10 regions in accordance with the USDA convention (USDA-ERS, 2022) to explore the effect of location.The regions included West Coast, Mountain, Northern Plains, Southern Plains, Lake States, Corn Belt, Delta States, Appalachian, Southeast, and Northeast.We then separated records from Wisconsin, Pennsylvania, and New York (the 3 states with the most lactation records in our data set) from the Lake States and the Northeast (Figure 2, Table 1).Some states on the West Coast and Mountain regions that are large dairy production states (e.g., California and Idaho; USDA-ERS, 2022) were not separated from their respective regions because they were not as well represented in our data set.
We accounted for correlated data associated with multiple lactations by the same animal as well as lactations within the same herd through the use of linear mixed-effects models, with animals and herds included as random effects.To eliminate herds that were not well represented in the data set, we further excluded 4,426 lactations from those with fewer than 5 lactations available, leaving 8,595,100 lactations (from a total of 21,056 herds) for further analyses.
To assess the effects of spatial, temporal, and management factors on the lactation curve parameters (a, b, c) and 305-d milk yield, we used the following mixed effects model: where y ijklmnpq was the response variable (lactation curve parameter a, b, c, or 305-d total milk yield) related to the ith lactation from the jth parity (3 groups, first, second, and third and later lactations), in the kth calving year (11 years, from 2006-2016), in the lth month of calving (12 mo, from January-December), in the mth region (13 levels, as described above), and at the nth milking frequency (2 levels, 2×/d or 3×/d) of the pth animal in the qth herd.The factors parity j , year k , month l , region m , and milkingfreq n referred to the fixed effects of parity, calving year, calving month, US region where the herd was located, and daily milking frequency, respectively.In addition, herd q (n = 21,056) and animal p (n = 5,375,200) were the random effects of herd and animal, respectively.Last, ε ijklmnpq represented the residual term, with an expectation of 0 and variance equal to σ 2 .Table 1 reports the number of lactations for each factor level.Eq. 1 was fit independently to each of the 3 parameter estimates from the selected subset of individual lactation curves (8,599,526 lactations).
In addition, we fit the following model separately for first lactations (3,324,803 lactations): The parameters in the model described in Eq. 2 were the same as those in the model for all lactations except that the effect of parity was replaced by the effect of age at calving.The corresponding variable, age at calving j , was a continuous variable with a median of 746 d, a mean value at 764 d, and a standard deviation of 84 d (Figure 3).This model included only herd as a random effect (n = 20,509).
for all-lactations models and for first-lactation models.
We conducted the linear mixed-effects model analyses using the lme4 (Bates et al., 2015) package in R statistical software.We investigated the effect of each factor on each response variable by examining the average fitted parameter estimate for each factor level and the ANOVA tables for the linear mixed-effects models.Additionally, we used likelihood ratio tests to test the significance of including random effects.Furthermore, we reported the coefficients of mixed-effects models using parity, year, month, region, and milking frequency specific lactation curve parameters.

RESULTS
Visual inspection of residual plots did not reveal any important deviations from homoscedasticity or normality in the mixed-effects model fits.The fixed effects of parity (or age at calving), year and month of calving, region, and milking frequency were significant (P < 0.05) in both sets of models (all lactations and first lactation only) for all variables (i.e., parameters a, b, and c, and 305-d milk yield), except for the effect of milking frequency on the first-lactation curve model of parameter b (Tables 2 and 3).Likelihood ratio tests demonstrated significance (P < 0.05) of the random effects of herd and animal for the all-lactations models and herd for first-lactation models.
The proportions of the total sum of squares for the scale parameter a were ranked as follows: lactation group (90.3%), month (7.7%), year (1.2%), milking frequency (0.8%), and region (0.1%).For the rate of change parameters b and c, the milking frequency accounted for the lowest proportion of the sum of squares, and the month of calving factor comprised the largest proportion of the sum of squares for parameter b, accounting for 59.3% of the total fixed-effect variation.Ranking of the proportion of the fixed-effect sum of squares for 305-d milk yield was similar to that for parameter a except that the year was ranked higher than the month (Table 2).The first-lactation mixedeffects models included the age at first calving, which explained the largest proportion (35.2%) of the variability in 305-d milk yield (Table 3).
Tables 4 and 5 summarize the estimates and standard errors for each coefficient in the mixed models.To reconstruct estimated lactation curves for specific categories, the mean of each parameter can be adjusted by adding the corresponding coefficients.In general, for all parameters, an increasing trend is observed as lactation number increases 2.47-kg increase of 305-d milk yield, and a decrease of 0.000032 of parameter b.Demonstratively, to obtain the parameter a for a second lactation that began in August 2010 in a PA herd milking 2×/d, one could extract the relevant coefficients from Table 4 and add them as follows: a = 19.9+ 2.16 + (−0.11) + 1.15 + (−0.74) = 21.4.
Figure 4 shows the lactation curves of mean parameter values from each parity and milking frequency.Because 2×/d milking had a lower peak level and a faster rise and decline, the time of peak production did not differ between milking frequencies.Moreover, the lactation shapes of first lactation had lower rate of increase and decrease than the others and intersected with the second-lactation curves.Figure 5 presents the upper and lower boundaries (±3 SD from the mean) of estimates for parameters a, b, and c across levels for the year, month, and region factors.The ranges of year and month parameters were modest; however, the regional ranges were substantial, particularly in regions with fewer data.New York and Pennsylvania exhibited similar parameter ranges and were closer to Wisconsin than the rest of the Northeast; the Northeast had smaller estimates for parameter b.Wisconsin, on the other hand, had a higher estimate for parameter a and a lower parameter b than the other Lake states.Delta states and the Southeast, Mountain, and West Coast had greater b and c parameters, with the widest range in parameter estimates indicating large regional variance.The West Coast had the highest parameter a and a large parameter b, which resulted in a greater lactation peak production among those 4 regions (Delta, Southeast, Mountain, and West Coast).Parameters from regions with a low standard deviation are more dependable in characterizing the region's general lactation curve.
Tables 6 and 7 contain estimates of variance components.For all 3 parameters, the majority (>70%) of variance remained unexplained in the residual variance.The proportion of total variance explained by the herd was similar to the proportion explained by the animal for the parameters a (8.0% and 7.1%) and b (9.2% and 6.9%); however, the proportion explained by the animal for parameter c was greater (18.0%) than the proportion explained by herd (6.7%).For 305-d milk yield, the herd effect (44.8%) explained a greater proportion of the variance than animal (25.8%) in the all-lactations model (Table 6).Herd explained 8.5%, 8.0%, 7.5%, and 46.5% of total variance in the first lactation models for parameters a, b, c, and 305-d milk yield, respectively (Table 7).

DISCUSSION
Our objective was to study lactation curves that follow typical and biologically reasonable lactation patterns, trying to avoid as much as possible problems related to data recording (e.g., discrepant observations, missing data, and minimum number of test-day observations per lactation) or biological issues such as disease.To accomplish this goal, we deliberately censored data to ensure all resulting lactation curves conformed with the expected shape of a lactation curve.Most lactations were censored due to missing data, recording errors, or highly abnormal lactations (e.g., heavily affected performance due to external factors such as disease).Such abnormal lactations were detected, for example, based on negative estimates for either the b or c parameters, which would result in inverted lactation curves that do not conform with expected biological lactation curves.We chose the Wood's model for this study because of its simplicity and flexibility, for which it has been selected by a variety of studies on milking performance (Silvestre et al., 2009;Cole et al., 2011) and many simulation analyses (Rotz et al., 1999;Calsamiglia et al., 2018).We also plan to use these results for simulation with the Ruminant Farm System Model (RuFaS; Hansen et al., 2021).

Fixed-Effects Factors
Lactation group parity accounted for 90%, 34%, 92%, and 93% of the variation in parameters a, b, c, and the 305-d milk yield, respectively.As a result of the considerable effect of lactation group on parameters a and c, even lactation groups with comparable 305-d production levels exhibited substantial differences in the shape of the estimated lactation curves, particularly during the decay period of lactation.According to our findings and those of Ehrlich (2011), it is critical to distinguish the lactation curve parameters for the second lactation from those of later lactations.As illustrated in Figure 3, the milk yield of first-lactation cows after 200 DIM for 3×/d (and 222 DIM for 2×/d) would be greater than that of the second-lactation cows.This information could be used, for example, in reproductive management decisions concerning the length of the open period for second-lactation cows or dry-off timing for first-lactation cows to benefit total farm productivity.
The fixed effect of calving month accounted for 8%, 59%, 7%, and 2% of the variation on parameters a, b, c, and the 305-d milk yield, respectively, in all-lactations models, and 67%, 86%, 87%, 7.4%, respectively, in first-lactation models.The calving month's considerable effect on parameter b indicated that the calving season influenced the rate of milk yield increase in the early days across all lactation groups.Also, the large effect of calving month on parameters in comparison to the 305-d milk yield suggested that even when the total milk yields were similar across calving months, the shape of the lactations, namely the rates of rise and decay, the initial and peak yield, or the peak time, were different.Thus, knowing the lactation curve parameters for a particular calving month would produce more precise milk yield projections to inform nutrition, reproduction, and facility management decisions.Specifically, parameter a was greater than its mean for all months from February to June; parameter b was greater than its mean value from August to February; and parameter c was greater than its mean from October to March.A greater parameter b in August to February denotes a faster post-parturient rise in production, which might suggest fewer transition management issues during these months' parturitions.The quotient of parameters b and c (b/c) determines the peak time (d) of lactation curves, and different estimates of these 2 parameters for different calving months result in changes in the expected peak time depending on the month of calving.For an August calving, the expected peak time is 78 DIM, which was the calving month with the longest expected time to peak production in our analysis.Conversely, calving in April resulted in the shortest expected time to peak production (68 DIM).Apart from the parameters, Torshizi (2016) and Dedkova and Nemcova (2003) found that winter calving led to a higher total milk production, which was also corroborated by our findings.Lactations starting during the hot season, namely from May to September, Parameters were estimated for Wood's lactation curve model as follows: y = at b e −ct + ε, where y is the observed daily milk yield (kg), t is the number of days postpartum (or DIM), parameter a is a scale factor for the production level, b indicates the growth rate in milk yield until peak production, c describes the rate of decline after the peak, and ε is the random residual (Wood, 1967).had lower-than-average estimates for the 305-d milk yield (e.g., July calving was associated with 3% milk yield decrease from the average), which might be due to negative effects from heat stress.Cows milked 3×/d had a higher milk production than 2×/d: an average increase of 2.2 kg/d and 2.5 kg/d for the first lactation and all-lactations 305-d milk yield model, respectively.This is consistent with previous studies; for example, Hart et al. (2013) reported an increase of 2.5 kg/d and 3.2 kg/d for first and second and later lactations when moving from milking 2×/d to 3×/d, respectively, and Erdman and Varner (1995) summarized the effect of milking 3×/d from the literature as 3.5 kg/d for all lactations.The magnitude of the change in 305-d yield due to milking frequency from our models (7.7% for all-lactations and 11.8% for first lactation) was less than the 15% reported by Smith et al. (2002) or the 12% for total and 14% for first lactation reported by Gisi et al. (1986).It is possible that the increased production over the last decade affected the base production level; in other words, the effect of increasing milking frequency may not be proportionate to the production as suggested by Erdman and Varner (1995).Also, Smith et al. (2002) reported results from an analysis of a DHI data set in which only 7.1% of the herds milked 3×/d in 2000, compared with our data set in which 44.5% of the herds milked 3×/d from 2006 to 2016.Thus, our study included a more balanced num-  Parameters were estimated for Wood's lactation curve model as follows: y = at b e −ct + ε, where y is the observed daily milk yield (kg), t is the number of days postpartum (or DIM), parameter a is a scale factor for the production level, b indicates the growth rate in milk yield until peak production, c describes the rate of decline after the peak, and ε is the random residual (Wood, 1967).ber of records from both milking frequencies, reducing the effect of other confounding factors associated with production such as geographic region and management.
The peak time and the overall shape of the lactation curves were not much different between 3×/d and 2×/d milkings (Figure 4).Capper and Cady (2020) reported that the national mean annual milk yield per cow increased from 9,164 kg/yr to 10,406 kg/yr (+13.6%)from 2007 to 2017.However, our data, before the mixed models (averaged individual nonlinear fitting results), showed an increase from 10,747 kg/yr to 11,812 kg/yr (+9.9%) over the 11yr period, suggesting a higher baseline production level but a smaller increase over a similar time period to that reported in Capper and Cady (2020).Between 2006 and 2016, the estimated overall mean 305-d milk yield in the mixed model outcomes increased from 10,108 kg/ yr to 10,577 kg/yr (+4.6%), reflecting a trend in milk production increase according to calving year.Over the years, the increasing milk yield trend (increasing parameter a and 305-d milk yield, decreasing parameter b and c, and later peak time) has been interpreted as a Parameters were estimated for Wood's lactation curve model as follows: y = at b e −ct + ε, where y is the observed daily milk yield (kg), t is the number of days postpartum (or DIM), parameter a is a scale factor for the production level, b indicates the growth rate in milk yield until peak production, c describes the rate of decline after the peak, and ε is the random residual (Wood, 1967).response to genetic selection (García-Ruiz et al., 2016) and management practice advancements (Brito et al., 2021), such as increasing milking frequency, as previously noted.
The geographic region affected milk production and lactation curve parameters (Carabaño et al., 1990;Hare et al., 2004;Cole et al., 2011).Cole et al. (2011) reported the lowest milk yields in the Southeast region, which is consistent with our findings, in which the Delta, Southern Plains, and Southeast regions (Southeast region in Cole et al., 2011) had the lowest 305-d milk yield.Although Cole et al. (2011) reported the highest milk yield in the Northeast, the West Coast was the region in our study with the highest 305-d milk yield.Furthermore, although all variables included in the mixed model were significantly affected, the proportion of variation explained by the region factor was small (<1.5%) for all models.This could be partially related to the substantial variation in total milk production within each region, which limited the explanatory power of region factor.Cole et al. (2011) also noted that it might be beneficial to define a region in a manner that captures the effect of environment on lactation curve parameters.Thus, future work could use Parameters were estimated for Wood's lactation curve model as follows: y = at b e −ct + ε, where y is the observed daily milk yield (kg), t is the number of days postpartum (or DIM), parameter a is a scale factor for the production level, b indicates the growth rate in milk yield until peak production, c describes the rate of decline after the peak, and ε is the random residual (Wood, 1967).alternative methods to define geographic factors based, for example, on average temperature and humidity or by longitude and latitude rather than state boundaries that do not necessarily correlate to environmental conditions experienced by the cow.
In the first lactation mixed model, age at first calving explained 16.3%, 1.4%, 1.4%, and 22.9% of the variation on the parameters a, b, c, and 305-d milk yield, respectively.The associations between age at first calving and lactation parameters were significant (P < 0.005), which is consistent with previous findings from Torshizi (2016).Increased age at first calving led to increased 305-d milk yield, increased values of parameters a and c, and decreased value of parameter b, similar to what Atashi et al. (2020) reported (with a few exceptions for parameters a and b).A positive relationship between 305-d milk yield and age at first calving was also consistent with previous studies (Ettema and Santos, 2004;Albarrán-Portillo and Pollott, 2011).Our model estimated that a heifer who is 1 d older at first calving resulted in a 2.47 kg of milk/lactation increase, which is comparable but higher than what was reported by Berry and Cromie (2009), who stated that a 1-mo increase in age at first calving resulted in a 55.5 kg of milk/lactation increase.Muir et al. (2004) found that younger heifers have increased lactation persistency, which is in alignment with our findings that parameter c increased with increasing age at first calving.According to Teke and Murat (2013), the increased 305-d milk yield related to the increased age at first calving can be explained by the additional time for mammary gland development.However, this relationship between age at first calving and 305-d milk yield had been reported in previous studies to be positive only when age at first calving happened between mo 18 and 26 (and negative thereafter; Nilforooshan and Edriss, 2004;Torshizi, 2016) or simply negative (Bewley et al., 2001).The median age at first calving in our data set (~24.5 mo)  is at the upper range of that reported by Nilforooshan and Edriss (2004), and 2 standard deviations above the mean (~30.5 mo) is well beyond that range.Thus, our results suggest that the increase in first-lactation milk yield due to delayed first calving may apply to older heifers than previously thought.

Random Effects Factors
In our study, including animal as a random variable improved model fit as measured by the Akaike information criterion and P-value of likelihood ratio tests.These results are consistent with the work of Piccardi et al. (2017), who included the random effect of ani-mal when they fit nonlinear mixed models of lactation curves directly and concluded that incorporating random deviations enables the model to estimate the expected production of an average cow more accurately.In alignment with the considerable difference in parameters between herds reported by Ehrlich (2013), and large diversity in parameter estimates from individual lactation curves reported by Macciotta et al. (2005), our results indicated including random effects of herd and animal helped explain random variation and enabled more precise estimation of fixed effects.In the all-lactation models, the animal accounted for a greater proportion of the variance than herd, whereas the herd was responsible for more of the variation in the 305-d milk yield than animal.The range of the proportion of total random variance explained by the herd effect was similar in the all-lactations and first-lactation models, indicating that herds affect animals from different parities in a similar manner.Further, the proportions of total random variance explained by the herd effect were large, which demonstrates the considerable variation between herds in our data set.

Applications
The lactation curve parameter estimates from our mixed models have a variety of potential applications; for example, they can serve as the expected lactation performance against which observed production can be compared for early detection on disease (Lucey et al., 1986) and identification of potential nutrition problems (Amaral-Phillips, 2015).Another direct application is through use of the parameters in dairy farm simulation models to estimate daily milk production, which informs estimates of BW change, and nutrient requirement calculations, diet formulation, and income prediction.Although milk production from dairy cattle increases and the shape of the lactation curve varies in response to many factors, most dairy farm simulation models available today use outdated lactation curve estimates that do not vary with geographic or management conditions.For example, the Integrated Farm System Model (Rotz et al., 2013) used lactation curve parameters from 1999 (Rotz et al., 1999), Calsamiglia et al. (2018) implemented lactation curve parameters from 1967 (Wood, 1967), andSilva-Villacorta et al. (2016) incorporated lactation curve parameters from 2006.Other lactation curve fitting analyses have been developed for specific farms or herds using more recent data (Murphy et al., 2014;Masía et al., 2020); however, due to the small number of animals included in these parameterizations, inferences from those works are limited to the conditions of the data used for model fitting.Furthermore, when characterizing lactation curves, including additional factors that can modify the lactation curve parameters enables model users to tailor the application to accurately represent their conditions of interest.
The findings of this study may contribute to a more accurate representation of current animal performance in a dairy farm system model, such as RuFaS (Li et al., 2018;Kebreab et al., 2019;Hansen et al., 2021).The RuFaS model is a daily time-step, process-based, whole-farm model that simulates individual animal production, growth, and life events daily.Accurate modeling of each animal's targeted daily milk yield is crucial for the model to reflect nutrient flows and overall farm performance.The results of the mixed models will be integrated into the RuFaS model to more precisely model lactation performance similar to the use of lactation curve parameter baseline estimates in other simulation models (Ben Abdelkrim et al., 2021).

CONCLUSIONS
We found that all temporal (calving year and month), spatial (location of the farm), and management (milking frequency, age at first calving, and parity) factors had a significant effect on 305-d milk yield and on the lactation curve shape, as defined by the 3 parameters from the Wood's Model.By including herd and animal as random effects in a mixed-model analysis of individual lactation parameter estimates, we were able to account for between-herd and animal variability to provide better predictions than reduced, fixed-effect models.Additionally, the estimates provided for each factor level of the fixed effects for age-at-first-calving, parity, year, month, region, and milking frequency support estimation of specific lactation curve parameters for future applications.The mixed models' estimated lactation curve parameters characterized well lactation curves in terms of temporal, spatial, and management factors.

ACKNOWLEDGMENTS
This study was supported by the Food and Agriculture Cyberinformatics and Tools grant no.2019-68017-

,
which represents the sum of estimated daily production during the first 305 d of lactation.

Figure 1 .
Figure 1.Histograms of fitted individual lactation curve (y = at b e −ct + ε, where y is the observed daily milk yield (kg), t is the number of days postpartum (or DIM), parameter a is a scale factor for the production level, b indicates the growth rate in milk yield until peak production, c describes the rate of decline after the peak, and ε is the random residual; Wood, 1967) parameters a, b, c, and 305-d milk yield a te dt b ct

Figure 2 .
Figure 2. Map showing different US regions in the study.
Figure 3. Histogram of number of lactations according to age at first calving.
3 305-d milk yield was calculated from fitted Wood's model parameters as follows: rows are expressed as a difference from the mean value.
Figure 4. Lactation curves plotted according to the estimated mean of the lactation curve parameters for each lactation group and milking frequency.Triangles indicate the peak.

Figure 5 .
Figure 5. Upper and lower boundaries (±3 SD from the mean) of lactation curve (y = at b e −ct + ε, where y is the observed daily milk yield (kg), t is the number of days postpartum (or DIM), parameter a is a scale factor for the production level, b indicates the growth rate in milk yield until peak production, c describes the rate of decline after the peak, and ε is the random residual; Wood, 1967) parameters [estimated from the mixed model for each year (top), month (middle), and region (bottom)].
Li et al.: TEMPORAL, GEOGRAPHIC, AND MANAGEMENT FACTORS ON LACTATION

Table 1 .
Li et al.: TEMPORAL, GEOGRAPHIC, AND MANAGEMENT FACTORS ON LACTATION Number of lactations in each level of factors in the linear mixed model for all-lactations models 1 Region: following the USDA convention according to https: / / www .ers.usda.gov/data -products/ dairy -data/ .

Table 2 .
Li et al.: TEMPORAL, GEOGRAPHIC, AND MANAGEMENT FACTORS ON LACTATION The residual sum of squares (Sum sq) and proportion of variance of the fixed effects in mixed models of the lactation curve parameters and 305-d milk yield for all-lactations models 2

Table 3 .
The residual sum of squares (Sum sq) and proportion of variance of the fixed effects in mixed models on the lactation curve parameters and 305-d milk yield for first-lactation models 2

Table 4 .
Li et al.: TEMPORAL, GEOGRAPHIC, AND MANAGEMENT FACTORS ON LACTATION Estimate (SE) of the coefficients of variables in the mixed models of lactation curve parameters and 305-d milk yield for all-lactations models

Table 5 .
Li et al.: TEMPORAL, GEOGRAPHIC, AND MANAGEMENT FACTORS ON LACTATION Estimate (SE) of the coefficients of variables in the mixed models of lactation curve parameters and 305-d milk yield for the firstlactation models

Table 6 .
Variance estimate and proportion of the total variance of random effects in mixed models on lactation curve parameters and 305-d milk yield for all-lactations models 1 P < 0.0001 for all likelihood ratio test of random effects.

Table 7 .
Li et al.: TEMPORAL, GEOGRAPHIC, AND MANAGEMENT FACTORS ON LACTATION Variance estimate and proportion of the total variance of random effects in mixed models on lactation curve parameters and 305-d milk yield for first-lactation models P < 0.0001 for all likelihood ratio test of random effects.29935/projectaccession no.1019780 from the USDA National Institute of Food and Agriculture (Washington, DC).The authors have not stated any conflicts of interest.