Bayes factor between Student t and Gaussian mixed models within an animal breeding context

The implementation of Student t mixed models in animal breeding has been suggested as a useful statistical tool to effectively mute the impact of preferential treatment or other sources of outliers in field data. Nevertheless, these additional sources of variation are undeclared and we do not know whether a Student t mixed model is required or if a standard, and less parameterized, Gaussian mixed model would be sufficient to serve the intended purpose. Within this context, our aim was to develop the Bayes factor between two nested models that only differed in a bounded variable in order to easily compare a Student t and a Gaussian mixed model. It is important to highlight that the Student t density converges to a Gaussian process when degrees of freedom tend to infinity. The twomodels can then be viewed as nested models that differ in terms of degrees of freedom. The Bayes factor can be easily calculated from the output of a Markov chain Monte Carlo sampling of the complex model (Student t mixed model). The performance of this Bayes factor was tested under simulation and on a real dataset, using the deviation information criterion (DIC) as the standard reference criterion. The two statistical tools showed similar trends along the parameter space, although the Bayes factor appeared to be the more conservative. There was considerable evidence favoring the Student t mixed model for data sets simulated under Student t processes with limited degrees of freedom, and moderate advantages associated with using the Gaussian mixed model when working with datasets simulated with 50 or more degrees of freedom. For the analysis of real data (weight of Pietrain pigs at six months), both the Bayes factor and DIC slightly favored the Student t mixed model, with there being a reduced incidence of outlier individuals in this population.


(Received 2 April 2007; accepted 19 December 2007)
Abstract -The implementation of Student t mixed models in animal breeding has been suggested as a useful statistical tool to effectively mute the impact of preferential treatment or other sources of outliers in field data. Nevertheless, these additional sources of variation are undeclared and we do not know whether a Student t mixed model is required or if a standard, and less parameterized, Gaussian mixed model would be sufficient to serve the intended purpose. Within this context, our aim was to develop the Bayes factor between two nested models that only differed in a bounded variable in order to easily compare a Student t and a Gaussian mixed model. It is important to highlight that the Student t density converges to a Gaussian process when degrees of freedom tend to infinity. The two models can then be viewed as nested models that differ in terms of degrees of freedom. The Bayes factor can be easily calculated from the output of a Markov chain Monte Carlo sampling of the complex model (Student t mixed model). The performance of this Bayes factor was tested under simulation and on a real dataset, using the deviation information criterion (DIC) as the standard reference criterion. The two statistical tools showed similar trends along the parameter space, although the Bayes factor appeared to be the more conservative. There was considerable evidence favoring the Student t mixed model for data sets simulated under Student t processes with limited degrees of freedom, and moderate advantages associated with using the Gaussian mixed model when working with datasets simulated with 50 or more degrees of freedom. For the analysis of real data (weight of Pietrain pigs at six months), both the Bayes factor and DIC slightly favored the Student t mixed model, with there being a reduced incidence of outlier individuals in this population.
Bayes factor / Gaussian distribution / mixed model / Student t distribution / preferential treatment 1

. INTRODUCTION
Genetic evaluations in animal breeding are generally performed using the mixed effects models pioneered by Henderson [9]. Usually, these models assume Gaussian distributions for most random effects, including the residuals, and in absence of contradictory evidence, it is practical to assume normality on the basis of both mathematical convenience and biological plausibility. Nevertheless, departures from normality are common in animal breeding, e.g. when more valuable animals receive preferential treatment [14,15]. This preferential treatment could be defined as any management practice that increases or decreases production and is applied to one or several animals, but not to their contemporaries [14]. Amongst others, these practices may include separate housing, better (or worse) or more (or less) feed, or better (or worse) sanitary attentions. Obviously, which animals or productive records receive preferential treatment is not known with any degree of certainty in real populations and this information loss could imply substantial bias in genetic evaluations [14,15]. Other potential causes of outliers or abnormal phenotypic records could be measurement errors, sickness, shortterm-changes in herd environment and mismanagement of data [11].
We generally lack a priori sufficient information relating to the presence or absence of preferential treatment in our livestock data sets. It has been recently demonstrated that the specification of heavy-tailed residual distributions (such as the Student t distribution) instead of the usual Gaussian process in best linear unbiased prediction (BLUP) models may effectively mute the impact of residual outliers, particularly in situations where the preferential treatment of some breed stock may be anticipated [16,21]. As a result, accurate statistical tests are required to compare the mathematical simplicity of the Gaussian mixed model with the improved goodness of fit (under preferential treatment or other unknown sources of outliers) of the Student t mixed model.
General statistical tools such as the deviance information criterion (DIC) [20] or other approaches to Bayes factors [6] have been used to make comparisons between Gaussian and Student t mixed models. However, they imply high computational demands because both the Gaussian and the Student t mixed model must be analysed to calculate the corresponding comparison parameter. Within this context, the Bayes factor developed by García-Cortés et al. [5] and Varona et al. [23] in the animal breeding context implies a substantial simplification because it compares two models that only differ in terms of a single bounded variable, and therefore only the analysis of the complex model is required. The Student t distribution converges with the Gaussian distribution when the number of degrees of freedom tends to infinity. This property can be exploited 396 to appropriately adapt Varona et al. [23] Bayes factor, generating a useful statistical tool for the analysis of field data, especially when used for genetic evaluation purposes. In this paper, we focused our efforts on describing the development of this Bayes factor to make comparisons between Gaussian and Student t processes, and we tested its performance on both simulated and real data sets, using DIC as the standard reference criterion.

