Capitolo 8 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 siamo ancora capaci di definire una graduatoria di merito tra i diversi livelli del trattamento sperimentale. Dopo l’ANOVA, è quindi logico chiedersi (in ordine crescente di importanza):

  1. Il generico trattamento A ha un effetto diverso dal trattamento B?
  2. Il trattamento A è migliore/uguale/peggiore di B?
  3. Quant’è la differenza tra A e B?

La prima domanda è abbastanza sciocca, in quanto due trattamenti non sono quasi mai esattamente uguali. Le altre due domande sono invece più rilevanti, specialmente la seconda, in quanto conoscere l’entità della differenza, oltre che la sua significatività, è fondamentale per comprendere anche la sua 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. Per rispondere alle domande precedenti, in genere, utilizziamo i contrasti lineari, che introdurremo con un esempio.

8.1 Esempio

Ammettiamo di aver effettuato una prova con un trattamento sperimentali caratterizzato da quattro livelli qualitativi (tesi di concimazione), con i risultati riportati di seguito:

yield <- c(20,21,23,22,19,20,12,15,13,19,18,16)
fert <- factor(rep(c("Minerale", "Minerale lento", 
          "Non concimato", "Organico"), each=3))
dataset <- data.frame(yield, fert)
rm(yield, fert)
dataset
##    yield           fert
## 1     20       Minerale
## 2     21       Minerale
## 3     23       Minerale
## 4     22 Minerale lento
## 5     19 Minerale lento
## 6     20 Minerale lento
## 7     12  Non concimato
## 8     15  Non concimato
## 9     13  Non concimato
## 10    19       Organico
## 11    18       Organico
## 12    16       Organico

L’ANOVA può essere eseguita facilmente con R, utilizzando l’ormai nota funzione ‘lm()’.

model <- lm(yield ~ fert, data=dataset)
anova(model)
## Analysis of Variance Table
## 
## Response: yield
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## fert       3 115.000  38.333  16.429 0.0008821 ***
## Residuals  8  18.667   2.333                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In questo caso, l’analisi della varianza ed il relativo test di F ci dicono che esiste una differenza significativa tra le medie, ma rimane il problema di classificare le soluzioni concimanti in ordine di efficacia. Chi è migliore di chi?

Per prima cosa, calcoliamo le medie dei trattamenti, utilizzando la funzione emmeans() del package emmeans:

library(emmeans)
medie <- emmeans(model, ~fert)
medie
##  fert           emmean    SE df lower.CL upper.CL
##  Minerale         21.3 0.882  8     19.3     23.4
##  Minerale lento   20.3 0.882  8     18.3     22.4
##  Non concimato    13.3 0.882  8     11.3     15.4
##  Organico         17.7 0.882  8     15.6     19.7
## 
## Confidence level used: 0.95

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

8.2 I contrasti

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

\[C = \frac{1}{3} \cdot 21.33 + \frac{1}{3} \cdot 20.33 - 1 \cdot 13.33 + \frac{1}{3} \cdot 17.67 = 6.446667\]

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

\[C = \frac{1}{3} + \frac{1}{3} - 1 + \frac{1}{3} = 0\]

Al contrario una combinazione lineare del tipo:

\[C1 = 1 \cdot 21.33 + 1 \cdot 20.33 - 1 \cdot 13.33 + 1 \cdot 17.67\]

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

Il primo contrasto C, ha un preciso significato biologico, in quanto stima la differenza tra il non fertilizzato e la media dei fertilizzati. Il risultato è 6.45 e, con esso, possiamo rispondere a tutte e tre le domande elencato all’inizio:

  1. La fertilizzazione (in media) ha un effetto diverso dal testimone non fertilizzato? Risposta: si, perché la differenza è non nulla
  2. La fertilizzazione in media è migliore/uguale/peggiore del testimone? Risposta: migliore, perché la differenza è positiva.
  3. Quant’è la differenza tra il fertilizzato e il non-fertilizzato? Risposta: è pari a 6.45

E’ evidente, tuttavia, che l’errore sperimentale produce incertezza, che sarebbe bene includere nei risultati, sotto forma di deviazione standard (errore standard) del contrasto, oppure come intervallo di confidenza.

8.2.1 Varianza del contrasto e test d’ipotesi

Ogni contrasto ha la sua varianza, ottenuta 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 è:

