L'analisi AMMI per le interazioni genotipo-ambiente

Il risultato degli esperimenti agronomici è sempre sotto la diretta influenza delle condizioni pedo-climatiche e, di conseguenza, le prove sperimentali vengono sempre ripetute in più ambienti (anni e/o località). Se, come di frequente capita, l’interazione tra il trattamento sperimentale in studio e l’ambiente è significativa, si viene ad introdurre un elemento di difficoltà per l’interpretazione e la sintesi dei dati, poichè diviene necessario discutere i risultati ambiente per ambiente, utilizzando numerose tabelle e/o grafici dai quali è difficile trarre conclusioni generali.

Esistono alcune tecniche esplorative di analisi multivariata che possono aiutarci a risolvere questo problema, come ad esempio l’analisi AMMI (Addittive Main effect Multiplicative Interaction), che permette di descrivere l’interazione tra fattori sperimentali in modo molto parsimonioso ed efficace, garantendo nel contempo un ottimo livello di leggibilità e sinteticità dei risultati.

Cerchiamo di scoprire questa tecnica con l’aiuto di un caso-studio pubblicato in Stagnari et al., 2007.

Un caso studio

Una prova di confronto varietale con 12 varietà di favino (in realtà si tratta di sei varietà in due epoche di semina, ma, per comodità, trascureremo questo secondo trattamento sperimentale) e quattro repliche è stata ripetuta in due località per tre anni (per un totale sei ambienti). Carichiamo il dataset e vediamo le medie.

rm(list=ls())
dataset <- read.csv("FabaBean.csv", header=T)
head(dataset)
##       Varieta Blocco Localita Prod
## 1    Chiaro_A      1    bad_1 4.36
## 2    Chiaro_P      1    bad_1 2.76
## 3 Collameno_A      1    bad_1 3.01
## 4 Collameno_P      1    bad_1 2.50
## 5    Palomb_A      1    bad_1 3.85
## 6    Palomb_P      1    bad_1 2.21

#Tabella pivot
library(reshape)
GEmedie <- cast(Varieta ~ Localita, data = dataset,
                 value = "Prod", fun=mean)
GEmedie
##        Varieta  bad_1  bad_2  bad_3  pap_1  pap_2  pap_3
## 1     Chiaro_A 4.1050 2.3400 4.1250 4.6325 2.4100 3.8500
## 2     Chiaro_P 2.5075 1.3325 4.2025 3.3225 1.4050 4.3175
## 3  Collameno_A 3.2500 2.1150 4.3825 3.8475 2.2325 4.0700
## 4  Collameno_P 1.9075 0.8475 3.8650 2.5200 0.9850 4.0525
## 5     Palomb_A 3.8400 2.0750 4.2050 5.0525 2.6850 4.6675
## 6     Palomb_P 2.2500 0.9725 3.2575 3.2700 0.8825 4.0125
## 7      Scuro_A 4.3700 2.1050 4.1525 4.8625 2.1275 4.2050
## 8      Scuro_P 3.0500 1.6375 3.9300 3.7200 1.7475 4.5125
## 9    Sicania_A 3.8300 1.9450 4.5050 3.9550 2.2350 4.2350
## 10   Sicania_P 3.2700 0.9900 3.7300 4.0475 0.8225 3.8950
## 11   Vesuvio_A 4.1375 2.0175 4.0275 4.5025 2.2650 4.3225
## 12   Vesuvio_P 2.1225 1.1800 3.5250 3.0950 0.9375 3.6275

#Medie varietali
Gmedie <- with(dataset, aggregate(Prod, list(Varieta), mean))
Gmedie
##        Group.1        x
## 1     Chiaro_A 3.577083
## 2     Chiaro_P 2.847917
## 3  Collameno_A 3.316250
## 4  Collameno_P 2.362917
## 5     Palomb_A 3.754167
## 6     Palomb_P 2.440833
## 7      Scuro_A 3.637083
## 8      Scuro_P 3.099583
## 9    Sicania_A 3.450833
## 10   Sicania_P 2.792500
## 11   Vesuvio_A 3.545417
## 12   Vesuvio_P 2.414583

