Attenuating effects of preferential treatment with Student-t mixed linear models: a simulation study

On a simule le traitement preferentiel de certaines vaches dans quatre troupeaux de selection utilisant la transplantation embryonnaire. La frequence et l'effet du traitement preferentiel ont dependu d'une fonction correlee a la valeur genetique vraie. On a compare trois modeles lineaires mixtes pour leur aptitude a prendre en compte le traitement preferentiel : le modele classique Gaussien, un modele avec des erreurs t-multivariates groupees par troupeau et un modele avec des erreurs t-distribuees independantes. Dans le modele ou les erreurs suivaient une distribution t, les parametres d'echelle et les degres de liberte ont ete consideres inconnus. Une analyse bayesienne a ete effectuee pour les trois modeles a partir de l'echantillonnage de Gibbs et les moyennes a posteriori ont ete utilisees pour en inferer au sujet de la variance genetique, des effets troupeau-annee, des valeurs genetiques et des reponses realisees a la selection. La performance des modeles a ete evaluee au travers des erreurs quadratiques moyennes. En l'absence de traitement preferentiel, les trois modeles ont eu une performance similaire. Quand le traitement preferentiel a ete frequent et d'effet important, le modele t-univariate a ete le meilleur et le modele Gaussien a ete clairement inadapte. Il apparait que des modeles lineaires robustes peuvent prendre en compte les traitements preferentiels mieux que les modeles lineaires mixtes Gaussiens classiques.


INTRODUCTION
Preferential treatment is any management practice that is applied nonrandomly to animals within a contemporary group [9]. For example, better housing and feeding, hormonal treatment, longer milking intervals on test day and feeding according to production are known to be applied selectively in dairy production. Preferential treatment occurs in dairy cattle, presumably to increase the economic value of a cow or the probability that it will be chosen as a bull-dam. Several studies (e.g. [17, 20!) have found that genetic evaluations for milk yield are inconsistent with expectations based on theory. This may be due to inadequate statistical assumptions or failure to account properly for selection or preferential treatment of cows. Preferential treatment is often suspected when no apparent reasons exist for such discrepancies. Kuhn et al. [9] simulated effects of preferential treatment on 'animal model' genetic evaluations. Mean squared error of prediction of breeding values increased as the extent of preferential treatment increased. Kuhn and Freeman [10] found that when the dam of a sire was treated preferentially, more than 30 daughters with untreated records were needed to offset the bias in prediction of breeding value caused by the dam's information. Bias increased as the proportion and number of daughters receiving preferential treatment increased. Bias decreased when all daughters given preferential treatment were in the same herd; this is so because the 'herd-year' effect in the model captures part of the preferential treatment administered in a particular herd-year.
In order to account for preferential treatment, Harbers et al. [7] included an environmental correlation between related females in a genetic evaluation model for a MOET (multiple ovulation and embryo transfer) scheme. This improved accuracy of cow evaluations when preferential treatment was mild. Weigel et al. [29] simulated different strategies of preferential treatment and found that it was not possible to detect it by monitoring within-herd variance; obviously, this parameter does not provide information about the probability that a cow within a herd is treated preferentially. Burnside and Meyer [3] simulated effects of bovine somatotropin (bST). Sire evaluations were least accurate when bST administration was targeted to the best producing cows.
In the context of prediction (e.g. !8]), a bias takes place when the expected values of the predictand and of the predi!tor differ. Evaluation of bias requires knowledge of the true model but, in practice, this is not available, so ad hoc assessments of bias have been suggested. Several studies [15,16,27,28] found upward 'biases' of cow's pedigree indexes for protein or milk yield in Finnish Ayrshire. It is unclear if this discrepancy is due to chance, but preferential treatment of dams of cows may be a culprit. On the other hand, Powell and Norman [19] found that pedigree indexes understated the first estimated breeding values of daughters of proven sires mated to lower producing dams. Little work has been undertaken on how to cope with preferential treatment in practice, at least from a statistical point of view. Kuhn and Freeman [11] studied power transformations of records but this was, at best, slightly effective in reducing bias due to preferential treatment. An alternative approach is to consider an error distribution with thicker tails than the normal, to allow for more variation. A commonly used one is the t-distribution, which is symmetric and leptokurtic. It has been advocated because of its simplicity [12], and because only one parameter (the degrees of freedom) is needed to describe robustness. A suitable robust distribution may be capable of attenuating the impact of outliers on data analysis. Many authors have employed statistical models with t-distributed residuals [4,12,13,25,31] in linear and non-linear regression models, with varying degrees of success. Use of the t-distribution in the context of mixed effects or hierarchical models is relatively recent !1, 2, 5, 6, 22-24, 26, 30].
Our objective was to assess frequentist properties of Bayesian point estimators obtained from mixed linear models where residuals were assumed to be either Gaussian or t-distributed. Milk production records obtained in herds in which some preferential treatment was practised were simulated. The analysis focused on mean squared error of estimation of genetic variance, herd-year effects, breeding values and genetic response to selection.

