5. Analysis of adaptation and identification of subregions

In Figure 5.1 (B), X_co = 251 mm (calculated over 9 cultivars included in the study). The crossover point can also be used for assigning new locations to subregions, depending on the site mean value of the covariate.

Concluding remarks and comparison of models

Sometimes, just one or two environmental covariates can explain a remarkable portion of GL interaction variation (Saeed and Francis, 1984; Biarnès-Dumoulin et al., 1996; Giauffret et al., 2000; see also case study in Chapter 8). However, careful selection of environmental variables, on the basis of common sense and their putative importance on adaptation patterns, is recommended prior to analysis - in order to both limit the calculation process and avoid the risk of multicollinearity problems. Such problems may occur when several highly correlated variables are considered simultaneously. Multicollinearity, which may lead to dangerously unstable predictions, may be eliminated by adopting a partial least squares regression model (Talbot and Wheelwright, 1989; Vargas et al., 1998). Vargas et al. (1999) report agreement between the results of this technique and those of factorial regression for a data set including a very large number of environmental variables; this suggests, however, that multicollinearity may be a minor problem in a wide range of situations. On the other hand, a reduction in the dimensions for the factorial regression model sought via the reduced rank factorial regression approach may sometimes prove unfeasible (van Eeuwijk, 1992). Although more complex, both these alternative regression techniques offer the additional advantage of a highly informative graphical output (Aastveit and Martens, 1986; van Eeuwijk, 1992).

Compared to joint regression and AMMI analysis, factorial regression allows for an explicit assessment of the relationships of environmental variables with GL effects. Our understanding of the reasons contributing to GL interaction (both in general and for individual genotypes - as depicted by estimates of β coefficients) is thus improved. The information generated facilitates the characterization of subregions and the spatial and temporal scaling-up of results (see Section 5.8), as well as the identification of genetic resources with definite response to key environmental factors. This technique, however, requires the availability of data for major environmental variables in the complete set of test environments, a condition that may severely limit its application. Conversely, with AMMI or joint regression modelling, locations with missing environmental data can be excluded from correlation or regression analyses assessing the relationships of these variables with GL effects, but used for modelling yield responses. In various reports (van Eeuwijk and Elgersma, 1993; Vargas et al., 1999; see also case study in Chapter 8), AMMI and factorial regression identify the same environmental factors as those related to GL interaction occurrence, producing similar results in terms of genotype adaptive responses and site similarity for GL effects, despite the higher requirements in terms of factorial regression input data.

Factorial regression models including genotypic covariates allow for proper statistical assessment of the relationships of adaptive responses with morphophysiological traits that are recorded in all environments (Denis, 1988; van Eeuwijk et al., 1996). However, correlation of genotype values for these traits (the average of observations from even a small number of sites) with values of β coefficients may provide indication of the traits involved in the response to each environmental covariate (see Section 6.3). Correlation results are expected to be in logical agreement with the covariate concerned (e.g. longer crop cycle positively correlated with β values for response to rainfall amount).

A word of caution is required with respect to the identification of causal factors for GL interaction. An environmental variable - either significant as a covariate in the factorial regression model or significantly related to site ordination on significant PC axes - may be correlated with an unmeasured causal factor rather than representing a causal factor itself. This may be evident (e.g. site altitude, underlying the effect of low temperatures at some growth stages) or obscure (e.g. soil pH, possibly hiding the effect of soil aluminium as the unmeasured, correlated causal factor). Similarly, putative adaptive traits as identified by the analyses (see Section 6.3) may be genetically correlated with the actual adaptive traits. Moreover, sampling errors may underemphasize the importance of environmental or genotypic variables (especially when modelling does not contemplate the use of these variables). However, false assumptions regarding causal factors do not necessarily have a very negative impact. For example, the exploitation of site altitude as an environmental variable for scaling up the results of subregion definition for breeding or variety recommendation (see Section 5.8) may be effective if the variable is highly correlated with the underlying causal factor (e.g. low temperatures). On the other hand, specific breeding for a subregion characterized as having low soil pH, based on artificial screening for tolerance to low pH, is largely ineffective if the real causal factor for GL interaction is soil aluminium. Indications of causal factors can often be verified with small ad hoc experiments.

When different analytical models have been applied to the same yield data (with possible different indications), it must be asked which approach is the most reliable. Model comparison is sometimes based on model R² (i.e. the portion of GL interaction SS accounted for by significant GL interaction parameters in the model). This criterion (which focuses on model accuracy with respect to the same data that generated the model itself) fails to recognize that model parsimony (expressed by low number of DF) is a desirable characteristic contributing to the ability to predict genotype performance in a novel, independent data set (Gauch, 1992^[19]). For example, compared to AMMI-1, the joint regression model is mathematically incapable of accounting for more GL interaction SS, while accounting for less DF. Comparison of models in terms of MS value takes into account accuracy and parsimony, but is only appropriate for models with one parameter (joint linear regression; AMMI-1; factorial regression with one covariate). Brancourt-Hulmel et al. (1997) propose a general measure of model quality provided by the ratio:

% GL interaction SS/% GL interaction DF

for SS and DF values accounted for by the significant parameters in the model.

An alternative criterion for model comparison may relate to the estimated variance of the significant parameters. The variances of PC axes are summed up because the principal components are uncorrelated (which means that they add independent pieces of information). Similarly, the variance of any additional environmental covariate may be added to that of the first covariate as it relates to partial regression SS, i.e. independent information. If Becker's (1984) procedure for estimation of the heterogeneity of genotype regressions variance component is extended to PC axes or environmental covariates, the following formulae may be used for ANOVA models 1, 2 or 3 in Tables 4.1 through 4.3 (with the same formulations used in the tables):

where the MS of a given component of the GL interaction is indicated by M_C and the estimated variance is indicated by S_C². For example, the variances of genotype regressions, PC 1 and PC 2 components of the GL interaction calculated from MS values of the ANOVA in Table 4.4 (according to the second formula) are 0.010, 0.007 and 0.005 (t/ha)², respectively, indicating the superiority of AMMI-2 (0.012 [t/ha]²) over the joint regression model. The criterion proposed by Brancourt-Hulmel et al. (1997) would rather favour the regression model (% GL interaction SS % GL interaction DF = 0.104/0.033 = 3.15) over AMMI-2 (0.380/0.176 = 2.16), despite the very low amount of explained GL interaction (R² = 10.4%).

Specific breeding is favoured by large genotype x subregion interaction and low GE interaction within subregions. Therefore, the MS ratio of these effects (calculated in the ANOVA including the factors: genotype, subregion and environment within subregions) may be used as an additional criterion for comparing different approaches for subregion definition in a breeding perspective. Especially in the absence of clear-cut indications, the estimation of yield gains predicted for each approach (see Section 8.2) provides an ultimate, fundamental criterion for comparison.

5.5 Pattern analysis

Adaptation patterns are herein investigated using a combination of ordination and classification techniques. There is ample literature providing method description (Williams, 1976b; Byth and Mungomery, 1981; DeLacy et al., 1996a) and documenting applications (Mungomery et al. 1974; Shorter et al., 1977; DeLacy et al., 1994, 1996c). The input data set is currently represented by a two-way matrix of genotype by location yield values (rather than GL interaction effects), averaged across years and/or experiment replicates. Yield data, however, are preliminarily standardized within-location to zero mean and unit standard deviation (by subtracting location mean yield and dividing by within-location standard deviation of genotype values). Standardization is recommended to remove the effects of site mean yield and heterogeneity of genotypic variance among locations from the assessment of site similarity for GL effects (Fox and Rosielle, 1982a; see also Section 5.6). Alternative techniques for classification of locations or genotypes, based on the analysis of a GL interaction data matrix, have been developed by Ramey and Rosielle (1983) and Lin and Butler (1988). An ordination technique, contemplating the execution of a principal components analysis on the genotype by location matrix of untransformed yields, is proposed by Yan et al. (2000), mainly as an alternative to AMMI analysis-derived procedures for the graphical display of the top-ranking genotype in each location.

