On Marker-Assisted Prediction of Genetic Value: Beyond the Ridge

Corresponding author: Daniel Gianola, 1675 Observatory Dr., Madison, WI 53706., gianola{at}calshp.cals.wisc.edu (E-mail)

Communicating editor: J. B. WALSH

	ABSTRACT

TOP ABSTRACT A MIXED-EFFECTS MODEL... EXTENDING THE HIERARCHY A BAYESIAN FORMULATION DISCUSSION LITERATURE CITED

Marked-assisted genetic improvement of agricultural species exploits statistical dependencies in the joint distribution of marker genotypes and quantitative traits. An issue is how molecular (e.g., dense marker maps) and phenotypic information (e.g., some measure of yield in plants) is to be used for predicting the genetic value of candidates for selection. Multiple regression, selection index techniques, best linear unbiased prediction, and ridge regression of phenotypes on marker genotypes have been suggested, as well as more elaborate methods. Here, phenotype-marker associations are modeled hierarchically via multilevel models including chromosomal effects, a spatial covariance of marked effects within chromosomes, background genetic variability, and family heterogeneity. Lorenz curves and Gini coefficients are suggested for assessing the inequality of the contribution of different marked effects to genetic variability. Classical and Bayesian methods are presented. The Bayesian approach includes a Markov chain Monte Carlo implementation. The generality and flexibility of the Bayesian method is illustrated when a Lorenz curve is to be inferred.

THE availability of a plethora of markers has led to consideration of the issue of the extent to which molecular information can be used to advantage in genetic improvement programs of agricultural species, such as maize or dairy cattle. There is a large body of literature on this matter (e.g., SOLLER and BECKMANN 1983 ; SMITH and SIMPSON 1986 ; LANDE and THOMPSON 1990 ); WHITTAKER 2001 gives a review. Marker-assisted selection can be effective as long as part of the genetic variance can be associated with segregating marker loci (LANDE and THOMPSON 1990 ).

The basic idea in marker-assisted selection is to exploit statistical dependencies (linkage disequilibrium) existing in the joint distribution of marker and quantitative trait loci (QTL) genotypes. For example, when two inbred lines are crossed, the disequilibrium is manifest in the F₂ generation. On the other hand, when there is linkage equilibrium at the population level, only the joint distribution of marker and QTL genotypes within a family is nontrivial (OLLIVIER 1998 ). GEORGES et al. 1995 have exploited high levels of within-family disequilibria in dairy cattle populations for QTL mapping. Most livestock populations have some disequilibrium due at least to chance (small effective size), as shown by FARNIR et al. 2000 .

Linkage disequilibrium between markers and QTL can be used for two main purposes: (1) to infer genomic location and effects of QTL affecting a trait and (2) to arrive at improved (in some statistical sense) predictions of genetic merit of candidates for selection in a breeding program. These two objectives may not be disjoint (e.g., FERNANDO and GROSSMAN 1989 ). However, genetic cartography of QTL is not a requirement for prediction of genetic merit or marker-assisted selection (LANGE and WHITTAKER 2001 ). In fact, remarkable advances have been made in prediction of breeding values in livestock since the introduction of the best linear unbiased predictor (BLUP; HENDERSON et al. 1959 ; HENDERSON 1973 ). This method uses only phenotypic and pedigree information, with the QTL viewed only in an abstract manner. On the other hand, it has been argued and shown, at least in simulations, that molecular information may enhance the accuracy of selection (e.g., WHITTAKER et al. 2000 ).

Our concern is with statistical models and methods for inferring genetic merit using molecular and phenotypic information. The objective is to describe phenotype-marker associations using multilevel hierarchical linear models. The setting is mainly as in WHITTAKER et al. 2000 and LANGE and WHITTAKER 2001 , i.e., situations where conditional independence of genetic sampling can be assumed, such as in an F₂ population derived from a cross between inbred lines. Model features include chromosome-specific effects, spatial associations of markers within chromosomes, existence of background genetic variability, and heterogeneity among families, if some such clustering exists. Classical and Bayesian methods are described. First, we present a mixed-effects model formulation and a BLUP implementation, with the dispersion components estimated by some likelihood-based procedure. An extension of the mixed-effects model is given subsequently. Finally, a Bayesian formulation is presented, including a Markov chain Monte Carlo procedure for drawing samples from target posterior distributions. Possible applications to outbred populations are discussed.

	A MIXED-EFFECTS MODEL FORMULATION

TOP ABSTRACT A MIXED-EFFECTS MODEL... EXTENDING THE HIERARCHY A BAYESIAN FORMULATION DISCUSSION LITERATURE CITED

Hierarchical representation:
Let the phenotypic value of individual i for a quantitative trait in an F₂ from a cross between inbred lines be described by the model

(1)

Here, the p x 1 vector ß contains some systematic effects representing, e.g., year of harvest, level of fertilization, or plant density; x^'_i is a known incidence vector relating ß to y_i; a_i is an unobserved genetic value; and e_i is an independently distributed random residual reflecting environmental variability or inadequacy of the model. It is assumed that the genetic value a_i results from an unknown number (K, say) of QTL acting additively, so that

where Q_ik is the genotype at biallelic (assumed for simplicity) locus k for individual i and _k is the per-allele effect. If a_i is random, (1) is a special case of the well-known mixed linear model (e.g., HENDERSON 1973 ). When K goes to infinity this becomes the classical infinitesimal model of quantitative genetics.

It is conceptually convenient to develop model (1) hierarchically. The first level of a Gaussian hierarchy is given by

(2)

where N (·) indicates a normal distribution and ²_e is the environmental variance. If the n environmental deviates are independently and identically distributed, this leads to the matrix representation

(3)

where X is an n x p incidence matrix assumed (without loss of generality) to have full-column rank, a = {a_i}, and I denotes an identity matrix, in this case n x n. Suppose next that individual i has been typed for marker genotypes at each of l loci; this is represented by the vector

with

Assume that all individuals have been typed for all markers (this is not realistic, but see the DISCUSSION). Then, the unobserved genetic value can be modeled as

(4)

where = {_k} is of order l x 1. We refer to _k as the "marked effect" of marker locus k on genetic value, while _i can be interpreted as some residual or "background" genetic effect not involved in the association between genetic value and the markers but, yet, having an effect on phenotype. The vector is the gradient or regression of the unobservable additive genetic value on the observable marker genotype, that is, = a_i/m_i. As noted by LANDE and THOMPSON 1990 , dominance can be introduced by expanding (4) as, e.g.,

where m^2'_i is a row vector with elements consisting of the squares of the corresponding entries of m_i. Interactions between marked effects for different marker loci can be modeled via cross-products between appropriate elements of m_i. For simplicity, additivity of marked effects is assumed throughout.

The second level of the hierarchy is represented by a distribution describing the uncertainty about genetic values, given the marked effects, that is the background genetic variability. We adopt the Gaussian model,

where ² is the background additive genetic variance. It is assumed (rightly or wrongly, depending on the context), given the marker genotypes and , that the "background" genetic effects _i of different individuals are mutually independent. This implies that either there is no family structure or that, conditionally on the marked effects , the family structure is not relevant. Family clustering is taken up in a later section. In matrix notation, and consistently with (4), the assumption of independence leads to

(5)

where M is the n x l matrix of known marker genotypes. Unless there is some prior knowledge about ²_e and ² (or some clustering of individuals, such as a family structure), the background effect _i must be lumped together with e_i, because of nonidentifiability. On the other hand, if the variances ² and ²_e are known a priori, it is possible to "predict" _i distinctly from e_i, in the same way that one can predict additive genetic and environmental effects via BLUP when dispersion parameters are known.

WHITTAKER et al. 2000 and WHITTAKER 2001 treat as a fixed parameter and employ ridge regression for estimation of this vector. From a Bayesian perspective (e.g., LINDLEY and SMITH 1972 ; ZELLNER and VANDAELE 1975 ), this is equivalent to regarding as having the distribution

(6)

with ² elicited in some manner. Assumption (6) is adopted as the third level of the hierarchy. In a frequentist setting, the assumption states that the marked effects are drawn at random from the multivariate normal distribution (6) in each conceptual repetition of a crossing experiment. In a Bayesian setting, this would be part of the prior ensemble of the model; this is discussed later. Regression coefficients for markers that are not adjacent to QTL are expected to be null under the assumption of no interference (ZENG 1993 ). However, without knowing the location of the QTL in relation to the markers (the usual situation), it is not obvious how such a prior consideration can be incorporated into the model. Since (6) is a prior distribution in the Bayesian sense, its influence on inferences can be tempered by measuring enough individuals. At this point, it suffices to say that can be treated as a random effect merely as a device for obtaining possibly improved (in some sense) predictions of genetic value. HAYES and GODDARD 2001 assumed that unobservable gene or chromosome effects followed a Gamma distribution, so these effects would be strictly positive even though estimated values can be negative. MEUWISSEN et al. 2001 used Gamma deviates for simulating gene effects, but then "tossed a coin" to determine their sign. Therefore, it is unclear what would be gained from assuming a Gamma distribution for the elements of . Note that (6) implies that the "marked effects" are independent and identically distributed, but this assumption is relaxed later on. The normality assumption in (6) is probably adequate enough and facilitates computation significantly.

