Deregressing estimated breeding values and weighting information for genomic regression analyses
- Dorian J Garrick^{1, 2}Email author,
- Jeremy F Taylor^{3} and
- Rohan L Fernando^{1}
https://doi.org/10.1186/1297-9686-41-55
© Garrick et al; licensee BioMed Central Ltd. 2009
Received: 2 July 2009
Accepted: 31 December 2009
Published: 31 December 2009
Abstract
Background
Genomic prediction of breeding values involves a so-called training analysis that predicts the influence of small genomic regions by regression of observed information on marker genotypes for a given population of individuals. Available observations may take the form of individual phenotypes, repeated observations, records on close family members such as progeny, estimated breeding values (EBV) or their deregressed counterparts from genetic evaluations. The literature indicates that researchers are inconsistent in their approach to using EBV or deregressed data, and as to using the appropriate methods for weighting some data sources to account for heterogeneous variance.
Methods
A logical approach to using information for genomic prediction is introduced, which demonstrates the appropriate weights for analyzing observations with heterogeneous variance and explains the need for and the manner in which EBV should have parent average effects removed, be deregressed and weighted.
Results
An appropriate deregression for genomic regression analyses is EBV/r^{2} where EBV excludes parent information and r^{2} is the reliability of that EBV. The appropriate weights for deregressed breeding values are neither the reliability nor the prediction error variance, two alternatives that have been used in published studies, but the ratio (1 - h^{2})/[(c + (1 - r^{2})/r^{2})h^{2}] where c > 0 is the fraction of genetic variance not explained by markers.
Conclusions
Phenotypic information on some individuals and deregressed data on others can be combined in genomic analyses using appropriate weighting.
Background
Genomic prediction [1] involves the use of marker genotypes to predict the genetic merit of animals in a target population based on estimates of regression of performance on high-density marker genotypes in a training population. Training populations might involve genotyped animals with alternative types of information including single or repeated measures of individual phenotypic performance, information on progeny, estimated breeding values (EBV) from genetic evaluations, or a pooled mixture of more than one of these information sources. In pooling information of different types, it is desirable to avoid any bias introduced by pooling and to account for heterogeneous variance so that the best use is made of available information.
Uncertainty as to whether or not EBV should be used directly or deregressed or replaced by measures such as daughter yield deviation (DYD) [2], and the manner in which information should be weighted, if at all, has been apparent for some time in literature related to discovering and fine-mapping quantitative trait loci (QTL). Typically in fixed effects models with uncorrelated residuals, observations would be weighted by the inverse of their variances. Morsci et al. [3] pointed out the counter intuitive behavior of using the reciprocal of the variance of breeding values as weights in characterization of QTL and followed the arguments of Rodriguez-Zas et al. [4] in using reliability as weights. Rodriguez-Zas et al. [4] did analyses that were limited by features of the chosen software so EBV/2 (i.e. predicted transmitting ability PTA) were multiplied by the square root of reliability and analyzed unweighted. Georges et al. [5] deregressed PTA to construct DYD and weighted these using the inverse of the variance of the DYD. Spelman et al. [6] had direct access to DYD and similarly weighted these by the inverse of their scaled variance, equivalent to using the inverse of reliability as weights. Other researchers have reported the use of PTA [7], standardized PTA [7, 8] or DYD weighted by respective reliabilities [8]. The uncertainty associated with using information for QTL discovery has recently been extended to genomic prediction. An Interbull survey [9] of methods being used in various countries for genomic prediction of dairy cattle reported that some researchers used deregressed proofs weighted with corresponding reliabilities, others used DYD weighted by effective daughter contributions, while yet others used EBV without any weighting. The objective of this paper is to present a logical argument for using deregressed information, appropriately weighted for analysis. For simplicity, we consider the residual variance from the perspective of an additive model but the deregression and weighting concepts extend to analyses that include dominance and epistasis.
Methods
An ideal model
where g is a vector of true genetic merit (i.e. breeding value BV) with var(g) = T , the scalar is the genetic variance and T can be constructed using the theory from combined linkage disequilibrium and linkage analyses [10], μ is an intercept, M is an incidence matrix whose columns are covariates for substitution, genotypic or haplotypic effects, a are effects to be estimated, var(Ma) = , G is a genomic relationship matrix [11–13], ε is the lack of fit, var(ε) = , hopefully small and will be 0 if BV could be perfectly estimated as a linear function of observed marker genotypes. In different settings, a might be defined as a vector of fixed effects [14] or a vector of random effects [1]. Even when a is fixed, Ma is random because M, which contains genotypes, is random. However, in genomic analyses M is treated as fixed because the analysis is conditional on the observed genotypes. The philosophical issues related to the randomness of M and a are discussed in detail by Gianola [15] but for our context it is sufficient to define var(Ma) = without explicitly specifying distributional properties of M or a.
Genotypes used as covariates in Ma are unlikely to capture all the variation in true genetic merit, either because they are not comprehensively covering the entire genome, or because linkage disequilibrium between markers and causal genes is not perfect. Knowledge of E is required in the analysis whether a is treated as a fixed (e.g. GLS) or random effect (e.g. BLUP). In practice with experiments that involve related animals, it is unreasonable to assume E has a simple form such as a diagonal matrix since that implies a zero covariance between lack of fit effects for different animals, however, it can be approximated using knowledge on the pedigree using the additive relationship matrix, A [16]. These lack of fit covariances can be accommodated by fitting a polygenic effect for each animal, in addition to the marker genotypes [17], or accounted for by explicitly modeling correlated residuals. For a non-inbred animal, , therefore and the proportion of the genetic variance not accounted for by the markers can be defined to be . The scalar c, will be close to 0 if markers account for most of the genetic variation and close to 1 if markers perform poorly.
A model using individual phenotypic records
with since cov(ε, e') = 0. This model can be fitted by explicitly including a random polygenic effect for ε, or by accounting for the non-diagonal variance-covariance structure of the residuals defined as var (ε + e). Including a polygenic term is not typically done in genomic prediction analyses [12, 18], and when undertaken does not seem to markedly alter the accuracy of genomic predictions [Habier D. Personal communication]. Assuming var (ε + e) is a scaled identity matrix facilitates the computing involved in fitting this model, as the relevant mixed model equations can be modified by multiplying the left- and right-hand sides by the unknown scale parameter as is typically done in single trait analyses. However, this is not an option if residuals are heterogeneous, for example, because they involve varying numbers of repeated observations.
A model using repeated records on the individual
which can be used for weighted regression analyses treating marker effects as fixed or random. When c = 0, the genetic effects can be perfectly explained by the model, and for n = 1, a single observation on the individual, the weight is 1 for any heritability. Scaling the weights is convenient because records with high information exceed 1 and the weights are trait independent which is useful when analysing multiple traits with identical heritability and information content.
Offspring averages as data
This expression can be used as weights in the fixed or random regression of full-sib progeny means on parent average marker genotypes.
Estimated breeding values as training data
There are at least two issues with this formulation of the problem, which may not be immediately apparent, and which both result from properties of BLUP. The first issue is that the addition of the prediction error term to the left- and right-hand side of (8) actually reduces rather than increases the variance, despite the fact that diagonal elements of must exceed 0, in contrast to the addition of non-genetic random residual effects in (3). That is , whereas var(g_{ i }) < var(y_{ i }), due to shrinkage properties of BLUP estimators [19]. Generally, but for BLUP so that implying . The reduction in variance of the training data comes about because prediction errors are negatively correlated with BV as can be readily shown since . This means that superior animals tend to be underevaluated (i.e. have negative prediction errors) whereas inferior animals tend to be overevaluated. This is a consequence of shrinkage estimation and prediction errors being uncorrelated with EBV, i.e. . In order to account for the covariance between the prediction errors and the BV, a model that accounted for such covariance would need to be fitted. Such models are computationally more demanding compared to models whereby the fitted effects and residuals are uncorrelated. The second issue resulting from the properties of BLUP, is that it is a shrinkage estimator, that shrinks observations towards the mean, the extent of shrinkage depending upon the amount of information. This is apparent if one considers the regression of phenotype on true genotype (i.e. BV) which is 1, whereas the regression of EBV on BV is equal to ≤ 1, where is the reliability of the EBV (for animal i) or squared correlation between BV and EBV. In the context of any marker locus, the contrast in EBV between genotypes at a particular locus is shrunk relative to the contrast that would be obtained if BV or phenotypes were used as data, with the shrinkage varying according to . We are, however, interested in estimating the effect of a marker on phenotype, but we get a lower value for the contrast if EBV with ≤ 1 are used as data, rather than using phenotypes. A further complication is that training data based on EBV typically comprise individuals with varying . This problem can be avoided by deregressing or unshrinking the EBV.
Deregressing estimated breeding values
for some matrix K chosen so that and is a constant. Since then this expression will be 0 when . For this value k_{ i }, , a constant for all animals regardless of their reliability. Accordingly, the deregression matrix is K = diagonal and the deregressed observations are . Note in passing that the nature of the deregression will depend upon the EBV base. Genetic evaluations are typically adjusted to a common base before publication, by addition or subtraction of some constant. The EBV should be deregressed after removing the post-analysis base adjustment or by explicitly accounting for the base in the deregression procedure [20]. To show the dependence of the deregression to the post-analysis base, supposes that EBV are adjusted to a base, b. Then a linear contrast in deregressed EBV without removing the base effect is unless . Marker effects are typically estimated as linear combinations of data, and will therefore be sensitive to the base adjustment.
A deregressed observation represents a single value that encapsulates all the information available on the individual and its relatives, as if it was a single observation with h^{2} = r^{2}. This can be shown by recognising that h^{2} is the regression of genotype on phenotype. Taking the deregressed observation to be the phenotype, . Training on deregressed EBV is therefore like training on phenotypes with varying h^{2}. Provided > h^{2}, training on deregressed EBV is equivalent to having a trait with higher heritability. However, as explained later, we recommend removing ancestral information from the deregressed EBV.
Weighting deregressed information
an expression analogous to (5) with n = 1 and h^{2} = . Note that the weight in (10) approaches as →1 in which case the weight tends to infinity as c→0. This is the same as would occur when the number of offspring p→∞, and p is used as a weight.
Removing parent average effects
then using the facts [19] that and leads to , and . Rearranging these equations, , and . The formula to derive the inverse of a 2 × 2 matrix applied to the coefficient matrix from (11) gives , and for .
Application of (15) provides the solution for that can be substituted in (14) to solve for , together enabling reconstruction of the coefficient matrix of (11).
Third, the right-hand side of (11) can be formed by multiplying the now known coefficient matrix by the known vector of EBV for PA and individual. The right-hand side on the individual, free of PA effects is The equation to obtain an estimate of EBV for animal i, free of its parent average, , based only on , is and the corresponding for use in constructing the weights in (10) is given by . The deregressed information is , which simplifies to and is analogous to an average. An iterative procedure using mixed model equations to simultaneously deregress all the sires in a pedigree, while jointly estimating the base adjustment and accounting for group effects was given by Jairath et al [20]. However, that method requires knowledge on the numbers of offspring of each sire.
Double counting of information from descendants
Genetic evaluation of animal populations results in EBV that are a weighted function of the parent average EBV, any information on the individual, adjusted for fixed effects, and a weighted function of the EBV of offspring, adjusted for the merit of the mates [2]. The previous section has argued for the removal of parent average effects in constructing information for genomic analyses. It could be argued that information from genotyped descendants should also be removed to avoid double counting. This can be achieved during the evaluation process, and is desirable in the absence of selection. If the genotyped descendants are a selected subset, the removal of their information will lead to biased information on the individual. Simulation suggests that the double counting of descendants performance has negligible impact on genomic predictions (results not shown).
Results
Weights for different information sources
Relative weights^{ a }for n phenotypic observations on the individual, p observations in twice the halfsib progeny mean with heritability 0.25 and repeatability 0.6, or deregressed EBV with reliability r^{2} for varying values of c, the proportion of genetic variation for which genotypes cannot account
c | |||||
---|---|---|---|---|---|
Information Source | 0.8 | 0.5 | 0.25 | 0.1 | |
Mean of n repeated records | n | ||||
1 | 0.79 | 0.86 | 0.92 | 0.97 | |
2 | 1.00 | 1.11 | 1.22 | 1.30 | |
5 | 1.19 | 1.35 | 1.52 | 1.65 | |
10 | 1.27 | 1.46 | 1.66 | 1.81 | |
2 × mean of p half-sib offspring | p | ||||
5 | 0.79 | 0.86 | 0.92 | 0.97 | |
10 | 1.30 | 1.50 | 1.71 | 1.88 | |
20 | 1.94 | 2.40 | 3.00 | 3.53 | |
Deregressed EBV with reliability r^{2} | r ^{2} | ||||
0.1 | 0.31 | 0.32 | 0.32 | 0.33 | |
0.2 | 0.63 | 0.67 | 0.71 | 0.73 | |
0.3 | 0.96 | 1.06 | 1.16 | 1.23 | |
0.4 | 1.30 | 1.50 | 1.71 | 1.88 | |
0.5 | 1.67 | 2.00 | 2.40 | 2.73 | |
0.6 | 2.05 | 2.57 | 3.27 | 3.91 | |
0.7 | 2.44 | 3.23 | 4.42 | 5.68 | |
0.8 | 2.86 | 4.00 | 6.00 | 8.57 | |
0.9 | 3.29 | 4.91 | 8.31 | 14.21 | |
1.0 | 3.75 | 6.00 | 12.00 | 30.00 |
Removing parent average effects
Suppose genomic training is to be undertaken for a trait using EBV available from national evaluations that have yet to be deregressed. Widely-used bulls have been genotyped and the EBV and r^{2} of those bulls are available, along with corresponding information on the sire and dam of each bull. Such a trio might have values of = 10, = 0.97; = 2, = 0.36; and = 15, = 0.68. Given h^{2} = 0.25, λ = 0.75/0.25 = 3, the PA information is , and . Using (15), with α = 5.97, δ = 0.523, then = 9.16 which substituted in (14) gives = 5.08.
Substituting these information contents into the coefficient matrix or left-hand side of (11) is with inverse . These values correspond to = 0.5 - 3 × 0.0558 = 0.33 and = 1.0 - 3 × 0.1066 = 0.68 the reported and confirming the equations used to determine the information content. The right-hand side of (11) can then be reconstructed by multiplying the coefficient matrix by the vector of EBV as . The element of interest is the right-hand side element corresponding to the individual, obtained as = -6 × 6 + 11.08 × 15 = 130. The deregressed information for use in subsequent analysis is obtained as and the corresponding reliability of this information free of PA effects is = 1.0 - 3/(5.08 + 3) = 0.63. The relevant scaled weight for use with the deregressed information on this individual assuming c = 0.5 can be found using (10) as . This implies that the deregressed information is 2.76 times more valuable than a single record on the individual.
Discussion
The relative value of alternative information sources varies according to c, the parameter that reflects the ability of the genotypic covariates to predict genetic merit. Genomic prediction models that fit well have small values for c and result in greater relative emphasis of reliable information than is the case when the genomic prediction model fits poorly and the residual variation is dominated by contributions from lack-of-fit. For example, the mean of 20 halfsib progeny has about 3.6 times the value of the mean of 5 progeny when c is 0.1, and 2.5 times the value when c is 0.8. Deregressed EBV with reliability 1.0 are 11 times as valuable as reliability 0.5 when c is 0.1 but only 3 times as valuable when c is 0.5. These results indicate that collecting genotypes and phenotypes on training animals with low to moderate reliability will be of more relative value to genomic predictions that account for only 50% genetic variation (i.e. correlation 0.7 between genomic prediction and real merit) than they will for genomic predictions that account for a high proportion of variance.
The impact of the assumed c is to influence the relative value of individuals with reliable information, such as progeny test results, in comparison to individuals with information from less reliable sources, such as individual records. The use of too large a value of c will result in overemphasis of less accurate information in relation to more accurate information. The use of too small a value of c will result in too little emphasis on less accurate records. The correct value of c will not be known prior to training analyses but can be estimated from validation analyses. Training analyses could then be repeated using the estimated value of c. Alternatively, sensitivity to c could be assessed by training using a range of values. The sensitivity to c varies according to the heterogeneity of information content in the training data.
In practice, information sources of phenotypic data on training individuals can vary more widely than the examples derived in this paper. For example, training individuals might have their own and a mix of half-and fullsib progeny observed. In such cases, a practical approach is to first set up the mixed model equations that would be appropriate to estimate breeding values on the training individuals and use these to solve for the deregressed information [2]. This approach could also be useful in circumstances where training individuals do not all have the appropriate phenotypes. Consider a situation where some individuals have carcass measurements while others have correlated observations such as live animal ultrasound measures. A bivariate analysis of these two traits could be used to produce a single deregressed value for the carcass trait for each animal that accounted for appropriately weighted ultrasound information.
Conclusions
The arguments put forward in this manuscript support the use of deregressed information, in agreement with practices adopted by many researchers [22]. The weighting factors proposed in this paper differ from any reported in the literature except when the parameter c = 0 in which cases the weights are effectively the same as those used by Georges et al. [5] and Spelman et al. [6]. In practice, the benefit of deregression and the subsequent weighting of alternative information sources will depend on the extent to which the number of repeat records, number of progeny and/or r^{2} varies among individuals in the training population.
Declarations
Acknowledgements
DJG and RLF are supported by the United States Department of Agriculture, National Research Initiative grant USDA-NRI-2009-03924 and by Hatch and State of Iowa funds through the Iowa Agricultural and Home Economic Experiment Station, Ames, IA.
Authors’ Affiliations
References
- Meuwissen THE, Hayes BJ, Goddard ME: Prediction of total genetic value using genome-wide dense marker maps. Genetics. 2001, 157: 1819-1829.PubMed CentralPubMedGoogle Scholar
- VanRaden PM, Wiggans GR: Derivation, calculation, and use of national animal model information. J Dairy Sci. 1991, 74 (8): 2737-2746. [http://www.hubmed.org/display.cgi?uids=1918547]View ArticlePubMedGoogle Scholar
- Morsci NMTJ, Schnabel RD: Association analysis of adinopectin and somatostatin polymorphisms on BTA1 with growth and carcass traits in Angus Association analysis of adinopectin and somatostatin polymorphisms on BTA1 with growth and carcass traits in Angus cattle. Anim Genet. 2006, 37: 554-562. 10.1111/j.1365-2052.2006.01528.x.View ArticlePubMedGoogle Scholar
- Rodriguez-Zas SL, Southey BR, Heyen DW, Lewin HA: Interval and composite interval mapping of somatic cell score, yield, and components of milk in dairy cattle. J Dairy Sci. 2002, 85 (11): 3081-3091.View ArticlePubMedGoogle Scholar
- Georges M, Nielsen D, Mackinnon M, Mishra A, Okimoto R, Pasquino AT, Sargeant LS, Sorensen A, Steele MR, Zhao X: Mapping quantitative trait loci controlling milk production in dairy cattle by exploiting progeny testing. Genetics. 1995, 139 (2): 907-920.PubMed CentralPubMedGoogle Scholar
- Spelman RJ, Coppieters W, Karim L, van Arendonk JA, Bovenhuis H: Quantitative trait loci analysis for five milk production traits on chromosome six in the Dutch Holstein-Friesian population. Genetics. 1996, 144 (4): 1799-1808.PubMed CentralPubMedGoogle Scholar
- Ashwell MS, Da Y, VanRaden PM, Rexroad CE, Miller RH: Detection of putative loci affecting conformational type traits in an elite population of United States Holsteins using microsatellite markers. J Dairy Sci. 1998, 81 (4): 1120-1125.View ArticlePubMedGoogle Scholar
- Van Tassell CP, Sonstegard TS, Ashwell MS: Mapping quantitative trait loci affecting dairy conformation to chromosome 27 in two Holstein grandsire families. J Dairy Sci. 2004, 87 (2): 450-457.View ArticlePubMedGoogle Scholar
- Loberg A, Durr JW: Interbull survey on the use of genomic information. Proc Interbull Intl Workshop. 2009Google Scholar
- Meuwissen THE, Goddard ME: Prediction of identity by descent probabilities from marker-haplotyes. Genet Sel Evol. 2001, 33: 605-634. 10.1186/1297-9686-33-6-605.PubMed CentralView ArticlePubMedGoogle Scholar
- Nejati-Javaremi A, Smith C, Gibson JP: Effect of total alleleic relationship on accuracy of evaluation and response to selection. J Anim Sci. 1997, 75: 1738-1745.PubMedGoogle Scholar
- VanRaden PM: Efficient methods to compute genomic predictions. J Dairy Sci. 2008, 91 (11): 4414-4423. 10.3168/jds.2007-0980.View ArticlePubMedGoogle Scholar
- Strandén I, Garrick DJ: Technical note: Derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit. J Dairy Sci. 2009, 92 (6): 2971-2975. 10.3168/jds.2008-1929. [http://www.hubmed.org/display.cgi?uids=19448030]View ArticlePubMedGoogle Scholar
- Falconer DS, Mackay TFC: Introduction to Quantitative Genetics. 1996, New York: Longman, Inc, fourthGoogle Scholar
- Gianola D, de los Campos G, Hill WG, Manfredi E, Fernando R: Additive genetic variability and the Bayesian alphabet. Genetics. 2009, 183: 347-363. 10.1534/genetics.109.103952.PubMed CentralView ArticlePubMedGoogle Scholar
- Van Vleck LD: Selection index and introduction to mixed model methods. Boca Raton: CRC 1993 chap. Genes identical by descent - the basis of genetic likeness, 49-Google Scholar
- Calus MPL, Meuwissen THE, de Roos APW, Veerkamp RF: Accuracy of genomic selection using different methods to define haplotypes. Genetics. 2008, 178: 553-561. 10.1534/genetics.107.080838.PubMed CentralView ArticlePubMedGoogle Scholar
- Weigel KA, de los Campos G, González-Recio O, Naya H, Wu XL, Long N, Rosa GJ, Gianola D: Predictive ability of direct genomic values for lifetime net merit of Holstein sires using selected subsets of single nucleotide polymorphism markers. J Dairy Sci. 2009, 92 (10): 5248-5257. 10.3168/jds.2009-2092.View ArticlePubMedGoogle Scholar
- Henderson CR: Best linear unbiased estimation and prediction under a selection model. Biometrics. 1975, 31: 423-449. 10.2307/2529430.View ArticlePubMedGoogle Scholar
- Jairath L, Dekkers JC, Schaeffer LR, Liu Z, Burnside EB, Kolstad B: Genetic evaluation for herd life in Canada. J Dairy Sci. 1998, 81 (2): 550-562.View ArticlePubMedGoogle Scholar
- Mrode R: BLUP univariate models with one random effect. In Linear Models for the Prediction of Animal Breeding Values. 2005, Cambridge: CABIView ArticleGoogle Scholar
- Thomsen H, Reinsch N, Xu N, Looft C, Grupe S, Kuhn C, Brockmann GA, Schwerin M, Leyhe-Horn B, Hiendleder S, Erhardt G, Medjugorac I, Russ I, Forster M, Brenig B, Reinhardt F, Reents R, Blumel J, Averdunk G, Kalm E: Comparison of estimated breeding values, daughter yield deviations and de-regressed proofs within a whole genome scan for QTL. J Anim Breed Genet. 2001, 118: 357-370. 10.1046/j.1439-0388.2001.00302.x.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.