AMMI model

The AMMI model (Gauch, 1988) is widely used to analyse the effect of GEI. This model includes two stages. First, an ANOVA is performed to obtain estimates for the additive main effects of environments and genotypes. Secondly, the residuals from the ANOVA are arranged in a matrix with genotypes in the rows and environments in the columns and a singular value decomposition (SVD) is applied in order to explore patterns related to GEI, still present in the residuals. The result of the first two multiplicative terms of the SVD is often presented in a biplot called GE and represents a two-rank approximation of GEI effects.

The rAMMI() function returns the GE biplot. Data in long format is required by this function, i.e. each row corresponds to one observation and each column to one variable (genotype, environment, repetition (if any) and the observed phenotype). If each genotype has been evaluated more than once at each environment, the phenotypic mean for each combination of genotype and environment is internally calculated and then the model is estimated. Extra variables that will not be used in the analysis may be present in the dataset. Missing values are not allowed (but can be imputated, see below).

The GE biplot for yan.winterwheat dataset is shown in Figure 1 along with the sentence used to obtain it. The first argument is the input dataset, then the names of the columns in which the necessary information to apply the technique is found and also the model to be obtained are indicated. Optionally, the percentage of GEI explained by the GE biplot can be added as a footnote with footnote = T, as well as a tittle with titles = T. In this example, BH93, KE93 and OA93 are the environments that contribute the most to the interaction as their vectors are the longest ones. The cultivars m12 and Kat present similar interaction patterns (their markers are close to each other in the biplot) and they are very different from Ann and Aug, for example. The closeness between the cultivar Dia and the environment BH93 indicates a strong positive association between them, which means that BH93 is a extremely favorable environment for that genotype. As OA93 and Luc markers are opposite, this environment is considerably unfavorable for that genotype. Finally, Cas and Reb are close to the origin, which means that they adapt equally to all environments.

rAMMI(yan.winterwheat, genotype = "gen", environment = "env", 
      response = "yield", type = "AMMI", footnote = F, titles = F)

Figure 1: GE biplot based on yield data of 1993 Ontario winter wheat performance trials. The 71.66% of GE variability is explained by the first two multiplicative terms. Cultivars are shown in lowercase and environments in uppercase.

The AMMI model, in its standard form, assumes that no outliers are present in the data. To overcome the problem of data contamination with outlying observations, Rodrigues et al. (2016) proposed five robust AMMI models, which can be obtained in two stages: (i) fitting a robust regression model with an M-Huber estimator (Huber, 1981) to replace the ANOVA model; and (ii) using a robust SVD or principal components analysis (PCA) procedure to replace the standard SVD. Until now, robust AMMI models were not available in any R package. All robust biplots proposed by Rodrigues et al. (2016) can be obtained using rAMMI(). The argument type can be used to specify the type of model to be fitted ("rAMMI", "hAMMI", "gAMMI", "lAMMI" or "ppAMMI"). Since the sample yan.winterwheat dataset does not present outliers, the conclusions obtained with robust biplots will not differ from those made with the classic biplot (Rodrigues et al., 2016). Thus, no interpretation of the robust biplots is presented in this tutorial.

Site Regression model

The Site Regression model (SREG, also called genotype plus genotype-by-environment model or GGE model) is another powerful tool for the analysis and interpretation of MET data in breeding programs. In this case, an ANOVA is performed to obtain estimates for the additive main effects of environments and a SVD is performed on the residuals matrix in order to explore patterns related to genotype (G) and GEI.

As rAMMI() function, GGEmodel() data needs to be presented in a long format and repetitions or extra variables in the dataset are allowed. All the combinations between genotypes and environments must be present.

GGE1 <- GGEmodel(yan.winterwheat, genotype = "gen", environment = "env", 
                 response = "yield")

The output from GGEmodel() is a list with the following elements:

• model: method for fitting the SREG model: "SREG","CovSREG","hSREG"or "ppSREG".
• coordgenotype: plot coordinates for all genotypes in each component.
• coordenviroment: plot coordinates for all environments in each component.
• eigenvalues: vector of eigenvalues for each component.
• vartotal: overall variance.
• varexpl: percentage of variance explained by each component.
• labelgen: genotype names.
• labelenv: environment names.
• axes: axis labels.
• Data: centered input data.
• SVP: SVP method.

The result of the first two multiplicative terms of the SVD is often presented in a GGE biplot (Yan et al., 2000), which represents a rank-two approximation of the G + GEI effects. Plant breeders have found GGE biplots as an useful tools for the analysis of mega-environment (Yan et al., 2001; Yan and Rajcan, 2002) and genotype and environment evaluation (Bhan et al., 2005; Kang et al., 2006; Yan et al., 2007). The GGE biplot addresses many issues relative to genotype and test environment evaluation. Considering the average performance of each genotype, this plot can be used to evaluate specific and general adaptation. In addition, environments can be visually grouped according to their ability to discriminate among genotypes and their representativeness of other test environments. GGE biplot reveals the which-won-where pattern and allows to recommend specific genotypes for each environment (Yan and Tinker, 2005).

Using the output from GGEmodel(), GGEPlot() builds several GGE biplots views, in which cultivars are shown in lowercase and environments in uppercase. The plot also displays the methods used for centering, scaling and SVD. Optionally, the percentage of G + GEI explained by the two axes can be added as a footnote with footnote = T, as well as a tittle with titles = T.

A basic biplot is produced with the option type="Biplot" (Figure 2). In this example the 78% of G + GE variability is explained by the fist two multiplicative terms. The angles between genotypes markers and environments vectors are considered to understand this plot. Thus, for example, Kat performs below the average in all environments, as it has an angle greater than 90\(°\) with all environments. On the other hand, Fun presents an above-average performance in all locations except OA93 and KE93, as indicated by the acute angles. The length of the environment vectors is a measure of the environment’s ability to discriminate between crops.

GGEPlot(GGE1, type = "Biplot", footnote = F, titles = F)

Figure 2: GGE biplot based on yield data of 1993 Ontario winter wheat performance trials. The scaling method used is symmetrical singular value partitioning (by default). The 78% of G + GE variability is explained by the first two multiplicative terms. Cultivars are shown in lowercase and environments in uppercase.

Breeders usually want to identify the most suitable cultivars for a particular environment of interest, i.e., OA93. To do this with GGE biplots, Yan and Kang (2003) suggest drawing a line that passes through the environment marker and the biplot origin, which may be referred to as the OA93 axis. The performance of the cultivars in this particular environment can be ranked projecting them onto this axis. This can be done by setting type = "Selected Environment" and providing the name of the environment (OA93) in selectedE (Figure 3). Thus, at OA93, the highest-yielding cultivar was Zav, and the lowest-yielding cultivar was Luc. The line that passes through the biplot origin and is perpendicular to the OA93 axis separates genotypes that yielded above and below the mean in this environment.

GGEPlot(GGE1, type = "Selected Environment", selectedE = "OA93", 
        footnote = F, titles = F)

Figure 3: comparison of cultivar performance in a selected environment (OA93). The scaling method used is symmetrical singular value partitioning (by default). The 78% of G + GE variability is explained by the first two multiplicative terms.

Another goal of plant breeders is to determine which is the most suitable environment for a genotype. Yan and Kang (2003) suggest plotting a line that passes through the origin and a cultivar marker, i.e., Kat. To obtain this GGE biplots view the argument type = "Selected Genotype" and selectedG = "Kat" must be indicated (Figure 4). Environments are classified along the genotype axis in the direction indicated by the arrow. The perpendicular axis separates the environments in which the cultivar presented a performance below or above the average. In this example, Kat presented a performance below the average in all the environments studied.

GGEPlot(GGE1, type = "Selected Genotype", selectedG = "Kat", 
        footnote = F, titles = F)

