Published online 6 May 2005
Published in Crop Sci 45:1045-1051 (2005)
© 2005 Crop Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
PLANT GENETIC RESOURCES
A Threshold Model for Multiyear Genebank Data Based on Different Rating Scales
K. Hartung and
H. P. Piepho*
Institut für Pflanzenbau und Grünland (430c), Universität Hohenheim, Fruwirthstrasse 23, D-70593 Stuttgart, Germany
* Corresponding author (piepho{at}uni-hohenheim.de)
 |
ABSTRACT
|
|---|
Genebanks routinely assess plant characteristics of accessions on an ordinal rating scale. A common problem with long-term data is that the rating scale is changed from time to time. This paper proposes a method for joint analysis of data from different rating scales, assuming a threshold model with a common latent scale for the different rating systems. The method may be used to derive mean scores on any one of the rating scales based on a long-term series of evaluation trials, as is illustrated using evaluation data on barley (Hordeum spp.). While the proposed method was motivated by data problems encountered in genebanks, it may be of equal interest in other areas of application, such as plant breeding, where rating scales change across time.
Abbreviations: BLUP, best linear unbiased prediction
|
INTRODUCTION
|
|---|
GENEBANKS EVALUATE a large number of accessions each year, mostly in unreplicated field trials. Many traits (e.g., resistance to diseases) are assessed on an ordinal rating scale comprising a small number of ordered categories, typically between five and nine. For example, a rating scale for a disease may have the following categories for degree of susceptibility: 1 = no disease, 2 = low, 3 = intermediate, 4 = high, 5 = very high.
Evaluation data covering many years of field trials raise the problem of how to combine different years into a single value per accession. For metric data (yield, thousand-kernel weight, etc.), the standard approach is to use an appropriate linear model for the series of trials and to estimate least squares means per accession (Piepho, 2003a). Rating data are often analyzed using the same type of model, despite the fact that ordinal data strictly do not meet the usual assumptions of homogeneity of variance, normality, and linearity/additivity. When replicate data (on a plant basis) are available per field plot, mean scores per plot often show no significant departures from the usual assumptions (Thöni, 1985; Schumacher and Thöni, 1990). Genebanks, however, assess scores on a plot basis, so results for replicate data per plot are not directly applicable.
Genebanks are frequently faced with the problem that rating scales change across years. This complicates the integration of multiyear data into a single score per accession. As of yet, sound statistical procedures for this problem are lacking. The present paper proposes a threshold model (McCullagh and Nelder, 1989; Piepho and Kalka, 2003) for rating data applicable to a series of evaluation trials. The standard threshold model is extended to jointly analyze data from different rating scales. The fitted model may be used to compute mean rating scores per accession on any one of the employed scales. Such scores are useful for classification of accessions (Franco et al., 2003) and for database queries by users of genetic resources. We exemplify the proposed procedure using evaluation data on barley of the genebank at the IPK (Institute of Plant Genetics and Crop Plant Research), at Gatersleben, Germany. The method is compared with a standard mixed linear model analysis which treats ratings as metric data.
|
THEORY
|
|---|
Metric data for accessions evaluated in several years may be analyzed based on the two-way model (Piepho, 2003a):
 | [1] |
where yij = value of ith accession in jth year; µj = mean of jth year; gi = main effect of ith accession; and eij = residual, comprising accession-by-year interaction and error.
We propose to assume years fixed and accessions random, so the effects gi and eij are assumed to follow normal distributions with zero mean and variances 
g2 and 2, respectively. The model may be used to estimate accession values vi = 
• + gi. The random assumption for accessions allows estimation of accession values by BLUP (best linear unbiased prediction; Searle et al., 1992), which has the desirable property of shrinkage, that is, relative to the least squares means, BLUPs of vi are shrunken toward the overall mean, the shrinkage factor depending on the two variance components. Shrinkage estimators have been shown to provide better estimates than standard least squares estimates that are based on fixed effects models (Hill and Rosenberger, 1985; Piepho, 1994, 1998). One might also assume years as random, arguing that environmental conditions in a particular year are not predictable (Searle et al., 1992). Usually, however, the magnitude of year effects is so large that accession estimates are little affected by whether years are fixed or random because the between-year information is negligible (Patterson, 1997).
Ordered categorical data (rating data) may be assumed to be generated according to a threshold model, which is based on an unobserved latent metric variable (McCullagh and Nelder, 1989). Observations in the ordered categories result from a subdivision of the latent scale into c categories. Denoting the latent variable as yij and assuming that threshold values 
k (k = 1, ..., c – 1) subdivide the latent scale into c categories, observed ratings zij are generated as follows:
 | [2] |