Conceptual population
Milk production records in a hypothetical 'adult' MOET nucleus scheme [18] were simulated. The scheme extended the simple hierarchical mating structure of Stranden et al. !21]. Our modification allowed bulls of the previous generation to mate current generation females. The nucleus consisted of 32 cows and eight bulls in every generation. In each generation, every nucleus cow produced (by multiple ovulation and embryo transfer to recipients) eight offspring, four females and four males. An animal could be selected only once into the nucleus as a parent and unselected animals were culled. The females were selected among those offspring to the nucleus that had completed a first lactation. Males were selected within those that had been born in the preceding generation. In practice, this would allow the bulls to have a progeny test outside the nucleus before selection. However, such progeny testing was not built in this simulation.
Thus, males within a full-sib family had the same estimated breeding value and three such males were randomly discarded. Each selected male was mated to four cows, chosen randomly from those that had been selected as replacements. per generation. With this scheme carried out for four generations, the data included 544 cows with records (32 in the base plus 32 x 4 x 4 = 512 female progeny) and 32 sires with daughters in production, i.e. a total of 576 animals.
A diagram of the simulated population is shown in figure 1.
Base generation cows were assigned to four herds in equal numbers, i.e. eight cows per herd. Female offspring of a cow remained in the same herd as her dam, whereas sires were used across herds. Breeding values of base animals were drawn at random from N(0,0.25) distribution. Records of the base animals were generated by adding a herd-year effect (independently, normally distributed) to a breeding value and to an independently drawn residual from N(0,0.75) distribution. Records in subsequent generations were simulated similarly, except that the breeding value of an individual was formed by averaging the breeding value of its parents and adding a N C 0 , 2 1L7 u 2(l F )J segregation residual, where a 2 is the additive genetic variance and F is the average inbreeding coefficient of the parents. The selection criterion in the breeding scheme was BLUP of breeding value with the true variance components. The statistical model included the herd-year as a fixed effect and animal as a random effect (but ignored preferential treatment, as discussed later) using all genetic relationships available up to the time of selection.

Preferential treatment
In practice, preferential treatment takes place in the course of a selection programme so this is the way that the present simulation proceeded. None of the base population cows were treated preferentially, so there were 512 cows eligible to receive preferential treatment. A scheme in which the preferential treatment assigned depended on the 'perceived' breeding value of an animal (e.g. based on a genetic evaluation available before the animal produces the record) was adopted. The records were generated as where y ij is the record of animal j made in herd-year i, h i is a herd-year effect, Uj is the breeding value of animal j, and e ij is an independent residual. The preferential treatment 02! was a stochastic effect taking the values: where !(.) is the standard normal cumulative distribution function, p min is a constant smaller than the herd-year effect h i , and u/j = !+(u!+v!) / afl + w is a 'value' function such that Ui rv A!(0,er!), v j -N(0, Q v), Cov(u_,,f j ) = 0, so Wj rv N(À, 1). In the preceding, 0' ; is the variance of breeding values and w j2 is an 'uncertainty' variance. The ratio O'! 2 describes the uncertainty the herd !u manager has about the true breeding value of animal j. For example, if the breeder is very uncertain about the breeding value of the animal, this ratio of variances should be high. Three values of the uncertainty were considered, The constant p min in equation (2) controls the range of production associated with preferential treatment. It was set equal to -5 Qh where a h is the standard deviation of herd-year effects. These were drawn from a normal distribution with mean zero and variance a 2 and two different values of the herd-year variance were considered: afl = au and 2 U 2 The constant A controls the proportion of cows to be preferentially treated. Normal distribution theory can be used to find a value of A such that a desired proportion of cows receives preferential treatment. The proportion of preferentially treated cows increases with .!, because Pr(u/j > 0) increases concomitantly. Three different prevalences of preferential treatment were considered: 1 out of 10, 1 out of 32, and 1 out of 64 cows. These correspond to A values of -1.2816, -1.8627 and -2.1539, respectively. It was intended to keep the proportion of preferentially treated animals roughly constant from generation to generation. To do so, it must be noted that selection is expected to increase mean breeding value and to reduce genetic variance over time. In order to account for these effects, the formula for w was changed to: where u is the mean breeding value of animals available for preferential treatment in the generation to which animal j belongs, and Su is the additive genetic variance for individuals born in that generation.
The probability distribution of the amount of preferential treatment (Di!) depends on the values of Qh and A as shown in the Appendix. The average amount of preferential treatment actually applied was assessed via a simulation of 1000 replicates of the MOET scheme. Mean increase (mean of 0) in production due to preferential treatment under varying prevalence of preferential treatment and amount of herd-year variance is in table I. As intended, production increased with prevalence of preferential treatment, and with or h* 2 Average value of preferential treatment was not affected by level of uncertainty j 2 2 . This is not shown in table 1, but it was expected because the distribution Q! of A zj does not depend on this ratio. Table I. Average increase in simulated lactation production due to preferential treatment as a function of herd-year variance (a h ) and of prevalence of preferential treatment (values in parenthesis are Monte Carlo standard errors from 1 000 replicates of the MOET scheme, au = additive genetic variance).