#Medie ambientali
Emedie <- with(dataset, aggregate(Prod, list(Localita), mean))
Emedie
##   Group.1        x
## 1   bad_1 3.220000
## 2   bad_2 1.629792
## 3   bad_3 3.992292
## 4   pap_1 3.902292
## 5   pap_2 1.727917
## 6   pap_3 4.147292

#Media generale
mu <- mean(dataset$Prod)

A questo dataset può essere adattato un modello lineare del tipo:

\[Y_{ij} = \mu + g_i + e_j + ge_{ij} + \epsilon_{ij}\]

dove le osservazioni sperimentali per ogni i-esima ‘varietà’ e per ognij-esimo ambiente possono essere descritte in modo additivo, sommando la media generale, l’effetto varietà, l’effetto ambiente e l’interazione genotipo x ambiente (GE). Il residuo \(\epsilon\) costituisce una misura dell’errore sperimentale.

Per un disegno bilanciato, gli effetti principali (varietà e ambiente: g ed e) sono dati dalle differenze tra le medie dei genotipi/ambienti e la media generale. Possiamo calcolarli in questo modo:

g <- Gmedie[,2] - mu
e <- Emedie[,2] - mu

g
##  [1]  0.47381 -0.25535  0.21299 -0.74035  0.65090
##  [6] -0.66243  0.53382 -0.00369  0.34757 -0.31076
## [11]  0.44215 -0.68868

e
## [1]  0.116736 -1.47347  0.88903  0.7990 -1.3753  1.04403

L’effetto di interazione (ge) è dato da:

\[ge_{ij} = Y_{ij} - \left( \mu + g_i + e_j \right)\]

Ad esempio, per Chiaro di Torre Lama (A) l’effetto del genotipo è 0.4738, mentre per il primo ambiente l’effetto è 0.1167. La produzione riscontrata per Chiaro di Torre Lama nel primo ambiente è 4.105, quindi l’interazione genotipo-ambiente è data:

$$ ge_{1,1} = 4.1050 - (3.1033 + 0.4738 + 0.1167) = 0.4112 $$

Questo valore positivo implica che quella particolare combinazione ‘varietà x ambiente’ è risultata particolarmente valida, più di quanto ci si sarebbe potuti aspettare guardando separatamente i singoli effetti principali g ed e. Altri esempi aiuteranno meglio a comprendere il concetto di interazione; nel caso di Chiaro di Torre Lama (A) nel terzo ambiente, abbiamo un interazione pari a:

$$ ge_{1,3} = 4.125 - (3.1033 + 0.4738 + 0.8890) = –0.3411 $$

che mostra come Chiaro di Torre Lama (A), pur avendo avuto una buona produzione in quell’ambiente (4.125), non è apparso particolarmente favorito. Viceversa, nello stesso ambiente, Collameno (P) ha mostrato un interazione:

$$ ge_{4,3} = 3.8650 - (3.1033 - 0.7403 + 0.8890) = 0.6130 $$

che dimostra come una varietà mediamente poco produttiva, in questo ambiente abbia trovato condizioni particolarmente favorevoli.

L’equazione sovrastante può essere scritta:

\[ge_{ij} = Y_{ij} - \left( \mu + e_j \right) - g_i = Y_{ij} - \mu_j - g_i\]

il che mostra che possiamo eseguire una centratura della matrice iniziale per le colonne (il che equivale a sottrarre da ogni dato la media dell’ambiente). La matrice risultante può essere ancora centrata per le righe, il che equivale a sottrarre da ogni dati l’effetto del genotipo, dato che, per ogni riga, abbiamo già rimosso la media generale e l’effetto dell’ambiente. Insomma la matrice risultante è doppiamente centrata (“doubly centered”):

GEint <- as.data.frame(t(scale( t(scale(GEmedie, center=T, scale=F)),
                                center=T, scale=F)))
