| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1 Department of Animal and Aquacultural Sciences, Norwegian University of Life Sciences, N-1432 Ås, Norway
2 GENO Breeding and AI Association, N-1432 Ås, Norway
3 Department of Dairy Science, University of Wisconsin, Madison 53706
Corresponding author: I. M. Andersen-Ranberg; e-mail: ina.ranberg{at}umb.no.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
2 statistic based on differences between observed and predicted outcomes for NR56 in a separate dataset. Comparison was also with respect to ranking of sires and correlations between sire posterior means (TL model) and PTA (LL model). We found very small differences in ability to predict NR56 between the 2 bivariate models, regardless of the dataset used. Correlations between sire posterior means (TL) and sire PTA (LL) and rank correlations between sire evaluations were all >0.98 in the 3 datasets. At present, the LL model is preferred for sire evaluations of NR56 and CFI in NRF. This is because the LL model is less computationally demanding and more robust with respect to the structure of the data than TL.
Key Words: female fertility • genetic parameter • model comparison • bivariate threshold model
Abbreviation key: CFI = interval from calving to first insemination, LL = linear-linear, MCMC = Markov chain Monte Carlo, NR56 = 56-d nonreturn rate, NRF = Norwegian Red, TL = threshold-linear.
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Roxström et al. (2001) and Andersen-Ranberg et al. (2005) found relatively high correlations between fertility in heifers and in first-lactation cows, varying from 0.54 to 0.67. However, because the estimated correlations are not 1, the 2 traits are probably not genetically the same trait (Thaller, 1997). Andersen-Ranberg (2005) concluded that the interval from calving to first insemination (CFI) should be included in the fertility index for NRF and reported low correlation between NR56 in first-lactation cows CFI in first-lactation cows. However, in a simulation study, Gates et al. (1999) concluded that genetic correlations involving categorical traits might be underestimated in linear sire models. This leads us to further investigation on the traits NR56 and CFI.
Estimates of heritability for fertility traits are typically low, ranging from 1 to 5% when linear models have been used for statistical analysis (Distl, 1982; Weigel and Rekaya, 2000; Roxström et al., 2001; Wall et al., 2003). A standard linear model assumes that the observed binary response follows a normal distribution, and variance component estimates obtained with such models may be misleading (Hoeschele et al., 1987). For example, heritability is frequency dependent, as the variance of the distribution depends on the mean (Dempster and Lerner, 1950), which implies that observed genetic change may be inconsistent with what would be expected from heritability estimates based on a linear model. Statistically, it is inappropriate to use standard linear model methodology for analyzing categorical response data (Gianola, 1982; Agresti, 1996).
From a genetic point of view, the threshold-liability (TL) model (Wright, 1934; Dempster and Lerner, 1950; Falconer and Mackay, 1996) has the theoretical appeal of producing estimates of parameters that are interpretable on an underlying continuous (liability) scale where gene substitutions are supposed to take place. This model was discussed by several authors (Thompson, 1979; Gianola 1982), and the threshold model was developed in practice in the early 1980s (Gianola and Foulley, 1983; Harville and Mee, 1984; Gilmour et al., 1985).
Because of advances in Markov chain Monte Carlo (MCMC) methods, such as Gibbs sampling, it has become possible to carry out exact (within the limits of Monte Carlo error) Bayesian analysis of large linear and nonlinear hierarchical models (Gilks et al., 1996; Wang, 1998; Sorensen and Gianola, 2002). The MCMC methods avoid the need for numerical integration by taking repeated samples from the posterior distributions of interest. A Bayesian MCMC implementation of the threshold model in a quantitative genetic context can be found in Sorensen et al. (1995).
Despite the theoretical appeal of the threshold model, there is an issue of how much difference it can make in a given animal breeding setting. For example, the simulation study of Meijering and Gianola (1985) compared sire evaluations obtained with linear and threshold models and found an advantage of the latter in situations with moderate or high heritability in the liability scale and low incidence (<25%) of the trait.
Studies with field data have reported very high correlations between sire PTA obtained from linear and threshold models (Weller et al., 1988; Heringstad et al., 2003). Conversely, Varona et al. (1999) found that a bivariate TL model for birth weight and calving ease had better predictions than the corresponding bivariate linear-linear (LL) model. It seems sensible to examine whether or not a threshold model can enhance genetic improvement of fertility of dairy cattle, as most fertility-related traits are recorded on a categorical scale.
In this study, heritability was estimated for liability to NR56 and for CFI, as well as their genetic correlation, in first-lactation NRF cows. A bivariate model was fitted, in which NR56 was regarded as a threshold trait, and CFI was treated as Gaussian. A main objective of this paper was to compare the bivariate sire TL with a bivariate LL model in terms of estimates of genetic parameters and of sire evaluations for NR56 and CFI. We also wanted to examine how differences between models are affected by varying progeny group sizes.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Four thousand herds were randomly chosen from DATA, yielding a data subset (TEST) with 41,467 records, which was used for model comparison (predictive ability) but not for parameter estimation. Three different datasets were created with decreasing mean progeny sizes. This was done because we wanted to investigate the effect of decreasing mean progeny group using the TL vs. LL model. The subset of records in DATA, other than those in TEST, was termed DATA1. Herds were randomly deleted from DATA1, which had 110,245 records, to make 2 smaller datasets (DATA2 and DATA3) with intended average sire progeny group sizes of 100 and 50 daughters, respectively. Hence, DATA3 was a subset of DATA2, which, in turn, was contained in DATA1. The actual mean progeny group sizes were 147.8, 102.7, and 56.5 daughters per sire in DATA1, DATA2, and DATA3, respectively.
Table 1
shows summary statistics for the 4 datasets. A total of 746, 743, and 742 NRF sires were represented in DATA1, DATA2, and DATA3, respectively (Table 1), whereas 742 were in TEST. A between-sire additive relationship matrix caused by sires and maternal grandsires was built by tracing back the pedigree through sires for as many generations as possible. This resulted in pedigree files with a total of 1120, 1117, and 1116 males for DATA1, DATA2, and DATA3, respectively (Table 1). Number of sires, cows per herd, mean NR56, and mean CFI were similar in all 4 datasets. Mean NR56 was 67.0% in DATA1 and DATA3, 67.1% in DATA2, and 67.3% in TEST. Similarly, CFI was 80.7 d in DATA1, and 80.5 d in DATA2, DATA3, and TEST (Table 1).
|
), usually assumed to follow a normal distribution (Wright, 1934; Gianola, 1982; Gianola and Foulley, 1983). The binary response, yNR56, takes the value 0 (return after first insemination) if NR56 is smaller than or equal to a conceptual threshold T, and 1 (nonreturn) otherwise:
 |
