Additive main effects and multiplicative interaction model

The AMMI model was composed of additive and multiplicative components where Y represents a vector of n independently distributed observations that can be predicted by the categorical variables for genotypes and environments. The additive component, with fixed main effects for genotype (α _i ) and environment (β _j ), was assumed to be a fixed effect, and inferences were restricted to the grey leaf spot response variables, disease incidence and severity (Searle et al. Reference Searle, Casella and McCulloch1992).

A least-squares method was used to estimate these effects via a two-way ANOVA using the means matrix for Y _(gxe). The multiplicative component was estimated by the singular value decomposition (SVD) of the residual matrix from the two-way ANOVA of the genotype–environment means Y _(gxe). This matrix will be denoted as R _(gxe). Generally, the SVD of a matrix A is defined as the product of an orthogonal matrix U by a diagonal matrix S and the transpose of the orthogonal matrix V; thus, $A_{mn} = U_{mn} S_{mn} V_{nn}^T $ . Required conditions are that U^TU = I and V^TV = I; the columns of U are orthonormal eigenvectors of AA^T ; the columns of V are orthonormal eigenvectors of A^TA; S is a diagonal matrix containing the square roots of eigenvalues from U or V in descending order; and A = VS²V^T (Ientilucci Reference Ientilucci2003).

Therefore, the model equation for the ith genotype in the jth environment in the rth block is (Gauch & Zobel Reference Gauch and Zobel1988):

$$Y_{ijr} = \underbrace {\mu + \alpha _i + \beta _j + \rho _{r(\,j)}} _{additive\; terms} + \underbrace {\sum\limits_{h = 1}^p {\lambda _h \gamma _{ih} \delta _{\,jh}}} _{multiplicative\; terms} + \,\varepsilon _{ijr} $$

where Y _ijr is the phenotypic trait (i.e., the proportion of plants affected by grey leaf spot) of genotype i in environment j for replicate r; μ is the grand mean; α _i is the fixed effect for genotype i, where i = 1, 2, …, g; β _j is the fixed effect for the environment j, where $j = {\rm 1},{\rm 2},..., e;\,\lambda _h $ is the singular value for the interaction principal component (IPC) axis k; γ _ih and δ _jh are the IPC scores (i.e., the left and right singular vectors) for genotype and environment, respectively, for axis k; ρ _r(j) is the effect of the rth block in the jth environment; r is the number of blocks; p is the rank of the R _(gxe) matrix that corresponds to the number of main effects from the interaction (PCI) retained by the residual matrix, $p\,{\rm =} \,min{\rm (}{\rm g}{\rm - 1}{\rm,} \;{\rm e}{\rm - 1)};\ {\rm (}\alpha \beta {\rm )}_{ij} = \sum\nolimits_{h = {\rm 1}}^p {\lambda _h \gamma _{ih} \delta _{jh}} $ is the specific interaction of the ith genotype with the jth environment and ε _ijr is the experimental error that is assumed to be independently and normally distributed with a mean of zero and variance $\sigma ^2, \;\varepsilon _{ijr} \sim N(0,{\rm \;} \sigma ^2 ){\rm.} $ The decomposition of the residual matrix into singular values (SVD) permits the partitioning of the least squares from the elements of the R _(gxe) matrix by reducing the number of axes, or K < p such that the model remains informative, where K is the number of axes or PC retained by the model, but with fewer degrees of freedom. This partition is:

$$\mathop \sum \limits_{h = 1}^p \lambda _h \gamma _{ih} \delta _{\,jh} = \mathop \sum \limits_{h = 1}^K \lambda _h \gamma _{ih} \delta _{\,jh} + \mathop \sum \limits_{h = 1 + K}^p \lambda _h \gamma _{ih} \delta _{\,jh}, $$

where $\sum\nolimits_{h = 1 + K}^p {\lambda _h \gamma _{ih} \delta _{jh} = \varphi _{ijr}} $ quantifies the disturbance and ϕ _ijr is the residual containing all of the multiplicative terms not included in the model.

Therefore, using the least-squares approximation to the R _(gxe) matrix by the first n components of the SVD, the reduced model ηijr is estimated by:

$$\hat Y_{ijr} = \underbrace {\hat \mu + \hat \alpha _i + \hat \beta _j + \hat \rho _{r(\,j)}} _{additive\; terms} + \underbrace {\sum\limits_{h = 1}^K {\lambda _h \gamma _{ih} \delta _{\,jh}}} _{multiplicative\; terms}$$

The PCA permits the components from the interaction to capture the decreasing proportion of the variation that is present in matrix GE, or $\lambda _{\rm 1}^{\rm 2} \ge \lambda _{\rm 2}^{\rm 2} \ge {\cdots} \ge \lambda _K^{\rm 2} $ . A sufficient number of components (K) to represent the target model can be identified using Gollob's test (Table 2) (Gollob Reference Gollob1968).

Table 2.

Complete variance analyses of the means according to Gollob (Reference Gollob 1968 )

Generalized linear models