print(GEint, digits=3)
##               bad_1    bad_2   bad_3   pap_1    pap_2   pap_3
## Chiaro_A     0.4112  0.23639 -0.3411  0.2564  0.20826 -0.7711
## Chiaro_P    -0.4572 -0.04194  0.4656 -0.3244 -0.06757  0.4256
## Collameno_A -0.1830  0.27222  0.1772 -0.2678  0.29160 -0.2903
## Collameno_P -0.5722 -0.04194  0.6131 -0.6419 -0.00257  0.6456
## Palomb_A    -0.0309 -0.20569 -0.4382  0.4993  0.30618 -0.1307
## Palomb_P    -0.3076  0.00514 -0.0724  0.0301 -0.18299  0.5276
## Scuro_A      0.6162 -0.05861 -0.3736  0.4264 -0.13424 -0.4761
## Scuro_P     -0.1663  0.01139 -0.0586 -0.1786  0.02326  0.3689
## Sicania_A    0.2624 -0.03236  0.1651 -0.2949  0.15951 -0.2599
## Sicania_P    0.3608 -0.32903  0.0485  0.4560 -0.59465  0.0585
## Vesuvio_A    0.4753 -0.05444 -0.4069  0.1581  0.09493 -0.2669
## Vesuvio_P   -0.4088  0.23889  0.2214 -0.1186 -0.10174  0.1689

L’ANOVA (considerando gli ambienti come un fattore fisso) è:

mod.aov <- aov(Prod ~ Varieta*Localita + Error(Localita:factor(Blocco)), data=dataset)
## Warning in aov(Prod ~ Varieta * Localita +
## Error(Localita:factor(Blocco)), : Error() model is singular
summary(mod.aov)
## 
## Error: Localita:factor(Blocco)
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Localita   5  316.6   63.31   168.6 1.83e-14 ***
## Residuals 18    6.8    0.38                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##                   Df Sum Sq Mean Sq F value Pr(>F)    
## Varieta           11  70.03   6.366  58.411 <2e-16 ***
## Varieta:Localita  55  30.97   0.563   5.167 <2e-16 ***
## Residuals        198  21.58   0.109                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Gli obiettivi e il principio di fondo

La matrice GEint, per la sua dimensionalità, non è facilmente esaminabile e rappresentabile graficamente e crea quindi alcune difficoltà di studio.

L’analisi AMMI (Zobel et al., 1988; Gollob, 1968) può aiutarci a superare queste difficoltà in quanto ci fornisce un metodo per ridurre la dimensionalità della matrice GEint, senza perdere una grossa porzione delle informazioni sperimentali in essa contenute. In pratica, l’analisi AMMI non è altro che una Analisi delle Componenti principali sulla matrice doppiamente centrata ed è quindi basata sulla sua decomposizione ai valori singolari. Ricordiamo che una qualunque matrice X con j colonne e k righe può essere decomposta secondo la seguente espressione:

\[X = U D V^T\]

dove \(U\) è la matrice dei primi n autovettori di \(XX^T\) (di dimensione j x n), \(V\) è la matrice dei primi n} autovettori di \(X^T X\) (di dimensione k x n), \(D\) è la matrice diagonale delle radici quadrate dei primi n autovalori di \(XX^T\) (che coincidono con i primi n autovalori di \(X^T X\)). Il teorema di Eckart-Young (1936; si veda Gollob, 1968) stabilisce che la decomposizione è perfetta se n è uguale al rango di X, mentre la miglior approssimazione si ottiene conservando solo le prime r (con r < n) colonne delle matrici U e V e i primi r termini di D.

I calcoli

