Comparison between estimation of breeding values and fixed effects using Bayesian and empirical BLUP estimation under selection on parents and missing pedigree information

Bayesian (via Gibbs sampling) and empirical BLUP (EBLUP) estimation of fixed effects and breeding values were compared by simulation. Combinations of two simulation models (with or without effect of contemporary group (CG)), three selection schemes (random, phenotypic and BLUP selection), two levels of heritability (0.20 and 0.50) and two levels of pedigree information (0% and 15% randomly missing) were considered. Populations consisted of 450 animals spread over six discrete generations. An infinitesimal additive genetic animal model was assumed while simulating data. EBLUP and Bayesian estimates of CG effects and breeding values were, in all situations, essentially the same with respect to Spearman's rank correlation between true and estimated values. Bias and mean square error (MSE) of EBLUP and Bayesian estimates of CG effects and breeding values showed the same pattern over the range of simulated scenarios. Methods were not biased by phenotypic and BLUP selection when pedigree information was complete, albeit MSE of estimated breeding values increased for situations where CG effects were present. Estimation of breeding values by Bayesian and EBLUP was similarly affected by joint effect of phenotypic or BLUP selection and randomly missing pedigree information. For both methods, bias and MSE of estimated breeding values and CG effects substantially increased across generations.


INTRODUCTION
Wang et al. [22] stated that one deficiency in the practical application of best linear unbiased estimation (BLUE) and best linear unbiased prediction (BLUP) is that errors of estimation of dispersion parameters are not taken into account when predicting breeding values. A two-stage estimation procedure (empirical BLUE/BLUP [5,8] (EBLUP)) is usually applied by first estimating variance components and then obtaining BLUE and BLUP of fixed and random effects, respectively, by replacing the parametric values of variance components by, usually, their restricted maximum likelihood (REML) estimates into the Mixed Model Equations (MME) [6]. Under random selection or absence of selection this EBLUP procedure converges in probability to BLUE and BLUP as the information in the data about variance components increases [22] and the distributions of variance components are symmetric and peaked [2,8]. The frequentist proprieties of EBLUP procedure under nonrandom selection are unknown [22].
The mean of the posterior distribution of breeding values can be viewed as a weighted average of BLUP predictions where the weighting function is the marginal posterior density of the heritability [5,15,16]. Estimation of breeding values by giving all weight to a REML estimate of heritability has been given theoretical justification [3]. When the information in the data about heritability is large enough, the marginal posterior distribution of this parameter should be symmetric and peaked. The modal value of the marginal posterior distribution should be a close approximation of its expected value. In this case, the posterior distribution of the breeding values can be approximated by replacing the unknown heritability by its REML estimate and an EBLUP procedure should yield a good approximation of the expected value of the marginal distribution of the breeding values.
Selection may increase the mean square error of the estimates of variance components [12] amplifying the uncertainty about genetic parameters. Gianola and Fernando [2], Wang et al. [21] and Sorensen et al. [15,16] advocated that Bayesian methods can fully take into account the uncertainty about dispersion parameters by considering the marginal posterior density of those parameters. Although the Bayesian methods provide an attractive theoretical framework for this problem, the practical benefits in prediction accuracy and precision are not clear. A comparison between sampling properties of EBLUP and Bayesian procedures under different scenarios including random and selected populations would be of interest.
The objectives of this study were to examine the effects of non-random selection on the parents (using phenotypic records or BLUP of breeding values) on the sampling properties of EBLUP and Bayesian estimates of breeding values assuming models with or without effects of contemporary groups, and to examine the impact of missing pedigree information on these two alternative methods.

