domenica 9 agosto 2009

ANOVA post-hoc tests #1

In questo post vediamo come applicare alcuni test per indagare quali gruppi hanno influito sulla significatività di una analisi della varianza ANOVA. Ricordo che una ANOVA si dice significativa quando p-value < 0.05, il che significa che tra i gruppi che si sono considerati, almeno una coppia di medie è significativamente differente. Il nostro obiettivo è quello di ricercare quali coppie. I test a disposizione sono numerosi; in questo e nei prossimi post ne vedremo alcuni

Consideriamo i seguenti 4 gruppi, ripresi da questo link. Immettiamo in R i nostri dati, e creiamo i vettori che ci serviranno per le funzioni che applicheremo.
La prima cosa da fare, è stabilire se questi seguono una distribuzione normale, al fine di applicare l'analisi della varianza parametrica. In questo post sono raccolti alcuni dei test a disposizione (qui utilizzerò solo lo Shapiro-Wilk normality test).
Quindi si verifica l'omogeneità delle varianze, così da poter procedere con una ANOVA omoschedastica. In questo post sono raccolti alcuni test che possiamo applicare (qui utilizzo solo il test di Bartlett).



a <- c(13.47, 10.21, 15.10, 14.65, 9.03, 15.14, 6.09, 3.43, 5.95, 10.72, 10.01, 8.17)
b <- c(4.02, 14.03, 3.09, 10.25, 6.25, 1.00, 7.50, 10.18, 6.07, 2.03, 4.17, 7.28)
c <- c(10.56, 14.61, 11.88, 11.43, 9.73, 12.37, 3.82, 13.04, 13.28, 15.28, 12.90, 10.98)
d <- c(7.74, 2.77, 10.29, 4.03, 10.23, 7.67, 7.93, 6.75, 5.60, 8.19, 2.02, 7.88)

dati <- c(a,b,c,d)
gruppi = factor(rep(letters[1:4], each = 12))

shapiro.test(a)

Shapiro-Wilk normality test

data: a
W = 0.939, p-value = 0.4855

shapiro.test(b)

Shapiro-Wilk normality test

data: b
W = 0.9619, p-value = 0.8103

shapiro.test(c)

Shapiro-Wilk normality test

data: c
W = 0.8519, p-value = 0.0388

shapiro.test(d)

Shapiro-Wilk normality test

data: d
W = 0.916, p-value = 0.2548

bartlett.test(dati, gruppi)

Bartlett test of homogeneity of variances

data: dati and gruppi
Bartlett's K-squared = 2.1573, df = 3, p-value = 0.5404


Il vettore c risulta non-normale. Teoricamente quindi non potremmo procedere con una ANOVA parametrica, ma dovremmo utilizzare il test di Kruskal-Wallis. Tuttavia, procederò con l'ANOVA parametrica ai fini dell'esercizio. Le varianze risultano omogenee.
A questo punto fittiamo il modello, e calcoliamo la tabella ANOVA:


model <- lm(dati ~ gruppi)
anova(model)

Analysis of Variance Table

Response: dati
Df Sum Sq Mean Sq F value Pr(>F)
gruppi 3 243.67 81.22 7.1736 0.0005046 ***
Residuals 44 498.19 11.32
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Il p-value è minore di 0.05, quindi l'ANOVA è significativa: non tutte le medie sono significativamente uguali.
Procediamo quindi con i test post-hoc. Quelli che qui esamineremo sono:

  • test LSD di Fisher

  • test HSD di Tukey






Procedura LSD di Fisher

La sigla LSD sta per "least significance difference". Disponendo di k medie, ognuna calcolata su n dati (numerosità campionaria di ciascun gruppo), si calcola la quantità LSD, secondo la formula:



Il valore t di Student è quello relativo al proprio livello di confidenza (qui useremo sempre alpha = 0.05), e i gradi di liberà della varianza residua (N-k = 44). In R tale valore si calcola così:


qt(0.975, 44)
[1] 2.015368


Questo t si moltiplica per la radice quadrata del rapporto tra la varianza residua (il cui valore è presente della tabella ANOVA: 11.32), fratto la numerosità campionaria del singolo gruppo (12). Il risultato è:


