Cumulative t-link threshold models for the genetic analysis of calving ease scores

In this study, a hierarchical threshold mixed model based on a cumulative t-link specification for the analysis of ordinal data or more, specifically, calving ease scores, was developed. The validation of this model and the Markov chain Monte Carlo (MCMC) algorithm was carried out on simulated data from normally and t4 (i.e. a t-distribution with four degrees of freedom) distributed populations using the deviance information criterion (DIC) and a pseudo Bayes factor (PBF) measure to validate recently proposed model choice criteria. The simulation study indicated that although inference on the degrees of freedom parameter is possible, MCMC mixing was problematic. Nevertheless, the DIC and PBF were validated to be satisfactory measures of model fit to data. A sire and maternal grandsire cumulative t-link model was applied to a calving ease dataset from 8847 Italian Piemontese first parity dams. The cumulative t-link model was shown to lead to posterior means of direct and maternal heritabilities (0.40 ± 0.06, 0.11 ± 0.04) and a direct maternal genetic correlation (-0.58 ± 0.15) that were not different from the corresponding posterior means of the heritabilities (0.42 ± 0.07, 0.14 ± 0.04) and the genetic correlation (-0.55 ± 0.14) inferred under the conventional cumulative probit link threshold model. Furthermore, the correlation (> 0.99) between posterior means of sire progeny merit from the two models suggested no meaningful rerankings. Nevertheless, the cumulative t-link model was decisively chosen as the better fitting model for this calving ease data using DIC and PBF.


INTRODUCTION
Data quality is an increasingly important issue for the genetic evaluation of livestock, both from a national and international perspective [13]. Breed associations and government agencies typically invoke arbitrary data quality control edits on continuously recorded production characters in order to minimize the impact of recording error, preferential treatment and/or injury/disease on predicted breeding values [5]. These edits are used in the belief that the data residuals should be normally distributed.
It has been recently demonstrated that the specification of residual distributions in linear mixed models that are heavier-tailed than normal densities may effectively mute the impact of residual outliers, particularly in situations where preferential treatment of some breedstock may be anticipated [41]. Based on the work of Lange et al. [24] and others, Stranden and Gianola [42] developed the corresponding hierarchical Bayesian models for animal breeding, using Markov chain Monte Carlo (MCMC) methods for inference. In their models, residuals are specified as either having independent (univariate) t-distributions or multivariate t-distributions within herd clusters. Outside of possibly longitudinal studies, the multivariate specification is of dubious merit [36,41,42] such that all of our subsequent discussion pertains to the univariate t-error specification only.
Auxiliary traits such as calving ease or milking speed are often subjectively scored on an ordinal scale. It might then be anticipated that data quality, including the presence of outliers, would be an issue of greater concern in these traits than more objectively measured production characters, particularly since record keeping is generally unsupervised, being the responsibility of the attending herdsperson. As one example of preferential treatment, a herdsperson may more quickly decide to assist or even surgically remove a calf from a highly valued dam. Luo et al. [25] has furthermore suggested that a decline in the diligence of data recording was partially responsible for their lower heritability estimates of calving ease relative to earlier estimates from the same Canadian Holstein population.
The cumulative probit link (CP) generalized linear mixed model, otherwise called the threshold model, is currently the most commonly used genetic evaluation model for calving ease [4,49]. MCMC methods are particularly well suited to this model since the augmentation of the joint posterior density with normally distributed underlying or latent liability variables facilitate implementations very similar to those developed for linear mixed effects models [2,39]. A cumulative t-link (CT) model has been proposed by Albert and Chib [2] for the analysis of ordinal categorical data, thereby providing greater modeling flexibility relative to the CP model. The CT model can be created by simply augmenting the joint posterior density with t-distributed rather than normally distributed underlying liability variables [18]. Since outliers on the observed categorical scale also correspond to outliers on the underlying liability scale [1], the CT model might be anticipated to be more robust to residual outliers relative to the CP model.
The objectives of this study were to validate MCMC inference of the CT generalized linear mixed (sire) model via a simulation study and to compare the fit of this model with the CP model for the quantitative genetic analysis of calving ease scores in Italian Piemontese cattle. In section 2, the CT model is constructed hierarchically. We then present a discussion of two model choice criteria that we believe are appropriate for the comparisons of the CP with the CT model in section 3. In section 4, we describe a simulation study that is used to validate posterior inference and model choice criteria for the CP and CT models, presenting the results of this study along with an application to Italian Piemontese calving ease data in section 5. We conclude with a discussion of these results in section 6.

