Getting started

In this section, we will use the data in data_ge. For more information, please, see ?data_ge. Other data sets can be used provided that the following columns are in the dataset: environment, genotype, block/replicate and response variable(s).

The first step is to inspect the data with the function inspect().

library(metan)
insp <- inspect(data_ge,
verbose = FALSE,
plot = TRUE)

print_table(insp)

Individual and joint ANOVA

It is suggested to check if genotype-vs-environment interaction is significant before proceeding with the AMMI analysis. A within-environment ANOVA considering a fixed-effect model is computed with the function anova_ind(). For each environment the Mean Squares for block, genotypes and error are shown. Estimated F-value and the probability error are also shown for block and genotype effects. Some measures of experimental precision are calculated, namely, coefficient of variation, $$CV = (\sqrt{MS_{res}}/Mean) \times 100$$; the heritability, $$h2 = (MS_{gen} - MS_{res})/MS_{gen}$$, and the accuracy of selection, $$As = \sqrt{h2}$$.

indiv <- anova_ind(data_ge, ENV, GEN, REP, GY, verbose = FALSE)
# Warning: verbose is deprecated. It will be defunct in the new release.
print_table(indiv$GY$individual)

The joint ANOVA is performed with the function anova_joint().

library(metan)
joint <- anova_joint(data_ge, ENV, GEN, REP, GY, verbose = FALSE)
print_table(joint$GY$anova)

The genotype-vs-environment interaction was highly significant. So we’ll proceed with the AMMI analysis.

The AMMI model

The estimate of the response variable for the ith genotype in the jth environment using The Additive Main Effect and Multiplicative interaction (AMMI) model, is given as follows:

${y_{ij}} = \mu + {\alpha_i} + {\tau_j} + \sum\limits_{k = 1}^p {{\lambda _k}{a_{ik}}} {t_{jk}} + {\rho _{ij}} + {\varepsilon _{ij}}$

where $${\lambda_k}$$ is the singular value for the k-th interaction principal component axis (IPCA); $$a_{ik}$$ is the i-th element of the k-th eigenvector; $$t_{jk}$$ is the jth element of the kth eigenvector. A residual $$\rho _{ij}$$ remains, if not all p IPCA are used, where $$p \le min(g - 1; e - 1)$$.

The AMMI model is fitted with the performs_ammi() function. The first argument is the data, in our example data_ge. The second argument (resp) is the response variable to be analyzed. The function allow a single variable (in this case GY) or a vector of response variables. The arguments (gen, env, and rep) are the name of the columns that contains the levels for genotypes, environments, and replications, respectively. The last argument (verbose) control if the code will run silently.

AMMI_model <- performs_ammi(data_ge,
env = ENV,
gen = GEN,
rep = REP,
resp = GY,
verbose = FALSE)
# New names:
# *  -> ...15

Note that using the arguments in the correct order, the model above may be fitted cleanly with:

AMMI_model <- performs_ammi(data_ge, ENV, GEN, REP, GY, verbose = FALSE)

The AMMI table

The following comand creates the well-known ANOVA table for the AMMI model. Note that since

print_table(AMMI_model$GY$ANOVA)

Nine interaction principal component axis (IPCA) were fitted and four were significant at 5% probability error. Based on this result, the AMMI4 model would be the best model to predict the yielding of the genotypes in the studied environments.

Estimating the response variable based on significant IPCA axes

The response variable of a two-way table (for example, the yield of m genotypes in n environments) may be estimated using the S3 method predict() applyed to an object of class waas. This estimation is based on the number of multiplicative terms declared in the function. If naxis = 1, the AMMI1 (with one multiplicative term) is used for estimating the response variable. If naxis = min(g - 1; e - 1), the AMMIF is fitted. A summary of all possible AMMI models is presented below.

