HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Cardoso, F. F. Articles by Tempelman, R. J. Search for Related Content PubMed PubMed Citation Articles by Cardoso, F. F. Articles by Tempelman, R. J.

Multiple-breed genetic inference using heavy-tailed structural models for heterogeneous residual variances¹

F. F. Cardoso^*, G. J. M. Rosa and R. J. Tempelman^,2

^* EMBRAPA Pecuária Sul (Brazilian Agricultural Research Corporation Southern Cattle and Sheep Center), Bagé, RS 96401-970, Brazil; and and Department of Animal Science, Michigan State University, East Lansing 48824

	Abstract

Top Abstract Introduction Materials and Methods Results and Discussion Implications Appendix Literature Cited

 
Multiple-breed genetic models recently have been demonstrated to account for the heterogenous genetic variances that exist between different beef cattle breed groups. We extend these models to allow for residual heteroskedasticity (heterogeneous residual variances), specified as a function of fixed effects (e.g., sex, breed proportion, breed group heterozygosity) and random effects such as contemporary groups (CG). We additionally specify the residual distributions to be either Gaussian or based on heavier-tailed alternatives such as the Student’s t or Slash densities. For each of these three residual densities using either homoskedastic (homogeneous variance) or heteroskedastic error specifications, we analyzed 22,717 postweaning gain records from a Nelore-Hereford population based on a Markov chain Monte Carlo animal model implementation. The heteroskedastic Student’s t error model (with estimated df = 7.33 ± 0.48) was clearly the best-fitting model based on a pseudo-Bayes factor criterion. Breed group heterozygosity and, to a lesser extent, calf sex seemed to be marginally important sources of residual heteroskedasticity. Specifically, the residual variance in F₁ animals was estimated to be 0.70 ± 0.16 times that for purebreds, whereas the male residual variance was estimated to be 1.13 ± 0.09 times that for females. The CG effects were important random sources of residual heteroskedasticity (i.e., the coefficient of variation of CG-specific residual variances was estimated to be 0.72 ± 0.06). Purebred Nelores were estimated to have a larger genetic variance (124.84 ± 21.75 kg²) compared with Herefords (40.89 ± 6.70 kg²) under the heteroskedastic Student’s t error model; however, the converse was observed from results based on a homoskedastic Student’s t error model (46.24 ± 10.90 kg² and 60.11 ± 8.54 kg², respectively). These results indicate that allowing for robustness to outliers and accounting for heteroskedasticity of residual variances has potentially important implications for variance component and genetic parameter estimates from data on multiple-breed populations.

Key Words: Bayesian Inference • Multiple Breeds • Outlier Robustness • Residual Heteroskedasticity • Structural Variance Models

	Introduction

Top Abstract Introduction Materials and Methods Results and Discussion Implications Appendix Literature Cited

 
Beef cattle performance is generally recorded across highly diverse production systems and environments, with data quality often compromised by frequent occurrences of recording error, preferential treatment, and the potential effect of injury or disease. Hierarchical Bayes models have provided a general framework to address problems arising from data characterized by different types of heterogeneity across various environments. For example, such models have been used in animal breeding to address the issue of heterogeneous residual variances, defined as residual heteroskedasticity, across various conditions or environments (Foulley et al., 1992; Gianola et al., 1992; San Cristobal et al., 1993), as well as robustness to outliers, perhaps due to preferential treatment (Stranden and Gianola, 1998, 1999). These two issues generally have been investigated separately; however, there is no conceptual difficulty in considering them jointly (Kizilkaya and Tempelman, 2003, 2005; Sorensen and Waagepeterson, 2003). As we shall demonstrate, these extensions can be further superimposed on a hierarchical model that accounts for heterogeneous genetic variances across breed groups in multiple-breed genetic evaluations (Birchmeier et al., 2002; Cardoso and Tempelman, 2004).

The objectives of this study were to 1) present alternative multiple-breed hierarchical Bayes models based on Gaussian or heavy-tailed specifications of residual heteroskedasticity, and 2) apply the proposed models to a dataset on postweaning gains of Nelore-Hereford cattle raised in diverse environments, so as to characterize the magnitude of potentially important sources of residual heteroskedasticity and the associated influence on breed group specific heritabilities.

	Materials and Methods

Top Abstract Introduction Materials and Methods Results and Discussion Implications Appendix Literature Cited

 
Construction of Hierarchical Bayes Model

The first stage of the model specifies the conditional sampling density of the n x 1 data vector y. We write the linear model for subject or element j of y as

[1]

where µ is a fixed reference or overall mean parameter, ß is a p x 1 vector of fixed effects (e.g., gender, age of dam, heterotic effects, etc.), u is a q x 1 vector of random effects, and x'_j, and z'_j are known row incidence vectors connecting y_j to ß and u, respectively. Finally, across all j = 1, 2, ..., n, where ²_ej represents the residual variance that is specific for y_j.

We write ²_ej to be a multiplicative function of fixed and random effects:

[2]

