Comparison of inference methods of genetic parameters with an application to body weight in broilers

REML (restricted maximum likelihood) has become the standard method of variance component estimation in animal breeding. Inference in Bayesian animal models is typically based upon Markov chain Monte Carlo (MCMC) methods, which are generally flexible but time-consuming. Recently, a new Bayesian computational method, integrated nested Laplace approximation (INLA), has been introduced for making fast non-sampling-based Bayesian inference for hierarchical latent Gaussian models. This paper is concerned with the comparison of estimates provided by three representative programs (ASReml, WinBUGS and the R package AnimalINLA) of the corresponding methods (REML, MCMC and INLA), with a view to their applicability for the typical animal breeder. Gaussian and binary as well as simulated data were used to assess the relative efficiency of the methods. Analysis of 2319 records of body weight at 35 days of age from a broiler line suggested a purely additive animal model, in which the heritability estimates ranged from 0.31 to 0.34 for the Gaussian trait and from 0.19 to 0.36 for the binary trait, depending on the estimation method. Although in need of further development, AnimalINLA seems a fast program for Bayesian modeling, particularly suitable for the inference of Gaussian traits, while WinBUGS appeared to successfully accommodate a complicated structure between the random effects. However, ASReml remains the best practical choice for the serious animal breeder.


Introduction
The restricted maximum likelihood (REML) method (Patterson and Thompson, 1971) for unbalanced mixed models has been extensively used in animal breeding and has become the standard method for the estimation of variance components.The Bayesian Markov chain Monte Carlo (MCMC) methods were introduced in quantitative genetics in the early 1990s (Wang et al., 1993;Sorensen et al., 1994), facilitated by the development of the Gibbs sampling procedure (Geman and Geman, 1984;Gelfand and Smith, 1990).The Gibbs sampler successively samples from conditional distributions of all parameters in a model in order to generate a random sample of the marginal posterior distribution, which is the target for Bayesian inference.MCMC methods represent the standard inference procedure for Bayesian animal models (Sorensen and Gianola, 2002), and through the years they have become an attractive alternative to REML.Recently, a non-samplingbased alternative to MCMC, the integrated nested Laplace approximations (INLAs), has been introduced (Rue et al., 2009).Using INLA, marginal posteriors for all parameters and random effects can be calculated.Because INLA is based on direct numerical integration instead of simulations, it is much faster than MCMC (Rue et al., 2009).Furthermore, Holand et al. (2013) have developed an R package (Anima-lINLA) making Bayesian animal models more accessible to animal breeders.
Published by Copernicus Publications on behalf of the Leibniz Institute for Farm Animal Biology.
Several programs are available for MCMC methods, but very few provide a flexible environment.WinBUGS (Lunn et al., 2000) is the most well-developed and general-purpose Bayesian software available to date.It has an interactive environment that enables the user to specify models that need to be compiled before starting the Gibbs sampling.Convergence diagnostics, model comparisons, e.g., via DIC (deviance information criterion), and other useful plots and diagnostics are available.Several distributions can be used for modeling the observations as well as priors, while full conditional distributions are automatically constructed and the appropriate MCMC algorithm for sampling is chosen (Lunn et al., 2000).In WinBUGS and in the context of animal breeding, an important issue is the importation of the animals' genetic relationship matrix.Methods proposed so far (Damgaard, 2007;Waldmann, 2009) either require prior transformation of the data using complex code or do not provide a generic procedure independent of the data structure.A good solution here is the use of the inverse of the numerator relationship matrix A −1 directly through the diagonal values of W −1 matrix, where A −1 = (T −1 ) W −1 T −1 (Henderson, 1976;Quaas, 1989), as suggested by Gorjanc (2010).Recently, Hallander et al. (2010) have developed a Bayesian method in WinBUGS based on the decomposition of the multivariate normal prior distribution into products of conditional univariate distributions, thus permitting the genetic evaluation of complex pedigree structures.In addition, more complicated covariance structures have been incorporated via Bayesian methods, allowing for the simultaneous estimation of both additive and dominance genetic effects (Waldmann et al., 2008;Mathew et al., 2012).
The primary goals of the present study were to apply and investigate the relative merits of three methods (REML, Gibbs sampling and INLA) in the context of animal breeding, using representative programs such as ASReml 3.0 (Gilmour et al., 2009), WinBUGS and AnimalINLA.For this purpose, both a Gaussian and a binary trait were explored and variance components and the genetic parameters along with breeding values across the three methods were estimated and compared.

