HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Su, G. Articles by Jensen, J. Search for Related Content PubMed PubMed Citation Articles by Su, G. Articles by Jensen, J.

Bayesian analysis of the linear reaction norm model with unknown covariates¹

Danish Institute of Agricultural Sciences, Department of Genetics and Biotechnology, DK-8830, Tjele, Denmark

	Abstract

Top Abstract INTRODUCTION MODEL AND METHODS RESULTS DISCUSSION IMPLICATIONS LITERATURE CITED

 
The reaction norm model is becoming a popular approach for the analysis of genotype x environment interactions. In a classical reaction norm model, the expression of a genotype in different environments is described as a linear function (a reaction norm) of an environmental gradient or value. An environmental value is typically defined as the mean performance of all genotypes in the environment, which is usually unknown. One approximation is to estimate the mean phenotypic performance in each environment and then treat these estimates as known covariates in the model. However, a more satisfactory alternative is to infer environmental values simultaneously with the other parameters of the model. This study describes a method and its Bayesian Markov Chain Monte Carlo implementation that makes this possible. Frequentist properties of the proposed method are tested in a simulation study. Estimates of parameters of interest agree well with the true values. Further, inferences about genetic parameters from the proposed method are similar to those derived from a reaction norm model using true environmental values. On the other hand, using phenotypic means as proxies for environmental values results in poor inferences.

Key Words: environmental sensitivity • environmental value • genotype x environment interaction • genetic parameter • Gibbs sampler • reaction norm model

	INTRODUCTION

Top Abstract INTRODUCTION MODEL AND METHODS RESULTS DISCUSSION IMPLICATIONS LITERATURE CITED

 
The reaction norm model (Falconer and Mackay, 1996) is attractive to describe genotype x environment interactions (G x E) partly because it can accommodate a very large number of environmental levels with few parameters. In its standard version, it requires that covariates are known (Karan et al., 1999; Ravagnolo and Misztal, 2000; Kingsolver et al., 2004). However, in animal breeding applications, one may postulate a linear relationship between the phenotypic expression of a given genotype and a particular environmental effect (e.g., herd effect). In this setup the covariate (i.e., herd effect) is unknown. One approximation reported in the literature is to compute the mean phenotypic performance in the appropriate environment and to use such an estimate in lieu of the unknown covariate in the model (Calus et al., 2002; Kolmodin et al., 2002; Calus and Veerkamp, 2003).

Including a function of the data as a covariable in the sampling model for the data is clearly unsatisfactory. Apart from the understatement of uncertainty due to treating phenotypic means as known parameters, one can imagine situations that would result in misleading representations of environmental values using this approach. An example would be the presence of a genetic trend. Because in the reaction norm model a breeding value is defined as a function of the environmental gradient, inappropriate inferences about environmental values may result in incorrect ranking based on predicted genetic values.

It is therefore important to find more adequate methods to account for unknown covariates in a reaction norm model. An appealing alternative is to infer environmental values simultaneously with the other parameters of the model. The objectives of this study were 1) to describe a method and its Bayesian Markov Chain Monte Carlo (MCMC) implementation that makes this possible, and 2) using a simulation study, to test the expectation that the proposed method leads to more satisfactory inferences about genetic parameters than the approximate method mentioned previously.

	MODEL AND METHODS

Top Abstract INTRODUCTION MODEL AND METHODS RESULTS DISCUSSION IMPLICATIONS LITERATURE CITED

 
Model
When genetic and nongenetic environmental sensitivities are taken into consideration, a reaction norm model can be written as

(1)

where y is the data vector (order n), b is the vector of fixed effects (order n_b), h is the vector of environmental values (order n_h), u₀ is the vector of intercepts (order n_u), u_h is the vector of slopes (order n_u) of reaction norms for nongenetic random effects (e.g., permanent effects), a₀ is the vector of intercepts (order n_g), a_h is the vector of slopes (order n_g) of additive genetic reaction norms, and e is the vector of residual effects (order n). X, E, Z_u, H_u, Z_a, and H_a are incidence matrices. When the covariate associated with the reaction norm is treated as unknown, the row k of matrices H_u and H_a has exactly one element equal to the effect of the environment (h_i or a function of h_i) where the observation is recorded, and the others equal to zero.

In principle h can be treated as a fixed or a random vector. Here it is treated as random to better meet identifiability requirements. In the present model identifiability is a complex topic. We limit ourselves to making the statement that the functions of the parameters that are estimated and reported later are identifiable.

The conditional distribution of y is assumed to be normal having the form

