Mixed effects linear models with t-distributions for quantitative genetic analysis: a Bayesian approach

Original article Mixed effects linear models with t-distributions for quantitative genetic analysis: a Bayesian approach Ismo Strandén Daniel Gianola a Department of Animal Sciences, University of Wisconsin, Madison, WI 53706, USA b Animal Production Research, Agricultural Research Centre - MTT, 31600 Jokioinen, Finland (Received 21 July 1998; accepted 27 November 1998) Abstract - A Bayesian approach for inferences about parameters of mixed effects linear models with t-distributions is presented, with emphasis on quantitative genetic applications. The implementation is via the Gibbs sampler. Data from a simulated multiple ovulation and embryo transfer scheme in dairy cattle breeding with nonrandom preferential treatment of some cows is used to illustrate the procedures. Extensions of the model are discussed. © Inra/Elsevier, Paris mixed effects models / Bayesian inference / robust estimation / Gibbs sampling / Student’s t-distribution Résumé - Modèles linéaires mixtes avec distributions de Student en génétique quantitative : approche bayésienne. On présente une approche bayésienne en vue de l’inférence concernant les paramètres de modèles linéaires mixtes avec des distributions de Student, en mettant l’accent sur les applications en génétique quantitative. L’application s’effectue grâce à l’échantillonnage de Gibbs. Des données provenant d’un schéma de sélection simulé utilisant le transfert embryonnaire chez les bovins laitiers en présence d’un traitement préférentiel de quelques vaches sont utilisées pour illustrer les procédures. Les extensions du modèle sont discutées. © Inra/Elsevier, Paris modèle mixte / inférence bayésienne / estimation robuste / échantillonnage de Gibbs / distribution t de Student * Correspondence and reprints E-mail: [email protected] 1. INTRODUCTION Mixed effects linear models are used widely in animal and plant breeding and in evolutionary genetics [27]. Their application to animal breeding was pioneered by Henderson [17, 19-21], primarily from the point of view of making inferences about candidates for genetic selection by best linear unbiased prediction (BLUP). Because BLUP relies on knowledge of the dispersion structure, estimation of variance and covariance components is central in practical implementation [14, 18, 29, 32]. Typically, the dispersion structure is estimated using a likelihood-based method and, then, inferences proceed as if these estimates were the true values (e.g. [8]). Although normality is not required by BLUP, it is precisely when normality holds that it can be viewed as an approximation to the best predictor [4, 8, 12, 19]. More recently, Bayesian methods have been advocated for the analysis of quantitative genetic data with mixed linear models [8, 9, 34, 39, 40], and the Bayesian solutions suggested employ Gaussian sampling models as well as normal priors for the random effects. It is of practical interest, therefore, to study statistical models that are less sensitive than Gaussian ones to departures from assumptions. For example, it is known in dairy cattle breeding that more valuable cows receive preferential treatment, and to the extent that such treatment cannot be accommodated in the model, this leads to bias in the prediction of breeding values [23, 24]. Another source of bias in inferences is an incorrect specification of the inheritance mechanism in the model. It is often postulated that the genotypic value for a quantitative trait is the result of the additive action of alleles at a practically infinite number of unlinked loci and, thus, normality results [4]. This assumption is refuted in an obvious manner when inbreeding depression is observed, or when unknown genes of major effect are segregating. However, in the absence of clearly contradictory evidence, normality is a practical assumption to make, as then the machinery of mixed effects linear models can be exploited. An appealing alternative is to fit linear models with robust distributions for the errors and for the random effects. One of such distributions is Student’s t, both in its univariate and multivariate forms. Several authors [2, 7, 26, 37, 38, 41, 42] have studied linear and non-linear regression problems with Student’s t-distributions, but there is a scarcity of literature on random effects models. West [41] described a one-way random effects layout with t-distributed errors and a heavy tailed prior for the random effects. Assuming that the ratio between residual variance and the variance of the random effects was known, he showed that this model could discount effects of outliers on inferences. Pinheiro et al. [30] described a robust version of the Gaussian mixed effects model of Laird and Ware [25] and used maximum likelihood. They hypothesized that the distribution of the residuals had the same degrees of freedom as that of the random effects, and, also, that random effects were independently distributed. The first assumption is unrealistic as it is hard to accept why two different random processes (the distributions of random effects and of the residuals) should be governed by the same degrees of freedom parameter. The second assumption is not tenable in genetics because random genetic effects of relatives may be correlated. In quantitative genetics the random effects or functions thereof are of central interest. For example, in animal breeding programs the objective is to increase a linear or non-linear merit function of genetic values which, ideally, takes into account the economics of production [16, 28, 33]. Here, it would seem natural to consider the conditional distribution of the random effects given the data, to draw inferences. There are two difficulties with this suggestion. First, it is not always possible to construct this conditional distribution. For example, if the random effects and the errors have independent t-distributions, the conditional distribution of interest is unknown. Second, this conditional distribution would not incorporate the uncertainty about the parameters, a well-known problem in animal breeding, which does not have a simple frequentist or likelihood-based solution (e.g. [10, 15]). If, on the other hand, the parameters (the fixed effects and the variance components) are of primary interest, the method of maximum likelihood has some important drawbacks. Inferences are valid asymptotically only, under regularity conditions, and finite sample results for mixed effects models are not available, which is particularly true for a model with t-distributions. In addition, some genetic models impose constraints such that the parameter space depends on the parameters themselves, so it would be naive to apply a regular asymptotic theory. For example, with a paternal half-sib family structure [6], the variance between families is bounded between 0 and one-third of the variance within families. Moreover, maximum likelihood estimation in the multi-parameter case has the notorious deficiency of not accounting well for nuisance parameters [3, 8, 13]. A Bayesian approach for drawing inferences about fixed and random effects, and about variance components of mixed linear models with t-distributed random and residual terms is described here. Section 2 presents the probability model, emphasizing a structure suitable for analysis of quantitative genetic data. Section 3 gives a Markov chain Monte Carlo implementation. A Bayesian analysis of a simulated animal breeding data set is presented in section 4. Potential applications and suggestions for additional research are in the concluding section of the paper. 2. THE UNIVARIATE MIXED EFFECTS LINEAR MODEL 2.1. Sampling model and likelihood function Consider the univariate linear model where y is an n x 1 vector of observations; X is a known, full rank, incidence matrix of order n x p for ’fixed’ effects; b is a p x 1 vector of unknown ’fixed’ effects; Z is a known incidence matrix of order n x q for additive genetic effects; u is a q x 1 vector of unknown additive genetic effects (random) and e is an n x 1 vector of random residual effects. Although only a single set of random effects is considered, the model and subsequent results can be extended in a straightforward manner. It is assumed that u and e are distributed independently. Suppose the data vector can be partitioned according to ’clusters’ induced by a common factor, such as herd or herd-year season of calving in a cattle breeding context. The model can then be presented as: where m is the number of ’clusters’ (e.g. herds). Here yi is the data vector for cluster i (i = 1, 2, ... , m), Xi and Zi and are the corresponding incidence matrices and ei is the residual vector pertaining to yi. Observations in each cluster will be modeled using a multivariate tdistribution such that, given b and u, data in the same herd are uncorrelated but not independent, whereas records in different clusters are (conditionally) independent. Let yi !b, u, 62 N t ni (X ib + Zi u, 1,,, o, e 2, ve ), where ni is the number of observations in cluster i (i = 1, 2, ... , m), or’ is a scale parameter and v. is the degrees of freedom. If ni = 1 for all i, the sampling model becomes univariate t. The conditional density of all observations, given the parameters, is Although the m distributions have the same v, and Qe parameters, these are not identical. In particular, note that E(y2 !b, u, Qe, ve) Xib + Zi u, and Var(y2!b, u, Qe, ve) = hzv!./(v! - 2), i = 1,2,...,m, so the mean vector is peculiar to each cluster. Homoscedasticity is assumed, but this restriction can be lifted without difficulty. When each cluster contains a single observation, the error distribution is the independent t-model of Lange et al. [26]; then, the observations are conditionally independent. When all observations are put in a single cluster, the multivariate t-model of Zellner [42] results; in this case, the degrees of freedom cannot be estimated. Each of the m terms in equation (3) can be obtained from the mixture of the normal distribution: with the mixing process being: where xv e is a chi-squared random variable on Ve degrees of freedom [26, 38, 41,42]. 2.2. Bayesian structure Formally, both b and u are location parameters of the conditional distribution in equation (3). The distinction between ’fixed’ and ’random’ is frequentist, but from a Bayesian perspective it corresponds to a situation where there is a differential amount of prior information on b and u [8, 13]. In particular, the Bayesian counterpart of a ’fixed’ effect is obtained by assigning a flat prior to b, so that the prior density of this vector would be: in Rp. This distribution is improper, but lower and upper limits can be assigned to each of the elements of b, as in Sorensen et al. [34], to make it proper. The prior distribution of additive genetic values u will be taken to be a multivariate t-distribution, and independent of that of b. From a quantitative genetics point of view this can be interpreted as an additive, multivariate normal model (as in [4]), but with a randomly varying additive genetic variance. Because the multivariate t-distribution has thicker tails than the normal, the proposed model is expected to be somewhat buffered against departures from the assumptions made in an additive genetic effects model, so ’genetic outliers’ stemming from nonadditivity or from major genes become, perhaps, less influential in the overall analysis. All properties of the multivariate normal distribution are preserved, e.g. any vector or scalar valued linear combination of additive genetic values has a multivariate t-distribution, the marginal distributions of all terms in u are t, and all conditional distributions are t as well. In particular, if the additive genetic values of parents and the segregation residual of an offspring are jointly distributed as multivariate t, the additive genetic value of the offspring has a univariate t-distribution with the same degrees of freedom. This implies that the coancestry properties of the usual Gaussian model are preserved. We then take as prior distribution: with density Above, q is the number of individuals included in u (some of which may not have data), A is a known matrix of additive relationships, au is a scale parameter and vu is the degrees of freedom parameter. Hence, Var(ulo, 2, v!) = A<r!tt/(ftt !2), which reduces to the variance-covariance matrix of additive genetic values of a Gaussian model when vu -! oo. The scale parameters or2 and or2 are taken to have independent scaled inverted chi-square distributions, with densities: respectively, for af > 0 and or > 0. Here, T’e ( Tu ) is a strictly positive ’degree of belief’ parameter, and Te (T u) can be thought of as a prior value of the scale parameter. These distributions have finite means and variances whenever the T parameters are larger than 2 and 4, respectively. In animal breeding research, it is common practice to assign improper flat priors to the variance components of a Gaussian linear model [8, 9, 39]. If uniform priors are to be used, it is advisable to restrict the range of values they can take, to avoid impropriety (often difficult to recognize, see [22]). Here, one can take Typically, the lower bounds are set to zero, whereas the upper bounds can be elicited from mechanistic considerations, or set up arbitrarily. Prior distributions for the degrees of freedom can be discrete as in Albert and Chib [1] and Besag et al. [2], or continuous as in Geweke [7], with the joint prior density taken as p(ve,v!) = p(v e )p(v u ). In the discrete setting, let fj, j = 1, 2, ... , d e, and wk, k = 1, 2, ... , d!, be sets of states for the residual and genetic values degrees of freedom, respectively. The independent prior distributions are: Because a multivariate t-distribution is assigned to the whole vector u, there is no information contained in the data about v,,. Therefore, equation (11) is recovered in the posterior analysis. There are at least two possibilities here: 1) to assign arbitrary values to Vu and examine how variation in these values affects inferences, or 2) to create clusters of genetic values by, e.g. half-sib or full-sib families, and then assume that clusters are mutually independent but with common degrees of freedom. Here the v. parameter would be estimable, but at the expense of ignoring genetic relationships other than those from half-sib or full-sib structures. Alternative (2) may be suitable for dairy cattle breeding (where most of the relationships are due to sires) or humans (where most families are nuclear). A third alternative would be to use (2), then find the mode of the posterior distribution of vu, and then use (1) as if this mode were the true value. In the following derivation, we adopt option (1). The joint prior density of all unknowns is then: with obvious modifications if equations (8) and (9) are used instead of equations (6) and (7). The joint posterior density is found by combining likelihood equation (3) and appropriate priors in equations (4)-(11), to obtain: where b E !p, u E K,j,o! > 0, Q! > 0 and ve E f f j , j = 1, 2, ... , ,de} if a discrete prior is employed. The hyper-parameters are 7e , TM! Te , T u and vu because we assume this last one to be known. Hereafter, we suppress the dependency on the hyper-parameters in the notation. 3. THE GIBBS SAMPLING SCHEME A Markov chain Monte Carlo method such as Gibbs sampling is facilitated using an augmented posterior distribution that results from mixture models. The t-distribution within each cluster in equation (3) is viewed as stemming from the mixture processes noted earlier. Likewise, the t-distribution in equation (5) can be arrived at by mixing the ulA, or 2,s2 - N(O, AU2 /82 ) process with s2 wu N X2. Iv,,. The augmented joint posterior density is m where s, = (se l , ... , sP m ) and N = ! n2. Integration of equation (14) with ! i=i respect to Se and s2 yields equation (13), so these posteriors are ’equivalent’. There is a connection here with the heterogeneous variance models for animal breeding given, e.g. in Gianola et al. [11] and in San Cristobal et al. [31]. These authors partitioned breeding values and residuals into clusters as well, each cluster having a specific variance that varied at random according to a scale inverted chi-square distribution with known parameters. The full conditional distributions required to instrument a Gibbs sampler are derived from equation (14). Results given in Wang et al. [40] are used. Denote C = !c2!!, i, j = 1, 2,...,p + q, , and r = {r j, i,j = 1, 2,...,p + q to be the coefficient matrix and right-hand side of Henderson’s mixed model equations, respectively, where p + is the number of unknowns (fixed and random effects), given the dispersion components Se , sfl and the scale parameters Qe and or 2 u The mixed model equations are: best linear unbiased estimator (BLUE) of b, and u is the best linear unbiased predictor (BLUP) of u. Collect the fixed and random effects into a’ = (b’, u’) _ (a,, a2 , ... , ap + q) . Let a’ i = (a l , a 2 , ... , ai-1, a i+l , ... , ap + q). The conditional posterior distribution of each of the elements of a is / P+9 B where 4i = c1 r! - y! c,j’aj L i, j = 1, 2, ... , p + q. This extends to blocks B j-1 / B j$i / of elements of a in a natural way. If ai is a sub-vector of a, the conditional distribution of ai given everything else, is multivariate normal with mean ai = Ci il ri - E Cija i for appropriate definitions of C ij , ri and a! as B j-1 / B j54, / matrices and vectors. The conditional posterior density of each of the se is in the form of a gamma density where s, _, is s, without S; i’ Equivalently, , &dquo;e e / Similarly, the conditional posterior density of s! also has the gamma density form The conditional posterior distribution of Qe is a scaled inverted chi-square distribution with form If a bounded uniform distribution is used as prior for Q e, its conditional posterior is the truncated distribution: The conditional posterior density of ou is: distributions, than those obtained with a mixed effects Gaussian linear model, the current paradigm in quantitative genetics [19, 21, 27]. Our approach was illustrated with simulated data from a dairy cattle breeding scheme, where cows were subject to fairly prevalent and strong preferential treatment. A univariate t-model for the errors led to more accurate inferences about additive genetic variance than either a herd-clustered t-model or a Gaussian sampling process. The posterior distributions of breeding values of some example animals were sharper in the univariate t-model. Our model and implementation can be extended in several respects. For example, if the degrees of freedom of the distribution of genetic values needs to be assessed, it is possible to cluster the genetic values into ’independent’ families and proceed as for the residual variance. However, such clustering would lead to a loss of accuracy in the specification of the genetic variance-covariance structure, because relationships between individuals in different clusters would not be taken into account. Another extension would be to take the degrees of freedom as continuous and use a rejection algorithm or a MetropolisHastings walk to draw samples, combined with the Gibbs sampler for the rest of the parameters of the model. Additional random effects, such as permanent environmental effects affecting all records of a cow, can be incorporated, e.g. by taking a univariate t-distribution as prior. Residuals can be clustered in different manners. For example, clustering errors by sire or full-sib families may cope with inadequate genetic assumptions, e.g. unknown major genes may be segregating. In addition, it is possible to allow for heterogeneous variance in the model without major difficulty. Other residual distributions such as the logistic or the slash may be considered as well. At present, it is not yet possible to apply these methods to the large data sets used for routine genetic evaluation in the dairy cattle breeding industry, where the models can have millions of individual breeding values. Hence, if it is established that models based on the t-distribution improve genetic evaluations, computationally simpler or faster methods should be developed. One possibility would be to employ Laplacian approximations to assess the mode of the joint distribution of the dispersion parameters and then use some form of conditional analysis to obtain point predictors of breeding values. In a model with tdistributed random effects, modal estimates of fixed and random effects, given the scale parameters and the degrees of freedom can be found using an iterative procedure [35]. This requires solving reweighted mixed model equations several times, as in threshold models. Approximate solutions after a couple of iterations may be adequate for practical purposes. This Laplacian-iteration approach would be counterpart to the standard REML-BLUP analysis in linear models under Gaussian assumptions. Finally, we would like to observe some interpretative differences between the Gaussian and the t-model from a quantitative genetic point of view. Heritability is the regression of genotype on phenotype, which is o’!/(<!+o’!) in a Gaussian model. In the t-models discussed, this regression is with the numerator (and the appropriate part of the denominator) being equal to o,2if genetic values are assumed to be Gaussian, instead of t-distributed, as was the case in our analysis. The posterior distribution of heritability can be estimated by forming a sample value for h2 from the corresponding draws of ol 2, v!, or2and ve. For our simulation, posterior mean estimates of heritability were 0.47 and 0.45 for the Gaussian and herd-clustered t-models, respectively, and 0.19 for the univariate t-model, this being closest to the input value. REFERENCES [1] Albert J.H., Chib S., Bayesian analysis of binary and polychotomous response data, J. Am. Stat. Assoc. 88 (1993) 669-679. [2] Besag J., Green P., Higdon D., Mengersen K., Bayesian computation and stochastic systems, Stat. Sci. 10 (1995) 3-66. [3] Box G.E.P., Tiao G.C., Bayesian Inference in Statistical Analysis, Wiley, New York, 1973. [4] Bulmer M.G., The Mathematical Theory of Quantitative Genetics, Clarendon Press, Oxford, 1980. [5] Casella G., George E.I., Explaining the Gibbs sampler, Am. Stat. 46 (1992) 167-174. [6] Falconer D.S., McKay T.F.C., An Introduction to Quantitative Genetics, 4th ed., Longman, New York, 1996. [7] Geweke J., Bayesian treatment of the independent Student-t linear model, J. Appl. Econometrics 8 (1993) S19-S40. [8] Gianola D., Fernando R.L., Bayesian methods in animal breeding theory, J. Anim. Sci. 63 (1986) 217-244. [9] Gianola D., Foulley J.L., Variance estimation from integrated likelihoods, Genet. Sel. Evol. 22 (1990) 403-417. [10] Gianola D., Foulley J.L., Fernando R.L., Prediction of breeding values when variances are not known, Genet. Sel. Evol. 18 (1986) 485-498. [11] Gianola D., Foulley J.L., Fernando R.L., Henderson C.R., Weigel K.A., Estimation of heterogeneous variances using empirical Bayes methods: theoretical considerations, J. Dairy Sci. 75 (1992) 2805-2823. [12] Goffinet B., Selection on selected records, Genet. Sel. Evol. 15 (1993) 91-97. [13] Harville D.A., Bayesian inference for variance components using only error contrasts, Biometrika 61 (1974) 383-385. [14] Harville D.A., Maximum likelihood approaches to variance component estimation and to related problems, J. Am. Stat. Assoc. 72 (1977) 320-340. [15] Harville D.A., BLUP (Best Linear Unbiased Prediction) and beyond, in: Gianola D., Hammond K. (Eds.), Advances in Statistical Methods for Genetic Improvement of Livestock, Springer-Verlag, Berlin, 1990, pp. 239-276. [16] Hazel L.N., The genetic basis for constructing selection indexes, Genetics 28 (1943) 476-490. [17] Henderson C.R., Specific and general combining ability, in: Gowan J.W. (Ed.), Heterosis, Iowa State College Press, Ames, IA, 1950, pp. 352-370. [18] Henderson C.R., Estimation of variance components and covariance components, Biometrics 9 (1953) 226-252. [19] Henderson C.R., Sire evaluation and genetic trends, in: Proceedings of the Animal Breeding and Genetics Symposium in honor of Dr J.L. Lush, Blacksburg, VA, August 1973, American Society of Animal Science, Champaign, IL, 1973, pp. 10-41. [20] Henderson C.R., Best linear unbiased estimation and prediction under a selection model, Biometrics 31 (1975) 423-447. [21] Henderson C.R., Applications of Linear Models in Animal Breeding, University of Guelph Press, Guelph, 1984. [22] Hobert J.P., Casella G., The effect of improper priors on Gibbs sampling in hierarchical linear models, J. Am. Stat. Assoc. 91 (1996) 1461-1473. [23] Kuhn M.T., Freeman A.E., Biases in predicted transmitting abilities of sires when daughters receive preferential treatment, J. Dairy Sci. 78 (1995) 2067-2072. [24] Kuhn M.T., Boettcher P.J., Freeman A.E., Potential biases in predicted transmitting abilities of females from preferential treatment, J. Dairy Sci. 77 (1994) 2428-2437. [25] Laird N.M., Ware J.H., Random effects models for longitudinal data, Biometrics 38 (1982) 963-974. [26] Lange K.L., Little R.J.A, Taylor J.M.G., Robust statistical modeling using the t distribution, J. Am. Stat. Assoc. 84 (1989) 881-896. [27] Lynch M., Walsh B., Genetics and Analysis of Quantitative Traits, Sinauer Associates, Sunderland, MA, 1997. [28] Meuwissen T.H.E., Goddard M., Selection of farm animals for non-linear traits and profits, Anim. Sci. 65 (1997) 1-8. [29] Patterson H.D., Thompson R., Recovery of interblock information when block sizes are unequal, Biometrika 58 (1971) 545-554. [30] Pinheiro J.C., Liu C., Wu Y., Robust estimation in linear mixed effects models using the multivariate-t distribution, Bell Labs Technical Report, 1997, 04/97. [31] San Cristobal M., Foulley J.L., Manfredi E., Inference about multiplicative heteroscedastic components of variance in a mixed linear Gaussian model with an application to beef cattle breeding, Genet. Sel. Evol. 25 (1993) 3-30. [32] Searle S.R., Casella G., McCulloch C., Variance Components, John Wiley, New York, 1992. [33] Smith H.F., A discriminant function for plant selection, Ann. Eugenics 7 (1936) 240-250. [34] Sorensen D.A., Wang C.S., Jensen J., Gianola D., Bayesian analysis of genetic change due to selection using Gibbs sampling, Genet. Sel. Evol. 26 (1994) 333360. [35] Stranden I., Robust mixed effects linear models with t-distributions and application to dairy cattle breeding, Ph.D. thesis, University of Wisconsin, Madison, WI, 1996. [36] Stranden L, Gianola D., Attenuating effects of preferential treatment with Student-t mixed linear models: a simulation study, Genet. Sel. Evol. 30 (1998) 565583. [37] Sutradhar B.C., Ali M.M., Estimation of the parameters of a regression model with a multivariate t error variable, Comm. Stat. Theory Meth. 15 (1986) 429-450. [38] Verdinelli I., Wasserman L., Bayesian analysis of outlier problems using the Gibbs sampler, Stat. Comput. 1 (1991) 105-117. [39] Wang C.S., Rutledge J.J., Gianola D., Marginal inferences about variance components in a mixed linear model using Gibbs sampling, Genet. Sel. Evol. 25 (1993) 41-62. [40] Wang C.S., Rutledge J.J., Gianola D., Bayesian analysis of mixed linear models via Gibbs sampling with an application to litter size in Iberian pigs, Genet. Sel. Evol. 26 (1994) 91-115. [41] West M., Outlier models and prior distributions in Bayesian linear regression, J. R. Stat. Soc. B Met. 46 (1984) 431-439. [42] Zellner A., Bayesian and non-Bayesian analysis of the regression model with multivariate Student-t error terms, J. Am. Stat. Assoc. 71 (1976) 400-405.