Statistical models and computations
Three linear statistical models were compared, both with and without preferential treatment incorporated in the simulation. The objective was to assess the relative ability of these models to handle perturbations caused by unknown preferential treatment. In all three models, the linear structure for the records included an unknown herd-year effect (treated as fixed computationally), the unknown breeding value of the cow and a residual, distributed according to an appropriate error distribution, as noted below. In the three models, a multivariate normal distribution N(0, Aa!), where A is a 576 x 576 relationship matrix, was used for the genetic effects, so there was no difference in this respect. The three models, differing only in the error distribution were the following. 1) G: a purely Gaussian model with errors Niid(0, Q e ).
2) t-l: errors were independently and identically distributed as univariate-t, t i (0, 0' e 2, v,). Here, the variance of the distribution is U2V ,I(V, e -2), where Q e is a scale parameter and v, are the unknown degrees of freedom.
3) t-H: within herd i (i = 1, 2, 3, 4), the error vector e i had the multivariatet distribution t ni (0, I ni U2 , v e ) where n i is the number of records in herd i. Here, Var(e i ) = I!o!fe/(fe -2). Although the errors are uncorrelated, they are not independent, this being a property of the multivariate t-distribution. Error vectors in different herds were mutually independent, however, but with the same or and v e parameters. We refer to this model as a 'herd-clustered' one.
The G model is the usual one; model t-1 discounts outliers y zj on a 'case' by 'case' basis, and model t-H discounts outlying vectors y z for the entire herd i. Because the 'value function' w! used to generate preferential treatment does not depend on the herd, there is no apparent reason why model t-H should outperform model t-1. It should be noted that as Ve -j oo, the two t-distributions tend towards the Gaussian one.
A Bayesian structure was adopted for inference. Prior distributions were the same for all three models. Herd effects were assigned a uniform prior and, as noted, a multivariate normal process was used as a prior distribution for the breeding values. The dispersion components Q u and Q e were assigned independent scaled inverted chi-square distributions with four degrees of freedom and mean equal to the true variance component, i.e. 0.75 for the residual variance and 0.25 for the genetic variance. In the t-models, the prior for a is for the scale of the distribution and not for the residual variance, which is a e 2v ,l(v, -2) as noted before. In the two models involving the t-distribution, the residual degrees of freedom parameter v, was considered unknown. Degrees of freedom values allowed in the herd-clustered t-model were 4, 10, 100 or 1 000, all equally likely, a priori. In the univariate t-distribution model, the space of v e was 4, 6, 8, 10, 12 or 14, all receiving equal prior probability. These values were chosen arbitrarily. It is possible to use a continuous prior for v, [23] but the discrete distribution employed here facilitated implementation. A Gibbs sampler was used to carry out the Bayesian computations employing the full conditional distributions described in Strandén [22]. Tests made in several simulations with varying starting values indicated that a burn-in period of 7 000 iterates with 70 000 Gibbs iterates thereafter (all samples kept) was enough to obtain sufficiently precise estimates of posterior means of the parameters. About 60 min of CPU time were required to perform 70 000 iterations, for any of the models, in an HP 9 000(3) computer.