Here, the term 2e functions as a reference parameter in Eq. [2], much like µ does in Eq. [1], albeit on a multiplicative rather than an additive scale. Secondly, = [1 2 ... o]' specifies regression parameters that potentially influence residual heteroskedasticity using information in covariates pj = [p1j p2j ... p oj ]' specified for subject j. Examples of elements in include not only overall or global regression parameters, but also classification specific regression parameters using, for example, geographically specific breed proportion or breed heterozygosity coefficients as separate covariates in pj. Now, [1(j) 2(j) ... r(j) ] specifies the "fixed" effects for levels of each of r different sets of classification factors (e.g., sex of calf) associated with subject j. For example, if the third level (i = 3) of a second factor is associated with subject j, then 2(j) specifies the effect of the third element of 2 = {2i}t2i = 1 specifying the t2 heteroskedastic effects for the second factor as further presented later. In other words, if = [1' 2' ... r']' is used to specify all fixed classification effects influencing residual heteroskedasticity, then being a row incidence vector connecting log() to based on the logarithm of Eq. [2]. Similarly, specifies the "random" effects for levels of each of s different sets of classification factors (e.g., contemporary groups) associated with subject j. Finally, wj represents an optional random "noise" term as useful for specifying heavy-tailed alternatives to a Gaussian residual distribution as illustrated later.

Equation [2] is similar to expressions considered by Foulley et al. (1992) and much subsequent work, except that previous presentations have been based extensively on the logarithm of Eq. [2], thereby rendering an expression in as a linear function of fixed and random effects. On that logarithmic scale, log(w^–1_j) then represents an additional source of random variability that has not been previously considered in heteroskedastic residual modeling applications. Note that random regression parameters additionally are possible in Eq. [2] for specifying, for example, subject-specific changes in residual variability over repeated measurements.

In the second stage, we specify our prior assumptions on all unknowns defined in Eq. [1] and [2]. For classical random effects, we adopt the typical structural prior:

[3]

In multiple-breed models, G() includes, at the very least, a function of more than one genetic (co)variance parameter for additive genetic effects as defined by Lo et al. (1993). That is, includes ²_Ab as the additive variance of breed b, b = 1, 2, ..., B, and ²_Sbb' as the variance due to the segregation between breed b and b', where b' b, and b' = 1, 2, ..., B. Note that ²_Sbb' can be interpreted as the additional genetic variability observed in the F₂ generation relative to the F₁ generation. Additional hierarchical modeling details on multiple-breed random effects are provided in Cardoso and Tempelman (2004). In addition, u may include nongenetic effects such as contemporary groups (CG) with variance ²_CG, which are then also included in . Finally, subjective multivariate normal or bounded uniform priors can be adopted on ß.

As noted above and by others (Sorensen and Gianola, 2002), the Bayesian distinction between classical fixed effects ß and random effects u is based on the specification of subjective vs. structural priors. That is, a structural prior density implies a pure distributional characterization of the elements of u such that elements of structurally define the hierarchical model and hence are inferred. Conversely, parameters defining subjective priors on ß are specified to be "known," either to ensure a proper joint posterior density or to incorporate knowledge from personal beliefs or previous studies. We have previously drawn very similar distinctions between these two classes of priors for mixed effects modeling of residual variances (Kizilkaya and Tempelman, 2005). Here, we adopt subjectively specified inverted gamma distributions for the fixed effects influencing residual heteroskedasticity as follows:

[4]

[5]

and

[6]

Here, _ki represents the fixed effect of level i of the t_k x 1 vector of fixed effects influencing residual heteroskedasticity as defined within each factor k = 1, 2, ..., r. Furthermore, the hyperparameters ^*_e, ß^*_e, ^*_m, ß^*_m, ^*_ki, and ß^*_ki are presumed to be known. To facilitate parameter identifiability in the absence of suitably informative prior specifications, such as Eq. [4], [5], and [6], or based on the use of bounded uniform priors, the effect _ktk, corresponding to the last level t_k within each factor k, is constrained (i.e., _ktk= 1, k = 1, 2, ..., r). This type of parameterization is not conceptually different from the corner parameterization (Clayton, 1996) typically invoked for elements of ß in regular linear models software and previously adapted by Kizilkaya and Tempelman (2005) in a simpler heteroskedastic error specification.

As with classical random effects u, structural rather than subjective inverted gamma priors are specified on random effects influencing residual heteroskedasticity (Kizilkaya and Tempelman, 2005):

[7]

Here, _li represents the random effect of level i within the t_l x 1 vector of random effects v_l of factor l, l = 1, 2, ..., s. It can be seen that E(_li) = 1 with (and hence only defined for _l > 2) for i = 1, 2, ..., t_l, such that the standard deviation and the CV are synonymous for random effects within each factor l = 1, 2, ..., s. Note that as _l , the influence of random factor l on residual heteroskedasticity diminishes. With the structural specification in Eq. [7], there is a borrowing of information across the t_l effects of v_l, very much analogous to what is done for classical random effects modeling of location parameters u as in Eq. [3] (Robinson, 1991).

Now, there are several possible distributional assumptions on the weights {w_j}, each yielding a different heavier-tailed specification on the distribution of the residuals e_j relative to the use of the Gaussian distribution (see also Lange and Sinsheimer, 1993; Rosa et al., 2003). We specifically consider two heavy-tailed alternatives to the Gaussian or normal distribution: the Student’s t and the Slash distributions. The Student’s t residual distribution has been recently noted as promising for mitigating the effects of preferential treatment and/or deviant observations (Stranden and Gianola, 1998; Pinheiro et al., 2001), whereas the Slash distribution was demonstrated to be a better fit relative to the Student’s t for modeling residual effects in other applications (Rosa et al., 2003).

The Gaussian distribution requires no distributional specification of w_j in Eq. [2], with w_j = 1 for all j = 1, 2, ...., n. A distributional specification on w_j represents one way of modeling the lack of fit of the marginal density of e_j to a Gaussian distribution. Let us momentarily consider the simpler homoskedastic error specification (i.e., ²_ej = ²_e for all j), such that conditionally using Eq. [2], . Given a distributional specification p(w_j|_e), conditional on _e, the marginal density p(e_j|²_e, _e) of e_j is determined by p(e_j|²_e, _e) = e_j|²_e,w_j)p(w_j|_e)dw_j. For example, consider a Gamma specification on w_j:

[8]

It then can be readily shown that the marginal density p(e_j|²_e, _e) is a Student’s t density with scale parameter ²_e and df _e, such that the marginal residual variance . Alternatively w_j may be specified as having the following density:

[9]

thereby leading to a Slash distribution specification for p(e_j|²_e, _e), with scale parameter ²_e and df _e such that the marginal variance of e_j is . It is possible to extend heavy-tailed specifications to heteroskedastic specifications using ²_ej from Eq. [2] rather than in p(e_j|²_ej) in tandem then with Eq. [8] or [9]. The resulting marginal residual densities would then be, respectively, heteroskedastic marginal Student’s t or Slash densities with subject-specific scale parameters as specified by the numerator of Eq. [2].

We recommend the same inverted gamma prior densities on breed-specific genetic variances and segregation variance and remaining variance components in as considered by Cardoso and Tempelman (2004). Similarly, arbitrarily vague or noninformative subjective prior densities could be specified on _l, l = 1, 2, ..., s, and on _e as well.

Markov Chain Monte Carlo Inference

Hierarchical models with several stages, such as those considered in this paper, are particularly amenable to Bayesian inference using Markov chain Monte Carlo (MCMC) methods. The full conditional densities of all unknown parameters/quantities or blocks thereof necessary for inference based on MCMC are derived in the Appendix to this paper.

Application to Analysis of Nelore-Hereford Data on Postweaning Gain

Data Source. Data analyzed in this study consisted of 22,717 postweaning gain (PWG) records of a beef cattle population comprising Herefords, Nelores, and their crosses under genetic evaluation in Brazil. The records were collected between 1974 and 2000 under a large-scale breeding program called "Conexão Delta G." Extensive details on the characteristics and structure of the data are provided in Cardoso and Tempelman (2004).

Candidate Models. The linear mixed model in Eq. [1] was based on very similar specifications as previously outlined in Cardoso and Tempelman (2004). The elements of ß in Eq. [1] included the classification effects of region and sex, the linear effect of the length of the postweaning test period, the linear and quadratic effects of age of dam and their interactions with the additive effects of dam breed, the Nelore breed effect (ß_A1), the Nelore-Hereford dominance or heterotic effect (ß_D12), and the Nelore-Hereford additive x additive epistatic effect (ß_AA12). The corresponding subject-specific covariate information for ß_A1, ß_D12, and ß_AA12are the proportion of Nelore alleles f^j₁, the heterozygosity coefficient f^j₁₂, and the epistatic coefficient f^j₁(1 – f^j₁) for subject j, respectively; further clarification of these terms is provided in Cardoso and Tempelman (2004). Additionally, interactions between the Nelore breed effect with the linear effect of postweaning test period length and the classification of sex were modeled in ß. Other interactions were not considered for reasons provided in Cardoso and Tempelman (2004).

The random effects u in Eq. [1] included the effects of 940 normally, independently, and identically distributed (NIID) CG effects with variance component ²_CG in addition to the additive genetic effects of 40,082 animals, including ancestors with no data of their own. As clearly outlined previously by Cardoso and Tempelman (2004), the (co)variances of the genetic effects depend on breed-specific genetic variances ²_A1 and ²_A2 for the Nelores and Herefords, respectively, as well as for the segregation variance ²_S12 between the two breeds.

Six different specifications of the distribution of the residuals were considered as based on a two-way factorial arrangement. The first factor was defined by the presence or absence of residual heteroskedasticity in a marginal sense (i.e., present based on the specification of Eq. [2] vs. absent based on the homoskedastic specification of ). The second factor was based on three alternative specifications for the marginal density of the residuals: Gaussian, Student’s t, or Slash. We labeled the six models M₁ through M₆.

M₁: Homoskedastic Gaussian Error. This model was previously used by Cardoso and Tempelman (2004), whereby residuals were assumed to be normally, independently, and identically distributed with common residual variance ²_e.

M₂: Homoskedastic Student’s t Error. In this model, the residuals are specified to have independently and identically distributed Student’s t densities with common scale parameter ²_e and common degrees of freedom _e. This specification is conditionally equivalent to , with w_j for j = 1, 2, ..., n, having independent

Gamma densities as indicated previously.

M₃: Homoskedastic Slash Error. Here, the residuals are specified to have independently and identically distributed Slash densities with common scale parameter ²_e and common degrees of freedom _e. This specification is conditionally equivalent to , where w_j for j = 1, 2, ..., n, are specified as independently distributed with density as in Eq. [9].

M₄: Heteroskedastic Gaussian Error. As with M₁, the residuals are specified to be independently and normally distributed, but with heterogeneous residual variances as determined by Eq. [2] and constant w_j = 1 for j = 1, 2, ..., n.

M₅: Heteroskedastic Student’s t Error. As with M₂, the residuals are specified to be independently distributed as independent Student’s t with common degrees of freedom _e but with heterogeneous scale parameter as specified by the numerator of Eq. [2]. Recall that this specification is conditionally equivalent to the distributional assumption for e_j in Eq. [2] followed by the distribution as specified in Eq. [8].

M₆: Heteroskedastic Slash Error. As with M₃, the residuals are specified to have independently distributed Slash densities with common degrees of freedom _e, but with heterogeneous scale parameters as specified by the numerator of Eq. [2]. Recall that this specification is conditionally equivalent to distributional assumption for each e_j in Eq. [2] in tandem with each w_j having the density specified in Eq. [9].