Data simulation
Data were generated using a stochastic procedure similar to that described in [10,14,19,20]. This simulation procedure was simple and fully discussed in the literature. The genetic model assumed a large number of unlinked loci contributing to the genetic variance of a single hypothetical metric trait. The base population consisted of 10 males and 40 females which were assumed to be unrelated, unselected, and randomly sampled from a conceptually infinite population. The base animals were mated at random (four females per male) to produce 40 males and 40 females of generation 1. Ten males were selected as parents for the next generation following one of three schemes, i.e., random selection, selection on the basis of highest phenotypes, and selection on the basis of highest estimated breeding values. The last two gave different degrees of selection for true merit.
Therefore, selection was only on males and generations were discrete. Six generations were simulated, including the base population. No attempt was made to control inbreeding.
The model for simulation of data was Y ij is the phenotypic observation of animal j in contemporary group (CG) i, b i is the effect of CG i, a ij is the additive genetic value of animal j in CG i, and e ij is the random residual term.
Values for e ij were independently drawn from a normal distribution with mean zero and variance σ 2 e . The additive genetic variance for the base population, before selection, was σ 2 a . Genetic values of base animals were independently drawn from a N(0, σ 2 a ). Genetic values of animals in later generations were simulated as a ij = 1 2 (a sj + a dj ) + m ij , where a sj and a dj are genetic values of the sire and dam of individual j, and m ij is the Mendelian sampling effect of individual j assumed to be independent of the genetic values of the sire and dam. The inbreeding coefficient (F) of the parents was taken into account, so that m ij was drawn from a N 0, 1 2 − 1 4 (F sj + F dj ) σ 2 a . Two models were used. The first model did not include CG effects, in which case b i was equal to 0 for every i. This model was denoted as RM for random simulation model. The second model, called mixed simulation model (MM), included CG effects that were simulated in the first replicate and kept constant for all replicates. Eight CG's were assigned per generation, four for males and four for females. Their effects (b i ) were drawn from a uniform distribution ranging from −5.5 to +5.5. Animals were assigned randomly to CG's within generation and sex in each replicate. Connectedness of CG's was guaranteed by requiring two sires to have progeny in all eight CG's within a generation, and guaranteed a minimum of two animals per CG.
Pedigree information was either complete or had 15% randomly chosen nonbase animals with both sire and dam declared missing. Low (0.2) and high (0.5) heritability values were used in the simulations. The sum of the genetic and residual variances was kept at 20.0. The genetic variance was either 4.0 or 10.0, and the residual variance was either 16.0 or 10.0, respectively. One hundred replicates were simulated for each combination of model, selection scheme, heritability level, and pedigree information, and each replicate included 400 animals with phenotypic records plus 50 base population animals without records.

Analyses
The operational model was defined to be the same as the true model used for simulation of a data set. An overall mean (µ) was included in the model for RM data sets because the phenotypic mean was unlikely to be zero in the selected populations. The univariate linear mixed model used to analyze the simulated data was: The distributional assumptions were: a ∼ N(0, Aσ 2 a ) and e ∼ N(0, Iσ 2 e ), where a is the vector of additive genetic effects and e is the vector of random residual effects, and A was the numerator relationship matrix that included base population animals and accounted for inbreeding.
REML estimates of variance components were obtained from the multiple trait derivative free programs of Boldman et al. [1]. The starting values of the variances were the true simulation values.

Bayesian analyses
Bayesian estimates were obtained via Gibbs sampling following Wang et al. [22] and Van Tassel et al. [20]. In addition to the previously mentioned assumptions about distributions, prior densities (PD) were assigned for all variance components and the location parameters µ and b. Two different priors were assumed for b: a flat improper prior, where p(b) ∝ constant (p() denotes a density function), indicating no prior knowledge about their effects (fixed b in a frequentist setting), or a proper prior, where b ∼ N(0, Iσ 2 b ). The overall mean µ was always assumed fixed, that is p(µ) ∝ constant. For the variance components independent scaled inverted chi-square distributions (χ −2 i ) were assumed: where ν i is a degree of belief parameter and s 2 i can be thought of as a prior value for the variance.
The joint posterior density of all unknowns (Θ, v) was where Θ = (µ, b , a ) are the location parameters, v = (σ 2 b , σ 2 a , σ 2 e ) are the variances, and ν = (ν b , ν a , ν e ), and s = (s 2 b , s 2 a , s 2 e ) are parameters describing the prior degrees of belief and prior variances, respectively. When a flat improper prior was assumed for b, (1) did not apply for σ 2 b and (2) did not involve the term related with σ 2 b . For ν i a prior value of 5 was used for all variances. This value was chosen so that the variance of the prior scaled inverted chi-square distribution was large but finite. Given the value for ν i , the prior values for s 2 i were specified such that expected values of the prior scaled inverted chi-square distribution were equal to the true values. These prior values used for ν i and s 2 i yielded prior coefficients of variation equal to 141.4% for any i and heritability level.

Gibbs sampling
The fully conditional posterior distributions for the location parameters were normal. Let Θ −i be Θ without its i-th element and v −i to be v without its i-th element, then where w ij is the element ij of the coefficient matrix and h i is the element i of the right-hand side of the MME.
The fully conditional posterior distribution of variance components were in the scaled inverted chi-square form. For σ 2 e it was with parameters For the other variance components (σ 2 i ) it was The previous fully conditional posterior distributions from (3) to (5) were used in the Gibbs sampling scheme. The starting values of the variances to obtain the first solution from MME were the true simulated values.
The Gibbs sampling loop was repeated 10 000 times. A burn-in period of 1 000 rounds was used and was based on previous analyses where the plots of all samples were subjectively evaluated for trend and variability.

Posterior parameter estimates
All samples after the burn-in period were used to estimate the posterior mean of the distribution of the location parameters. Therefore, breeding values and CG effects were evaluated at their posterior mean value.