The threshold T and the residual variance of liability (
e12) are not identifiable when the response is binary (Harville and Mee, 1984), so these parameters were arbitrarily set to 0 and 1, respectively. A generalized linear mixed model with a probit link (Chang, 2002) was used for NR56.
Let Pj(yNR56 = 1| ß, h, s) be the probability of the event that cow j falls in category 1 (nonreturn), given values of some parameters to be defined later. The threshold probit model postulates that the probability of nonreturn to insemination can be written as
 |
where ß, h, s are location vectors of the model for NR56, as described later; xj', zjh', and zjs' are incidence row vectors; and 
(.) is the standard normal cumulative distribution function.
Bivariate TL Model
A bivariate TL model was fit in which NR56 was treated as binary-threshold trait and CFI was assumed Gaussian, as developed by Foulley et al. (1983). The bivariate sire model for joint analysis of the underlying liability to NR56 and observed value of CFI was
 |
 |
where NR56 is the vector of all unobserved liabilities to NR56 and yCFI is the vector of observed CFI of daughters. The vector ß = [ß', ßCFI']' included effects of age at first calving in weeks and of month-year at first calving, specific to each trait. Cows in age or month-year classes with <50 records were grouped together with nearby groups to circumvent the "extreme category problem" (all responses equal to 0 or 1). There were 82, 75, and 69 age classes in DATA1, DATA2, and DATA3, respectively; 82 month-year classes in DATA1, and 79 month-year classes in DATA2 and DATA3. Further, h = [h', hCFI']' is a vector of random herd effects, and s = [s', sCFI']' is a vector of sire transmitting abilities for the 2 traits. Because the herds in Norway are small, the random herd effects were included as single herd effect without year. The incidence matrices X, Zh, and Zs were the same for NR56 and CFI, as each cow had records on each of the 2 traits as noted previously. The vectors e1 and e2 are residuals on liability to NR56 and on CFI, respectively. Residuals were assumed to be correlated within cows, and multivariate normally distributed as
 |