Classification and ordination procedures

In its classification mode, pattern analysis implies the grouping of locations through a hierarchical cluster analysis based on site similarity in the Euclidean space represented by the standardized yields of genotypes. Classification usually also concerns genotypes, grouped together on the basis of similarity for GL effects and main effects, i.e. the adaptive response (as the standardization adopted does not remove the yield differences among genotypes, while compensating for a scale effect of individual locations on these differences [Yau, 1991]). A squared Euclidean distance as the dissimilarity measure and Ward's clustering method are normally recommended (DeLacy et al., 1996a; Cooper et al., 1996a). For site classification, a possible truncation criterion may be provided by the lack of significant GL interaction within a group of locations. However, alternative criteria are also available (DeLacy et al., 1996a).

Pattern analysis by the recommended procedure has the theoretical advantage of providing a grouping of sites that reflects the opportunities for exploiting indirect selection responses to locations, as these responses are maximized within potential subregions and minimized across subregions (Cooper et al., 1996a). Another advantage of this method over previous classification techniques is the rapidity of the analysis (no prior modelling of adaptive responses is required). A third advantage is its usefulness in the analysis of largely unbalanced data sets, using a procedure described by DeLacy et al. (1996b) (summarized below).

Cluster analysis following the modelling of GL interaction effects has the theoretical advantage of assessing site similarity after separating the pattern (retained in the significant GL interaction terms, such as PC axes or environmental covariates) from the noise portion of GL interaction variation. This advantage - confirmed in investigations by Crossa et al. (1991) and Smith and Gauch (1992) in relation to AMMI analysis - may produce in comparison with pattern analysis results: a poorer fit to the data observed (also in terms of predicted yield gain on the basis of the suggested definition of subregions); but a better prediction of actual, future responses of selected material. For pattern analysis, Crossa et al. (1991) therefore suggested using the expected yields of genotype-location combinations according to the selected AMMI model, rather than the observed yield values. This may be extended to the use of expected yields as provided by appropriate modelling through joint regression or factorial regression. Without ruling out this possibility, clustering locations according to the site characteristics that are directly associated with the occurrence of GL interaction (i.e. the mean yield, the score on significant PC axes or the value of significant environmental covariates: see Sections 5.2, 5.3 and 5.4) probably provides a faster and straightforward option for cluster analysis to complement joint regression, AMMI and factorial regression techniques.

Ordination of genotypes and locations in pattern analysis is performed by principal components analysis (Kroonenberg, 1995) or principal coordinates analysis (Gower, 1966; Mungomery et al., 1974) of the standardized data. It aims mainly to complement the cluster analysis results (e.g. by highlighting the response of individual sites or genotypes, or the extent of variation within the groups of sites and genotypes). The lack of integration with ANOVA results does not provide the user with an objective criterion for testing the significance of principal components and excluding the allegedly noise portion of GL effects from modelling of adaptation patterns. Indeed, the current modelling represents the standardized values observed, usually by means of a "performance plot" showing the mean performance of each genotype group averaged across sites of each location group (for the groups that are selected by cluster analysis). In addition, data standardization may introduce distortion to the modelled responses of genotypes (see Section 5.6). These characteristics make the previous analytical techniques preferable to pattern analysis when analysis is aimed at the targeting of genotypes. Therefore, pattern analysis is here envisaged mainly for site classification aiming at the definition of adaptation strategies. However, the opportunity offered by this technique for classifying genotypes according to their adaptive response may also prove useful in various breeding contexts (e.g. selection of breeding lines or identification of potential parent material), especially when the appreciation of these responses by graphs (such as Figs 5.1 or 5.4) is hindered by the large number of evaluated entries or the multidimensional adaptation pattern. Finally, contribution analysis (Shorter et al., 1977; DeLacy et al., 1996a) is an additional tool for obtaining insight into the relative contribution of individual test sites to GL interaction and genotype classification carried out by Ward's method.

Analysis of largely unbalanced data sets

The availability of largely unbalanced data sets composed of several historical data sets, each one relating to a specific year (or crop cycle) and having several test locations but few (or no) tested genotypes in common with the other data sets, is relatively frequent in multi-environment testing of breeding lines or commercial cultivars (especially where there is a high turnover of plant material). The combined analysis of this information for location classification may be realized using a procedure that requires different steps (DeLacy et al., 1996b):

• estimation of the phenotypic correlations among test locations for genotype original yields in each individual data set;

• averaging across data sets of the correlations for each pair of locations; and

• transformation of the similarity matrix (as provided by correlation coefficients) into a dissimilarity matrix of squared Euclidean distances, inputting it (rather than the genotype by location matrix of standardized yields) into the cluster analysis.

Examples of this procedure are provided by DeLacy et al. (1994 and 1996c). Correlations based on fewer than four genotypes should not be considered. Averaging of correlations relative to the same pair of locations should be performed on values previously submitted to the z (inverse tanh) transformation, back-transforming to r the average z value (using tables provided in statistical textbooks: e.g. Steel and Torrie, 1960; Dagnelie, 1975a^[20]). A weighted average should be used for correlations based on a variable number of genotypes, using the following formula (expressed for z values):

z = ∑ (n_i - 3) z_i/∑ (n_i - 3)

where z is the weighted average, and z_i and n_i are the z value and the associated number of genotypes for the correlation, respectively, in the data set i. For example, the weighted average of the three phenotypic correlations r₁ = 0 .50, r₂ = 0.80 and r₃ = 0.90, with the respective number of genotypes n₁ = 16, n₂ = 10 and n₃ = 15, can be obtained through the z transformation as:

z = [(13 × 0.55) + (7 × 1.10) + (12 × 1.47)]/(13 + 7 + 12) = 1.01

which, once back-transformed, provides the average phenotypic correlation r_p = 0.77 for insertion in the similarity matrix. Of course, pairs of locations may differ for number of correlations contributing to the average r_p value, since some sites may be absent in some data sets. At least one individual correlation is needed for each pair of locations to allow both sites in the analysis. A higher number of individual correlations (e.g. 2) and more than four genotypes per correlation (e.g. 6-10) may be set as minimum values to increase the reliability of results. The transformation of the matrix of average correlations into the matrix of squared Euclidean distances (D values) among locations can be performed by transforming each r_p value according to Gower's (1966) relationship (DeLacy et al., 1996a):

D = 2 (1 - r_p)

The dissimilarity matrix relates to distances as they would be estimated on the basis of genotype yields standardized within location. A principal coordinates analysis is required for the possible ordination of locations in this case. By adopting appropriate, freely available software (see Section 5.9), it is possible to input the correlation matrix without any further calculation for execution of the cluster analysis. In this case, no truncation criterion based on the level of within-group GL interaction is available (as no combined ANOVA can be performed). Opportunities for subregion characterization and scaling-up of pattern analysis results are considered in Section 5.8.

5.6 Data transformation

Heterogeneity of genotypic variance

