Capitolo 9 Contrasti e confronti multipli

La scomposizione della varianza (ANOVA fisheriana) rappresenta frequentemente il primo passo nell’elaborazione statistica dei dati sperimentali. Essa consente di quantificare l’errore sperimentale e ci permette di sapere se l’effetto del trattamento (nel suo complesso) è risultato significativo. Tuttavia, con la sola ANOVA non possiamo ancora definire una graduatoria di merito tra i diversi trattamenti sperimentale e capire cosa è stato meglio di cosa. Per questo motivo, dopo l’ANOVA, l’analisi dei dati prosegue, di solito, con il confronto tra medie o tra gruppi di medie, per capire se vi siano differenze significative tra un trattamento e l’altro. Ovviamente, La significatività di queste differenze è solo uno degli aspetti che più ci interessano, in quanto è anche molto importante valutare la loro rilevanza biologica. Infatti una differenza potrebbe essere significativa, ma irrilevante da un punto di vista agronomico o, viceversa, essa potrebbe essere non significativa, ma estremamente rilevante, quindi meritevole di ulteriori approfondimenti scientifici.

Questa parte del lavoro viene usualmente eseguita utilizzando i contrasti lineari, che introdurremo utilizzando lo stesso dataset che abbiamo già iniziato ad analizzare nei due capitoli precedenti (‘mixture.csv’).

9.1 Esempio 1

Torniamo al nostro esperimento in cui abbiamo confrontato due erbicidi e la loro miscela con il testimone non trattato, in un disegno sperimentale completamente randomizzato con quattro repliche. Ricarichiamo il dataset ed eseguiamo l’ANOVA, come abbiamo visto nel Capitolo 7. Ricordiamo anche che, nel Capitolo 8, l’analisi grafica dei residui ha confermato che, per questo dataset, non vi sono problemi con le assunzioni di base per i modelli lineari.

repo <- "https://www.casaonofri.it/_datasets/"
file <- "mixture.csv"
pathData <- paste(repo, file, sep = "")
dataset <- read.csv(pathData, header = T)
model <- lm(Weight ~ Treat, data=dataset)
anova(model)
## Analysis of Variance Table
## 
## Response: Weight
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Treat      3 1089.53  363.18  23.663 2.509e-05 ***
## Residuals 12  184.18   15.35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

L’analisi della varianza ed il relativo test di F ci dicono che esiste una qualche differenza significativa tra i trattamenti, ma rimane il problema di classificare le soluzioni diserbanti in ordine di efficacia. Per prima cosa, calcoliamo le medie dei trattamenti, utilizzando la funzione emmeans() del package ‘emmeans’:

library(emmeans)
medie <- emmeans(model, ~Treat)
medie
##  Treat           emmean   SE df lower.CL upper.CL
##  Metribuzin__348   9.18 1.96 12     4.91     13.4
##  Mixture_378       5.13 1.96 12     0.86      9.4
##  Rimsulfuron_30   16.86 1.96 12    12.59     21.1
##  Unweeded         26.77 1.96 12    22.50     31.0
## 
## Confidence level used: 0.95

Vediamo che l’output di R riporta anche gli errori standard delle medie (SEM) e gli intervalli di confidenza.

9.2 I contrasti

Si definisce CONTRASTO una combinazione lineare delle medie dei trattamenti, in modo che i coefficienti sommati diano 0. Ad esempio, considerando i parametri del modello precedente, una combinazione lineare del tipo:

k=139.18+135.13+1316.86126.77=16.385

è un contrasto, in quanto la somma dei coefficienti è:

13+13+131=0

Al contrario una combinazione lineare del tipo:

k2=19.18+15.13+116.86126.77

non è un contrasto, perché la somma dei coefficienti non è zero (1 + 1 + 1 - 1 = 2).

