5. Analysis of adaptation and identification of subregions

5.1 Objectives of the analysis

Techniques for analysis of adaptation have been developed for three main objectives:

• definition of adaptation strategies for breeding programmes;

• variety recommendation; and

• identification of optimal test or selection locations.

Different objectives may co-exist in the analysis of one data set, but partially different approaches may be required. Differences between the first and second objectives concerning the relevant GL interaction effects and the underlying criteria for subregion definition are outlined in Chapter 2 and summarized in Table 5.1. Different procedures are required for subregion definition within each analytical technique; consequently, results for site classification generally differ between the two objectives. In particular, location classification based on the similarity between relevant GL effects is essential for the preliminary definition of distinct subregions for breeding. On the other hand, accurate modelling of best-yielding material across locations is required for the definition of subregions for genotype recommendation. Data transformation to compensate for heterogeneity of genotypic variance among sites is useful for definition of adaptation strategies (see Section 5.6). The indications for this objective in terms of relevant GL effects and data transformation are also useful for identifying optimal test or selection locations (see Section 6.2), although subregion definition may not be implied in this case.

TABLE 5.1 - Major differences between two possible objectives for analysis of adaptation

^a Related to lack of genetic correlation among locations.
^b Complemented by comparison of wide vs. specific adaptation.
^c Against heterogeneity of genotypic variance among sites.

Scope for modelling

Modelling GL interaction effects through joint linear regression, AMMI (additive main effects and multiplicative interaction) factorial regression etc. is justified, as a simplified description of genotype responses helps target genotypes and understand adaptation patterns. Modelling is also important for improving the prediction of future responses of genotypes in comparison with observed data, because the portion of uncontrolled error variation in the GL interaction (so-called "noise") tends to be removed from the meaningful "pattern" portion retained in the model (Gauch, 1990, 1992). The higher predictive ability of expected (i.e. modelled) yields relative to observed yields for genotype-location combinations can be justified by the fact that model parameters (hence, all plot values used for estimation of these parameters), rather than only plot values for the specific genotype by location cell mean, concur to the estimation. In fact, genotype ranks on individual sites based on observed data can be quite different from those expected according to the selected AMMI model (Crossa et al., 1990; Gauch and Zobel, 1997; Ebdon and Gauch, 2002). Differences often concern also the best-yielding genotypes on the sites; these "winning" genotypes become a smaller subset based on expected data (thereby facilitating cultivar recommendation), because noise effects allowing lower-yielding material to occasionally win have been removed (Gauch and Zobel, 1997). For trials repeated in time, the noise portion relates mainly to the error term for GL interaction, i.e. non-repeatable GL interaction or GY interaction within locations. On the basis of results from an independent data set (Annicchiarico et al., unpublished results), when recommendation of the two top-yielding cultivars in each location was based on yield response modelled by factorial regression, joint regression or AMMI analysis (as reported in the case study in Section 8.3), yield gains were on average 4 to 5 percent higher than when recommendation was based on observed yields for each site. The elimination of noise effects from genotype responses - of crucial importance for making variety recommendation - may also be contemplated when grouping sites on the basis of their similarity for these responses in the context of subregion definition for breeding.

Analytical approaches

The following techniques - namely, joint regression model, AMMI models, factorial regression models including environmental covariates and pattern analysis - are a subset of those which may be adopted for analysis of adaptation. Major techniques not described here include:

• shifted multiplicative model (Cornelius et al., 1992; Crossa et al., 1993, 1995);

• application of canonical variates analysis (Seif et al., 1979; Calinski et al., 1987);

• factorial regression model contemplating genotypic as well as environmental covariates (Denis, 1988; van Eeuwijk et al., 1996);

• reduced rank factorial regression (van Eeuwijk, 1992, 1995);

• partial least squares regression (Aastveit and Martens, 1986; Talbot and Wheelwright, 1989; Vargas et al., 1998); and

• a recent method based on a Bayesian approach for variety recommendation (Theobald et al., 2002).

Ordinary crop simulation models offer sometimes excellent potential (e.g. Hammer and Vanderlip, 1989), but their adoption is generally limited by the difficulty of reliably representing the GL effects by empirically-estimated genetic coefficients - also in view of the large amount of experimental data required (Hunt et al., 1993; Russell et al., 1993). Recent models incorporating gene action may be considered an integration rather than an alternative to current models for analysis of adaptation, mainly for assessing the value of single traits and combinations of traits in the context of a given adaptation target (Hoogenboom et al., 1997; Chapman et al., 2002). However, some models may also support decisions on adaptation strategies, by simulating results of multi-environment trials and even the outcome of different selection strategies (Cooper et al., 1999).

5.2 Joint linear regression modelling and complementary analyses

Joint linear regression analysis was developed by Yates and Cochran (1938). Slightly different models have since been proposed by Finlay and Wilkinson (1963), Eberhart and Russell (1966), Perkins and Jinks (1968) and others (as reviewed by DeLacy et al., 1996a). According to Perkins and Jinks' model (applied to GL rather than GE interaction analysis), the GL_ij effects (relative to the genotype _i and the location j) are modelled as a function of the location main effect (L_j) or the location mean value (m_j), which represents an indicator of the ecological potential of the site for the crop (positive L_j or high m_j = high potential; negative L_j or low m_j = low potential):

GL_ij = β_i L_j + d_ij = α_i + β_i m_j + d_ij

where β_i is the regression coefficient of the genotype _i and d_ij is the deviation from the model, i.e. the residual GL interaction; in the second expression an intercept value α_i (equal to -m β_i, where m = grand mean) is also present. The β coefficients, with a mean value equal to zero, and the genotype mean yields are the relevant estimated parameters of genotype adaptation, since the genotype merit on a given site depends on its mean yield and the expected GL interaction effect (which varies according to the β coefficient). Largely positive and negative β values, if associated with relatively high mean yield, result in specific adaptation to high-yielding (favourable) and low-yielding (unfavourable) sites, respectively. Conversely, β values around zero indicate a lack of specific adaptation (and wide adaptation, if combined with high mean yield). No definite indication on genotype adaptive responses can be inferred solely from β values: for example, an entry that is highly susceptible to a disease occurring on high-yielding sites would show a negative β value even if it was only mid-ranking for yield on low-yielding sites.

Estimation and test of model parameters

In the joint regression analysis table, the GL interaction SS (relative to ANOVA models in Tables 4.1 through 4.3) is partitioned into two components relating, respectively, to:

• heterogeneity of genotype regressions; and

• deviations from regressions.

The SS proportion accounted for by heterogeneity (model R²) may be obtained by summing up the model SS across separate regression analyses of GL effects performed for each individual genotype, and calculating the proportion in the total SS in the regressions. The complement to one of this value indicates the proportion of GL variation accounted for by deviations from regressions. As an alternative, SS values for genotype regressions and deviations from regressions can be provided by model and error SS, respectively, in the analysis of covariance of all GL effects as a function of genotype main effect and genotype × site mean yield interaction (which also provides estimates of β_i and α_i values). Appropriate software can simplify the analysis (see Section 5.9). For integration with ANOVA results in a joint regression analysis table, the SS for the two terms obtained from analysis of a genotype by location cell mean basis must be multiplied by N' = no. test years (or crop cycles) x no. experiment replicates. They may also be calculated as the respective proportion of GL interaction variation multiplied by the ANOVA GL interaction SS. The GL interaction DF can be partitioned into:

• a (g - 1) portion for heterogeneity of regressions among g genotypes; and

• a remainder for deviations from regressions.

MS values for these terms can then be calculated and reported in a joint regression analysis table (e.g. Table 4.4.). The genotype regressions term should preferably be tested for significance using an F ratio with deviations from regressions MS as the error term (Cruz Medina, 1992). The deviations from regressions MS are then tested for significance using the appropriate error term for overall GL interaction in the ANOVA (GLY interaction in Table 4.4). Likewise, it is preferable to test each β coefficient for difference to zero (i.e. presence of interaction with sites of the genotype) using the deviation from regression of the individual genotype as the error term (Freeman, 1973), as provided by separate regression analyses for the individual genotypes (testing of β values by analysis of covariance refers, conversely, to an average value of deviations).

Finlay and Wilkinson (1963) assessed MS values for genotype regressions and deviations from regressions with the same procedure as Perkins and Jinks' (1968) and described the response to site mean yield of the genotype _i by the coefficient b_i as equal to:

b_i = β_i + 1

Therefore, b values close to unity are indication of the negligible interaction with sites of the genotype; testing for difference to unity is equal to testing β for difference to zero.

Eberhart and Russell's (1966) b coefficients for genotype regressions are the same as Finlay and Wilkinson's. However they proposed a different and less convenient layout of sources of variation in the analysis (Freeman, 1973; Westcott, 1986).

Modelled yield responses and cultivar recommendations

According to Perkins and Jinks' model, the expected yield response (R_ij) of the genotype _i at the location j can be expressed as a function of L_j or m_j as follows:

R_ij = m + G_i + L_j + β_i L_j = G_i + m_j + α_i + β_i m_j

where m = grand mean and G_i = genotype main effect. The same response can be described according to Finlay and Wilkinson's b coefficient as follows:

R_ij = m + G_i + Lj + (b_i - 1) Lj = a_i + b_i mj

[5.1]

The latter simple expression is particularly useful. The intercept values a_i (different from previous α_i values) are equal to:

a_i = G_i - m β_i = m_i - m b_i

where m_i = genotype mean yield. Therefore, using results generated by genotype regressions for GL effects (needed for calculation of SS and MS values and relative to β_i and α_i parameters), a_i and b_i may be estimated, and b_i may be tested for difference to unity analysis. As usual in regression analysis, the inference on expected yields is valid in the range of observed site mean yields. Non-significant parameters should also be retained in the model as they provide the best estimation of yield responses.

One last expression of genotype responses that can be useful (see Section 7.2) and will largely be considered in AMMI and factorial regression modelling, is response in terms of nominal yields, i.e. expected yields from which the main effect of location (which has no effect on ranking of genotypes on each site) has been eliminated. The nominal yield (N_ij) of genotype i in location j can be estimated as follows from the genotype mean value (m_i) and Perkins and Jinks' regression parameters, and the site mean value:

N_ij = m_i + α_i + β_i m_j

[5.2]

Adaptation patterns of best-yielding genotypes across the relevant range of site mean yields in Annicchiarico and Mariani (1996) are reported in Figure 5.1 (A). For recommendation purposes, the responses suggest the presence of three mega-environments:

• high-yielding (site mean yield > 4.5 t/ha), in which 'Karel' is the top-ranking but 'W 4267' and 'M 104' yield almost as well;

• medium-yielding (mean yield 3-4.5 t/ha), in which 'Messapia' is the best-ranking genotype; and

• low-yielding (mean yield < 3 t/ha), in which 'D 3415' appears preferable.

Figure 5.1 - Genotype yields modelled by joint regression (A) and nominal yields modelled by factorial regression on rainfall amount (B), for durum wheat varieties across Italian locations

Source: Annicchiarico and Mariani, 1996

Recommending (for the sake of biodiversity) at least two varieties per subregion would imply the addition of a second best-ranking genotype in each case, with a possible slight change in the yield boundary of subregions (e.g. 'D 3415' and 'Messapia' could be recommended in a low-yielding subregion with mean yield < 3.3 t/ha).

It is possible to test the significance between pairs of genotypes at given values of location mean yield (e.g. Pecetti and Annicchiarico, 1993). However, such tests are time-consuming without powerful statistical software. An outline, rather than a precise assessment, of the extent of significant differences between genotypes at variable levels of site mean yield can be provided by the mean value of Dunnett's one-tailed critical difference. Adopting formula [4.2] (Section 4.2), the appropriate error MS (M_err) for trials repeated in time is here represented by the average GY interaction within locations - in the absence of significant deviations from regressions (i.e. when joint regression is a fully reliable model). This term is readily available from ANOVAs in which the time factor is nested into location (Table 4.3) and is easily obtained by pooling GY and GLY interactions (dividing the sum of SS by the sum of DF) for ANOVAs with time and location crossed factors (Table 4.2). If the deviations from regressions term is significant, it can be pooled with the GY interaction within locations, giving rise to a pooled error term comprising all GE interaction variation not accounted for by heterogeneity of genotype regressions. This leads to a certain loss of protection for the tests. With respect to data in Figure 5.1 (A), by pooling SS and DF values of GY and GLY interactions and deviations from regressions (Table 1 in Annicchiarico and Mariani, 1996), M_err = 0.60. For P < 0.20 and p = 8 (cultivars are nine), t' ≈ 1.71 (from values in Table 4.6, which provide a good approximation for the pooled error DF = 128), whereas N (relating to comparisons in individual locations) can be equalled to: no. years × no. replicates = 9. The mean critical difference is:

d ≈ 1.71 √ (2 × 0.60/9) = 0.62 t/ha

This value, superimposed on yield responses in Figure 5.1 (A), would show that the genotype 'Messapia' is hardly ever significantly less yielding than any other genotype within the range of site mean yields. For ANOVAs with no time factor (Table 4.1), an appropriate error MS may be based on the the pooled experimental error (pooled with the deviations from regressions term, if significant). If there is proportionality between site mean yield and the appropriate error (i.e. GY interaction or pooled experimental error within individual locations), the mean critical difference is too large for comparisons on low-yielding sites and too small in high-yielding locations. Possible solutions for increasing the reliability of genotype comparison are discussed in Section 5.6.

Location classification for breeding