Statistical background for Student t mixed models
Take as a starting point a standard linear model [9] such as where y is the vector with n phenotypic data, X, W, Z are the incidence matrices of systematic (b), permanent environmental (p) and additive genetic effects (a), respectively, and e is the vector of residuals. The probability density of phenotypic data can be modeled under a multivariate Student t distribution with m degrees of freedom (with m being equal to or greater than 2): where x i , w i and z i are the ith row of X, W and Z, respectively, y i is the ith scalar element of y, r 2 e is the residual variance and C(.) is the standard gamma function with the argument as defined within parentheses. For small values of m, the Student t distribution shows a Gaussian-like pattern with increased probability in tails, whereas this distribution converges to a Gaussian distribution when m tends to infinity [16]. For mathematical convenience, we can define d = 2/m (0 d < 1) and then, the conditional density (2) reduces to a normal density when d = 0 (as is, m tends to infinity).
Following Strandén and Gianola [21], the previous model can be extended to an alternative parameterization if the data vector is partitioned according to Student t versus Gaussian mixed models m 'clusters' typified by a common factor (e.g. animal, maternal environment, herd-year-season at birth), with the previous linear model defined as: where p y j b; p; a; r 2 e ; s j À Á is a multivariate normal distribution weighted by s 2 j , p y j b; p; a; r 2 e ; s 2 I n j being an identity matrix with dimensions n j · n j , and the conditional distribution of the mixing parameter (s 2 j ) is a Gamma density with it having an expectation of 1 when d = 0 [4,21].

Bayes factor between Student t and Gaussian linear models
The Bayes factor developed by Verdinelli and Wasserman [25], and applied to the animal breeding context by García-Cortés et al. [5] and Varona et al. [23], contrasts nested linear models that only differ in terms of a bounded variable. We adapted this methodology to compare a Student t mixed linear model with its simplification to the Gaussian mixed linear model when m tends to infinity or, for mathematical convenience, d = 2/m = 0. Within this context, the posterior distribution of all the parameters of a Student t mixed model can be stated in two ways, with a pure Student t Bayesian likelihood (Model T 1 ): or with a Gaussian · Gamma Bayesian likelihood (Model T 2 ): where A is the numerator relationship matrix between individuals. Following in part Varona et al. [23], the prior distribution assumed for the bounded variable (d) was assumed The permanent environmental and the additive genetic effects were assumed to be drawn from multivariate normal distributions, with I p being an identity matrix with dimensions equal to the number of elements of p. The prior distributions for the remaining parameters of the model were defined as: Student t versus Gaussian mixed models where k 1 , k 2 , k 3 and k 4 are four values that were small enough to ensure a flat distribution over the parameter space [23].
The joint posterior distribution of all the parameters in the alternative Gaussian mixed model (Model G) was proportional to where the Bayesian likelihood was defined as multivariate normal, and the prior distributions were identical to the prior distributions of Model T 1 (or Model T 2 ). The Bayes factor between Model T 1 (or Model T 2 ) and Model G (BF T/G ) can be easily calculated from the Markov chain Monte Carlo sampler output of the complex model (Student t mixed model). Under Model T 1 , the conditional posterior distribution of all the parameters in the model did not reduce to wellknown distributions and generic sampling processes such as Metropolis-Hastings [8] are required. Simplicity was gained under the alternative Model T 2 during the sampling process. In this case, sampling from all the parameters in Model T 2 can be performed using a Gibbs sampler [7], with the exception of d, which requires a Metropolis-Hastings step [8]. Following García-Cortés et al. [5] and Varona et al. [23], the posterior density p T (d = 0|y) suffices to obtain BF T/G , because p T (d = 0) = 1 (see equation (9)). Alternatively, The BF T/G can be obtained by averaging the full conditional densities of each cycle at d = 0 using the Rao-Blackwell argument [26]. At this point, computational simplicity is gained with Model T 1 (or a normal density for d = 0), whereas Model T 2 tends to computationally unquantifiable extreme probabilities when d is close to zero. A BF T/G greater than 1 indicates that the Student t mixed model is more suitable, whereas a BF T/G smaller than 1 indicates that the Gaussian mixed model is more suitable.
From the standard definition of the Bayes factor [13], where PO T/G is the posterior odds between models, PrO T/G is the prior odds between models, and p T and p G are the a priori probabilities for Student t mixed model and Gaussian mixed model, respectively. In the standard development of the Bayes factor described above, we assumed that prior odds were 1 and p T and p G were both 0.5. Nevertheless, we could modify prior odds depending on our a priori knowledge, e.g. Student t mixed model is a more parameterized model and it could be easily penalized with a smaller-than-1 prior odds. Posterior odds can be viewed as the weighted value of the Bayes factor, conditional to our a priori degree of belief.

