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é
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.
RispondiEliminaFammi sapere!
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?
RispondiEliminanon 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è bilaterale, devi usare alfa mezzi
Eliminase 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?
RispondiEliminaScusate, mi potete aiutare su come inserire le lettere? grazie
RispondiEliminaGreat blog post.
RispondiElimina