MODEL CONSTRUCTION
Suppose that elements of the n × 1 data vector can take values in any one of C mutually exclusive ordered categories. The classical CP model for ordinal data [17] can be written as follows: where j = 1, 2, . . . , C denotes the index for categories. Also, Φ(.) denotes the standard normal cumulative distribution function, β and u are the vectors of unknown fixed and random effects, and τ = [τ 0 τ 1 . . . τ C ] is a vector of unknown threshold parameters satisfying τ 1 < τ 2 . . . < τ C with τ o = −∞ and τ C = +∞. Furthermore, x i and z i are known incidence row vectors. Latent liability variables (L = {L i } n i=1 ) can be introduced to alternatively define the same specification as in equation (1) but in two hierarchical stages: for i = 1, 2, . . . , n. Here 1(.) denotes an indicator function, which is equal to 1 when the expression in the function is true and is equal to 0 otherwise. As shown by Albert and Chib [2], and in an animal breeding context by Sorensen et al. [39], this model augmentation using L facilitates a tractable MCMC implementation.
The CT model is a simple generalization of (1), that is, for j = 1, 2, . . . , C where F v represents the cumulative density function of a standard Student t-distribution with degrees of freedom v. Note that as v → ∞, (3) → (1) such that the standard CP model is simply a special case of the CT model. Like the CP model, the CT model can also be represented as a twostage specification, with the first stage as in equation (2a) but the second stage specified as: , (4) i.e., L i is Student t-distributed with location parameter µ i = x i β + z i u, scale parameter σ 2 e > 0 and degrees of freedom v > 2 for i = 1, 2, . . . , n. In turn, equation (4) can be represented by a two-stage scale mixture of normals: Note that (5b) specifies a Gamma density with parameters v/2 and v/2, thereby having an expectation of 1. The remaining stages of our hierarchical model are characteristic of animal breeding models. We write β ∼ p(β), (6) where p(β) is a subjective prior, typically specified to be flat or vaguely informative. Furthermore, the random effects are typically characterized by a structural multivariate prior specification: Here G(ϕ) is a variance-covariance matrix that is a function of several unknown variance components or variance-covariance matrices in ϕ, depending on whether or not there are multiple sets of random effects and/or specified covariances between these sets; an example of the latter is the covariance between additive and maternal genetic effects. Furthermore, flat priors, inverted Gamma densities, inverted Wishart densities or products thereof may be specified for the prior density p(ϕ) on ϕ, depending, again, on the number of sets of random effects and whether there are any covariances thereof [21]. Finally, a prior is required for the degrees of freedom parameter v to ensure a proper joint posterior density. We use the prior: which is consistent with a vaguely informative Uniform(0,1) prior on 1/(1+v).
As with the CP models, there are identifiability issues involving elements of τ with σ 2 e such that constraints are necessary. The origin and scale are arbitrary so that, as done by others (e.g. [17]), τ 1 is set here to zero and σ 2 e to 1. We chose this parameterization such that inference on σ 2 e is not subsequently considered in this paper.
Presuming that the elements of Y are conditionally independent given β and u, we can write the joint posterior density of all unknown parameters and latent variables (L) as follows: (9) where λ = {λ i } n i=1 . An MCMC inference strategy involves determining and generating random variables from the full conditional densities (FCD) of each parameter or blocks thereof. Many of the FCD can be directly derived using results from Sorensen et al. [39] jointly with the results from Stranden and Gianola [42]. Let θ = [β u ] . It can be readily shown that the FCD of θ is multivariate normal: where The generation of individual elements θ j , j = 1, 2, . . . , p + q or blocks thereof of θ from their respective FCD is straightforward using the strategy presented by Wang et al. [48]. The FCD of individual elements of L and τ are straightforward to generate from, using results from Sorensen et al. [39]. We, however, prefered the Metropolis-Hastings and method of composition joint update of L and τ presented by Cowles [11]. She demonstrated and we have further noted in our previous applications [23] that the resulting MCMC mixing properties using this joint update are vastly superior to using separate Gibbs updates on individual elements of L and τ as outlined by Sorensen et al. [39]. A lucid exposition on Cowles' update is also provided by Johnson and Albert [22].
If some partitions of ϕ form a variance-covariance matrix, then their respective FCD can be readily shown to be inverted-Wishart [21] whereas if other partitions of ϕ involve scalar variance but no covariance components, then the FCD of each component can be shown to be inverted-gamma.
The FCD of λ i can be shown to be: (12) that is, the kernel in (12) specifies that the distribution to be Here λ −i denotes all elements of λ = {λ i } except for λ i , i = 1, 2, . . . , n.
Finally, the FCD of v can be shown to be: given the specification for p(v) in (8). Equation (13) is not a recognizable density such that a Metropolis-Hastings update is required. We utilized a random walk implementation [10] of Metropolis-Hastings sampling; specifically, a normal density with expectation equal to the parameter value from the previous MCMC cycle was used as the proposal density for drawing from the FCD of κ = log(v), using equation (13) and the necessary Jacobian for this transformation. The Metropolis-Hastings acceptance ratio was tuned to intermediate rates (40-50%) during the MCMC burn-in period to optimize MCMC mixing [10], adapting the tuning strategy of Müller [32]. Since the variance of a t-density is not defined for v ≤ 2, we truncate the sample from (13) such that v > 2, or equivalently κ > log(2), consistent with work by previous investigators ( [42,47]).