Data description
Data on body weight (BW) at 35 days of age from a broiler line were made available by Aviagen Ltd.Given that, in the Windows version of AnimalINLA 1.1, limitations in the size of the data set exist, a small data set was randomly selected, consisting of 2319 records.This comprised 1171 males and 1148 females in 40 hatch weeks, while the pedigree included a total of 2456 animals.All sires (n = 32) and dams (n = 105) were assumed to be non-inbred and nonrelated.To make results directly comparable, all phenotypic values were standardized to the standard normal distribution via y = y 0 − ȳ σ y 0 , where y ∼ N(0, 1) is the standardized BW, y 0 the original phenotypic values of BW, ȳ the mean BW in the population and σ y 0 the standard deviation of BW.A preliminary analysis of variance showed that the statistically significant (P <0.05) fixed effects included hatch week and sex.Hence, these fixed effects were included in all models.In this data set, each dam was mated with two sires producing from 2-57 offspring with records (average full-sib family size: 16), while sires were mated with two to seven dams and produced 2-97 offspring (average half-sib family size: 56).Such a structure enabled the inclusion of maternal environmental effects (c 2 ) through proper modeling.The latter are modifications of the offspring phenotype caused by the environment provided by the mother and consider any influence of a dam on its progeny, excluding the effects of directly transmitted genes.
A binary response trait was also constructed, using the original BW values and a threshold at the highest 20 % phenotypic values.Thus, the new variable y B followed the Bernoulli distribution, with values 0 and 1 denoting low and high weight, respectively.In this data set, only the gender of the animals was statistically significant (P <0.05) and was thus included in analyses as the only fixed effect.

Gaussian trait
Three animal models were considered for BW.Model M 1 was a purely additive animal model, while model M 2 allowed for the inclusion of maternal environmental effects and model M 3 was as model M 2 but with a covariance σ uc between additive genetic and maternal environmental effects.In summary, the models in matrix notation were as follows: where y = n × 1 is the vector of observations (n: number of records, 2319), b = p × 1 is the vector of fixed effects (p: number of fixed effects classes, 42), u = q × 1 is the vector of direct additive genetic effects (q: number of additive effects, 2456), c = k × 1 is the vector of maternal environmental effects (k: number of dams with offspring, 105) and e = n × 1 is the vector of residuals; X, Z and Z c denote the incidence matrices relating the observations to the corresponding fixed and random effects.The vector of direct genetic effects was assumed to follow the normal distribution: u ∼ N 0 n , σ 2 u A , where 0 n denotes a n × 1 vector of 0s, σ 2 u denotes the direct genetic variance and A denotes the additive genetic relationship matrix.The maternal environmental effects were assumed to follow a normal distribution given by c ∼ N 0 k , σ 2 c I k , where I k is an identity matrix of order k and σ 2 c is the maternal environmental variance.Finally, residuals for the two traits were assumed to be normal as follows: e ∼ N 0 n , σ 2 e I n , where σ 2 e is the residual variance.
From a Bayesian perspective, the data y are assumed to be y|b, u, σ 2 e ∼ N (Xb + Zu, σ 2 e I n ) and y|b, u, c, σ 2 e ∼ N(Xb + Zu + Z c c, σ 2 e I n ) for models M 1 and M 2 , respectively.The vector of the data y for model M 3 was assumed to be y|b, u, c, r, σ 2 e ∼ N(Xb + Zu + Z c c, σ 2 e I n ), where the correlation was r = cov(u,c) σ u σ c .The vector of b (p × 1) for all three models was partitioned into two sub-vectors, denoting hatch (h) and sex (s).It was assumed that both sub-vectors followed univariate normal, according to h|σ 2 h ∼ N(0, σ 2 h )I and s|σ 2 s ∼ N(0, σ 2 s )I.Gelman (2006) investigated the statistical properties of different priors on variance components and found that a uniform prior on the standard deviation is a reasonable choice in a number of situations.Therefore, vague uniform priors were utilized for the standard deviation of the additive genetic effects σ u ∼ U (0, 100) as well as for the c 2 effects σ c ∼ U (0, 100).The inverse gamma distribution (0.001, 0.001) for the residual variance σ 2 e or the uniform distribution σ e ∼ U (0, 100) for the residual standard deviation were utilized in order to account for the effect of the priors on the estimations.Both approaches gave indifferent results.The same priors were used in AnimalINLA and in WinBUGS to attain comparability.Inferences were made by REML and by estimating the marginal posterior distribution using either Gibbs sampling or INLA.Estimates of heritability (h 2 ) as well as c 2 were calculated as ratios of the estimates of direct additive genetic (σ 2 u ) and maternal environmental (σ 2 c ) variances, respectively, to the phenotypic variance (σ 2 p ).The phenotypic variance accounts for the sum of all variance components, according to the model.
For measuring the mixing and efficiency of the MCMC samples, the effective sample size (ESS) was used.The ESS of the posterior samples of each parameter corresponds to the number of independent samples having the same estimation accuracy as the dependent MCMC samples and is given by Waagepetersen et al. (2008) , where K the total number of correlated MCMC samples and ρ k is the Markov chain lag-k autocorrelation.