Statistical methods have been developed to fit linear models to the continuous latent variable yij when the only data available are the categories into which the observations fell (McCullagh and Nelder, 1989). Here, we use the linear predictor in Eq. [1]. The threshold values 1 < 2 < ... < c–2 < c–1 are unknown parameters which need to be estimated along with the parameters for the linear model. Identifiability of all parameters requires two restrictions to be imposed, for example, 2 = 1 and • = 0. The need for constraints arises essentially because the latent scale is an unobserved construct. Thus, the same probabilities in the different categories can be obtained by different sets of parameter values. For example, we may shift both the thresholds and all year means in the same direction by the same amount without altering the probabilities. To resolve this indeterminacy, a restriction needs to be imposed on the location of either the thresholds or of the year means, such as • = 0. Similarly, the spread of year effects, thresholds, and random effects may be changed simultaneously without altering the probabilities. Thus, the variance on the latent scale needs to be fixed, and one way of doing this is to set 2 = 1. Assuming normality for gi and eij on the latent scale, we have that conditionally on
 | [3] |
the ordinal ratings zij follow a multinomial distribution with probabilities
 | [4] |
where 
(.) denotes the standard normal distribution function, and 0 = –
and c = +.
When, after some years, the rating system is changed, this may be accommodated by assuming two different sets of thresholds for the two periods and by assuming that, on the latent scale, the same linear model applies in both periods. This idea works even when the number of categories is not the same in each period (Piepho and Kalka, 2003). The linear model on the latent scale is
 | [5] |