where 
is the residual (co)variance matrix. As stated, the residual variance of the liability to NR56 was set equal to 1, e22 is the residual variance of CFI, and e12 is the residual covariance between liability to NR56 and CFI.
A Bayesian approach using a Gibbs sampling algorithm (Sorensen et al., 1995; Sorensen and Gianola, 2002) was used for inferring parameters of interest.
Prior distributions.
The following uniform prior distribution was assigned to each of the age and month-year effects:
 |
where ßk is a given scalar element of ß. The prior distribution of herd effects was assumed to be multivariate normal:
 |
where
is the 2 x 2 (co)variance matrix between herd effects on the 2 traits and I is an H x H identity matrix; H is the number of herds appropriate to each of the 3 datasets (Table 1
). A priori effects of different herds were assumed to be independent, but effects of the same herd on the 2 traits were correlated. The vector of sire effects was assigned the multivariate normal prior distribution:
 |
where 
is the 2 x 2 (co)variance matrix between sire transmitting abilities, A is the Q x Q additive relationship matrix between sires, and Q is the number of males in the pedigree file appropriate to the dataset used (Table 1). The prior distributions of the H0 and G0 matrices were assumed to be independent scaled inverse Wishart with densities
 |
where 
H and G are the degrees of freedom (=2) and 
are scale matrices. The prior distribution of e22 was scaled inverse 2 distribution, with 2 degrees of freedom (e2) and scale parameter equal to 4 (Se22). The parameter space of the residual correlation between the 2 traits is the interval (–1, 1). This implies that the residual covariance (e12) can take any value in
. The following uniform bounded prior was adopted for e12:
 |
Posterior density.
The joint posterior density of all parameters had the form
 |
where p(y| ß, h, s, e22 , e12) is the conditional density of all NR56 (binary) and CFI observations. Pairs of records from different cows were assumed to be conditionally independent.
Draws from posterior conditional distributions of the parameters were obtained using a Gibbs sampler, following augmentation of the joint posterior densities with the unobserved liabilities to NR56 (Sorensen et al., 1995). Details of the Gibbs sampling scheme are provided in Chang (2002).
Convergence diagnostics.
Chain lengths of 50,000, 70,000, and 150,000 samples were run for DATA1, DATA2, and DATA3, respectively, with 10,000 samples discarded as burn in. The sample size and length of burn in were based on the convergence diagnostics method of Raftery and Lewis (1992) and on visual inspections of trace plots.
Sire evaluation.
Genetic evaluations (posterior means) of sire effects for NR56 were computed in the liability scale. These posterior means were transformed from the underlying liability scale to the probability (0 to 1) scale using
 |
where pi is the probability of NR56 for daughters of sire i, (.) is the standard normal cumulative distribution function, µ is the probit corresponding to the mean liability to NR56 for the population, and si is the posterior mean of sire transmitting ability of liability to NR56.
Bivariate Linear Model Analysis
For comparison purposes, the same datasets were analyzed with a bivariate linear model procedure in which NR56 was treated as Gaussian instead of binary. The explanatory structure was as described previously, the modification being that NR56 was replaced by yNR56, a vector of observations for NR56 (0 = return; 1 = nonreturn). (Co)variance components were estimated by bivariate REML, and sire PTA were predicted with BLUP (REML) using the DMU package (Jensen and Madsen, 1994).
Model Comparison
The 2 bivariate models were compared with respect to ranking of sires by calculating simple and rank correlations between sire posterior means (TL model) and PTA (LL model) within each of DATA1, DATA2, and DATA3.
We also evaluated the ability of TL and LL models to predict the probability of NR56 = 0 (or NR56 = 1) in the TEST dataset. The ratio between the number of daughters with NR56 = 0 and the number of daughters in each sire group was calculated in DATA1, DATA2, and DATA3 to arrive at sire-specific observed rates of return. Subsequently, logistic regressions of these empirical rates on either TL or LL model sire evaluations obtained from DATA1, DATA2 or DATA3 were calculated. For the TL model, the sire evaluations in the liability scale were transformed to the observable scale. The GENMOD procedure (SAS, 1999) was employed to compute logistic regression, and each empirical rate was weighted by the number of daughters of the appropriate sire. The logistic regression gives the fitted probability of a daughter not returning to insemination, conditionally on the sire evaluation from the appropriate model. As in Caraviello et al. (2004), the expected number of daughters not returning to insemination (or returning) was calculated in TEST for each sire by multiplying the total number daughters of the sire by the fitted probability of NR56 = 0 (or by the probability of NR56) from the logistic regression. The following 2 statistic (Caraviello et al., 2004) was computed for each sire:
 |