Frequentist comparison
Each replicate of the simulation consisted of a data set generated as per the scheme in figure 1 under the appropriate assumptions of preferential treatment. A Bayesian analysis of the data set according to each of the three models was carried out in each replicate. Mean squared errors of posterior mean estimates were computed, over replicates, for: a) genetic variance, b) herd-year effects, and c) breeding values. Mean squared errors were also computed for three classes of breeding values: sires, cows who had been preferentially treated and cows without preferential treatment. d) An additional end-point of interest was mean squared error of estimated response to selection, assessed by predicting breeding values using posterior means from the three models contrasted. 'True' response was the mean difference in true breeding value (due to selection using BLUP) between animals born in the last generation and those born in the first generation. Differences in mean squared errors between models should reflect the relative accuracy of estimation of genetic trend.
A 'pilot run' [14] was conducted to assess the number of replicates needed to attain enough precision for a parameter of interest. The approximate number of replications required to achieve an absolute precision r for the confidence interval given a pilot run of n replicates was found using: where t i -1 , 1 -a/2 is the value of a t-distribution with i -1 degrees of freedom at the 100(1 -a) percentile ('confidence'). Our pilot study consisted of carrying n = 20 replicates for each of the three models. The number of replications required to achieve 0.05 precision with 95 % confidence for the genetic variance was less than 60 for most cases. Hence, it was decided that all cases would be replicated 60 times. Absolute precision was recalculated after 60 replicates, and a further 40 replicates were made for the schemes involving 1/10 prevalence of preferential treatment. One scheme 92 2 = 3 , -& m d a s h ; 2 = 100 ) required an ( a2 a2 1) u u 100 / additional 40 replicates to achieve the required precision. Table II indicates the schemes and number of replicates performed.. Because of its heavy computing requirements, the analysis was performed using a network of machines administered by Professor Miron Livny of the Department of Computer Science, University of Wisconsin at Madison. This cluster was accessed using the Condor system, which allows running jobs simultaneously at many computers while the data and program reside in one computer. Each replicate of each model was a process to be executed in this network of computers. There were between 10 and 15 computers available at any time, giving at least a 10-fold increase in computing power compared to using only the HP9000(3).

Absence of preferential treatment
The objective here was to examine possible losses in efficiency due to using the two t-distribution models when there is no preferential treatment and the Gaussian assumption holds throughout. Averages and mean squared errors of estimates of additive genetic variance are given in table III. The posterior means of afl for each of the three models were practically unbiased, in light of the Monte Carlo variation. However, the mean squared error was larger for the two t-models than for the Gaussian one. Hence, if the Gaussian assumption holds, posterior means of additive genetic variance for the t-models are less accurate than those from the G-model. The increase in mean squared error over the Gaussian model was about 5-6 % for the t-H model, and 7-18 % for the t-1 model.
Tables IV and Vgive the posterior distributions of the degrees of freedom for the two t-models in the absence of preferential treatment. The analysis carried out with the herd-clustered t-model clearly favoured a model with Gaussian errors, as indicated by a posterior probability of about 90 % for the degrees of freedom being larger than 10. Also, the univariate t-model assigned the highest posterior probability, about 40 %, to the largest value of the degrees of freedom (v e = 14) considered. The posterior distributions were not sharp, this being a function of the low informational content the data have about v e . However, both analyses favoured the larger values of v e or, equivalently, the Gaussian assumption for the errors. For example, in the herd-clustered t-model, the posterior odds ratio of v, = 1000 relative to v, = 4 was 17.7 and 29:7 for In summary, in the absence of preferential treatment and with the Gaussian assumption holding throughout, the t-models were less accurate for estimation of Q!, but were as competitive as the Gaussian model for estimation of breeding values and of genetic trend. 3 figure 2. Differences between models were clearest when preferential treatment was more prevalent (1/10) and when the herd-year variance was high (this affects the distribution of 0). Also, differences between models were largest when uncertainty about true breeding values was low, so the value function is a high correlate of breeding value. There was little difference between the G and the t-H models, but the univariate t-model had the best performance when prevalence of preferential treatment was medium (1 out of 32 cows) or high (1 out of 10 cows). The univariate t-model was robust to variation in the uncertainty parameter; this was not the case for the G and the t-H models, whose performance was hampered under severe forms of preferential treatment.