MODEL COMPARISON
Model choice is an important issue that has not received considerable attention in animal breeding until only very recently (e.g. [20,35]). Likelihood ratio tests have been used to compare differences in fit between various models and their reduced subsets; however, these tests do not facilitate more general model comparisons. The Bayes factor has a strong theoretical justification as a general model choice criterion; however algorithms for Bayes factor computations are either computationally intensive (e.g. [9]) or numerically unstable [33]. Furthermore, as Gelfand and Ghosh [15] indicate, Bayes factors lack clear interpretation in the case of improper priors which are particularly frequent specifications in animal breeding hierarchical models. The Akaike information criterion or Schwarz Bayesian criterion are analytical measures that provide an asymptotic representation of Bayes factors and reflect a compromise between goodness of fit and number of parameters. Since the total number of parameters and latent variables often exceeds the number of observations in animal breeding (e.g. animal model) analysis, the effective number of parameters in hierarchical models is not always so obvious. The MCMC sample average of the posterior log likelihoods, or data sampling log densities, may be used as a means for comparing different models [12]; however, as Speigelhalter et al. [40] indicate, it is not always so obvious how to proceed when these densities are similar but the number of parameters and/or the numbers of hierarchical stages of the candidate models vary. Speigelhalter et al. [40] proposed the deviance information criterion (DIC) for comparing alternative constructions of hierarchical models. The DIC is based on the posterior distribution of the deviance statistic, which is −2 times the sampling distribution of the data as specified in the first stage of a hierarchical model. However, it may not be obvious how to specify the data sampling stage in a hierarchical model. For example, the data sampling stage for the CT model may be specified in one way as: (14) given the specifications of (2a) and (5a) or it may be specified more marginally using (3). We prefer a more marginalized or heavier-tailed first stage specification such as (3) for CT and (1) for CP, potentially leading to a more stable implementation with justification provided by Satagopan et al. [38] but with their context being the stabilization of the Bayes factor estimator of Newton and Raftery [33].
The DIC is computed as the sum of average Bayesian deviance (D) plus the "effective number of parameters"(p D ) with respect to a model, such that smaller DIC values indicate better fit to the data. Let G denote the number of cycles after convergence in an MCMC chain. Furthermore, we represent all unknown parameters in the marginalized first stage specification by ϑ = (β, u, τ, v) with ϑ excluding v = ∞ in the CP model. Then, for the CT model, the average Bayesian deviance can be estimated using (3) bȳ where the superscript [g] denotes the MCMC cycle g, g = 1, 2, . . . , G for the sampled value of the corresponding parameter. Furthermore, p D can be estimated as Here the bar notation (e.g.θ) denotes the corresponding posterior mean vector.
We alternatively considered the conditional predictive ordinate (CPO) as the basis for model choice [14]. Defined for observation i, we write the CPO as: using (1) for the CP model (Model M 1 ) and Here y −i denotes all observations other than y i . The log marginal likelihood (LML) of the data for a certain model, say M k , can then be estimated as: A pseudo Bayes factor (PBF) between two models, say Model M 1 and Model M 2 , can be determined by computing the antilog of their LML difference, that is, (15) Under the assumption of equal prior model probabilities, PBF 1,2 can be interpreted as a surrogate Bayes factor measure [14] and hence the approximate posterior odds of Model 1 relative to Model 2.