Member of AMMI family Espected response of the i-th genotype in the jth environment
AMMI0 $$\hat{y}_{ij} = \bar{y}_{i.} + \bar{y}_{.j} - \bar{y}_{..}$$
AMMI1 $$\hat{y}_{ij} = \bar{y}_{i.} + \bar{y}_{.j} - \bar{y}_{..} +\lambda_1 a_{i1}t_{j1}$$
AMMI2 $$\hat{y}_{ij} = \bar{y}_{i.} + \bar{y}_{.j} - \bar{y}_{..} +\lambda_1 a_{i1}t_{j1}+\lambda_2 a_{i2}t_{j2}$$
AMMIF $$\hat{y}_{ij} = \bar{y}_{i.} + \bar{y}_{.j} - \bar{y}_{..} +\lambda_1 a_{i1}t_{j1}+\lambda_2 a_{i2}t_{j2}+...+\lambda_p a_{ip}t_{jp}$$

Procedures based on postdictive success, such as Gollobs’s test (Gollob 1968) or predictive success, such as cross-validation (Piepho 1994) should be used to define the number of IPCA used for estimating the response variable in AMMI analysis. This package provides both. The waas() function compute traditional AMMI analysis showing the number of significant axes according to Gollobs’s test. On the other hand, cv_ammif() function provides cross-validation of AMMI-model family, considering a completely randomized design (CRD) or a randomized complete block design (RCBD).

predicted <- predict(AMMI_model, naxis = 4)

AMMI model for unbalanced data

Singular Value Decomposition requires a complete two-way table (i.e., all genotypes in all environments). Sometimes (for several reasons), a complete two-way table cannot be obtained in a multi-environment trial. metan offers an option to impute the missing cells of the two-way table using Expectation-Maximization algorithms. If an incomplete two-way table is identified in performs_ammi() a warning is issued, impute_missing_val() is called internally and the missing value(s) is(are) imputed using a low-rank Singular Value Decomposition approximation estimated by the Expectation-Maximization algorithm. The algorithm will (i) initialize all NA values to the column means; (ii) compute the first axis terms of the SVD of the completed matrix; (iii) replace the previously missing values with their approximations from the SVD; (iv) iterate steps 2 through 3 until convergence or a maximum number of iterations be achieved.

As an example we will run the AMMI model by omiting H2 from E1 in data_ge2.

miss_val <-
data_ge2 %>%
remove_rows(4:6) %>%
droplevels()

mod_miss <-
performs_ammi(miss_val, ENV, GEN, REP, PH)
# ----------------------------------------------
# Convergence information
# ----------------------------------------------
# Number of iterations: 23
# Final RMSE: 6.007683e-11
# Number of axis: 1
# Convergence: TRUE
# ----------------------------------------------
# Warning: Data imputation used to fill the GxE matrix
# New names:
# *  -> ...9
# variable PH
# ---------------------------------------------------------------------------
# AMMI analysis table
# ---------------------------------------------------------------------------
#     Source  Df Sum Sq Mean Sq F value   Pr(>F) Percent Accumul
#        ENV   3  6.952  2.3173 117.463 5.92e-07       .       .
#   REP(ENV)   8  0.158  0.0197   0.863 5.51e-01       .       .
#        GEN  12  2.470  0.2058   9.000 3.03e-11       .       .
#    GEN:ENV  35  5.286  0.1510   6.603 1.39e-13       .       .
#        PC1  14  4.413  0.3152  13.780 0.00e+00    83.4    83.4
#        PC2  12  0.588  0.0490   2.140 2.12e-02    11.1    94.5
#        PC3  10  0.290  0.0290   1.270 2.59e-01     5.5     100
#  Residuals  94  2.150  0.0229      NA       NA       .       .
#      Total 152 17.016  0.1119      NA       NA    <NA>    <NA>
# ---------------------------------------------------------------------------
#
# All variables with significant (p < 0.05) genotype-vs-environment interaction
# Done!
p1 <- plot_scores(mod_miss, type = 2, title = FALSE)

