Multi-environment trials data analysis: An efficient biplot analysis approach


 The analysis of multi-environment trials (MET) data has a long history in plant breeding and agricultural research, with the earliest approaches being based on ANOVA methods. ANOVA-based biplot analysis has been used for a long time in analyzing MET data, and advances have been made employing different modeling approaches. This paper presents MET data analysis using mixed model approaches, and compares three methods of biplot analysis, namely genotype main effects plus genotype by environment interaction (GGE) analysis, factor analytic multiplicative mixed (FAMM) model analysis, and combined model analysis. Ten grain yield datasets from the national variety trial series conducted by the Ethiopian institute of agricultural research were used for this study. Our results revealed that spatial and FA model provide a significant improvement in analyzing MET data. This was demonstrated with evidence of heritability measure. We demonstrated that biplot analysis based on the approached of combined model analysis provides a substantial increase in the total percentage of genotype by environment (G×E) variance explained by the first two multiplicative components for both types of balanced and unbalanced datasets. Thus, by estimating the G×E mean values with the best linear unbiased predictions using spatial+FA (FAMM model analysis), and thereby conducting biplot analysis based on the combined model analysis, plant breeding and trial evaluation programs can have a more robust platform for evaluation of crop cultivars with greater confidence in discriminating superior cultivars across a range of environments.