Il primo contrasto k, ha un preciso significato biologico, in quanto stima la differenza tra il non diserbato e la media dei diserbati. Il risultato è -16.39, il che farebbe pensare che, nei vasetti diserbati vi siano, mediamente, meno infestanti che non in quelli non trattati. Tuttavia, noi non siamo interessati ai soli vasetti impiegati in prova, ma vogliamo trarre conclusioni di carattere generale. Sappiamo già che l’errore sperimentale produce incertezza sulla stima delle medie delle popolazioni, che si propaga al contrasto, per il quale dovremmo cercare di calcolare la deviazione standard e l’intervallo di confidenza.

La varianza di un contrasto può essere calcolata come combinazione lineare di varianze, attraverso il metodo di propagazione degli errori. Considerando che le medie sono, usualmente, indipendenti, la varianza di un contrasto tra medie è:

var(Aμ1+Bμ2)=(Aσμ1)2+(Bσμ2)2

dove A e B sono i coefficienti del contrasto, μ1 e μ2 sono due medie e σμ1 e σμ2 sono gli errori standard di μ1 e μ2. Nel nostro caso, la varianza del contrasto è:

var(k)=(13)23.83703+(13)23.83703+(13)23.83703+(1)23.83703=5.11604

mentre la deviazione standard (cioè l’errore standard) del contrasto è pari a:

SE(k)=5.11604=2.261866

Insomma, per il contrasto k abbiamo una stima puntuale (-16.4) e un errore standard (2.26), che può essere utilizzato per calcolare l’intervallo di confidenza del contrasto, secondo le metodiche usuali che abbiamo gia visto in un capitolo precedente. Ovviamente, se l’intervallo di confidenza di un contrasto contiene lo 0, siamo portati, intuitivamente, a concludere che il contrasto stesso non è significativo.

In modo più formale, potremmo chiederci: “E’ possibile che il contrasto, nella realtà, sia uguale a 0?”. Ovviamente è possibile: il nostro è solo un campione e, se ripetessimo il campionamento, potremmo ottenere valori di k totalmente diversi da quello effettivamente osservato. Poniamo l’ipotesi nulla in questi termini:

H0:κ=0

dove κ è il valore ‘vero’ del contrasto, per le popolazioni che hanno generato i dati. Scriviamo la statistica:

T=kES(k)=16.3852.261866=7.244

Se l’ipotesi nulla è vera, che probabilità abbiamo di osservare T = -7.244? Abbia già visto che il rapporto tra una stima ed il suo errore standard ha una distribuzione t di Student, con un numero di gradi di libertà pari a quelli del residuo dell’ANOVA. Di conseguenza possiamo saggiare l’ipotesi nulla che il contrasto è uguale a 0 calcolando la probabilità di trovare un valore di T pari o superiore (in valore assoluto) a quello da noi ottenuto. Nell’esempio sottostante abbiamo moltiplicato la probabilità trovata per 2, dato che si tratta di un test a due code:

Tval <- -16.385 / 2.261866
2 * pt(abs(Tval), 12, lower.tail = F)
## [1] 1.023001e-05

Possiamo quindi concludere che la differenza tra la media dei vasetti trattati e quella dei vasetti non trattati è pari a -16.34 g e ci sono elementi sufficienti per ritenere che essa sia diversa da 0.

9.3 I contrasti con R

Nel caso in esempio, si potrebbero pianificare quattro contrasti (incluso quello già discusso):

  1. trattato vs. non trattato (in media)
  2. miscela vs. erbicidi singoli (in media)
  3. miscela vs. rimsulfuron
  4. miscela vs. metribuzin

Per eseguire i contrasti con R, dobbiamo, in primo luogo, definire i vettori dei coefficienti. Per il primo contrasto, abbiamo già visto che questo vettore è:

k1 <- c(1/3, 1/3, 1/3, -1)

Per gli altri tre contrasti, i coefficienti sono:

k2 <- c(1/2, -1, 1/2, 0)
k3 <- c(0, -1, 1, 0)
k4 <- c(1, -1, 0, 0)

Una volta definiti i coefficienti, possiamo utilizzare il package ‘emmeans’, con la funzione contrast(), passandole, come argomento, l’oggetto ‘medie’, ottenuto con la funzione emmeans() (vedi sopra), ed una lista contenente i quattro vettori dei coefficienti (k1, k2, k3, k4).

