Optimized Animal Models for the Genetic Evaluation of Conformation Traits, Milking Ease, and Milking Temperament in Dairy Gir Cattle

TL;DR

This study improves genetic evaluation models for conformation and milking traits in Brazilian Dairy Gir cattle to boost productivity and economic returns.

Contribution

The study identifies optimal statistical models for genetic evaluation of conformation and milking traits in Dairy Gir cattle.

Findings

Models fitting only significant fixed effects and treating contemporary groups as random effects showed better performance.

Linear models outperformed threshold models for most traits.

Using parsimonious models reduces estimation error and improves genetic evaluation accuracy.

Abstract

Enhancing the genetic evaluation of important traits can help increase productivity in dairy animals. This study focused on two groups of traits evaluated in Brazilian Dairy Gir cattle. Conformation traits, which describe the animal’s body and are linked to key production and health traits. Milking traits, which reflect the animals’ capacity to be milked and are especially important for zebu cattle, known for their responsiveness to the milking process. These traits are easy to measure, have moderate heritability, and are assessed early in the first lactation. The first step is to improve the statistical models used to estimate genetic parameters and breeding values. In this study, different models were tested to identify which provided the best goodness-of-fit, leading to the most suitable approach for evaluating conformation and milking traits. As a result, fitting only significant fixed effects in linear models and treating contemporary groups as random effects can be more advantageous for the genetic evaluation of conformation and milking traits in Dairy Gir cattle in Brazil. Implementing these findings in the genetic evaluation of the Dairy Gir breed can enhance genetic gains and ultimately improve farmers’ economic returns.

This study aimed to evaluate four different models for the genetic evaluation of sixteen conformation traits, milking ease, and milking temperament in Dairy Gir cattle. The models vary based on whether they include only statistically significant fixed effects or all recorded effects, along with contemporary groups (CGs) treated as fixed or random effects. Categorical traits were also analyzed using a threshold model. The adjusted R-squared (Radj2) was used to compare the goodness-of-fit of the linear models. Spearman’s rank correlation and the average accuracy of bull estimated breeding values (EBVs) with at least 20 phenotyped daughters were compared. Models fitting CG as a random effect performed better based on their Radj2 values and had a greater average accuracy of EBVs for most traits. Spearman’s rank correlation coefficients indicated low to medium bull EBV re-ranking between most of the models. The linear models performed better than threshold models for almost all traits. When possible, more parsimonious linear models fitting only significant fixed effects should be used to reduce the standard error of estimation. Additionally, fitting CGs as a random effect seems more beneficial for the genetic evaluation of conformation and milking traits in Dairy Gir cattle in Brazil.

Full text

1. Introduction

In the Brazilian Dairy Gir National Breeding Program, eighteen linear traits are routinely evaluated using scores (from 1 to 9), centimeters, or angles. These traits include 16 body conformation traits, milking ease, and milking temperament. Conformation traits can be important indicators of herd life and general health [1,2]. Studies have shown that rump conformation traits impact calving ease, uterine and vaginal prolapse, and udder fixation [3,4]. Mammary system traits, mainly the anterior udder, udder depth, and teat length, are related to mastitis and susceptibility to trauma [5]. Traits such as milking ease and temperament, which evaluate animal behaviour during milking, are essential for effective daily dairy management [6]. Zebu breeds, such as Gir cattle, are generally more reactive than European breeds. This reactivity can lead to injuries for themselves and other animals, cause accidents, and prolong the milking process [7,8].

The estimation of breeding values and genetic parameters for these traits relies on optimized statistical models that explain the variability and produce reliable estimates. Some conformation and milking traits are categorical traits that are evaluated using scores. Hence, evaluating whether threshold models could fit these traits better than linear models is crucial. Theoretically, threshold models are preferable to linear models for categorical data [9], as linear models assume that residuals are normally distributed, which is typically not the case with score-based traits [10].

In dairy cattle, genetic evaluation models typically compare animal performance within groups, with the assumption that animals in the same group are raised in similar environments. These contemporary groups (CGs) can be defined by the combination of herd–year–season, herd–year, or their interactions [11]. Since 1973, when Henderson introduced the best linear unbiased prediction (BLUP) for the sire model evaluation [12], CGs have been treated as a fixed effect to avoid bias in situations where good sires are used in well-managed herds. However, this problem was mitigated with the implementation of animal models in evaluation systems, as CGs are genetically connected not only through sires but also through dams [13,14].