The three-stage hierarchy can be condensed by inserting (5) into (3), so that the model describing the phenotypic values can be written as

(7)

(8)

where = {_i} and e = {e_i}. This can be viewed as a frequentist mixed-effects model, where ß is a fixed location parameter and and are random terms. Unless additional assumptions are made or some knowledge about the partition of variance in the population is available, _i and e_i are "confounded." However, it is conceptually useful to maintain these two vectors as distinct. The marginal (frequentist) distribution of the phenotypes induced by model (8) is the normal process

(9)

where MM'² is the variance-covariance matrix of marked genotypic values, conditionally on the marker genotypes M observed in the experiment. In scalar notation, the "total variance" of a single observation is

Hence,

(10)

is interpretable as the fraction of genetic variance attributable to marked effects, in an experiment repeated over and over with the marker genotypes fixed across replications. The variance "due to" the association with the markers is ^l_k=1m²_ik², which depends nontrivially on m_i, the specific marker genotype of individual i. On the other hand, since the marker genotypes vary at random over replications,

(11)

where the expectation and covariance matrix are taken over the distribution of marker genotypes in the population. Knowledge of the distribution of marker genotypes is needed for evaluation of (11).

Best prediction and best linear prediction:
Under standard assumptions, with ß and the variance components ², ², and ²_e known, the joint distribution of and , given the phenotypes and the marker genotypes, is the multivariate normal process:

(12)

Here,

(13)

where = ²_e/², = ²_e/², and

(14)

The best linear predictor [BLP; also the best predictor (BP) under normality; HENDERSON 1973] of the unobserved total genetic values a = M + is

(15)

and the variance-covariance matrix of the prediction error (under normality this is also the covariance matrix of the conditional distribution of a) is

(16)

The best predictor has the smallest possible mean-squared error of prediction. Hence, it would be difficult to improve upon this, provided that the model is reasonable, normality holds, and parameters are known. Thus, (13–16) provide an alternative to a ridge regression approach to prediction. Generalization to multiple traits measured in different individuals is straightforward, but this is not dealt with here. Since, given M, the BLP or BP is unbiased (e.g., HENDERSON 1973 ), it follows automatically that it is also unbiased unconditionally. However, unconditionally,

with the expectation taken with respect to the joint distribution of marker genotypes in the entire population. A drawback of BP or BLP is that it is unrealistic to assume that ß, ², ², and ²_e are known without error.

Best linear unbiased prediction:
An obvious improvement is to use BLUP. BLUP takes into account uncertainty about ß, which is not the case of BP or BLP above, where ß is treated as known. Under normality, BLUP(a) can be interpreted as the mean of the conditional distribution of the predictand a = M + , given a vector of "error contrasts," denoted as w. For example, take w = y - X, where is either the ordinary least-squares or the generalized least-squares estimator of ß. In such a setting, BLUP is the best predictor under normality, but only in the class of linear translation invariant predictors (SEARLE 1974 ; GIANOLA and GOFFINET 1982 ). It is well known that for (8) and (9),

(17)

