IntroductionTop

The main objectives of multi-environmental trials (METs) are assessing cultivar adaptation patterns under different environmental conditions (locations) and investigate genotype – environment (G×E) interactions. The common approach used to analyze multi-environmental trials is the classical analysis of variance (ANOVA). However, the classical ANOVA models are restricted by assumptions regarding variances of experimental effects and experimental error (Crossa et al., 2010; Hu et al., 2013). ANOVA models assume homogeneity of variances and do not include covariances for random effects, which is often unrealistic in case of multi-environmental trials (So & Edwards, 2009; Hu & Spilke, 2011).

Lately, linear mixed models (LMMs) have become recognized and widely used for analyzing METs data (Smith et al., 2005; Yang, 2010). These models are especially useful in the analysis of incomplete datasets. In LMMs, an asset is a chance to use more complex variance-covariance structures for describing the G×E interaction effects. The most realistic one is unstructured (UN) variance-covariance matrix (Meyer, 2009). This variance-covariance matrix is the most liberal. Each parameter (variance or covariance) in the matrix is different and is estimated uniquely from the data (Littell et al., 2006; Gilmour et al., 2009; Hu & Spilke, 2011). The disadvantage of the UN structure in comparison to other structures is the large number of parameters to estimate. This can contribute to complexity of numerical calculations. On the other end of the spectrum, we have the compound symmetry (CS) structure, which is the most restrictive variance-covariance matrix for LMM. This structure assumes equal variances and equal covariances.

Recently, researchers more and more often suggest using factor-analytic (FA) structure for the METs data (Piepho, 1997, 1998; Kelly et al., 2007; Meyer, 2009; Stefanova & Buirchell, 2010; Burgueño et al., 2011). That variance-covariance matrix also offers the modeler a flexibility comparable to the UN structure but with a lower number of parameters to estimate. The FA model can be regarded as the mixed model equivalent of the additive main effects and multiplicative interaction (AMMI) model with similar fixed-effects as genotype main effects and G×E interaction effects (GGE). In contrast to AMMI model, which is based on principal components analysis via singular value decomposition, the factor-analytic model is based on factor analysis with a Cholesky factorization (Smith et al., 2001; Burgueño et al., 2011). The factor-analytic model depends on the decomposition of an unstructured variance-covariance matrix.

In case of the analysis of the METs data with LMMs, it is always recommended to carry out a comparison of the models with different variance-covariance structures first (So & Edwards, 2009; Burgueño et al., 2011). The second step is to use the selected model for proper analysis and evaluation of the cultivars. Comparison of the models involves assessment of the goodness of fit and the prediction accuracy. Several studies have evaluated the models with different variance-covariance structures (including FA). Piepho (1998) used and evaluated LMMs with FA structure starting the analysis with one component and ending with seven. He recommended using the two-component FA structure. The model was later used by Kelly et al. (2007) and Burgueño et al. (2011). Hu & Spilke (2011) compared LMMs for five variance-covariance structures, including the FA structure with one component, and concluded that the FA variance-covariance structure is more accurate than the classical ANOVA and numerically less complicated than the UN structure.

The LMMs analysis for METs data are mainly used to assess the yield adaptability of cultivars, through the estimation of the adjusted means conducted as linear combinations of appropriate BLUP (best linear unbiased prediction) and BLUE (best linear unbiased estimation) effects. Adjusted means of yield for genotype and environment combinations obtained from LMM are then used to analyze the G×E interaction effects by AMMI model, GGE plot or pattern analysis. Improperly chosen variance-covariance structures may lead to the biased estimation of adjusted means. Based on the incorrectly estimated adjusted means, G×E interaction analysis may result in incorrect conclusions. So far, other studies were focused on measures of the goodness of fit or the accuracy of prediction for LMMs with different variance-covariance structures. However, it may also be important to assess the G×E interactions in models with different variance-covariance structures. Therefore, in this work we focused on adaptive response of cultivars on the environments modeled by the LMMs with different variance-covariance structures. We also tried to identify possible limitations of inference when using an inadequate variance-covariance structure.

Material and methodsTop

Linear mixed models for METs

We analysed raw data for grain yield using a two-stage approach (Smith et al., 2005; Möhring & Piepho, 2009; Welham et al., 2010; Piepho et al., 2012). In the first stage, raw data were analysed for each trial (each location-year combination) separately using the ANOVA mixed-model for a split-block design. In the analysis, we assumed cultivars to be fixed effects and blocks to be random (Spilke et al., 2005). The least squares (LS) means for cultivars were estimated. Then, the LS means were combined across trials and years to obtain an unbalanced G×L×Y data table. In the second stage, we performed the combined analysis of the means in the unbalanced G×L×Y table based on the linear mixed model:

where X_ijk is the LS mean of yield for the 3-factorial combination of the k-th year, the i-th (i=1,2, …, I) location, the j-th cultivar; m is the general mean; L_i is the fixed main effect of the i-th location; G(L)_ij is the random effect of the j-th cultivar within the i-th location; Y_k is the random main effect of the k-th year; YL_ik is the random interaction effect of the k-th year and the i-th location; YG(L) _ijk is the random interaction effect of the j-th cultivar within the i-th location and the k-th year; and is the random residual effect of the mean experimental error.