K <- list(k1 = k1, K2 = k2, k3 = k3, k4 = k4)
contrast(medie, method = K, adjust="none")
##  contrast estimate   SE df t.ratio p.value
##  k1         -16.39 2.26 12  -7.244  <.0001
##  K2           7.89 2.40 12   3.289  0.0065
##  k3          11.73 2.77 12   4.235  0.0012
##  k4           4.05 2.77 12   1.461  0.1697

Il test mostra che il diserbo ha avuto, in media, un effetto significativo, che la miscela è più efficace dei trattamenti singoli, che la miscela è più efficace di rimsulfuron da solo, ma non è significativamente più efficace di metribuzin da solo.

9.4 I confronti multipli a coppie (pairwise comparisons)

Non sempre le prove sperimentali consentono di saggiare pochi contrasti pre-stabiliti, ma spesso è necessario confrontare tutte le possibili coppie di trattamenti (pairwise comparisons). In questo caso dovremmo definire un contrasto per ogni coppia di medie, anche se l’impiego del package ‘emmeans’ agevola, non di poco, il lavoro.

In particolare, possiamo immaginare due situazioni di riferimento: tutti contro tutti (confronti tipo “Tukey”) e tutti verso uno (confronti tipo “Dunnett”). Questo secondo tipo di contrasto può essere interessante, quando sia importante confrantare tutti i trattamenti verso un riferimento prescelto, ad esempio la miscela tra i due erbicidi.

Nel quadro sottostante mostriamo un confronto tipo Tukey (tutti contro tutti), eseguito utilizzando la funzione contrast() (come sopra) e passando il valore ‘pairwise’ all’argomento ‘method’. Vediamo che ci sono sei contrasti a coppie, tutti significativi, meno uno.

#Confronti multipli a coppie
contrast(medie, adjust="none", method="pairwise")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.1697
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0168
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  <.0001
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0012
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0038

Per i confronti del tipo ‘tutti verso uno’ è possibile utilizzare la stessa funzione, assegnando però il valore ‘dunnett’ (invece che ‘pairwise’) all’ argomento ‘method’.

contrast(medie, adjust="none", method="dunnett")
##  contrast                         estimate   SE df t.ratio p.value
##  Mixture_378 - Metribuzin__348       -4.05 2.77 12  -1.461  0.1697
##  Rimsulfuron_30 - Metribuzin__348     7.68 2.77 12   2.774  0.0168
##  Unweeded - Metribuzin__348          17.60 2.77 12   6.352  <.0001

Così facendo vediamo che R confronta tutte le tesi con metribuzin, che è il primo livello in ordine alfabetico, mentre noi volevamo confrontare tutte le tesi con la miscela. Per ottenere questo risultato basta aggiungere l’argomento ‘ref’ ed assgnare il valore ‘2’, considerando che la miscela è la seconda tesi in ordine alfabetico:

contrast(medie, adjust="none", method="dunnett", ref = 2)
##  contrast                      estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378     4.05 2.77 12   1.461  0.1697
##  Rimsulfuron_30 - Mixture_378     11.73 2.77 12   4.235  0.0012
##  Unweeded - Mixture_378           21.64 2.77 12   7.813  <.0001

Il risultato delle elaborazioni sovrastanti mostra i contrasti con il loro errore standard e potrebbe essere interessante calcolare anche l’intervallo di confidenza per le differenze stimate. Ciò può esser fatto facilmente assegnando il risultato della funzione contrast() ad una variabile ed esplorando quest’ultima con il metodo confint().

# Confidence intervals
con <- contrast(medie, adjust="none", 
                method="pairwise", ref = 2)
confint(con)
##  contrast                         estimate   SE df lower.CL upper.CL
##  Metribuzin__348 - Mixture_378        4.05 2.77 12    -1.99    10.08
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12   -13.72    -1.65
##  Metribuzin__348 - Unweeded         -17.60 2.77 12   -23.63   -11.56
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12   -17.77    -5.70
##  Mixture_378 - Unweeded             -21.64 2.77 12   -27.68   -15.61
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12   -15.95    -3.88
## 
## Confidence level used: 0.95