Data transformation in relation to heterogeneity of genotypic variance among environments (and in particular, locations) only concerns analyses aimed at the definition of adaptation strategies and yield stability targets. With special reference to pattern analysis of a genotype by environment matrix of yield data, Fox and Rosielle (1982a) proposed standardizing genotype values within environments, in particular to counterbalance the trend of high-yielding environments - supposedly characterized by higher variance of genotype values - having greater weight in the assessment of environment similarity for GE effects. Yau (1991) then showed that positive correlation of within-environment phenotypic standard deviation of genotype yields with environment mean yield is common in data sets characterized by large variation in mean yields. This correlation may occur because the variation in genotype yields is limited by the zero value in lowest-yielding environments. These indications can be extended to cluster analysis of locations based on site similarity for GL interaction effects, as it is substantially the case also for cluster analysis of PC scores of sites following AMMI modelling. In the presence of positive correlation between site mean yield (m_loc) and within-location phenotypic standard deviation of genotype yields averaged across replicates (s_p), higher-yielding sites are characterized by larger residuals from additivity (i.e. GL effects) and tend to show larger Euclidean distances from other locations. Conversely, low-yielding locations tend to be more similar than they really are in relation to genotype ranking. This is exemplified for the yield response of four hypothetical genotypes in four locations. Starting from responses in Data set A (s_p constant), the effect of proportionality of s_p with m_loc is superimposed, resulting in Data set B (Table 5.2). This effect results in larger Euclidean distances of the highest-yielding site 'd' from the other sites, in comparison with distances for s_p constant, with reference to the relevant dissimilarity matrix (i.e. that for GL effects) (Table 5.3). As a result, the highest-yielding site tends to cluster alone, while the remaining sites cluster together (data not reported). It can also be seen from the matrices of Euclidean distances relative to Data set A that the assessment of site similarity for GL effects based on original yields is biased by an effect of site mean yield (which inflates the distances between the sites with contrasting mean yield) even in the absence of proportionality between s_p and m_loc (Table 5.3). The within-location standardization of genotype values in Data set B (by subtracting site mean yield, and dividing by the within-site standard deviation of genotype values) can compensate for the proportionality of s_p with m_loc, providing Euclidean distances between sites (Data set E in Table 5.3) that coincide (in terms of relative differences) with those of Data set A. The distances assessed on genotype yields or GL effects coincide with standardized data (Table 5.3). Also two alternative data transformations, namely, the logarithmic one (Data set C in Table 5.2) and the division of genotype yields by site mean yield, i.e. the relative yield (Data set D in Table 5.2), provide Euclidean distances based on GL effects that coincide (in relative terms) with those of Data set A (Table 5.3).

TABLE 5.2 - Mean yield of four genotypes in four locations (a, b, c and d) and across locations, mean yield of locations (m_loc) and within-location phenotypic standard deviation of genotype mean yields (s_p)

^a s_p constant.
^b s_p proportional to m_loc.
^c Log₁₀ transformation of Data set B.
^d Relative yield of Data set B, as ratio of m_loc.

TABLE 5.3 - Squared Euclidean distance between four locations (a, b, c and d) based on genotype means in each site (G means) or genotype × location interaction (GL) effects

^a Within-location phenotypic standard deviation of genotype mean yields, sp constant.
^b sp proportional to location mean yield, m_loc.
^c Log₁₀ transformation of Data set B.
^d Relative yield of Data set B, as ratio of m_loc.
^e Standardization of Data set B to m_loc = 0 and s_p = 1.
Note: See Table 5.2 for G means of Data sets A, B, C and D.

Concern for the heterogeneity of genotypic variance among locations is justified by its irrelevance to the assessment of adaptation strategies, as it relates to a scale effect implying no change of relative merit between genotypes (Cooper et al., 1996a). Heterogeneity may arise from the relationship between s_p and m_loc, as well as from the presence of the aforementioned main crossover point. The crossover point - already highlighted for unidimensional modelling of genotype responses by joint regression (Fig. 5.1 [A]), AMMI-1 (Fig. 5.4 [A]) or factorial regression (Fig. 5.1 [B]) - is expected to also be present in more complex models, including two or more GL interaction PC axes or environmental covariates (although estimation is far more complex in these cases). The contribution of the main crossover point to heterogeneity of genotypic variance can be accepted because it is intrinsic to adaptation patterns and, in addition, can be very useful for site classification. Its elimination by data standardization (which eliminates any cause of heterogeneity of genotypic variance) has the additional disadvantage that it may distort the modelled adaptive responses of genotypes - as recognized by McLaren (1996) and exemplified in Figure 5.8 with reference to responses to site mean yield of two hypothetical genotypes. In this situation, where both the proportionality of s_p with m_loc and the main crossover point contribute to heterogeneity of genotypic variance, a logarithmic transformation could remove the former while maintaining the latter. However, data standardization can be recommended for site classification by pattern analysis (in which no modelling of genotype responses is preliminarily performed).

Figure 5.8 - Yield response of two hypothetical genotypes modelled as a function of site mean yield for original, log-transformed and within-location standardized data under the assumption of within-location phenotypic standard deviation of genotype yields proportional to site mean yield

Assessing the extent of heterogeneity of genotypic variance among environments prior to other analyses (as proposed in Section 2.6) also provides information on the extent of heterogeneity among locations. Large heterogeneity among environments may also be troublesome for assessing yield stability targets for breeding, because it inflates the estimation of relevant GE interaction components of variance by effects that do not relate to change in the relative response of genotypes between environments. If the results point to the transformation of data (see Section 2.6), one solution could be to remove mainly the relationship of s_p with m_loc, while maintaining the main crossover point. A decision may be based on the regression of the within-location phenotypic variance of genotype mean yields (s_p²) on m_loc, with both terms expressed on a logarithmic scale. A regression slope:

• b ≈ 2 (implying the relationship s_p ≈ k m_loc) suggests a logarithmic transformation;
• b ≈ 1 (implying the relationship: s_p² ≈ k m_loc) suggests a square root transformation; and
• b ≈ 0 (implying no relationship of s_p² with m_loc) discourages any transformation.

This approach is similar to that suggested for coping with heterogeneity of experimental errors (see Section 4.4). It also explains why a logarithmic transformation could effectively remove the effect of proportionality of s_p with m_loc (i.e. b ≈ 2) in Tables 5.2 and 5.3 (Data set C) and Figure 5.8. However, the relationship b ≈ 1 may also occur, justifying a different transformation.

Results for four data sets (Table 4.7) suggest that a relationship of s_p with m_loc may occur where there is wide variation in site mean yield including very low yields. In particular, significant correlation is only found in the Algerian data set, showing the largest difference in relative terms between extreme values of site mean yield. The regression value (b = 1.92) definitely suggests a logarithmic transformation, which can also help stabilize error variances (although low priority indications favour a square root transformation). Modelling of genotype adaptive responses as a function of the site score on the first GL interaction PC axis is shown in Figure 5.9 for original and log-transformed yield values of this data set. Data transformation removed the clear trend of higher-yielding sites towards larger variation in genotype mean yields, while safeguarding the occurrence of the main crossover point (CP). It also determined a change in the site classification based on the CP criterion (which includes three - rather than only one - sites in the smaller group) or several other criteria, and could considerably reduce the relative size of the heterogeneity of genotypic variance among environments (see Section 8.2).

Figure 5.9 - Nominal yields of genotypes modelled as a function of the first GL interaction PC axis of locations for original and log-transformed data, with indication of scaled PC 1 scores of sites (black triangles) and estimation of the main crossover point (CP)

Note: Data relate to the case study in Chapter 8.

The importance currently attributed to the relationship of s_p with m_loc for decisions on data transformation in the presence of large heterogeneity of genotypic variance among environments justifies the previous indication (Section 2.6) to verify, in terms of correlation, this relationship as a faster alternative to estimating variance components relative to heterogeneity of genotypic variance and lack of genetic correlation among environments. Indeed, the lack of correlation between s_p and m_loc implies either a relatively low level of heterogeneity of genotypic variance relative to the lack of genetic correlation component, or (rarely and on the basis of experience and common sense) the opposite situation (in which hardly any data transformation can be suggested - besides the data standardization for pattern analysis).