of genotype by environment (G×E) variance explained by the first two 1 Introduction Plant breeding and evaluation programs conduct a series of field trials, often called multi-environment trials (MET), to evaluate the performance and stability of individual crop varieties across a range of agroecological environments. MET data analysis has a long history with the earlier statistical methods of ANOVA, and has seen advancement in the generation of accurate and informative results using efficient statistical models (Piepho et al, 2021;Yan et al, 2007;Kelly et al, 2009;Smith et al, 2005).
Additive main effects and multiplicative interaction (AMMI) and genotype main effects plus genotype by environment interaction (GGE) analysis are the most popular methods for MET data analysis, as they can provide biplot techniques and measures of stability in addition to ANOVA results (Rodrigues et al, 2016;Yan et al, 2000). The GGE and AMMI analysis have been used for a long time to provide insight into the nature of genotype × environment (G×E) interaction. These methods use the approach of a two-stage data analysis, in which genotype means are first estimated separately for each trial using ANOVA methods, and then combined to form the data for an overall analysis. Yan et al (2007) demonstrated that the GGE biplot tool is superior to AMMI biplots in mega-environment analysis and genotype evaluation because it explains more GGE. Compared to AMMI biplot analysis, the GGE biplot methodology has become an increasingly popular data visualization tool for MET data analysis, as it allows visualizing the which-won-where pattern of the MET data, the interrelationship among test environments, and the ranking of genotypes based on both mean performance and stability (Yan and Tinker, 2006).
The other more efficient methodology for the analysis of MET data is based on mixed model approaches. Recently, mixed model approaches have become predominant as they provide a flexible framework in which incomplete data can be easily handled and non-genetic variances between and within environments (spatial variation within environments and error variance heterogeneity across environments) can be appropriately modeled (Smith et al, 2001a). The error variation within environments can be modeled using the approach of Gilmour et al (1997), as it appropriately models the three patterns of spatial trends associated with the field of trials, namely local, extraneous, and global trends. Smith et al (2001bSmith et al ( , 2005 extended the GGE analysis by employing factor analytic multiplicative mixed (FAMM) models. These authors assume random genotype effects and fixed environment effects and use the approach of one-stage analysis, in which the models for residual effects are estimated simultaneously with models for G×E effects. A key feature of the factor analytic (FA) model for MET data is its importance for the estimation of the associated variance structure for G×E effects, as it provides a good and parsimonious approximation to the unstructured form and is generally more computationally robust (Kelly et al, 2007). Additionally, the best linear unbiased predictions (BLUPs) of the G×E effects and estimations of loadings and scores can be obtained, and more importantly, biplot analysis can be done in order to better understand the G×E interaction patterns.
Some studies have demonstrated the application of the combined model analysis (AMMI or GGE analysis based on the BLUPs estimated from the FAMM model analysis). Sinebo and Taye (2011) conducted GGE analysis based on the BLUPs versus the cell means of G×E two-way table data with genotypes in the rows and environments (location-year combinations) in the columns. The BLUPs were estimated from RCBD-based mixed-effect model analysis treating genotypes as random. As a result, these authors demonstrated that the GGE biplot analysis based on the BLUPs gives more insight into the nature of G×E interaction than the actual cell mean-based analysis. Santos and Marza (2020) also showed that the combined model analysis can accumulate greater variance than the conventional way of analysis in the first two principal components. By using one dataset of six sites/environments, Zhang et al (2020) demonstrated that the GGE biplot analysis based on the raw versus BLUPs data gives different results in mega-environment delineation, test-environment evaluation, and genotype evaluation.
The biplot analysis is an important aspect in plant breeding and evaluation programs in analyzing MET data, as it directly provides breeders with informative results. The inefficiency of ANOVA methods in analyzing incomplete G×E datasets motivate the use of mixed model approaches. The aim of this study, therefore, was to present MET data analysis using mixed model approaches, and to compare the GGE, FAMM model, and combined model analysis for the efficiency of biplot analysis. In this study, the combined model analysis focuses on the GGE analysis based on the BLUPs of G×E mean values estimated from FAMM model analysis Ten MET datasets from the national variety trial (NVT) series conducted by the Ethiopian institute of agricultural research (EIAR) were used in this study. In the FAMM model analysis, the error variation within environments is spatially modeled using the approach of Gilmour et al (1997), and we use FA models to model the variance structure for G×E effects. A measure of heritability is employed as a measure of improvement in the MET data analysis results.
The percentage of G×E variance explained by the first two multiplicative components is used as a measure of efficiency or goodness of fit for the biplot analysis. Data analysis is done using ASReml and GGEBiplots statistical R packages.
The paper is arranged as follows. In section 2 we describe the MET datasets used in this study. ANOVA-based models, mixed models, and the heritability formula are introduced. The data analysis procedure is described. In Section 3 we write the results from MET data analysis using mixed model approaches, biplot analysis , and the discussion follows in Section 4. Finally, we conclude our study in Section 5.

Data
Ten MET datasets from the national variety trial (NVT) series conducted by the EIAR were used in this study. All trials were carried out by randomized complete block design (RCBD) laid out in a rectangular (row × column) array of plots, and sown in different parts of the country by the research programs of common bean, chickpea, wheat, and maize from 2014 to 2020. The entries have a minimum of two replicates in each trial. In this study, trial and environment are used interchangeably, where environment/trial is year by site combination. Similarly, we use entry, genotype, and variety interchangeably. The trait which we use for illustration is the harvested grain yield, expressed in tonnes per hectare. Concurrence of entries within and between years was high. The name of each dataset was designated using the first letter of the trial series name, two letters for the crop name, and the last two digits of years for the duration of the trial series (Table 1). In this study, both types of datasets, complete and incomplete, were considered. Thus, the datasets, namely LCB14-16, ACB14-16, SCB14-16, and BWT19-20 are complete (all entries grown in all trials), and the concurrence of entries just equals the total number of entries. The total number of environments for each dataset ranges from eight to sixteen.

ANOVA based models
The base-line statistical model for MET data analysis can be written as where Y ijk is grain yield of the k th replicate of entry i in environment j (i = 1, 2, . . . , m, j = 1, 2, . . . , s, k = 1, 2, . . . , n), η ij is the empirical/least-square mean effect of entry i in environment j, µ is an overall mean effect, α i is the main effect for genotype i, β j is the main effect for environment j, γ ij is the interaction effect for genotype i in environment j, and ϵ ijk is the random error effect for replicate k of genotype i in environment j, assumed to be NID(0, σ 2 ). The methods follow the approaches of two stage data analysis, in which the two-way table means η ij are estimated first using ANOVA methods, and then the G×E analysis using GGE or AMMI model. The models for the second stage analysis can be written as where l = 1, 2, . . . , c, λ l is the singular value of the l th multiplicative/principal component (PC), with c ≤ min(m − 1, s), τ il is the eigenvector of genotype i for PC l, θ jl is the eigenvector of environment j for PC l, and δ ij is the residual associated with genotype i in environment j, assumed to be NID(0, σ 2 /r) where r is the number of replications within an environment. The models are subject to the constraints λ 1 ≥ λ 2 , ..., λ c ≥ 0 and orthogonality constraints on the τ il scores, that is similar constraints on the θ ij scores by replacing symbols (i, m, τ ) with (j, s, θ). AMMI analysis uses the model in equation 2 whereas GGE analysis uses the model in equation 3

Mixed models
A general form of mixed model for the n × 1 vector y of individual plot yields combined across trials can be written as where η is the ms × 1 vector of G×E effects, µ is the overall mean, the vectors 1 ms and 1 s denote the unit vectors of length ms and s, respectively. The operator, ⊗, denotes the Kronecker product. The matrix I m denotes the t × t identity matrix, β e is the s × 1 vector of fixed environment main effects, and γ g is an ms × 1 vector of random genotype effects for individual environments (ordered as genotype within environment). The n × ms matrix, M, is an indicator matrix that may contain columns whose elements are all zero if the corresponding genotype by environment combination is not present in the data. The vectors β p and γ p are vectors of fixed and random effects, respectively, with associated design matrices X p (assumed to have full column rank) and Z p . These vectors represent environment-specific effects that are peripheral to the genotype effects for each trial. β p and γ p together, therefore, include all effects associated with the block structure of each environment and spatial terms within environments for extraneous and global variation. And, the vector e is the n × 1 vector of plot error effects combined across trials. Some statistical assumptions are made about the random terms of the general mixed models. The random effects in equations 4 are, therefore, assumed to follow a Gaussian distribution with zero mean and variance matrix, namely Smith et al (2001bSmith et al ( , 2005 extended the model for η by employing FAMM models. These authors assume random genotype effects and fixed environment effects and use the approach of one-stage analysis in which the models for residual effects are estimated simultaneously with models for the G×E effects. The Smith model, therefore, which is more of an analogue of the GGE model, can be written as where Λ e is the s×c matrix of environment loadings, f v is the associated mc×1 vector of genotype scores, c is the number of components (multiplicative terms) included in the analysis, and δ is the ms × 1 vector of residual for η effects. It is assumed that β e comprises vector of fixed environment main effects, and f v and δ are random effects with where ψ is a diagonal s×s matrix of specific variances. Here, we do have FAMM models from Smith et al (2001b) modeling approach that includes models in equation 4 and 5 for ϵ and η, respectively. The combined model approach follows the model specification in equation 3, where the value for η ij is the BLUP estimated from FAMM model analysis.
Our data for this approach, therefore, would be empirical BLUPs for the mean value of entry i in environment j.

Heritability formula
Following the approach of Cullis et al (2006), the heritability (H 2 j ) value for the j th trial was calculated from a generalized formula that takes unbalanced data into account, as where A j is the average pairwise prediction error variance of genotype effects for the j th environment, and σ 2 gj is the genetic variance at environment j.

Data analysis procedure
In this study, data analysis was conducted using mixed model approaches of spatial techniques for plot error effects and FA models for G×E effects. In the first step, a separate analysis for each trial was conducted to account for non-genetic effects through spatial models using the approach of Gilmour et al (1997). The second step analysis, modeling the G×E effects, was then implemented using model fitting procedures demonstrated by De Faveri (2013) and Smith (1999). For G×E analysis, we fitted a diag model which is a combined form of individual trial models constructed in step one analysis using the diag function built into the ASReml package. The model forms the basis of a sequence of models to be fitted for the G×E analysis, and it helps to organize the trial-specific models in combined form and to confirm the presence of genetic variance in each trial. If any trial is found to have no genetic variance then it would be excluded from the MET data analysis. FA models were then considered through maintaining the spatial models, as specified in the diag mode. This was simply done by replacing the dig function with the fa function of the ASReml. The adequacy of the FA models of several k orders was formally tested, as it is fitted within a mixed model framework. A model with k factors denoted FA-k, is nested within a model with k + 1 factor. A model with FA-1 was compared with a model with FA-2, a model with FA-2 was compared with a model with FA-3, and so on, and a plausible model was retained from the comparison. Residual maximum likelihood ratio test (REMLRT) was used to compare such models. Where REMLRT was not plausible, we used the total percentage of the G×E variance explained by the factor components of the model to compare these models. The GGE biplot analysis was done based on the G×E observed means table. The G×E observed means table was constructed for only entries grown in all trials for the reason that the GGE biplot analysis works only for a complete and/ or balanced dataset. This was done with the help of the GGEModel function built into the GGEBiplots R package. The biplot analysis based on the FAMM model analysis, conducted for both types of balanced and unbalanced datasets, was done using the summary.fa function of ASReml. It computes the loadings, factor scores, and the total percentage of G×E variance accounted for by the FA regression. The third biplot analysis was just conducted based on the combined model analysis. The two-way G×E means table for this type of biplot analysis was BLUPs estimated from FAMM model analysis.

Analysis using mixed model approaches
The results of the FA model comparisons for each dataset are presented in Table 2, which includes residual log-likelihoods (LR), residual maximum likelihood ratio tests (REMLRT), and the total percentage of the G×E variance explained by the factor components of the model (%var). The datasets LCB14-16, SCB14-16, LCB19-20, DCP14-16, and DCP16-18 were fitted with a model FA-2. Except for the dataset LCB14-16, an FA-2 model could explain more than 75% of the G×E variance with the first two factor components, implying that the precision of the biplot analysis results may be trusted. Despite the fact that the REMLRT indicates that an FA-2 is the final model for the dataset LCB14-16, the percentage of G×E variance explained by the two factor components was less than 50%, which is a modest percentage. An FA-3 model was used to fit the datasets ACB14-16, SCB19-20, IHMZ20, and BWT19-20, although the dataset ACB14-16 was not successfully modeled, with the twofactor component accounting for around 28% of the G×E variation. We also found an FA-4 as a final model for the dataset BCB15-18, with the two factor components explaining around 60% of the G×E variance. We found low parentage of pairs of environments with absolute correlation greater than or equal to 0.5 (see Table reftab3) for datasets ACB14-16, LCB14-16, and BCB15-18, where G×E variance was not well explained by the two factor components, or/and REMLRT could not show significant improvement for each fitted factor component. Figure 1 presents the comparisons of the mixed model approaches, namely non-spatial, spatial, and spatial+FA (FAMM model analysis) in analyzing each dataset using average heritability (H 2 ). Non-spatial analysis conducted as the combined ANOVA with random genotype and replication factors, whereas spatial analysis comprises additional factors for variations associated with the trial field. Spatial+FA extended the analysis for modelling the G×E interaction effects with FA models incorporating spatial information. The spatial analysis improves average heritability, and more improvement is resulted from Spatial+FA. This suggests that by detecting non-genetic variation associated with the agricultural field experiments using spatial models and parsimoniously leveraging the information stored in the G×E interaction using FA models, we could enhance precision in MET data analysis.  combined model analysis results in a substantial improvement in the %varPC2 for balanced and unbalanced datasets. We noted that the biplot analysis based on the BLUPs estimated from FAMM model analysis is more efficient than the biplot analysis based on the GGE and FAMM model approaches, and this suggests that the interpretations in mega-environment delineation, testenvironment evaluation, and genotype evaluation would be plausible on the biplots. In the GGE biplot analysis, the %varPC2 is relatively large for balanced datasets compared to that of obtained using the FAMM model analysis. This is also noted for the unbalanced dataset LCB19-20, but it has a high parentage of pairs of environments with an absolute correlation greater than or equal to 0.5 compared to the other datasets (Table 3). Figure 2 also shows Multi-environment trials data analysis   Fig. 2 The %varPC2 (total percentage of G×E variance explained by the first two multiplicative components) using three methods of biplot analysis: GGE, FAMM, and combined model analysis

Discussion
The MET data analysis being based on the ANOVA methods has a long history and has been used for the evaluation of crop varieties across a range of growing environments. The GGE and AMMI analysis are the commonly used and popular ANOVA-based methods for the MET data analysis in crop breeding research programs. The GGE biplot methodology has become increasingly popular data visualization tool in analyzing MET data, as it allows visualizing the which-won-where pattern of the MET data, the interrelationship among test environments, and the ranking of genotypes based on both mean performance and stability (Zhang et al, 2020). The ANOVA-based model approaches, on the other hand, have been criticized with three frequently mentioned drawbacks. The first is that MET data analysis being based on the ANOVA methods requires balanced /complete genotype by environment means table to conduct the analysis. It is difficult to get a complete dataset for all trials, and data can also be lost during data collection, and even a variety may not be grown in all trials due to seed limitation. Second, the efficiency of ANOVA methods would be unsatisfactory when the field of trials is found to be characterized by large plot variability. The third one is that it is a very rare case to get no significant statistical tests for the distributional assumptions surrounding the residual error.
Mixed models are suitable in analyzing MET data (including incomplete MET data), and to go for extended analysis through relaxing the distributional assumptions surrounding the residual error. This was demonstrated by modeling spatial variation, popularized by the approach of Gilmour et al (1997), focusing on three scales; global, extraneous, and local trends, and by further modeling G×E interaction using FA models, popularized by Smith et al (2001b). As shown in Figure 1, the spatial analysis improves estimations for genetic parameters when compared to non-spatial analysis. Modeling G×E interaction with FA models in combination with models for spatial variations resultes in a substantial improvement in MET data analysis results. The FA models are found to be useful not only for adequately estimating/predicting G×E interaction effects but also for estimating the variance of G×E and conducting biplot analysis. Correlated environments can be established based on the estimated G×E variance and breeders can choose genotypes based on BLUPs averaged across correlated environments.
In this study, the MET data analysis under a linear mixed model was conducted based on the one-stage approach, in which the models for residual effects are estimated simultaneously with models for G×E effects, and this is more efficient and outperforms the two-stage analysis like in GGE or AMMI analysis. The results in Table 3 show that the G×E effects can be modeled using a few factor components for datasets of correlated environments. For datasets with a strong G×E pattern, we may be forced to fit many factor components. The variance component analysis for important genetic merits, such as genetic variance and heritability was adequately estimated from the FAMM model analysis.
The biplot analysis is a popular and very useful technique for generating an important summary of results exploited from the analysis of MET data. In this study, three methods, namely GGE, FAMM model, and combined model analysis were compared for the efficiency biplot analysis. The biplot analysis based on GGE analysis is popular and has a well-established computing program, such as the GGEBiplots R package developed by Yan and Kang (2002) and is routinely used in most plant breeding programs. Yan et al (2007) revealed that the GGE biplot tool is superior to AMMI biplots in megaenvironment analysis and genotype evaluation through more GGE. Thus we abandoned AMMI biplot analysis and focused on GGE biplot analysis. The analysis, however, necessitates a complete G×E means table estimated using the methods of ANOVA, which would be impossible to achieve with an unbalanced dataset. The problem with GGE biplot analysis is that the G×E means table may not be adequately estimated because the spatial variability associated with the field cannot be accounted for and the assumptions surrounding the residual error may be unfeasible. As compared to the biplot analysis based on the FAMM model analysis, the %varPC2 was large for balanced datasets (Figure 2). Unlike the GGE analysis, FAMM model analysis is not computationally simple, and we may not have a platform like the GGEBiplots R package to directly submit the MET datasets for biplot analysis. However, the G×E effects and it's variance could be adequately estimated by modeling nongenetic variation and the G×E interaction for both balanced and unbalanced datasets using the ASReml facilities. As a result, the precision for estimating G×E effects and the type of dataset may not be an issue for the biplot analysis using the FAMM model.
Most of the studies follow the two-stage approach for the biplot analysis by estimating the genotype by environment means table from fixed effect analysis of the individual trials. Peixouto et al (2016) and Sinebo and Taye (2011) conducted the GGE biplot analysis through fixing genotype as random for both separate and combined analysis. The genotype by environment means table were BLUPs estimated from the individual trial data analysis. Smith et al (2001a) criticized this approach, claiming that it results in a doubleshrinkage in the estimation of G×E effects and that it must be unshrunk back to fixed-effect estimates.
Our biplot analysis based on the combined model approach differs in that we used fixed-effect analysis for the G×E analysis despite the fact that the G×E means value were BLUPs estimated from fitting FAMM model, which addresses the problem of double-shrinkage in the estimation of G×E effects. Zhang et al (2020) demonstrated that the biplot analysis based on the combined model approach is more reasonable due to the stronger analytical ability of the first two multiplicative components as compared to the biplot analysis based on GGE analysis. But they didn't show the efficiency of such analysis as compared to the biplot analysis directly computed from the FAMM model analysis. Gogel et al (2018) considered the approach of one-stage analysis to be the gold standard analysis of MET data, and they demonstrated that it outperforms the two-stage analysis. Two-stage analysis, an approximation of one-stage analysis, is a more convenient approach particularly when we are faced with a computational facility and an incomplete data problem. Yan and Tinker (2006) well documented the GGE biplot analysis, which is based on two-stage analysis, and they made it easy for use by developing a graphical computational tool. Our findings here, GGE analysis based on the combined model approach considered the advantages of the approach of one-stage analysis (uses BLUPs estimated from fitting FAMM model) and the Yan and Kang (2002) GGE biplot methodology, and this would bring efficiency in the biplot analysis.

Conclusions
The MET data cannot be always balanced and/or complete, and the ANOVAbased approaches may not handle its analysis. The linear mixed model provides a powerful framework in which unbalanced and/or incomplete data can be easily handled. The non-genetic and G×E effects can be modeled using mixed model approaches that enhance the precision in the analysis results.
The spatial+FA (FAMM model analysis) provides a significant improvement in the MET data analysis results, and this is demonstrated with the evidence of heritability measure. The improvement in the analysis is due to the fact that each individual trial is spatially analyzed, and the G×E effects are modeled using FA models. The G×E patterns can be explained with a few factor components for datasets that have a high proportion of correlated environments.
Biplot analysis has evolved into an important method to conduct MET data analysis. The biplot analysis based on combined model analysis provides substantial increases in the total percentages of the G×E variance explained by the first two multiplicative components. The combined model analysis is conducted using GGE biplot methodology based on the BLUPs data estimated from FAMM model analysis.
By estimating the GÖE mean values with the BLUPs estimated from FAMM model analysis, and then conducting the biplot analysis based on the combined model analysis, plant breeding and trial evaluation programs can have a more robust platform for evaluation of crop cultivars with greater confidence in discriminating superior cultivars across a range of environments.