If genotype adaptive responses depend substantially on the location mean yield, breeding for specific adaptation is possible in a few subregions characterized by contrasting levels of site mean yield. Classification of test locations may derive from multiple comparisons between sites in the context of the ANOVA, or from a cluster analysis including site mean yield as the only variable. Of the different clustering methods available (Everitt, 1980), a hierarchical clustering strategy - using either Ward's incremental sum of squares or the average linkage (i.e. group average) methods and holding a squared Euclidean distance as the dissimilarity measure - is recommended (Williams, 1976b; DeLacy et al., 1996a). The P level for significance of the GL interaction within group of locations (here related to variation for mean yield, as assumed by the model) can provide a simple and convenient truncation criterion for definition of groups of locations (Ghaderi et al., 1980). Going backwards from the last fusion stage to the formation of groups, an ANOVA performed for each newly-formed group on data from the member locations would prevent the further subdivision of the group when GL interaction ceases to be significant at a specified P level (P < 0.05 or P < 0.01). An application of this procedure is exemplified in the next section. Ward's clustering method, which minimizes the within-group SS at each clustering stage, is particularly appropriate for adoption in combination with this truncation criterion.

An alternative criterion of location classification for breeding purposes in the context of joint regression analysis was proposed by Singh et al. (1999). It contemplates the estimation of the value of site mean yield for which crossover interactions between genotypes (whether between high- or low-yielding entries) reach the highest frequency. This main crossover point serves as a cut-off for defining two groups of locations. If Finlay and Wilkinson's regression coefficient = b_i and the intercept value = a_i for the genotype _i, the main crossover point (X_co) may be estimated as:

X_co = - ∑ a_i (b_i - 1)/ ∑ (b_i - 1)²

[5.3]

While this criterion distinguishes only two provisional subregions, the criteria based on cluster analysis imply two or more. Another more complex technique, described by Crossa and Cornelius (1997), involves classifying sites into three or more groups on the basis of crossover GL effects.

Concluding remarks

For both breeding and genotype recommendation, the assignment of new locations to subregions depends on the level of long-term site mean yield. Similarly, test locations can be re-assigned to subregions on the basis of their long-term yield.

Due to its simplicity, the joint regression model has been the most popular approach for analysis of adaptation (Becker and Léon, 1988; Romagosa and Fox, 1993). The method has some statistical limitations. Caution should be applied for low numbers of genotypes and locations, especially when extreme values of site mean yield are represented by just one location (Westcott, 1986; Crossa, 1990). But the main drawback is its limited ability to describe genotype adaptive responses in two notable situations:

• when these responses are not ecologically simple, for example, in the presence of a complex of variably occurring environmental constraints (implying multidimensional adaptation patterns, owing to the different adaptive genes or sets of genes that are possibly involved in the plant responses); and

• when GL effects and location mean yield are affected mainly by a different environmental factor (Annicchiarico, 1997a; Brancourt-Hulmel et al., 1997).

Furthermore, non-linear genotype responses to site mean yield cannot be taken into account, unless more complex analytical models are adopted (e.g. Pooni and Jinks, 1980; Mariani et al., 1983). The possible criteria for comparing different analysis of adaptation models in terms of their ability to describe GL interaction are discussed at the end of Section 5.4. In all cases, an adequate model R² should be relatively high (e.g. > 30%).

5.3 AMMI modelling and complementary analyses

Work by Williams (1952), Gollob (1968), Mandel (1971) and Bradu and Gabriel (1978) has made an important contribution to the development of AMMI models. Their application to agricultural research was proposed by Kempton (1984) and Zobel et al. (1988), but their use became widespread thanks to the comprehensive monograph by Gauch (1992). The potential field of application of this technique in agriculture goes far beyond the study of GE interactions (Gauch and Zobel, 1996a).

In analysis of adaptation, AMMI analysis requires at first the estimation of genotype and location main effects by ANOVA. Residuals from additivity of these effects (i.e. GL effects) are then partitioned into:

• the multiplicative term of the model, of which the estimated parameters relate to the statistically significant axes of a double-centred principal components analysis performed on the GL interaction matrix; and

• a deviation from the model term:

  GL_ij = ∑ u_in v_jn l_n + d_ij = ∑ (u_in √ l_n) (v_jn √ l_n) + d_ij

where u_in and v_jn are eigenvectors (scaled as unit vectors, i.e. ∑ u_in² =∑ v_jn² = 1) of the genotype _i and the location j, respectively, and l_n is the singular value (i.e. the square root of the latent root or eigenvalue) for the principal component (PC) axis n; and d_ij is the deviation from the model.

The proportion of GL interaction variation accounted for by each PC axis is equal to the relative size of its eigenvalue. The further scaling of eigenvectors through multiplication by √l_n allows for a straightforward estimation of the GL_ij effects expected on the PC axis n by multiplication of the scaled genotype and location scores on that axis.

There are several possible AMMI models characterized by a number of significant PC axes ranging from zero (AMMI-0, i.e. additive model) to a minimum between (g - 1) and (l - 1), where g = number of genotypes and l = number of locations. The full model (AMMI-F), with the highest number of PC axes, provides a perfect fit between expected and observed data. Models including one (AMMI-1) or two (AMMI-2) PC axes are usually the most appropriate where there is significant GL interaction. Due to their simplicity, they provide a notable reduction of dimensionality for the adaptation patterns relative to observed data.

While principal components analysis is usually executed on the correlation matrix (Dagnelie, 1975b; Joliffe, 1986), for AMMI modelling it is executed on the covariance matrix. Furthermore, two (not one) analyses are performed simultaneously: in the analysis the genotypes are individuals (rows) and the locations original variables (columns); in the other, vice versa.

For example, Figure 5.2 shows the scaled scores of three genotypes (left) and three locations (right) in the space of the first two PC axes, for the AMMI analysis of data reported in Table 4.5:

• The first axis (PC 1) maximizes the variation of GL effects in one dimension (i.e. it minimizes the sum of squared projections of the points off that axis).

• The second axis (PC 2) maximizes the residual variation in a second dimension that must be perpendicular to PC 1 (correlation is therefore zero between PC 1 and PC 2 scores, for both genotypes and locations).

Figure 5.2 - AMMI analysis of the genotype-location data matrix in Table 4.4

Note: Scaled scores of genotypes (G) and locations (L) in the space of the first two PC axes.

For a larger GL interaction matrix, the cloud of points on the graph relating to genotypes or locations could be represented by an ellipse of concentration, where the longer axis is proportional to the variation (as standard deviation) for PC 1 and the shorter axis to the variation for PC 2 scores (Fig. 5.2). Any additional PC axis should be perpendicular to the PC axes already in the model. PC 1 represents mainly a contrast of: i) genotype 3 with the other genotypes; and ii) location 3 with the other locations (Fig. 5.2).

Genotype 3, which has the same sign as location 3, is expected to show a large positive GL effect in this location according to this PC axis (as estimated by the product of the genotype and location PC 1 scores). If PC 2 (representing mainly a contrast of genotype 1 with the other genotypes, and of location 1 with the other locations) is also kept in the model, the slight negative GL interaction effect of genotype 3 in location 3 (as estimated by the product of the respective PC 2 scores) should also be taken into account. A large negative GL effect could be estimated for genotype 2 in location 3 and for genotype 3 in location 2, owing to scores with different signs on both PC axes (Fig. 5.2). These effects reflect the observed GL effects reported in Table 4.5.