Heterogeneity of error for genotype comparison on specific sites

Genotype comparison at specific site values in the analysis of adaptation using the mean value of Dunnett's one-tailed critical difference (see Sections 5.2 to 5.4) is not reliable when site ordination for mean yield (joint regression), PC scores (AMMI) or environmental covariates (factorial regression) is correlated with the appropriate error for the comparison on each site (i.e. the within-location GY interaction for trials repeated in time, and the pooled experimental error for trials not repeated in time). This may occur, in particular, when site ordination is coincident or strictly correlated with site mean yield, especially in the presence of large variation in site mean yield including very low values. For example, it would take place in the joint regression analysis of the Algerian data set in Table 4.7, for which the within-site GY interaction mean square (M_gy) varies as a function of site mean yield (m_loc); on the contrary it would not occur for any of the other data sets, where no such relationship can be found. In situations of concern, transforming data may increase the reliability of genotype comparison. For trials not repeated in time, the data transformation discussed in Section 4.4, which removes the effect of site mean yield on the heterogeneity of experimental errors, is beneficial also in this context. For trials repeated in time, the regression of the M_gy values as a function of m_loc on a logarithmic scale may indicate a transformation for removing the effect of site mean yield on the heterogeneity of within-site GY interaction, using the same criterion already discussed in relation to heterogeneity of genotypic variance among sites. This suggests a logarithmic transformation (b = 1.67) for the Algerian data set modelled by:

• joint regression;

• AMMI analysis, if mean yield and PC 1 of sites are strictly correlated; or

• factorial regression, if one or more covariates are closely associated with site mean yield.

Data transformation may be envisaged only when the statistical comparison of genotypes on some sites is of specific interest, because of the difficulties arising in various aspects of the analysis (see Section 4.4) and the fact that it rarely has a significant effect on genotype merit and the crossover points of each site. Its potential adoption for the description of adaptive responses in terms of yield reliability may be adequate only if site mean yield also affects the average yield stability in time of genotypes (see Section 7.2).

5.7 Unbalanced data sets

Situations and possible options

The analysis of largely unbalanced data comprising several historical data sets aimed at the definition of adaptation strategies can be performed by pattern analysis (see Section 5.5). The same analysis may also give indications for genotype recommendation, using site groups (defined on the basis of site similarity for GL interaction effects relative to the whole germplasm sample) as possibly different recommendation domains, and identifying best-yielding genotypes across environments of each group of sites with procedures which take account of the high degree of data imbalance (i.e. least squares or, preferably, BLUP-based genotype means). Site groups characterized by specific best-yielding material define the different subregions. For largely unbalanced data sets, an alternative procedure for the definition of genotype recommendations may be obtained from multiple regression models of original yields as a function of the available site covariates developed independently for each genotype (see Section 5.4).

Methods for modelling genotype adaptive responses in the presence of some empty genotype by location cell means are available for AMMI analysis (Gauch and Zobel, 1990; Gauch, 1992^[21]) and conceivable for ordinary factorial regression (Piepho et al., 1998). However, because of their complexity, it is necessary to adopt specific, powerful software (see Section 5.9). When at least one plot value can be available for each genotype by location cell mean, the imbalance in the number of observations per cell mean (determined by missing plot values, variable number of replicates across trials or variable number of test years) may represent a minor problem for the analyses. As mentioned earlier (Section 4.1), estimating missing plot values or analysing genotype by environment cell means is recommended, whenever possible, to eliminate data imbalance. In the remaining cases (e.g. variable number of test years per location in trials repeated in time), the least squares estimation of genotype by location cell means (i.e. genotype means across environments of each location) is recommended. If this is not possible, the analysis of adaptation should be based on a balanced data set (i.e. genotype by location cell means based on same number of test years), once less represented locations have been eliminated. This is sometimes preferable also for defining a more reliable model for analysis of adaptation, through the elimination of cell means estimated with lower accuracy. This occurred in the case study in Chapter 8, where ANOVA and analysis of adaptation were carried out on a balanced data set including 14 locations with two test years each. Three locations with only one test year were then attributed to subregions by different methods depending on the classification criterion.

Balanced data sets simplify the analysis and offer better opportunities for scaling up results. Even in the presence of intense turnover of commercial cultivars, it is preferable to test a fixed set of cultivars for recommendation during a given time lag (e.g. two years for annuals), and then update the set of material for the next evaluation cycle.

Assignment to subregions of test sites excluded from analysis

Test locations including a subset of genotypes or a smaller number of test years, may be assigned to subregions on the basis of genetic correlation (r_g) analysis. Falconer (1989^[22]) pointed out that the same character (in this case, yield) evaluated on a set of genotypes in two environments (or locations) could be considered as two distinct characters, controlled: almost entirely (high r_g value); partly (low positive r_g); or not at all (r_g ≈ 0), by the same genes or sets of genes. Negative r_g values may also occur - an indication that genes favouring adaptation to one environment (or location) are unfavourable for adaptation to another. The genetic correlation of genotype values between two locations j and j' may be estimated (Burdon, 1977) as:

r_g = rp/(hj hj')

where rp is the phenotypic correlation of genotype values on site j with those on site j', and hj and hj' are square roots of broad sense heritability (h²) values estimated on a genotype mean basis for sites j and j', respectively. Each h² value can be estimated from components of variance through ANOVAs performed for each individual location:

h² = s_g²/(s_g² + s_e²/r)

h² = s_g²/(s_g² + s_gy²/y + s_e²/y r)

where s_g², s_gy² and s_e² are, respectively, genotype, GY interaction and experimental error components of variance, r is the number of experiment replicates and y the number of test years (or crop cycles) for the location. The first formula is appropriate for one test year on the site, and the second formula for two or more test years. Variance components can be estimated (as in Section 4.3) for sites including one environment. For sites with two or more environments (corresponding to different years), they are estimated usin_g formulae for Model 1 in Table 4.1.

These formulae may be used to identify the location included in the analysis of adaptation which is most similar in terms of genotype responses (i.e. with the highest r_g value) to a given location not included in the analysis. Conclusions concerning the former location (e.g. in terms of recommended material) may thus be extended to the latter. It is preferable that r_g values are calculated with reference to the same test year(s) for the sites.

Alternatively, the same approach may be used to identify the mega-environment m (among those issued from analysis of adaptation) most similar (in terms of genotype responses) to the site j not included in the analysis and to which the site j can therefore be assigned. In this case:

r_g = r_p/(h_j h_m)

h_m² = s_g²/(s_g² + s_ge²/e +s_e²/er)

where r_p is the phenotypic correlation of genotype mean values across the subregion m with genotype values on site j; h_m² is the broad sense heritability value for the same subregion; s_g², s_ge² and s_e² are, respectively, genotype, GE interaction and pooled error components of variance, estimated by formulae for Model 1 (Table 4.1) from the ANOVA including the subset of environments (as individual trials) pertaining to the subregion; e is the number of environments and r the number of experiment replicates in the ANOVA; and h_j is estimated on the basis of earlier formulae applied tosite j.

The above procedures based on genetic correlations are particularly useful when subregion definition derives from application of AMMI or pattern analysis. When subregions are defined by joint regression or one-covariate factorial regression modelling, an alternative and simpler procedure is to assign test sites (with only a subset of evaluated material or test years) to the subregion where the range of site mean yields for a common set of genotypes and test years (joint regression), or the range of site values for the relevant covariate across the same years (factorial regression), includes the value of the site. For multicovariate factorial regression, the assignment to subregions may also derive from cluster analysis of all test sites for values of relevant covariates recorded in the same years.

5.8 Characterization of subregions and scaling-up of results

Scopes and data requirements

The definition of subregions, initially limited to individual test locations, should ideally be extended to cover the entire target region. The scaling-up of results is more objective if it is able to exploit information on environmental variables associated with the occurrence of GL interaction. Hereafter, values of individual sites for these variables are divided into:

• test year data (averaged across test years, for climatic or other variables; specific to the experiment site, for soil or crop management variables); and

• mean values (long-term, for climatic or biotic variables; average values for the area including site, for other variables).

The use of long-term rather than test year data for climatic or biotic variables may also allow for a temporal scaling-up, in which test locations are re-assigned to subregions. Likewise, the use of mean values rather than test year data for soil or crop management factors may imply a re-assignment of test sites to subregions.

Different procedures can be used depending on:

• the aim of subregion definition (cultivar recommendation or breeding);

• the analytical method adopted for subregion definition; and

• the availability of data for environmental variables associated with the occurrence of GL interaction (no data; only test year data, for all or most test locations; also site mean values, for test locations as well as other sites; GIS data).

Techniques allowing for accurate scaling-up of results may also require extensive calculation. Their adoption for the final definition of subregions for breeding may be viewed as a refinement of the initial definition of subregions based on test locations when, for these subregions, there is evidence of some advantage for a specific adaptation strategy. In this context, the scaling-up procedure can improve the estimation of the relative size of each candidate subregion, relative to the estimation based on the relative frequency of test sites assigned to each subregion; this information is then used in formulae for comparing adaptation strategies (see Sections 6.1 and 8.2).

Joint regression analysis-based results

For subregion definition based on joint regression modelling (for whatever purpose), locations can be assigned (or re-assigned) to subregions depending on their long-term mean yield. However, critical levels of site mean yield for definition of subregions relate to the average performance of several tested genotypes, of which only some can be expected to be widely grown in the region. Therefore, these critical levels may be redefined in relation to mean yield of one or a few widely grown control genotypes (raising or lowering the original levels, depending on whether the response of the control genotype[s] is higher or lower than the average response of all tested material). Ideally, the control genotype(s) should not possess a markedly specific adaptation, to avoid a bias in the estimation of site mean yield at very low or high yield levels. Scaling up the results of joint regression modelling may prove more difficult to realize than anticipated, due to the difficulty in obtaining a reliable estimation of the long-term yield of locations.

AMMI analysis- and pattern analysis-based results

When no information on environmental variables is available or no correlation is found between test year data of these variables and scores on significant GL interaction PC axes of sites, subregion characterization following AMMI modelling can rely only on the supposedly close relationship between geographic proximity of locations and site similarity relating to either best-ranking material (genotype recommendation) or overall GL effects (definition of adaptation strategy). On a geographic map of the target region, new sites may be attributed to the subregion to which the nearest test site belongs. The reliability of the procedure increases as the density of test sites in the region increases.

The knowledge of correlations between environmental variables and ordination on significant PC axes of sites can be exploited to characterize subregions and make their geographic definition less arbitrary. Correlations may also be assessed when environmental data are only available for a subset of test sites; it is here that AMMI modelling is of special interest, because other models requiring the observation of environmental variables in the complete set of test environments cannot be used. In the simplest case, subregions can be characterized by mean values of relevant environmental variables across the test sites they include, and each new site can be assigned to the most similar subregion on the basis of the known, general level of the variables in the area of the site. This simple procedure (used for the provisional definition of subregions for breeding - see Fig. 5.7) has proved quite useful in the exploitation of GL interaction effects by growing specifically adapted material (Annicchiarico, 1998).

The availability (for relevant environmental variables) of mean values for test sites and many other locations may permit a real up-scaling of AMMI analysis results, which may be achieved in two steps:

• Firstly, a relationship is established between adaptation patterns and test year data of environmental variables recorded on test sites.

• Secondly, the relationship is exploited to predict future responses or define subregions on the basis of mean values (i.e. the most probable values) of the environmental variables on new sites or test sites.

The recommended procedure differs according to the purpose of the analysis and the criterion adopted for site classification.

For cultivar recommendation, an equation for estimating the scaled PC scores of sites as a function of environmental variables recorded in test years is searched for by means of a stepwise multiple regression analysis for each significant PC axis. The main drawback of this analysis is represented by the risk of multicollinearity (Aastveit and Martens, 1986). To limit the risk, the simultaneous assessment of several highly correlated environmental variables should be avoided, while forward selection of variables is preferable to backward elimination for model selection (Dagnelie, 1975b^[23]; Aastveit and Martens, 1986). The stepwise method, if available, is generally preferred. The procedure is considered sufficiently reliable if, for each PC axis, a large part of variation (e.g. R² = 60%) can be explained by the selected multiple regression equation. If so, from an equation of the following type for a given axis (e.g. PC 1):

v_j' = a + ∑ b_n V_jn

where vj' = scaled score of the site j; a = intercept value; b_n = partial regression coefficient of the significant (e.g. P < 0.10) environmental variable n; and V_jn = value on the site j of the variable n, the scaled scores for new sites or test sites can be estimated as a function of site mean values of the environmental variables. These scores can then be inputted in formulae, such as [5.4] or [5.5], in order to estimate nominal yields of the genotypes in each location. More simply, graphical expressions, such as Figure 5.4 (A) for AMMI-1 and Figure 5.5 (A) for AMMI-2, can be exploited by reading, for each location, the best-ranking genotype(s) expected on the basis of the estimated site score(s). Finally, locations with the same best-yielding material are grouped together to form one subregion. The opportunity to interface GIS data (rather than data for a number of individual sites) allows for a very fine-tuned definition of subregions on a geographical map following the scaling-up of AMMI (or factorial regression) results (e.g. for subregions having a distinct pair of top-ranking cultivars - Annicchiarico et al., 2002c).

The above procedure may also be used to calculate and report the expected top-yielding genotype(s) for given levels of relevant environmental variables. For example, Figure 5.10 shows the expected pair of winning cultivars for each combination of 10 rainfall by 10 winter temperature levels of sites for the case study data in Chapter 8. These variables are included in the selected multiple regression model for estimating the scaled score of sites on PC 1 (the only significant PC axis). For any combination of levels of environmental variables, nominal yields of genotypes are estimated with formula [5.4], using the site score on PC 1 (v_j1') predicted by the regression equation for the hypothetical site j possessing the given climatic characteristics:

N_ij = m_i + u_i1' [-2.025 + (0.00176 × rainfall) + (0.1612 × winter mean temp.)]

Figure 5.10 - Expected pair of top-ranking cultivars for yield at each combination of 10 annual rainfall by 10 mean winter temperature values of locations, according to an AMMI-1 model in which site score on PC 1 is predicted by multiple regression as a function of the two environmental variables (source of data: case study in Chapter 8)

The graphical expression in Figure 5.10 is similar to that in Figure 8.3, although the current procedure exploits only indirectly the information on environmental variables. Both types of graph can allow for spatial and temporal scaling-up of results by defining the top-ranking, recommended cultivars as a function of the site mean value of environmental variables at a given site (for an accurate definition, a denser grid of points for these variables is recommended). Since recommendation domains of locations can simply be read by users on the basis of the environmental variables, subregions do not need to be represented on a geographical map here. The graphs, however, can only be produced for AMMI models (with one or more PC axes) in which no more than two environmental variables are involved in the multiple regression equations used for prediction of site PC scores, or for factorial regression models with no more than two covariates. Moreover, approximating site values of environmental variables in the graphs may introduce some degree of inaccuracy in the recommendations. These limitations could be overcome by the development of a simple "Decision-aid System" that outputs the expected nominal yields of the cultivars (possibly together with a Dunnett's one-tailed critical difference value) as a function of the inputted site mean value of environmental variables.

The scaling-up of AMMI analysis results relative to subregions for breeding, defined according to the main crossover criterion, can be based largely on the same procedure. If a reliable estimation of site scores on PC 1 as a function of environmental variables proves possible on the ground of a stepwise multiple regression analysis performed on test year data, the equation may be exploited to calculate the expected PC 1 score of new sites or test sites for which mean values of the relevant variables are available. Sites are assigned to either subregion on the basis of the estimated PC 1 score relative to the main crossover point acting as the cut-off value.

When site classification for breeding purposes is derived from cluster analysis of site scores on significant GL interaction PC axes, an alternative procedure may be used based on discriminant analysis of groups of test locations as a function of test year data of relevant environmental variables. The same procedure may be used for subregion definition based on pattern analysis. The case described below refers to two subregions (A and B). It is probably the most common in breeding and the only case in which multiple regression analysis can also be used. Multiple regression analysis has the advantage of being more widespread and easily available in computer packages than formal discriminant analysis. Site groups should be coded as a dummy variable assuming two possible values: Y_A for each site of subregion A and Y_B for each site of subregion B (Dagnelie, 1975b^[24]; Afifi and Clark, 1984^[25]), equal to:

Y_A = N_B/(N_A + N_B)

Y_B = - N_A/(N_A + N_B)

where N_A and N_B are the number of test sites classified in subregions A and B, respectively.

The dummy variable is used as the dependent variable in a stepwise multiple regression analysis that holds test year data of environmental variables as the independent variables. The risk of multicollinearity and the above recommendations concerning the selection of variables for assessment and inclusion in the model, also apply here. The best model including significant (e.g. P < 0.10) variables should explain a sizeable portion of variation for the dummy variable (e.g. R² > 40-50%) to be considered reliable enough for discrimination. If so, the expected values (Y_j) of the dummy variable for the site j, calculated on the basis of test year data of significant environmental variables:

Y_j = a + ∑ b_n V_jn

where a = intercept value, bn = partial regression coefficient for the environmental variable n, and V_jn = value on the site j of the variable n, can be used for estimating the dividing point C for discrimination of the two groups of sites (subregions):

C = (Z_A + Z_B)/2

where Z_A = mean of Y_j values for test sites of subregion A, and Z_B = mean of the Y_j values for test sites of subregion B. The same multiple regression equation is then used for predicting Y_j values of new sites or test sites on the basis on site mean values of relevant environmental variables. Each site is assigned to subregion A if Y_j > C, and to subregion B if Y_j < C. From the multiple regression R² value, it is possible to estimate the correct probability of site misclassification (Afifi and Clark, 1984^[26]). Formal discriminant analysis is required when the procedure involves three or more subregions for site classification.

Factorial regression analysis-based results

Factorial regression modelling and complementary analyses make scaling-up of results more straightforward. For cultivar recommendation, nominal yields of genotypes can generally be estimated for new sites or test sites as a function of site mean values of the significant covariates in the model through formula [5.7]. For models with just one or two covariates, graphical expressions (such as Fig. 5.1 [B] and Fig. 8.3) can considerably simplify the task: the best-yielding material is read as a function of the covariate value(s) of the sites, and locations with the same winning germplasm form a recommendation domain. Developing a simple Decision-aid System (in which site mean values of covariates are inputted) may be of special interest for defining recommendations based on genotype yields expected from multiple regression models developed independently for each genotype (as may be the case when analysing largely unbalanced data sets).

For the definition of subregions for breeding based on the main crossover criterion, new sites or test sites can be allocated to either subregion depending on the mean value of the covariate. Finally, subregion definition by cluster analysis implies the scaling-up of results (spatially and temporally), if based on site mean values of significant covariates and including new sites besides the test sites.

Limitations and verifications

Similar procedures may be used for scaling up the results obtained through other analytical methods. Techniques incorporating environmental covariates for modelling (e.g. the partial least squares regression) facilitate the extension of results over the region compared with those that do not (e.g. the shifted multiplicative model); however, their adoption requires the availability of environmental data for all trials.

It is likely that the mean value of relevant variables (mean yield, for joint regression; environmental variables, for other methods) for some new sites exceeds the range of test year values used for modelling by joint regression or factorial regression, or for multiple regression analyses applied to AMMI analysis results. The scaling-up of results is less reliable for these sites. The difficulty in identifying actual causal environmental factors (see Section 5.4), combined with the error in the estimation of model parameters, is likely to determine a rather high error for results of any up-scaling procedure (where the error size can be properly assessed only on the basis of new data). For test sites, the error level associated with the procedure may even offset the theoretical benefit of the temporal scaling-up, leading to predictions (e.g. for best-yielding genotypes) which are no more accurate than those obtained through modelling of experimental data (Annicchiarico et al., unpublished results). However, these disadvantages are largely compensated for by extending results to a number of non-test locations in the region.

Information from new trials on old or new test sites including a subset of genotypes could be used for verifying the results predicted for the locations (Gauch, 1992^[27]; Annicchiarico, 1998). In particular, for areas scarcely represented by test sites in the analysis of adaptation, the subregion definition could be verified and, if necessary, modified by means of small experiments including:

• only genotypes recommended in at least part of the target region (for variety recommendation); and

• a few genotypes characterized by contrasting adaptation pattern and a specific ranking order in each subregion (for specific breeding).

The latter case, in particular, may be considered an application of the concept of probe (Fox and Rosielle, 1982b) or indicator (Bramel-Cox et al., 1991) genotypes, by which a synthetical characterization of environments is searched for indirectly by the response of plant material. For example, the three cultivars, 'La Rocca', 'Prosementi' and 'Orchesienne', in Figure 5.4 (A) could be used for verifying the subregion definition in both respects (variety recommendation and specific breeding); on the other hand, the first and the third varieties may suffice for verifying the membership of sites to either of two subregions for breeding defined according to the main crossover point (equal to X_co = -0.19, it almost coincides with the crossover point of the two varieties).

Scaling-up of yield reliability responses for cultivar recommendation

Modelling of genotype responses to locations is herein envisaged in relation to the most probable value expected for site parameters, exploiting directly (for factorial regression) or indirectly (for AMMI analysis) the available information on relevant environmental variables. Within-site temporal variation may also be taken into account for recommendation of more stable-yielding material as depicted from expected responses in unfavourable years. However, assessing expected yields of genotypes on each site across a range of possible values for climatic or biotic variables can be very complex, especially when more variables are involved. Moreover, this approach requires the modelling of GE rather than GL effects (with implications for model complexity - see Table 2.1), in order to not neglect the effect of variables associated markedly with GY interaction within sites (i.e. yield stability) and limitedly with GL interaction (as may be the case for climatic or biotic variables characterized by wide year-to-year and limited site-to-site variation).

For the frequent case of trials repeated in time, a simpler alternative is represented by the synthetical assessment of yield stability of genotypes from the observed temporal variation of yield responses within locations. Should this information prove important (i.e. if genotypes differ in yield stability), it may then be integrated with that on adaptive responses to site mean values of environmental variables. In particular, substituting a yield reliability value (for a given level of risk aversion) for the mean value of genotypes in formulae used for the calculation of nominal yields (e.g. equations [5.4], [5.5] or [5.7]) allows for modelling the yield reliability rather than the mean yield of genotypes as a function of site characteristics (see Section 7.2). The described procedures for spatial and temporal scaling-up of results for cultivar recommendation can then be applied to crossover interactions between top-ranking material for yield reliability estimated as a function of site mean values of relevant environmental variables (see Section 8.3).

5.9 Computer software

IRRISTAT software

Most of the techniques considered for analysis of adaptation can be applied through ordinary, relatively inexpensive statistical software in combination with substantial worksheet and manual calculation, following the techniques described in previous sections. However, the adoption of IRRISTAT can enlarge the set of applicable techniques and considerably reduce the amount of calculation (worksheets are also possible with IRRISTAT). In addition, various graphs of possible interest can be outputted by the software. Indications on the use of IRRISTAT provided here or in other sections are exemplified in the analysis of the case study (Chapter 8).

IRRISTAT modules of specific interest for analysis of adaptation have been developed for investigating genotype responses to generic environments. Their adoption may require some modification to the recommended procedure when trials are repeated in time, in which environments correspond to location-year combinations rather than individual locations. Joint regression and AMMI modelling can be performed by the G × E Interaction module. The input data file may conveniently be the one of genotype by location cell means that can be outputted by the Single Site Analysis module for trials not repeated in time, and by the ANOVA module (used for the combined ANOVA) for trials repeated in time. In the former case, a weighted analysis (in which the cell means are weighted proportionally to the error MS of sites: see Section 4.4) can also be available for analysis of adaptation. Some missing genotype by location cell means can be accepted by the program, which estimates their value by a proper algorithm. Joint regression analysis can be requested by indicating the yield response variable as the Site Index in the Stability Analysis and AMMI Model window. The number of PC axes for the AMMI model can be specified in the same window. Usually, several PC axes are indicated in a former AMMI analysis that is aimed at model selection, whereas the number of PC axes for the selected AMMI model is specified in a second analysis. The output comprises mean yields and main effects of genotypes and locations, as well as the matrix of GL effects. For joint regression analysis, b parameters of genotypes are estimated and tested for difference to unity (using the deviation from regression of the individual genotype as the error term). The intercept values of genotypes are not printed out but they can be estimated for each genotype i as:

a_i = (m_i - m b_i)

using the outputted information relative to slope (b_i) and mean value (m_i) of the genotype and the grand mean (m). The output also provides the DF and SS values for heterogeneity of genotype regressions (termed "treatment x site regression") and deviations from regressions, as well as F test results for the former component using the latter as the error term. For integration with ANOVA results (such as in Table 4.4), SS values must be multiplied by N' = (no. test years [or crop cycles] × no. experiment replicates), before calculating correct MS values. Then, the deviations from regressions MS can be tested on the appropriate error MS. The estimated parameters of genotype regressions may be used for producing a graph such as Figure 5.1 (A), or for estimating the main crossover point (according to equation [5.3]). The expected yields of genotypes on each site are stored in an output data file that may also be used for producing graphs of genotype regressions. For AMMI analysis, the output includes the scaled PC scores for genotypes and locations and, for the requested AMMI model, the DF and SS values of PC axes and residual GL interaction and the expected yields of genotypes on test sites. Mean yields and PC scores of genotypes and locations are also stored in another output data file. The outputted SS values for PC axes and residual GL interaction need to be multiplied by N' before calculating MS values and testing PC axes on the appropriate error MS in the combined ANOVA. The expected yield values for the selected AMMI model are sufficient for the identification of the best-performing entries on each test site. By subtracting the site main effect from these values the nominal yields of genotypes on test sites may be calculated. Alternatively, nominal yields can be estimated through worksheet calculation from outputted information on mean yields of genotypes and scaled PC scores of genotypes and locations (according to equations [5.4], [5.5] etc.). AMMI-1 nominal yields for the two sites with extreme PC 1 scores can be used for producing a graph such as Figure 5.4 (A), or for estimating the main crossover point (according to formula [5.6]). For AMMI-2 models, the output includes a graph that shows the genotype with the highest expected yield in the range of PC 1 and PC 2 site scores (like Figure 5.5 (A), with PC 1 and PC 2 on abscissa and ordinate axes, respectively). Other graphs can also be obtained showing mean values and scaled scores on PC 1, or scaled scores on PC 1 and PC 2, for genotypes, locations or both. The clarity of biplots is enhanced if the first two characters of the genotype and location classification variables, which are used for representation of items, are unique. Finally, the singular value (l_n) for each PC axis n is also outputted by the program. Its square root may be multiplied by the scaled PC score of sites to obtain the original PC score values to be inputted in the cluster analysis of locations.

Although no specific IRRISTAT module has been developed for factorial regression modelling, this method can be applied through some procedures that imply the analysis of GL interaction effects. A new variable including these values can be defined in the file of genotype by location cell means according to formula [4.1], after creating two more variables that assign genotype mean yield and site mean yield values to each genotype-location combination through IRRISTAT's logical operators. The calculated GL effects may be checked by comparison with those outputted by the G × E Interaction module. An additional variable should be created for each environmental covariate, assigning site values to each genotype-location combination through logical operators. Each of the possible one-covariate models can be assessed by the G × E Interaction module, specifying the GL effects as the response variable and the considered environmental variable as the Site Index. The best of these models can easily be identified as the one with the highest R² value (ratio of heterogeneity of genotype regressions SS to total GL interaction SS) in the output. Again, the outputted SS values for genotype regressions and deviations from regressions must be multiplied by N' in order to be able to calculate MS values and test them on the appropriate error MS in the combined ANOVA. In the output, regression slopes correspond to β values and are conveniently tested for difference to zero with respect to deviation from the model of the individual genotype. If a one-covariate model was ultimately selected, the intercept value α_i for each genotype _i could be estimated from the relevant β_i value and the mean value of the covariate (V) by the formula: α_i = -V β_i. Nominal yields of genotypes on test sites, on the other hand, could be estimated according to formula [5.7], either using worksheets, or by adding the genotype mean value to the expected GL interaction effect stored in an output data file for each genotype-location combination. Nominal yields for the sites with extreme values of the covariate could then be used for producing a graph such as 5.1 (B), or for estimating the main crossover point (according to formula [5.8]). The assessment of two-covariate models (all including the best single covariate) is more time-consuming. For each model, it is necessary to:

• perform a multiple regression analysis of GL effects as a function of the two covariates for each genotype, using the Regression module with the Multiple data selection option; and

• sum up SS values for regression and residual terms over the analyses, calculating the R² value (ratio of genotype regressions SS to total SS).

In the F test for the additional covariate (conveniently limited to the model with highest R²), it is necessary to multiply SS values for regressions and residual terms by N' before calculating MS values and testing them on the appropriate error MS in the combined ANOVA. The possible assessment of models with three or more covariates can be carried out using the same procedure, with the Regression module. Estimates of genotype parameters α_i and β_in in equation [5.7] can be found, together with testing results of β_in values for difference to zero, in the output of multiple regression analyses for the individual genotypes. Nominal yields of genotypes on test sites can be estimated: by using the same equation in a worksheet; or, more easily, by requiring the Regression module to store the expected GL interaction values for each genotype in a previously created variable in the input data file, and then adding these values to the genotype mean value for each genotype-location combination. Likewise, a graphical expression of winning genotypes for a two-covariate model (as exemplified in Fig. 8.3) can be obtained after estimating nominal yields of genotypes for different pairs of covariate values according to equation [5.7] using worksheets.

Pattern analysis (of special interest here in its classification mode) is performed using IRRISTAT's Pattern Analysis module. The input data file is that for genotype by location cell means already considered for AMMI and joint regression analysis. The type of data and classification method are specified in the Pattern Model window. The default option, contemplating the within-location standardization of genotype values and Ward's incremental SS method, respectively, is generally recommended. Other clustering methods, including average linkage (Group average), are also available. In the same window, the classification can be requested for: locations (considered as columns of the genotype by location data matrix); genotypes (considered as rows); or both. Other options of major interest concern:

• the output of dissimilarity matrices for genotypes and locations;

• the execution of a contribution analysis and the printing of mean values - for the specified number of genotype and location groups and the type of data (i.e. the standardized); and

• the scale for the dendrogram axis (with Fusion level representing more faithfully the degree of similarity, and Fusion number often providing a clearer output of classification results).

Dendrograms that summarize the classification results for sites and genotypes can be printed out, together with additional information on fusion stages. ANOVAs including a variable number of location and genotype groups are performed, together with contribution analysis, for the specified type of data. In its ordination mode, this module can produce results, plots and biplots relative to a principal components analysis performed on the standardized genotype by location mean values. This module also allows for classification and ordination based on GL interaction effects (requested as Mean polish type of data). Ordination results coincide with those of AMMI analysis in this case, while site classification results may differ from those of cluster analysis performed on significant GL interaction PC scores, because site similarity is currently assessed with regard to all GL effects (including the noise portion of GL interaction variation).

Cluster analysis based on original or scaled PC scores of sites (complementing AMMI analysis) or environmental variables (complementing factorial regression analysis) can be performed through IRRISTAT's Pattern Analysis module, albeit with some difficulty. In particular, site classification should be requested using each significant PC or environmental covariate variable as if it were a genotype, i.e. using a matrix of v (variables) by l (locations) as if a matrix of g (genotypes) by l (locations). Therefore, the input data file should contain: one variable that identifies locations; a second variable reporting the name of the PC axis or covariate in the analysis; and a third that includes the PC score or covariate value for the specified observation (i.e. the location and PC or covariate name). Raw (i.e. non-standardized) type of data and the appropriate classification method (e.g. Ward's) should be indicated in the Pattern Model window. The dendrogram summarizing the site classification is the only output of interest. Unfortunately, no output is produced when only one PC axis or environmental variable is included in the analysis (because a matrix of at least two genotypes/variables is expected by the program in view of its conventional use). However, a simple artifice may allow the analysis to be done in this situation. It is sufficient to include a second PC axis or covariate as a dummy variable of which the randomly assigned values vary within a negligible range of variation (e.g. between -0.00005 and 0.00005) compared with PC 1 or the environmental variable. This dummy variable has practically no influence on cluster analysis results. The same artifice can be used for cluster analysis based on mean yield of locations (complementing joint regression analysis).

As anticipated (Section 4.5), IRRISTAT's Regression module also allows for the performance of a stepwise multiple regression analysis, which may be useful in various instances (see Sections 5.4 and 5.8). In addition, IRRISTAT can be used for the extensive worksheet calculation that may be needed for up-scaling spatially and temporally the genotype adaptive responses according to relevant equations (e.g. [5.1], [5.4], [5.5], [5.7]). In particular, these equations may be calculated for all relevant genotype-location combinations (if necessary, with only a subset of better-performing genotypes) included in a large IRRISTAT data file. It is possible to assign to each genotype-location combination through logical operators estimated adaptation parameters of genotypes (mean yield, together with joint regression coefficient and intercept value, scaled scores on significant PC axes, or factorial regression coefficients and intercept value) and long-term values of mean yield or environmental variables of location that may be needed directly (for joint regression and factorial regression models) or for estimating site scores on significant PC axes (for AMMI models). Likewise, extensive worksheet calculations may be used for classifying new sites or test sites into subregions for specific breeding based on long-term values of the environmental variables that are required for predicting either the site score on PC 1 (according to the main crossover criterion), or the site value of the dummy variable for group discrimination (following the classification by AMMI + cluster analysis or by pattern analysis).

Other software

Various computer programs have been developed mainly or exclusively for the analysis of GE interaction effects. Their use in the analysis of GL interaction effects in trials repeated in time mostly requires the input of genotype by location cell means and, when needed, the manual calculation of F tests for components of the GL interaction using the appropriate error MS. Genotype by location by year cell means may sometimes be inputted as if they were genotype by location by replicate (i.e. plot) values.

GEBEI is a program (available free of charge) for pattern analysis developed by Watson et al. (1996) (University of Queensland). Most of its utilities have been included in IRRISTAT's Pattern Analysis module. GEBEI, however, also allows for inputting a similarity matrix of phenotypic correlation coefficients for locations, which is transformed by the program into a dissimilarity matrix of squared Euclidean distances prior to site (or genotype) classification or ordination. This opportunity can be of special interest for classifying locations based on information from historical data sets (see Section 5.5).

EIGAOV (Cornelius et al., 1996) is a program written by Dr P.L. Cornelius (University of Kentucky) for modelling GE interaction effects by joint regression, AMMI and shifted multiplicative models. It provides estimates of model parameters and tests them for significance by different criteria. In this context, it also allows for inputting the value of an appropriate error term obtained from a combined ANOVA (and expressed on a cell mean basis). Expected values for the selected model(s) can be printed and stored in an output data file. The main constraint is the low user-friendliness compared to other specific software. At present, EIGAOV is available free of charge for public institutions.

MATMODEL Version 2.1 (Gauch and Furnas, 1991; Gauch, 1998) is written by Dr H.G. Gauch (Cornell University) for modelling GE effects by joint regression and AMMI analysis. It can estimate missing cell means by an appropriate algorithm (Gauch and Zobel, 1990). In comparison with IRRISTAT's G × E Interaction module, MATMODEL offers additional opportunities for: i) testing PC axes by a cross-validation procedure (Gauch, 1992^[28]); ii) obtaining (for the selected AMMI model and in combination with the freely-available program, AMMIWINS, which utilizes its output) a concise summary of subregion definition according to the winning genotype criterion and the yield gain expected from growing the top-ranking genotype in each subregion relative to the top-ranking genotype over the region. The use of the program for trials repeated in time is discussed by Annicchiarico (1997b).

Finally, INTERA Version 3.3 (Decoux and Denis, 1991) is a program developed by INRA's Laboratory of Biometrics at Versailles that allows for a comprehensive investigation of GE interaction by joint regression, factorial regression and AMMI analysis (called "Mandel's model"). The program, which can estimate missing cell means, includes a method for classifying locations and genotypes, as well as procedures for performing ANOVAs with location or genotype group factors and producing graphs. Extensive documentation (also in English) is available for the user, who must however have a relatively good background in statistical methods.

The GENSTAT software provides users with specific procedures for joint regression, factorial regression and partial least squares regression modelling. In addition, a procedure for AMMI analysis (written by M. Talbot, K. Brown and M.F. Smith) is distributed informally and will be included in the next version of the software (R.W. Payne, personal communication, 2001). These procedures relate to analysis of GE interaction effects and may, therefore, require some modification to the recommended instructions when analysing GL interaction effects for trials repeated in time. Other powerful statistical software offers fewer or no specific procedures for data analysis. For SAS users, however, the analysis of GL effects by joint regression and AMMI models, as well as additional information (e.g. estimation of heterogeneity of genotypic variance and lack of genetic correlation components of GE interaction), can be obtained for balanced data sets by the program (intended as a set of SAS commands and procedures) developed by Hussein et al. (2000). Instructions are given by Denis et al. (1997) for factorial regression and by Piepho (1999) for joint regression analysis of GE effects also for unbalanced data sets. In addition, instructions for GE interaction analysis are given by Crossa et al. (1993) for the shifted multiplicative model using the SAS software, and by van Eeuwijk et al. (1995) for the reduced rank factorial regression model using GENSTAT.