Simulation study
A simulation study was used to validate the CT model and the utility of the DIC and the PBF for model choice between CP and CT. Three replicated datasets were generated from each of two different populations as characterized by the distribution of the liability residuals. Population I had a residual density that was standard Student-t distributed with scale parameter σ 2 e = 1 and degrees of freedom v = 4 whereas Population II had a residual density that was standard normal. All datasets were generated based on a simple random effects (sire) model with a null mean. Liability data for 50 progeny from each of 50 unrelated sires was generated by summing independently drawn sire effects from N(0, σ 2 s = 0.10) with independently drawn residuals from N(0, σ 2 e = 1.00) for a total of 2500 records. These underlying liabilities were mapped to ordinal data with four categories based on the threshold parameter values of τ 1 = −0.50, τ 2 = 1.00, and τ 3 = 2.00 for all populations. Ordinal data from each replicated dataset was analyzed using both CP and CT sire models. For the purpose of parameter identifiability, we invoked the restrictions σ 2 e = 1 and τ 1 = −0.50. As a positive control, the underlying liability data for each replicate was analyzed using both normal and t distributed error mixed linear models. For the t-distributed error model, the MCMC procedure adapted was similar to that presented in Stranden and Gianola [42], except that the degrees of freedom parameter (v > 2) was inferred as a continuous (rather than discrete) parameter, using the Metropolis-Hastings update as presented earlier. Graphical inspection of the chains based on preliminary analyses was used to determine a common length of burn-in period. For each replicated data set within each population, a burn-in period of 20 000 cycles was seen to be sufficiently large upon which random draws from each of an additional 100 000 MCMC cycles were subsequently saved. Furthermore, DIC and LML values were computed for each model on each replicated dataset to validate those measures as model choice criteria. For the direct mixed linear model analysis of liability data, DIC and LML measures were based on normal and t-error data sampling densities for their respective models, similar to that implemented for the robust regression example in Speigelhalter et al. [40]. In all cases, flat unbounded priors were invoked on the variance components and on the fixed effects and the vaguely informative prior in (8) was used for v. Furthermore, the effective number of independent samples (ESS) for each parameter was determined using the initial positive sequence estimator of Geyer [16] as adapted by Sorensen et al. [39].