Given the importance of these traits in selecting productive and durable animals, understanding how poor statistical models can negatively influence genetic evaluations is crucial. Therefore, this study aimed to determine which models would perform better in the genetic evaluation of conformation traits, milking ease, and milking temperament in Dairy Gir cattle by testing the significance of fixed environmental effects to be included in the model and assessing the inclusion of CG as random effect in the model.

2. Materials and Methods

2.1. Data Description

The dataset used in this study was provided by Embrapa Dairy Cattle (Juiz de Fora, MG, Brazil), which is part of the Brazilian Dairy Gir National Breeding Program. The records were collected by the Brazilian Association of Dairy Gir Breeders and the Brazilian Association of Zebu Breeders between 1992 and 2018 from herds located in the southeastern region of Brazil.

The traits evaluated in this study included 16 conformation traits, milking ease, and milking temperament (Table 1). Stature (STA), heart girth (HG), body length (BL), navel length (NL), rump length (RL), pin width (PW), hook width (HW), teat length (TL), and teat diameter (TD) were measured in centimeters. Foot angle (FA), rear legs–side view (LSV), rear legs–rear view (LRV), fore udder attachment (FUA), rear udder width (RUW), udder depth (UD), milking ease (ME), and milking temperament (TEM) were recorded using a visual score (9-point scale). Rump angle (RA) was recorded in degrees using a slope inclinometer protractor or, in some cases, approximated from a visual linear score. Details on the traits were previously described by Panetto et al. [15].

The raw dataset contained 10,956 linear-type evaluations performed in the early stage of lactation on 7913 Gir cows. The number of evaluations per cow ranged from one to four, with 6153 cows being evaluated only once. The mean age of the cows at the evaluation was 1820 ± 788 days. Cows were born between 1976 and 2015 and calved between 1991 and 2018. The pedigree file used in the analysis for all traits consisted of 18,402 animals up to 15 generations.

2.2. Model Description

As the dataset included repeated measures, four alternative animal repeatability models were used (M1, M2, M3, and M4), as follows:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the observation vector for the linear conformation traits, milking ease or milking temperament; $[eqn]$ is the incidence matrix for fixed effects; $[eqn]$ is the incidence matrix for the animal additive genetic effects; $[eqn]$ is the incidence matrix for the permanent environment effects of the animal; $[eqn]$ is the incidence matrix for the CG effect; $[eqn]$ is the vector of fixed environmental effects; $[eqn]$ is the vector of random additive genetic effects; $[eqn]$ is the vector of random animal permanent environment effects; $[eqn]$ is the vector of random CG (defined by animals in the same herd and year of evaluation) effects; and $[eqn]$ is the vector of random residual effects. For models M1 and M2, CG was fit as a fixed effect in $[eqn]$ .

For the models above, $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ were assumed, where $[eqn]$ is the additive relationship matrix; $[eqn]$ is the additive genetic variance; $[eqn]$ is an identity matrix; $[eqn]$ is the permanent environmental variance; $[eqn]$ is the contemporary group variance; and $[eqn]$ is the residual variance. The genetic, permanent environment, contemporary group and residual effects were assumed to be independent.

Model M2 was the most complete model for all traits with respect to potential environmental effects, including all available environmental effects as fixed effects, except for the permanent environmental effect. Model M4 was similar to model M2, but with the CG effect treated as a random effect to exploit the use of the inter CG information and allow for the use of small-sized CG in the genetic evaluation. Model M1 was similar to model M2, but included just fixed environmental effects that were statistically significant to assess the use of more parsimonious models. Model M3 was similar to model M1, included only fixed environmental effects that were statistically significant to assess the use of more parsimonious models, but always included CG as a random effect. More details of the models are given below.

2.3. Models M2 and M4