We used models with homogeneous and heterogeneous variance-covariance structures for cultivar within environment G(L) effects. The random effects of G(L)_ij were modelled with CS, UN and FA structures. The FA structure for the presented linear mixed model was fitted with different number of components, from one component, FA1, to seventeen components, FA17. The linear mixed models with different variance-covariance structures were fitted using ASReml-R (Gilmour et al., 2009).

The popular method used to choose the best variance-covariance structure is based on information criteria, such as Akaike’s information criterion (AIC) or Bayesian information criterion (BIC) (So & Edwards, 2009; Hu & Spilke, 2011). It depends on evaluating the predictive accuracy of the effects (Piepho, 1998; Kelly et al., 2007). In our study, the classical evaluation of the variance-covariance structures was based on the AIC. This is a commonly used criterion to choose the best variance-covariance structure in linear mixed model for METs. The AIC was obtained using the following formula:

where LogL is the logarithm of maximum restricted likelihood of the model; and p is the number of estimated variance-covariance parameters in the model. The smaller AIC values, the better fit of the variance-covariance structure. The expression “–2 × LogL” is called “deviance” and is also used in calculating other information criteria (e.g. BIC) and likelihood ratio tests.

In the cultivar selection and recommendation process, the cultivar ranking (Roostaei et al., 2014) and its compatibility across locations (Hu, 2015), are very important. This is why in this research we analyzed the Spearman rank correlation. We calculated the correlation coefficients for cultivars effects G and cultivars within the locations G(L) effects estimators from the models with difference variance-covariance structures. Inconsistency of cultivars ranking between the models informs us which model leads to an inaccurate evaluation of cultivars adaptability patterns.

For the description of the cultivars adaptation to the agricultural environments, we calculated adjusted means for cultivar and location G(L) combination in all compared models with different variance-covariance structures. The adjusted means were calculated as combinations of appropriate BLUPs or/and BLUEs effects, using an algorithm described by Welham et al. (2004). The modeled adjusted means for all considered cultivars across 18 locations using different variance-covariance structures are presented in Figs. 1-2. Additionally, in the results we indicated two winter wheat cultivars characterized by wide environmental adaptation ('Mulan' and 'Smaragd'), bred in Germany.

Figure 1. The winter wheat cultivars adaptive yield response patterns across eighteen environments modeled with compound symmetry (a) and unstructured (b) variance-covariance structure. Environmental means, dashed line; 'Mulan' cultivar, blue line; 'Smaragd' cultivar, red line.

Figure 2. The winter wheat cultivars adaptive yield response patterns across eighteen environments modeled with FA variance-covariance structure: a) FA1, b) FA2, c) FA3, d) FA4. Environmental means, dashed line; 'Mulan' cultivar, blue line; 'Smaragd' cultivar, red line.

ResultsTop

The lowest AIC value was observed for the model with FA2 variance-covariance structure for G(L) effects (Table 1). The FA structure with three components was the second best. The next lowest levels of AIC were observed in FA structures with one, four and five components. In contrast, the highest value of AIC was found for an UN variance-covariance matrix, indicating interior fit. The model with this variance-covariance structure has the lowest value of the logarithm of the restricted likelihood ratio LogL. For the model with compound symmetry structure we obtained the highest value of LogL. The LogL value is not in constant decrease with increasing number of components for models with FA variance-covariance structure.

Table 1. Akaike information criterion and restricted log-likelihood values for models with different variance-covariance structures fitted to winter wheat dataset.

We observed significant and positive correlation for cultivar effects between all variance-covariance structures (Table 2). The results with the UN structure were highly correlated with ones obtained using the FA4 structure. They showed that the rankings of cultivar effects were compatible. The correlation calculated among G and G(L) effects using the FA structure with smaller number of components and the UN structure were weak. The results acquired for the CS structure were moderately correlated with other variance-covariance structures. This indicates that the cultivar rankings were not the same. In case of cultivars within location effects we observed similar patterns. However, the values of Spearman correlation coefficients were slightly lower than for cultivars effects. The use of CS model or FA with smaller number of components models can contribute to erroneous selection of superior cultivars.

Table 2. Spearman rank correlations matrix for main cultivars effects (lower triangular matrix) and cultivar within the locations effects (upper triangular matrix) estimates using difference variance-covariance structures.

In the Fig. 1a, we present cultivars yield response across 18 agricultural environments (locations) modeled with the CS variance-covariance structure. The yields of cvs. 'Mulan' and 'Smaragd' modeled using this structure were on the first place in all tested locations. The modeled locations range (max-min) of cultivars’ grain yield stayed on a similar level, amounting to about 0.80 t/ha, in all tested locations. The minimum yield was 0.63 t/ha for location L6 and maximum 0.98 t/ha for L18 (Table 3).