9.5 Display a lettere

I risultati di un confronto multiplo a coppie possono essere presentati anche con un display a lettere, nel quale le medie seguite da lettere diverse sono significativamente diverse per un livello di probabilità di errore minore di quello dato.

Realizzare un display a lettere manualmente è piuttosto facile, utilizzando la seguente procedura:

  1. ordinare le medie in senso crescente/decrescente
  2. partire dalla prima media e aggiungere la lettera A a tutte quelle che non sono significativamente diverse
  3. passare alla seconda media e aggiungere la lettera B a tutte quelle che non sono significativamente diverse
  4. procedere analogamente con tutte le medie successive, finche non vengono più individuate differenze significative.

Con R si può sfruttare il package ‘emmeans’, utilizzando i comandi sottostanti.

multcomp::cld(medie, adjust="none", Letters=LETTERS)
##  Treat           emmean   SE df lower.CL upper.CL .group
##  Mixture_378       5.13 1.96 12     0.86      9.4  A    
##  Metribuzin__348   9.18 1.96 12     4.91     13.4  A    
##  Rimsulfuron_30   16.86 1.96 12    12.59     21.1   B   
##  Unweeded         26.77 1.96 12    22.50     31.0    C  
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

9.6 Tassi di errore per confronto e per esperimento

Operando nel modo anzidetto, ogni contrasto/confronto ha una probabilità di errore del 5% (αC=0.05). Se i contrasti/confronti sono più di uno (‘famiglia’ di n contrasti), la probabilità di sbagliarne almeno uno (maximum experimentwise error rate: αE) è data da:

αE=1(1αC)n

Bisogna premettere che l’anzidetta formula vale quando i contrasti sono totalmente indipendenti tra loro, cosa che quasi mai avviene, dato che, anche in un semplice modello ANOVA, i contrasti condividono la stessa varianza d’errore e sono quindi più o meno correlati tra di loro. Con contrasti non indipendenti la formula anzidetta fornisce una sovrastima di αE (per questo si parla di maximum experimentwise error rate).

Il numero di confronti a coppie per esperimento può essere anche molto elevato: se ho k medie il numero dei confronti possibili è pari a k(k1)/2. Di conseguenza, la probabilità di errore per esperimento (αE) può essere molto più alta del valore αC prefissato per confronto.

Ad esempio, se ho 15 medie, ho (1514)/2=105 confronti possibili. Se uso αC=0.05 per ogni confronto, la probabilità di sbagliarne almeno uno è pari (in caso di confronti indipendenti) a 1(10.05)105=0.995. Sostanzialmente, vi è pressoché la certezza che in questo esperimento qualcosa sia sbagliato!

Per questo motivo, quando si elaborano i dati di un esperimento nel quale è necessario fare molti contrasti, o confronti, o, più in generale, molti test d’ipotesi simultanei, si potrebbe voler esprimere un giudizio globale (simultaneo) sull’intera famiglia di contrasti/confronti, minimizzando la possibilità che anche solo uno o pochi di essi siano sbagliati. In particolare, ciò potrebbe capitare quando:

  1. vogliamo trovare i migliori di k trattamenti, senza correre rischi di escluderne erroneamente qualcuno. In questa situazione, facendo ogni confronto con il 5% di probabilità di errore, la probabilità di escludere erroneamente anche solo un trattamento dal lotto dei migliori è molto più alta di quella prefissata, perché basta sbagliare anche uno solo dei k - 1 confronti con il migliore.
  2. Abbiamo utilizzato un display a lettere e intendiamo affermare che ‘i trattamenti seguiti da lettere diverse sono significativamente diversi’. In questo caso, stiamo tirando una conclusione basata sull’intera famiglia di confronti e non possiamo lecitamente riportare la probabilità di errore di un singolo confronto.

9.7 Aggiustamento per la molteplicità