Empirical BLUE/BLUP analyses
The MME were used to predict breeding values and to estimate CG effects. The true variances were replaced by the REML estimates. The models were the same as used for the Bayesian analyses.

Criteria for comparing methods
Methods were compared based on their biases, mean square errors (MSE), and Spearman's rank correlations of predicted breeding values and estimated CG effects with respect to their true values. The rank correlation was used as an attempt to measure the ability of each method in properly ranking animals and environmental effects.
Bias and MSE were defined, respectively, as the average deviation and the average squared deviation of predicted breeding values from their corresponding true values or of estimated contrasts of CG effects from their corresponding true values: Where q is the number of animals or the number of CG's, andˆrefers to the predicted or estimated value of the parameter ω.
Because an overall mean was included in all analyses, the effects of CG were not estimable when they were treated as fixed effects. Thus, the estimable contrasts between each level of CG effect and the first level were used to calculate the rank correlation, bias, and MSE for all analyses.
The differences in biases, MSE and rank correlations between methods were tested by a paired t-test [9,13] at the 5% significance level. For bias, the paired t-test was not performed when the biases of both methods were not significantly different from zero by a t-test.

Spearman's rank correlations
The results presented in Table I and Table II (for low and high heritability, respectively) showed that there was no difference between Bayesian and EBLUP estimation regarding the overall rank correlation of breeding values and of estimable contrasts of CG effects with their true values for any combination of simulation model, selection scheme, true heritability (h 2 ), and level of pedigree information (PI). Rank correlations were also calculated within each generation (data not shown) and there were no differences between the two procedures across all simulated scenarios.
Bayesian and EBLUP estimation yielded rank correlations between true and predicted breeding values that were equally decreased by randomly missing PI and by both phenotypic and BLUP selection. The joint effect of selection and missing PI produced the smallest rank correlations for both RM and MM data sets.
For all analyses, regardless of the true heritability, the rank correlations between Bayesian and EBLUP estimates of breeding values and of contrasts of CG effects were higher than 0.998 (data not shown). Table I Selection and missing PI did not affect rank correlations between true and estimated contrasts of CG effects. The insensitivity of rank correlations between true and estimated contrasts of fixed CG effects to missing PI [4,7] and to phenotypic selection, which was characteristically not translation invariant to the fixed effects [7], was not expected. The simulation procedure may have facilitated the estimation of CG effects for several reasons. First, animals were assigned randomly to CG in each generation. Thus, great differences in genetic mean among CG's from the same generation were not expected. Larger differences may be found with real data. Second, sires were selected across CG's, but within each discrete generation. Finally, the average number of animals within CG levels (10) was large enough to allow estimation of their effects with reasonable accuracy [18]. With real data, some CG's are often smaller and especially the variability of CG size is usually much larger.
Use of proper informative priors for CG effects and their variance or considering CG effects as random in EBLUP analyses had negligible effect on rank correlations of breeding values.

Biases
Table III presents the empirical mean over 100 replicates of the biases in each generation and over generations for high h 2 . Bayesian and EBLUP estimation showed the same pattern regarding the biases of both predicted breeding values and estimated contrasts of CG effects. The small differences between biases of the two methods were (Tab. II), however, often significant (p < 0.05). For low h 2 , similar results were found.
Phenotypic and BLUP selection did not cause bias on predicted breeding values from Bayesian and EBLUP analyses when pedigree information was complete (Tabs. I and II).
Nonrandom selection in conjunction with 15% randomly missing PI had large impact on biases of estimates from both procedures for RM and MM data sets (Tab. III, analyses 5 and 6, and 11 and 12, respectively). In these cases, biases in breeding values increased negatively and consistently as generation number increased. For both, phenotypic and BLUP selected populations, the bias in the last generation was around 29% and 34% of the true additive genetic mean of the population at this generation for RM and MM data sets, respectively. For the case of full PI, the same figures were 2% and 1%. For low h 2 , the biases were 23% and 28%, and 1% and 2% for the cases of missing and full PI, respectively.
In the MM data sets, the increase of bias in the estimated contrasts of CG effects was in the opposite direction (positive) to that of the breeding values. When changes in the expectations of genetic values are not modeled through a complete additive relationship matrix in an animal model or the use of a genetic  grouping strategy [24], solutions for the genetic effects might be confounded by fixed effects, generating bias and increased MSE. The use of Bayesian procedure did not lessen the effects of not accounting for missing PI when non random selection was applied. These results reinforce the importance and need for properly account for missing PI regardless of the procedure used for estimation. Rodriguez et al. [11] gave one example of an application of Bayesian analyses with genetic groups in the model.