The ordination of genotypes and locations according to their scaled scores on the first two GL interaction PC axes is sometimes reported in a single graph (known as biplot) to help appreciate location and genotype similarity for GL and allow for the graphical estimation of these effects. In particular, the closeness between pairs of locations or pairs of genotypes in the biplot is proportional to their similarity for GL effects. Genotypes represented by a point near the origin of axes reveal limited GL interaction. On the contrary, material that is far from the origin has a much better response to locations which are far from the origin and for which the angle formed between the genotype point, the origin and the location point is small. Locations far from the origin show rather peculiar responses of genotypes. For example, the scores on the first two PC axes of 18 genotypes and the mean values of scores for four groups of locations (A, B, C and D) are reproduced in Figure 5.3 (Annicchiarico and Perenzin, 1994). PC scores for the 31 locations included in this study (which would make the current representation rather confusing) were indicated in another graph in which results of site classification by cluster analysis were also superimposed (as in Fig. 5.5 [B] for another set of locations). While highlighting the pattern and extent of interaction effects between genotypes and subregions (or locations), this graph does not allow for a definite appreciation of genotype adaptation, as it provides no information on the other relevant parameter in this context: genotype mean yield. When there is only one significant GL interaction PC axis, the allocation of mean values on the abscissa axis and PC 1 scores on the ordinate axis of the biplot helps to appreciate both determinants of genotype performance (Gauch, 1992^[12]; Gauch and Zobel, 1997). However, other graphs (proposed by these authors and described later on in this section) can provide a straightforward representation of genotype adaptive responses, thereby facilitating the identification of recommendation domains on the basis of best-yielding material.

Figure 5.3 - Scores on the first two GL interaction PC axes of 18 bread wheat genotypes (numbers) and four Italian subregions (letters)

Source: Annicchiarico and Perenzin, 1994.

Estimation and test of model parameters

For the definition of an AMMI analysis table, the GL interaction SS in the ANOVA is divided into portions relating to each significant PC axis and to a residual term. The SS for each PC can be obtained as the proportion of GL interaction variation accounted for by its eigenvalue (which is unaffected by the definition of genotypes or locations as individuals in the analysis) multiplied by the ANOVA GL interaction SS. According to Gollob (1968) the DF for the PC axis n can be calculated as:

DF = g + l - 1 - 2 n

where g = number of genotypes and l = number of locations. The GL interaction SS and DF not accounted for by significant PC axes are pooled in the residual term. An MS can then be calculated for each PC and the residual. The procedure may be verified for data in Table 4.4. The proportion of variation explained by eigenvalues is 0.22 for PC 1 and 0.16 for PC 2, while the model R² is 0.22 + 0.16 = 0.38 = 38%. For example:

• PC 2 SS is 0.16 × 587.9 = 94.1

• PC 2 DF is 18 + 31 - 1 - 4 = 44

Statistical testing of PC axes has mainly been investigated for analysis of a GE interaction matrix (applicable to GL interaction analysis for trials not repeated in time). An early F test devised by Gollob (1968) is too liberal (biased towards too many significant results) both in theory and on the basis of simulation results (Mandel, 1971; Cornelius, 1993). Cornelius (1993) proposed the F_GH2 test, using approximations published in an earlier paper (Cornelius, 1980), as well as the F_GH1 test and some simulation tests that are expected to produce results very similar to the F_GH2 test while requiring more extensive calculation. In addition, Cornelius et al. (1992) described an alternative test: the F_R ratio. Simpler to calculate, F_R is more robust in the presence of heterogeneous within-site experimental errors than the F_GH2 test (Piepho, 1995). Gauch (1988, 1992^[13]) proposed a further test criterion based on a cross-validation procedure. The available plot values (replicates) for each genotype-environment combination are randomly split into modelling and validation data, and the possible AMMI models (from AMMI-0 to AMMI-F) are compared in terms of predictive ability (i.e. the ability to explain variation in data not used in constructing the model itself). Following assessment, the selected AMMI model is reparametrized using information from all data. For trials in randomized complete block designs, Piepho (1994a) recommends splitting complete blocks within environments rather than single replicates, in order to avoid the bias towards too few significant PC axes caused by an error term inflated by the block effect. Nevertheless, this test tends to be too conservative because modelling relies on a subset of the available data, which is less accurate than modelling all data (as performed after model choice) and is therefore somewhat biased towards simpler, lower-ranking AMMI models (Cornelius, 1993). Finally, Gauch (1992^[14]) suggests a simple guideline for model selection based on the estimation of the noise level expected in the GE interaction matrix.

The application of most of the above tests to AMMI analysis of GL interaction for trials repeated in time is discussed by Annicchiarico (1997b). For cross-validations, two procedures are proposed, aimed at eliminating undue sources of variation from the relevant error term in the context of ANOVA models in Table 4.2. For the application of the remaining criteria it is necessary to utilize the appropriate error term for GL interaction - i.e. GLY interaction (Table 4.2) or GY interaction within locations (Table 4.3) - rather than the pooled error term (which is applicable in the absence of repetition in time). Although the preferred testing procedures may vary according to the size of the GL matrix and the number of test years (Annicchiarico, 1997b), the adoption of the F_R test can be recommended in a wide range of situations (also because of its simplicity and robustness).

The F_R test verifies the significance of the residual GL interaction variation in each AMMI model, beginning with AMMI-0. By an ordinary F ratio, the MS of the residual is tested on the same MS used as the error term for the GL interaction in the ANOVA. A significant result implies the addition of one more PC to the model (i.e. the significance of the newly-added PC). The test for PC 1, in which the residual is represented by the entire GL interaction, coincides with the ANOVA F test for the interaction. For example, consider the test for PC 2 in Table 4.4. For the residual GL interaction from AMMI-1:

• SS is 587.9 - 129.3 = 458.6.

• DF are 510 - 46 = 464.

• MS is 458.6/464 = 0.99.

• F ratio (using the appropriate error term) = 0.99/0.78 = 1.27 (with 464, 1020 DF) - higher than the critical F value at P < 0.01 (≈ 1.20).

Therefore, PC 2 is also admitted in the model. The main drawback of the F_R test is that no PC 1 can be declared significant in the absence of significant GL interaction, even if it accounts for most of the GL interaction SS. In such a rare event, reported by Zobel et al. (1988) for one data set, the signal-to-noise ratio criterion is prone to the same drawback; more complex tests, such as F_GH2, should therefore be envisaged.

AMMI analysis is aided by adopting customized statistical software (including IRRISTAT, see Section 5.9). When such software is not available, two separate principal components analyses must be executed on the GL interaction matrix. The analyses consider genotypes or locations in turn as individuals and original variables. Location and genotype scaled scores on the PC axis n can be obtained from their respective eigenvectors issued by the two analyses, multiplied by the square root of the relevant singular value (√ l_n). Alternatively, the scaled scores can be obtained from the original scores of locations and genotypes on PC n, dividing by √ l_n. In both cases, ln can be calculated from the eigenvalue ln², which is equal to the GL interaction SS accounted for by PC n in the AMMI analysis, expressed on a genotype by location cell mean basis by dividing by N' (where N' = no. test years [or crop cycles] x no. experiment replicates). Considering 10 genotypes grown in 20 locations in trials with 3 replicates conducted over 4 years, if the ANOVA GL interaction SS is 256.9, and the proportion of GL variation explained by PC 1 is 0.274:

l₁² = (0.274 × 256.9)/(4 × 3) = 5.866

√ l_n = (l₁²)^0.25 =1.556

location (or genotype) PC 1 scaled score = location (or genotype) eigenvector × 1.556

location (or genotype) PC 1 scaled score = location (or genotype) original score/1.556

The last of these equations may also be used to obtain original scores of locations, of possible interest for site classification performed by cluster analysis (discussed below), when using specific software providing scaled PC score values.

The performance of two independent principal components analyses the additional problem of a potentially wrong sign of the genotype or the location eigenvectors, assigned independently in the two analyses. The polarity of eigenvectors and PC scores is arbitrary in AMMI analysis, but pairs of genotype and location scores must possess the same sign (positive or negative) to produce a positive GL interaction effect, and different signs to produce a negative GL effect. Therefore, signs of eigenvectors and PC scores may need to be inverted for genotypes (or locations), following the comparison between modelled and observed GL interaction effects. Several genotype-location pairs possessing the highest scores (in absolute terms) on the PC axis reveal the sign of the observed GL effect. The correct sign provides the best match between modelled and observed effects.

Modelled yield responses and cultivar recommendations

The expected yield response (R_ij) of the genotype _i in the location j is (according to the selected AMMI model):

R_ij = m + G_i + L_j + ∑ (u_in √ l_n) (v_jn √ l_n) =m + G_i + L_j + ∑ (u_in' v_jn')

where m is the grand mean, G_i and Lj are the main effects for the genotype _i and the location j, u_in' = (u_in √ l_n) indicates the scaled score on the axis n for the genotype and v_jn' = (v_jn √ l_n) that for the location. Since changes in genotype rank across locations only depend on the multiplicative term, the adaptive responses can conveniently be represented as a function of the scaled scores of locations on the statistically significant PC axes. The location main effect has no influence on genotype value and it complicates the graphic representation of adaptation patterns. Gauch (1992) and Gauch and Zobel (1997) therefore proposed to eliminate this effect from modelling by representing the expected yield responses in terms of nominal yields. This has the additional advantage of simplifying the calculations. For an AMMI-1 model, the nominal yield (N_ij) of the genotype _i in a given location j can be estimated from the genotype mean value (m_i) and the scaled scores of genotype and location on PC 1:

N_ij = m + G_i + (u_i1 √ l_l) (v_j1 √ l_l) =m_i + (u_i1' v_j1')

[5.4]

The nominal yield responses of genotypes may be represented by straight lines as a function of location PC 1 scores reported in abscissa. This is shown in Figure 5.4 (A) for seven lucerne varieties grown in ten locations in northern Italy (the lowland sites in Annicchiarico [1992]). For each genotype, nominal yields are simply calculated for the two locations with extreme PC 1 score values (sites 'p' and 'o') and the two values connected by a straight line. On the basis of adaptive patterns:

• 'La Rocca', 'Prospera' and (to a lesser extent) 'Robot' can be recommended for the area represented by sites 'p' and 'l';

• 'Prosementi' and (to a lesser extent) 'Garisenda' and 'Robot' are preferable in a vast area, including sites 'z', 'a', 'b', 's' and 'e'; and

• 'Orchesienne', 'Europe' and (to a lesser extent) 'Prosementi' are the best-performing in the area comprising sites 'u', 'o' and 'c'.

Figure 5.4 - Nominal yields of lucerne varieties modelled as a function of the score on the first GL interaction PC axis of ten Italian locations (A), or the first GE interaction PC axis of four artificial environments (B)

Note: Site classification based on cluster analysis of site scores on PC 1 (see Fig. 5.7 for geographic position of coded locations).

Source: Annicchiarico, 2002.

These subregions are similar to those summarized in Figure 5.4 (A) representing the cluster analysis site classification based on GL effects for all cultivars (discussed below also in geographical terms). Incidentally, in subsequent trials, recommendations allowed for higher yields (6% in the first and 3% in the third subregion) than those widely recommended (on the basis of average yields of genotypes across sites) (Annicchiarico, 1998). The extension of results to new locations is less straightforward for AMMI modelling than for joint regression, because it cannot be based on a definite site characteristic, such as mean yield. In general, the response of a new site is expected to be similar to that of the closest test site, on the basis of a supposedly close relationship between similarity for GL effects and geographical proximity. However, additional information on environmental factors related to GL interaction may greatly facilitate the characterization of subregions and the spatial and temporal scaling-up of results (see Section 5.8, in relation to subregion definition for genotype recommendation and breeding).

Figure 5.5 - Scores on the first two GL interaction PC axes of 20 Italian locations (letters), definition of four subregions for bread wheat variety recommendation by grouping sites with the same expected winning genotype (A), and classification of locations into five groups based on cluster analysis of site scores on PC 1 and PC 2 (B)

An indication of the extent of statistical differences between varieties on varying sites score on PC 1 in a graph like Figure 5.4 (A) can be provided by the mean value of Dunnett's one-tailed critical difference calculated in a similar way to joint regression. An appropriate error term is represented: by the average GY interaction within locations (obtained, for ANOVAs with the time factor crossed with location, by pooling SS and DF values of GY and GLY interactions) for trials repeated in time, and by the pooled experiment error for trials not repeated in time. The latter applies to Figure 5.4 (A), where: M_err = 4.09; t' ≈ 1.79 (Table 4.6) for P < 0.20 and p = 10 (eleven cultivars overall, of which seven are shown); N = no. replicates = 4; and mean critical difference d ≈ 1.79 √ (2 × 4.09/4) = 2.56 t/ha. This value, superimposed on yield responses in Figure 5.4 (A), would slightly enlarge the range of recommended material. The critical difference is potentially valid also for comparison of nominal yields in AMMI models of higher complexity.

For specific recommendation on each site also, indications based on genotype comparison at the conventional P < 0.05 rate of Type 1 error may have negative implications due to concurrently high Type 2 errors. For example, recommending the AMMI-modelled best-yielding set of varieties based on Dunnett's one-tailed critical difference at P < 0.05 (average of five varieties) has on average provided over 3 percent lower yields than merely recommending the top-ranking variety (for data sets in Table 3.1, using 4 years for modelling and 2 other years for verification of yield gains - Annicchiarico, unpublished results).

For an AMMI-2 model, nominal yields can be estimated from the genotype mean value (m_i) and the product of genotype and location scaled scores on PC 1 and PC 2:

N_ij = m_i + (u_i1' v_j1') + (u_i2' v_j2')

[5.5]