In tutte le condizioni analoghe a quelle più sopra accennate si pone il problema di applicare un aggiustamento per la molteplicità, in modo da rispettare un certo livello prestabilito di protezione per esperimento (e non per confronto). La prima possibilità che abbiamo è quella di aggiustare il P-level per ogni confronto, in modo da diminuire la probabilità di errore per l’intera famiglia di sei confronti. Utilizzando la formula che lega la probabilità d’errore per esperimento (αE) alla probabilità d’errore per confronto (αC; vedi sopra), possiamo prendere la sesta colonna del dataframe ‘con’, quella che contiene i P-levels non corretti, e trasformarla come segue:

alphaC <- as.data.frame(con)[,6]
1 - (1 - alphaC)^6
## [1] 6.722991e-01 9.683462e-02 2.190543e-04 6.923077e-03 2.869757e-05
## [6] 2.255183e-02

Più facilmente, possiamo arrivare allo stesso risultato con il package ‘emmeans’:

contrast(medie, method = "pairwise", adjust = "sidak")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.6723
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0968
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0069
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0226
## 
## P value adjustment: sidak method for 6 tests

Vediamo che il secondo confronto, che era significativo, non è più tale adottando la correzione di Sidak.

Un’alternativa più nota (e semplice) è quella di utilizzare la diseguaglianza di Bonferroni:

αE=αCk

Quest’ultima è un po’ più conservativa della precedente, nel senso che fornisce un P-level aggiustato leggermente più alto dell’altra.

alphaC * 6
## [1] 1.018071e+00 1.009900e-01 2.190743e-04 6.943132e-03 2.869792e-05
## [6] 2.276671e-02

Oppure possiamo utilizzare la funzione contrast():

contrast(medie, method = "pairwise", adjust = "bonferroni")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  1.0000
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.1010
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0069
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0228
## 
## P value adjustment: bonferroni method for 6 tests

Esistono altre procedure di aggiustamento del P-level (metodi di Holm, Hochberg, Hommel), ma nessuna di queste tiene conto della correlazione eventualmente esistente tra i contrasti e tutte quindi sono da definirsi più o meno ‘conservative’.

Oltre che aggiustare il P-level, possiamo anche utilizzare altre procedure di aggiustamento, basate sulla distribuzione multivariata di t e in grado di tener conto dell’eventuale correlazione esistente tra i contrasti stessi. Una di queste procedure viene impiegata di default nella funzione contrast() nel package ‘emmeans’:

#Confronti multipli a coppie, basati sul t multivariato
contrast(medie, method="pairwise")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.4885
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0698
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0055
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0173
## 
## P value adjustment: tukey method for comparing a family of 4 estimates

Possiamo notare che i P-levels sono leggermente più bassi di quelli ottenuti con Bonferroni, che conferma quindi di essere una procedura molto conservativa, mentre l’impiego del t multivariato consente di rispettare esattamente il tasso di errore ‘per esperimento’.

Ovviamente la correzione per la molteplicità ed il conseguente innalzamento del P-level sono fortemente dipendenti dal numero di contrasti effettuati; se utilizziamo un set di confronti tipo ‘dunnett’ invece che tipo ‘pairwise’, avremo meno confronti e quindi una correzione più ‘leggera’.

#Confronti multipli a coppie, basati sul t multivariato
contrast(medie, method="dunnett")
##  contrast                         estimate   SE df t.ratio p.value
##  Mixture_378 - Metribuzin__348       -4.05 2.77 12  -1.461  0.3711
##  Rimsulfuron_30 - Metribuzin__348     7.68 2.77 12   2.774  0.0442
##  Unweeded - Metribuzin__348          17.60 2.77 12   6.352  0.0001
## 
## P value adjustment: dunnettx method for 3 tests

9.8 Esempio 2: confronti multipli con dati trasformati