where yipj is the latent value of the ith accession in the jth year within the pth period (p = 1, 2), and µpj is mean of the jth year within in the pth period. Equation [5] is essentially the same as Eq. [1], except that we have added the index p to represent the stratification of years into two periods with different rating scales. This alternative expression is convenient for specifying constraints to be imposed on the parameters, as discussed below. The kth threshold for the pth rating system is denoted as kp (k = 1, ..., cp – 1, where cp is the number of categories of the pth rating scale). For identifiability, we require 2 = 1, 1• = 0, and 2• = 0. In contrast to Eq. [1], we now require two constraints on the year means, one for each period. This is necessary because, without further assumptions on the thresholds, the indeterminacy problem occurs independently for each of the two periods. The restriction that means across years are the same for both periods is a strong assumption which cannot strictly be valid since random variation is expected among years. A more realistic model would consider the year factor as random, thus enabling the weaker assumption that µpj has the same expected value in both periods. This assumption is not made here, however, since it would render the estimation by maximum likelihood computationally very demanding. The reason is that Maximum Likelihood estimation in mixed effects threshold models requires integrating random effects out of the likelihood using numerical methods such as Gaussian quadrature. The computational demand remains relatively low if independent blocks of data (subjects) can be identified so that integration can be done per subject and the log-likelihood values per subject added subsequently. With fixed year effects, subjects correspond to accessions, and there is a single random effect per subject. By contrast, when there are random year effects and random genotype effects, there are no independent blocks of data and there is essentially just a single subject. Thus, a single integral needs to be computed across many effects rather than many small integrals across a single random effect each. The former integration task is computationally much more demanding than the latter. The model with fixed years may be seen as an approximation to the more useful random-years model. The more years there are per period, the less problematic is the restriction on the fixed year means, because similarity of period means across years will increase with the number of years.
A more favorable situation occurs when the two rating scales are anchored in at least one common threshold. In this case, the assumption of common year means in both periods is no longer necessary, and we may simplify the restriction on year effects as •• = 0. For example, the Gatersleben genebank used two different rating scales for assessing disease severity of mildew on barley in the periods from 1949 to 1991 and from 1993 to 2002, respectively (Fig. 1). On both scales, the first category is defined as no disease. In this particular case, it seems reasonable to assume that the evaluator's perception as to whether or not a plant is diseased will not be affected by the definition of the neighboring category on either scale. Thus, there is justification to assume that the first threshold is identical on both scales, as is illustrated in Fig. 1. It should be stressed, however, that identical verbal definition of a category alone is no sufficient justification for the assumption of a common threshold. For example, in the highest category, disease severity is described as very high on both scales. The neighboring categories are high in the first period and high to very high in the second period, respectively. The difference in definition of the neighboring category makes the assumption of a common threshold doubtful. For example, perception of the scale might be such that a plant appearing as very highly diseased in the first period is judged only as highly to very highly diseased in the second period due to the finer subdivision of the scale. In the mildew example, one common threshold could be assumed because both scales had the category no disease in common. More generally, whenever two scales share the category not present, the assumption of a common threshold is tenable, so that the simple restriction •• = 0 for year means can be imposed. In our experience, many rating scales share the not present category.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 1. Illustration of threshold model based on two different rating systems for mildew in barley assuming a common latent scale. Position of thresholds does not correspond exactly to parameter estimates.
|
|
Estimation of the threshold model by maximum likelihood requires an integration of the likelihood over the random effects in the linear predictor. Since the integral cannot be solved analytically, numerical methods are needed. We use adaptive Gaussian quadrature (Pinheiro and Bates, 1995) for this task, as available in the NLMIXED procedure of SAS (SAS Institute, 1999; Piepho, 2002). Gaussian quadrature essentially replaces the integral with a weighted sum. The accuracy of the approximation depends on the number of quadrature points. Adaptive quadrature selects the number of points so that the relative change between log likelihood values is small (<0.0001 by the default setting of NLMIXED). In our experience, the adaptive procedure works very well in practice. Usually we feel more comfortable, however, by checking that increasing the number of quadrature points above that determined adaptively by NLMIXED does, in fact, lead to only marginal changes in the parameter estimates. The Laplace approximation to the likelihood (Vonesh and Chinchilli, 1997) is implemented by setting the number of quadrature points to unity. This method is computationally less demanding than adaptive quadrature. Estimators of vi are obtained by computing the expected value of vi, given the data and the estimated model parameters, for example, the empirical Bayes estimates (Pinheiro and Bates, 1995; SAS Institute, 1999, p. 2476). The implied score on any of the two rating scales is computed by determining the category into which the estimate of vi falls. These estimated scores may be used for characterizing accessions. While less informative than vi itself, users of genetic resources often prefer estimated scores for reporting in database queries and for simple classifications of accessions according to different traits.
|
EXAMPLE
|
|---|
Mildew on barley has been evaluated by the Gatersleben genebank since 1949. Here, we analyze the first of three mildew evaluation dates for winter barley. We use data covering all years, except for 1949, 1962, 1982, and 1992, when there were no or only a few accessions per year. There were two different rating scales. The old one used until 1991 has scores from 0 to 5 and the one used since then has scores from 1 to 9. To harmonize the two scales, we shifted the old rating scale by adding one to the old scale, transforming the 0-to-5 scale to a 1-to-6 scale (Table 1, Fig. 1). There is a total of 2880 accessions, of which 1079 are represented only in years 1993 to 2002 using the new rating scale. Generally, the accession x year classification is very sparse. The mean number of winter barley accessions per year evaluated between 1950 and 2002 was 168, with a minimum of 16 and a maximum of 849. A total of 960 accessions were tested in only 1 yr; 465 accessions were tested in 2 yr; and 1457 accessions were tested in >2 yr, among these, one has been evaluated in 19 yr. Details are described in Fig. 2.
View this table:
[in this window]
[in a new window]
|
Table 1. Ad hoc conversion of scores on old rating scale (1949–1991) to scores comparable to new rating scale (1993–2002), as used by the Gatersleben genebank.
|
|
View larger version (35K):
[in this window]
[in a new window]
|
Fig. 2. Histogram for number of evaluation years per accession in winter barley between 1950 and 2002 (Gatersleben genebank).
|
|
Since both scales are anchored in the same first category (Fig. 1), we estimated the model with a common first threshold for both rating scales (Table 2). Parameter starting values for adaptive Gaussian quadrature were determined from a fit of the threshold model without a random genetic effect. The genetic variance was initialized by a mixed model fit of (1) based on observed rating scores, standardizing the residual variance to unity. SAS code used for fitting the model is described in the Appendix. Using the default settings, the NLMIXED procedure adaptively chose a single quadrature point for integrating the likelihood, which is equivalent to the Laplace approximation. Despite very good starting values, computing time on a standard Pentium PC was about 37 h. To check accuracy of the Laplace approximation, we ran the procedure with three quadrature points. The resulting parameter estimates were very similar (Table 2), as were the resulting empirical Bayes estimates (Pearson's correlation r = 0.99999992, Spearman's rank correlation rs = 0.999998). The rank correlation of estimated mean scores on the new scale by both methods was 1.000 exactly, so the rank orders were identical. For comparison, rating scores were analyzed directly based on Eq. [1] using BLUP, assuming random accession effects and fixed year effects. In accordance with current practice at the Gatersleben genebank, the scores on the old scale were transformed to scores corresponding to the new rating scale by an ad hoc rule described in Table 1. The scatter plot of empirical Bayes estimates (Laplace approximation) based on the threshold model versus BLUPs shows a slight bend, indicating some nonlinearity (Fig. 3). The nonlinearity may be partly due to the ad hoc scale change used to translate the old to the new scale. In addition, nonlinearity indicates that the threshold model is not well approximated by a linear model for the observed scores, which supports the use of a threshold model. The association between both estimates is rather high (Pearson's r = 0.983, Spearman's rank correlation rs = 0.985). Table 3 shows the frequency cross-tabulation of estimated mean scores for the old and new rating scales. The tabulation shows some discrepancy compared with the ad hoc one-to-one transformation currently used by the Gatersleben genebank (Table 1). Specifically, each of the categories 2, 3, and 4 on the old scale falls into more than one category on the new scale, while the ad hoc transformation places each category on the old scale into a single category on the new scale, suggesting that the ad hoc transformation does not truthfully reflect the transition of categories.
View this table:
[in this window]
[in a new window]
|
Table 2. Estimates (SEs) of genetic variance and thresholds for two different rating scales using the Laplace approximation and adaptive Gaussian quadrature with three quadrature points (winter barley).
|
|
View larger version (21K):
[in this window]
[in a new window]
|
Fig. 3. Scatter plot of empirical Bayes estimates based on threshold model (fitted by adaptive Gaussian quadrature) vs. best linear unbiased predictions (BLUPs), treating rating scores as metric data. Some circles repsesent more than one observation.
|
|
View this table:
[in this window]
[in a new window]
|
Table 3. Frequency table of estimated mean scores for 2880 winter barley accessions on two different rating scales (see Fig. 1) based on a threshold model estimated by adaptive Gaussian quadrature with a single quadrature point (Laplace approximation).
|
|
For comparison, we also fitted the threshold model without constraints on the thresholds, but with the constraint that both periods had the same mean 
. For each of the two periods, the resulting threshold parameter estimates are shifted by the same amount (about 0.04 for the first and –0.17 for the second period; see Table 2). The difference in the two shifts is identical to the difference in period means under the model with restriction on the first thresholds (results not shown). Also, the genetic variance estimate, as well as the log-likelihood values are identical, as are the rank orders of the empirical Bayes estimates. This shows that we have, in fact, just fitted a reparameterized version of the model with one common threshold: The strong assumption on the thresholds has been replaced by a strong assumption on the year means in both periods. The notable shift in the threshold parameters indicates that the assumption of common period means is really a strong one. This is particularly true for the second period, which has just 10 yr. With many more years in the second period, one would not expect the shift of thresholds in the second period to be quite as large. At the same time, the results demonstrate that the assumption of two common period means is equally restrictive as the assumption of one common threshold. Unless the two scales are truly connected by a common category or the number of years is large in each period, there is no justification for making either assumption. It is concluded that genebank managers and breeders should be very careful when redefining rating scales. They should make sure the new scale is connected to the old in such a way that at least one common threshold can be assumed.
|
DISCUSSION
|
|---|
This paper has shown that a threshold model can be used to analyze ordinal genebank data from different rating scales. Our example concerned the integration of two rating scales. Extension to more than two scales is straightforward. Empirical Bayes estimates based on the threshold model yielded an accession rank order which was quite similar to the ranking by BLUPs based on observed scores (rs = 0.985), though there were a considerable number of rank changes (only 16 out of 2880 genotypes had the same rank by both methods). This suggests that a mixed model analysis which treats ordinal scores as metric data will yield meaningful, though somewhat biased results, and that some gain in efficiency is possible by using a threshold model. The proposed procedure rests critically on the assumption that the two rating scales are anchored in a common threshold. To be able to use the procedure, care should be exercised when redefining rating scales. It will generally be advantageous to have partial agreement between categories on both scales.
The threshold model may also be used when there is a metric scale, such as a percentage scale (Piepho, 2002), underlying the observed ratings. In this case, the thresholds need not be estimated, but follow directly from the definition of the rating scale. In fact, the metric scale is no longer a latent scale in this case. A critical point then is the distributional assumption on the percentage scale. While, for a latent scale, we can always assume there is a monotone transformation to normality on the latent scale, this assumption is not generally tenable when the metric scale is not latent. Assuming some transformation of the percentage scale to normality can be found (Piepho, 2003b), the thresholds may be transformed accordingly. Different transformations may be compared via the log-likelihood function. This option will be investigated in the future.
We have assumed accessions as a random factor for the threshold model. Assuming accessions as fixed would not be feasible due to the large number of parameters to be estimated and the resulting breakdown of asymptotic theory (McCulloch and Searle, 2001). To obtain good estimates of the genetic variance, one needs to make sure that a sufficient number of accessions is evaluated in at least 2 yr. In the example, roughly two-thirds of the accessions were tested in at least 2 yr, which provides plenty of information on the genetic variance.
We have not employed a model separating accession-by-year interaction from experimental error. This would require replicate data for at least some accessions and/or a geostatistical model for within-trial variation (Hartung et al., 2005). Moreover, computing time would increase dramatically due to the need to integrate not only the genetic effect gi, but also the random accession-by-year interaction out of the likelihood. Apart from computational cost, it may be conjectured that the resulting genetic effect estimates will not differ largely from those based on Eq. [1], because the accession-by-year interaction variance is usually more important than within-trial variation.
In the mildew example, the Laplace approximation (a single quadrature point) worked very well and yielded parameter estimates virtually identical to those obtained with three quadrature points, which is mainly due to the large number of accessions tested each year. Adaptive quadrature, as implemented in the NLMIXED procedure of the SAS system, will automatically determine the number of quadrature points. It is our experience, however, that occasionally parameter estimates change when the number of quadrature points is increased above that determined by the default setting (this was not the case for the mildew data). Thus, it is recommended that the number of quadrature points is increased until the change in parameter estimates becomes negligible.
While our paper was motivated by data problems encountered frequently in genebanks, the methods we propose may be of interest in other areas as well. For example, the change-of-rating-scale problem is quite common in large breeding programs involving different trial sites or long-term data. In addition, different breeders may evaluate the same materials using different rating scales. If the datasets have a structure similar to the one we found at the Gatersleben genebank, our procedures may be applicable with minor modifications. Typically, however, datasets from plant breeding programs will differ from genebank datasets in several important respects. For example, plant breeding trials usually have several replicates per genotype and trial and involve some form of incomplete blocking. An efficient analysis will need to model this design structure at the trial level, requiring several random effects. Moreover, genotypes will have a specific pedigree structure, giving rise to genetic correlations, which can be conveniently modeled by additional random effects (Piepho and Pillen, 2004). While it is usually straightforward to set up an appropriate mixed effects linear predictor for a given data structure, fitting several random effects simultaneously using Gaussian quadrature may be computationally rather prohibitive. An investigation of the approach presented in this paper as applied to complex data structures commonly found in plant breeding applications would be an interesting topic for future work.
|
APPENDIX
|
|---|
The SAS code presented in Fig. A1 was used to fit the threshold Eq. [5] with a common first threshold to the barley data.
Some explanations of the code follow. Starting values for all parameters are provided in the dataset ‘p_start’. The dataset mildewdat contains dummy variables y51 to y2001 for the different years. The dummy pertaining to a year in question is set equal to one, while all other dummies are set equal to zero. For year 1950, all dummies were set equal to –1 to impose the constraint •• = 0. Year mean parameters are coded year51 to year2002, while the genetic variance is coded as S2g. The dataset mildewdat contains a variable accession identifying accession. It is important that the data be ordered by this variable, because it is used as the subject variable in the random statement. Each time the variable changes its value, subsequent records are assumed to be from a new subject with a new realization of the random effect gi. The ordinal rating is stored in the variable mildew. The variables n1 to n9 are set to zero, except that n[i] = 1 when mildew = i. Thresholds are reparameterized as a sum of increments D_ij, which are constrained to be nonnegative in the bounds statement. This ensures that threshold estimates are ordered correctly. Thresholds are obtained from D_ij using ESTIMATE statements. The variable year codes different years. This variable is needed to define the likelihood differently for the two periods using IF statements. Multinomial probabilities under the threshold model are computed using the PROBNORM function. The PREDICT statement produced empirical Bayes estimates of the random effects, which are stored in a dataset g.
|
ACKNOWLEDGMENTS
|
|---|
Thanks are due to the Institute of Plant Genetics and Crop Plant Research (IPK), Gatersleben, for supplying the data. We also thank Dr. H. Knüpffer (IPK) for drawing our attention to the problem of changing rating scales. Moreover, we appreciate constructive comments by two referees, who helped improve the paper. This research is supported by the German Research Foundation (DFG), grant number PI 377/5.
Received for publication May 13, 2004.
|
REFERENCES
|
|---|
- Franco, J., J. Crossa, S. Taba, and H. Shands. 2003. A multivariate method for classifying cultivars and studying group x environment x trait interaction. Crop Sci. 43:1249–1258.[Abstract/Free Full Text]
- Hartung, K., H.P. Piepho, and H. Knüpffer. 2005. Analysing gene bank evaluation data using geostatistical methods. Genet. Res. Crop Evol. (In press).
- Hill, R.R., Jr., and J.L. Rosenberger. 1985. Methods for combining data from germplasm evaluation trials. Crop Sci. 25:467–470.[Abstract/Free Full Text]
- McCullagh, P., and J. Nelder. 1989. Generalized linear models. 2nd ed. Chapman and Hall, London.
- McCulloch, C.E., and S.R. Searle. 2001. Generalized, linear, and mixed models. John Wiley & Sons, New York.
- Patterson, H.D. 1997. Analysis of series of variety trials. p. 139–161. In R.A. Kempton and P.N. Fox (ed.) Statistical methods for plant variety evaluation. Chapman and Hall, London.
- Piepho, H.P. 1994. Best linear unbiased prediction (BLUP) for regional yield trials: A comparison to additive main effects multiplicative interaction (AMMI) analysis. Theor. Appl. Genet. 89:647–654.
- Piepho, H.P. 1998. Empirical best linear unbiased prediction in cultivar trials using factor analytic variance–covariance structures. Theor. Appl. Genet. 97:195–201.[CrossRef][ISI]
- Piepho, H.P. 2002. Auswertung von Bonituren mit der SAS Prozedur NLMIXED. Z. Agrar. 10:30–41.
- Piepho, H.P. 2003a. Model-based mean adjustment in germplasm evaluation data. Genet. Res. Crop Evol. 50:281–290.
- Piepho, H.P. 2003b. The folded exponential transformation for proportions. Statistician 52:575–589.
- Piepho, H.P., and E. Kalka. 2003. Threshold models with fixed and random effects for ordered categorical data. Food Qual. Pref. 14:343–357.[CrossRef]
- Piepho, H.P., and K. Pillen. 2004. Mixed modelling for QTL x environment interaction analysis. Euphytica 137:147–153.[CrossRef]
- Pinheiro, J.C., and D.M. Bates. 1995. Approximations to the log-likelihood function in the nonlinear mixed effects model. J. Comput. Graph. Stat. 4:12–35.
- SAS Institute. 1999. SAS/STAT user's guide. v. 8. SAS Inst., Cary, NC.
- Schumacher, E., and H. Thöni. 1990. Auswertung von Boniturwerten. Agrarinformatik 18:51–62.
- Searle, S.R., G. Casella, and C.E. McCulloch. 1992. Variance components. John Wiley & Sons, New York.
- Thöni, H. 1985. Auswertung von Bonituren: Ein empirischer Methodenvergleich. EDV Med. Biol. 16:108–114.
- Vonesh, E.F., and V.M. Chinchilli. 1997. Linear and nonlinear models for the analysis of repeated measurements. Marcel Dekker, New York.
Related articles in Crop Science:
- THIS ISSUE IN CROP SCIENCE
Crop Science 2005 45: xiii.
[Full Text]
TOP
ABSTRACT
INTRODUCTION