\[\textrm{var}(A \mu_1 + B \mu_2) = (A \cdot \sigma_{\mu_1} )^{2} + (B \cdot \sigma_{\mu_2} ) ^ 2\]

dove A e B sono i coefficienti del contrasto, \(\mu_1\) e \(\mu_2\) sono due medie e \(\sigma_{\mu_1}\) e \(\sigma_{\mu_2}\) sono gli errori standard di \(\mu_1\) e \(\mu_2\). Nel nostro caso, la varianza del contrasto è:

\[var(C) = \left( \frac{1}{3} \cdot 0.882 \right)^2 + \left( \frac{1}{3} \cdot 0.882 \right)^2 + \left( 1 \cdot 0.882 \right)^2 + \left( \frac{1}{3} \cdot 0.882 \right)^2 = 1.037\]

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

\[ ES(C) = \sqrt{1.037} = 1.0183\]

Insomma, per il contrasto C abbiamo una stima puntuale (6.45) e un errore standard (1.0183). Utilizzando questo errore standard possiamo calcolare l’intervallo di confidenza del contrasto, secondo le metodiche usuali che abbiamo gia visto in un capitolo precedente. Gli intervalli di confidenza dei contrasti sono trattati in modo più dettagliato in fondo al capitolo.

A questo punto 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 C totalmente diversi da quello effettivamente osservato. Poniamo l’ipotesi nulla in questi termini:

\[H_0: K = 0\]

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

\[T = \frac{C}{ES(C)} = \frac{6.45}{1.0183} = 6.331\]

Se l’ipotesi nulla è vera, che probabilità abbiamo di osservare T = 6.331? In generale, 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:

t <- 6.4467 / 1.018277
2 * pt(t, 8, lower.tail = F)
## [1] 0.0002251554

Con questo abbiamo risposto a tutte e tre le nostre domande di ricerca:

  1. Fertilizzare non produce lo stesso effetto che non fertilizzare;
  2. Fertilizzare permette di produrre di più che non fertilizzare
  3. la differenza è pari a 6.45 q/ha e ci sono elementi sufficienti per ritenere che essa sia diversa da 0.

8.2.2 I contrasti con R

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

  1. non concimato vs concimato (in media)
  2. concime organico vs. concimi minerali (in media)
  3. minerale tradizionale vs. lento rilascio.

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 è:

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

Per gli altri due contrasti, i coefficienti sono:

m2 <- c(0.5, 0.5, 0, -1)
m3 <- 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 tre vettori dei coefficienti (m1, m2, m3), ai quali si può assegnare un nome, ad esempio, C1, C2 e C3.

m1 <- c(1/3, 1/3,  -1, 1/3)
m2 <- c(0.5, 0.5, 0, -1)
m3 <- c(1, -1, 0, 0)
contrast(medie, method=list(C1=m1, C2=m2, C3=m3), 
           adjust="none")
##  contrast estimate   SE df t.ratio p.value
##  C1           6.44 1.02  8 6.328   0.0002 
##  C2           3.17 1.08  8 2.932   0.0189 
##  C3           1.00 1.25  8 0.802   0.4458

Il test mostra che la concimazione ha, in media, un effetto significativo, che il concime organico differisce significativamente dai minerali (in media) e che il concime a lento rilascio non è significativamente diverso da quello normale.

8.3 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, nel nostro caso, per verificare quale dei concimi consenta un incremento di produzione rispetto al testimone non concimato, qualora non sia necessario confrontare tra di loro le tesi concimanti.

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

#Confronti multipli a coppie
contrast(medie, adjust="none", method="pairwise")
##  contrast                       estimate   SE df t.ratio p.value
##  Minerale - Minerale lento          1.00 1.25  8  0.802  0.4458 
##  Minerale - Non concimato           8.00 1.25  8  6.414  0.0002 
##  Minerale - Organico                3.67 1.25  8  2.940  0.0187 
##  Minerale lento - Non concimato     7.00 1.25  8  5.612  0.0005 
##  Minerale lento - Organico          2.67 1.25  8  2.138  0.0650 
##  Non concimato - Organico          -4.33 1.25  8 -3.474  0.0084

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
##  Minerale lento - Minerale    -1.00 1.25  8 -0.802  0.4458 
##  Non concimato - Minerale     -8.00 1.25  8 -6.414  0.0002 
##  Organico - Minerale          -3.67 1.25  8 -2.940  0.0187