Un altro aspetto di cui tener conto è relativo al fatto che, nel caso in cui si fosse reso necessario qualche tipo di trasformazione stabilizzante, il test di confronto multiplo deve opportunamente adeguato. Per capire questo aspetto, possiamo riprendere in mano il dataset ‘insects.csv’ già esplorato nel capitolo precedente e relativo alle quindici piante trattate con tre insetticidi (cinque piante per insetticida), su ciascuna delle quali, alcune settimane dopo il trattamento, sono stati contati gli insetti presenti e vitali. Ricarichiamo il dataset e, tenendo conto dell’esigenza emerse nel capitolo precedente, analizziamo i dati previa trasformazione logaritmica. Quest’ultima viene eseguita direttamente nel corso del ‘model fitting’ e non nel dataset, per un motivo che sarà chiaro in seguito.

fileName <- "https://www.casaonofri.it/_datasets/insects.csv"
dataset <- read.csv(fileName, header = T)
dataset$Insecticide <- factor(dataset$Insecticide)
mod <- lm(log(Count) ~ Insecticide, data = dataset)

A questo punto, possiamo calcolare le medie marginali attese, che, tuttavia, sono su scala logaritmica, come indicato nel codice sottostante.

medie <- emmeans(mod, ~Insecticide)
medie
##  Insecticide emmean    SE df lower.CL upper.CL
##  T1            6.34 0.178 12     5.96     6.73
##  T2            5.81 0.178 12     5.43     6.20
##  T3            3.95 0.178 12     3.56     4.34
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Presentare le medie su scala logaritmica, in molti casi, potrebbe non essere di immediata o facile lettura. Per questo, potremmo pensare di retro-trasformare le medie, utilizzando la trasformazione inversa di quella logaritmica, cioè l’antilogaritmo. Ad esempio, per la prima media:

e6.34=566.7963

In questo modo la nostra unità di misura ridiventa quella originale, anche se il valore ottenuto non coincide con la media dei dati originali; in effetti la trasformazione è non lineare e la media dei logaritmi non può coincidere con il logaritmo della media. Possiamo osservare che la media del trattamento A, sulla scala originale, è:

mean(dataset[dataset$Insecticide == "T1","Count"])
## [1] 589.8

e risulta più alta della media retro-trasformata. In realtà, se è vero che i logaritmi sono normalmente distribuiti, la media dei logaritmi (6.34) dovrebbe essere uguale alla mediana (ricordiamo che in una distribuzione normale media e mediana coincidono). La mediana è il valore centrale; dato che la retro-trasformazione è monotona, il valore centrale resta centrale, anche se io retro-trasformo. Quindi la media retro-trasformata è uno stimatore della mediana della popolazione originale, non della sua media. Questo non è uno svantaggio: infatti il QQ-plot suggerisce un’asimmetria positiva (vedi capitolo precedente) cosa che è confermata dal fatto che la mediana è minore della media. Se la distribuzione dei dati è asimmetrica, la mediana è un indicatore di tendenza centrale migliore della media, perché meno sensibile ai valori estremi, che sono più frequenti in caso di asimmetria.

Il problema è che, se vogliamo utilizzare la media retro-trasformata, dobbiamo trovare un valore di errore standard per questo stima, con il quale esprimere la sua incertezza. In realtà, anche l’errore standard può essere retro-trasformato, con una tecnica detta metodo ‘delta’, che costituisce un estensione della legge di propagazione degli errori per le trasformazioni non-lineari. È inutile andare nel dettaglio; diciamo solo che la funzione emmeans() rende semplicissima l’implementazione del metodo delta, con il comando seguente:

retroMedie <- emmeans(mod, ~Insecticide, type = "response")
retroMedie
##  Insecticide response     SE df lower.CL upper.CL
##  T1             568.6 101.01 12    386.1    837.3
##  T2             335.1  59.54 12    227.6    493.5
##  T3              51.9   9.22 12     35.2     76.4
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Con questo abbiamo tutto quello che ci serve: stime ed errori standard, che, ovviamente, sono diversi per le diverse medie retro-trasformate, coerentemente con la mancanza di omoscedasticità. Il test di confronto multiplo è, pertanto:

library(multcomp)
cld(retroMedie, Letters = LETTERS)
##  Insecticide response     SE df lower.CL upper.CL .group
##  T3              51.9   9.22 12     35.2     76.4  A    
##  T2             335.1  59.54 12    227.6    493.5   B   
##  T1             568.6 101.01 12    386.1    837.3   B   
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

