Bayesian estimation of lactation curve’s parameters in Iranian dairy cows

The graphical representation of daily milk yield during lactating period is lactation curve. There are many advantages of evaluation of the lactation curve in dairy cows, such as designing suitable breeding and management strategies for dairy cattle, genetic evaluation of dairy cows, and the prediction of the total milk yield of cow To explain the flow of milk production over the course of lactation in dairy cows, various mathematical models have been developed (Wilmink 1987; Wood 1967).There are many earlier studies that have successfully applied lactation curves in modeling the milk yield-DIM of Iranian Holstein. In all of these research the classic or Frequentist statistical methods have been used for parameter estimation and statistical inferences. The parameters of these models have not been estimated for Iranian dairy cows using Bayesian method. In the Bayesian approach, the parameters of the model are random variables, and inference is made on parameters using their posterior distributions (in some cases without the assumption of normality of the studied data) and need for huge data (Iqbal et al. 2019).This study was Bayesian estimation of the parameters of Wood, Milkbot, Gompertz, Dijkstra, Cobby and Le Du, Von Bertalanffy, Brody and Logistic mathematical models for the lactation curve. This was done using 30618, 30685 and 30627 days in milk (DIM) records of milk yield, fat percentage and protein percentage, respectively. These records were related to the days of 5 to 305 days of the first lactation period of Iranian Holstein cows with 3685 cows out of 350 herds. These data have been collected by National Breeding Center and improvement of animal production of Iran in registered herds of Holstein cows.

Data were initially adjusted for significant fixed effects that were herd-test-date (HTD) and age at calving. In order to be accurate and due to the high computational needs of Bayesian method, in this study, 30618, 30685 and 30627 records of milk production test day, fat percentage and protein percentage related to the first lactation period belonged to 3685 Holstein cows out of 350, respectively. For test-day records of milk yields, an outlier control assessment was conducted. 99.73 % of the observations accounted for data within the μ ± 3 standard deviations range. Records outside this range have been taken into consideration outliers (Junior et al. 2018). Bayesian inference was used to estimate the subsequent distribution of both unknown parameters of each lactation curve model, using a mixed nonlinear model to which the random lactation effect of each cow was added to account for the individual lactation curve of each cow. The parameters of lactation curves were estimated using the cow's test day records in the MCMC procedure of SAS software. To consider the individual lactation curve of each cow, the effect of each cow was used as a random effect in all nonlinear models. For sampling of the posterior distribution of parameters, the Monte Carlo Markov chain sampling algorithm was obtained by considering the burn-in period, sampling interval (thin) and number of cycles of 150,000, 100 iterations and 400,000 cycles for each, respectively. The effective sample size and Geweke detection test was used to evaluate the sampled amount and calculate the convergence indices. Comparisons of the models were made based on the Deviance Information Criterion (DIC). To determine the significance of the difference between the two models, if the difference between the two nonlinear models was less than 10, the DIC index was considered non-significant and if it was 10 or more, it was considered significant (Iqbal et al., 2019).

By examining the statistics of functions and convergence indices, Brody function for milk production and Wood function for fat percentage and protein percentage showed the most suitable model and better fit. The Convergence has been achieved according to Geweke test (Geweke 1992) statistic P-value. In the protein percentage fitting curve, Wood model was selected despite the significance of parameter a in the Geweke test (P <0.05) due to the appropriate trace plot. Effective sample size and the simulation sample size (4000) indicate convergence in all parameters. The differences between DIC values were found greater than 10 points for most of the cases, indicating high significant difference between the fitted models. In the results, the estimated parameters of the Brody function for milk production were 37.219, 0.544 and 0.084 for a, b and c, respectively. For milk fat percentage and milk protein percentage, the estimated parameters of Wood function were 4.29 and 3.53 for a, -0.08 and -0.04 for b, -0.0008 and -0.0004 for c, respectively. For the first 30 days, a modified gamma function gave the best fit for the first lactation (Sherchland et al. 1995).The Wood curve was superior to the Gompertz function in fitting the data and, hence, it was used for biological inference (Hansen et al. 2012).In another study by Bangar and Verma (2017) to compare four nonlinear models, Wood, quadratic model, mixed logarithmic model and Wilmink to fit the shape of lactation curves of milk production and its production traits in Gir crossbred cows in India, Wood model as the best model for Fitting of milk production data and production traits were introduced.

Bayesian method can be used in modeling complex nonlinear functions for lactation curves, especially with a small amount of data. Biological interpretation of these parameters makes it possible to use these estimates in a selection index to genetically modify the lactation curve. This shows that Brody and Wood functions, respectively, as the best model in predicting milk production and production traits in Iranian Holstein cows in the first lactation period.