Binary trait
Initially, a simple animal model was fitted via REML, considering y B as a normally distributed trait.Subsequently, a generalized linear model (McCullagh and Nelder, 1994) was used for the analysis of the binary variable.In this analysis, the observed binary variable y B is related to an underlying unobservable continuous variable λ, such that the observed binary response (y B ) is the result of the following relationship: where τ is fixed and y Bi corresponds to observation i.Several link functions (logit, probit, cloglog) can be applied to link the binary variable to the underlying scale (Gilmour et al., 2009).In our study, the logit function was used: , where µ is the probability of success and λ the vector of linear predictors of the unobserved variable on the underlying scale.An animal model was assumed for λ such that λ = Xb + Zu + e.A uniform prior was assumed here for the standard deviation of the additive genetic effects on the underlying scale σ u ∼ U (0, 100).On the logit scale σ 2 e = π 2 3 ≈ 3.29, and heritability is thus estimated as Gilmour et al., 2009).
In order to investigate the relative merits of the three approaches, data for both the Gaussian and the binomial case were simulated and models were applied accordingly.

Simulation study
The initial analysis of data revealed a marginal importance of the c 2 effects and a possible covariance between u and c.To further test the behavior of the three programs under a scenario of two correlated random effects with a marginal contribution by one of them, a simulation study was conducted, emulating the pedigree structure and the variance components of the real data.In total, 20 sires and 70 dams were used in the pedigree, and 2240 progeny with records were simulated.Each sire was assumed to mate to seven dams, while each dam produced offspring with two different sires.All sires and dams were assumed to be non-inbred and non-related.Each full-sib family consisted of 16 offspring.The direct genetic effect for founder i (1, ..., 90) was drawn as u i ∼ N 0, σ 2 u , while the maternal environmental effect of dam j (1, ..., 70) was c j ∼ N 0, σ 2 c , with σ 2 u = 7 and σ 2 c = 3.Two scenarios were explored regarding the correlation between the direct genetic and the c 2 effects (r uc ): (a) r uc = −0.2(low) and (b) r uc = −0.8(high).The direct genetic effects of offspring i (1,... ,2240) were calculated by u i = 1 2 u j + u k +ms, where u j and u k denote direct genetic effects of dam and sire, respectively, while ms i represented the Mendelian sampling deviation drawn conditional upon the c 2 effects: ms . The total phenotypic variance was estimated according to σ 2 p = σ 2 u + σ 2 c + σ 2 e .The residuals were sampled as e i ∼ N 0, σ 2 e , where σ 2 e = 32, thus resulting in σ 2 p = 42, h 2 = 0.17 and c 2 = 0.07.In total, 30 samples from each scenario were generated.These samples were then analyzed via models M 1 -M 2 (ASReml and AnimalINLA) and M 2 -M 3 (WinBUGS).The mean squared error (MSE) was employed to quantify the performance of the predictors throughout, along with the coverage of interval estimates.The MSE was computed as follows: , where θ stands for the true and θi for the estimated parameter, θi − θ corresponds to bias, and N = 30 is the number of samples.