For each of the three heteroskedastic error models M₄ through M₆, the specific model chosen for residual heteroskedasticity was based on a single fixed classification effect for sex, representing the male effect by ₁, thereby requiring the identifiability constraint ₂ = 1 for females. The model also included fixed regression parameters ₁ and ₂, using information on Nelore breed proportion (p_1j = f^j₁) and heterozygosity coefficient (p_2j = f^j₁₂), respectively, for subjects j = 1, 2, ..., n. Finally, the effects v = [_{1 2} ... ₉₄₀]' of 940 CG within a single random effects factor also were used to model residual heteroskedasticity with a CV of as based on Eq. [7]. For example, the residual variance for a particular F₁ male animal j (i.e., f^j₁ = 0.50, f^j₁₂ = 1.00) in CG 2 based on the Gaussian heteroskedastic error model would be determined as ²_{e 1}(₁)^0.50(₂)^1.00 ₂, whereas a female within the same breed group and CG would have a specified residual variance of ²_e (₁)^0.50(₂)^1.00 ₂.

Prior Specifications. For all six models, we used the same linear mixed model in Eq. [1] and the same subjective or structural prior specifications on all parameters with the natural exceptions on w^j, j = 1, 2, ..., n and _e, where applicable. Bounded uniform priors were placed on all identifiable elements of ß, , and . Conjugate but highly dispersed inverted gamma prior densities were specified for variance components in , specifically in a manner similar to that for ²_e in Eq. [5]. Specifically, the corresponding hyperparameters for ²_CG were *_CG = 2.5 and ß^*_CG = 2125. Similarly, for both ²_A1 and ²_A2, the hyperparameters were ^*_A = 2.5 and ß^*_A = 200, whereas for ²_S12, the hyperparameters were ^*_S = 2.5 and ß^*_S12 = 25. These specifications also were used previously for the same data in Cardoso and Tempelman (2004). For the reference parameter ²_e , we specified ^*_e = 2.5 in Eq. [5], with ß^*_e chosen such that the prior mean on the marginal residual variance was the same for each of the three residual density types (i.e., Gaussian, Student’s t or Slash). Now ß*_e was chosen to be based on a REML estimate of 350 for ²_e using the PWG data. Accordingly, ß^*_e = 2.5 x 350 for both Gaussian error models (M₁ and M₄), for both Student’s t error models (M₂ and M₅), and for both Slash error models (M₃ and M₆), where E(_e) denotes the prior mean of _e within the respective model. These specifications ensure the same prior mode or mean on the marginal residual variance for each model.

Finally, we specified a gamma prior on but based on hyperparameters ^*= 0.03 and ß^*= 0.01, leading to a prior mean of E() = 3 but with very large prior variance of var() = 300. We similarly specified gamma prior densities for _e using ^*_e= 0.04 and ß^*_e= 0.01 for Student’s t error models (M₂ and M₅), and ^*_e= 0.015 and ß^*_e= 0.01 for the Slash error models (M₃ and M₆), such that the respective prior means were E(_e) = 4 and E(_e) = 1.5; these two prior means corresponded to a similar degree of kurtosis for these two distributions as based on a simulation study.

Definition of Genetic Parameters. The additive genetic variance 2Ag of any individual belonging to breed group g was obtained by 2Ag = fg12A1 + fg22A2 + 2(fs1fs2 + fd1fd2)2S12 (Lo et al., 1993), where fib indicates the allelic proportion of breed b in breed group i. Here b = 1 and b = 2 specify Nelore and Hereford, respectively, and superscript i specifies either the breed group of the individual itself (i = g), the breed group of that individual’s sire (i = s), or the breed group of the individual’s dam (i = d). Within each of the three homoskedastic error models (M1 through M3), the marginal residual variance (2E) was the same for each breed group, being , respectively. For Model M4, the marginal residual variance specific for breed group g (2Eg) was determined using the expression 2Eg = 2e11 fg11 fg122, noting that E(l) = 1 across all CG. Similarly, the marginal residual variance was specified to be . Here, 1 = 0.3197 represents the proportion of calves with records that were male such that we report inference on residual variance weighted according to the calf sex proportions in the data structure. Alternatively, setting 1 = 1 would specify 2 Eg with respect to male calves and = 0 specifies 2 Eg relative to female calves within breed group g. The marginal phenotypic variance for subjects in breed group g was then specified as 2 Pg = 2 Eg + 2 Ag, such that the additive heritability was h2Ag = 2 Ag/ 2 Pg for the same breed group. Note, as in our previous work (Cardoso and Tempelman, 2004), we deliberately omit 2CG in the determination of 2Pg to facilitate comparisons of our heritability estimates with intraherd estimates reported in the animal breeding literature (e.g., Koots et al., 1994), in which CG are more commonly modeled as fixed effects.

Definition of Overall Marginal Residual Variance: For the three homoskedastic error models (M₁ through M₃), the overall marginal residual variance ²_E is naturally defined to be identical for all breed groups using , respectively. Similarly, an overall marginal residual variance ²_E as weighted by the existing calf sex proportion (₁ = 0.3197), proportion of calves that were Nelore (₁ = 0.2281), and average heterozygosity coefficient (₁₂ = 0.4347) was determined for each heteroskedastic error model. Hence, for models M₄ through M₆, ²_E was defined by , respectively.

MCMC Implementation. The length of the MCMC chain for PWG was 200,000 cycles after 15,000 cycles of burn-in for all six models. Samples were saved every 10 cycles after burn-in, such that 20,000 samples were saved for posterior inference. Means, modes, and SD of the parameters were obtained from their respective marginal posterior densities. Furthermore, 95% credible sets or posterior probability intervals (95% PPI), as based on the 2.5th and 97.5th percentiles of each posterior density, were determined as interval estimates. For each dispersion parameter, the initial monotone sequence approach (Geyer, 1992) was used to calculate effective sample sizes (ESS) based on Sorensen et al. (1995), which estimates the number of independent samples with information content equivalent to that contained within the 20,000 postburn-in dependent samples.

