ANOVA in R: aov(), la tabella ANOVA e il test di Tukey HSD

Perché l'ANOVA e non una pila di t-test

Un t-test confronta due medie. Con tre o più gruppi viene la tentazione di fare un t-test su ogni coppia, ed è proprio questa la trappola. Ogni test eseguito al livello 0.05 comporta un rischio di falso positivo del 5%, e i rischi si accumulano: con 4 gruppi sono 6 test a coppie e circa il 26% di probabilità di ottenere almeno un risultato "significativo" spurio; con 5 gruppi (10 test), circa il 40%. Finiresti per fabbricare scoperte dal puro rumore.

L'ANOVA (ANalysis Of VAriance, analisi della varianza) risolve il problema ponendo una domanda con un p-value: le medie di tutti i gruppi sono uguali, oppure almeno una differisce? Nonostante il nome, confronta delle medie: lo fa semplicemente analizzando la varianza. Se le medie dei gruppi sono più disperse di quanto il rumore interno ai gruppi possa spiegare, sta succedendo qualcosa di reale.

ANOVA a una via con aov()

PlantGrowth è fatto apposta: pesi di piante in una condizione di controllo e in due condizioni di trattamento. Stima il modello con aov() e poi, cosa importante, stampa la tabella con summary():

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

Leggi la formula come "weight dipende da group?". La variabile di raggruppamento deve essere un factor: PlantGrowth$group lo è già, ma se i tuoi gruppi sono codificati come numeri (dosi, ID di lotto), avvolgili: aov(y ~ factor(dose), ...). Altrimenti aov() stima in silenzio una retta di regressione attraverso i codici dei gruppi invece di confrontare le medie dei gruppi: modello sbagliato, nessun messaggio di errore.

Leggere la tabella ANOVA

La tabella ha due righe, il fattore e i residui, e cinque colonne. Cella per cella:

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df: i gradi di libertà. La riga del fattore riceve gruppi − 1 (3 gruppi → 2); la riga dei residui riceve osservazioni − gruppi (30 − 3 = 27). Un rapido controllo che R abbia visto il disegno sperimentale che intendevi.
• Sum Sq: la variabilità, divisa in due mucchi. La riga group è il mucchio tra i gruppi: quanto le medie dei gruppi si allontanano dalla media generale. Residuals è il mucchio entro i gruppi: quanto le piante variano attorno alla media del proprio gruppo. Insieme sommano la variabilità totale dei dati.
• Mean Sq: ogni Sum Sq diviso per i suoi Df, che trasforma i mucchi di variabilità in tassi confrontabili per grado di libertà. Il Mean Sq dei residui (0.389) è il livello di rumore.
• F value: il rapporto tra il Mean Sq di group e il Mean Sq di Residuals (1.8832 / 0.3886 ≈ 4.85). Se tutte le medie dei gruppi fossero davvero uguali, questo rapporto starebbe intorno a 1. Più F cresce, più le differenze tra i gruppi superano ciò che il rumore può spiegare.
• Pr(>F): il p-value, cioè la probabilità di ottenere un F così grande se le tre medie vere fossero identiche. Qui vale 0.016, abbastanza piccolo al livello convenzionale di 0.05 per rifiutare "tutte uguali". Come sempre, non è la probabilità che l'ipotesi nulla sia vera, e non dice nulla su quali gruppi differiscano né di quanto.

Quest'ultimo punto è il limite cruciale: un F significativo dice "c'è qualche differenza, da qualche parte". Niente di più.

Quali gruppi differiscono? TukeyHSD()

La domanda successiva richiede un test post-hoc. L'Honest Significant Difference di Tukey confronta ogni coppia mantenendo il tasso di errore a livello di famiglia al 5% su tutti i confronti insieme:

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

Una riga per coppia, quattro numeri per riga:

• diff: la differenza stimata tra le medie (il secondo gruppo nominato meno il primo).
• lwr, upr: l'intervallo di confidenza al 95% a livello di famiglia per quella differenza.
• p adj: il p-value, già corretto per il fatto di fare tre confronti.