mod_comp <-
data_ge2 %>%
performs_ammi(ENV, GEN, REP, PH)
# New names:
# *  -> ...9
# variable PH
# ---------------------------------------------------------------------------
# AMMI analysis table
# ---------------------------------------------------------------------------
#     Source  Df Sum Sq Mean Sq F value   Pr(>F) Percent Accumul
#        ENV   3  7.719  2.5728 127.913 4.25e-07       .       .
#   REP(ENV)   8  0.161  0.0201   0.897 5.22e-01       .       .
#        GEN  12  1.865  0.1554   6.929 6.89e-09       .       .
#    GEN:ENV  36  5.397  0.1499   6.686 5.01e-14       .       .
#        PC1  14  4.466  0.3190  14.230 0.00e+00    82.8    82.8
#        PC2  12  0.653  0.0545   2.430 8.40e-03    12.1    94.9
#        PC3  10  0.277  0.0277   1.240 2.76e-01     5.1     100
#  Residuals  96  2.153  0.0224      NA       NA       .       .
#      Total 155 17.294  0.1116      NA       NA    <NA>    <NA>
# ---------------------------------------------------------------------------
#
# All variables with significant (p < 0.05) genotype-vs-environment interaction
# Done!
p2 <- plot_scores(mod_comp, type = 2, title = FALSE)

arrange_ggplot(p1, p2, labels = c("Missing data", "Complete data"))

Rendering engine

This vignette was built with pkgdown. All tables were produced with the package DT using the following function.

library(DT) # Used to make the tables
# Function to make HTML tables
print_table <- function(table, rownames = FALSE, ...){
datatable(table, rownames = rownames, extensions = 'Buttons',
options = list(scrollX = TRUE, dom = '<<t>Bp>', buttons = c('copy', 'excel', 'pdf', 'print')), ...) %>%
formatSignif(columns = c(1:ncol(table)), digits = 3)}`

References

Farshadfar, E. 2008. “Incorporation of AMMI stability value and grain yield in a single non-parametric index (GSI) in bread wheat.” Pakistan Journal of Biological Sciences 11 (14): 1791–6. http://www.ncbi.nlm.nih.gov/pubmed/18817218.

Gollob, H. F. 1968. “A statistical model which combines features of factor analytic and analysis of variance techniques.” Psychometrika 33 (1): 73–115. https://doi.org/10.1007/BF02289676.

Olivoto, T., A. D. C Lúcio, J. A. G. Da silva, V. S. Marchioro, V. Q. de Souza, and E. Jost. 2019. “Mean performance and stability in multi-environment trials I: Combining features of AMMI and BLUP techniques.” Agronomy Journal 111 (6): 2949–60. https://doi.org/10.2134/agronj2019.03.0220.

Olivoto, T., A. D. C Lúcio, J. A. G. Da silva, B. G. Sari, and M. I. Diel. 2019. “Mean performance and stability in multi-environment trials II: Selection based on multiple traits.” Agronomy Journal 111 (6): 2961–9. https://doi.org/10.2134/agronj2019.03.0221.

Piepho, H. P. 1994. “Best Linear Unbiased Prediction (Blup) for Regional Yield Trials: A Comparison to Additive Main Effects and Multiplicative Interaction (Ammi) Analysis.” Theor. Appl. Genet. 89 (5): 647–54. https://doi.org/10.1007/BF00222462.

Purchase, J. L., Hesta Hatting, and C. S. van Deventer. 2000. “Genotype × environment interaction of winter wheat ( Triticum aestivum L.) in South Africa: II. Stability analysis of yield performance.” South African Journal of Plant and Soil 17 (3): 101–7. https://doi.org/10.1080/02571862.2000.10634878.

Sneller, C. H., L. Kilgore-Norquest, and D. Dombek. 1997. “Repeatability of Yield Stability Statistics in Soybean.” Crop Science 37 (2): 383–90. https://doi.org/10.2135/cropsci1997.0011183X003700020013x.

Zali, H., E. Farshadfar, S. H. Sabaghpour, and R. Karimizadeh. 2012. “Evaluation of genotype × environment interaction in chickpea using measures of stability from AMMI model.” Annals of Biological Research 3 (7): 3126–36. http://eprints.icrisat.ac.in/7173/.