where, for P = I - X(X'X)^-1X',

(18)

gives = BLUP() and = BLUP(). Further,

(19)

The BLUP of a = M + is

and the variance-covariance matrix of the prediction errors is

As noted, BLUP of the unobservable genetic values (marked or background genetic effects or any linear combination thereof) takes into account the uncertainty about ß, although given , , and ²_e. Since BLUP is conditionally unbiased, it follows that it is also so unconditionally (averaged over marker genotypes). The unconditional covariance matrix of the prediction errors is the average of (19) taken over the distribution of marker genotypes.

The needed variance components:
Sensible values of the dispersion parameters must be specified for implementing BLUP. As stated, unless there is prior or external knowledge, it is not possible to separate the variance of the background genetic effects (²) from that of the environmental influences (²_e). The model can be rewritten as

(20)

where e* N(0, I^2*_e). Methods for obtaining maximum-likelihood estimates of ² and of (² + ²_e) for (20) are well known, e.g., SEARLE et al. 1992 . In this setting only the marked part of the genetic value can be predicted. Here, the conditional distribution used for inferring would be

(21)

where

is the best predictor (given the error contrasts) or BLUP(). The covariance matrix of the conditional distribution (21) is

with this being the same as the covariance matrix of the prediction errors - , given the marker genotypes M. The BLUP of the marked quantitative trait genotype is M, and the prediction error dispersion matrix is MC^*M'^2*_e. Hence, M could be used as a criterion for genetic evaluation in marker-assisted selection ignoring background effects. If reasonable maximum-likelihood (ML) or restricted maximum-likelihood (REML) estimates of ² and of ^2*_e are available, these can be treated as "true" values in (21).

Suppose that estimates of the additive genetic (²_a) and of the environmental variance (²_e) are available (or of heritability ² and the phenotypic variance ²_y) from an analysis ignoring marker information. Then, the variance of the background genetic effects (assuming additivity) could be estimated as ² = ^2*_e - ²_e = ^2*_e - (1 - ²)²_y (hoping that this value will be positive). Then, form

where ² is the ML or REML estimate of ², and proceed with the BLUP implementation in (17). The resulting predictor is an empirical BLUP having properties that depend largely on the accuracy of the variance component estimates and on the adequacy of model in (7) and (9).

Differences with ridge regression:
A main difference with the ridge regression approach of WHITTAKER et al. 2000 and WHITTAKER 2001 is how the value of (or of ^* = ^2*_e/²) is obtained. Here, the formalism of the random-effects model (with a justification from a Bayesian viewpoint) is favored over an ad hoc procedure for doing shrinkage and tempering colinearity in a fixed-effects model, such as ridge regression. LINDLEY and SMITH 1972 , ZELLNER and VANDAELE 1975 , and BOX 1980 discuss ridge regression from a Bayesian perspective. Further, a random-effects treatment of a fixed-effects model can lead to estimators and predictors with a superior frequentist performance (e.g., ZELLNER and VANDAELE 1975 ; GIANOLA 1990 ; WEIGEL et al. 1991 ).

An advantage of a random-effects treatment is the flexibility to accommodate additional model features. For example, the ridge regression estimator of WHITTAKER 2001 implies that all marked effects in are independent (a priori). Evidence of coexpression of genes in at least the same chromosome (CARON et al. 2001 ) indicates that the assumption that marked sections or segments within a chromosome have independent effects is not always tenable. Adjacent QTL may coexpress more (that is, be less independent) than QTL further apart. Thus, some spatial covariance structure along the chromosome might be in order. Ridge regression overstates the prior precision of the model, by virtue of implicitly assuming independence. A possible consequence is that "more genotype is marked" than warranted. Also, there may be chromosomes that are "hotbeds" of QTL, whereas other chromosomes may be arid, leading to between-chromosome variation. This heterogeneity can be accommodated by introducing chromosome effects in the model, with markers in different chromosomes having (possibly) distinct distributions. These extensions are dealt with in the following section.

	EXTENDING THE HIERARCHY

TOP ABSTRACT A MIXED-EFFECTS MODEL... EXTENDING THE HIERARCHY A BAYESIAN FORMULATION DISCUSSION LITERATURE CITED

Model without background genetic effects:
For simplicity, assume first that the state of prior knowledge does not allow disentangling background genetic effects from environmental deviates. The first tier is then as in (20), so that

(22)

The second tier is given by the distribution of the marked genetic effects. Here, instead of assuming that |² N(0, I²), a "one-way layout" is adopted, to partition the variability of marked effects into between- and within-chromosome components. Arrange markers sequentially according to their order within each chromosome and assume a dense marker map. The model for the second tier is

(23)

where c = {c_i} is an n_c x 1 vector of "marked chromosomal effects" (n_c is the number of pairs of chromosomes), T is a known incidence matrix of order l x n_c relating marked effects to chromosomes, and v = {v_ij} is an l x 1 vector of marked "within-chromosome" deviates. In scalar notation, _ij = c_i + v_ij is the "marked quantitative effect of the jth marker at the ith chromosome." It is assumed that

(24)

Here, V = Diag(V_i) is taken to be a block-diagonal matrix with n_c blocks, where V_i is the l_i x l_i variance-covariance matrix of marked within-chromosome deviates in chromosome i, and l_i is the number of markers in chromosome i. The block diagonality of V in (24) implies that marked within-chromosome effects are independent across chromosomes. However, a within-chromosome dependence will be introduced. The matrix V_i can be structured in several different manners, and some forms of modeling such dependence are discussed below. In (23) and (24) the n_c x 1 vector of chromosome effects is assumed to possess the distribution

(25)

where ²_c is the variance between chromosome effects. Some of the possible forms of V_i (among many possible ones) are considered next:

1. Within-chromosome deviates are correlated according to a first-order autoregressive process. Partition with
   

   

   Then, if markers are equally spaced,
   

   

   where ²_vi is the variance between deviates in chromosome i (i = 1, 2, ... , n_c) and is a parameter taking values between -1 and 1. This covariance structure satisfies the assumption that adjacent within-chromosome deviates are more strongly correlated than those further apart. A form of relaxing the requirement that markers be equally spaced is immediately below.

2. Within-chromosome deviates are correlated according to a Gaussian decay model (e.g., VERBEKE and MOLENBERGHS 1997 ). The covariance structure here is
   

   

   where
   

   

   is the correlation between marked within-chromosome deviates k and t in chromosome i. Here, d_i,kt is the distance (e.g., in physical units such as kilobases) between markers k and t in chromosome i, and _i > 0 is a chromosome-specific parameter. When the distance between markers 0, _i,kt 1; on the other hand, when d_i,kt , the correlation between within-chromosome deviates is null. The parameter _i governs the rate at which the correlation decreases. When _i is close to 0, such correlation falls rapidly; when _i increases, the drop is more gentle. It may be convenient to assume that _i = for all chromosomes.

3. Suppose that the number of markers is constant across chromosomes (with l markers and n_c chromosomes there would be z = l/n_c markers per chromosome). If the markers are equally spaced and if the variance of within-chromosome deviates is constant across chromosomes, one may pose the Toeplitz covariance specification
   

   

   (VERBEKE and MOLENBERGHS 1997 ) for i = 1, 2, ... , n_c. Here, the correlation between within-chromosome deviates k and t is
   

   

Note that _kt here, even though a correlation coefficient, is not the same as _i,kt in the preceding section. Using (23) in (22) the model can be expressed as

(26)

Under Gaussian assumptions, the joint distribution of c and v, given a vector or linearly independent error contrasts w, is

(27)

Here, = BLUP(c) and = BLUP(v) are computed as

(28)

and

(29)

The BLUP of the M is . Since the predictor is conditionally unbiased (given M), it is also unbiased unconditionally. The unconditional covariance matrix of the prediction errors is the average of (29) taken with respect to the distribution of marker genotypes.

Here, it is possible to predict marked effects on a chromosomal basis. Consider model (26) and recall that = Tc + v. Suppose that a hypothetical species has four pairs of chromosomes. For the n individuals assayed, the incidence matrix M can be written as

where m^'_ij are the codes for the marker genotypes of individual i at chromosome j. Likewise, the "true" marked effects on a chromosome-by-chromosome basis can be written as

where c_j is an effect peculiar to all markers in chromosome pair j and 1_lj is a vector of 1's of order l_j. Then, after solving (28), one has the following relationships:

.

For many reasonable covariance structures, all dispersion parameters are identifiable in this model and can be estimated by tailoring some suitable algorithm for maximum likelihood. VERBEKE and MOLENBERGHS 1997 , VERBEKE and MOLENBERGHS 2000 present examples of applications of standard statistical software for restricted maximum-likelihood estimation.

Model with background genetic effects:
The model is now an expanded version of (26):

(30)

The best linear unbiased estimator of ß and the BLUP of c, v, and can be obtained by solving the system:

(31)

As before, if an estimate of the total genetic variance is available, one can set ² as ² = ^2* - (1 - ²) ²_y. Some estimates (e.g., REML) of ²_c and of the parameters defining the structure of V_i can be obtained when estimating ^2*_e from the model where background genetic effects are ignored or from the full model but fixing the value of ². The BLUP of the "total" genetic value of the entire collection of individuals in the sample is then

(32)

and the covariance matrix of the prediction errors is

(33)

Differential contribution of marked effects to genetic variability:
A representation of the degree of inequality of a frequency distribution is given by the Lorenz curve (e.g., GASTWIRTH 1971 , GASTWIRTH 1972 ). This is used for measuring income inequality in economics and was applied first in quantitative genetics by URIOSTE et al. 2001 in an assessment of heterogeneity of variance of milk yield in cattle among herds. We place the cumulative proportion of the markers 1/l, 2/l, ... , i/l, ... , (l - 1)/l in the abscissa and plot this against the cumulative proportion of a measure of the marked "genetic variability." A hypothetical example is in Fig 1. A straight line with a slope of 45° indicates perfect equality. For example, if 20% of the markers explain 20% of the variability, 40% of the markers account for 40%, and so on, a straight line results. On the other hand, if a few markers account for a large portion of the variability, the curve bends downward because the distribution of effects is unequal. A measure of the inequality of contributions to variability is the Gini coefficient or twice the "inequality area," that is, the area between the Lorenz curve and the straight line. Since the square in Fig 1 has an area equal to 1, the maximum value of the inequality area is 0.5. Hence, the Gini coefficient takes values between 0 (all marked effects contributing equally to variability) and 1 (a single marked effect explains all variability).

View larger version (20K): In this window In a new window Download PPT slide   Figure 1. Lorenz curve depicting the relative contribution of marked effects to marked genetic variability. Ordinate, cumulative proportion of marked genetic variance explained by marked effects; abscissa, cumulative proportion of markers.

Let M_unique be a set of rows of M defining u unique marker genotypes (u may be much smaller than n), m^'_unique,j be the jth row of M_unique, and m_unique,jk be the entry for marker locus k in m^'_unique,j. Define the corresponding "unique marked effect" as

(34)

Arbitrarily, define "total marked genetic variability" as the sum of the squared unique marked genetic effects, so

(35)

This definition does not take into account the frequency of the unique marker genotypes in the population, but this can be accounted for easily by weighting each µ_j appropriately. Let now µ²_[j] be the jth ordered value of µ²_j, sorted in an increasing order. The ordinate values in the Lorenz curve, i.e., the cumulative proportion of the "total marked genetic variation," are calculated as

(36)

where i/u is the cumulative proportion of unique marker genotypes. The Lorenz curve results from plotting L(i/u) against i/u; observe that L(0) = 0 and L(1) = 1. It can be shown (AGUILAR-GUTIERREZ 2000 ; URIOSTE et al. 2001 ) that the Gini coefficient is approximately equal to

(37)

A main difficulty in estimating the Lorenz curve and the Gini coefficient is that nonlinear functions of unknown quantities are involved. Simple method of moment estimators and are obtained by replacing µ²_j in (36) and (37) by

with obtained from (31). This statistic does not take into account the uncertainty associated with .

Family heterogeneity:
Suppose that the phenotypic and molecular marker information can be partitioned into F nonoverlapping or "independent" clusters representing some familial aggregation. For example, XU 1998 discussed the use of "families" of line crosses for QTL mapping. This author pointed out that single-family methods present the risk that if two lines involved in the cross are not segregating at a QTL, then the latter would not be detected, irrespective of the number of individuals scored in a backcross or F₂ population. He allowed for multiple families in the model, but treated their effects as fixed. Irrespective of whether these effects are viewed as fixed or random, they constitute statistical nuisances from the point of view of QTL mapping (XU 1998 ). Here, a different point of view is taken. Allowance is made for heterogeneity of family effects on the marker-phenotype relationships. The stylized situation described in standard treatments of quantitative genetics is adopted; for example, individuals may be clustered into full-sib or half-sib families. We emphasize that our approach does not include recombination rates between markers and QTL as model parameters. A rigorous treatment for outbred populations requires extremely involved algebra (e.g., WANG et al. 1998 ). The purpose here is to bypass these complex procedures via straightforward linear models having a much simpler covariance structure.

Let the vector of phenotypic records be partitioned as y = [y^'₁, y^'₂, ... , y^'_F]^', where y_i contains f_i observations pertaining to family i (i = 1, 2, ... , F). Assume now that heterogeneity in the marker-phenotype association exists. In the presence of family-specific marked effects the entire vector can be written as

where _i is an l x 1 set of effects peculiar to family i. The model for phenotypic values, allowing now for family-specific background effects, _F, becomes

(38)

where X = [X^'₁, X^'₂, ... , X^'_F]^' and X_i has as many rows as there are observations for family i (f_i) and p columns. Further, M = Diag(M_i). The incidence matrix M_i has f_i rows and l columns, the number of markers. The F x 1 vector _f = {_i} contains F family background effects related to the phenotypic values via the known incidence matrix F (consisting of vectors of ones, 1_i, in appropriate locations and of zeroes elsewhere); w is an n x 1 residual vector of independently distributed within-family deviations, with common variance ²_w and where n = ^F_i=1f_i.

In the absence of marker information, the phenotypic variance can be partitioned into between- and within-family components. For example, with half-sib families, the between-family variance contains one-quarter of the additive genetic variance (supposing additive inheritance), and the within-family component includes the environmental variance plus three-quarters of the additive genetic variance. In the setting of (38) it is assumed that _f and w are independently distributed vectors, with distributions _f|²_f N(0, I²_f) and w|²_w N(0, I²_w), respectively. This partitioning reassigns the unmarked additive genetic variance into between (²_f) and within-family (²_w) components. Now, the _i marked effects can be treated as random regressions following the distribution

(39)

Here, ₀ is an l x 1 parameter that is common to all families (the "fixed" part of the regression) and B is the l x l variance-covariance matrix between the regression coefficients. This implies that _i = ₀ + _i, where _i N(0, B) is a vector of deviations from the common regression ₀ that is specific to family i. In a simpler hierarchical model it could be postulated, for example, that _i|² N(0, I_i²) for i = 1, 2, ... , F. The model for the more general specification becomes

(40)

where ß*' = [ß', ^'₀], X* = [X M₀], M₀ = [M^'₁ M^'₂ · M^'_F]', and ' = [₁, ^'₂, ... , _F] is a vector of family-specific deviations from the overall regression in (39).

Given the dispersion parameters B, ²_f, and ²_w, the best linear unbiased estimator of ß* and the BLUP of each of the family-specific deviations _i are found by solving the system

(41)

where ^F_i=1 is the "direct sum" of matrices notation (e.g., SEARLE 1982 ), denoting a block-diagonal concatenation of the F matrices M^'_iM_i + B^-12_w, with all these matrices being l x l. Note that

and

The variance-covariance matrix of the prediction errors is given by the inverse of the coefficient matrix in (41) multiplied by ²_w. Unless there are many families (e.g., independent half-sib groups in dairy cattle) it will be difficult to obtain reliable maximum-likelihood estimates of B.

The hierarchy can be extended as in (23) and (24) to include chromosome and within-chromosome effects that would now be family specific. Put

(42)

where _i is l x 1, c_i is an n_c x 1 vector of chromosome effects peculiar to family i, T_i is an l x n_c incidence matrix, and v_i is an l x1 vector of within-chromosome effects for family i. It is convenient to assume that the within-chromosome marked effects are independent across families and chromosomes, but possibly dependent within chromosomes. Here, the variance-covariance matrix V in (24) would be Fl x Fl and would take the form

where V_i = Var(v_i) is an l x l variance-covariance matrix that is specific to family i. Now, partition

where v_ij is an l_j (number of markers in chromosome j) x 1 vector of within-chromosome deviates for family i in chromosome j. Then, under the assumption of independence of deviates across chromosomes,

With this,

and V_ij can be assigned any of the spatial structures discussed earlier. Then, the distribution of all family-specific deviations from regression would be

(43)

The hierarchy is completed by assigning a distribution to the family-specific chromosome effects c_i. An extension of assumption (25) is

(44)

where K = {k_ij} is an n_c x n_c matrix of covariances between chromosome effects of the same family, with this dispersion structure assumed homogeneous across families. For example, element k₄₅ of K would be the covariance between the effects of chromosomes 4 and 5 within the same family. As a side note, observe that unconditionally to c_i, (43) and (44) imply that _i N(0, T_iKT^'_i + V_i). Hence, assumption (39) is equivalent to (43) and (44) if and only if T_iKT^'_i + V_i = B for all i.

Using (42) in (40) leads to the following model for phenotypes,

(45)

where M = Diag(M_i) and T = Diag(T_i). The BLUP of all family-specific chromosome and within-chromosome effects is obtained by solving

(46)

The variance-covariance matrix of the prediction errors is the inverse of the coefficient matrix in (46) times ²_w.

Lande-Thompson regression model revisited:
LANDE and THOMPSON 1990 suggested splitting the molecular scores into between- and within-family components, so that the components could be combined "optimally" for the purpose of genetic improvement. In our context and for simplicity, consider the following variant of their model, using (38). Write the marker codes for family i as M_i = _i + (M_i - _i), where _i is an f_i x l matrix with columns equal to the mean values of the corresponding columns of M_i, with all its rows being equal. A model of interest might be

(47)

where _B is the regression of phenotypes on mean family molecular scores, M_W is a matrix of within-family deviations from the mean scores for the appropriate families, and _W is the vector of regressions of phenotypes on the within-family deviations. Assume that _B N(0, I_F²_F) and _W N(0, I_F²_W) where ²_F and ²_W are components of variance, and suppose that the two random vectors are independently distributed. The best linear unbiased predictor of the between- and within-family regressions can be calculated by solving

(48)

The variance-covariance matrix of the prediction errors can be obtained from the inverse of the coefficient matrix in (48), times ²_w. Model (47) can be extended by allowing the within-family regressions to be family-specific. In this extended model, in addition to _B, there would be _Wi (i = 1, 2, ... , F) within-family regressions.

	A BAYESIAN FORMULATION

TOP ABSTRACT A MIXED-EFFECTS MODEL... EXTENDING THE HIERARCHY A BAYESIAN FORMULATION DISCUSSION LITERATURE CITED

We adopt now a Bayesian point of view and describe a Markov chain Monte Carlo (MCMC) implementation. The focus is on model (40) and its hierarchical expansion in (42–44), leading to (45) and (46) in the mixed model and BLUP treatment. This is the most general and richly parameterized specification among those discussed so far. The developments follow the typical hierarchical structure of Bayesian multilevel models (LINDLEY and SMITH 1972 ).

Joint posterior distribution:
Start with (40) as sampling model (data-generating process) or first level of the Bayesian hierarchy. Conditionally on ß*, , and _f, it is assumed that the phenotypic values are drawn independently from the distribution

(49)

It is convenient to express the joint density of all phenotypes, given ß*, , _f, and ²_w as the product of F densities corresponding to the contributions to information made by the phenotypic values in each of the families. Thus,

(50)

The notation N(y_i|X^*_iß^* + M_ii + 1_ii, I²_w) indicates a normal density or distribution with y_i as random variable, X^*_iß* + M_ii + 1_ii as mean vector, and I²_w as covariance matrix; a similar notation is used if the distribution involves a scalar.

The second level assigns a prior distribution to all unknown parameters of the first tier. It is assumed that the joint prior density can be written as

(51)

Above, ß^*₀ is the known mean of the prior distribution of ß* and s²_ß* is a known scalar that tunes the degree of vagueness of the prior for this parameter. Different degrees of vagueness may be assigned to the two components of ß*; thus, there would be two distinct tuning parameters, s²_ß and s²₀. The notation V() means that the variance-covariance matrix of v (the within-chromosome deviates) depends on some parameter vector . For example, if within-chromosome deviates are correlated according to an autoregressive process having as parameter and with chromosome-specific variances ²_v1, ²_v2, ... , ²_vnc, then = [, ²_v1, ²_v2, ... , ²_vnc]'. The autoregressive process is assumed hereafter, to illustrate one of the possible specifications. In (51), s²_w and _w are known parameter values (hyperparameters) of a scaled inverted chi-square distribution (e.g., GELMAN et al. 1995 ) assigned to ²_w. Now, as stated in the developments following (42) and leading to (43), given the chromosome effects c, the marked effects of different families are taken to be mutually independent. Hence, (51) can be put as

(52)

The third level of the hierarchy consists of the prior distribution assigned to c, , and ²_f, the parameters of the second tier. It is assumed that the corresponding density takes the form

(53)

In the preceding, it is assumed that (a) the chromosome effects of different families are mutually independent, although correlated within a family with an n_c x n_c covariance matrix K, as in (44); (b) the parameter has a prior distribution indexed by some known parameter vector and is independently distributed (a priori) of all within-chromosome variances ²_vi; (c) the n_c variances ²_vi are independently and identically distributed as scaled inverted chi square with known hyperparameters s²_v and _v; and (d) the variance between families ²_f follows a scaled inverted chi-square distribution with known parameters s²_f and _f.

The fourth and final level of the hierarchical model is the prior distribution assigned to K, the covariance matrix between effects of the same family on different chromosomes. It is assumed that K follows an inverted Wishart process of order n_c (the order of K) with density p(K|_cS_c, _c), where _cS_c is a known scale matrix and _c is a known positive parameter usually referred to as the "degrees of freedom" of the distribution (e.g., GELMAN et al. 1995 ). Our parameterization is such that, a priori, E(K|_cS_c, _c) = _cS_c/(_c - n_c - 1). For a very large value of _c, S_c then approximates the mean value of the prior distribution. In less structured models, e.g., independence of chromosome effects within a family, K would take a simpler form, in which case the prior distribution would be modified appropriately.

Some of the prior distributions given above are difficult or impossible to elicit. In such a situation, one may consider using some of the standard default improper priors, carry out a technically involved reference prior analysis (BERNARDO 1979 ; BERNARDO and SMITH 1994 ), or resort to maximum entropy fits; see CANTET et al. 1992 and SORENSEN and GIANOLA 2002 for an elementary discussion of the latter in a genetics context. Collecting (50), (52), (53), p(K|_cS_c, _c) and letting H denote the set of all known hyperparameters, the joint posterior density of the uncertain parameters is, after rearrangement,

(54)

Markov chain Monte Carlo scheme:
At least in theory, all marginal or joint distributions of sets of parameters of interest can be derived from (54) by effecting the necessary integrations, to take uncertainty about nuisance parameters into account properly. Since the joint posterior distribution is not in a form amenable to analytical treatment, MCMC (e.g., GILKS et al. 1996 ) can be used for simulating draws from the posterior. Features of any target posterior distribution can be estimated from the simulated samples. Monte Carlo error can be made negligible by sheer computing force; however, suitable reparameterizations may enhance the efficiency of the simulation.

The hybrid algorithm proposed uses the conditional posterior distributions as candidate-generating densities, whenever possible. Save for the presence of p(|) (whose structure has not been specified yet), the joint posterior process with density (54) is in a normal-inverse Gamma or normal-inverse Wishart form. All fully conditional distributions (except that of ) can be identified and sampled with relative ease. This defines a Gibbs sampler for drawing from most posterior distributions of interest. Since the conditional process is not recognizable for , a Metropolis procedure is tailored. The most relevant expressions needed for implementing this hybrid algorithm are presented below. The notation [parameter or parameters|ELSE] is used throughout to denote a fully conditional posterior distribution, that is, the posterior distribution of a scalar or vector parameter given all other parameters, the data, and H. For details, see WANG et al. 1993 , WANG et al. 1994 , SORENSEN et al. 1994 , and SORENSEN and GIANOLA 2002 .

Conditional distribution of ß*: This is arrived at by retaining in (54) only the terms involving ß*. The resulting conditional density is

This can be identified as the density of the multivariate Gaussian distribution,

(55)

where

and

Obtaining samples from (55) is straightforward, especially if the order of ß* is not too large. When the prior distribution of ß* is diffuse (s²_ß* ), this step involves essentially least-squares type computations, after making an offset of the data vector.

Conditional distributions of family-specific marked effects : Inspection of (54) reveals that all family-specific marked regressions _i are mutually independent, given all other parameters. The density of the conditional posterior distribution of _i (i = 1, 2, ... , F) is

The distribution can be shown to be the l-variate normal process

(56)

where

and

All family-specific marked effects can be drawn on a family-by-family basis by sampling from the multivariate normal distribution given in (56).

Conditional distributions of family-specific background effects : From (54), it follows that

so that family-specific background effects are conditionally independent, given all other parameters. In particular, for the ith family one has

This can be shown to be the kernel of the density of the univariate normal distribution,

(57)

where

and

Note that 1^'_i1_i = f_i is the number of individuals in family i and that ²_f/(²_f + ²_w/f_i) is proportional to the "heritability of a family mean" (whenever the variance of background effects between families is proportional to the additive genetic variance), in the terminology of standard quantitative genetics. Sampling from distribution (57) is straightforward.

Conditional distributions of chromosome effects: The density of the conditional distribution of the vector of chromosome effects c is, from (54),

Thus, family-specific chromosome effect vectors c_i (i = 1, 2, ... , F) are conditionally independent of each other. After algebra, as in GIANOLA and FERNANDO 1986 or GIANOLA et al. 1990 , it can be shown that the ith conditional distribution is the n_c-variate normal process,

(58)

where

and

For instance, in a species with 30 pairs of chromosomes such as cattle, a draw would need to be made from a 30-variate normal distribution for each of the F families represented in the data set.

Between- and within-family variances of background effects: The density of the conditional posterior distribution of the between-family variance of background effects, ²_f, is

Writing the F Gaussian densities and the density of the scaled inverted chi-square prior process of ²_f explicitly, one arrives at

This indicates that the conditional posterior distribution of ²_f is scaled inverted chi square (e.g., WANG et al. 1993 ; GELMAN et al. 1995 ), that is,

(59)

To draw samples from (59), one extracts a random deviate from a central chi-square distribution on F + _f degrees of freedom, takes its reciprocal, and multiplies the inverted deviate by (^F_i=1²_i + s²_ff).

Likewise, the density of the conditional posterior distribution of the within-family variance is

Writing the densities explicitly and rearranging one arrives at

This is the distribution of the scaled inverted chi-square random variable

(60)

The sampling process is similar to that described in connection with (59).

Autoregressive parameter : Note in (54) that the only terms involving are the prior densities of the family-specific marked effects _i and the prior density p(|). Hence,

Now, recall that V_i() is an l x l covariance matrix such that

where V_ij() is the l_j x l_j covariance matrix of within-chromosome deviates for family i in chromosome j. Specifically,

(61)

Note that this covariance matrix is chromosome specific. Using this, the conditional posterior density above has the form

(62)

When viewed as a function of , (62) is not in a recognizable form, irrespective of the form of the prior density p(|) adopted for the autoregressive parameter. Effecting a draw from this distribution requires a more involved procedure, such as a single-site Metropolis or Metropolis-Hastings algorithm (METROPOLIS et al. 1953 ; TANNER 1993 ; GELMAN et al. 1995 ; GILKS et al. 1996 ; SORENSEN and GIANOLA 2002 ). Perhaps the simplest is to adopt a Metropolis step. Here, a proposal distribution is needed for generating candidate values of . As in FISHER 1921 , and for the purpose of facilitating the generation of proposals, consider transforming the correlation coefficient into as

with the Jacobian of the transformation being

Then, the conditional posterior density of is

(63)

where = [, ²_v1, ²_v2, ... , ²_vnc]^'. Suppose now that the state of at iteration t of the algorithm is ^[t] and that all other parameters in its conditional posterior distribution have been updated. This implies that the state of is

We update via a Metropolis jump. Here, a candidate value * must be sampled from a symmetric candidate generating distribution with density Q(^*|^[t]); symmetry means that Q(^*|^[t]) = Q(^[t]|^*) for all pairs (^*, ^[t]) and for every t (GELMAN et al. 1995 ). A distribution meeting this requirement is a t-process; the t-distribution has thicker tails than the normal, so "extreme" values appear at higher density, specially when the degrees of freedom are low (4–10, say). A possible Metropolis jump is as follows.

First, sample a deviate r from a Gamma(d.f./2, d.f./2) where d.f. denotes the degrees of freedom of the t-distribution.

Second, following BOX and TIAO 1973 , sample a candidate * from

Third, using (63), compute the Metropolis ratio:

The integration constant cancels out in the numerator and denominator, so is calculated by direct evaluation of ratios involving (63).

Fourth, draw a deviate U from a Uniform(0, 1) distribution, and set

Invert the candidate as

and set

The preceding completes the updating of the autoregressive parameter.

Variances of within-chromosome deviates ²_vi: From (54),

It is convenient to rewrite the model for marked effects presented in (42) into a chromosome-specific basis. Write

where ^*_j is a vector of order Fl_j x 1 containing the marked effects of the l_j markers in chromosome j for the F families, and v^*_j is a vector of within-chromosome deviates. An assumption made earlier is that within-chromosome deviates are conditionally mutually independent across chromosomes and families, but correlated within chromosomes according to the covariance matrix in (61). Since the preceding is merely a rearrangement of (42), it follows that

where the identity matrix is F x F. Hence,

indicating that the within-chromosome variances ²_vj have independent fully conditional posterior distributions. In particular, and writing the densities explicitly, one has for ²_vj,

(64)

This is the density of the scale-inverted chi-square variate:

(65)

Thus, the within-chromosome variance can be sampled independently, chromosome by chromosome.

Covariance matrix of chromosome effects K: Retaining in (54) only those terms involving K leads to

Now, we write explicitly the F normal densities and the inverse Wishart density p(K|S_c, v_c), to obtain

(66)

This indicates that the fully conditional posterior distribution is inverse Wishart of order n_c, degrees of freedom equal to F + v_c, scale matrix ^F_i=1c_ic^'_i + S_cv_c, and posterior expectation

(given all other parameters). Procedures for sampling from inverse Wishart distributions in a quantitative genetics context are in, e.g., JENSEN et al. 1994 , KORSGAARD et al. 1999 , and SORENSEN and GIANOLA 2002 .

The MCMC scheme in a nutshell:
A discussion of how Markov chain Monte Carlo schemes can be tuned, run, and monitored for convergence can be found in GELMAN et al. 1995 , GILKS et al. 1996 , COWLES and CARLIN 1996 , ROBERT and CASELLA 1999 , and SORENSEN and GIANOLA 2002 . A possible implementation follows.

1. Set starting values that are hopefully "inside" of the target joint posterior distribution.

2. Sample systematically from each of the following conditional posterior distributions, perhaps changing the order of visitation at each round: (55) for ß*; (56), (57), and (58) for the family-specific marked effects, background effects, and chromosome effects, respectively; (59) and (60) for the between- and within-family variances of the background effects; (62) for the autoregressive parameter, after incorporating the Metropolis step; (65) for the n_c variances of within-chromosome deviates; and (66) for K, the covariance matrix of chromosome effects. This constitutes a single completed iteration of the scheme. At each subsequent loop, update the values of the appropriate conditioning values with the new draws.

3. Repeat the iteration as many times as needed to ensure (a) that draws can be reasonably claimed as belonging to the posterior distribution and (b) that posterior features are estimated with a sufficiently small Monte Carlo error.

4. If necessary, discard early iterations, in what is called the "burn-in" period.

5. Using the remaining samples, estimate any feature of the posterior distribution, e.g., a posterior mean or variance, a posterior density, or the distribution of an order statistic.

To illustrate the flexibility of the approach, consider inferring the Lorenz curve and the Gini coefficient in (36) and (37), respectively, and contrast this with the crude methods discussed in Differential contribution of marked effects to genetic variability. Let ^(k) be a draw from the posterior distribution of the marked effects in the simplest model. Then, as in (34), µ^(k)_j = m^'_unique,j^(k) is a draw from the posterior distribution of µ_j and S(^(k)) is a sample from the posterior distribution of S() in (35). Next, order all the µ^2(k)_j, so that

is a sample from posterior distribution of the ordered squared "unique" marked effects. Then

is a sample from the posterior distribution of the Lorenz curve in (36). Likewise,

is a sample from the posterior distribution of the Gini coefficient in (37), measuring unequal contribution of marked effects to genetic variability. The Bayesian analysis produces an entire description of the uncertainty about the Lorenz curve and the Gini coefficient.

Draws from the posterior distribution of the "total" genetic value (in the simplest model) would be obtained as M^(k) + ^(k), and posterior means, variances, and covariances can be computed directly from such samples.

	DISCUSSION

TOP ABSTRACT A MIXED-EFFECTS MODEL... EXTENDING THE HIERARCHY A BAYESIAN FORMULATION DISCUSSION LITERATURE CITED

Several methods for incorporating molecular information into predictors of genetic merit of candidates for selection in improvement programs have been proposed in recent years. For example, LANDE and THOMPSON 1990 suggested using multiple regression (in a least-squares sense) of the phenotypic value of an individual on the number of copies of alleles at polymorphic loci; in situations where the linkage map is known, separate regressions could be calculated for each chromosome. They argued that, given a sufficient number of markers, nearly all of the additive genetic variance could be accounted for by the regression model. These authors calculated a "molecular score" (the fitted value of the least-squares regression) and combined the molecular score with the phenotype for the quantitative trait via standard selection index techniques. LANDE and THOMPSON 1990 also introduced a familial structure and discussed how selection could be optimized (in some sense) by combining family means and within-family deviations for both the molecular score and the phenotypic values via selection index theory (SMITH 1936 ; HAZEL 1943 ). There are some problems with this approach. First, the method breaks down when formulated in a vectorial manner. Even in the trivial case where candidates are independent and identically distributed (essentially their setting), the covariance matrix of the molecular scores is singular since the distribution of fitted values in regression is defined only in a p-dimensional space, rather than in n dimensions (the number of individuals with molecular scores). Hence, a selection index approach coined in the more general framework of best linear prediction (HENDERSON 1973 ) leads to an infinite number of solutions. Hence, care must be exercised to ensure that the functions of genetic merit to be predicted are indeed predictable (HARVILLE 1976 ).

Another difficulty with the least-squares method arises when the number of markers is almost as large as the number of individuals typed. For example, >1.5 million single-nucleotide polymorphisms (the current desideratum of genetic marker) have been identified in the human genome and their positions located at an average spacing of 2 x 10^-3 cM (HARTL and JONES 2002 ). If advantage is to be taken of this map density, this leads to a situation where l >> n. Hence, colinearity of predictor variables creates an insurmountable parameter identification problem unless some dimension-reduction technique, such as singular value decomposition, or some ad hoc marker model selection procedure is employed. Also, it is well known in statistics that least-squares is an inadmissible estimation procedure. STEIN 1955 showed for the orthonormal linear regression model that whenever l > 2, the least-squares estimator is dominated (in the mean-squared error sense) throughout the parameter space by a nonlinear, shrinkage type, estimator. This result, which torpedoes best linear unbiased estimation and maximum-likelihood estimation under normality, has been ignored in quantitative genetics (GIANOLA 1990 ), although it is widely acknowledged in econometrics (THEIL 1971 ; JUDGE et al. 1985 ). On these grounds, not much should be expected from least-squares derivatives in marker-assisted prediction of genetic value, as found by MEUWISSEN et al. 2001 .

WHITTAKER et al. 2000 and WHITTAKER 2001 noted that having "too many" markers in the regression model produces serious colinearity, causing unstable least-squares estimates and a poor prediction of the molecular score. The method leads to apparent paradoxes, as follows. At least in theory, the regression coefficients of phenotype on markers flanking a QTL must have the same sign (WHITTAKER et al. 1996 ). However, the estimated regression coefficients can differ in sign even when the true parameters agree in sign (HWANG and NETTLETON 2002 ). The latter authors point out that sign inconsistency is greatest when flanking markers are close together and, thus, highly correlated. As an alternative to selection of subsets of markers to include in the model, WHITTAKER et al. 2000 and WHITTAKER 2001 suggested employing ridge regression (HOERL and KENNARD 1970 ) for estimation. This ad hoc procedure shrinks the least-squares estimates toward zero, improves the condition of the coefficient matrix in least-squares, and often produces a reduction in mean-squared error of estimation. Here, one encounters the typical dilemma: the shrinkage estimator may be more precise but less accurate, since it is biased. Whether or not this reduces mean-squared error depends on a bias vs. variance trade-off, and the performance of ridge relative to least-squares regression is a function of the true value of the regression vector, which is obviously unknown. One never knows how much shrinkage should be effected, although there are some informal recipes for tuning the shrinkage parameter. On the other hand, ridge regression makes more sense when viewed from a Bayesian perspective (LINDLEY and SMITH 1972 ; ZELLNER and VANDAELE 1975 ; GIANOLA and FERNANDO 1986 ). Ridge regression is equivalent to adopting a normal prior for the regression vector centered at zero and with a prior covariance structure equal to an identity matrix times a scalar, the variance of the prior distribution. A small variance means large prior precision and more shrinkage (a stronger belief that the regression is null); a large variance implies diffuse prior knowledge of the regression vector. The scalar prior covariance matrix implies that the individual regression coefficients of phenotypes on markers are viewed as independently and identically distributed a priori. Shrinking all marked effects toward a common prior mean may be questionable. Also, assuming that different chromosome sections have independent expression levels is contradicted by evidence from gene expression trials (CARON et al. 2001 ).

In practice, all individuals cannot be typed for all markers. Further, if typing costs continue to be abated, the number of predictor variables (markers) may grow faster than the number of observations. This situation is encountered routinely in animal breeding, and it can be solved by positing random-effect models. The BLUP methods used in animal breeding can be viewed as based on conditional or posterior distributions (GIANOLA and FERNANDO 1986 ), illustrating why the dimensionality and identification problems arising in least-squares are circumvented. MEUWISSEN et al. 2001 noted this and compared least-squares with BLUP and with two Bayesian methods for marker-assisted prediction. In their simulation, effects of 50,000 marker haplotypes needed to be inferred from 2200 observations. Since this is not possible via least-squares, they calculated regressions on markers found in chromosome segments of 1 cM each and evaluated a log-likelihood (assuming normality) for every segment. Then, "significant" segments were selected using an arbitrary procedure on the basis of the said likelihoods. Eventually, the authors combined all such segments into a single multiple regression model. In the BLUP analysis, all chromosome effects were regarded as independent and identically distributed random effects with known variance. Their first Bayesian model was as in BLUP, but assigning independent scaled inverted chi-square distributions to the variances of the chromosome segments. The second Bayesian model used a mixture prior distribution for the variance of the chromosome segment variances. Their main conclusions were: (1) least-squares had a low accuracy of prediction of genetic value, (2) the BLUP and Bayesian methods differed by little, and (3) the Bayesian methods assigning prior distributions to the variance associated with chromosome segments gave the most accurate predictions.

The classical and Bayesian methods outlined here do not accommodate complex genealogical structures such as those arising in animal breeding. To illustrate, consider the hypothetical situation in which the markers and the QTL are identical, so that there are no background effects if gene action is additive. Hence, in the notation of A MIXED-EFFECTS MODEL FORMULATION, a_i = ^l_k=1m_ikk, and, therefore, one would expect a priori that the covariance between the genetic value of any two relatives i and i' is equal to

In the absence of knowledge of genotype frequencies, the term E(m^'_im^'_i) cannot be assessed. In a fully Bayesian method one needs to assign a prior to the genotypic frequencies and take the MCMC method even further (e.g., HOESCHELE 2001 ). For this reason, all analyses presented in our article treat the matrix M as nonstochastic. Concerning the background genetic effects, one can introduce a full relationship matrix under additive inheritance.

Model (38) introduces a family structure, and (47) expands LANDE and THOMPSON 1990 to a matrix domain, but with a random-effects treatment. To illustrate a Bayesian analysis of the latter, consider the following terms: (1) _iB, the f_i x 1 vector of unobserved (contrary to least-squares fitted values) mean molecular scores for all members of family i (all its elements are equal); (2) (M_i - _i)_W, the within-family contributions to the molecular scores; (3) _i or background genetic effect of family i (e.g., the additive genetic contribution, net of the marked part of the variability); and (4) w_i, a vector of within-family deviations (also net of the marked part of the variability). The mean of the posterior distribution of the total genetic merit of the individuals in family i would be obtained by averaging _i^(k)_B + (M_i - _i)^(k)_W + 1_i^(k)_i over the MCMC samples. Similarly, the posterior mean of the within-family deviations can be estimated by averaging w^(k)_i = y_i - X_iß^(k) - M_i^(k)_i - 1_i^(k)_i. LANDE and THOMPSON 1990 consider "selection optimization" by defining total merit of a candidate for selection as a linear combination of family and within-family information. For instance, one could write the merit vector

where h_F and h_F are some nonstochastic relative weights assigned to the unobserved molecular and quantitative trait family means, respectively, and h_w, h_w are corresponding weights assigned to the unobserved within-family deviations. The means, variances, covariances, and any feature of a posterior distribution involving the unobserved total merits T_i can be estimated from the MCMC samples. The procedure takes the ideas of LANDE and THOMPSON 1990 to a higher level of generality, with the uncertainty about nuisance parameters taken into account fully in the usual Bayesian manner (BOX and TIAO 1973 ).

We have assumed throughout that all individuals have been typed for all markers, but this is seldom the case, as noted earlier. In the Bayesian context missing markers can be dealt with automatically via an augmentation of the posterior distribution. Subsequently, MCMC is used to make imputations for the unobserved part of the molecular information, at least if missingness is at random in the sense of RUBIN 1976 . See HOESCHELE 2001 for technical details.

It remains to be seen to what extent the proposed procedures hold in practice or in simulations what they seem to promise on theoretical grounds. In a simulation, WHITTAKER et al. 2000 found that a simple ridge regression procedure enhanced response to selection over plain phenotypic selection or selection based on least-squares fits of the molecular scores. The results of MEUWISSEN et al. 2001 point in the same direction. This is encouraging. Another issue is the computer intensiveness of the required calculations. The BLUP-REML analysis can be carried out with standard software for large data sets. On the other hand, the MCMC procedures require a higher level of computation; however, they are more flexible and lead to richer inference (e.g., estimation of the Lorenz curve). At any rate, computational burden is unimportant relative to the cost of genotyping individuals. Also, the Bayesian analysis requires elicitation of some complex priors distributions and this is difficult. In practice, a possibility is to use some of the standard default noninformative priors (e.g., SORENSEN and GIANOLA 2002 ) although at the risk of causing impropriety. Another option is to carry out the Bayesian analysis with different priors. If inferences are invariant with respect to the set of priors used, this would suggest that the data are "informative enough" and that the prior is overwhelmed by the information contributed by the data.

	ACKNOWLEDGMENTS

The authors thank Rohan Fernando, Davorka Gulisija, and Guilherme Rosa for useful comments. Research was supported by the Wisconsin Agriculture Experiment Station and by grants from the NRICGP/U.S. Department of Agriculture (99-35205-8162) and National Science Foundation (DEB-0089742).

Manuscript received May 3, 2002; Accepted for publication September 27, 2002.

	LITERATURE CITED

TOP ABSTRACT A MIXED-EFFECTS MODEL... EXTENDING THE HIERARCHY A BAYESIAN FORMULATION DISCUSSION LITERATURE CITED

AGUILAR-GUTIERREZ, X., 2000 Desigualdad y Pobreza en México, son inevitables? Instituto de Investigaciones Económicas, Universidad Nacional Autónoma de México, Mexico City.

BERNARDO, J. M., 1979  Reference posterior distributions for Bayesian inference. J. R. Stat. Soc. B 41:113-147.

BERNARDO, J. M., and A. F. M. SMITH, 1994 Bayesian Theory. John Wiley & Sons, New York.

BOX, G. E. P., 1980  Sampling and Bayes' inference in scientific modelling and robustness. J. R. Stat. Soc. A 143:383-430.

BOX, G. E. P., and G. C. TIAO, 1973 Bayesian Inference in Statistical Analysis. Addison-Wesley, Reading, UK.

CANTET, R. J. C., R. L. FERNANDO, and D. GIANOLA, 1992  Bayesian inference about dispersion parameters of univariate mixed models with maternal effects: theoretical considerations. Genet. Sel. Evol. 24:107-135.

CARON, H., B. VAN SCHAIK, M. VAN DER MEE, F. BAAS, and G. RIGGINS et al., 2001  The human transcriptome map: clustering of highly expressed genes in chromosomal domains. Science 291:1289-1292.[Abstract/Free Full Text]

COWLES, M. K. and B. P. CARLIN, 1996  Markov chain Monte Carlo convergence diagnostics: a comparative review. J. Am. Stat. Assoc. 91:883-904.

FARNIR, F., W. COPPIETERS, J. J. ARRANZ, P. BERZI, and N. CAMBISANO et al., 2000  Extensive genome-wide linkage disequilibrium in dairy cattle. Genome Res. 10:220-227.[Abstract/Free Full Text]

FERNANDO, R. L. and M. GROSSMAN, 1989  Marker-assisted selection using best linear unbiased prediction. Genet. Sel. Evol. 33:209-229.

FISHER, R. A., 1921  On the probable error of a coefficient of correlation deduced from a small sample. Metron 4:3-32.

GASTWIRTH, J. L., 1971  A general definition of Lorenz curve. Econometrica 39:1037-1039.

GASTWIRTH, J. L., 1972  The estimation of the Lorenz curve. Rev. of Econ. Stat. 54:306-316.

GELMAN, A., J. B. CARLIN, H. S. STERN and D. B. RUBIN, 1995 Bayesian Data Analysis. Chapman & Hall, London.

GEORGES, M. D., D. NIELSEN, M. MACKINNON, A. MISHRA, and R. OKIMOTO et al., 1995  Mapping quantitative trait loci controlling milk production in dairy cattle by exploiting progeny testing. Genetics 139:907-920.[Abstract]

GIANOLA, D., 1990 Can BLUP and REML be improved upon?, pp. 445-449 in Proceedings of the 4th World Congress on Genetics Applied to Livestock Production, Vol. XIII. Joyce Darling, Penicuik, UK.

GIANOLA, D. and R. L. FERNANDO, 1986  Bayesian methods in animal breeding theory. J. Anim. Sci. 63:217-244.[Abstract/Free Full Text]

GIANOLA, D. and B. GOFFINET, 1982  Sire evaluation with best linear unbiased predictors. Biometrics 38:1085-1088.

GIANOLA, D., S. IM and F. W. MACEDO, 1990 A framework for prediction of breeding value, pp. 210–238 in Advances in Statistical Methods for Genetic Improvement of Livestock, edited by D. GIANOLA and K. HAMMOND. Springer-Verlag, Heidelberg.

GILKS, W. R., S. RICHARDSON and D. J. SPIEGELHALTER, 1996 Markov Chain Monte Carlo in Practice. Chapman & Hall, London.

HARTL, D. L., and E. W. JONES, 2002 Essential Genetics: A Genomics Perspective. Jones & Bartlett, Sudbury, MA.

HARVILLE, D. A., 1976  Extension of the Gauss-Markov theorem to include the estimation of random effects. Ann. Stat. 4:384-395.

HAYES, B. and M. E. GODDARD, 2001  The distribution of the effects of genes affecting quantitative traits in livestock. Genet. Sel. Evol. 33:209-229.[Medline]

HAZEL, L. N., 1943  The genetic basis for constructing selection indexes. Genetics 28:476-490.[Free Full Text]

HENDERSON, C. R., 1973 Sire evaluation and genetic trends, pp. 10–41 in Proceedings of the Animal Breeding and Genetics Symposium in Honor of Dr. Jay L. Lush. American Society of Animal Science and American Dairy Science Association, Champaign, IL.

HENDERSON, C. R., O. KEMPTHORNE, S. R. SEARLE, and C. M. VON KROSIGK, 1959  The estimation of environmental and genetic trends from records subject to culling. Biometrics 15:192-218.

HOERL, A. E. and R. W. KENNARD, 1970  Ridge regression: biased estimation and applications for non-orthogonal problems. Technometrics 12:55-82.

HOESCHELE, I., 2001 Mapping quantitative trait loci in outbred pedigrees, pp. 599–644 in Handbook of Statistical Genetics, edited by D. J. BALDING M. BISHOP and C. CANNINGS. John Wiley & Sons, Chichester, UK.

HWANG, J. T. and D. NETTLETON, 2002  Investigating the probability of sign inconsistency in the regression coefficients of markers flanking quantitative trait loci. Genetics 160:1697-1705.[Abstract/Free Full Text]

JENSEN, J., C. S. WANG, D. A. SORENSEN, and D. GIANOLA, 1994  Bayesian inference on variance and covariance components for traits influenced by maternal and direct genetic effects using the Gibbs sampler. Acta Agric. Scand. 44:193-201.

JUDGE, G. G., W. E. GRIFFITHS, R. CARTER, H. LUTKEPOHL and T. C. LEE, 1985 The Theory and Practice of Econometrics. John Wiley & Sons, New York.

KORSGAARD, I. R., A. H. ANDERSEN, and D. SORENSEN, 1999  A useful reparameterisation to obtain samples from conditional inverse Wishart distributions. Genet. Sel. Evol. 31:177-181.

LANDE, R. and R. THOMPSON, 1990  Efficiency of marker-assisted selection in the improvement of quantitative traits. Genetics 124:743-756.[Abstract]

LANGE, C. and J. C. WHITTAKER, 2001  On prediction of genetic values in marker-assisted selection. Genetics 159:1375-1381.[Abstract/Free Full Text]

LINDLEY, D. V. and A. F. M. SMITH, 1972  Bayes estimates for the linear model. J. R. Stat. Soc. B 34:1-41.

METROPOLIS, N., A. W. ROSENBLUTH, M. N. ROSENBLUTH, A. H. TELLER, and E. TELLER, 1953  Equations of state calculation by fast computing machines. J. Chem. Phys. 21:1087-1091.

MEUWISSEN, T. H. E., B. J. HAYES, and M. E. GODDARD, 2001  Prediction of total genetic value using genome-wide dense marker maps. Genetics 157:1819-1829.[Abstract/Free Full Text]

OLLIVIER, L., 1998  The accuracy of marker-assisted selection for quantitative traits in populations in linkage equilibrium. Genetics 148:1367-1372.[Abstract/Free Full Text]

ROBERT, C. P., and G. CASELLA, 1999 Monte Carlo Statistical Methods. Springer-Verlag, New York.

RUBIN, D. B., 1976  Inference and missing data. Biometrika 63:581-592.[Abstract/Free Full Text]

SEARLE, S. R., 1974 Prediction, mixed models and variance components, pp. 229–266 in Reliability and Biometry, edited by F. PROSCHAN and R. J. SERFLING. Society for Industrial and Applied Mathematics, Philadelphia.

SEARLE, S. R., 1982 Matrix Algebra Useful for Statistics. John Wiley & Sons, New York.

SEARLE, S. R., G. CASELLA and C. E. MCCULLOCH, 1992 Variance Components. John Wiley & Sons, New York.

SMITH, C. and S. P. SIMPSON, 1986  The use of genetic polymorphisms in livestock improvement. J. Anim. Breed. Genet. 103:205-217.

SMITH, H. F., 1936  A discriminant function for plant selection. Ann. Eugen. 7:240-250.

SOLLER, M. and J. S. BECKMANN, 1983  Genetic polymorphisms in varietal identification and genetic improvement. Theor. Appl. Genet. 67:25-33.

SORENSEN, D., and D. GIANOLA, 2002 Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics. Springer-Verlag, New York.

SORENSEN, D., C. S. WANG, J. JENSEN, and D. GIANOLA, 1994  Bayesian analysis of genetic change due to selection using Gibbs sampling. Genet. Sel. Evol. 26:333-360.

STEIN, C., 1955 Inadmissibility of the usual estimator for the mean of a multivariate normal distribution, pp. 192-206 in Proceedings of the Third Berkeley Symposium on Mathematics, Statistics, and Probability. University of California Press, Berkeley.

TANNER, M. A., 1993 Tools for Statistical Inference. Springer-Verlag, New York.

THEIL, H., 1971 Principles of Econometrics. John Wiley & Sons, New York.

URIOSTE, J. I., D. GIANOLA, R. REKAYA, W. F. FIKSE, and K. A. WEIGEL, 2001  Evaluation of extent and amount of heterogeneous variance for milk yield in Uruguayan Holsteins. Anim. Sci. 72:259-268.

VERBEKE, G., and G. MOLENBERGHS, 1997 Linear Mixed Models in Practice: A SAS Oriented Approach. Springer-Verlag, New York.

VERBEKE, G., and G. MOLENBERGHS, 2000 Linear Mixed Models for Longitudinal Data. Springer, New York.

WANG, C. S., J. J. RUTLEDGE, and D. GIANOLA, 1993  Marginal inferences about variance components in a mixed linear model using Gibbs sampling. Genet. Sel. Evol. 25:41-62.

WANG, C. S., J. J. RUTLEDGE, and D. GIANOLA, 1994  Bayesian analysis of mixed linear models via Gibbs sampling with an application to litter size in Iberian pigs. Genet. Sel. Evol. 26:91-115.

WANG, T., R. L. FERNANDO, and M. GROSSMAN, 1998  Genetic evaluation by best linear unbiased prediction using marker and trait information in a multibreed population. Genetics 148:507-515.[Abstract/Free Full Text]

WEIGEL, K. A., D. GIANOLA, R. J. TEMPELMAN, C. A. MATOS, and I. H. C. CHEN et al., 1991  Improving estimates of fixed effects in a mixed linear model. J. Dairy Sci. 74:3174-3182.[Abstract]

WHITTAKER, J. C., 2001 Marker-assisted selection and introgression, pp. 673–693 in Handbook of Statistical Genetics, edited by D. J. BALDING, M. BISHOP and C. CANNINGS. John Wiley & Sons, Chichester, UK.

WHITTAKER, J. C., R. THOMPSON, and P. M. VISSCHER, 1996  On the mapping of QTL by regression of phenotype and marker-type. Heredity 77:23-32.

WHITTAKER, J. C., R. THOMPSON, and M. C. DENHAM, 2000  Marker-assisted selection using ridge regression. Genet. Res. 75:249-252.[Medline]

XU, S., 1998  Mapping quantitative trait loci using multiple families of line crosses. Genetics 148:517-524.[Abstract/Free Full Text]

ZELLNER, A., and W. VANDAELE, 1975 Bayes-Stein estimators for k-means, regression and simultaneous equation models, pp. 627–653 in Studies in Bayesian Econometrics and Statistics, edited by S. FIENBERG and A. ZELLNER. North-Holland, Amsterdam.

ZENG, Z-B., 1993  Theoretical basis for separation of multiple linked gene effects in mapping quantitative trait loci. Proc. Natl. Acad. Sci. USA 90:10972-10976.[Abstract/Free Full Text]

This article has been cited by other articles:

				H. P. Piepho Ridge Regression and Extensions for Genomewide Selection in Maize Crop Sci., June 26, 2009; 49(4): 1165 - 1176. [Abstract] [Full Text] [PDF]

				G. de los Campos, H. Naya, D. Gianola, J. Crossa, A. Legarra, E. Manfredi, K. Weigel, and J. M. Cotes Predicting Quantitative Traits With Regression Models for Dense Molecular Markers and Pedigree Genetics, May 1, 2009; 182(1): 375 - 385. [Abstract] [Full Text] [PDF]

				N. Yi and S. Xu Bayesian LASSO for Quantitative Trait Loci Mapping Genetics, June 1, 2008; 179(2): 1045 - 1055. [Abstract] [Full Text] [PDF]

				D. Gianola and J. B. C. H. M. van Kaam Reproducing Kernel Hilbert Spaces Regression Methods for Genomic Assisted Prediction of Quantitative Traits Genetics, April 1, 2008; 178(4): 2289 - 2303. [Abstract] [Full Text] [PDF]

				O. Gonzalez-Recio, D. Gianola, N. Long, K. A. Weigel, G. J. M. Rosa, and S. Avendano Nonparametric Methods for Incorporating Genomic Information Into Genetic Evaluations: An Application to Mortality in Broilers Genetics, April 1, 2008; 178(4): 2305 - 2313. [Abstract] [Full Text] [PDF]

				D. Gianola, R. L. Fernando, and A. Stella Genomic-Assisted Prediction of Genetic Value With Semiparametric Procedures Genetics, July 1, 2006; 173(3): 1761 - 1776. [Abstract] [Full Text] [PDF]

				H. Wang, Y.-M. Zhang, X. Li, G. L. Masinde, S. Mohan, D. J. Baylink, and S. Xu Bayesian Shrinkage Estimation of Quantitative Trait Loci Parameters Genetics, May 1, 2005; 170(1): 465 - 480. [Abstract] [Full Text] [PDF]

				C. Cervantes-Martinez and J. S. Brown A Haplotype-Based Method for QTL Mapping of F1 Populations in Outbred Plant Species Crop Sci., September 1, 2004; 44(5): 1572 - 1583. [Abstract] [Full Text] [PDF]

				M. Perez-Enciso, M. A. Toro, M. Tenenhaus, and D. Gianola Combining Gene Expression and Molecular Marker Information for Mapping Complex Trait Genes: A Simulation Study Genetics, August 1, 2003; 164(4): 1597 - 1606. [Abstract] [Full Text] [PDF]