The fixed effects considered in M2 and M4 included all recorded environmental effects, regardless of their statistical significance. For M2, CG (herd (H), year of evaluation (YE), and their interaction), season of evaluation (S), conformation trait evaluator (E), diet (D), and age of the dam as covariate (linear, L, and quadratic, Q, effects) were included as fixed effects. The random effects in M2 included animal and permanent environment effects. The fixed effects included in M4 were S, E, D, and the age of the dam as a covariate (with L and Q effects). The effects of the animal, the permanent environment, and the contemporary group were included as random effects.

For the eight traits measured with scores, the two models above were run as both linear models and threshold models. However, only linear models were fit for all other traits. See Section 2.5 for more details.

2.4. Models M1 and M3

A general linear model (GLM) was applied, using the R software [16] as a preliminary step to define which fixed effects would be included in models M1 and M3. Fixed effects were selected based on statistical significance (F test, p < 0.05) and biological relevance (Table 2). For M1, the fixed effects tested were H, YE, the interaction between H and YE, S, E, D, and the age of the dam as a covariate (L and Q effects). For M3, the same effects were tested, except for H, YE, and their interaction, which were fit as a random CG effect.

The random effects considered in M1 were animal and permanent environment effects. For M3, in addition to the CG, random animal and permanent environment effects were also included.

For eight traits measured with scores, the two models above were run as both linear models and threshold models. Therefore, only linear models were fit for all other traits. See Section 2.5 for more details.

2.5. Data Editing

Data that were outside of possible biological ranges for traits evaluated in cm and degrees were removed. Evaluators (between 13 and 23 evaluators were removed according to each trait) and CGs (between 542 and 570 CGs were removed according to each trait) with fewer than five records were removed from the analysis. The connectedness between CGs was assessed using the AMC program [17], in which the degree of connectedness is measured through the existence of genetic links attributed to a common ancestor. The program used ten or more genetic links to consider a CG connected, and any unconnected CG was removed (between 4 and 8 CGs were removed according to each trait). The descriptive statistics of this final dataset are provided in Table 1.

Due to differences in the models for BL, LSV, LRV, and TD, in which CG was not included either as a fixed effect or as a random effect in M1, because the CG effect was not significant for these traits, a second dataset was created, which was not edited based on CG size and connectedness. Disconnected CGs or those with fewer than five records within each group were not removed for these traits in M1 (Table 3). This dataset (dataset 2) was used exclusively for the analysis of the traits BL, LSV, LRV, and TD under M1. The remaining trait for M1 and all traits under M2, M3, and M4 were evaluated using dataset 1. Despite the differences in the number of animals used in both datasets, the means and standard deviations were similar across the datasets (Table 1 and Table 3).

2.6. Data Distribution and Transformation

Skewness and kurtosis tests were performed in R software (version 4.5.2) [16] using the R package “moments” (version 0.14.1) [18] to verify distribution properties of the categorical traits. Several data transformations were tested when substantial deviations from symmetry were observed, including log transformation, square root transformation, cube root transformation, Yeo–Johnson transformation, and ordered quantile normalization transformation. These transformations were performed using base R or the “fitdistrplus” (version 1.2-4) [19] and “bestNormalize” (version 1.9.2) [20] packages. The ordered quantile normalization method [20] was chosen to transform the score records because it yielded the best skewness and kurtosis values after the transformation. The transformed data were used for the analyses of the linear models. For the threshold models, the original data were used for the analyses.

2.7. Genetic Parameter Estimates

Variance components were estimated using the average information restricted maximum likelihood (AIREML) method [21] under a single-trait repeatability animal linear model and threshold model. All analyses were conducted using ASREML 4.1 software [22].

For the multinomial threshold analyses, the !MULTINOMIAL qualifier and the !LOGIT function were used in ASREML 4.1 software. The logit link function is as follows:

$$[eqn]$$

where $[eqn]$ is the mean on the data scale and $[eqn]$ is the linear predictor on the underlying scale. On the logit scale, the residual variance is as follows:

$$[eqn]$$

Heritability ( $[eqn]$ ) and repeatability ( $[eqn]$ ) estimates for the threshold models were estimated as follows:

$$[eqn]$$

where $[eqn]$ is the additive genetic variance and $[eqn]$ is the permanent environmental variance [22]. The heritability ( $[eqn]$ ) and repeatability ( $[eqn]$ ) of the linear models were calculated as follows:

$$[eqn]$$

where $[eqn]$ is the residual variance.

For models that include CG as a random effect (M3 and M4), the denominator of the $[eqn]$ and $[eqn]$ equations also include the CG variance ( $[eqn]$ ).

2.8. Model Evaluation Criteria

First, models with all fixed effects included were compared to models that included only significant effects (M1 versus M2). Then, the inclusion of CG as fixed or random effects was compared (M1 versus M3 or M2 versus M4). In the first comparison, the structure of the fixed effects changes between models, and in the second, the structure of both the fixed and random effects changes at the same time, making the use of an appropriate parameter to compare the models difficult. As noted by Verbyla [23], models that vary in the structure of their fixed effects cannot have their restricted (residual) likelihood used to compare the models. Instead, the full and maximum likelihoods should be used. Because ASREML 4.1 software uses and outputs the restricted maximum likelihood rather than the maximum likelihood of the models, in this study, we chose to use the adjusted R-squared ( $[eqn]$ ) value, which is a measurement that adjusts the goodness-of-fit to the number of fixed effects considered in the model. $[eqn]$ was calculated for each model and trait individually as follows:

$$[eqn]$$

where $[eqn]$ is the number of records; $[eqn]$ is the correlation between $[eqn]$ and $[eqn]$ squared, calculated as $[eqn]$ ; and r( $[eqn]$ ) is the sum of the degrees of freedom. A higher value of $[eqn]$ indicates that the model explains more of the phenotype. Even thought, $[eqn]$ is not the best criterion to compare mixed models, as it does not take into account the partition of the variance across random effects, it would be a better alternative than using restricted likelihood based criteria in the case of alternative mixed models that vary in their fixed effects.

2.9. Impacts of the Models on the Animal Evaluation

The average accuracy of the estimated breeding value (EBV) of bulls with at least 20 phenotyped daughters for each trait was calculated for all the models. The accuracy of the EBV for each bull was estimated as follows:

$$[eqn]$$

where $[eqn]$ is the standard error of prediction of the EBV of the ith animal as estimated by ASREML; $[eqn]$ is the inbreeding coefficient of the ith animal; and $[eqn]$ is the population additive genetic variance [24]. A paired Student’s t-test was used to test differences (p < 0.05) between the average EBV accuracy of bulls with at least 20 phenotyped daughters between the alternative models for each trait.

Additionally, Spearman’s rank correlation coefficient was computed between the EBV for the same bulls used to calculate the average EBV accuracies to assess the impact of using different models in the ranking of the bulls used as parents in the next generations. The comparisons were based on linear model M1 versus M2, M1 versus M3, M2 versus M4, linear model M3 versus threshold model M3, and linear model M4 versus threshold model M4.

3. Results and Discussion

Optimizing models to be used in the genetic evaluation of important economic traits is necessary since they can improve the estimation of genetic parameters, breeding values, and accuracies of predictions. An optimized model can impact the selection of animals to be the sires and dams of the next generations and, consequently, can lead to greater genetic gain and increased profitability. This study compared different models for the genetic evaluation of conformation traits, milking ease, and milking temperament in Dairy Gir cattle to identify the most suitable model. The models differed in their fixed effect structure and the inclusion of contemporary groups as fixed or random effects. Additionally, linear and multinomial threshold models were compared for the genetic evaluation of the categorical traits.

The results and discussion that follow present genetic parameter estimates across all traits and models, with a comparison of goodness of fit, average EBV accuracies, and bulls’ EBV rankings of linear models for groups of traits, followed by a comparison between linear and threshold models for categorical traits, and a summary of results and limitations of this study.

3.1. Genetic Parameter Estimates

The estimates of additive genetic, permanent environmental, and residual variances for all the evaluated models are shown in Table 4. Heritability and repeatability estimates are presented in Table 5 for both linear and threshold models. The estimated standard errors for the variance, heritability and repeatability estimates are presented in Table S1 and Table S2 (Supplementary Material), respectively.

The estimates of genetic, permanent environmental, and residual variances were similar between the linear models for almost all the traits (Table 4). Body length, rear legs–side view, rear legs–rear view, and teat diameter had a smaller estimated residual variance for linear models M2, M3, and M4 than for linear model M1. This result may be because for these traits, CG was not included in the linear model M1 but it was included in the linear models M2, M3, and M4 as a fixed or random effect, which could have reduce the residual variance.