Mean square errors
Table IV presents the empirical mean over 100 replicates of the MSE in each generation and over generations for high h 2 . The mean MSE over generations of predicted breeding values and estimated contrasts of CG effects were usually smaller (p < 0.05) for Bayesian than for EBLUP analyses (Tab. II) although the differences were small. Selection associated with missing PI greatly increased the MSE of predicted breeding values and estimated contrasts of CG effects from both Bayesian and EBLUP analyses. Similar results were found for low h 2 .
As shown in Tables I and II, phenotypic and BLUP selection did not cause bias on predicted breeding values from Bayesian and EBLUP analyses when pedigree information was complete, but increased MSE, when MM data sets were analyzed (analyses 7 vs. 9 and 10).
Weigel et al. [23] investigated the improvement of fixed effect estimates in a mixed linear model and concluded that it was possible to improve upon unbiased estimators in a mean squared error sense by allowing bias. In agreement with Weigel et al. [23], treating CG's as random (Tabs. III and IV), shrunk their solutions towards zero, created some bias on CG estimates, but reduced the MSE of CG and, in less extent, of breeding value estimates in nonrandomly selected populations with full PI (analyses 9 and 10 vs. 15 and 16). With missing PI (analyses 11 and 12 vs. 17 and 18), the reduction in the MSE of breeding value estimates was more accentuate. With full PI, treating CG's as random introduced a small bias in the breeding values estimates (Analyses 7, 9, and 10 vs. 13, 15, and 16).

General discussion
The asymmetry of the marginal posterior distribution of σ 2 a , when there was random selection and full pedigree information, as illustrated in Figure 1 (analysis 7) for low and high h 2 , suggests that the simulated data sets did not have a high degree of resolution concerning inferences about genetic parameters. Sorensen et al. [16] argued that this fact is taken into account when computing the marginal posterior distribution of breeding values. This is in marked contrast with the estimation of breeding values that is obtained    using EBLUP, which assumes the h 2 known and gives 100% weight to an estimate of this h 2 . In this study, however, Bayesian estimation did not differ from EBLUP estimation regarding rank correlations with true values for both estimated contrasts of CG effects and predicted breeding values. With respect to bias and MSE, EBLUP and Bayesian estimation showed the same pattern over the range of simulated scenarios and exhibited only small differences in their values. The small differences in biases and MSE of the estimates from the two methods could be speculated to be due to the influence of the vague prior densities on variance components used in the Bayesian analyses. If this is the case, those small differences should disappear when larger, more informative data files were analyzed, because the likelihood function of the data would dominate the prior information.
Markov Chain Monte Carlo (MCMC) error of Bayesian estimates was indirectly assessed for the variance components on the basis of effective chain sizes, which ranged from 145 to 385 for all variances. The effective chain sizes were reasonably large to yield acceptable MCMC errors on the posterior means.
The great appeal of the Bayesian analyses via Gibbs sampling is that it yields Monte Carlo estimates of the full marginal posterior distribution of all parameters of interest, for instance breeding values, from which the probabilities that the parameter lies between specified values can be computed [17,20]. This is particularly interesting when asymptotic normality of the posterior distributions is difficult to justify, which can be the case with selected populations [16] when the variance components are not known. In this case, the uncertainty about the variance components is accounted for in the Bayesian probability intervals of predicted breeding values [22].
There are also situations where the infinitesimal model is not a sound approximation and, therefore, normality does not hold after cycles of selection. Bayesian analyses could be more flexible to incorporate more appropriate or robust distributions.
Bayesian analyses via Gibbs sampling are becoming more and more feasible as computer power increases and as better algorithms are developed. The applicability of Bayesian methods for genetic evaluation is already possible routinely for moderately sized problems.

CONCLUSIONS
Bayesian and EBLUP estimation did not differ over the range of simulated situations in this study with respect to Spearman's rank correlations between true and predicted breeding values and between true and estimated contrasts of CG effects. Hence, the two methods showed the same ability to rank animals and environmental CG effects. The sample properties, Bias and MSE, of Bayesian and of EBLUP estimates showed the same pattern over the range of simulated scenarios. The bias and MSE of Bayesian estimates were often less than of EBLUP estimates, but the differences were small and likely due to the vague prior information on variance components used in the Bayesian analyses.
Phenotypic and BLUP selection did not cause bias in predicted breeding values by Bayesian or EBLUP when pedigree information was complete, but caused small increases in MSE, when MM data sets were analyzed.
Bayesian and EBLUP prediction of breeding values were similarly affected by the joint effect of phenotypic or BLUP selection and randomly missing pedigree information. For both methods, bias and MSE of predicted breeding values and estimated contrasts of CG's substantially increased across generations, because the change in the expectation of breeding values was not accounted for in the model.