Figure 4: comparison of the performance of cultivar Luc in different environments. The scaling method used is symmetrical singular value partitioning (by default). The 78% of G + GE variability is explained by the first two multiplicative terms.

It is also possible to compare two cultivars, i.e. Kat and Cas, linking them with a line and a segment perpendicular to it. To obtain this GGE biplots view the argument type = "Comparison of Genotype" and the genotypes to be compared selectedG1 = "Kat" and selectedG2 = "Cas" must be indicated (Figure 5). Cas was more yielding than Kat in all environments as they all are in the same side of the perpendicular line as Cas.

GGEPlot(GGE1, type = "Comparison of Genotype", 
        selectedG1 = "Kat", selectedG2 = "Cas", 
        footnote = F, titles = F, axis_expand = 1.5)

Figure 5: comparison of the cultivars Kat and Cas. The scaling method used is symmetrical singular value partitioning (by default). The 78% of G + GE variability is explained by the first two multiplicative terms. Cultivars are shown in lowercase and environments in uppercase.

The polygonal view of the GGE biplots provides an effective way to visualize the which-won-where pattern of MET data (Figure 6). Cultivars in the vertices of the polygon (Fun,Zav, Ena, Kat and Luc) are those with the longest vectors, in their respective directions, which is a measure of the ability to respond to environments. The vertex cultivars are, therefore, among the most responsive cultivars; all other cultivars are less responsive in their respective directions.

The dotted lines are perpendicular to the polygon sides and divide the biplot into mega-environments, each of which has a vertex cultivar, which is the one with the highest yield (phenotype) in all environments found in it. OA93 and KE93 are in the same sector, separated from the rest of the biplot by two perpendicular lines, and Zav is the highest-yielding cultivar in this sector. Fun is the highest-yielding cultivar in its sector, which contains seven environments, namely, EA93, BH93, HW93, ID93, WP93, NN93, and RN93. No environments fell in the sectors with Ena, Kat, and Luc as vertex cultivars. This indicates that these vertex cultivars were not the best in any of the test environments. Moreover, these cultivars were the poorest in some or all of the environments.

GGEPlot(GGE1, type = "Which Won Where/What", footnote = F,
        titles = F, axis_expand = 1.5)

Figure 6: polygon view of the GGE biplot, showing which cultivars presented highest yield in each environment. The scaling method used is symmetrical singular value partitioning (by default). The 78% of G + GE variability is explained by the first two multiplicative terms. Cultivars are shown in lowercase and environments in uppercase.

Selecting cultivars within each mega-environments is an issue among plant breeders. Figure 6 clearly suggests that Zav is the best cultivar for OA93 and KE93, and Fun is the best cultivar for the other locations. However, breeders do not select a single cultivar in each megaenvironment. Instead, they evaluate all cultivars in order to get an idea of their performance (yield and stability).

In the GGE biplot it is also possible to visualize mean yield and stability of genotypes in yield units per se (Figure 7 and 8). The GGE biplot based on genotype-focused scaling, obtained indicating SVP = "row" in GGEmodel(), provides an useful way to visualize both mean performance and stability of the tested genotypes. This is because the unit of both axes for the genotypes is the original unit of the data.

Visualization of the mean and stability of genotypes is achieved by drawing an average environment coordinate (AEC). For example, Figure 7 shows the AEC for the mega-environment composed of he environments BH93, EA93, HW93, ID93, NN93, RN93, WP93. The abscissa represents the G effect, thus, the cultivars are ranked along the AEC abscissa. Cultivar Fun was clearly the highest-yielding cultivar, on average, in this mega-environment, followed by Cas and Har,and Kat was the poorest. The AEC ordinate approximate the GEI associated with each genotype, which is a measure of the variability or instability of the genotype. Rub and Dia are more variable and less stable than other cultivars, by the contrary, Cas, Zav, Reb, Del, Ari, and Kar, were more stable.