The graphical definition of mega-environments based on best-yielding material is complicated by the difficulty of representing the genotype responses as a function of site scores on two PC axes simultaneously. After representing the site scores as points in the bidimensional space of PC 1 and PC 2, a third dimension is required to represent the yield levels. However, the "winning" genotype expected for each location can be determined by calculating AMMI-2 nominal yields, and locations with the same best-yielding genotype can be grouped into the same subregion. For example, the definition of subregions and best-yielding material on the basis of repeated testing of 10 bread wheat varieties in 20 Italian locations (Annicchiarico et al., unpublished results) is reported in Figure 5.5 (A). Genotype '7' may be recommended on 12 sites, and cultivars '2', '3' and '4' on a small number of sites. Despite the small number, a specific recommendation (giving rise to a four mega-environment scenario) may be justified. The straight lines - subdividing the 20 coded locations into four subregions and defining boundaries in which the winning genotypes of two adjacent subregions have the same expected yield - can be calculated using appropriate software (see Section 5.9). The graph can also help in assigning new locations to subregions (see Section 5.8). However, it is possible to group sites merely on the basis of similarity of winning genotypes. If genotype recommendation considers more than one cultivar per subregion, sites can be grouped on the basis of the same two (or more) best-yielding genotypes. For estimation of nominal yields in AMMI-3 models, the product of genotype and location scaled scores on PC 3 is added to the terms in equation [5.5]. Estimated parameters for higher-order PC axes are also included for AMMI models of increasing complexity.

The present approach for modelling GL effects and defining subregions for genotype recommendation from trials repeated in time differs from that proposed by Gauch (1992^[15]), because it is based on the mean response of locations for GL effects. The results are thus simplified (e.g. the 20 points for locations in Fig. 5.5 [A] would become 80 for results expressed for individual environments). The average inconsistency of locations for GL effects (i.e. GLY interaction) is taken into account by the error term adopted for assessing the GL interaction effects. For a relatively low number of test sites, the results of an additional AMMI analysis performed on GE effects may reveal whether or not a location is characterized by unusual variation for PC scores between its environments, leading to inconsistency of the winning genotype between years (or crop cycles) (Fig. 6.5 in Gauch, 1992^[16]). However, a reliable assessment of the frequency of different winning genotypes in each location (i.e. the information required for recommendation in this context) may be difficult when based on the limited number of test years (< 3) normally available. This approach is difficult to apply when more than one recommended genotype is considered for each site.

Location classification for breeding

The provisional definition of subregions for breeding requires that site classification be based on similarity for overall GL effects (Table 5.1). The performance of cluster analysis, using as variables (i.e. dimensions in an Euclidean space) the scaled or unscaled scores of sites on the significant PC axes, means that the noise portion of GL interaction - which tends to be pooled into the non-significant, residual GL interaction term - may be substantially eliminated from the assessment. This procedure is consistent with one of the major objectives of ordinary principal components analysis, namely, summarizing the information provided by many correlated original variables (each one corresponding to the set of GL effects of genotypes in a given location) into a few uncorrelated PC variables, thus simplifying the multivariate comparison of the individuals (in this case, locations). In this context, the assessment of site similarity based on original rather than scaled PC scores of sites (e.g. Annicchiarico, 1992; Annicchiarico and Perenzin, 1994) may be preferable when there are two or more PC axes (while scaling does not affect assessment based on PC 1 scores only). Ward's incremental SS or the average linkage methods (with a squared Euclidean distance as the dissimilarity measure) can again be recommended for cluster analysis. When applied to PC axes, the dissimilarity measure for cluster analysis requires no prior standardization of variables, because the variation in site score on each PC axis is proportional to the importance of each variable (Dagnelie, 1975b^[17]; Afifi and Clark, 1984^[18]). As mentioned above, Ward's method is particularly suited to the truncation criterion based on the significance of the within-group GL interaction in ANOVAs.

For example, consider the cluster analysis applied to PC 1 scores of nine locations (all except site 'e') in Figure 5.4 (A). Results of site grouping for the final fusion stages are reported in Figure 5.6, together with significance levels for the within-group GL interaction according to ANOVA. P levels < 0.01 require further subdivision of groups, whereas higher P levels are here acceptable for group definition (as GL interaction is tested on a relatively small error term, i.e. the pooled experimental error). Three groups of locations may be defined in the analysis. They are superimposed (as subregions A, B and C) on the ordination of sites in Figure 5.4 (A). In the report by Annicchiarico (1992), location 'e' was actually represented by two environments clustered, respectively, in subregions B and C, implying a border position between these subregions. The results of site classification, combined with information on environmental factors related to the occurrence of GL interaction (Annicchiarico, 1992, 2000), have allowed for the geographical definition of potential subregions reported in Figure 5.7 - confirmed by AMMI and the cluster analysis results obtained for an independent set of subsequent trials (Annicchiarico, 2000). Predicted yield gains from different selection strategies, as well as opportunities for selection in artificial environments, are further discussed in Section 6.1.

Results of location classification for the bread wheat example based on: cluster analysis of PC 1 and PC 2 site scores as original variables; an average linkage method; and P > 0.05 for the within-group GL interaction as a truncation criterion for definition of groups of locations, are given in Figure 5.5 (B). In this case (unlike the lucerne example) the preliminary definition of mega-environments for breeding is distinctly different from that for variety recommendation (Fig. 5.5 [A]). The provisional subregions 1 and 2 are large enough to justify specific breeding (Fig. 5.5 [B]), whereas subregions 3 to 5 (represented by one location each) are too small and can be merged with subregion 1 (the closest). The site classification is geographically meaningful: the provisional subregion 2 groups locations in southern Italy, while the larger subregion 1 groups sites in northern and central Italy.

Figure 5.6 - Cluster analysis of test locations performed on site scores on the first GL interaction PC axis, using the lack of significant (P < 0.01) GL interaction within group of locations as the truncation criterion for definition of groups

Note: * = P < 0.05; *** = P < 0.001 (see Fig. 5.4 for PC scores and Fig. 5.7 for geographic position of locations).

Figure 5.7 - Geographic position of test locations for lucerne in northern Italy, and provisional definition of three subregions for a specific adaptation strategy based on AMMI analysis complemented by cluster analysis and additional information on relevant environmental variables of sites

Source: Annicchiarico, 2000.

For AMMI-1 models, an alternative criterion of site classification for breeding purposes can be provided by estimating the main crossover point for nominal yields of genotypes modelled as a function of the site scaled score on PC 1, to be used as a cut-off for the definition of two subregions. The procedure proposed is an extension to AMMI analysis of the procedure developed by Singh et al. (1999) for joint regression. The main crossover point X_co in formula [5.3] relates to genotype slopes which vary around a zero mean (b_i - 1 = β_i) and are different for intercept value a_i (i.e. nominal yields that have zero grand mean and are expressed as a function of site mean yield). A similar formula can be obtained for AMMI nominal yields expressed as in Figure 5.4, after performing a double change of origin that causes the zero value for the abscissa and the ordinate axes to coincide with the lowest site score on PC 1 and the grand mean, respectively.

The slope of the genotype _i as a function of PC 1 can be expressed as:

β_i = (N_i2 - N_i1)/(v₂₁' - v₁₁')

where j = 1 is the location with the lowest PC 1 score; j = 2 is the location with the highest PC 1 score; v₁₁' = (v₁₁ √ l₁) and v₂₁' = (v₂₁ √ l₁) are the scaled scores on PC 1 for these locations, as expressed in the graph; and N_i1 and N_i2 are the nominal yields of the genotype _i on the site with lowest and highest PC 1 score, respectively (these values are already calculated for the graphical expression of adaptive responses).