Purtroppo vediamo che R confronta tutte le tesi con ‘Minerale’ (primo livello in ordine alfabetico), mentre noi volevamo confrontare tutte le tesi con ‘non concimato’. Per ottenere questo risultato dobbiamo ricodificare il vettore ‘fert’, in modo da definire ‘Non concimato’ come livello di riferimento. Per far questo si utilizza la funzione ‘relevel()’:

dataset$fert <- relevel(dataset$fert, ref="Non concimato")

A questo punto dobbiamo ripetere integralmente le analisi:

model <- lm(yield ~ fert, data=dataset)
medie <- emmeans(model, ~fert)
contrast(medie, adjust="none", method="dunnett")
##  contrast                       estimate   SE df t.ratio p.value
##  Minerale - Non concimato           8.00 1.25  8 6.414   0.0002 
##  Minerale lento - Non concimato     7.00 1.25  8 5.612   0.0005 
##  Organico - Non concimato           4.33 1.25  8 3.474   0.0084

8.4 Display a lettere

I risultati di un confronto multiplo a coppie (pairwise) 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 letter 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

Con R si può sfruttare il package ‘emmeans’, con la sintassi sottostante.

multcomp::cld(medie, adjust="none", Letters=LETTERS)
##  fert           emmean    SE df lower.CL upper.CL .group
##  Non concimato    13.3 0.882  8     11.3     15.4  A    
##  Organico         17.7 0.882  8     15.6     19.7   B   
##  Minerale lento   20.3 0.882  8     18.3     22.4   BC  
##  Minerale         21.3 0.882  8     19.3     23.4    C  
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05

8.5 Problemi di molteplicità: tassi di errore per confronto e per esperimento

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