La retrotrasformazione è immediata in quanto la funzione emmeans() è stata in grado di riconoscere la trasformazione stabilizzante applicata all’interno del ‘model fitting’. La situazione diviene più complessa se la trasformazione è applicata prima del ‘model fitting’ oppure se non viene automaticamente riconosciuta. Ad esempio, una trasformazione inversa non viene automaticamente riconosciuta e, di conseguenza, la retrotrasformazione non viene effettuata.

mod2 <- lm(1/Count ~ Insecticide, data = dataset)
emmeans(mod2, ~Insecticide, type = "response")
##  Insecticide  emmean      SE df lower.CL upper.CL
##  T1          0.00182 0.00279 12 -0.00426  0.00789
##  T2          0.00317 0.00279 12 -0.00291  0.00924
##  T3          0.02120 0.00279 12  0.01513  0.02727
## 
## Confidence level used: 0.95

In questo caso, dobbiamo eseguire la trasformazione prima del ‘model fitting’ e poi cambiare la griglia di riferimento per il modello (update(ref_grid(mod2), ....)) specificando la trasformazione effettuata (tran = "inverse"). La griglia così modificata viene passata al posto del modello alla funzione emmeans(), come indicato di seguito.

dataset$InvCount <- 1/dataset$Count
mod3 <- lm(InvCount ~ Insecticide, data = dataset)
updGrid <- update(ref_grid(mod3), tran = "inverse")
emmeans(updGrid, ~Insecticide, type = "response")
##  Insecticide response    SE df lower.CL upper.CL
##  T1             550.9 845.9 12    126.8      Inf
##  T2             315.8 278.0 12    108.2      Inf
##  T3              47.2   6.2 12     36.7       66
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the inverse scale

Questo metodo si può utilizzare con numerose funzioni, come, ad esempio “identity”, “1/mu^2”, “inverse”, “reciprocal”, “log10”, “log2”, “asin.sqrt” e “asinh.sqrt”.

Un approccio analogo può essere utilizzato se vogliamo effettuare la trasformazione di Box-Cox nella sua forma semplificata (W=Yλ), utilizzando un valore λ diverso da quelli ‘semplici’ (cioè diverso da 1, 0, 0.5 o -1). Ad esempio, se vogliamo utilizzare λ=0.25, dobbiamo operare in modo simile a quello descritto in precedenza, utilizzando la funzione make.tran(), come evidenziato nel box seguente.

dataset$lCount <- dataset$Count^0.25
mod4 <- lm(lCount ~ Insecticide, data = dataset)
updGrid <- update(ref_grid(mod4), tran = make.tran("power", 0.25))
emmeans(updGrid, ~Insecticide, type = "response")
##  Insecticide response   SE df lower.CL upper.CL
##  T1             573.5 78.7 12    420.3    765.3
##  T2             340.5 53.2 12    238.5    472.1
##  T3              53.1 13.2 12     29.6     88.3
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the mu^0.25 scale

9.9 E le classiche procedure di confronto multiplo?

Il confronto multiplo tradizionale è basato sul calcolo di una differenza critica minima, da utilizzare come base per il confronto tra due medie. In pratica, due medie sono considerate significativamente diverse quando la loro differenza supera la differenza critica. Nella letteratura scientifica si trovano molte procedure di confronto multiplo, tra le quali segnaliamo:

  1. Minima Differenza Significativa (MDS o LSD)
  2. Honest Significant Difference di Tukey
  3. Test di Dunnett
  4. Test di Duncan
  5. Test di Newman-Keuls
  6. Test di confronto multiplo di Tukey

Il primo metodo corrisponde esattamente a quello che abbiamo utilizzato all’inizio, per fare confronti multipli ‘tutti contro tutti’, senza correzione per la molteplicità. Il secondo e il terzo metodo corrispondono rispettivamente al test di confronto ‘tutti verso tutti’ e ‘uno verso tutti’ indicati in precedenza, con correzione per la molteplicità.