The 2 statistics were summed across sires for each of the 6 predictive situations (3 datasets x 2 models of sire evaluation; TL or LL models). The model producing the smallest sum was considered to have a better predictive ability of NR56.
|
RESULTS AND DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
). Monte Carlo estimates of posterior means (and of standard deviations) of dispersion parameters for the TL model are presented in Table 2. Corresponding REML estimates from the LL model analysis are presented in Table 3; Figure 1 shows posterior distribution of dispersion parameters for the 2 traits. The posterior distribution of the within-herd heritability of liability to NR56 in DATA1 and DATA2 was reasonably sharp and symmetric, as shown in Figure 1, with mean 0.028 and standard deviation equal to 0.004 (Table 2). The posterior distribution of heritability of liability to NR56 from DATA3 (Figure 1) was also symmetric but with higher posterior mean, 0.038 (SD = 0.007) (Table 2), which reflected the higher degree of uncertainty stemming from a smaller dataset. These point estimates of heritability of liability to NR56 are more than twice as large as those found with the LL model in the 3 datasets (Table 3). In large California Holstein herds, Weigel and Rekaya (2000) also found slightly larger heritability (1.6%) of NR56 with a univariate threshold model than with a linear model (1.4%). On the other hand, Weigel and Rekaya (2000) found, in Minnesota herds, that heritability was lower with the threshold model (3.4%) than with the linear model (4.1%), contrary to what theory for a pure random effects model would lead one to expect (Dempster and Lerner, 1950; Gianola, 1982). Other studies have found larger heritability of binary traits with threshold models than with linear models (Weller and Ron, 1992; Boichard and Manfredi, 1994; Varona et al., 1999). Estimates of heritability of NR56 from the LL model (Table 3) agreed with a previous study in NRF by Andersen-Ranberg et al. (2005) in which the point estimate was <1%.
|
|
|
and 3). For the TL model, the posterior SD was 0.003 in DATA1 and DATA2, and it was 0.005 in DATA3. The estimates of heritability are somewhat smaller than the 0.03 found by Andersen-Ranberg et al. (2005), using a larger data-set. Other studies have found heritablities in the same range, from 3 to 4.1% (Hoekstra, 1994; Pryce et al., 1997; Dechow et al., 2001). The posterior distribution of within-herd heritability of CFI was symmetric and clearly sharper in DATA1 and DATA2 than in DATA3 (Figure 1).
All point estimates of the genetic correlation between NR56 and CFI were near zero (Tables 2 and 3). However, precision was low, especially in the smallest data-set, as indicated by the wide posterior distributions in Figure 1; values ranging from –0.2 to 0.2 had appreciable posterior density. Credibility intervals include the estimate of value of the genetic correlation found by Andersen-Ranberg et al. (2005) and Hoekstra et al. (1994) of 0.08 and 0.13, respectively.
Table 4 shows the 2 statistics calculated from predicted (using DATA1 through DATA3) and observed (TEST dataset) NR56 for the TL and LL models; the model having the smallest 2 sum was interpreted as having a better ability of predicting NR56. There were very small differences in this respect for all 3 datasets. The TL model was slightly better in DATA2, but the LL model had slightly lower sums in the other 2 data-sets. Thus, the 2 comparison of the 2 models was not conclusive. The very low heritability implies that the ability of a sire posterior mean or a PTA of predicting a future NR56 must be low, as there is little variation in values of the explanatory variable; hence, the observed results could be due to sampling. Varona et al. (1999) compared a TL with a LL model for calving ease and birth weight by using a mean squared error statistic and found the TL model performed better. The advantages of using a threshold model depends on the frequency and the heritability of the trait considered (Hoeschele, 1988). Even though subclass sizes were larger in our study than in practice, we were not able to find clear differences between methods. This is in agreement with Matos et al. (1997) who found, in a study of reproductive traits, similar predictive ability for threshold and linear models. In contrast, Luo et al. (2001) compared estimation errors and variances of estimates in a simulation study linear vs. threshold model and found that linear model biased the estimates of genetic parameters for categorical traits.
|
. All correlations were >0.98 for both traits in all 3 datasets. In particular, correlations for CFI were perfect or nearly so, which was expected, as both models treated CFI in the same Gaussian manner. In NRF 7 through 10, new elite sires are selected each year, and Table 6 shows the top 10 sires for NR56, as in the LL model, and their ranking in the TL model for DATA1 through DATA3. There were few differences in ranking of sires for DATA1 and DATA2. However, in DATA3, sires ranked as 12, 25, and 23 with the TL model were included among the top 10 sires in the LL analysis. This illustrates that model choice has consequences on sire selection decisions, especially if there are few daughters per sire.
|
|
|
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
This study did not indicate any difference in predictive ability between TL and LL models, even though the TL specification seems to allocate a larger portion of the variance to the additive genetic component. From a practical point of view, a case cannot be made for using the TL model, because it is computationally more intensive and more prone to numerical difficulties. Research on alternative methods of utilizing existing data and addition of proper fertility traits is needed to further enhance the accuracy of genetic evaluations for fertility.
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Received for publication December 10, 2004. Accepted for publication February 11, 2005.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
This article has been cited by other articles:
|
|
 |