When a distribution is non-normal, GLMs expand the possibilities for statistical modelling. These models allow fitting of n random variables y _i , where i = 1, 2,…, n, that are independently distributed with mean μ _i and an exponential probability density function. These random variables are associated with the explanatory variables x _j, j = 1, 2, …, p by means of a link function g(μ _i ) designated as the linear predictor (η _i ) that is monotonic and differentiable such that:

$$\eta _{\rm i} = g{\rm (}\mu _i {\rm )} = \sum\limits_{{\rm j} = {\rm 1}}^{\rm p} {x_{{\rm ij}} \psi _{\rm j}} $$

where ψ _j represents the coefficients of the linear predictor.

The maximum-likelihood method is the most useful method for estimating the vector of the unknown parameters of the linear predictor ψ _j . Cordeiro & Demétrio (Reference Cordeiro and Demétrio2008) explained that the robust and fast GLM algorithm rarely fails to converge. However, when this does happen, the fitting procedure must be restarted using the current estimate as the starting value for another model.

The deviance function derived from the likelihood ratio statistic tests the significance of the coefficients of the linear predictor. Therefore, in a sequence of k nested models (which have the same probability distribution and link function, but the linear component M ₀ is a special case of the general linear component M ₁) (Dobson Reference Dobson2002), tests of significance are performed using an analysis of deviance (ANODEV) table. Thus the deviance function from the GLMs is analogous to the residual squared sums from least squares. Standardized Pearson residuals, standardized deviance residuals, and Cook's distance measures were used to diagnoses in the quasi-likelihood models.

Quasi-likelihood models

Quasi-likelihood has been used due to the characteristics of the data and the model. Although the GLM represents a great advance in statistical modelling because it allows the fitting of a large number of models, in some instances the choice of an exponential model is not adequate (McCullagh & Nelder Reference McCullagh and Nelder1989), so Wedderburn (Reference Wedderburn1974) proposed quasi-likelihood estimation. Assuming that Var(μ) is a known function of the mean, and φ is the dispersion parameter, the quasi-likelihood function for every observation is

$$Q_i = Q_i (y;\mu ) = \mathop{\int}\limits_{y_i}^{\mu _i} {\displaystyle{{y_i - t} \over {\phi {\cdot} Var{\rm (}t{\rm )}}}dt} {\rm,} \quad y_i \le t \le \mu _i $$

Inference in quasi-likelihood models is similar to that in GLM because quasi-likelihood estimates maximize Q or solve the following system of equations:

$$\sum\limits_{i = 1}^n {\displaystyle{{{\rm (}y_i - \mu _i {\rm )}} \over {\phi V{\rm (}\mu _i {\rm )}}}} \displaystyle{{\partial \mu _i} \over {\partial \psi _j}} = 0{\rm,} \;\,j = 1{\rm,} \ldots {\rm,} p$$

and

$$\sum\limits_{i = 1}^n {\displaystyle{{{\rm (}y_i - \mu _i {\rm )}x_{ij}} \over {\phi V{\rm (}\mu _i {\rm )}}}} \displaystyle{{\partial \mu _i} \over {\partial \eta _i}} = 0$$

In the first system, $\mu _i = g^{ - 1} {\rm (}\eta _i {\rm )} = g^{ - 1} {\rm (}z_i^T \psi _j {\rm )} = h{\rm (}x_i^T \psi _j {\rm )}$ and this expression are based on the GLM theory. The dispersion φ is estimated using the method of moments on the residual vector ${\rm (}Y - \hat \mu {\rm )}$ :

$$\hat \phi = \displaystyle{1 \over {n - p}}\sum\limits_{i = 1}^n {\displaystyle{{{\rm (}y_i - \hat \mu _i {\rm )}^2} \over {V{\rm (}\hat \mu _i {\rm )}}}} = \displaystyle{{\chi ^2} \over {n - p}}$$

where χ ² is Pearson's generalized chi-square statistic for goodness-of-fit, n is the number of observations, and p is the number of parameters (ψ _j ). The generalized function can be estimated in a similar manner to the deviance function, using the difference between the quasi-likelihood logarithm of the current and the saturated models:

$$\eqalign{D{\rm (}y{\rm ;}\hat \mu {\rm )} = & 2\phi {\rm \{} Q{\rm (}y{\rm ;}y{\rm )} - Q{\rm (}\hat \mu {\rm ;}y{\rm )\}} \cr = & - 2\phi {\rm \{} Q{\rm (}\hat \mu {\rm ;}y{\rm )} - Q{\rm (}y{\rm ;}y{\rm )\}}} $$

Because the contribution from the saturated model is zero, then:

$$\eqalign{D{\rm (}y{\rm ;}\hat \mu {\rm )} & = - 2\phi Q = - 2\phi \sum {\mathop\int\limits_{y_i} ^{\mu _i} {\displaystyle{{y_i - t} \over {\phi V{\rm (}t{\rm )}}}}} dt \cr & = - 2\sum {\mathop\int\limits_{y_i} ^{\mu _i} {\displaystyle{{y_i - t} \over {V{\rm (}t{\rm )}}}}} dt}$$