Model evaluation criteria
According to the method applied, the model comparison was based on four evaluation criteria: the Akaike information criterion (AIC; Akaike, 1973), the Bayesian information criterion (BIC; Schwarz 1978), the conditional Akaike information criterion (cAIC; Vaida and Blanchard, 2005) and the DIC (Spiegelhalter et al., 2002).All criteria are based upon the computation of the deviance (D): D = −2 log(p(y| θ )) = −2 log L, where θ denotes the p × 1 vector of the model parameters and p(y| θ) denotes the likelihood of the data y evaluated at the maximum likelihood estimate θ .While likelihood ratio tests (LRTs) suggest the direct comparison of logLs between the various nested models, AIC, BIC and cAIC suggest penalizing the deviance by appropriate complexity terms.However, the determination of the number of the model parameters is nontrivial when random effects are of interest and are being estimated using methods such as BLUP.For such cases the AIC is shown to be asymptotically biased (Crainiceanu and Ruppert, 2004).An asymptotically unbiased criterion is the cAIC, defined by Vaida and Blanchard (2005) as cAIC = −2 log L i + 2ρ, where ρ are the effective degrees of freedom (Hodges and Sargent, 2001), given by the trace of the hat matrix H that maps the vector of observed values to the vector of the fitted values.In all criteria, models with smallest values are to be preferred, denoting a better balance between complexity and fit.

Gaussian trait
Table 1 summarizes the estimated variance components and genetic parameters of BW, along with likelihoods, ρ and the model evaluation criteria.With regard to the Bayesian methods, posterior means and posterior medians were very close for all parameters of interest.The closeness of mean, median and mode was also suggested by visual inspection of the posterior densities, which displayed unimodality.Therefore, only the posterior means are presented.For our data to achieve convergence via WinBUGS, a burn-in of 10 000 iterations, a total number of 1 000 000 samples and a thinning interval of 20 were necessary.The latter was concluded on graphical inspection of the trace and autocorrelation plots, yielding a sample of 50 000 iterations.Such runs took approximately 14 to 16 h, depending on modeling assumptions.Heritability for BW ranged from 0.15 to 0.34, while c 2 ac-counted for 0-0.08 of the total phenotypic variance, depending on the model and the method applied.All evaluation criteria, regardless of the method considered, concur in the choice of a purely additive animal model without the inclusion of the c 2 effects.With M 1 , heritability estimates ranged slightly among the methods, from 0.31 (ASReml) to 0.34 (AnimalINLA), while 95 % confidence and credible intervals between ASReml and the Bayesian programs always coincided.The ESS of all parameters estimated via model M 1 and WinBUGS exhibited the highest values (higher than 7000) among models, indicating best MCMC mixing properties.
Under model M 2 , REML-based estimates were significantly different than those obtained from the two Bayesian approaches.In this case, REML heritability was seriously underestimated (0.15) when contrasted with MCMC and INLA methods (0.31 and 0.32, respectively).Furthermore, while c 2 was 0.07(±0.03) in REML, no detectable variance due to c 2 was identified with the Bayesian methods.As a result, the sum of the additive and the c 2 effects given as a proportion of the phenotypic variance was significantly lower in REML (0.22) when compared to Bayesian methods (0.31-0.32).Such a paradox may arise from covariances between the various random effects.To test for such a hypothesis, we fitted model M 3 that accounted for a covariance between the additive genetic and the maternal environmental effects.
This could be effectively modeled only via the WinBUGS software.Under model M 3 , h 2 and c 2 estimates were comparable (0.17 and 0.08, respectively) to ASReml estimates (for model M 2 ), while the covariance in question was not statistically significant (Table 1).A negative additive genetic maternal environmental correlation was detected (−0.20), although with large standard error (0.30) that did not allow for firm conclusions.
To further quantify the implications of model and method evaluation on selection decisions, Pearson as well as rank correlations of animals' EBVs and the percentage of common animals selected were calculated across the models and methods applied (results not shown).The correlations in question were extremely high (0.97-0.99) when the focus was on the whole population and/or a proportion of the best 20 % of animals.During this phase, an additional advantage of the WinBUGS software was its ability to estimate (via the rank tool) the uncertainty associated with the ranking of the individuals from the posterior distributions of the EBVs. Figure 1 presents 12 selected examples from the posterior distribution of the EBV ranks, with four animals each from the top, middle and low end of the spectrum.These ranks were based upon the whole posterior density and properly accounted for characteristics such as the variance and skewness of the posterior.Both, a 95 % rank interval as well as the median rank are provided, thus presenting an easy and flexible way of animal selection.The large uncertainty associated with selecting among similar animals is also illustrated.Here, rank correlations were remarkably high, ranging from 0.96 to 0.99 among all methods and models considered.Further-more, standard errors of the EBVs and solutions for the fixed effects were comparable among the methods, with no statistically significant differences.All models and methods suggested the same animals, resulting in correlations between the estimated breeding values that ranged from 0.96 to 0.99.