In R i calcoli si possono eseguire banalmente con la funzione rda() del package vegan. Usualmente, per l’analisi AMMI si preferisce uno scaling simmetrico, anche se questa non è una regola:

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.4-3
ammi <- rda(GEint, scale=F)
summary(ammi, scaling=3)
## 
## Call:
## rda(X = GEint, scale = F) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          0.7039          1
## Unconstrained  0.7039          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5
## Eigenvalue            0.4809 0.1263 0.05869 0.02382 0.01422
## Proportion Explained  0.6832 0.1794 0.08338 0.03385 0.02019
## Cumulative Proportion 0.6832 0.8626 0.94596 0.97981 1.00000
## 
## Scaling 3 for species and site scores
## * Both sites and species are scaled proportional to eigenvalues
##  on all dimensions
## * General scaling constant of scores:  1.66812 
## 
## 
## Species scores
## 
##            PC1      PC2      PC3     PC4     PC5
## bad_1 -0.83079  0.09477 -0.46711 -0.3173 -0.1507
## bad_2  0.04402 -0.41802  0.07037  0.3706 -0.4025
## bad_3  0.67043 -0.12977 -0.52477  0.1705  0.2977
## pap_1 -0.66137  0.51268  0.28858  0.3143  0.2210
## pap_2 -0.06863 -0.62703  0.41985 -0.2941  0.2079
## pap_3  0.84634  0.56736  0.21307 -0.2440 -0.1733
## 
## 
## Site scores (weighted sums of species scores)
## 
##                  PC1       PC2        PC3      PC4      PC5
## Chiaro_A    -0.60711 -0.383733  0.0008682  0.20801 -0.06349
## Chiaro_P     0.55193  0.026531 -0.0809519  0.04480  0.16410
## Collameno_A  0.08445 -0.542186 -0.0063099  0.17593  0.05688
## Collameno_P  0.80677 -0.065753 -0.1321600 -0.17195  0.07930
## Palomb_A    -0.32131  0.110117  0.5908136 -0.08279  0.38856
## Palomb_P     0.28105  0.345909  0.2816496  0.04235 -0.25302
## Scuro_A     -0.62639  0.139186 -0.1626008  0.01712 -0.08010
## Scuro_P      0.22961  0.076556  0.1817985 -0.20708 -0.24155
## Sicania_A   -0.06287 -0.323857 -0.3547204 -0.27987  0.09016
## Sicania_P   -0.21433  0.683297 -0.4016608  0.14811  0.15052
## Vesuvio_A   -0.43787 -0.007914  0.0202503 -0.29981 -0.17690
## Vesuvio_P    0.31606 -0.058153  0.0630236  0.40517 -0.11447

G <- summary(ammi, scaling=3)$specie
E <- summary(ammi, scaling=3)$sites
print(G%*%t(E), digits=3)
##       Chiaro_A Chiaro_P Colla_A Colla_P Palomb_A Palomb_P Scuro_A
## bad_1    0.411  -0.4572  -0.183  -0.572  -0.0309 -0.30757  0.6162
## bad_2    0.236  -0.0419   0.272  -0.042  -0.2057  0.00514 -0.0586
## bad_3   -0.341   0.4656   0.177   0.613  -0.4382 -0.07236 -0.3736
## pap_1    0.256  -0.3244  -0.268  -0.642   0.4993  0.03014  0.4264
## pap_2    0.208  -0.0676   0.292  -0.003   0.3062 -0.18299 -0.1342
## pap_3   -0.771   0.4256  -0.290   0.646  -0.1307  0.52764 -0.4761
##       Scuro_P Sicania_A Sicania_P Vesuvio_A Vesuvio_P
## bad_1 -0.1663    0.2624    0.3608    0.4753    -0.409
## bad_2  0.0114   -0.0324   -0.3290   -0.0544     0.239
## bad_3 -0.0586    0.1651    0.0485   -0.4069     0.221
## pap_1 -0.1786   -0.2949    0.4560    0.1581    -0.119
## pap_2  0.0233    0.1595   -0.5947    0.0949    -0.102
## pap_3  0.3689   -0.2599    0.0585   -0.2669     0.169

Le matrici G ed E (rispettivamente scores dei genotipi e degli ambienti) moltiplicate tra di loro restituiscono la matrice GEint. Ci si può a questo punto chiedere: si, ma cosa abbiamo ottenuto fin ora? In realtà nulla, abbiamo solo trovato una nuova forma algebrica per esprimere la matrice dell’interazione.