Simulation studies
The Bayes factor methodology developed above was validated through simulation. Seven different scenarios were analyzed following a Student t residual process, with degrees of freedom equal to 5 (d = 0.4), 10 (d = 0.2), 20 (d = 0.1), 50 (d = 0.04), 100 (d = 0.02), 200 (d = 0.01) and 300 (d = 0.007), respectively. Twenty-five replicates were simulated for each case and each replicate included five non-overlapping generations with 200 individuals (10 sires and 190 dams) and random mating. Following Model T 2 , each individual had a phenotypic record and was assigned its own independent cluster. Data were generated from a normal density N ðl; Ir Ã e Þ weighted by a cluster-characteristic value drawn from equation (6). Note that l included a unique systematic effect (10 levels randomly assigned with equal probability and sampled from a uniform distribution between 0 and 1) and a normally distributed additive genetic effect generated Student t versus Gaussian mixed models under standard rules [1]. Residual and additive genetic variances were equal to 1 and 0.5, respectively. This simulation process generated seven different scenarios with 25 data sets which were analyzed twice, through the previously described Bayes factor and through a standard Gaussian model (Model G). For each analysis, a single chain was launched that contained 100 000 rounds, after discarding the first 10 000 rounds as burn-in [19]. Comparisons between the two models were performed through three approaches: (a) Bayes factor between nested models, (b) DIC [20], and (c) correlation coefficient between simulated and predicted breeding values (q a,ã ). Note that DIC is based on the posterior distribution of the deviance statistic [20], which is À2 times the sampling distribution of the data as specified in formula (2) or as the conjugated distribution of (5) and (6)