Binary trait
The estimated variance components and genetic parameters of y B for a purely additive animal model across the three methods are presented in Table 2.A model incorporating c 2 effects was also fitted; however, convergence was not achieved under any method applied.In ASReml, heritability on the observed scale (h 2 o ) was estimated to be as high as 0.10, while the respective estimate on the underlying scale was significantly higher (h 2 U = 0.19).Using the classical formula (Dempster and Lerner, 1950), the ratio between the two estimates would be , where p is the level of incidence and z x p is the ordinate of a standard normal curve cutting off an area equal to p.For p = 0.2 (as in here) the ratio is ( ) in full agreement with our results.Estimates from AnimalINLA were comparable to those of AS-Reml (h 2 U = 0.21).Interestingly, the WinBUGS heritability estimate was significantly higher (up to 0.36), exceeding the original h 2 .Differences were also detected on the 95 % confidence or credible intervals of the point estimates of the additive variance as well as the heritability on the underlying scale.More specifically, the 95 % credible interval of h 2 U given by WinBUGS was in the region of (0.21, 0.56), that of AnimalINLA was in (0.13, 0.30) and finally that of ASReml was in (0.09, 0.29).The ESS of all parameters estimated via WinBUGS were 1293 and 1436 for h 2 and additive genetic variance, respectively.
As in the case of the Gaussian trait, rank correlations across the three methods remained high, ranging from 0.92 to 0.99 (results not shown).In addition, the proportion of common animals selected among the three methods exceeded 93 %, suggesting minor implications of method usage on selection decisions.