\[\alpha_E = 1 - (1 - \alpha_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 \(\alpha_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 \cdot (k-1)/2\). Di conseguenza, la probabilità di errore per esperimento (\(\alpha_E\)) può essere molto più alta del valore \(\alpha_C\) prefissato per confronto.

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

8.5.1 Correzione per la molteplicità

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. Vediamo alcuni esempi di quando questo potrebbe capitare.

  1. Non vogliamo correre rischi di escludere erroneamente alcun trattamento dal lotto dei migliori. Infatti, poniamo di voler trovare i migliori di k trattamenti, intendendo con ciò quelli che non sono significativamente inferiori a nessun altro. 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.

In tutte le condizioni analoghe a quelle più sopra accennate si pone il problema di aggiustare il p-level di ogni contrasto in modo da rispettare un certo livello prestabilito di protezione per esperimento (e non per confronto). Per far questo, si utilizzano metodi che, invece che essere basati sulla distribuzione di t di Student, estendono quest’ultima al caso multivariato. Non darò dettagli, in quanto i calcoli non sono eseguibili a mano; tuttavia è importante sottolineare che, in R, la correzione per la molteplicità viene eseguita di default, basta rimuovere dai comandi l’argomento ‘correct=“none”’.

contrast(medie, method="pairwise")
##  contrast                       estimate   SE df t.ratio p.value
##  Non concimato - Minerale          -8.00 1.25  8 -6.414  0.0009 
##  Non concimato - Minerale lento    -7.00 1.25  8 -5.612  0.0022 
##  Non concimato - Organico          -4.33 1.25  8 -3.474  0.0342 
##  Minerale - Minerale lento          1.00 1.25  8  0.802  0.8518 
##  Minerale - Organico                3.67 1.25  8  2.940  0.0724 
##  Minerale lento - Organico          2.67 1.25  8  2.138  0.2203 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
contrast(medie, method="dunnett")
##  contrast                       estimate   SE df t.ratio p.value
##  Minerale - Non concimato           8.00 1.25  8 6.414   0.0006 
##  Minerale lento - Non concimato     7.00 1.25  8 5.612   0.0014 
##  Organico - Non concimato           4.33 1.25  8 3.474   0.0219 
## 
## P value adjustment: dunnettx method for 3 tests

Possiamo notare che i p-values sono più alti di quelli che abbiamo ottenuto senza correzione; in questo caso stiamo rispettando il nostro livello di protezione a livello di esperimento, non a livello di singolo confronto.

8.6 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 in fondo al capitolo.

8.7 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).

8.8 Per approfondire un po’…

8.8.1 Intervallo di confidenza di un contrasto

In precedenza, abbiamo calcolato la significatività dei contrasti con un test t di Student. Questo approccio all’analisi dei dati, estremamente comune nella pratica sperimentale, è da taluni ritenuto sub-optimale. Ad esempio, secondo il famoso statistico John Tukey, testare la significatività di un contrasto è sciocco3 almeno per due motivi:

  1. la domanda non è realistica: due trattamenti diversi o due gruppi di trattamenti diversi non possono che dare risultati diversi, magari in modo impercettibile, ma pur sempre diversi;
  2. l’eventuale rifiuto dell’ipotesi nulla non ci da nessuna informazione sulla rilevanza biologica della differenza, che è indipendente dalla sua significatività.

Pertanto, sempre secondo Tukey, è molto più rilevante parlare di effect size, cioè di ampiezza dell’effetto, da quantificare tramite un intervallo di confidenza. Le formule sono quelle usuali, tramite i quantili della distribuzione t di Student, con un numero di gradi di libertà pari a quello del residuo ANOVA. Ad esempio, per il primo contrasto, l’intervallo di confidenza è:

limSup <- 6.4467 + qt(0.975, 8) * 1.018277
limInf <- 6.4467 - qt(0.975, 8) * 1.018277
limInf; limSup
## [1] 4.098549
## [1] 8.794851

Con R, possiamo utilizzare la funzione ‘confint()’, passandole l’oggetto ‘medie’, ottenuto come output della funzione ‘emmeans()’:

confint(contrast(medie, method=list(C1=m1, C2=m2, C3=m3), 
           adjust="none"))
##  contrast estimate   SE df lower.CL upper.CL
##  C1         -2.889 1.02  8    -5.24   -0.541
##  C2         -0.333 1.08  8    -2.82    2.157
##  C3         -8.000 1.25  8   -10.88   -5.124
## 
## Confidence level used: 0.95

L’uso degli intervalli di confidenza può essere preferibile al test d’ipotesi formale perché ci fa vedere l’entità degli effetti; ad esempio, possiamo vedere che la differenza tra il concime tradizionale e quello a lento rilascio (contrasto C3), anche se non significativa, potrebbe essere rilevante da un punto di vista agronomico (3.88 q/ha). In altri casi, differenze significative potrebbero risultare biologicamente irrilevanti.

Insomma, l’analisi statistica moderna sta un po’ spostando l’attenzione dal P-level e dal test d’ipotesi formale, in favore di una più attenta valutazione dell’ ‘effect size’ e della rilevanza biologica delle differenze. Si tratta di una tendenza che va tenuta nella debita considerazione.

8.8.2 Correzione per la molteplicità: qualche dettaglio ulteriore

In questo capitolo, abbiamo mostrato come è possibile, con R, introdurre la correzione per la molteplicità, ma pensiamo che sia necessario fornire alcune informazioni aggiuntive. Prendiamo, per esempio, il dataset ‘mixture’, che abbiamo già utilizzato nel capitolo 7, che riguarda un esperimento nel quale sono poste a confronto quattro soluzioni erbicide per il diserbo del pomodoro. Le quattro medie, con gli errori standard, sono riproposte qui sotto.

library(aomisc)
library(emmeans)
data(mixture)
mod <- lm(Weight ~ Treat, data = mixture)
medie <- emmeans(mod, ~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

I confronti multipli possibili sono sei (non molti quindi…) e sono eseguiti senza alcuna correzione per la molteplicità, con la funzione ‘contrast()’.

confronti <- contrast(medie, method = "pairwise", 
                      adjust = "none")
confronti
##  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

Immaginiamo che sia necessario adottare una correzione per la molteplicità, anche se il numero di confronti non è elevatissimo. La prima possibilità è quella di aggiustare il P-level per ogni confronto, in modo da diminuire la probabilità di errore per l’intera famiglia di sei confronti. In questo capitolo, abbiamo visto che la probabilità d’errore per esperimento (\(\alpha_E\)) dipende dalla probabilità d’errore per confronto (\(\alpha_C\)), secondo la formula seguente:

\[\alpha_E = 1 - (1 - \alpha_C)^n\]

dove \(n\) è il numero dei confronti. Pertanto, possiamo correggere il P-level con la formula anzidetta (metodo di Sidak). Prendiamo quindi la sesta colonna del dataframe ‘confronti’, quella che contiene i P-levels non corretti e la trasformiamo:

alphaC <- as.data.frame(confronti)[,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, 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:

\[\alpha_E = \alpha_C \cdot k\]

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

O meglio:

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

Sono possibili 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’.

Invece che aggiustare il P-level con uno dei metodi indicati più sopra è possibile considerare che, nel caso di contrasti e/o confronti, ogni singolo test d’ipotesi consiste in un rapporto tra una stima e il suo errore standard. Se l’ipotesi nulla di non-significatività è vera, la sampling distribution è la distribuzione t di Student univariata. Di conseguenza, l’intera famiglia di confronti/contrasti segue la distribuzione di t multivariato, con una matrice di correlazione che è deducibile dal contesto, come indicato da Bretz et al., (2011), pag. 73.

In altre parole, noto che sia il valore di t di ogni contrasto/confronto, posso desumere la relativa probabilità dalla distribuzione di t multivariata, invece che da quella univariata. Ovviamente il calcolo manuale è complesso e dovremo affidarci al software, come abbiamo già esemplificato più sopra:

#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’.

Questo tipo di correzione, così come le altre, tiene conto del numero di confronti che vengono effettuati. Se ad esempio, volessimo eseguire un set di confronti tipo ‘dunnett’, cioè tutti verso uno, avremmo 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

Nell’ottica esposta nel precedente sotto-capitolo, che prevede l’uso preferenziale degli intervalli di confidenza al posto del test d’ipotesi, è molto più interessante creare degli intervalli di confidenza familywise. Nel caso più semplice dei confronti a coppie nell’ANOVA per disegni ortogonali (bilanciati), si può utilizzare al posto del valore \(t_{\alpha/2, \nu}\) il valore ottenuto dalla distribuzione t multivariata, che, per nostra fortuna, si può facilmente desumere dalle tabelle dello ‘Studentised Range’, in funzione del numero di trattamenti in prova. Ad esempio, si può consultare questo link, da dove desumiamo che lo Studentised Range per 4 medie e 12 gradi di libertà dell’errore è 4.1985. Di conseguenza, se consideriamo ancora il secondo dei sei confronti a coppie illustrati in precedenza (Metribuzin__348 vs. Rimsulfuron_30; SE = 2.7702), l’intervallo di confidenza non corretto sarebbe:

limSup <- -7.6850 + qt(0.975, 12) * 2.770209
limInf <- -7.6850 + qt(0.025, 12) * 2.770209
limInf; limSup
## [1] -13.72077
## [1] -1.649233

Mentre l’intervallo di confidenza corretto sarebbe :

limSup <- -7.6850 + 1/sqrt(2) * 4.1985 * 2.770209
limInf <- -7.6850 - 1/sqrt(2) * 4.1985 * 2.770209
limInf; limSup
## [1] -15.90916
## [1] 0.5391627

Possiamo osservare che lo Studentised Range viene diviso per \(\sqrt{2}\). Con R possiamo ottenere lo stesso risultato:

confint(contrast(medie, method="pairwise"))
##  contrast                         estimate   SE df lower.CL upper.CL
##  Metribuzin__348 - Mixture_378        4.05 2.77 12    -4.18   12.272
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12   -15.91    0.539
##  Metribuzin__348 - Unweeded         -17.60 2.77 12   -25.82   -9.373
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12   -19.96   -3.508
##  Mixture_378 - Unweeded             -21.64 2.77 12   -29.87  -13.421
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12   -18.14   -1.688
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 4 estimates

Nel caso dei confronti tutti contro uno (tipo Dunnet), l’intervallo di confidenza può essere analogamente calcolato con la distribuzione t-multivariato. Le tabelle da consultare in questo caso sono diverse, perché, a parità di numero di medie, il numero di confronti è inferiore. Tuttavia, preferiamo utilizzare R:

confint(contrast(medie, method="dunnett"))
##  contrast                         estimate   SE df lower.CL upper.CL
##  Mixture_378 - Metribuzin__348       -4.05 2.77 12  -11.541     3.45
##  Rimsulfuron_30 - Metribuzin__348     7.68 2.77 12    0.192    15.18
##  Unweeded - Metribuzin__348          17.60 2.77 12   10.104    25.09
## 
## Confidence level used: 0.95 
## Conf-level adjustment: dunnettx method for 3 estimates

Sono possibili altre procedure di correzione più avanzate (Shaffer, Westfall), che tuttavia sono valide in presenza di alcune assunzioni aggiuntive e debbono quindi essere valutate con attenzione.

Per chiudere questo trattazione dei contrasti/confronti, vogliamo puntualizzare qualcosa che non sarà sfuggito ai lettori più esperti, cioè la differenza con le procedure di confronto multiplo tradizionale, basate sulla definizione di una differenza critica.

Esistono diversi metodi per calcolare una differenza critica; il più diffuso (Minima Differenza significativa di Fisher o Least Significant Difference, abbreviato MDS o LSD) impiega la distribuzione t di Student:

\[ MDS = SED \times t_{\alpha/2; \nu}\]

Dove \(\alpha\) è il tasso di errore per confronto e \(\nu\) è il numero di gradi di libertà del residuo. Per dataset bilanciati, l’uso della MDS fornisce esattamente gli stessi risultati della funzione ‘contrast()’ utilizzata più sopra, con l’argomento ‘adjust = “none”’.

Un’altra differenza critica molto utilizzata è la Honest Significant Difference di Tukey (per i confronti a coppie), che utilizza, invece della distribuzione t univariata, la distribuzione t multivariata. Questa procedura fornisce risultati analoghi a quelli della funzione ‘contrast()’, con aggiustamento per la molteplicità.

La differenza critica di Dunnett consente invece di confrontare tutte le medie con un testimone (o con il migliore/peggiore dei trattamenti), fornendo risultati analoghi a quelli ottenuti con la funzione ‘contrasts(method = “dunnet”)’.

Esistono almeno altre tre procedure classiche di confronto multiplo, che elenchiamo di seguito:

  1. Test di Duncan;
  2. Test di Newman-Keuls;
  3. Test di confronto multiplo di Tukey.

In genere queste procedure sono sconsigliabili, per i seguenti motivi:

  1. sono basate su differenze critiche multiple (crescenti al crescere della distanza dei trattamenti in graduatoria) e quindi non consentono la definizione di un’intervallo di confidenza. Di conseguenza, tra le domande ‘biologiche’ alle quali si cerca la risposta con i confronti multipli (si veda all’inizio) sono in grado di rispondere solo alla prima e non alla seconda e alla terza (non consentono il calcolo di un intervallo di confidenza).
  2. Non danno protezione ne’ per un tasso di errore per confronto ne’ per esperimento, ma rimangono a metà strada, in modo imprecisato (quindi il p level non è effettivamente noto, né a livello di singolo confronto né a livello di intero esperimento.

8.8.3 Altre letture

Il riferimento definitivo è:

Hsu, J., 1996. Multiple comparisons. Theory and methods. Chapman and Hall.

Altre letture:

  1. Bretz, F., T. Hothorn, and P. Westfall. 2002, On Multiple Comparison Procedures in R. R News 2: 14-17.
  2. Calinsky, T. and L. C. A. Corsten. 1985, Clustering means in ANOVA by simultaneus testing. Biometrics, 41, 39-48.
  3. 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.
  4. 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.
  5. Carmer, S. G. 1976, Optimal significance levels for application of the least significant difference in crop performance trials. Crop Science, 16, 95-99.
  6. Chew, V. 1976, Comparing treatment means: a compendium. Hortscience, 11(4), 348-357.
  7. Cousens, R. 1988, Misinterpretetion of results in weed research through inappropriate use of statistics. Weed Research, 28, 281-289.
  8. Edwards, A. W. F. and L. L. Cavalli-Sforza. 1965, A method for cluster analysis. Biometrics, 21, 362-375.
  9. Gates, C. E. and J. D. Bilbro. 1978, Illustration of a cluster analysis method for mean separation. Agronomy journal, 70, 462-465.
  10. O’Neill, R. and G. B. Wetherill. 1971, The present state of multiple comparison methods. Journ. Roy. Stat. Soc. , 2, 218-249.
  11. Petersen, R. G. 1977, Use and misuse of multiple comparison procedurs. Agronomy Journal, 69, 205-208.
  12. Scott, A. J. and M. Knott. 1974, A cluster analysis method for grouping means in the analysis of variance. Biometrics, 30, 507-512.
  13. Willavize, S. A., S. G. Carmer, and W. M. Walker. 1980, Evaluation of cluster analysis for comparing treatment means. Agronomy journal, 72,317-320.

  1. All we know about the world teaches us that the effects of A and B are always different-in some decimal place-for any A and B. Thus asking “are the effects different?” is foolish.