In questo post continuiamo a vedere come utilizzare R per i test post-hoc di una ANOVA significativa. Raccomando la lettura del
, dove si è parlato del test LSD di Fisher, e del test HSD di Tukey. Qui utilizzerò gli stessi dati, per esaminare insieme i seguenti test:
Test t di Bonferroni
Il test per i confronti multipli simultanei proposto dallo statistico italiano Bonferroni, si basa sul confronto tra il valore t di Bonferroni, e una quantità che viene calcolata, di volta in volta, per ogni coppia di gruppi che si esamina.
Il t di Bonferroni si ricerca sulle tavole t di Student, in corrispondenza dei gradi di libertà della varianza d'errore (N-k), e del valore alpha ottenuto dividendo l'alpha di scelta (solitamente 0.05) per il doppio del numero di confronti da effettuare. Nel nostro caso abbiamo 4 gruppi; il numero di confronti da effettuare è quindi 6, che moltiplicato per 2 fa 12. In R tale t si calcola in questo modo. Supponiamo di avere una lista, con i k gruppi; il doppio del numero di confronti da effettuare è semplicemente k*(k-1). Quindi
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 <- list(g1=a, g2=b, g3=c, g4=d)
confronti2 <- length(dati) * (length(dati) - 1)
alpha.bonferr <- 0.05 / confronti2
abs(qt(alpha.bonferr, 44))
[1] 2.762815
Questo valore va confrontato con il risultato della seguente formula. Siano r ed s due generici gruppi in esame:
}}$$ "$$t_{bonferroni}=\frac{|\overline{x_r}-\overline{x_s}|}{\sqrt{S_e^2(\frac{1}{n_r}+\frac{1}{n_s})}}$$")
con 
= varianza entro gruppi.
Se il t di Bonferroni calcolato è maggiore di quello tabulato (nel nostro caso 2.76), allora la differenza è significativa (le medie sono significativamente differenti); altrimenti sono uguali.
Procedendo manualmente, quindi, si crea una matrice triangolare con i t di Bonferroni calcolati, e si individuano quelli maggiori (in valore assoluto) del t-tabulato, ossia le differenze significative (qui le indico con un asterisco). Con i nostri dati avremmo:
In R non è disponibile una funzione per eseguire direttamente il test di Bonferroni. Tuttavia sono riuscito a trovare in rete una bozza della funzione, che ho modificato per svolgere i calcoli visti finora. Ecco il codice, quindi, per eseguire il test di Bonferroni in R:
bonferroni.test = function(x, y)
{
# x is a list of vectors
# y is the alpha
num <- 0
denom <- 0
grandmean.num <- 0
sb2.num <- 0
nn <- numeric(length(x))
ss2 <- numeric(length(x))
for(i in 1:length(x)) {
nn[i] <- length(x[[i]])
ss2[i] <- (sum((x[[i]] - mean(x[[i]]))^2))/(nn[i] - 1
)
num <- ((nn[i] - 1) * ss2[i]) + num
denom <- nn[i] + denom
grandmean.num <- (nn[i] * mean(x[[i]])) +
grandmean.num
}
sw2 <- num/(denom - length(x))
grandmean <- grandmean.num/denom
for(j in 1:length(x)) {
sb2.num <- (nn[j] * ((mean(x[[j]]) - grandmean)^2)) +
sb2.num
}
sb2 <- sb2.num/(length(x) - 1)
FF <- sb2/sw2
p.val <- 1 - pf(FF, length(x) - 1, denom - length(x))
effort <- numeric(length(x) * (length(x) - 1))
effort1 <- character(length(x) * (length(x) - 1))
tt <- matrix(effort, ncol = length(x))
pp <- matrix(effort, ncol = length(x))
qq <- matrix(effort1, ncol = length(x))
alpha.star <- y/((length(x) * (length(x) - 1)))
deg.free <- denom - length(x)
for(k in 1:(length(x) - 1)) {
for(l in 2:length(x)) {
if(k >= l) {
tt[k, l] <- 0
}
else {
tt[k, l] <- abs((mean(x[[k]]) - mean(x[[l
]]))/sqrt(sw2 * ((1/nn[k]) + (1/nn[
l]))))
alpha <- 1 - (alpha.star)
t.tab <- qt(alpha, deg.free)
pp[k, k] <- 0
if(tt[k, l] == 0) {
qq[k, l] <- c("Zero")
}
else if(tt[k, l] < t.tab) {
qq[k, l] <- c("Not")
}
else if(tt[k, l] > t.tab) {
qq[k, l] <- c("Signif")
}
}
}
}
results <- list(grandmean, sw2, sb2, FF, p.val, tt,
alpha.star, qq)
results <- list(Grand.mean = grandmean, Within.var = sw2,
Between.var = sb2, F.value = FF, p.value = p.val,
bonferroni.t = tt, alpha.star = alpha.star,
bonferroni.sig = qq)
results
}
L'input è costituito dalla lista dei gruppi (x), e dal valore di alpha da cui vogliamo partire (qui userò alpha = 0.05). Applichiamo subito la funzione, e studiamo l'output:
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 <- list(g1=a, g2=b, g3=c, g4=d)
bonferroni.test(dati, 0.05)
$Grand.mean
[1] 8.725417
$Within.var
[1] 11.32240
$Between.var
[1] 81.22261
$F.value
[1] 7.173622
$p.value
[1] 0.0005046445
$bonferroni.t
[,1] [,2] [,3] [,4]
[1,] 0 2.796572 1.086478 2.4793041
[2,] 0 0.000000 3.883050 0.3172684
[3,] 0 0.000000 0.000000 3.5657816
$alpha.star
[1] 0.004166667
$bonferroni.sig
[,1] [,2] [,3] [,4]
[1,] "" "Signif" "Not" "Not"
[2,] "" "" "Signif" "Not"
[3,] "" "" "" "Signif"
Per comodità la funzione fornisce in output anche le varianze (tra gruppi ed entro gruppi), il risultato della statistica F di una ANOVA parametrica, e il suo corrispettivo p-value (per controllare se è utile eseguire i confronti multipli oppure no).
Quindi viene presentata la matrice triangolare con i valori t di Bonferroni calcolati come sopra. E' in pratica la stessa tabella che abbiamo scritto manualmente. Quindi viene riportato il valore alpha (0.05 / 12), ed infine una matrice triangolare contenente la significatività del test.
Risultano quindi significative (ossia le medie sono differenti) i confronti tra i gruppi 1-2, 2-3, 3-4. Osservate che il test HSD di Tukey (applicato in questo post) ha dato lo stesso risultato, ossia ha individuato come significative le stesse coppie.
Test di Scheffé per i confronti multipli
Il test implementato da Scheffé è molto simile al testo di Bonferroni. Per ogni coppia di dati si calcola la seguente quantità:
}}$$ "$$S=\frac{|\overline{x_r}-\overline{x_s}|}{\sqrt{S_e^2\cdot(\frac{1}{n_r}+\frac{1}{n_s})}}$$")
con = varianza entro gruppi.
Questo valore (che è lo stesso della procedura di Bonferroni) viene confrontato con il valore critico:
F_{(alpha, k-1, N-k)}}$$ "$$S_\alpha=\sqrt{(k-1)F_{(alpha, k-1, N-k)}}$$")
Se la quantità calcolata supera il valore critico, allora la differenza è significativa.
Il test è quindi molto simile al test di Bonferroni, eccetto il valore critico di confronto. Consideriamo il consueto dataset, e calcoliamo il valore critico manualmente:
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 <- list(g1=a, g2=b, g3=c, g4=d)
N <- length(a) + length(b) + length(c) + length(d)
alpha <- 1 - .05
F.cv <- sqrt( (length(dati)-1) * qf(alpha, (length(dati)-1), (N-length(dati))) )
> F.cv
[1] 2.906785
Proviamo ora a calcolare la quantità S per la coppia di gruppi a-b:
num <- abs( mean(a) - mean(b) )
dati.tot <- c(a,b,c,d)
gruppi = factor(rep(letters[1:4], each = 12))
within.var <- anova(lm(dati.tot~gruppi))[2,3]
den <- sqrt( within.var*(1/length(a)+1/length(b)) )
S <- num/den
> S
[1] 2.796572
Per la coppia a-b, S < S-critico, quindi la media di questi due gruppi è significativamente uguale.
Possiamo procedere in modo analogo per le altre coppie. Oppure utilizziamo la funzione sottostante, che è sostanzialmente uguale a quella del test di Bonferroni visto poco sopra, salvo le modifiche del test di Scheffé. Questo è il codice per eseguire il test di Scheffé in R:
scheffe.test = function(x, y)
{
# x is a list of vectors
# y is the alpha
num <- 0
denom <- 0
grandmean.num <- 0
sb2.num <- 0
nn <- numeric(length(x))
ss2 <- numeric(length(x))
for(i in 1:length(x)) {
nn[i] <- length(x[[i]])
ss2[i] <- (sum((x[[i]] - mean(x[[i]]))^2))/(nn[i] - 1
)
num <- ((nn[i] - 1) * ss2[i]) + num
denom <- nn[i] + denom
grandmean.num <- (nn[i] * mean(x[[i]])) +
grandmean.num
}
sw2 <- num/(denom - length(x))
grandmean <- grandmean.num/denom
for(j in 1:length(x)) {
sb2.num <- (nn[j] * ((mean(x[[j]]) - grandmean)^2)) +
sb2.num
}
sb2 <- sb2.num/(length(x) - 1)
FF <- sb2/sw2
p.val <- 1 - pf(FF, length(x) - 1, denom - length(x))
effort <- numeric(length(x) * (length(x) - 1))
effort1 <- character(length(x) * (length(x) - 1))
tt <- matrix(effort, ncol = length(x))
pp <- matrix(effort, ncol = length(x))
qq <- matrix(effort1, ncol = length(x))
alpha.star <- y/((length(x) * (length(x) - 1)))
deg.free <- denom - length(x)
for(k in 1:(length(x) - 1)) {
for(l in 2:length(x)) {
if(k >= l) {
tt[k, l] <- 0
}
else {
tt[k, l] <- abs((mean(x[[k]]) - mean(x[[l
]]))/sqrt(sw2 * ((1/nn[k]) + (1/nn[
l]))))
alpha <- 1 - y
f.tab <- sqrt( (length(x)-1) * qf(alpha, (length(x)-1), (denom-length(x))))
pp[k, k] <- 0
if(tt[k, l] == 0) {
qq[k, l] <- c("Zero")
}
else if(tt[k, l] < f.tab) {
qq[k, l] <- c("Not")
}
else if(tt[k, l] > f.tab) {
qq[k, l] <- c("Signif")
}
}
}
}
results <- list(grandmean, sw2, sb2, FF, p.val, tt,
alpha.star, qq)
results <- list(Within.var = sw2,
Between.var = sb2, F.value = FF, p.value = p.val,
Scheffé.S = tt, alpha = y, S.critical.value = f.tab,
Scheffé.sig = qq)
results
}
Adesso studiamo l'output della funzione:
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 <- list(g1=a, g2=b, g3=c, g4=d)
scheffe.test(dati, .05)
$Within.var
[1] 11.32240
$Between.var
[1] 81.22261
$F.value
[1] 7.173622
$p.value
[1] 0.0005046445
$Scheffé.S
[,1] [,2] [,3] [,4]
[1,] 0 2.796572 1.086478 2.4793041
[2,] 0 0.000000 3.883050 0.3172684
[3,] 0 0.000000 0.000000 3.5657816
$alpha
[1] 0.05
$F.critical.value
[1] 2.906785
$Scheffé.sig
[,1] [,2] [,3] [,4]
[1,] "" "Not" "Not" "Not"
[2,] "" "" "Signif" "Not"
[3,] "" "" "" "Signif"
Nell'ordine ci vengono forniti:
- la varianza entro gruppi
- la varianza tra gruppi
- il valore della statistica F dell'ANOVA e il suo p-value
- la quantità S calcolata per ogni coppia
- il valore alpha prescelto
- il valore critico di confronto
- la significativà (una coppia è significativa quando le medie sono differenti)
In base al test di Scheffé sono significative le coppie 2-3 e 3-4.