The estimates of additive genetic variance for teat diameter and body length were higher for linear model M1 than for M2, M3, and M4. However, estimates of permanent environmental variance were higher for linear models M2, M3, and M4 than for linear model M1. The linear M2, M3, and M4 models included CG (herd–year evaluation) for both traits, either as a fixed or random effect, whereas linear model M1 included only the year of evaluation. The difference in the genetic additive variance is reflected in a higher heritability estimate for linear model M1 for these traits. However, the repeatability estimates were similar between linear models M1 and M2 (Table 5).

Small differences in the repeatability estimates for the hook width between linear models M1 and M2 and linear models M3 and M4 were observed, which could be attributed to the inclusion of the evaluator as a fixed effect in M2 and M4.

The estimates of genetic additive and permanent environment variance were very similar for the threshold models M1 and M2 for all the traits, except for the rear legs–side view, in which M1 showed smaller estimates of genetic additive and permanent environment variance (0.47 and 1.68, respectively) than did M2 (0.71 and 2.50, respectively) (Table 4). These differences may be attributed to the fixed effects (CG and quadratic effect of the dam’s age) added to the second model. When comparing threshold models M1 versus M3 and M2 versus M4, the estimates of genetic additive and permanent environment variance were smaller for M3 and M4 for almost all traits, except for rear legs–side view and rear udder width, which is expected since one more random effect (CG) was included in threshold models M3 and M4.

The estimates of genetic additive variance (Table 4) and heritability (Table 5) for rear legs–rear view were low and not reliable for all four models when both linear and threshold models were used, which may be due to the inability to record this trait correctly in the field or the poor definition of the trait. Traits related to feet and legs are highly influenced by the environment and management [25,26], which explains the low heritability estimates. Heritability estimates in the literature for rear legs–rear view are scarce; however, Duru et al. [27] reported a heritability of 0.11 (±0.18), and Ptak et al. [28] reported a heritability equal to 0.09 (±0.02) for this trait.

When threshold models M1 and M2 were compared, heritability and repeatability estimates were similar between models (Table 5). As expected, thresholds M3 and M4, which included CG as a random effect, yielded estimates of heritability and repeatability smaller than those of M1 and M2, which included CG as a fixed effect (Table 5). This result is expected because the inclusion of CG variance in the estimated phenotypic variance leads to a higher denominator in the fraction, resulting in a smaller heritability and repeatability estimate, which should be interpreted as estimates across CGs in contrast to the estimates within CG, which are obtained when CG is treated as a fixed effect and, therefore, does not contribute to the estimated phenotypic variance.

Heritability estimates for conformation traits, milking ease, and milking temperament vary in the literature, which could be due to several factors, including the statistical model, breed studied, fixed and random effects considered in the models, and data editing, among others [29]. However, the heritability estimates in this study were generally similar to the results in the literature [27,30,31,32,33,34,35,36].

3.2. Heart Girth, Rump Length, Teat Length, and Navel Length

For heart girth, rump length, teat length, and navel length, all tested fixed effects were significant (p < 0.05), resulting in the same models for M1 and M2 (Table 2) and the same values of $[eqn]$ . Consequently, M3 and M4 were also the same models with the same $[eqn]$ values. When the linear models with CGs fitted as a fixed effect (M1 and M2) were compared with the linear models with CG fitted as a random effect (M3 and M4), the linear models M3 and M4 showed a better fit based on their higher $[eqn]$ values (Table 6). For these traits, the average EBV accuracies were the same between M3 and M4 (Table 7). The differences in average EBV accuracy of bulls with at least 20 daughters with records between M1 and M3 or M2 and M4 were statistically significant (p < 0.05) for all compared traits (Table S3 in the Supplementary Material). The EBV Spearman rank correlation coefficients between M1 and M3 (or M2 and M4) for these traits ranged between 0.98 and 0.99 (Table 8), which indicates that changing the CG from fixed to random effect had little impact on the ranking of bulls for these traits.

3.3. Stature, Pin Width, Rump Angle, Udder Depth, Milking Ease, and Milking Temperament