Italian Piemontese calving ease data
First parity calving ease scores recorded on Italian Piemontese cattle from January, 1989 to July, 1998 by ANABORAPI (Associazione nazionale allevatori bovini di razza Piemontese, Strada Trinità 32a, 12061 Carrù, Italy) were used for this study. In order to limit computing demands, only the 66 herds that were represented by at least 100 records over that nine-year period were considered for the demonstration of the proposed methods in this paper, leaving a total of 8847 records. Calving ease was coded into five categories by breeders and subsequently recorded by technicians who visited the breeders monthly. The five ordered categories are: (1) unassisted delivery; (2) assisted easy calving (3) assisted difficult calving (4) caesarean section and (5) foetotomy. Since the incidence of foetotomy was less than 0.5%, the last two ordinal categories were combined, leaving a total of four mutually exclusive categories. The general frequencies of first parity calving ease scores in the data set were 951 (10.75%) for unassisted delivery; 5514 (62.32%) for assisted easy calving; 1316 (14.88%) for assisted difficult calving; and 1066 (12.05%) for caesarean section and foetotomy.
The effects of dam age, sex of the calf, and their interaction were considered by combining eight different age groups (20 to 23, 23 to 25, 25 to 27, 27 to 29, 29 to 31, 31 to 33, 33 to 35, and 35 to 38 months) with the sex of the calf for a total of 16 nominal subclasses. A total of 1212 herd-year-season (HYS) contemporary subclasses were created from combinations of herd, year, and two different seasons (from November to April and from May to October) as in Carnier et al. [7] who also analyzed calving ease data from this same population. The sire pedigree file was further pruned by striking out identifications of sires having no daughters with calving ease records and appearing only once as either a sire or a maternal grandsire of a sire having daughters with records in the data file. Pruning results in no loss of pedigree information on parameter estimation yet is effective in reducing the number of random effects and hence computing demands. The number of sires remaining in the pedigree file after pruning was 1929.
As in Kizilkaya et al. [23], the CP and CT models used for the analysis of calving ease data included the fixed effects of age of dam classifications, sex of calf and their interaction in β, the random effects of independent herd-yearseason effects in h, random sire effects in s and random maternal grandsire effects in m. We assume: , with σ 2 s denoting the sire variance, σ 2 m denoting the maternal grandsire variance, σ sm denoting the sire-maternal grandsire covariance, and σ 2 h denoting the HYS variance. Furthermore, ⊗ denotes the Kronecker (direct) product and A is the numerator additive relationship matrix between sires due to identified male ancestors [19]. Also, h is assumed to be independent of s and m. Flat unbounded priors were placed on all fixed effects and variance components. Based on the poor mixing results from the simulation study and the increased relative computing demands for this larger data set, v was not inferred upon. That is, since the simulation study was much simpler in design compared to the calving ease dataset, any attempt to infer upon v would prove even more difficult. To provide a stark contrast to the CP model, v was then held constant to 4 in the CT model. MCMC inference was based on the execution of three different chains for each model. For each chain in the CP model, a total of 5000 cycles of the burnin period followed by saving samples from each of 100 000 additional cycles was executed based on the experiences of Kizilkaya et al. [23]. Because of initially anticipated slower mixing, the corresponding burn-in period for each chain in the CT model was 10 000 cycles followed by saving each of 200 000 additional cycles. To facilitate diagnosis of sufficient MCMC convergence, the starting values on variance components for each chain within a model were widely discrepant, with one chain starting at the posterior mean of all (co)variance components based on the analysis of Kizilkaya et al. [23], another chain starting at the posterior mean minus 3 posterior standard deviations for each (co)variance component and the final chain starting at the posterior mean plus 3 posterior standard deviations for each (co) variance component.
As with the simulation study, the ESS for each inferred parameter was determined. Furthermore, key genetic parameters, specifically direct heritability (h 2 d ), maternal heritability (h 2 m ) and the direct-maternal genetic correlation (r dm ) were inferred upon in the calving ease data using the functions of G o as presented by Kizilkaya et al. [23] and Luo et al. [25], for example. The only difference in the computation of heritabilities between the CP and the CT model was that the marginal residual variance for the underlying liabilities was not σ 2 e in CT, as it is in CP, but is equal to v v − 2 σ 2 e [42]. Posterior means and the standard deviation of elements of s were also compared between the CP and the CT model. Table I summarizes inferences on v based on the replicated datasets from the two populations, comparing the CP versus CT models for the analysis of ordinal categorical data and comparing the Gaussian linear mixed model versus the t-error linear mixed model for the analysis of the matched latent or underlying normal liabilities, as if they were directly observed. Inference on v Each population replicated for 3 datasets each.