Model Choice Criterion. We considered the pseudo-Bayes Factor (PBF) (Gelfand, 1996) as a means for model choice. The PBF involves the evaluation of the first stage density in Eq. [1] for each MCMC cycle; let p(y_j |y_(–j), M_f) be the conditional predictive ordinate (CPO) for observation y_j under Model M_f, f = 1, 2, ..., 6. The CPO is intended to be a cross-validation density, suggesting what values of y_j are likely when Model M_f is fitted to all other observations y_(–j) except y_j. Writing = [ß' u' ²_e ' v' w']' where an MCMC approximation for the CPO is estimated by a harmonic mean (Gelfand, 1996), with ⁽ⁱ⁾ denoting the postburn-in MCMC sample for i = 1, 2, ..., G = 20,000.

An approximate log marginal likelihood (LML) over all observations can be obtained by the following: n

Finally, for comparing, say, models M₁ and M₂, the corresponding PBF is determined to be PBF_1:2 = exp(LML₁ – LML₂).

	Results and Discussion

Top Abstract Introduction Materials and Methods Results and Discussion Implications Appendix Literature Cited

 
Model Choice

The log marginal likelihood functions (LML₁ through LML₆) for models M₁ through M₆ were computed to be –100,563, –99,623, –99,711, –99,455, –99,052, and –99,098, respectively. Based on the PBF as derived from these LML, the heteroskedastic t-error model provided the best fit to the PWG data among the six different residual specifications considered. In fact, this model had a PBF of 9.328 x 10¹⁹ compared with the next best-fitting model, the heteroskedastic slash error model. Furthermore, it can be definitively concluded that the conventional homoskedastic Gaussian error model was the poorest choice for this data set, given that its PBF relative to comparisons with the other five models always approached zero. It can be further noted from the LML values that heteroskedasticity was a more important specification relative to heavy-tailedness for the distribution of the residuals, particularly as the heteroskedastic Gaussian error model fitted the data much better than either non-Gaussian homoskedastic error models. Nevertheless, the Student’s t distribution always was a better fit to the data compared with the Slash distribution, whether those comparisons were made within homoskedastic or heteroskedastic error specifications.

Robustness to Outliers

Based on the homoskedastic Gaussian error model, there were 45 standardized residuals lying outside the range of ± 4.0 SD as estimated by the square root of the posterior mean of ²_e. Moreover, the estimated kurtosis of this standardized residual distribution was 2.72, indicating that the distribution was somewhat leptokurtic. These features potentially explain our results of a substantially better fit of non-Gaussian error models to the data. The estimated skewness of the estimated residuals was 0.34, being of moderate magnitude; however, this skewness may be an artifact of heteroskedasticity not considered in a homoskedastic error model. Bayesian hierarchical modeling procedures that allow for skewness in residuals have been developed (Fernandez and Steel, 1998) and implemented in animal breeding (Von Rohr and Hoeschele, 2002); our model hierarchy could be extended accordingly.

The posterior distributions of the robustness parameters (_e) were reasonably symmetric for all non-Gaussian error models (data not shown). The posterior means ± SD of _e were 7.33 ± 0.48 and 2.20 ± 0.09 for the heteroskedastic Student’s t and Slash error models, respectively, whereas corresponding estimates for the homoskedastic Student’s t and Slash error models were 4.79 ± 0.21 and 1.66 ± 0.05, respectively. Note that higher values of _e were associated with heteroskedastic error relative to homoskedastic error specifications. This association was somewhat expected as residuals that otherwise may seem to be substantial outliers under homoskedastic error specifications, and thereby contributing to lower estimates of _e, may not be as such in subclasses with large variances as determined using heteroskedastic error models.

Assessment of Heteroskedasticity Sources

Fixed Effects and Regression Parameters. Because the Slash error specification was uniformly poorer than the Student’s t error specification for fit to the PWG data using both homoskedastic and heteroskedastic error specifications, we will not discuss results involving the Slash error distribution further. Posterior inferences on fixed effects for residual heterogeneity obtained by the two remaining heteroskedastic error models are presented in Table 1.

There was no obvious evidence that residual variability depended on breed proportion based on the wide 95% PPI for ₁ including unity in Table 2. The effects of heterozygosity (₂) on residual heteroskedasticity, however, was more convincing. That is, it was determined that Pr(₂ > 1|y) = 0.0449 for the heteroskedastic Student’s t error model, indicating that as heterozygosity (i.e., heterosis) increases, residual variability decreases. This result is consistent with the theory that heterozygosity acts as homeostatic buffer to minimize environmental variation (Lerner, 1954).

It is instructive to note that the posterior mean ± posterior SD of (0.84 ± 0.07) under the heteroskedastic Gaussian error model was slightly larger than that (0.72 ± 0.06) estimated for the heteroskedastic Student’s t error model. This result was somewhat expected as the Student’s t error specification should attenuate the effects of residual outliers, which would otherwise inflate certain CG-specific variance estimates and therefore estimates of under a Gaussian error specification.

We did not explicitly consider region as a fixed-effects factor for modeling residual heteroskedasticity. Because CG are nested within region, regional differences in residual variability can be investigated by examining elements of the posterior mean of v, or estimated scaling factors, by region. The means ± the SD of these estimates were 0.93 ± 0.39, 0.97 ± 0.46, and 1.02 ± 0.57, respectively, for Regions 1, 2, and 3. Hence, there was no evidence of significant regional differences in average residual variances. The variability of these estimated scaling factors, however, tended to increase from Regions 1 to 2 to 3, such that the number of highly variable CG was greatest in Region 3. This result, however, may be due in part to the larger number of CG in this region (621 in Region 3 vs. 198 in Region 1 and 121 in Region 2), thereby increasing the probability of observing extremely highly or lowly variable CG. At any rate, regional differences in the degree of residual heteroskedasticity across CG may be possible. Our model could be readily extended to allow regional differences in the degree of residual heteroskedasticity with region-specific parameters (e.g., _region, region = 1, 2, and 3) if needed.

Regional differences in residual heteroskedasticity also could be due to differences in the level of production and environmental quality as variances tend to be proportional to mean responses for growth traits in beef cattle (Koots et al., 1994). Farms belonging to Region 1 were generally located in poorer environments compared with Regions 2 and 3 (i.e., the mean PWG were 82.4 ± 21.2 kg, 105.5 ± 39.5 kg, and 107.7 ± 39.3 kg for Regions 1, 2, and 3, respectively). Larger posterior means of v tended to be associated with larger CG-specific posterior means of u, with an estimated correlation coefficient between these two variables of 0.40 (i.e., CG with above-average performances also tended to be more variable). A positive correlation estimate of a similar nature, but for random genetic rather than for CG effects, has been reported for litter size in sheep (San Cristobal-Gaudy et al., 2001), whereas, conversely, a negative correlation estimate was observed for litter size in pigs by Sorensen and Waagespeterson (2003). Both studies modeled this relationship by specifying u and log(v) for the same factor as multivariate normal with a genetic correlation between the two classes of random effects.

Variance Components and Heritabilities

Despite the same variance-covariance specifications for u as based on across the six different models employed to analyze PWG, inference on the genetic components of changed considerably depending on the residual distributional specification. Posterior inferences of elements of and for the marginal residual variance (²_E) based on homoskedastic error and heteroskedastic error models are presented in Tables 2 and 3, respectively.

Among the genetic variance components, the segregation variance ²_S12 was the least affected by the different model specifications, as its posterior mean and mode were quite similar across the four models (Tables 2 and 3). Conversely, breed specific genetic variances were widely affected by the heteroskedastic versus homoskedastic error specifications. Based on point and interval estimates, Herefords seemed to have a substantially larger genetic variance than Nelores (i.e., ²_A2 ²_A1 ), under both homoskedastic error models (Table 2), whereas the opposite (i.e., ²_A1 ²_A2 ) was observed using heteroskedastic error models (Table 3).

Despite a seemingly nonsignificant influence of breed proportion on genetic variance under both heteroskedastic error models (Table 3), the effect of heteroskedasticity on breed group-specific phenotypic variances (²_Pg) was nevertheless quite appreciable, as seen from Table 4. The wide 95% PPI for ²_Pg of the Nelore breed obtained using both heteroskedastic error models indicated poor precision for inferring upon this parameter. This result was anticipated because purebred Nelores were only represented by parents without records, such that information to estimate that breed group-specific genetic and residual variance derived primarily from their crossbred progeny. Greater uncertainty on the Nelore ²_Pg using either heteroskedastic error model then better reflects this data structure limitation compared with the relatively sharp 95% PPI obtained using the homoskedastic error models also observed from Table 4. Conversely, for breed groups with substantial amounts of data, corresponding phenotypic variance estimates should be rather constant across candidate models. This was, to some extent, observed on Herefords and Nelores (Brafords) in Table 4. However, the largest discrepancy in posterior means of phenotypic variance between models within breed groups with records was 9.6% for F₁ cattle between the two Student’s t error models; this result is most likely due to the significant effect of ₂ on the residual variance of F₁ animals in the heteroskedastic error model.

View this table: [in this window] [in a new window]   Table 4. Posterior mean (PMEAN), posterior SD (PSD), posterior mode (PMODE), 95% posterior probability intervals (PPI), and effective sample size (ESS) of phenotypic variances (2Pg) as estimated for postweaning gain of different genotypes (g), obtained by four different models

View larger version (26K): [in this window] [in a new window]   Figure 1. Posterior density of heritabilities of postweaning gain for different breed composition groups: Nelore, Hereford, F1, and advanced Nelore (A38), obtained by a) Gaussian homoskedastic, b) Student’s t homoskedastic, c) Gaussian heteroskedastic, and d) Student’s t heteroskedastic error models.