Posterior distribution of the degrees of freedom
Posterior probabilities of the degrees of freedom under the herd-clustered t-model were often higher for the larger values of v e , thus supporting the Gaussian model, especially when preferential treatment was uncommon, or uncertainty was high. Only under medium (1/32) or high (1/10) prevalence of preferential treatment and a high herd-year variance the largest values of the degrees of freedom did not have the highest posterior probability. However, this depended on the level of uncertainty and on the amount of herd-year variance. For example, when prevalence of preferential treatment was 1/10 and with ( T2 = 3a', low values of the degrees of freedom had higher posterior probabilities when uncertainty was low; however, as uncertainty increased, the posterior distribution tended to favour larger values of the degrees of freedom. Posterior probabilities of v, for the univariate t-model are given in figure 3. Here, posterior distributions tended to be flat. Higher probabilities were assigned to the largest values of the degrees of freedom only when preferential treatment was rare and the herd-year variance low. As in the herd-clustered t-model, high uncertainty often resulted in higher probabilities assigned to the highest degrees of freedom values, as one would expect. However, other degrees of freedom values also received relatively high probabilities. When preferential treatment was prevalent (1/10) and the herd-year variance was large, the posterior distribution was sharp, with a modal value of v e = 4 at all levels of uncertainty. This points away from a Gaussian distribution of the residuals. With a small data set such as the one in this simulated MOET scheme, one should not expect the posterior distribution of the degrees of freedom parameter to be highly peaked. Nevertheless, the univariate t-model recognised the non-Gaussian situation even when prevalence was rare (1/64), provided that the variance between herds was relatively large. This is because the expected value of the preferential treatment, E(Di!), increased with or2 h, as illustrated in Average of mean squared error of estimates of herd-year effects was similar for the three models except when preferential treatment was prevalent or herdyear variance was high, but it was always smallest for the univariate t-model (figure l!). When preferential treatment was common (1/10), the univariate t-model clearly had the smallest mean squared error at each level of uncertainty and value of o,2 It * The average of mean squared error of estimates of all breeding values is shown in figure 5. This criterion was about the same with all models except when preferential treatment was common and the herd-year variance high. Here, when uncertainty was high, there were no differences between the models, but at low levels of uncertainty, the univariate t-model was markedly superior. The picture for mean squared errors of estimates of sire and cow breeding values and for genetic response was similar to that for of all breeding values, so the figures are not presented. In all cases, differences between models were clear, favouring the univariate t-model when preferential treatment was more prevalent (1/10). The same was true for preferentially treated cows (figure 6), but mean squared errors were larger than for breeding values of cows that were not treated preferentially. The univariate t-model had a similar or slightly worse performance than the Gaussian or herd-clustered models when preferential treatment was rare or mildly prevalent, but it was superior when such treatment was common. In particular, at the lowest level of uncertainty and at the highest herd-year variance, the univariate t-model gave predictions of breeding value of preferentially treated cows that had a mean squared error of about a third of that observed with the Gaussian model. In this situation, the herd-clustered . model improved estimates somewhat relative to the Gaussian model.

CONCLUSIONS
In the absence of preferential treatment, the t-models were as good as the Gaussian model for estimating breeding values and response to selection. When preferential treatment was mildly prevalent (1/32) the models performed similarly. However, when preferential treatment was common (1/10) and especially when the herd-year variance was large relative to the additive genetic variance, the univariate t-model was clearly the best, at least in terms of mean squared error. Under preferential treatment, the posterior distribution of the degrees of freedom in the univariate t-model pointed away from the correctness of the Gaussian assumption. The univariate t-model was quite robust to variation in the simulation parameters, but it is unknown whether this robustness holds across different forms of preferential treatment.
This simulation could not differentiate clearly between the Gaussian and the herd-clustered t-models, although the latter was always slightly better under preferential treatment. A reason for the lack of difference between these two models may be the low number of herds in the simulation. With a few clusters (herds) the statistical information about the degrees of freedom is low, so the posterior distribution of this parameter cannot be estimated accurately.
In conclusion, it appears that the univariate t-model can attenuate adverse effects of preferential treatment as applied here. It leads to better inferences about breeding values and genetic trends than those obtained with the Gaussian model, especially when preferential treatment is prevalent, at least under the conditions of the study. If, on the other hand, preferential treatment is nonexistent, or the assumption of a Gaussian distribution of the residuals seems to be true, there is little loss in efficiency from using a robust model, such as the univariate t. It is encouraging that a symmetric error distribution, such as Student t, improved upon the Gaussian one under a single-tailed form of preferential treatment as in equation (2). This suggests that a robust asymmetric distribution may do even better, but perhaps at the expense of conceptual and computational simplicity.