Analysis of weight at six months in Pietrain pigs
After editing, 2330 records of live weight at six months in Pietrain pigs were analyzed, with an average weight (± SE) of 102.9 (± 0.265) kg. These pigs were randomly chosen from 641 litters from successive generations grouped in 135 batches during the fattening period, and their records were collected between years 2003 and 2006 in a purebred Pietrain farm registered in the reference Spanish Databank (BDporcÒ, http://www.bdporc.irta.es). At the beginning of the fattening period (two months of age), batches were created with pigs from different litters in order to homogenize piglet weight, and these groups were maintained up to slaughter (six months of age). Pigs were reared under standard farm management during the suckling and fattening periods. Pedigree expanded up to five generations and comprised 2601 individuals, with 109 boars and 337 dams with known progeny.
The operational model included the additive genetic effect of each individual, the permanent environmental effect characterized by the batch during the fattening period, and three systematic sources of variation: sex (male or female), year · season with 11 levels, and age at weighing (180.0 ± 0.3 days) treated as a covariate. Data were analyzed by applying the Bayes factor described above and assuming a different cluster for each pig with phenotypic data. To easily 402 compare this method with a standard Gaussian model, data were also analyzed under Model G. The empirical correlation between estimated breeding values (posterior mean) was calculated in the two models and, as for the simulated data sets, DIC was calculated for Model T and Model G. Each Gibbs sampler ran with a single chain of 450 000 rounds after discarding the first 50 000 iterations as burn-in [19].

Simulated datasets
Summarized results of the 25 replicates for each simulated Student t process (5, 10, 20, 50, 100, 200 and 300 degrees of freedom) are shown in Table I. Estimates for additive genetic variance showed coherent behavior with average estimates slightly greater than 0.5. Average residual variance estimated using the Student t mixed model clearly agreed with the simulated value. Nevertheless, residual variance was clearly over-estimated for simulations with few degrees of freedom in which a Gaussian mixed model was applied, showing higher standard errors in data sets with few degrees of freedom. Simulations with 5 degrees of freedom showed the highest average residual variance under the Gaussian mixed model (1.664 ± 0.038), whereas the average residual variance was reduced to 1.222 ± 0.025 for replicates with 10 degrees of freedom, and converged to one for datasets with 300 degrees of freedom (showing a standard error smaller than 0.020). Under the Student t mixed model, average estimates of degrees of freedom fitted with true values without any noticeable bias, although precision decreased with larger degrees of freedom (Tab. I). Substantial discrepancies were observed between the two models in terms of predicted breeding values in extreme heavy-tailed simulations. Although the correlation coefficients between predicted breeding values in the Student t and Gaussian mixed models increased quickly in line with the degrees of freedom, the empirical correlation in replicates with 5 degrees of freedom was very small (0.377 ± 0.030) and average correlations greater than 0.9 were observed in simulations with 100 or more degrees of freedom (Tab. I).
Empirical correlations between simulated and predicted breeding values increased with degrees of freedom in both the Student t and Gaussian mixed models, although the Student t mixed model reached higher correlations when simulated degrees of freedom were small. As seen in Table II, simulations under extremely heavy-tailed processes (5 degrees of freedom) showed average correlations of 0.420 and 0.377 for Student t and Gaussian mixed models, respectively, suggesting substantial bias for genetic evaluations performed with Student t versus Gaussian mixed models

404
standard Gaussian models when normality did not hold. Differences between the two models quickly decreased with increasing degrees of freedom and were almost negligible for m = 20 and higher values. As was expected, both the Bayes factor and DIC clearly favored Student t mixed models in simulated scenarios with few degrees of freedom and showed similar behavior throughout the analyzed framework (Fig. 1A). Datasets simulated under a residual Student t process with 5 degrees of freedom reached an average Bayes factor favoring the Student t mixed model of 4.6 · 10 92 , and the average difference between DIC was also huge (À 490). Although both comparison criteria decreased when degrees of freedom increased, our results suggest that the Student t mixed model was preferable instead of the Gaussian model up to 20 degrees of freedom, and that even for simulations with 50 degrees of freedom, the superiority of the Student t model remained almost total (Tab. III). Substantial discrepancies between the Bayes factor and DIC appeared in the last two scenarios (200 and 300 degrees of freedom). While the average Bayes factor slightly favored the Gaussian model, the DIC continued to produce smaller estimates for the Student t model, although with only a minimal difference in the last scenario (Tab. II). This suggests that the Bayes factor was more conservative, favoring the less parameterized model. This hypothesis was confirmed in Table III where the Bayes factor supported the Gaussian mixed model in 0, 0, 0, 1, 11, 18 and 19 data sets (for simulations with 5, 10, 20,    406 50, 100, 200 and 300 degrees of freedom), whereas DIC favored the Gaussian mixed model in 0, 0, 0, 0, 3, 5 and 10 data sets, respectively (Fig. 1B).

Analysis of Pietrain pig weight at six months
Analysis under a Student t linear mixed model placed the highest posterior density region at 95% (HPD95) for the d parameter between 0.004 and 0.009, with the modal value at 0.005 (Tab. IV). Degrees of freedom (m) determine  (Fig. 2). The posterior mean of the weights (s 2 i ) for the Student t mixed model ranged from 0.924 to 1.026, although  The Bayes factor favored the Student t model rather than the Gaussian model (BF T/G = 2.532), although the small size of this value is not worth more than a bare mention according to Jeffreys' [12] scale of evidence. He/she suggested that differences could be minimal, although DIC reported substantial discrepancies between the Student t mixed model (16 445) and the Gaussian mixed model (16 450); this difference was smaller than the average difference between DIC in simulations with 200 degrees of freedom. On the contrary, the empirical correlation coefficient between predicted breeding values in each model was 0.999 (this was 0.993 if only breeding animals were considered) and showed that, although the posterior probabilities of both models were slightly (Bayes factor) or substantially (DIC) different, genetic evaluations performed with a Student t or a Gaussian mixed model provided an almost identical genetic ranking for this data set.

DISCUSSION
The Bayes factor as originally proposed by Verdinelli and Wasserman [25] has been applied to various models used in animal breeding and the genetics research field. Although the method was initially developed to test for the genetic background of linear traits [5] and the location of quantitative trait loci (QTL) [23], this Bayes factor has been recently modified to discriminate between linked and pleiotropic QTL [24], to test for the genetic background of threshold traits [2,18], and to compare different structures of random genetic groups [3], all with encouraging results. Previous research highlighted its computational stability [5]. Although two nested models were compared, only analysis of the complex model was required, which allowed a substantial reduction in computational requirements. Within this context, these properties highlight an interesting Bayesian tool that could be extended to compare a wide range of models. Analysis of simulated data sets showed correct fitting of residual variances under a Student t mixed model. On the contrary, Gaussian mixed models suffered from substantial over-estimations of residual variance when phenotypic records were drawn from a Student t process that involved relatively few degrees of freedom. It confirmed previous research [11] and this phenomenon could be directly related to heavy-tailed characteristics of Student t density. Nevertheless, it is important to highlight the fact that in each simulation scenario, the average estimate of additive genetic variance across 25 replicates was placed close to the true value (Tab. I) in both the Student t and Gaussian mixed models. A slight overestimation was suggested for average additive genetic variance, but this is commonly observed in Bayesian animal models in which modal estimates of r 2 a are unbiased, whereas the mean tends towards slight or moderate overestimation [22]. This stability in the estimation of genetic variance contrasted with the substantial discrepancies observed between the two models in terms of empirical correlation between predicted breeding values (between the two models), as well as empirical correlations between simulated and predicted breeding values (within a given model). For small degrees of freedom, the Student t mixed model reached higher correlations and verified the results previously reported by Strandén and Gianola [21] and Jamrozik et al. [11], who suggested that Gaussian models had only a limited ability to mute the effects of outliers in genetic evaluations. As expected, the Bayes factor confirmed the superiority of the Student t mixed models when data were generated from a Student t process with few degrees of freedom, and its performance was comparable to that provided by DIC (Fig. 1A). Our study shows that, although differences between Student t and Gaussian mixed models may be small when degrees of freedom tend towards large values, Student t mixed models could be preferable within a wide range of degrees of freedom. This confirmed previous results reported by Jamrozik et al. [11]. Indeed, the Bayes factor provided substantial evidence favoring the Student t model for six data sets in the replicate with 300 degrees of freedom. Slight discrepancies appeared between the Bayes factor and DIC in simulations with 50 or more degrees of freedom, where the former 410 could be viewed as a more conservative method. As suggested in Figure 1B and Table II, both the Bayes factor and DIC favored the Gaussian mixed model in 16 cases, whereas they produced contradictory results in a further 20 cases. This does not invalidate our methodology, although it suggests that for large degrees of freedom, the more parameterized model (Student t) would tend to be more penalized using the Bayes factor.
The analysis of Pietrain pig weight at six months provided an example of the Bayes factor applied to real data. Estimated variance components and their ratios were similar for both the Student t and the Gaussian mixed models. This was not surprising given the mode of degrees of freedom in the Student t mixed model (191.12; Tab. IV). Heritability for weight at six months was around 0.25, which was similar to the values reported by Hovernier et al. [10] and Noguera et al. [17]. There are no comparable results in the animal breeding literature for Student t mixed models with high values of degrees of freedom, whereas Jamrozik et al. [11] and Chang et al. [4] reported extremely heavy-tailed distributions (2.4 and 8.5 degrees of freedom, respectively) for mastitis in dairy cattle. This suggests a clearly higher incidence of unknown sources of outliers in dairy cattle than in fattening pigs, as well as a greater probability of the influence of preferential treatment. As seen in Table V, the weight parameters (s 2 j ) of the Student t mixed model were placed around 1, showing less dispersion than in other research involving a smaller estimate of degrees of freedom [4,11]. Nevertheless, some outlier individuals have been observed (e.g. 20 individuals with s 2 j lower than 0.95) and this result confirms the slight advantage of the Student t model over the Gaussian model suggested by the Bayes factor and DIC. Moreover, whereas the Bayes factor showed a small value that was close to 1, DIC showed substantial differences close to 5 units. This confirmed previous results obtained under simulation, in which the Bayes factor was suggested as being the more conservative. Nevertheless, given the extreme correlation obtained between the predicted breeding values in the two models and the computational simplicity provided by Model G, it would seem preferable to use a standard Gaussian mixed model for genetic evaluations. The high correlation between breeding values predicted by the two models confirmed the results obtained under simulation. Nevertheless, these results are conditional to our data set and they can not be taken as a general rule for genetic evaluations in swine.
In conclusion, a straightforward Bayes factor was developed to compare Gaussian and Student t mixed linear models, proposing two alternative parameterizations of the Student t model. This methodology could be of special interest to check preferential treatment or other sources of outlier data, which is a common phenomena in livestock data. Student t versus Gaussian mixed models