We previously reported substantial rerankings of animals for genetic evaluations between statistical models that either do or do not account for heterogeneous genetic variances across multiple-breed groups (Cardoso and Tempelman, 2004). We observed rankings to be even more discrepant between homoskedastic error and heteroskedastic error models, in agreement with Hill (1984), even when heterogeneous genetic variances are modeled across breed groups in both models (data not shown).

	Implications

Top Abstract Introduction Materials and Methods Results and Discussion Implications Appendix Literature Cited

 
Our proposed genetic evaluation models help identify factors (e.g., sex, contemporary group, breed group) that are responsible for high vs. low variability of economically important traits, while at the same time providing results that are robust to data outliers. The implications of the use of these models relative to conventional genetic evaluation models that assume constant variability across environments were demonstrated to be potentially large for breed group-specific heritability estimates in a multiple-breed cattle population. Our proposed models will facilitate the identification of management systems and contemporary groups with superior product uniformity, an increasingly important beef marketing issue. Substantial rerankings in expected progeny differences for individual animals should be expected between genetic evaluation systems that either do or do not account for heterogeneous genetic and residual variances across environments.

	Appendix

Top Abstract Introduction Materials and Methods Results and Discussion Implications Appendix Literature Cited

 
The use of MCMC procedures requires the specification of full conditional densities (FCD) of each unknown parameter or blocks thereof. We briefly characterize each such FCD here and briefly describe strategies for sampling from those FCD that are not recognizable. A more extensive presentation of the derivation of the FCD for a heteroskedastic error model with only one fixed factor and one random factor is presented in Kizilkaya and Tempelman (2005).