Ovviamente, noi non siamo interessati alla decomposizione completa di cui sopra, ma alla decomposizione parziale, ottenuta estraendo le prime r componenti principali. Se utilizzando la prima componente principale (cioè la prima colonna di G ed E, che chiameremo G1 ed E1, possiamo calcolare una matrice d’interazione “semplificata” (approssimata):

PC <- 1
G1 <- G[,1:PC]
E1 <- E[,1:PC]
print ( G1 %*% t(E1), digits=2)
##      Chiaro_A Chiaro_P Colla_A Colla_P Palomb_A Palomb_P Scuro_A
## [1,]    0.504   -0.459 -0.0702  -0.670    0.267   -0.233   0.520
## [2,]   -0.027    0.024  0.0037   0.036   -0.014    0.012  -0.028
## [3,]   -0.407    0.370  0.0566   0.541   -0.215    0.188  -0.420
## [4,]    0.402   -0.365 -0.0559  -0.534    0.213   -0.186   0.414
## [5,]    0.042   -0.038 -0.0058  -0.055    0.022   -0.019   0.043
## [6,]   -0.514    0.467  0.0715   0.683   -0.272    0.238  -0.530
##      Scuro_P Sicania_A Sicania_P Vesuvio_A Vesuvio_P
## [1,]  -0.191    0.0522    0.1781     0.364    -0.263
## [2,]   0.010   -0.0028   -0.0094    -0.019     0.014
## [3,]   0.154   -0.0421   -0.1437    -0.294     0.212
## [4,]  -0.152    0.0416    0.1418     0.290    -0.209
## [5,]  -0.016    0.0043    0.0147     0.030    -0.022
## [6,]   0.194   -0.0532   -0.1814    -0.371     0.267

Gli elementi contenuti in questa matrice hanno una devianza pari a 5.29, che costituisce il 68.3% della devianza degli elementi in GEint.

sum(unlist(GEint^2))
## [1] 7.742996
sum(unlist (G1 %*% t(E1))^2)
## [1] 5.290174

Allo stesso modo, se utilizziamo le prime due componenti principali possiamo osservare che utilizzando anche la PC2, la matrice risultante ha una devianza pari 5.29 + 1.388. Di conseguenza, l’inclusione della seconda PC consente di spiegare un ulteriore 17.9% della devianza dell’interazione.

Quante componenti principali includere nell’analisi AMMI?

Come nella PCA, la scelta del numero di componenti principali da considerare può essere semplicemente basata sulla percentuale di devianza spiegata. Ciò rientra perfettamente nella logica di esplorazione dei dati (Exploratory Data Analysis) che ci spinge ad eseguire l’analisi AMMI.

Tuttavia, se ci si vuole muovere in una logica più tradizionale di test d’ipotesi, si può considerare che i gradi di libertà associati ad ogni PC sono dati da (Zobel, 1988):

$$ df = j + k - 1 + 2n $$

dove \textit{j} è il numero delle varietà, \textit{k} è il numero degli ambienti ed \textit{n} è il numero delle PC considerate. Quindi in questo caso la PC1 ha 15 gradi di libertà, mentre la PC2 ne ha 13. Inoltre, le devianze che abbiamo appena calcolato vanno moltiplicate per 4, dato che abbiamo quattro repliche.

Se si dispone delle repliche (e quindi di una stima dell’errore ‘puro’), le informazioni precedenti possono essere utilizzate nell’ANOVA per verificare quante PC inserire, tenendo conto che la prova aveva quattro repliche. In questo caso, non sono necessarie altre PC oltre le prime 2.

Secondo molti autori il test di F riportato nella tabella precedente è troppo poco conservativo e sono state proposte altre procedure più complesse (Cornelius, 1993).

Perchè l’analisi AMMI è utile? Il Biplot

A questo punto, l’utilità dell’analisi AMMI dovrebbe essere chiara: se la prima (o al limite le prime due) PC spiegano una gran parte della devianza dell’interazione, quest’ultima può essere descritta utilizzando una sola colonna (o al limite due) di G ed E. Di conseguenza G$_1$ ed E$_1$ possono essere facilmente plottati su un grafico.

In particolare, si può utilizzare un biplot comq quello in figura 1, dove possiamo rappresentare sulle ascisse la produzione media di ogni varietà e di ogni ambiente (gli effetti principali in t/ha, nel nostro caso) e sulle ordinate gli scores di ciascuna varietà ed ambiente (che a loro volta rappresentano l’interazione ‘genotipo x ambiente’. Il biplot diviene quindi uno strumento riassuntivo molto efficace, in quanto rappresenta un’altissima quota della devianza totale della prova (in questo caso il 91.6%).

La lettura del grafico è molto semplice: oltre al livello produttivo, le varietà e gli ambienti che hanno factor scores vicini a 0 non interagiscono (infatti il prodotto degli scores, cioè l’interazione nella sua forma semplificata, sarebbe vicino a 0).

Allo stesso modo, le varietà e gli ambienti che hanno factor scores dello stesso segno mostrano un interazione positiva (perchè tale è il prodotto di grandezze di ugual segno), tanto più alta quanto più alti sono i factor scores.

Viceversa per le varietà e gli ambienti che hanno factor scores di segno opposto interagiscono negativamente, in modo tanto più marcato quanto più alti sono i factor scores (in valore assoluto).

Figura 1. Biplot 1

Si può ancora osservare che il gruppo di varietà indicate con A sone mediamente più produttive di quelle indicate con P. Nel primo gruppo, Palombino, Vesuvio, Scuro e Chiaro di Torre Lama sono state le più produttive nel 2001 in entrambe le località (PC score per l’ambiente vicino a quello della varietà), mentre Collameno e Sicania si sono rivelate in media meno produttive, ma anche meno influenzate dall’ambiente (scores quasi pari a 0).

Nel secondo gruppo (P), Scuro di Torre Lama, Sicania e Chiaro di Torre Lama si sono rivelate le più produttive; nel caso di Sicania, le prestazioni sono state basse nel 2003 in entrambe le località (score per l’ambiente molto lontano dallo score di varietà) mentre Chiaro di torre Lama non è andata particolarmente bene nel 2001.

Per quanto riguarda gli ambienti, il 2003 è stato l’anno più favorevole in entrambe le località per le varietà del gruppo P, mentre il 2001 è stato favorevole per le varietà del gruppo A. Al contrario, il 2002 è stato un anno poco produttivo per tutte le varietà.

Insomma, un solo grafico può darci una quantità notevole di informazioni, se solo lo si riesce a leggere in modo appropriato!

Se volessimo concentrarci sulla sola interazione, trascurando gli effetti principali, potremmo usare un biplot analogo a quelli utilizzati nel PCA

biplot(ammi, scaling = 3)

Bibliografia essenziale

1. Annichiarico, P. (1997). Addittive main effects and multiplicative interaction (AMMI) analysis of genotype-location interaction in variety trials repeated over years. Theorethycal applied genetics, 94, 1072–1077.
2. Annichiarico, P. (1997). Joint regression analysis of genotype-environment interactions for cereals in Italy. Euphytica, 94, 53–62.
3. Annichiarico, P. (1997). Nuove metodologie per lo studio dell’interazione trattamento-ambiente e per la definizione della raccomandazione tecnica. Rivista di Agronomia,31, 817–823.
4. Annichiarico, P. (2000). Analisi della sperimentazione ripetuta in più ambienti. L’ Informatore agrario, 29–34.
5. Annichiarico, P., Bozzo, F., Parente, G., Gusmeroli, F., Mair, V., Marguerettaz, O., and Orlandi, D. (1995). Analysis of adaptation of grass/legume mixtures to Italian alpine and subalpine zones through an addittive main effects and multiplicative interaction model. Grass and Forage Science, 50, 405–413.
6. Ariyo, O. J. (1998). Use of additive main effects and multiplicative interaction model to analyse multilocation soybean varietal trials. J. Genet. and Breed, 129–134.
7. Cornelius, P. L. (1993). Statistical tests and retention of terms in the Addittive Main Effects and Multiplicative interaction model for cultivar trials. Crop Science, 33,1186–1193.
8. Crossa, J. (1990). Statistical Analyses of multilocation trials. Advances in Agronomy, 44, 55–85.
9. Gollob, H. F. (1968). A statistical model which combines features of factor analytic and analysis of variance techniques.Psychometrika, 33, 73–114.
10. Stagnari F., Onofri A., Jemison J., Monotti M. (2006). Multivariate analyses to discriminate the behaviour of faba bean (Vicia faba L. var. minor) varieties as affected by sowing time in cool, low rainfall Mediterranean environments. Agronomy For Sustainable Development, 27, 387–397.
11. Zobel, R. W., Wright, M.J., and Gauch, H. G. (1988). Statistical analysis of a yield trial. Agronomy Journal, 388–393.