Simulation study
b posterior mean ± standard deviation. c 95% equal-tailed posterior probability interval based on 2.5th and 97.5th percentiles. d the effective number of independent samples using the initial positive sequence estimator of Geyer [16].
was surprisingly sharp and seemingly unbiased for the t-error mixed model analysis of liability data from Population I (v = 4), with 95% equal-tailed posterior probability intervals (PPI) not exceeding 1.5 in width; furthermore, the corresponding ESS were relatively large indicating stable MCMC inference. We define a 95% PPI to be the range of the posterior density falling between the 2.5th and 97.5th percentile. Conversely, inference on v based on the t-error mixed model analysis of liability data from Population II (v = ∞) indicated extremely wide 95% PPI and posterior means exceeding 100, indicating stronger evidence of Gaussian distributed versus t-distributed residuals for data from that population. Furthermore, ESS were generally very small (∼ 20) indicating that inference on v was rather unreliable for data from Population II, at least given the specified MCMC sampling scheme. That is, five times as many MCMC samples would be needed to attain a minimum of 100 as advocated by previous investigators [6,44].
Inference on v in ordinal data under the CT model was also interesting. In Population I, the 95% PPI correctly concentrated on low values for v although the PPI were understandably wider than for the corresponding analyses of liability data under t-error linear mixed models. Also, the ESS on v were considerably smaller (< 100) for the CT model analysis of ordinal data than for the corresponding matched linear model analyses of liability data, such that acceptably accurate inference on v would require substantially more sampling than what we considered in this study. In replicated ordinal data from Population II, the 95% PPI on v were wide and concentrated on high values of v, consistent with what was expected. Furthermore, as with the t-error mixed model analysis of liability data, MCMC mixing on v using the CT model on ordinal data was seen to be particularly problematic in Population II as manifested by the small ESS. Table II summarizes inference on σ 2 s based on the replicated datasets from the two populations, comparing the CP versus CT models for the analysis of ordinal categorical data and comparing the Gaussian linear mixed model versus the t-error linear mixed model for the analysis of the underlying liabilities. In the analyses of liability data from replicated datasets from both populations, the 95% PPI were in good agreement with σ 2 s = 0.10; furthermore, very large ESS indicated very good MCMC mixing. MCMC mixing on σ 2 s was understandably slower in the analysis of ordinal categorical data, particularly in replicated data from Population II using CT, due to the generally high posterior sampling correlation between σ 2 s and v. Because of the problem of MCMC mixing of v in ordinal data, the MCMC chains were rerun with v held constant (v = 4) for more reliable DIC and LML comparisons between models. Both DIC and LML are reported in Table III for each replicated dataset within each population for the CP versus the CT model analyses of ordinal data and for the Gaussian linear mixed model versus the 0.08 ± 0.02 0.05-0.13 42 514 0.12 ± 0.03 0.07-0.19 187 a Population I specified with t residual distribution with 4 degrees of freedom; Population II specified with Gaussian residual distribution. Each population replicated for 3 datasets. b posterior mean ± standard deviation. c 95% equal-tailed posterior probability interval based on 2.5th and 97.5th percentiles. d the effective number of independent samples using the initial positive sequence estimator of Geyer [16]. t-error linear mixed model analyses of liability data. Speigelhalter et al. [40] suggested that a DIC difference exceeding 7 was a substantial indication of an important difference in model fit. Raftery [34] suggested Bayes factors exceeding 12 or, equivalently, LML differences exceeding 2.5, as being important.
Given that, the model choices based on DIC and PBF for the linear mixed model analyses of liability data were always resoundingly in favor of the correct model. However, in the comparison between CP and CT models for ordinal data analysis, the correct (CT) model was decisively chosen in only one of the replicates of Population I (v = 4) and, similarly, the correct (CP) model was decisively chosen in only one of the three replicates of Population II (v = ∞), with all other comparisons being indecisive (i.e. DIC differences < 7 and LML differences < 2.5). These results should not be too surprising since the much lower information content of ordinal data relative to underlying continuous liability data would make it that much more difficult to clearly resolve model choice. It is interesting to note, however, that LML and DIC measures led to identical conclusions with respect to decisive or indecisive model choice. In particular, −2LML and DIC values were virtually identical to each other.

Genetic parameter inference
Sire and maternal grandsire CP and CT models were used for the analysis of calving ease scores in Italian Piemontese cattle. Because of the MCMC mixing problems encountered in inferring upon v, this parameter was held constant to v = 4. Posterior inferences on key genetic parameters are summarized in Table IV and are based on the combined results from each of the three separate chains. The posterior mean, median and modal estimates (not shown) of the two heritabilities, and the genetic correlation using the MCMC algorithms were similar to each other across models, implying that the posterior densities were symmetric and unimodal. This property was further manifested by the fact that the 95% PPI are closely matched by the posterior mean ±2 standard deviations. In this study, the total ESS for dispersion parameters across the three chains ranged from 1420 to 9305, indicating sufficient MCMC mixing under both models. Table IV shows that the ESS from the CT model were found to be almost double those of the CP model, attributable to the twice as large post-burn-in period for the CT model. Considering v as known improved mixing of, particularly, genetic parameters, in the CT model relative to joint inference with v (results not shown). Although the n × 1 auxiliary variable vector λ is additionally included in the CT model, this augmentation does not appear to adversely impact ESS and hence mixing of key genetic parameters relative to the CP model.
In this study, the CT model produced posterior means of genetic variance components that were nearly twice as large as those estimated using the CP model. Furthermore, the marginal residual variance is v v − 2 σ 2 e in the CT model such that seemingly twice as much residual variance is inferred under a CT model (with v = 4) than under a CP model (with v = ∞). Although these d the total effective number of independent samples across the three MCMC chains using the initial positive sequence estimator of Geyer [16]; interval values [ ] refer to a range in ESS across the three chains. overlapped substantially between the two models. Furthermore, the posterior density of r dm was very similar between CT and CP, with most of the density being between −0.2 and −0.8. In order to compare the CP and CT models for fit to the calving ease data, LML and DIC values, broken down into its componentsD and P D , are reported in Table V. Since it is difficult to quantify the degree of Monte Carlo error on DIC [50], we report DIC and LML values separately for each of the three chains under each model. It can be seen that there are relatively inconsequential differences in the measures from one chain to the next within each model relative to chains between models, thereby indicating minutely small Monte Carlo errors on the DIC difference between the two models. Both model choice criteria were overwhelmingly in favor of the CT model with v = 4. As anticipated, the model complexity, as measured by P D , is higher for the CT model; however, the complexity penalty was strongly counteracted by a smaller mean devianceD, thereby resulting in a smaller DIC favoring choice of the CT model.