Simulation study
Descriptive statistics of the simulated data and the estimators across models and methods are given in Table 3.Average values of the simulated data were equal to the true ones (h 2 = 0.17 and c 2 = 0.07).Note that during simulations, c 2 was statistically significant.Using model M 1 under either ASReml or AnimalINLA always resulted in inflated predictions for the true parameters.More specifically, the estimated heritability ranged from 0.35 to 0.51, with a tendency for inflated estimates particularly in AnimalINLA and under the strongly negative-r uc scenario for both software packages (ASReml and AnimalINLA).Overestimation of the heritability was due to both higher estimates of the additive genetic variance and lower estimates of the residual variances.
Estimates under model M 2 were in close proximity to the true values only in the case of ASReml and the low-r uc scenario (h 2 = 0.15, c 2 = 0.07).Slightly higher estimates for h 2 and c 2 were observed in ASReml in the high-r uc scenario (h 2 = 0.21, c 2 = 0.08).Under AnimalINLA, the respective h 2 estimator was seriously inflated (h 2 = 0.34) due to overestimation of the additive genetic effects and failure to account for the c 2 effects.This trend was more evident in the strong-vs.the low-r uc scenario.The WinBUGS estimates for Model M 2 under the high-r uc scenario were slightly better than those obtained by AnimalINLA.Finally, model M 3 was fitted via WinBUGS for the high (r uc = −0.8)scenario.In this case, a statistically significant r uc was detected (as high as −0.60), but h 2 and c 2 were systematically overestimated.Only minor differences were observed in the mean estimates using WinBUGS and two prior distributions for the residual variance.In Table 3, results are derived from the uniform distribution case.
The MSEs across models and methods are presented in Table 4. Irrespectively of the method and/or model, MSEs were lower in the low-vs.the high-correlation scenario.Furthermore, better estimates (in terms of MSEs) were attained in ASReml using M 2 model under the low correlation.Lowest MSEs were observed under model M 2 in ASReml and highest under model M 1 in AnimalINLA.Interestingly, lowest MSEs were attained even under the strongly negative-r uc scenario using model M 2 in ASReml.The WinBUGS software, although able to account for the specific correlation, exhibited the highest MSE of σ 2 e when the prior distribution chosen was inverse gamma (0.001, 0.001), with an analogous effect on the estimators of h 2 and c 2 .In contrast to the real data, WinBUGS estimates of the simulated data exhibited better performance when the prior utilized for σ e was the uniform distribution (MSE 44.75 vs. 215.21for the inverse gamma prior for σ 2 e ).All other parameters (σ 2 u and σ 2 c ) estimated via model M 3 in WinBUGS had relatively low MSE.
The coverage of interval estimates for the three models and the respective methods of analysis are shown in Table 5.To construct Bayesian 95 % credible intervals, the quantiles of the relevant posterior distributions (as estimated by MCMC and INLA) were used.ASReml's intervals were constructed based on asymptotic normality of the maximum likelihood using θi ± 1.96 • se( θ ), where se denotes the estimated standard error of the parameter.In the case of low r uc , the best coverages were given by ASReml and model M 2 , with narrower intervals than the Bayesian methods.In contrast, Win-BUGS exhibited the best coverage performance in the case of the high r uc , at the expense of wider intervals.Anima-  lINLA experienced difficulty in attaining nominal coverage of interval estimates when model M 1 was assumed as well as under the strongly negative-r uc scenario.In addition, DIC via WinBUGS favored the true model that incorporated the r uc in 76.67 % of the samples.