Thus, the quasi-deviance function does not depend on the dispersion parameter φ. The quasi-deviance function $D{\rm (}y{\rm ;}\hat \mu {\rm )}/\phi $ is compared with the percentiles of the χ ² distribution with (n–p) degrees of freedom, although the null distribution of $\phi ^{ - 1} D{\rm (}y{\rm ;}\hat \mu {\rm )}$ is not usually known (Paula Reference Paula2004).

Generalized additive main effects and multiplicative interaction

The GAMMI theory requires the same basic assumptions as the GLMs. The response variable Y must be independently distributed and have a known exponential family distribution, and p associated explanatory variables X _j, where j = 1, 2, …, p are determined by a link function g(μ _i ) that designates a linear predictor η _i that is monotonic and differentiable such that:

$$\eta _i = g(\mu _i ) = \mathop \sum \limits_{\,j = 1}^p X_{ij} \psi _j $$

These linear predictors have been useful to estimate the mean severity of grey leaf spot. Because the link function was the logit in which η = g(μ) = log(μ/1 − μ), the mean proportion of disease was estimated by the relationship

$$g^{ - 1} (\eta _{ijr} ) = \displaystyle{{exp(\mu + \alpha _i + \beta _j + \rho _{r(\,j)} + \sum\limits_{h = 1}^K {\lambda _h \gamma _{ih} \delta _{\,jh} )}} \over {1 + exp(\mu + \alpha _i + \beta _j + \rho _{r(\,j)} + \sum\limits_{h = 1}^K {\lambda _h \gamma _{ih} \delta _{\,jh} )}}} $$

where K is the number of axes considered.

The GAMMI model is applied using van Eeuwijk's algorithm adapted from Sumertajaya (Reference Sumertajaya2007) in R software version 3.0·2 (R Development Core Team 2013) using the R package gnm (generalized nonlinear models) (Turner & Firth Reference Turner and Firth2009). This algorithm uses iterative alternating generalized regression of rows and columns to estimate the parameters. The first step in determining the appropriate model is to identify the distribution and handling of the experimental data. An error plot should be used to visualize whether the data have, for example, a Poisson or binomial distribution instead of a normal distribution. The second step is to fit the GAMMI model, in which each regression includes a GLM class that is arrived at iteratively. This algorithm involves convergence in row regression, in column regression, and in alternating regression (Hadi et al. Reference Hadi, Mattjik and Sumertajaya2010). If the model converges, then ANODEV may then be performed. Finally, the data matrix is represented as a biplot. Figure 1 shows the algorithm necessary for applying the GAMMI model.

Fig. 1.

van Eeuwijk's algorithm for modelling GAMMI, adapted from Sumertajaya (Reference Sumertajaya2007). *Analysis of Deviance (ANODEV).

To determine the number of axes or the number of multiplicative terms in a GAMMI model, a generalization of the AMMI method via the tests described below may be used. The F test does not require a special table and is easy to calculate. The statistic used is F = (Dev. restricted/D.F. sv restricted)− $({\rm Dev}{\rm. \; full}/ {\rm D}{\rm. F}{\rm. \; full})/\hat \phi $ , which approximates the F(_{D.F.source of variation}; _D.F.error) distribution. Where: ${\rm Dev}{\rm. :deviance},\,\hat \phi $ is the dispersion parameter from quasi-likelihood estimation, D.F..sv: degrees of freedom from source of variation that is being tested.

The test proposed by Gollob (Reference Gollob1968) allocates (g − 1)(e − 1) − (2k − 1) = g + e − 1 − 2k degrees of freedom to the eigenvalues associated with the kth axis, where k = 1, 2,…, n and n = minimum (g–1, e–1), which corresponds to the difference between the number of parameters to be estimated and the number of factors applied. Thus, the mean deviance is tested against the estimated error.

Stability is the maintenance or predictability of the response variable in various environments (Annicchiarico et al. Reference Annicchiarico, Bellah and Chiari2005; Cruz et al. Reference Cruz, Regazzi and Carneiro2006). For the incidence or severity of disease, a genotype is considered to be stable when its disease severity percentage is low and constant with respect to environmental variation under both specific and broad conditions. Stability is estimated by analysing the magnitude and sign of the biplot scores corresponding to the selected GAMMI model. Genotypes and environments with low (near zero) scores are considered stable, which is expected for genotypes and environments that have a small contribution to the overall interaction (Duarte & Vencovsky Reference Duarte and Vencovsky1999).

The adaptability of a genotype indicates its ability to take advantage of environmental effects to ensure a high level of productivity. Adaptability is predicted as a function of the responses for each combination of genotype and environment in the model selected by GAMMI IPCAk (axis k: axis of interaction PCA). The correlation between cultivars and the environment is based on the angles between vectors determined by coordinates of the interaction (axis 1, axis 2) and the vertex. The cosine of the two vectors indicates the level of correlation between two corresponding variables (Rencher Reference Rencher2002). Therefore, a small angle indicates highly positively correlated variables, perpendicular vectors indicate non-correlated variables, and an angle greater than 90° indicates a negative correlation.