The main crossover point expressed in the same PC 1 score units as in the graph can be estimated as:

X_co = (-∑ a_i β_i/∑ β_i²) + v₁₁' = (-∑ a_i β_i/∑ β_i²) - | v₁₁'|

[5.6]

where a_i = (N_i1 - m) is the intercept value for the genotype i.

For the lucerne example in Figure 5.4 (A), X_co = -0.19 (calculated over 11 cultivars). Therefore, the coded locations 'p', 'l', 'z', 'a' and 'b' would be assigned to one subregion and the remaining sites to another subregion. While only two groups of locations can be defined by this approach, three or more groups may be defined by the shifted multiplicative model (Crossa et al., 1993, 1995) based on crossover GL effects.

Concluding remarks

AMMI analysis does not provide a direct explanation for the occurrence of GL interaction, since environmental variables are not included in the model. However, values of these variables for sites where they are available can be related to site scores on significant PC axes by correlation or regression techniques. It is thus possible to distinguish the environmental factors likely to play a major role in this context (Annicchiarico and Perenzin, 1994; van Eeuwijk, 1995; Bidinger et al., 1996). Such external information - which may be incorporated into biplots of PC 1 and PC 2 axes (van Eeuwijk, 1995) - can contribute to subregion characterization and scaling-up of results (see Section 5.8). Likewise, morphophysiological traits of genotypes (even if recorded in a small number of sites), averaged across observations for each entry, can be used for correlation analysis with estimated parameters of genotype adaptation (i.e. significant PC scores and mean yields) in order to highlight characters associated with specific or wide adaptation (see Section 6.3). Correlation results for sites and genotypes are expected to be consistent with each other. For example, if the ordination of sites on PC 1 has high positive correlation with rainfall, the ordination of genotypes on the same axis probably reveals positive correlation with lateness of flowering (or crop cycle) and/or negative correlation with physiological indicators of drought resistance, as mechanisms of drought escape and resistance may produce specific adaptation to low rainfall (i.e. low PC 1 score) sites. In joint regression analysis, environmental variables and adaptive traits associated with the occurence of GL interaction may be revealed, respectively, by correlation of site mean yield with environmental variables and by correlation of regression slope with morphophysiological traits of genotypes.

5.4 Factorial regression modelling and complementary analyses

Early models of genotype regression on environmental factors are proposed by Shukla (1972a) for one environmental covariate, and by Hardwick and Wood (1972) for two or more covariates in a multiple regression approach. Applications of these methods, mostly targeted to analysis of GE rather than GL effects, are reported by Wood (1976), Saeed and Francis (1984), Kang and Gorman (1989), Piepho et al. (1998) and others. One development of this technique has involved the simultaneous use of environmental and genotypic covariates, as proposed by Denis (1980 and 1988) and applied by Biarnès-Dumoulin et al. (1996), Balfourier et al. (1997) and others. Regressions are usually linear, although quadratic terms may also be included as additional covariates in the multiple regression model (e.g. Saeed and Francis, 1984).

Factorial regression is considered herein for models including only environmental covariates. The interaction effects GL_ij are modelled as a function of the mean value (V_jn) on the site j of the environmental variable n. If β_in represents the regression coefficient of the genotype _i on the covariate n:

GL_ij = α_i + ∑ β_in V_jn + d_ij

where α_i is the intercept value for the genotype and d_ij is the deviation from the model. For one environmental variable equal to site mean yield, the model coincides with Perkins and Jinks' joint linear regression. Positive β values are indication of a good response to sites with a high level of the environmental variable concerned, while negative values indicate a good response to sites with a low level of variable. Specific adaptation to these sites implies relatively high values of the other adaptation parameter: genotype mean yield. Qualitative covariates (e.g. soil type) may be used, but information must be incorporated through a set of dummy variables (Piepho et al., 1998).

Estimation and test of model parameters

The factorial regression model is constructed with the progressive addition of the most important covariates (Denis, 1988). In the absence of powerful or specific statistical software (see Section 5.9), the best one-covariate model may be identified on the basis of the proportion of GL interaction SS accounted for, i.e. the model R², using procedures similar to those described for joint linear regression. There are two possible procedures:

• the model SS relating to regressions of GL effects for individual genotypes are summed up and expressed as the proportion of the total SS in the regressions; or

• SS values for genotype regressions and deviations from the model (or residual GL interaction) are obtained as the model and error SS, respectively, in an analysis of covariance of all GL effects as a function of the genotype factor and the interaction of genotype with site mean value of the covariate.

Testing for statistical significance of the covariate (i.e. heterogeneity of genotype responses to the environmental variable) can be limited to the best model. For integration with ANOVA results in a factorial regression table, the SS for the covariate and the residual GL interaction need be multiplied by N' = no. test years (or crop cycles) × no. experiment replicates. Alternatively, they can be calculated by multiplying the model R² (for the SS for the covariate) and its complement to one (for the residual GL interaction) by the ANOVA GL interaction SS. The GL interaction DF is partitioned into a (g - 1) portion for the covariate and a remainder for the residual. The MS for both terms can be tested on the appropriate error term for the GL interaction in the ANOVA (Saeed and Francis, 1984). An alternative error for the covariate term can be represented by the variation for the same error pooled with that of the residual GL interaction (Denis, 1988).