data <- yan.winterwheat[yan.winterwheat$env %in% c("BH93", "EA93","HW93", "ID93",
                                                   "NN93", "RN93", "WP93"), ]
data <- droplevels(data)
GGE2 <- GGEmodel(data, genotype = "gen", environment = "env", 
                 response = "yield", SVP = "row")

GGEPlot(GGE2, type = "Mean vs. Stability", footnote = F, titles = F, sizeEnv = 0)

Figure 7: average environment view of the GGE biplot based on genotype-focused scaling, showing mean yield and stability of genotypes.

Figure 8 compares the cultivars to the “ideal” one with the highest yield and absolute stability. This ideal cultivar is represented by a small circle and is used as a reference, as it rarely exists. The distance between cultivars and the ideal one can be used as a measure of convenience. Concentric circles help to visualize these distances. In the example, Fun is the closest one to the ideal crop, and therefore the most desirable one, followed by Cas and Hay, which in turn are followed by Rum, Ham, Rub, Zav, Del and Reb, etc.

GGEPlot(GGE2, type = "Ranking Genotypes", footnote = F, titles = F, sizeEnv = 0)

Figure 8: Classification of genotypes with respect to the ideal genotype. Genotype-focused scaling is used.

Although METs are performed to study cultivars, they are equally useful for the analysis of the environments. This includes several aspects: (i) evaluating whether the target region belongs to one or more megaenvironments; (ii) identifying better test environments; (iii) detecting redundant environments that do not provide additional information on cultivars; and (iv) determining environments that can be used for indirect selection. To obtain GGE biplots for comparing environments the environment-focused scaling should be used as is most informative of interrelationships among them (Figure 9 and 10). This is obtained indicatig SVP = "column" in GGEmodel().

In Figure 9 environments are connected to the origin through vectors, allowing us to understand the interrelationships between them. The coefficient of correlation between two environments it is approximated by the cosine of the angle formed by the respective vectors. In this example the relation between the environments for the mega-environment with BH93, EA93, HW93, ID93, NN93, RN93 and WP93 is considered. The angle between the vectors for the environments NN93 and WP93 is approximately 10\(º\); therefore, they are closely related; while RN93 and BH93 present a weak negative correlation since the angle is slightly greater than 90\(º\). The cosine of the angles does not translate precisely into coefficients of correlation, since the biplot does not explain all the variability in the dataset. However, they are informative enough to understand the interrelationship between test environments.

GGE3 <- GGEmodel(data, genotype = "gen", environment = "env", 
                 response = "yield", SVP = "column")
GGEPlot(GGE3, type = "Relationship Among Environments", footnote = F, titles = F)

Figure 9: Relationship between environments. Environment-focused scaling is used.

Discrimination ability as well as representativeness with respect to the target environment are fundamental measures for an environment. An ideal test environment should be both discriminating and representative. If it does not have the ability to discriminate, it does not provide information on cultivars and is therefore of no use. At the same time, if it is not representative, not only does it lack usefulness but it can also provide biased information on the evaluated cultivars.

To visualize these measurements, an average environment coordinate is defined and the center of a set of concentric circles represents the ideal environment. Figure 10 shows the GGE biplots view for the mega-environment with BH93, EA93, HW93, ID93, NN93, RN93 and WP93. The angle between the vector of an environment and the AEC provides a measure of representativeness. Therefore, EA93 and ID93 are the most representative, while RN93 and BH93 are the least representative of the average environment, when the mega-environment is analyzed. On the other hand, an environment to be discriminative must be close to the ideal environment. HW93 is the closest to ideal environment and therefore the most desirable of the mega-environment, followed by EA93 and ID93. By the contrary, RN93 and BH93 were the least desirable test environments of this mega-environment.

GGEPlot(GGE3, type = "Ranking Environments", footnote = F, titles = F)

Figure 10: classification of environments with respect to the ideal environment. Environment-focused scaling is used.