Discussion
The theoretical aspects and advantages of REML and MCMC methods for fitting hierarchical multilevel models, such as the animal model, have been extensively explored elsewhere, either with a statistical focus (Browne and Draper, 2006) or from an animal breeder's perspective (van Tassel et al., 1995;van Tassel and van Vleck, 1996).However, this is the first study applying REML and MCMC methods along with another Bayesian approach, i.e., INLA, within the context of poultry breeding.Our main concerns were the practical aspects of the applicability of three available typical software programs for the standard animal breeder.Given that both the size and the structure of data sets may have an impact on the performance of the analytical approach (Blasco, 2001), no general inference can be made based on the present results.
In the present study, an attempt to compare coverage intervals derived from Bayesian and REML approaches was pursued.However, there are two main differences between credible and confidence intervals.While a credible interval incorporates information from the prior distribution into the estimate, confidence intervals are based solely on the data, treating the parameter as fixed and the interval itself as random.Credible intervals are different from confidence intervals essentially because credible intervals are probability intervals; i.e., they say that the true value should be within the interval with a determined probability.Confidence intervals do not say that the true value is within the limits with a determined probability.In conceptual repetitions of an experiment, different confidence intervals can be obtained; 95 % of these intervals contain the true value.Thus, we treat the in-terval as containing it, knowing that, in the long run, we will be wrong 5 % of times.Although different in philosophy, the comparison between these types of intervals may be useful within the context of a study such as ours.
From a frequentist's point of view, the standard method entails the use of the REML and BLUP methods.In the present study, ASReml (Gilmour et al., 2009) software was employed.The software is stable and fast and can handle many different models, data structures and thousands of data records.In addition, the necessary files are not especially complicated to construct, while a valuable manual, containing a lot of information and numerous examples, is available for the animal breeder.For binary trait modeling, a variety of link functions (logit, probit, cloglog) can be chosen.
An obvious obstacle when using commercial programs is their limited flexibility, i.e., the inability to model complex structures between (random) effects.A good example here was the presence of negative correlation between u and c effects which could not be appropriately accommodated within the context of a typical REML package.This covariance is typically ignored (assumed to be 0), but this need not be always the case.Although in need of a more concise biological explanation, scenarios relate the negative correlation between the u and the c effects to maternally transmitted immunoglobulins, antioxidants (particularly carotenoids and vitamin E) and yolk androgens.While yolk androgens correlate positively with offspring growth (Schwabl, 1996;Groothuis et al., 2005;Müller et al., 2007), they suppress the immune system (Ketterson and Nolan, 1999;Groothuis et al., 2005) and may promote oxidative stress (von Schantz et al., 1999) in the offspring.On the other hand, maternally transmitted immunoglobulins (Buechler et al., 2002;Boulinier and Staszewski, 2008) and carotenoids (Surai and Speake, 1998) may enhance immune function, but at the expense of offspring growth.
Modeling the covariance in question was made possible only via WinBUGS.This is a very valuable feature when testing assumptions of the standard animal model with regard to possible correlation structures between the various random effects.This program allows for the application of a large group of competing models and Bayesian model evaluation criteria (Sorensen and Gianola, 2002).A further important attribute of WinBUGS is the rank tool, which can simultaneously incorporate the uncertainty associated with ranking the individuals, thus assisting in animal selection.In theory, REML and INLA would probably struggle if the likelihood was very flat, whereas MCMC methods should be able to cope (Blasco, 2001).Such scenarios could be important for practical breeding purposes and might be properly encountered by MCMC methods.Bayesian methods, such as MCMC implemented in WinBUGS, can be especially useful in complex situations at the cost of being computationally expensive and time-consuming.For our data, approximately 14 to 16 h were needed to achieve convergence, depending on modeling assumptions.
The AnimalINLA has proved to be a remarkably timeefficient experience.It took less than 10 s to produce the required posterior distributions, while providing comparable estimates with the other packages.Although computationally efficient, the current version of this R package (Anima-lINLA 1.1) could not accommodate more than 4000 records in the animal model, probably due to compatibility problems with Windows.Although time-efficient, AnimalINLA has displayed certain problems in terms of bias and accuracy, particularly for a binary trait.The latter has also been confirmed by Holand et al. (2013) and is supported by a more detailed investigation of simulated data.Finally, it is not as flexible in modeling as the WinBUGS and the documentation is still under development.
In conclusion, WinBUGS can be of great assistance to the animal breeder because of its flexibility in modeling complex models while unraveling existent data structures that the usual REML-based packages neglect.Within the animal breeding context, its applicability remains rather limited since only small to moderate data sets or populations can be handled in a time-efficient manner.Furthermore, the choice of the priors should be made with caution, particularly when the posteriors may vary with priors.The AnimalINLA software appears to be a promising future perspective for the animal breeder dedicated to the Bayesian paradigm since it is remarkably fast.It seems, however, to be a package still under development.Our own experience on large data sets has shown that ASReml can effectively handle analyses for up to 200 000 records and related pedigree structures fast (< 1 h) and mostly independent of initial values (Maniatis et al., 2013).Furthermore, as the simulation results have shown, even when a large covariance between random effects is neglected, it may provide estimates of the parameters in question with relatively small bias and error.Given all the above, ASReml remains the best practical choice for the serious animal breeder among the software packages examined.
Edited by: K. Wimmers Reviewed by: two anonymous referees

Figure 1 .
Figure 1.Distribution of ranking for 12 representative animals, based on the EBVs estimated by WinBUGS.Four animals each were taken from the top, middle and low end of the spectrum.u[i] refers to the EBV of i animal; rank 1, ..., 2456.

Table 1 .
Estimates of variance components and genetic parameters for body weight (BW) at 35 days of age.

Table 2 .
Estimates of variance components and genetic parameters for the binary transformed BW.

Table 3 .
True values and descriptive statistics of the estimators under two levels of additive genetic maternal environmental correlation.
: ratio of the maternal environmental variance to the phenotypic variance; σ uc : additive genetic maternal environmental covariance; r uc : additive genetic maternal environmental correlation; in parentheses: standard deviations; in square brackets: range[min, max].

Table 4 .
Mean squared errors of the variance components and the genetic parameters under two levels of additive genetic maternal environmental correlation.

Table 5 .
Actual coverage of nominal 95 % intervals of estimated variance components and genetic parameters.