After identifying the best one-covariate model and verifying the significance of its covariate term, the best two-covariate model is found by comparing (on the basis of their R²) the possible multiple linear regression models including the best single covariate (previously identified). Model SS and R² and residual GL interaction SS are obtained: by summing the results of multiple regression analyses of GL effects executed for individual genotypes; or through an analysis of covariance of all GL effects including, in addition to the genotype factor, a genotype × covariate site mean value interaction term for each of the two environmental variables. Testing of the second covariate (limited to the best two-covariate model) relates to the portion of GL interaction SS explained by its addition, i.e. the partial regression SS, calculated as the difference in model SS between the two-covariate and the best one-covariate model (multiplied by N' = no. test years [or crop cycles] x no. experiment replicates, for inclusion in the factorial regression table). Alternatively, the partial regression SS may be calculated as the difference in R² between the two-covariate and the one-covariate model, multiplied by the ANOVA GL interaction SS. The DF accounted for by the second (and additional) covariate(s) are always (g - 1). If the MS of the added covariate (tested as for the first covariate) is significant, the best three-covariate model can be searched for, and the significance of the last entered variable verified (and so forth, for models of increasing complexity, following a forward selection strategy). Otherwise, the best one-covariate model is retained (if its covariate is significant).

The estimation and testing of parameters of genotype adaptive response, obtained either from analysis of the single genotypes or from analysis of covariance (in which α_i and β_in values are estimated by parameters for, respectively, the genotype and the interaction of genotype with site mean value of the covariate n), can conveniently be considered only in relation to the selected factorial regression model. Testing of each β_in value for difference to zero in relation to deviations from the model of the genotype i (as obtained by individual analysis) is somewhat more reliable than testing in relation to an average value of deviations (as provided by analysis of covariance).

Factorial regression performed on original yields rather than GL effects (a procedure so far not considered) would inflate the GL interaction component accounting for deviations from the model. This source of variation would include the lack of fit not only of individual genotype regressions but of site mean yield (the latter would be responsible for the increase or decrease of all observed genotype yields relative to regressions, whenever the mean yield response of a site to the environmental variable(s) deviated from a linear response). However, this procedure may be useful for predicting genotype responses from information provided by historical, largely unbalanced data sets. In this case, estimating genotype main effects may prove difficult. Therefore, mean yield and GL interaction effects could be estimated (separately for each genotype) by means of a stepwise multiple regression analysis of original yields as a function of the available site covariates. It is very unlikely, however, that the same environmental variables will be selected for all genotypes (thereby increasing the complexity of modelling).

Modelled yield responses and cultivar recommendations

The expected yield of the genotype _i in the location j is, according to the selected factorial regression model:

R_ij = m + G_i + L_j + α_i + ∑ β_in V_jn

where m is the grand mean and G_i and L_j are the main effects for the genotype i and the location j, respectively. Modelling is reliable in the range of observed covariate values, and estimates of parameters should be kept in the model even when they are not statistically significant. As for AMMI analysis, the expected yield responses can be expressed in terms of nominal yields (N_ij) to simplify their calculation and, possibly, their graphical expression.

N_ij = m_i + α_i + ∑ β_in V_jn

[5.7]

where m_i is the mean value of the genotype i.

In the graphical expression there are problems and opportunities similar to those for AMMI analysis. When genotype adaptive responses can be modelled in one dimension as a function of the only covariate in the model, nominal yields can be represented by straight lines (connecting their values - calculated only for the extreme covariate levels in the graph), as reported in Figure 5.1 (B) for response to site mean rainfall (from January to June) of the same cultivars already modelled by joint regression (Fig. 5.1 [A]). Straight lines are not adequate for modelling actual yields (by addition of site main effects), unless the response to rainfall of site mean yield is perfectly linear. The comparison of joint regression and factorial regression results reveals some inconsistency (as is normal for different models). In particular, 'D 3415' is never top-yielding in the range of site mean rainfall (Fig. 5.1 [B]). Specific genotype recommendation may be envisaged for two subregions:

• relatively low mean rainfall (< 400 mm), where 'Messapia' is the preferable cultivar; and

• relatively high mean rainfall (> 400 mm), where 'Karel' and 'W 4267' are preferable.

Multiple comparison among genotypes at different rainfall values may be performed by means of Dunnett's one-tailed critical difference, as described for joint regression (the error term for the test - when the residual GL interaction is not significant - is the GY interaction within locations in the presence of repetition in time of the trials, and the pooled experiment error in the absence of repetition in time). The critical difference may also be calculated for models with two or more environmental covariates.

The graphical expression of genotype responses for a two-covariate scenario is more complex and requires a bidimensional representation of best-yielding genotypes as a function of site values for the covariates. Figure 8.3 in Section 8.3 provides an example of outputted information with the expected pair of top-yielding genotypes for each combination of 10 rainfall by 10 winter temperature levels of sites. Combinations (represented as cells) with the same winning genotypes can be grouped together. The graph requires the preliminary estimation of nominal yields of the genotypes (usually a subset of better-yielding entries, to avoid extensive calculations) for each pair of covariate values according to equation [5.7]. A denser grid of points in Figure 8.3 (e.g. formed by the combination of 20 rainfall by 20 temperature values) allows for a more fine-tuned definition of winning genotypes. The spatial and temporal scaling-up of results is straightforward by the factorial regression approach, depending on the site mean value of the significant covariates. Subregions are defined, grouping locations with same best-yielding material.

Piepho et al. (1998) propose an alternative procedure for genotype comparison at different covariate levels, considering possible differences between components of the error term as they may occur for specific pairs of compared genotypes. Although more appropriate, this approach requires more complex calculation, as error terms are estimated for each considered pair and different factorial regression models may need to be adopted, depending on the value of the site(s) of specific interest.

Location classification for breeding

For preliminary definition of subregions for breeding, factorial regression analysis may be complemented by hierachical cluster analysis performed on site values of the significant environmental covariates. Analysis may also include non-test sites in the target region for which data on the relevant environmental variables are available; results are thus scaled up to new sites. The temporal scaling-up of results can be obtained by inputting long-term rather than test year data of relevant climatic or biotic variables. Earlier recommendations concerning the clustering method and the dissimilarity measure may also be applied. However, the inclusion of new sites or long-term data makes the truncation criterion based on the level of within-group GL interaction more difficult to apply and less reliable. An average linkage clustering strategy (more useful for detecting true discontinuities among locations, although less useful for minimizing the within-group GL interaction - DeLacy et al. [1996a]) is the most appropriate in this context.

Refinement of the cluster analysis technique that usually requires a modest increase in the calculations is represented by the assignment to each environmental variable of a weight proportional to its importance in the factorial regression model, as suggested (in a slightly different context, i.e. modelling of site mean yield) by Brown et al. (1983). The weight is proportional to the partial regression SS of each environmental variable in the selected factorial regression model. The SS value may be calculated as the difference in SS between the regression model with the considered covariate and the regression model without it (all other significant covariates being included). The assessment requires: no additional multiple regression analysis for a two-covariate selected model; one additional multiple regression analysis (excluding the best single covariate) for a three-covariate selected model; etc. For the application of weights, the standardization of each variable to a standard deviation value proportional to its weight is recommended. For example:

• considering three covariates V₁, V₂ and V₃ with, respectively:

  - partial regression SS of 50, 30 and 20
  - estimated standard deviation of s1, s2 and s3

• and assigning a unity weight to the least important variable:

  - the respective weights are 2.5, 1.5 and 1; and
  - the standardization for a given site j is obtained by the variable transformations:

  V'_j1 = V_j1 2.5/s₁

  V'_j2 = V_j2 1.5/s₂

  V'_j3 = V_j3/s₃

For factorial regression models including only one covariate, the estimation of the main crossover point for genotype nominal yields modelled as a function of the covariate value can provide an alternative criterion for the preliminary definition of two subregions for breeding based on the test locations. Extending Singh et al. (1999)'s approach to factorial regression - as was done for AMMI modelling - nominal yields as expressed in Figure 5.1 (B) can be submitted to a double change of origin, whereby the zero value for the abscissa and the ordinate axes coincides with the lowest site mean value of the covariate and the grand mean, respectively. The slope of the genotype _i as a function of the covariate V1 is coincident with that expressed by the b_i1 value in formula [5.7]. The main crossover point expressed in the same covariate units as in the graph can be estimated as:

X_co = ( - Σ a_i β_i / Σ β_i²) + V₁₁

[5.8]

where V₁₁ is the lowest site value for the covariate; a_i = (N_i1 - m) is the intercept value for the genotype _i in the current formulation (which is different from α_i in formula [5.7]); and N_i1 is the nominal yield of the genotype at the value V₁₁ of the covariate.