where R is the matrix (order n) of random residual covariances. Without loss of generality, it is assumed that residuals are homoscedastic and independent of each other, so that R = I²_e, where I is the identity matrix and ²_e is the residual variance.

Prior Distribution of Location Parameters
The prior distribution of vector b is assumed to be improper uniform. The random vectors h, (u₀', u_h')', and (a₀', a_h')' are assumed to have normal, mutually independent prior distributions.

The prior distributions of the variance of h_i (²_h) and the residual variance (²_e) are assumed to be scaled inverse ² distributions. The (co) variances of u_0i and and the (co) variances of a_0i and are assumed to follow inverse Wishart distributions.

Joint Posterior Distribution of All the Parameters
Let be the vector of all location parameters except h; i.e., = (b', u₀', u_h', a₀', a_h')'. Then the joint posterior distribution of all the parameters is

(2)

Fully Conditional Posterior Distribution of the Location Parameters
The fully conditional posterior distribution of can be directly derived from (2) by extracting terms involving . This results in

(3)

Further, assuming h is known, define

(4)

Because p(y | , h, ²_e) = p(y | , h, ²_e), the fully conditional posterior distribution of is

(5)

Using results in Lindley and Smith (1972), Gianola and Fernando (1986), and Sorensen and Gianola (2002), it is easy to show that the posterior distribution of location parameters (), given dispersion parameters and h, is multivariate normal. Now, write the mixed model equations associated with (4) as C = r. Then | U₀, G₀, ²_e, h, y is normally distributed with mean equal to and variance equal to C^–1, i.e.,

(6)

Let _i denote an arbitrary element (or set of elements) of , and let _–i denote the vector with _i excluded. From standard multivariate normal theory, it can readily be established that, if the distribution of | U₀, G₀, ²_e, h, y is normal, the distribution of _i | _–i, U₀, G₀, ²_e, h, y is

(7)

Fully Conditional Posterior Distribution of h
From (2), the density of the fully conditional posterior distribution of h is

(8)

Based on (1), an observation y can be described as

Therefore, an alternative formulation of the reaction-norm model (1) is

(9)

where E* is the coefficient matrix obtained by replacing the nonzero element in the row k of matrix E with (1 + z_uk'u_h + z_ak'a_h).

Assuming is known, define

(10)

Because p(y | , h, ²_e) = p(y_h | , h, ²_e), the fully conditional posterior distribution of h is

(11)

Now, write the mixed model equations associated with (10) as C_hĥ = r_h. Then

(12)

and for the element i, the fully conditional posterior distribution is

(13)

Fully Conditional Posterior Distribution of Dispersion Parameters
The fully conditional posterior distribution of dispersion parameters is deduced from (2). Let H be the vector of all of the location parameters, and let W = (X: E: Z_u:H_u: Z_a: H_a). For the residual variance one obtains

(14)

which is recognized as a scaled inverse ² distribution with degrees of freedom _e + n and scale parameter [(y – W) (y – W) + _es²_e]/(_e + n), where _e is the degrees of freedom for the prior distribution of ²_e, and n is the number of observations.

The fully conditional posterior distribution of the variance of environmental values is

(15)

which is again recognized as a scaled inverse ² distribution with degrees of freedom _h + n_h and scale parameter (h'h + _hs²_h)/(_h + n_h), where _h is the degrees of freedom for the prior distribution of ²_h, and n_h is the order of h.

The fully conditional posterior distribution of the covariance matrix of the reaction norm of the nongenetic random effects is

(16)

This is an inverse Wishart distribution of dimension k_u = 2, with degrees of freedom _u + n_u and scale matrix (S²_u + V^–1_u)^–1, where V_u is the scale matrix, _u is the degrees of freedom for the prior distribution of U₀, and n_u is the order of u₀ or u_h.

The fully conditional posterior distribution of the covariance matrix of the additive genetic reaction norm is

(17)

which is an inverse Wishart distribution of dimension k_g = 2 with degrees of freedom _g + n_g and scale matrix (S²_g + V^–1_g)^–1, where V_g is the scale matrix, _g is the degrees of freedom for the prior distribution of G₀, and n_g is the order of a₀ or a_h.

Implementation of the Gibbs Sampler
The Gibbs sampler is a Monte Carlo method for obtaining samples from joint or marginal posterior distributions of all parameters in the model, by repeated sampling from fully conditional posterior distributions. The algorithm for the proposed model is as follows:

1. Compute y_h, E*, C_h, and r_h and draw h_i from N(C^–1_h(i,i)(r_hi – C_h(i,–i)h_–i), C^–1_h(i,i)).
   
2. Compute y, C, and r and draw _i from N(C^–1_{(i, i)}(r_i – C_{(i, –i)–i}),C^–1_{(i, i)}).
   
3. Sample a new from Inv – X²(_h + n_h, (h'h + _hs²_h)/(_h + n_h)).
   
4. Sample a new U₀ from Inv – W₂((S²_u + V^–1_u)^–1, _u + n_u).
   
5. Sample a new G₀ from Inv – W₂((S²_g + V^–1_g)^–1, _g + n_g).
   
6. Sample a new ²_e from Inv – X²[_e + n, ((y – W)'(y – W) + _es²_e)/(_e + n)].
   
7. Repeat (1) through (6) until enough samples are available.
   

The Gibbs sampling algorithm was implemented using the DMU package (Madsen and Jensen, 2004). This package uses the iteration on data technique to avoid storing Cand C_h.

Simulation Studies
Data Generation.
The proposed method was evaluated using simulated data. Observations were generated with the model y = 1µ + Eh + Za₀ + Ha_h + e, where h was the vector of environmental values (herd-year effects), a₀ was the vector of levels, a_h was the vector of slopes of additive genetic reaction norms, and e was the vector of random residuals. Vectors h, (a₀', a_h')' and e were assumed to be mutually independent and were sampled from

To simplify the simulation procedure, but without loss of generality, data were generated according to a structure mimicking swine. Thus, 5 generations (5 yr) of data were simulated and distributed over 50 herds. In each generation, 50 sires were mated to 1,000 dams, and each dam produced 5 offspring with records. Sires and dams were chosen randomly. Sires were used across herds, and each sire was mated to 20 dams from 5 herds. Dams were used within herds. Consequently, there were 100 individuals from 5 sires and 20 dams in each herd each generation.

The parameters used in the simulation were: ²_h = 80, ²_a0 = 100, ²_ah = 1, r_a0,ah = 0.5, and ²_e = 300. This corresponds to a G x E variance as follows: Var(a_hh) = ²_ah·²_h = 80; and a marginal variance of a datum (phenotypic variance across herd-years) as follows: ²_P = ²_h + + ²_ah²_h + ²_e = 560.

Statistical Analysis.
The simulated data were analyzed using the following models:

The model with unknown covariate of the reaction norm, treating herd-years as random effects (the proposed approach),

(M1)

the model using true herd-year effects as covariate (H^t) of the reaction norm and including herd-years as fixed effects,

(M2)

the model using phenotypic means of herd-years as proxies for the unknown covariates (H^m) of the reaction norm and including herd-years as fixed effects,

(M3)

Note that in models M2 and M3, the covariates of the reaction norm (H^t, H^m) are not necessarily equivalent to the corresponding elements of h, whereas in M1, the nonzero elements of H are equivalent to those of h.

The additive genetic variance () and heritability (h²_a) in a particular herd-year were calculated as ²_a = ²_a0 + ²_ahh² + 2_a0ahh and Because the covariate features in the additive genetic and phenotypic variances, for ease of comparison of heritabilities among methods (Figure 1) the covariate was expressed in units of the appropriate standard deviation (h* = h/_h). On defining ²_ah* = ²_ah²_h and _a0ah* = _a0ahh, then, ²_a = ²_a0 + ²_ahh² + 2_a0ahh = ²_a0 + ²_ah*h^*2 + 2_a0ah*h^*, where, using M1, ²_h was the empirical variance of the estimated herd-year effect; using M2, ²_h was the variance of true herd-year effects; and using M3, ²_h was the variance of herd-year averages.

View larger version (17K): [in this window] [in a new window]   Figure 1. Heritability as a function of herd-year value (covariate of the reaction norm expressed in units of standard deviation). True = the simulated heritability; M1 = the heritability obtained from the proposed method; M2 = the heritability obtained using true herd-year effects as covariates of the reaction norm; M3 = the heritability obtained using phenotypic means of herd-years as covariates of the reaction norm. Heritability was calculated as the average of the posterior means over 20 replicates.

	RESULTS

Top Abstract INTRODUCTION MODEL AND METHODS RESULTS DISCUSSION IMPLICATIONS LITERATURE CITED

As shown in Table 1, the proposed method (M1) yielded estimates (based on posterior means) of variance components with no detectable bias, whereas using herd-year averages as proxies for herd-year effects (M3) resulted in biased estimates. Averaged over the 20 replicates, the variance components estimated from the proposed method and from the model using true herd-year effects as covariates in the reaction norm (M2) resulted in similar inferences. These estimates agreed well with the true values. On the other hand, using herd-year averages as covariates in the reaction norm caused an overestimation of the variance component associated with level (²_a0) and an underestimation of the variance components associated with the slope (²_ah, and _a0,ah). These biases were statistically significant. After scaling herd-year averages using the ratio of the standard deviation of herd-year averages to the standard deviation of true herd-year effects, the estimates were 0.78 for ²_ah and 4.28 for _a0,ah, whereas the realized values in the simulated data were equal to 1.01 and 5.11, respectively. The sampling standard deviation of the estimates of ²_ah was largest using M1, lowest using M2, and intermediate using M3, whereas the standard deviation of the estimates of ²_a0 was largest using M3. Mean squared errors favored M1 over M3 in all cases.

	DISCUSSION

Top Abstract INTRODUCTION MODEL AND METHODS RESULTS DISCUSSION IMPLICATIONS LITERATURE CITED

 
In the present work we describe a method to infer unknown environmental values simultaneously with other parameters in a reaction norm model. Using computer simulation, this method is compared with an approximation traditionally implemented in the literature, whereby the unknown environmental value is replaced by the average of the observations in the appropriate environment. It is shown that the proposed method leads to better inferences than those derived from the approximate method.

The approximate procedure has some shortcomings. The variance among phenotypic means of production environments includes a genetic component. This results in an overestimation of the variation of environmental values. Even in a random mating population, as is the case in the current study, the variance among phenotypic means was 35% larger than the variance of true herd-year effects. Further, the correlation between herd-year means and true herd-year effects was 0.901, obviously lower than the correlation of 0.970 obtained using the method proposed in this paper. Because in the reaction norm model a breeding value is defined as a function of the environmental gradient, improper estimates of environmental values may result in incorrect ranking based on predicted genetic values.

The approximated method also results in biased estimation of variance components. Thus, the variance component associated with the slope (²_ah) was underestimated by 42%, and that associated with level (²_a0) was overestimated by 11%. The underestimation of the variance of the slope (²_ah) and the covariance between the level and the slope (_a0,ah) was partly due to the fact that the variation of herd-year averages was larger than that of true herd-year effects. After scaling herd-year averages using the ratio of the standard deviation of herd-year averages to the standard deviation of true herd-year effects, the estimates were 0.78 for ²_ah and 4.28 for _a0,ah, whereas the realized values in the simulated data were equal to 1.01 and 5.11, respectively. Thus, the bias persisted despite the fact that the scaling involved the use of the standard deviation of true herd-year effects. Underestimation of the variance of the slope understates the importance of G x E interaction. Further, inadequate inferences about production environmental effects and the resulting biased estimates of other genetic parameters reduce the accuracy of prediction of breeding values. This would be expected to have an unfavorable effect on genetic progress by selection.

The amount and sign of the bias associated with the approximate method depend on the data structure. To illustrate, an additional study was carried out with data simulated from the same sampling model as reported earlier, but with the difference that from generation 1 onward individuals were selected on the basis of their predicted additive genetic values for level. The results showed that the correlation between herd-year averages and true herd-year effects was approximately 0.80 and the variance of herd-year averages was approximately 5 times larger than the variance of the true herd-year effects. Using the herd-year average as a covariate of the reaction norm, ²_a0 was overestimated by 50%, whereas ²_ah was underestimated by 88%.

Many approximations and ad-hoc procedures have been reported in previous studies to account for unknown covariates in reaction norm models. In a study of production and fertility traits in dairy cattle, Kolmodin et al. (2002) estimated herd-year values using herd-year means computed from data that had been preadjusted for fixed effects other than herd-years. In addition, herd-year values were estimated using herd-year means that were computed from data including animals with records in the appropriate herd-year, whereas dispersion parameters and breeding values were inferred from data that only included individuals whose sires were to be evaluated. The adequacy of this approximation could not be tested because it was applied using real (as opposed to simulated) data. Kolmodin et al. (2002) made a plea in their conclusion for the development of alternative procedures that avoid using functions of the data in the sampling model for the data.

Calus et al. (2004) proposed to estimate environmental values via an iterative procedure whereby the estimated environmental effect in a given iteration replaces the value of the covariate in the next. Using computer simulation, the authors observed a negligible reduction in bias of estimates of variance components using this approach when compared with the standard procedure of replacing covariates with phenotypic averages. They suggested replacing environmental values with estimates of herd effects obtained from a large number of animals per herd, instead of from herd-years, at the cost of losing information on G x E interaction.

The overall picture that emerges is that the conventional approximations do not always produce reliable results, and it is difficult to decide a priori how they will behave in any given data set/modeling scenario. In contrast, the method that we propose here avoids ad-hoc constructs, is theoretically coherent, easy to implement, and leads to adequate inferences. An important caveat associated with the reaction norm model with unknown covariates is that of identifiability of parameters in the likelihood. This is a technically elaborate problem that is presently under investigation and hopefully will be reported elsewhere.

	IMPLICATIONS

Top Abstract INTRODUCTION MODEL AND METHODS RESULTS DISCUSSION IMPLICATIONS LITERATURE CITED

 
The reaction norm model is becoming a popular approach for the analysis of genotype x environment interactions because it can deal with a very large number of environmental levels with few parameters. A substantial variance of the slope indicates genotype x environment interaction. A correlation between the level and the slope far from plus 1 will cause reranking of animals in different environments. Unknown effects of environments (environmental values) are commonly used as an environmental gradient in reaction norm models. Using phenotypic means of appropriate environments as proxies for environmental values could lead to biased estimates of genetic parameters and breeding values. The method proposed here estimates environmental values simultaneously with the other parameters in the reaction norm model. The method was tested using simulated data and was shown to lead to estimates of parameters with no detectable bias and with smaller mean squared errors than those obtained using the conventional approximations.

	Footnotes

 
¹ The research was supported by the Research Project DARCOF II: Research in organic farming 2000–2005—Organic dairy production systems.

² Corresponding author: guosheng.su{at}agrsci.dk

Received for publication September 14, 2005. Accepted for publication January 31, 2006.

	LITERATURE CITED

Top Abstract INTRODUCTION MODEL AND METHODS RESULTS DISCUSSION IMPLICATIONS LITERATURE CITED

 


Calus, M. P. L., A. F. Groen, and G. de Jong. 2002. Genotype x environment interaction for protein yield in Dutch dairy cattle as quantified by different models. J. Dairy Sci. 85:3115–3123.[Abstract/Free Full Text]

Calus, M. P. L., and R. F. Veerkamp. 2003. Estimation of environmental sensitivity of genetic merit for milk production traits using a random regression model. J. Dairy Sci. 86:3756–3764.[Abstract/Free Full Text]

Calus, M. P. L., P. Bijma, and R. F. Veerkamp. 2004. Effects of data structure on the estimation of covariance functions to describe genotype by environment interaction in a reaction norm model. Genet. Sel. Evol. 36:489–507.[CrossRef][Medline]

Falconer, D. S., and T. F. C. Mackay. 1996. Introduction to Quantitative Genetics, 4th ed. Longman Group, Essex, UK.

Gianola, D., and R. F. Fernando. 1986. Bayesian methods in animal breeding. J. Anim. Sci. 63:217–244.[Abstract/Free Full Text]

Karan, D., B. Moreteau, and J. R. David. 1999. Growth temperature and reaction norms of morphometrical traits in a tropical drosophilid: Zaprionus indianus. Heredity 83:398–407.[Medline]

Kingsolver, J. R., G. J. Ragland, and J. G. Shlichta. 2004. Quantitative genetics of continuous reaction norms: Thermal sensitivity of caterpillar growth rates. Evolution 58:1521–1529.[CrossRef][Medline]

Kolmodin, R., E. Strandberg, P. Madsen, J. Jensen, and H. Jorjani. 2002. Genotype by environment interaction in Nordic dairy cattle studied using reaction norms. Acta Agric. Scand., Sect. Anim. Sci. 52:11–24.

Lindley, D. V., and A. F. M. Smith. 1972. Bayes estimates for the linear model. J. R. Stat. Soc. Sect. B 34:1–41.

Madsen, P., and J. Jensen. 2004. A User’s Guide to DMU, version 6, release 4.5. Danish Inst. Agric. Sci., Tjele, Denmark.

Ravagnolo, O., and I. Misztal. 2000. Genetic components of heat stress in dairy cattle, parameter estimation. J. Dairy Sci. 83:2126–2130.[Abstract]

Sorensen, D., and D. Gianola. 2002. Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics. Springer-Verlag, New York, NY.

This article has been cited by other articles:

				M. M. Shariati, G. Su, P. Madsen, and D. Sorensen Analysis of Milk Production Traits in Early Lactation Using a Reaction Norm Model with Unknown Covariates J Dairy Sci, December 1, 2007; 90(12): 5759 - 5766. [Abstract] [Full Text] [PDF]

				M. Piles, M. Garcia-Tomas, O. Rafel, J. Ramon, N. Ibanez-Escriche, and L. Varona Individual efficiency for the use of feed resources in rabbits J Anim Sci, November 1, 2007; 85(11): 2846 - 2853. [Abstract] [Full Text] [PDF]

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Su, G. Articles by Jensen, J. Search for Related Content PubMed PubMed Citation Articles by Su, G. Articles by Jensen, J.