Inferences on sire effects
Posterior means of elements of s were determined to be corresponding point estimates of progeny differences (EPD) under both CT and CP models.
The relationships between these estimates are shown to be strongly linear with no hint of substantial reranking as indicated by a Pearson correlation of 0.99. The CT model had greater spread in EPD's compared to the CP model, as further manifested by a least-squares estimated slope of 1.38 of CT on CP EPD's. This is not too surprising since a larger additive genetic (sire) variance was inferred in the CT model such that posterior means of elements of s should be more dispersed in the CT model relative to the CP model. However, as discussed later, this is not practically important since the variance of L is defined only proportionately to the marginal residual variability which is also larger in the CT model.
Posterior standard deviations of elements of s are analogous to standard errors of prediction (SEP) in mixed effects model analysis and can be used to derive approximate reliabilities of EPD's [49]. That is, the SEP's were simply determined as the standard deviation of the MCMC samples of elements of s. In this case, the correlation between the estimated standard errors was near unity. The estimated slope of CT on CP SEP's was nearly 1.41 indicating that the posterior standard errors were on average 41% larger under the CT model than under the CP model; nevertheless, for the purpose of reporting reliabilities, these standard errors need to be considered relative to σ 2 s which was concentrated on higher values under the CT model.

DISCUSSION
Given the recent momentum in using heavy-tailed residual specifications for the analysis of production data in animal breeding [36,41,42,47] a hierarchical threshold (CT) mixed model based on a cumulative t-link specification was developed, validated by simulation and applied to a small calving ease dataset from Italian Piemontese cattle. The simulation study indicated that inference on v is possible in a CT model; however, it appears that either a more suitable MCMC strategy is needed or many more samples are required compared to that considered in our study to ensure a more reliable inference on v. Until this issue is satisfactorily resolved, we advocate fixing v at some arbitrarily low value in a CT model analysis. We chose v = 4 since this value minimally assures defined first, second, and third moments while providing a liability variable distribution that is maximally heavy-tailed. One can then use model choice criteria such as DIC to assess whether or not the CT model is a better data fit than the CP model. We further note that for the case where v is fixed, that MCMC mixing was not negatively affected using the CT model, even though our data augmentation (of λ) implementation might be of concern to those who might prefer Metropolis-Hastings sampling on all parameters [47] instead of introducing augmented variables. More recently, it has been demonstrated that data augmentation can be strategically used to enhance MCMC mixing; the strategies discussed by van Dyk and Meng [45] may facilitate more favorable mixing on v and hence deserve further consideration in CT model applications to animal breeding. One reviewer pointed out that several Metropolis-Hastings updates on v per each MCMC cycle, as utilized by Bink et al. [6], and Uirmari et al. [44] may facilitate improvements in mixing. Of particular note (results not reported) was that due to the implementation of the algorithm of Cowles [11], MCMC mixing of τ was not seen to have been the most limiting in terms of ESS in contrast to previous animal breeding implementations [39,46].
Our point estimates for heritabilities are substantially higher than corresponding threshold model estimates for calving ease reported by Manfredi et al. [27,28], McGuirk et al. [30,31], Varona et al. [46], Luo et al. [25], and Bennett and Gregory [3]. Nevertheless, our inference on a strongly negative direct-maternal genetic correlation is in agreement with previous work on calving ease using threshold models [3,25,46] and linear mixed models using data from the same source [7]. Hence, it appears from our results, in agreement with other studies, that selection of sires for calving ease of their progeny as calves should result in antagonistic effects in the ability of their daughters to calve easily as dams in successive generations. What was the most surprising is that the 95% PPI for h 2 d are greater and do not overlap with corresponding PPI in Kizilkaya et al. [23] who used a larger data set on first parity records from herds with greater than 50 rather than 100 records over the nine year period from the same data source. This result may be indicative of heterogeneity in genetic and residual variance due to the size of the herd or other confounding factors (e.g. region); this is an area for further research that our group has started with respect to residual variability, extending the work of San Cristobal-Gaudy et al. [37].
Two Bayesian model choice criteria, DIC and LML, were used to choose between the CP and CT models. In a simulation study, it was demonstrated that both DIC and LML were able to decisively choose the correct model in some cases whereas, in the remaining cases, these measures were too similar between the two models to allow a definitive choice. In the analysis of calving ease scores in Italian Piemontese cattle, the CT model was overwhelmingly chosen as the best fitting model by both model choice criteria. Nevertheless, in the examination of EPD's there were no real tangible differences between the two models in terms of sire genetic rankings.
Our study involved sire models, where calf records are connected directly to sires with genetic relationships identified only through known paternal ancestry. Sire models are different from animal models where each record is connected directly to a calf identification with all known genetic (maternal and paternal) relationships explicitly modeled. For animals (e.g. dams), with effectively less data information, we would anticipate greater differences in predicted genetic merit between CT and CP animal models than those given in our study.
However, a sire model has been seen to be more stable than an animal model in CP implementations [29], particularly when paternal half-sib relationships are dominant as in our data. Furthermore, the fact that a sire and maternal grandsire model only accounts directly for a portion of the additive genetic variance and of the maternal genetic variance implies that a t-error assumption is essentially placed on a composite source of error, i.e. the sum of the residual variance and remaining genetic variation, attributable to unknown dams and to Mendelian sampling. From the perspective of using heavy tailed densities to mute residual outliers, this is significant since deviant dam and/or Mendelian sampling effects may be muted as well in a CT sire model specification.
Presently, it does not appear to be feasible to apply MCMC methods to the very large datasets used for routine genetic evaluation of livestock by breed associations and national recording organizations. Kizilkaya et al. [23] recently demonstrated, however, very little differences in predicted genetic merit and standard error of prediction in a CP sire model between inference provided by MCMC and by approximate empirical Bayes procedures currently utilized by the industry (e.g. [4]). Empirical Bayes procedures are based on using the joint posterior modes of the sire effects as the EPD's, conditional on estimated variance components as if they were known with certainty. Based on results from Gianola and Foulley [17], the CT model can be readily implemented using empirical Bayes methods since the probability density function and the cumulative density function of a standard Student t-distribution could be substituted for the corresponding Gaussian functions needed to derive the necessary scoring equations used to compute the joint posterior modes. Furthermore, the degrees of freedom parameter, v, could be jointly estimated with the variance components using the Laplace method and tested for statistical significance using marginal likelihood ratio tests [43]. Empirical Bayes implementation of the CT sire model and the comparison of results with MCMC inference under the CT sire model deserve further consideration. Unfortunately, however, these comparisons may not necessarily apply to CP or CT animal models since joint modal estimates of EPD's in the CP animal model can be badly biased [29].
In this study, we have considered only two cumulative link models for the analysis of calving ease; conceptually, there are many others including those proposed by Chen and Dey [8]. In fact, Albert and Chib [2] demonstrate that the cumulative logistic link model is roughly equivalent to the CT model with v = 8. Other models can be contrived by considering alternative heavytailed distributions for the underlying variables, such as those considered by Rosa [36]. Some of the resulting models may be shown to have better fit to calving ease data than demonstrated with the CT model in our paper. Specifications based on skewness [47] may also have merit.
The substantial residual and genetic correlations between birth weight and calving difficulty imply that genetic evaluations of calving ease would substantially benefit from a bivariate threshold/linear multiple trait analysis with birth weight [26,46]. Further work on jointly providing modeling flexibility with t-distributed specifications on both traits is needed given our results and those already presented by Stranden and Gianola [41,42] and Rosa [36].