qt(0.975, 44) * sqrt(2*11.32/12)
[1] 2.768


A questo punto si calcolano le differenze tra le medie di tutti i gruppi. Quindi la media del gruppo a meno la media del gruppo b; media-a meno media-c, e così via, prese in valore assoluto. Questo risultato si confronta con il valore LSD. Se la differenza è maggiore di LSD allora quella coppia è significativa (le medie sono significativamente differenti); se è minore, allora le medie sono uguali.
Manualmente si può procedere così:


> abs(mean(a) - mean(b))
[1] 3.841667

> abs(mean(a) - mean(c))
[1] 1.4925

#e così via


Con molti gruppi, il lavoro può essere laborioso. Ad ogni modo nella libreria agricolae (da scaricare e installare), la procedura LSD è a disposizione. La funzione LSD.test() riceve in input:
- i dati
- i gruppi
- i gradi di libertà della varianza residua (che possiamo calcolare in R così: df <- df.residual(model) )
- la varianza residua (che possiamo calcolare in R così: MSerror <- deviance(model)/df )

Applichiamo subito la funzione ai nostri dati:


library(agricolae)
model <- lm(dati ~ gruppi)
df<-df.residual(model)
MSerror<-deviance(model)/df
comparison <- LSD.test(dati, gruppi, df, MSerror, group=TRUE)

Study:

LSD t Test for dati
......
Alpha 0.050000
Error Degrees of Freedom 44.000000
Error Mean Square 11.322399
Critical Value of t 2.015368

Treatment Means
gruppi dati std.err replication
1 a 10.164167 1.1196506 12
2 b 6.322500 1.0982832 12
3 c 11.656667 0.8498728 12
4 d 6.758333 0.7694191 12

Least Significant Difference 2.768521
Means with the same letter are not significantly different.

Groups, Treatments and means
a c 11.65667
a a 10.16417
b d 6.758333
b b 6.3225


Specificando group=TRUE, possiamo leggere il valore LSD calcolato (2.768). Se vogliamo confrontare tutte le differenze con LSD, dobbiamo specificare group=FALSE.


comparison <- LSD.test(dati, gruppi, df, MSerror, group=FALSE)

Study:

LSD t Test for dati
......
Alpha 0.050000
Error Degrees of Freedom 44.000000
Error Mean Square 11.322399
Critical Value of t 2.015368

Treatment Means
gruppi dati std.err replication
1 a 10.164167 1.1196506 12
2 b 6.322500 1.0982832 12
3 c 11.656667 0.8498728 12
4 d 6.758333 0.7694191 12

Comparison between treatments means

tr.i tr.j diff pvalue
1 1 2 3.8416667 0.0076
2 1 3 1.4925000 0.2832
3 1 4 3.4058333 0.0171
4 2 3 5.3341667 0.0003
5 2 4 0.4358333 0.7525
6 3 4 4.8983333 0.0009


L'ultima parte dell'output è quello che ci serve in questo momento. Sono state calcolate tutte le differenze (primo gruppo meno secondo gruppo; primo gruppo meno terzo gruppo, e così via), ed è stato calcolato il relativo p-value di confronto con LSD.
- p-value > 0.05: differenza non significativa (le medie dei due gruppi sono uguali)
- p-value < 0.05: differenza significativa (le medie dei due gruppi non sono uguali)
E così si può ricavare quali sono le coppie che hanno reso significativa la ANOVA.

Il test LSD modificato da Winer propone di utilizzare la correzione di Bonferroni per calcolare la probabilità alpha del comparisonwise (in pratica alpha viene diviso per il doppio del numero di confronti da effettuare simultaneamente; quindi cambia il valore del t di Student). Se scegliamo di tener conto di questa modifica, in R dobbiamo semplicemente richiedere di aggiustare il p-value secondo il metodo di Bonferroni:


comparison <- LSD.test(dati, gruppi, df, MSerror, group=FALSE, p.adj="bonferroni")

Study:

LSD t Test for dati
P value adjustment method: bonferroni
......
Alpha 0.050000
Error Degrees of Freedom 44.000000
Error Mean Square 11.322399
Critical Value of t 2.762815

Treatment Means
gruppi dati std.err replication
1 a 10.164167 1.1196506 12
2 b 6.322500 1.0982832 12
3 c 11.656667 0.8498728 12
4 d 6.758333 0.7694191 12

Comparison between treatments means

tr.i tr.j diff pvalue
1 1 2 3.8416667 0.0458
2 1 3 1.4925000 1.0000
3 1 4 3.4058333 0.1024
4 2 3 5.3341667 0.0021
5 2 4 0.4358333 1.0000
6 3 4 4.8983333 0.0053


Con la correzione di Bonferroni risultano significativamente differenti le coppie 1-2, 2-3, 3-4.




Test HSD di Tukey

L'acronimo HSD sta per "honestly significant difference". I principi generali visti per il test LSD vengono applicati anche qui. L'unica differenze è il calcolo della quantità LSD, che qui viene indicata con HSD:



Il valore Q (Q-studentizzato), dipende dal livello di confidenza alpha (0.05), dal numero di medie a confronto (k) e dai gradi di libertà della varianza d'errore (N-k). Si ricerca quindi in tavole tabulate il valore, che in questo caso è pari a 3.789. Quindi la quantità HSD ha valore 3.668. Le differenze tra le medie di tutti i gruppi vanno confrontate con il valore HSD (come per il test LSD).
Sempre nel pacchetto agricolae è già disponibile la funzione che svolge tutti i calcoli. L'input è uguale alla funzione per il test LSD


comparison <- HSD.test(dati, gruppi, df, MSerror, group=T)

##cancello parte dell'output, perchè ci serve solo
##per osservare il valore HSD:

Honestly Significant Difference 3.667801

##...

comparison <- HSD.test(dati, gruppi, df, MSerror, group=FALSE)

Study:

HSD Test for dati
......
Alpha 0.050000
Error Degrees of Freedom 44.000000
Error Mean Square 11.322399
Critical Value of Studentized Range 3.775958

Treatment Means
gruppi dati std.err replication
1 a 10.164167 1.1196506 12
2 b 6.322500 1.0982832 12
3 c 11.656667 0.8498728 12
4 d 6.758333 0.7694191 12

Comparison between treatments means

tr.i tr.j diff pvalue
1 1 2 3.8416667 0.0369
2 1 3 1.4925000 0.6995
3 1 4 3.4058333 0.0773
4 2 3 5.3341667 0.0019
5 2 4 0.4358333 0.9888
6 3 4 4.8983333 0.0048


Anche in questo caso sono significative (ossia significativamente diverse) le differenze che in valore assoluto sono superiori al valore HSD, e che nell'output del programma R presentano p-value < 0.05. Con il test HSD di Tukey risultano significativamente differenti le medie dei gruppi 1-2, 2-3, 3-4 (le stesse coppie del test LSD con la modifica di Winer).

In questo post vengono trattati il test t di Bonferroni e il test S di Scheffé

6 commenti:

  1. Bel lavoro, grazie mille. Solo una cosa: non sono d' accordo col fatto che tu testi la normalitá sui dati stessi. Bisonga testare la normalitá dei residui ottenuti dal fit del modello lineare.
    Fammi sapere!

    RispondiElimina
  2. scusate ma come si calcolano come è scritto nel blog i gradi di liberà della varianza residua (N-k = 44) . mi spiego meglio cosa sono N e k nella formula?

    RispondiElimina
  3. non ho capito perchè, se l'alpha è 0,05, nella formula per il calcolo del valore critico di t viene usato 0,975 e non 0,95. Io avrei usato il primo se l'alpha fosse stato 0,025. Grazie.

    RispondiElimina
  4. se avessi un plot di dati con piu vetteori e dove ogni vettore "viaggia su piu giorni" esiste qualcosa che mi resituisca il p-value giornaliero tra i gruppi in cui voglio dividere i miei vettori?

    RispondiElimina
  5. Scusate, mi potete aiutare su come inserire le lettere? grazie

    RispondiElimina