Per PlantGrowth: trt2-trt1 mostra una differenza di circa 0.87 con p adj ≈ 0.012 e un intervallo lontano dallo zero, quindi il trattamento 2 fa crescere di più del trattamento 1. Entrambe le righe trattamento contro controllo (trt1-ctrl, trt2-ctrl) hanno intervalli che comprendono lo zero e p adj ben sopra 0.05: in questo campione nessuno dei due trattamenti si distingue dal controllo. La regola pratica è la stessa degli intervalli di confidenza in ogni contesto: l'intervallo esclude lo zero ⇔ p adj sotto 0.05.

Nota come l'F complessivo abbia detto "esiste una differenza" mentre Tukey la individua esattamente in una coppia. La struttura in due passi è tutto il senso del metodo: un test complessivo onesto, poi un lavoro da detective a coppie correttamente aggiustato.

ANOVA a due vie: due fattori e la loro interazione

Con due variabili di raggruppamento, una sola chiamata le testa entrambe e verifica anche se interagiscono. ToothGrowth incrocia il tipo di integratore (supp) con la dose (0.5, 1, 2 mg, numerica, quindi serve factor()):

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) si espande in tre effetti, ciascuno con una riga della tabella:

• supp: in media sulle dosi, il tipo di integratore conta? (Sì: p ≈ 0.0002.)
• factor(dose): in media sugli integratori, la dose conta? (Decisamente: p è minuscolo.)
• supp:factor(dose): l'interazione. L'effetto dell'integratore cambia a seconda della dose? Qui p ≈ 0.022, quindi sì. Il succo d'arancia batte l'acido ascorbico alle dosi basse, ma il divario si chiude alla dose di 2 mg.

Un'interazione significativa è un'etichetta di avvertimento sugli effetti principali: "il tipo di integratore conta" è vero solo in media, e la media nasconde una storia che dipende dalla dose. Quando l'interazione è significativa, descrivi le combinazioni (con un TukeyHSD(fit2, "supp:factor(dose)") o una tabella delle medie di gruppo tramite i riepiloghi per gruppo) invece di riportare solo gli effetti principali. Usa + al posto di * solo quando vuoi deliberatamente un modello senza interazione.

Le assunzioni e l'alternativa basata sui ranghi

La matematica dell'ANOVA si appoggia a tre assunzioni, in ordine decrescente di negoziabilità:

• Indipendenza: le osservazioni non si influenzano a vicenda. Niente rimedia a una violazione a posteriori; è una proprietà del disegno dello studio.
• Varianze uguali tra i gruppi: il Mean Sq dei residui è un'unica stima aggregata del rumore, quindi i gruppi dovrebbero essere più o meno ugualmente rumorosi. Controllo in una riga: bartlett.test(weight ~ group, data = PlantGrowth) (un p-value grande significa nessuna evidenza di varianze diverse).
• Residui approssimativamente normali: conta meno di tutto con disegni bilanciati e gruppi di dimensioni decenti; dai un'occhiata a un istogramma o a un boxplot dei residui tramite residuals(fit).

Quando i dati sono molto asimmetrici, ordinali o pieni di valori anomali, l'alternativa basata sui ranghi all'ANOVA a una via è kruskal.test(weight ~ group, data = PlantGrowth), il fratello a più gruppi di wilcox.test(). È una riga sola e scambia un po' di potenza con la robustezza.

Cosa ti porti a casa

• L'ANOVA pone una sola domanda sulle medie di 3 o più gruppi con un solo p-value: è il rimedio alla molteplicità dei t-test.
• fit <- aov(y ~ group, data = df); summary(fit), e la variabile di raggruppamento deve essere un factor.
• Il valore F è il segnale tra i gruppi diviso per il rumore entro i gruppi; un Pr(>F) piccolo significa "qualche differenza da qualche parte", e niente di più.
• TukeyHSD(fit) trova quali coppie differiscono, con la correzione per i confronti multipli già inclusa.
• a * b stima due fattori più la loro interazione; un'interazione significativa significa che gli effetti principali da soli non raccontano tutta la storia.
• Verifica l'uguaglianza delle varianze (bartlett.test) e la normalità dei residui; kruskal.test() è l'alternativa basata sui ranghi.

Prossimo passo: dal confronto tra medie di gruppo alla modellazione di una relazione, con la regressione lineare e lm().