Genetic improvement aimed at increasing milk production in warm tropical regions is based, among other strategies, on crosses between adapted breeds (BI; Bos indicus) and specialized breeds from temperate and cold climes (B7; Bos taurus). Different degrees of crossbreeding provide different advantages¹. Beginning in the 1960s, crossbreeding began between the Holstein and Zebu breeds in Cuba with the aim of developing genotypes apt for milk production and adapted to adverse tropical conditions. Two synthetic breeds resulted^2,3: Siboney (5/8 Holstein + 3/8 Zebu) and Mambí (3/4 Holstein + 1/4 Zebu). The Siboney breed was created as part of the National Genetic Improvement Program to produce genotypes less demanding than Holstein that could attain adequate production potential under tropical conditions with feeding based mainly on grazing^4,5.

The lactation curve (LC) is used to analyze dairy cow milk production in three periods with four main components: 1) initial production; 2) increased production or ascending phase; 3) maximum production or peak phase; and 4) descending rate or reduced production, called persistence^6,7,8. Nonlinear models (NLM) have been used to characterize and analyze LCs. In this approach, regression coefficients are interpreted in association with LC stages, or help to derive other indicators that describe the LC^9,10. Adjusting NLM has been done using nonlinear regression analysis involving iterative methods such as those of Gauss, Marquardt and Newton, among others¹¹. Statistically analyzing milk production through lactation measurements involves repeated measurements over time. This generates a variance/covariance structure that can be analyzed using mixed models^12,13, with positive effects in NLM fit analyses. Nonlinear regression analyses, including random effects related to the regression coefficients, are run using a different iterative method (e.g. Firo, Gauss, Hardy and Isamp)^14,15 than used in NLMs with only fixed effects.

Understanding LC is important for dairy cattle management, feed and genetic improvement programs¹⁶. Analysis of LC is done at 305 d lactation in intensive or corral production systems with specialized breeds such as Holstein. However, the environmental effects and handling conditions in tropical production systems shorten cow productive life, and as a consequence the definition, analysis and characterization of lactation is done over shorter intervals, from 210 to 240 d^17,18,19. In an effort to characterize and analyze the Siboney cattle LC, the present study objectives were to a) select a NLM that fits the Siboney LC; b) compare nonlinear regression procedures for the only fixed effects NLM versus the mixed NLM; and c) use the NLM with the best fit to estimate LC components or indicators.

The data analyzed in this study were from the data base of the Research Center for Genetic Improvement of Tropical Livestock (Centro de Investigaciones para el Mejoramiento Animal de la Ganadería Tropical), a division of the Ministry of Agriculture of the Republic of Cuba. A total of 15,324 daily milk production (kg) records were analyzed, in an interval of 1 to 240 d, for first lactations in 2,809 Siboney cows that had given birth between January 2000 and December 2012. These data were collected from 28 herds located on the island of Cuba (20 and 23° N, 74 and 85° W) and managed using a grazing system. The island experiences two clear seasons, a rainy season from May to October during which 70 to 80 % of rainfall occurs, and a dry season from November to April. Mean annual temperature is 23.1 °C and average relative humidity ranges from 60 to 70 % in the day to 89 to 90 % at night²⁰.

Five NLMs were evaluated for their effectiveness in analyzing and characterizing LC^7,9,21:

Where: y_t is milk production (kg) on day t (1 to 240); β₁, β₂ and β₃ are each model's regression coefficients; and ε_ij are random effects, associated with errors or residuals, of j-th production on the i-th day of lactation. Four LC indicators were estimated using each model's β values: lactation initial production (IP); days to maximum production (DPMAX); peak or maximum production (PMAX); and total or accumulated production (PTOTAL) at 240 d lactation ^9,10,22.

Two analyses were run with the NLMIXED procedure in the SAS program, using the Gauss iterative method^13,23: 1) a NLM without random effects related to the regression coefficients and only the random effect (e) of the residuals (ANA1); and 2) a mixed NLM in which the random effect of each model's β₁ was added to the residuals (ANA2). The non-correlated random effects of ε_ij and β₁ fit the assumptions of a normal distribution, with a mean equal to zero and variances of σ²_e and σ²_β1, respectively. The random effects of the regression coefficients β₂ and β₃ were not included because their variance components were very low magnitude and tended towards zero. Selection of the best fit model was done based on six criteria^24,25,26:

Where pmi is the production measured on the i-th day of lactation; epi is the estimated production for the i-th day of lactation; n is the total number of observations; ssr is the sum of squares of residuals; tss is the total sum of squares; Log lik is the likelihood function logarithm; and k is the number of model parameters. The APE analyzes the relationship between measured and estimated production for the i-th day of lactation, and, based on the combination of "+" and "-" signs, the NLM is shown to either over- (+) or underestimate (-) the predictions. For APE, PEV, AIC and BIC, the model with the lowest value is considered that with the best fit. In contrast, for R² the model with the highest value is considered to have the best fit. The DW statistic analyzes autocorrelations in the errors, using three categories: a) 2<DW<4 indicates a negative autocorrelation; b) 0<DW<2 indicates no autocorrelation; and c) DW≤0 indicates a positive autocorrelation.

No differences in APE, PEV and DW were present between the two analyses run to select the model with the best fit (Table 1). Based on APE and the DW statistic, all the NLM predictions tended to underestimate the analyzed information, with no autocorrelation between residuals. The mixed NLM exhibited the best fit based on the three results in Table 1: a) the residual variances of ANA2 were lower (up to 82 % in the SIK and WOD models) than those estimated in ANA1; b) consequently, the ANA2 R² values (97 % on average) were higher than those produced by ANA1 (87 % in all), in an interval of 9 to 12 %; and c) the AIC and BIC criteria congointly exhibited the best fit in ANA2. As a function of the AIC, BIC and R² criteria, the best fit within ANA2 was with the WOD model, followed by the SIK and BRO models. With the ANA2 results, estimations were made for IP, DPMAX, PMAX and PTOTAL (Table 1).

The WOD NLM is widely seen to have the best statistical fit when using different selection and hierarchization criteria; it is the most popular model to describe dairy cow LC. Its regression coefficients are associated with different LC components and allow derivation of other partial or total production indicators over time^7,8,27. In one study, the WOD NLM was found to be the most apt of sixteen different NLM in LC for six cross genotypes using different BT:BI proportions²¹. This model was also found to be apt, with an ideal fit, for modeling milk production in beef cattle in different sampling schemes^28,29.

In the NLM, β₁ is associated with lactation initial production, β₂ with the ascending phase or peak production, and β₃ with descending rate or persistence. The combination of signs in β₂ and β₃ showed the WOD and WIL models to fit four unimodal curve types³⁰: 1) standard or typical curve, WOD + β₂,- β₃ and WIL - β₂,- β₃; 2) continuously descending curve (atypical), WOD -β₂,-β₃ and WIL +β₂,-β₃; 3) inverted standard curve (U), WOD - β₂,+ β₃ and WIL + β₂,+ β₃; and 4) continuously increasing curve: WOD + β₂,+ β₃ and WIL - β₂,+ β₃. The standard LC based on the WIL model and the continuously increasing LC based on the WOD model exhibited similar estimations for PMAX and PTOTAL (Figure 1), but differed notably in their estimations of IP and DPMAX (Table 1). In contrast to the high β₁ estimates they produced, both the COB and BRO models generated underestimates for initial lactation points such as IP (1.8 kg) and DPMAX (20 d). This may indicate that these models have certain problems estimating the early portions of the LC.

Figure 1 Siboney cattle lactation curves adjusted to 240 days based on the Wiltmink (a; WIL) and Wood (b; WOD) models 

Average PMAX was 7.8 kg and average PTOTAL was 1,660 kg, levels higher than previously reported for Siboney cattle^5,31,32, and other BI x BT genotypes¹⁹. However, they are lower than reported in Mambí cattle^33,34.

Dairy cattle production systems under warm tropical conditions depend on factors linked to the animal, environment, feeding and production technology. Variability in production at the animal level can be attributed to genotype as a function of the BI:BT proportion in a given individual^35,36. At the technology level, milk production is affected by the cow/calf interaction and nursing management¹⁸. In tropical systems, feeding is largely based on grazing of forage grasses, occasionally complemented by forage legumes³⁷. Factors involved in the environment level include prevalence of internal and external parasites, stress caused by high temperatures and rainfall, and forage quality and availability^1,38. In the present results, the Wood model had the best fit for characterizing the lactation curve in Siboney cattle at 240 d. The curve was one of continuous increase with a peak production of 7.71 kg at 40 d and a cumulative total production of 1,653 kg. The nonlinear mixed models exhibited the best fits based on the Akaike and Bayesian information criteria, and reduced residual variance, with increases of 11 % in the determination coefficient.