Stature, pin width, rump angle, udder depth, milking ease, and milking temperament showed the same $[eqn]$ values between linear models M1 and M2 (Table 6), which could indicate that differences in the structure of fixed effects for these traits were not sufficient to affect their $[eqn]$ values. The $[eqn]$ values between linear models M3 and M4 for these traits were the same. However, when the linear models M1 and M2 were compared with M3 and M4, respectively, the models that included CG as a random effect showed higher values of $[eqn]$ , indicating a better fit of the linear models M3 and M4 with CG as a random effect. Although the $[eqn]$ values were the same between M3 and M4, fitting nonsignificant fixed effects in a model can lead to increased standard error of the estimates. When the average EBV accuracies for these traits were compared between linear models M1 versus M3 and M2 versus M4, linear M3 and M4 presented higher averages for all traits, except udder depth and milking temperament, which presented similar values (Table 7). The differences in average EBV accuracy for bulls with at least 20 daughters with records between M1 and M3 and between M2 and M4 were statistically significant (p < 0.05) for all compared traits, except for udder depth (Table S3 in the Supplementary Material). Similarly, the EBV Spearman rank correlation coefficient between linear M1 and M2 was greater than 0.99 (Table 8) for these traits, indicating that differences in the fixed effects fitted in the model led to little or no reranking between both models, also supporting the use of more parsimonious models. The EBV Spearman rank correlation coefficients between linear models M1 and M3 and between M2 and M4 were similar and ranged from 0.97 to 0.99 (Table 8), indicating that just a few bulls were reranked, especially for pin width and udder depth, when CG was considered a random effect.

3.4. Rear Legs–Side and –Rear Views

For the rear legs–side and –rear views, linear model M1 showed higher values of $[eqn]$ than the linear model M2. When linear model M1 was compared to linear model M3, linear model M3 showed higher values of $[eqn]$ (Table 6). This result could indicate that considering CG as a random effect could have led to better model performance for the rear legs–side and –rear views. The average EBV accuracies for these traits were greater for linear model M1 (rear legs–side view) and M4 (rear legs–rear view) (Table 7). A higher average EBV accuracy for linear model M1 for the rear legs–side view could be explained by a slightly higher additive genetic variance estimated for this model. Differences in average EBV accuracy for bulls with at least 20 daughters with records were statistically significant (p < 0.05) for rear legs–side and –rear views between M1 and M2, M1 and M3, and M2 and M4 (Table S3 in the Supplementary Material). The EBV Spearman rank correlation coefficient between linear models M1 and M2 indicated greater reranking when the fixed effects changed between models for rear legs–side view (0.91) than for rear legs–rear view (0.68). The removal of CG from the models, due to its lack of significance, impacted the bulls’ rankings, as well as the genetic parameter estimates, as discussed previously. A large reranking was also observed when linear models M1 and M3 were compared for rear legs–rear view (0.80, Table 8) and again a much smaller re-ranking was found for rear legs–side view (0.94). Changing the CG from fixed to random effect impacted much less the ranking (0.95 and 0.99) when comparing linear models M2 and M4 for rear legs–rear and –side views, respectively.

3.5. Body Length, Hook Width, Foot Angle, Fore Udder Attachment, Rear Udder Width, and Teat Diameter

The body length, hook width, foot angle, fore udder attachment, rear udder width, and teat diameter showed higher values of $[eqn]$ (Table 6) for linear model M2 than for linear model M1. When the linear model M2 was compared to the linear model M4, the linear model M4 presented higher values of $[eqn]$ , indicating a better fit of this model for these traits. The average EBV accuracies were higher for linear model M4 for foot angle and hook width. For body length and teat diameter, linear model M1 had the highest average EBV accuracies, and for fore udder attachment and rear udder width, linear model M2 showed the highest average EBV accuracies (Table 7). The highest average EBV accuracies found for linear models M1 and M2 for body length and teat diameter and for fore udder attachment and rear udder width could be explained by higher additive genetic variance estimated for these models and traits, respectively. Differences in average EBV accuracy for bulls with at least 20 daughters with records were statistically significant (p < 0.05) for all compared traits between M1 and M2, M1 and M3, and M2 and M4 (Table S3 in the Supplementary Material). Removing CG as an effect from linear model M1 for body length and teat diameter may explain the moderately high reranking observed between linear models M1 and M2 (0.88 and 0.89) and linear models M1 and M3 (0.92) (Table 8). Little reranking was observed between linear models M2 and M4 (0.98 to 0.95) for body length, hook width, foot angle, fore udder attachment, rear udder width, and teat diameter.