Non è necessario dettagliare gli altri test, in quanto, seppur siano ancora molto utilizzati, vengono ormai ritenuti obsoleti e sconsigliabili, da parecchi ricercatori. Chi vuole, trova altre informazioni nella letteratura indicata in fondo al capitolo.

9.10 Consigli pratici

La cosa fondamentale è muoversi in coerenza con le finalità dell’esperimento. Si consiglia di:

  1. Quando è possibile, pianificare gli esperimenti in modo da ottenere le risposte cercate con pochi contrasti di interesse. In questo modo il problema della molteplicità è minimizzato.
  2. Non usare mai contrasti con serie di dati quantitative. In questo caso la regressione è l’approccio corretto e ne parleremo in un prossimo capitolo. In generale, utilizzare i contrasti solo se sono coerenti con la logica dell’esperimento.
  3. Pianificare esattamente il numero di contrasti necessari ed eseguirli, fornendo il valore del contrasto e il suo errore standard.
  4. Decidere è necessario aggiustare il p-level (e gli intervalli di confidenza) per la molteplicità (tasso di errore comparisonwise o experimentwise).
  5. Se si decide di aggiustare il p-level, considerare che le procedure di Bonferroni o Sidak possono essere eccessivamente protette. Preferire quindi le procedure di aggiustamento basate sulla distribuzione t multivariata, il che, a livello di confronto multiplo con dati bilanciati, è equivalente ad utilizzate la Tukey HSD o il test di Dunnett.
  6. Evitare le procedure di Duncan e Newmann-Keuls: non danno il livello di protezione cercato e, inoltre, non sono basate su una differenza critica costante (quindi sono difficili da discutere).

9.11 Altre letture

  1. Hsu, J., 1996. Multiple comparisons. Theory and methods. Chapman and Hall.
  2. Bretz, F., T. Hothorn, and P. Westfall. 2002, On Multiple Comparison Procedures in R. R News 2: 14-17.
  3. Calinsky, T. and L. C. A. Corsten. 1985, Clustering means in ANOVA by simultaneus testing. Biometrics, 41, 39-48.
  4. Cargnelutti, A. F., L. Storck, L. A. Dal’Col, L. Pisaroglo de Carvalho, and P. Machado dos Santos. 2003, A precisao experimental relacionada ao uso de bordaduras nas extremidades das fileiras em ensaios de milho. Ciencia Rural 33: 607-614.
  5. Carmer, S. G. and M. R. Swanson. 1971, Detection of difference between mens: a Monte Carlo study of five pairwise multiple comparison procedures. Agronomy Journal, 63, 940-945.
  6. Carmer, S. G. 1976, Optimal significance levels for application of the least significant difference in crop performance trials. Crop Science, 16, 95-99.
  7. Chew, V. 1976, Comparing treatment means: a compendium. Hortscience, 11(4), 348-357.
  8. Cousens, R. 1988, Misinterpretetion of results in weed research through inappropriate use of statistics. Weed Research, 28, 281-289.
  9. Edwards, A. W. F. and L. L. Cavalli-Sforza. 1965, A method for cluster analysis. Biometrics, 21, 362-375.
  10. Gates, C. E. and J. D. Bilbro. 1978, Illustration of a cluster analysis method for mean separation. Agronomy journal, 70, 462-465.
  11. O’Neill, R. and G. B. Wetherill. 1971, The present state of multiple comparison methods. Journ. Roy. Stat. Soc. , 2, 218-249.
  12. Petersen, R. G. 1977, Use and misuse of multiple comparison procedurs. Agronomy Journal, 69, 205-208.
  13. Scott, A. J. and M. Knott. 1974, A cluster analysis method for grouping means in the analysis of variance. Biometrics, 30, 507-512.
  14. Weisberg, S., 2005. Applied linear regression, 3rd ed. John Wiley & Sons Inc. (per il metodo ‘delta’)
  15. Willavize, S. A., S. G. Carmer, and W. M. Walker. 1980, Evaluation of cluster analysis for comparing treatment means. Agronomy journal, 72,317-320.