|
  C. Sun, P. Madsen, U. S. Nielsen, Y. Zhang, M. S. Lund, and G. Su Comparison between a sire model and an animal model for genetic evaluation of fertility traits in Danish Holstein population J Dairy Sci, August 1, 2009; 92(8): 4063 - 4071. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Konig, Y. M. Chang, U. U. v. Borstel, D. Gianola, and H. Simianer Genetic and Phenotypic Relationships Among Milk Urea Nitrogen, Fertility, and Milk Yield in Holstein Cows J Dairy Sci, November 1, 2008; 91(11): 4372 - 4382. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. Holtsmark, B. Heringstad, P. Madsen, and J. Odegard Genetic Relationship Between Culling, Milk Production, Fertility, and Health Traits in Norwegian Red Cows J Dairy Sci, October 1, 2008; 91(10): 4006 - 4012. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Konig, X. L. Wu, D. Gianola, B. Heringstad, and H. Simianer Exploration of Relationships Between Claw Disorders and Milk Yield in Holstein Cows via Recursive Linear and Threshold Models J Dairy Sci, January 1, 2008; 91(1): 395 - 406. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. T. Kuhn, J. L. Hutchison, and G. R. Wiggans Characterization of Holstein Heifer Fertility in the United States J Dairy Sci, December 1, 2006; 89(12): 4907 - 4920. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
B. Heringstad, I. M. Andersen-Ranberg, Y. M. Chang, and D. Gianola Short communication: Genetic analysis of nonreturn rate and mastitis in first-lactation Norwegian Red cows. J Dairy Sci, November 1, 2006; 89(11): 4420 - 4423. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
B. Heringstad, Y. M. Chang, I. M. Andersen-Ranberg, and D. Gianola Genetic analysis of number of mastitis cases and number of services to conception using a censored threshold model. J Dairy Sci, October 1, 2006; 89(10): 4042 - 4048. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
Y. M. Chang, I. M. Andersen-Ranberg, B. Heringstad, D. Gianola, and G. Klemetsdal Bivariate Analysis of Number of Services to Conception and Days Open in Norwegian Red Using a Censored Threshold-Linear Model J Dairy Sci, February 1, 2006; 89(2): 772 - 778. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