3.6. Linear and Threshold Models

A comparison of the $[eqn]$ values between the linear and threshold models is not possible since they are in different scales. However, an evaluation of the impacts of both models on average EBV accuracies and the bull rankings is possible. When comparing the linear model with higher $[eqn]$ values (M3: rear legs–side and –rear views; or M4: foot angle, fore udder attachment, rear udder width, udder depth, milking ease, and milking temperament) to their respective thresholds for the categorical traits, the average EBV accuracies were higher for the linear models (M3 and M4) than for the threshold models (M3 and M4), except for milking ease (Table 7). The differences in average EBV accuracy for bulls with at least 20 daughters with records between the linear and threshold models were statistically significant (p < 0.05) for all compared traits, except for milking temperament for linear model M3 versus threshold model M3 and milking ease for linear model M4 versus threshold model M4 (Table S3 in the Supplementary Material). The EBV Spearman rank correlation coefficients were greater than 0.94 for all traits (Table 8) indicating a low degree of reranking between linear and threshold models. The similar results found between linear and threshold models can also be attributed to the transformation of the categorical data to fit normality when linear models were applied. The use of non-normal distributed data can lead to biased genetic parameter estimates and improper animal ranking.

Meijering [37] reported no advantage of using threshold models over linear models, although threshold models were theoretically a better option. Vanderick et al. [38] observed that the threshold model had better goodness of fit than linear models. However, according to the authors, no clear advantage was found in terms of predictive ability, and linear models would be more suitable and practical for application in genetic evaluation. Weller et al. [39] also identified several advantages of using threshold models over linear models, and even though threshold models estimate larger variance components, the rank correlation estimated between the two models was greater than 0.90. Threshold models can theoretically be the best fit for categorical traits. However, this study revealed that the lower average EBV accuracies and minimal changes in the animals’ rankings for most of the studied traits may not justify the implementation of these models in a genetic evaluation.

3.7. Summary and Study Limitations

In this study, models including all fixed effects (M2) performed better or the same as models that included only statistically significant fixed effects (M1) based on their $[eqn]$ values for most traits, except rear legs–side and –rear views. The inclusion or exclusion of fixed effects in a model needs to be performed carefully. The inclusion of nonsignificant fixed effects can lead to overparameterization of the statistical models and increase the standard error of the estimates. However, when important effects are not included in the model can lead to biases in the evaluation. For traits with the same $[eqn]$ values between linear models M1 and M2, using a model that includes only significant fixed effects may be recommended to reduce the standard error of the estimates.

Models incorporating CG as a random effect showed better fitting across all assessed traits in this study, as indicated by their higher $[eqn]$ values. As suggested by Schaeffer [13], when an animal model is used, CG should be treated as a random effect, whereas CG was commonly treated as a fixed effect in the past based on Henderson’s selection bias theory under a sire model. Fitting CG as a random effect can result in other benefits. For instance, when CG is treated as a random effect, the loss of information from the removal of CGs with a small number of individuals (e.g., less than 5) would be avoided [40].

A limitation of this study is the use of $[eqn]$ to compare linear mixed models. This limitation could be resolved if the models’ maximum likelihood values, instead of restricted maximum likelihood values, were available, which could be used to compare models that differ in their fixed effects.

4. Conclusions

This study provides valuable information for the optimization of models for future applications in the genetic evaluation of conformation traits, milking ease, and milking temperament in Dairy Gir cattle in Brazil. When linear and threshold models were compared, no advantages were observed for the threshold model. The inclusion of only significant environmental fixed effects in the model should be considered instead of all recorded environmental effects for some traits to avoid overparameterization and increased standard errors of estimation. Finally, contemporary group should be included as a random effect in the animal models since this modelling provided better model fitting and higher EBV accuracies and little bull EBV re-ranking.