Table 3. Minimal, maximum and range values for the cultivars adjusted yield means (t/ha) for all tested environments modeled by differences variance-covariance structures.

The cultivar yield response across environments modeled by linear mixed model with unstructured variance-covariance matrix is presented in the Fig. 1b. Contrary to cultivars response on environments modeled by the CS structure, the UN structure diversifies the ranges of yield values in tested environments. For this variance-covariance structure we observed that one of the environments had exceptionally greater diversity of cultivars yield than others. The modeled yields of cvs. 'Mulan' and 'Smaragd' were above the mean in all environments but, in most cases, they were not the highest. The modeled ranges of cultivars yields in different locations were from 0.69 t/ha for L6 to 2.44 t/ha or L5. There are environments where the range of grain yield values is three times greater than in other tested environments (Fig. 1b).

Fig. 2 presents reaction of cultivars to the tested environments modeled with an FA structure for the first four components. Particularly for the model with FA1 we observe that some environments are characterized by a very small difference in cultivars yield. The range of yield values is around only 0.20 t/ha, e.g. for location L10 their yield difference was 0.16 t/ha. For this model in environment L5 cultivars yield difference was equal to 2.30 t/ha, minimum yield of 5.50 t/ha, and the maximum one –7.84 t/ha (Table 3). In most locations the cultivars yields on average differed by about 0.40 t/ha. For FA1 structure, 'Mulan' and 'Smaragd' cultivars had yield below environmental mean in some locations. It might suggest that these cultivars were characterized by narrow adaptation.

With increasing number of components for the FA models, the cultivars adaptive yield response patterns across eighteen environments become slightly more similar to the results of the model with UN structure. For FA4, the modeled range of cultivar yields in different locations were from 0.65 t/ha for location L_6 to 2.37 t/ha for L_5. These results were compatible with the UN structure. For the FA model with more than four components the cultivars response to environments did not change much (Table 3, figure not presented).

DiscussionTop

The adaptive yield response of winter wheat cultivars on environments modeled by the CS structure suggests a lack of interaction. In Fig. 1 we observed almost parallel lines of cultivars reaction on different environments, so the yield diversity of the tested cultivars was the same across environments according to the model. The results are far away from the reality. It is very rare that all environments have equal discrimination power for cultivars. In reality there are always differences in cultivar yield reactions to different environments. Thus, evaluation of cultivar adaptability to environments using the CS model is incomplete and may lead to erroneous conclusions. Especially, the ranking for cultivars and cultivars within the locations from the CS structure was not consistent with those from other variance-covariance structures. Other empirical studies also suggest poor performance of the models with CS structure (So & Edwards, 2009; Hu & Spilke, 2011).

A completely different situation occurs when using the UN structure. With the help of this model we can discern environments where cultivars yields were very different from the ones with low discriminatory power. When the profiles for the cultivars reaction to the environment are not parallel, this indicates significant G×E interaction. Using this structure allows us to evaluate both qualitative and quantitative G×E interactions. The UN structure has a number of numerical and computational problems, inter alia predictions may be characterized by lower accuracy. Nevertheless, cultivars yield reaction to environments described by the UN structure more realistically describes G×E interactions in comparison to the CS structure.

Analysis with the FA structure with an optimal number of components gives the same results as the UN structure. But it is free of most of the UN structure problems. The model for PVTS data, which also describes well the reaction of varieties to the environments is the FA4 model. The models with FA4 and UN structures have consistent ranking of cultivars mean and cultivar within the locations effects. However, we have to remember that models with a number of components too low can contribute to an inappropriate description of cultivars response patterns. A larger number of components allows to perform the real description of the cultivars’ response on the environments. The FA models with more than 4 components did not increase the value of Spearman correlations among G or G(L) effects. This means that the ranking of cultivars effects among these models was compatible (similar). To select the optimal number of components in the FA structure we may use the information criterion, e.g. AIC. However, as shown in this paper, this solution is not perfect. In our case, AIC indicated as the best-fitting an FA structure with 2 components, i.e. FA2. Looking at the Fig. 2 we cannot say that the cultivars adaptability is properly characterized by the FA2 model. By contrast, a model with four components FA4, was on the 3rd place according to the AIC criterion. Thus, the use of AIC can give us only a suggestion on how many components for FA structure should be used, requiring further examination. We can pre-select models using the AIC criterion and then plot their results against each other for further interpretation. Piepho (1998) and Kelly et al. (2007) suggested that the optimum number of components is two. However, the datasets used in their studies were lower in the number of levels within each factor than in our dataset.

In order to fully describe the cultivars adaptive patterns we should use the variance-covariance structure for linear mixed models without restrictive assumptions. A good example of that structure is the unstructured variance-covariance matrix. In contrast to the model with compound symmetry structure, it does not limit the inference on cultivars adaptive patterns. The factor-analytic structure is also a good tool to describe the cultivars’ yield reaction to the environments. In our case, for the winter wheat datasets coming from the Polish Post-registration Variety Testing System, the reasonably low optimal number of components is 4 but each dataset requires separate determination of the optimal number of the factor-analytic components.