We write for ²_ej defined in Eq. [2]. With flat bounded priors on identifiable elements of ß, the FCD of can be seen to be

(A1)

where

provided that X is full rank and

with 1 being a n x 1 unitary vector. Univariate sampling strategies for drawing from Eq. [A1] are provided in Wang et al. (1994).

Using the prior density in Eq. [4], the FCD for ²_e can be demonstrated to have an inverted gamma density with parameters . Similarly, other parameters defining residual heteroskedasticity have inverted gamma FCD. For example, the FCD for _ki or the ith ordered element of _k that specifies effects for factor k, k = 1, 2, ..., s can be shown to be inverted gamma with parameters is an indicator variable such that ^[k]_ij = 1 if the subscript address for the element of interest having row address i within _k is the same as that associated with observation j; otherwise ^[k]_ij = 0. Furthermore, indicates the number of occurrences that _ki was associated with a particular observation, and specifies the product of all elements of = [ ₁' ₂' ... ']', except _ki, that connect with y_j. Similarly, the FCD of _li the ith ordered effect in v_l can be shown to be inverted gamma with parameters is an indicator variable, such that ^[l]_ij = 1 if the subscript address for the element of interest having row address i within v_l is the same as that associated with observation j; otherwise ^[l]_ij = 0. In addition, is defined analogously to , but for effects in v = [v₁' v₂'... v_r']', connected to y_j and not including _li.

The FCD of _m, m = 1, 2, ..., p is not as readily recognizable as the other components of Eq. [2]:

(A2)

where _–m denotes all other elements of _m except _m. Hence, we sampled from Eq. [A2] using a Metropolis-Hastings (MH) step. This step was based on a random walk algorithm (Chib and Greenberg, 1995) using an inverted gamma distribution as the proposal density with scale parameter equal to the last sample value of _m times the shape parameter of the proposal density. This shape parameter was tuned during the first half of the burning period, such that the acceptance of proposed values was intermediate for optimal MCMC mixing (Chib and Greenberg, 1995).

The FCD for the weights w_j depend on the choice of p(w_j | ). Adopting the specification in Eq. [8], thereby leading to the Student’s t specification, these FCD correspond to a series of gamma densities with parameters represents the numerator of Eq. [2]. On the other hand, adopting Eq. [9] leads to a marginal Slash residual density, the corresponding FCD are truncated gamma densities (0 < w_j < 1), with parameters .

The FCD of _l also is not recognizable. That is, it can be shown (see also Kizilkaya and Tempelman, 2005) that

(A3)

where _–l refers to all elements of except _l and p(_l) is an arbitrarily defined prior. Here, Eq. [A3] does not have recognizable form and also requires MH sampling. In this case, we sampled _l = log(_l) using a random walk scheme with a Gaussian proposal density centered at the value of _l in the previous cycle and with the variance of the proposed density tuned at the midpoint of MCMC burn-in for optimal MCMC mixing (Chib and Greenberg, 1995).

The FCD for the robustness parameter _e depends on the choices of p(w_j |_e) and, of course, on p(_e). Under the Student’s t specification (i.e., p(w_j | _e) as in Eq. [8]), the FCD of _e is given by

(A4)

which does not have a recognizable form regardless of the specification on p(_e) but can be sampled using a MH algorithm similar to that adapted for Eq. [A2] and [A3].

Alternatively, adopting the Slash specification with p(w_j |_e) as provided in Eq. [9] and a conveniently chosen conjugate Gamma prior (see also Rosa et al., 2003), the FCD of _e also can be shown to be gamma with parameters .

The reader interested in the structure of and strategies for sampling from the the FCD for various elements of in a multiple-breed animal model is referred to the Appendix of Cardoso and Tempelman (2004).

	Footnotes

 
¹ Based on the first author’s Ph.D. dissertation, presented to Michigan State Univ. Funded by CAPES—Brasília/Brazil and the Michigan Agric. Exp. Stn. (Project MICL 01822). We are grateful to I. Misztal for making available Sparsem90 and Fspak90, and to the Brazilian beef cattle genetic evaluation alliance, "Conexão Delta G," for providing the Nelore-Hereford data.

² Correspondence: 1205 Anthony Hall (phone: 517-355-8445; fax: 517-353-1699; e-mail: tempelma{at}msu.edu).

Received for publication November 15, 2004. Accepted for publication March 21, 2005.

	Literature Cited

Top Abstract Introduction Materials and Methods Results and Discussion Implications Appendix Literature Cited

 


Birchmeier, A. N., R. J. C. Cantet, R. L. Fernando, C. A. Morris, F. Holgado, A. Jara, and M. S. Cristal. 2002. Estimation of segregation variance for birth weight in beef cattle. Livest. Prod. Sci. 76:27–35.

Cardoso, F. F., and R. J. Tempelman. 2004. Hierarchical Bayes multiple-breed inference with an application to genetic evaluation of a Nelore-Hereford population. J. Anim. Sci. 82:1589–1601.[Abstract/Free Full Text]

Chib, S., and E. Greenberg. 1995. Understanding the Metropolis-Hastings algorithm. Am. Stat. 49:327–335.

Clayton, D. G. 1996. Generalized linear mixed models. Pages 275–302 in Markov Chain Monte Carlo in Practice. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, ed. Chapman & Hall, London, U.K.

Fernandez, C., and M. F. J. Steel. 1998. On Bayesian modeling of fat tails and skewness. J. Am. Stat. Assoc. 93:359–371.

Foulley, J. L., M. S. Cristobal, D. Gianola, and S. Im. 1992. Marginal likelihood and Bayesian approaches to the analysis of heterogeneous residual variances in mixed linear gaussian models. Comput. Stat. Data. Anal. 13:291–305.

Garrick, D. J., E. J. Pollak, R. L. Quaas, and L. D. Van Vleck. 1989. Variance heterogeneity in direct and maternal weight traits by sex and percent purebred for simmental-sired calves. J. Anim. Sci. 67:2515–2528.

Gelfand, A. E. 1996. Model determination using sampling-based methods. Pages 145–161 in Markov Chain Monte Carlo in Practice. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, ed. Chapman & Hall, London, U.K.

Geyer, C. J. 1992. Practical Markov chain Monte Carlo. Stat. Sci. 7:473–511.

Gianola, D., J. L. Foulley, R. L. Fernando, C. R. Henderson, and K. A. Weigel. 1992. Estimation of heterogeneous variances using empirical Bayes methods—Theoretical considerations. J. Dairy Sci. 75:2805–2823.[Abstract]

Hill, W. G. 1984. On selection among groups with heterogeneous variance. Anim. Prod. 39:473–477.

Kizilkaya, K., and R. J. Tempelman. 2003. An assessment of heteroskedastic t-error linear mixed models for the analysis of field data collected from diverse environments. Page 62 in Proc. Int. Biometric Soc. ENAR Spring Mtg., Tampa, FL.

Kizilkaya, K., and R. J. Tempelman. 2005. A general approach to mixed effects modeling of residual variances in generalized linear mixed models. Genet. Sel. Evol. 37:31–56.[Medline]

Koots, K. R., J. P. Gibson, C. Smith, and J. W. Wilton. 1994. Analyses of published genetic parameter estimates for beef cattle production traits. 1. Heritability. Anim. Breed. Abstr. 62:309–338.

Lange, K. L., and J. S. Sinsheimer. 1993. Normal/independent distributions and their applications in robust regression. J. Comput. Graph. Stat. 2:175–198.

Lerner, I. M. 1954. Genetic Homeostasis. Oliver and Boyd, Edinburgh, U.K.

Lo, L. L., R. L. Fernando, and M. Grossman. 1993. Covariance between relatives in multibreed populations—Additive-model. Theor. Appl. Genet. 87:423–430.

Pinheiro, J. C., C. H. Liu, and N. Y. Wu. 2001. Efficient algorithms for robust estimation in linear mixed effects models using the multivariate t distribution. J. Comput. Graph. Stat. 10:249–276.

Robinson G. K. 1991. That BLUP is a good thing—The estimation of random effects. Stat. Sci. 6:15–51.

Rodriguez-Almeida, F. A., L. D. Van Vleck, L. V. Cundiff, and S. D. Kachman. 1995. Heterogeneity of variance by sire breed, sex, and dam breed in 200-day and 365-day weights of beef-cattle from a top cross experiment. J. Anim. Sci. 73:2579–2588.[Abstract]

Rosa, G. J. M., C. R. Padovani, and D. Gianola. 2003. Robust linear mixed models with normal/independent distributions and Bayesian MCMC implementation. Biometrical J. 45:573–590.

San Cristobal, M., J. L. Foulley, and E. Manfredi. 1993. Inference about multiplicative heteroskedastic components of variance in a mixed linear Gaussian model with an application to beef-cattle breeding. Genet. Sel. Evol. 25:3–30.

San Cristobal-Gaudy, M., L. Bodin, J.-M. Elsen, and C. Chevalet. 2001. Genetic components of litter size variability in sheep. Genet. Sel. Evol. 33:249–271.[Medline]

Sorensen, D. A., S. Andersen, D. Gianola, and I. Korsgaard. 1995. Bayesian inference in threshold models using Gibbs sampling. Genet. Sel. Evol. 27:229–249.

Sorensen, D. A., and D. Gianola. 2002. Likelihood, Bayesian and MCMC methods in Quantitative Genetics. Springer-Verlag, New York, Inc., New York, NY.

Sorensen D. A., and R. Waagepetersen. 2003. Normal linear models with genetically structured residual heterogeneity: A case study. Genet. Res. 82:207–222.[Medline]

Stranden, I., and D. Gianola. 1998. Attenuating effects of preferential treatment with Student-t mixed linear models: A simulation study. Genet. Sel. Evol. 30:565–583.

Stranden, I., and D. Gianola. 1999. Mixed effects linear models with t-distributions for quantitative genetic analysis: A Bayesian approach. Genet. Sel. Evol. 31:25–42.

von Rohr, P., and I. Hoeschele. 2002. Bayesian QTL mapping using skewed Student-t distributions. Genet. Sel. Evol. 34:1–21.[Medline]

Wang, C. S., J. J. Rutledge, and D. Gianola. 1994. Bayesian-analysis of mixed linear-models via Gibbs sampling with an application to litter size in Iberian pigs. Genet. Sel. Evol. 26:91–115.

This article has been cited by other articles:

				F. F. Cardoso, G. J. M. Rosa, and R. J. Tempelman Accounting for outliers and heteroskedasticity in multibreed genetic evaluations of postweaning gain of Nelore-Hereford cattle J Anim Sci, April 1, 2007; 85(4): 909 - 918. [Abstract] [Full Text] [PDF]

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Cardoso, F. F. Articles by Tempelman, R. J. Search for Related Content PubMed PubMed Citation Articles by Cardoso, F. F. Articles by Tempelman, R. J.