sabato 22 agosto 2009

Calcolo della tabella ANOVA partendo da dati aggregati

Può capitare di dover eseguire un test per verificare l'uguaglianze delle medie tra più gruppi, di cui però non disponiamo dei dati grezzi (il valore di ogni singolo dato), ma solamente della media, della deviazione standard e della numerosità campionaria di ciascun gruppo. Questo capita frequentemente nel momento in cui vogliamo fare ulteriori analisi sui dati pubblicati in articoli scientifici, che ci forniscono i dati descrittivi dei gruppi, e non i gruppi in sè.
Si possono ricavare facilmente le formule per ottenere la tabella ANOVA partendo dai dati descrittivi (in questo pdf sono presenti le formule per eseguire la one-way ANOVA, e le ANOVA bi- e tridimensionale partendo dai dati descrittivi). Supponiamo di avere i seguenti 4 gruppi:


a <- c(7.4, 6.6, 6.7, 6.1, 6.5, 7.2)
b <- c(7.1, 7.3, 6.8, 6.9, 7)
c <- c(6.8, 6.3, 6.4, 6.7, 6.5, 6.8)
d <- c(6.4, 6.9, 7.6, 6.8, 7.3)


Da questi dati calcoliamo: numerosità campionaria, media e deviazione standard, e raccogliamo tutti i dati in una tabella di riepilogo:


> mean(a)
[1] 6.75
> mean(b)
[1] 7.02
> mean(c)
[1] 6.583333
> mean(d)
[1] 7
>
> length(a)
[1] 6
> length(b)
[1] 5
> length(c)
[1] 6
> length(d)
[1] 5
>
> sd(a)
[1] 0.4764452
> sd(b)
[1] 0.1923538
> sd(c)
[1] 0.2136976
> sd(d)
[1] 0.4636809




Adesso dimentichiamo i dati campionari, e facciamo una ANOVA parametrica sui dati di quest'ultima tabella. Le formule sono:



La GrandMean è la media di tutti i gruppi; la SSB è la devianza tra gruppi, e la SSE è la devianza entro gruppi (o d'errore). Le altre formule sono quelle consuete, che già conoscevamo:



I df sono i gradi di liberta; MSB e MSE sono la varianza tra gruppi e d'errore; F è il valore della statistica test.

Traduciamo questi calcoli in R, e procediamo:


nA <- 6; nB <- 5; nC <- 6; nD <- 5
meanA <- 6.75; meanB <- 7.02; meanC <- 6.58; meanD <- 7
sdA <- .476; sdB <- .192; sdC <- .214; sdD <- .464

GrandMean <- (meanA * nA + meanB * nB + meanC * nC + meanD * nD) / (nA + nB + nC + nD)
> GrandMean
[1] 6.821818

SSB <- (((meanA - GrandMean)^2)*nA) + (((meanB - GrandMean)^2)*nB) + (((meanC - GrandMean)^2)*nC) + (((meanD - GrandMean)^2)*nD)
> SSB
[1] 0.7369273

SSE <- (sdA^2 * (nA-1)) + (sdB^2 * (nB-1)) + (sdC^2 * (nC-1)) + (sdD^2 * (nD-1))
> SSE
[1] 2.3705

dfSSB <- 4 - 1
dfSSE <- (nA+nB+nC+nD) - 4

MSB <- SSB / dfSSB
> MSB
[1] 0.2456424

MSE <- SSE / dfSSE
> MSE
[1] 0.1316944

F <- MSB / MSE
> F
[1] 1.865245


Possiamo calcolare il p-value di una distribuzione F, per prendere la nostra decisione statistica, in questo modo:


p.value <- 1 - pf(F, dfSSB, dfSSE)
> p.value
[1] 0.1716690


Raccogliamo tutti i dati nella classica tabella ANOVA:



Essendo p-value > 0.05, accettiamo l'ipotesi nulla: i 4 gruppi hanno medie uguali.

Adesso vediamo i valori che avremmo ottenuto, se avessimo utilizzato i dati campionari:


a <- c(7.4, 6.6, 6.7, 6.1, 6.5, 7.2)
b <- c(7.1, 7.3, 6.8, 6.9, 7)
c <- c(6.8, 6.3, 6.4, 6.7, 6.5, 6.8)
d <- c(6.4, 6.9, 7.6, 6.8, 7.3)

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

anova(lm(dati ~ gruppi))

Analysis of Variance Table

Response: dati
Df Sum Sq Mean Sq F value Pr(>F)
gruppi 3 0.72730 0.24243 1.8402 0.1760
Residuals 18 2.37133 0.13174


Potete confrontare i valori contenuti nella tabella, dopo averla resa leggibile in Latex, con la funzione xtable, disponibile nel pacchetto xtable (converte le tabelle di R in codice LaTeX):

Confrontando i valori, essi risultano molto simili, approssimazioni a parte.

giovedì 20 agosto 2009

ANOVA post-hoc tests #2

In questo post continuiamo a vedere come utilizzare R per i test post-hoc di una ANOVA significativa. Raccomando la lettura del primo post in riguardo allo stesso argomento, 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

  • test S di Scheffé






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:



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



con = varianza entro gruppi.
Questo valore (che è lo stesso della procedura di Bonferroni) viene confrontato con il valore critico:



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.

mercoledì 19 agosto 2009

Regressione logistica ordinale

In vari post precedenti (tra cui questi: [1^][2^][3^]), abbiamo discusso su come eseguire un modello di regressione logistica qualosa la variabile dipendente (Y) sia di tipo dicotomico. Se la variabile dipendente può invece assumente più di 2 valori (ossia la variabile dipendente è policotomica), si ricorre alla regressione logistica ordinale, in inglese ordered logistic regression. Esistono numerosi approcci a questo problema, e i principali modelli che sono stati formulati sono: cumulative logit o proportional odds model (Walker and Duncan 1967; McCullagh 1980), il continuation ratio model (Feinberg 1980), il constrained and unconstrained partial proportional odds models (Peterson and Harrell 1990), ed il adjacent-category logistic model (Agresti 1984). In questo articolo vengono esaminati alcuni di questi: LINK.
Eseguire in R una ordered logistic regression è alquanto difficoltoso; mancano difatti pacchetti dedicati esclusivamente all'analisi di dati categoriali (nominali). Inoltre, a quanto ne so, l'unico modello che è stato implementato è il proportional-odds model, ed è questo che utilizzerò negli esempi che ora vedremo.

Utilizzerò il dataset disponibile in questo pdf, ottima dispensa che ho trovato in rete sulla regressione logistica.
La seguente tabella di contingenza presenta i dati raccolti da un sondaggio universitario sottoposto ad alcuni studenti. L'obiettivo che ci proponiamo è quello di realizzare un modello che spieghi la variabile dipendente Y ordinale "giudizio complessivo), in base a "chiarezza" (variabile indipendente X1) e "interesse" (variabile indipendente X2).
Disponiamo quindi di:
- variabile dipendente policotomica
- 2 variabili indipendenti policotomiche



Per prima cosa dobbiamo cercare un modo per importare i dati di questa tabella in R. Io utilizzo il seguente schema: creo un file.txt sostituendo ai descrittori delle variabili ("decisamente no"..."decisamente si"), i numeri, da 1 a 4; si creano 3 colonne (X1, X2, Y) più una quarta colonna contenente le frequenze osservate. Il contenuto del file è il seguente:


chiarezza interesse giudizio frequenza
1 1 1 10
1 1 2 2
1 1 3 3
1 1 4 0
1 2 1 6
1 2 2 2
1 2 3 0
1 2 4 0
1 3 1 4
1 3 2 10
1 3 3 1
1 3 4 0
1 4 1 2
1 4 2 1
1 4 3 3
1 4 4 1
2 1 1 10
2 1 2 6
2 1 3 3
2 1 4 2
2 2 1 3
2 2 2 6
2 2 3 4
2 2 4 1
2 3 1 4
2 3 2 32
2 3 3 31
2 3 4 1
2 4 1 6
2 4 2 18
2 4 3 20
2 4 4 6
3 1 1 10
3 1 2 10
3 1 3 21
3 1 4 2
3 2 1 1
3 2 2 27
3 2 3 73
3 2 4 4
3 3 1 1
3 3 2 41
3 3 3 375
3 3 4 47
3 4 1 3
3 4 2 12
3 4 3 195
3 4 4 138
4 1 1 2
4 1 2 2
4 1 3 19
4 1 4 10
4 2 1 1
4 2 2 12
4 2 3 48
4 2 4 26
4 3 1 0
4 3 2 14
4 3 3 288
4 3 4 240
4 4 1 1
4 4 2 7
4 4 3 182
4 4 4 1224


Ogni elemento della riga deve essere spaziato (potremmo anche utilizzare la virgola). Salvato il file dataset.txt, importiamo i dati in R con il seguente comando:


dati <- read.table("/datasetR/dataset.txt", header=T)


All'interno di dati si trova ora la tabella, ripresa dal file txt di cui specifichiamo l'indirizzo (attenzione ad usare lo slash /, anzichè il backslash \, altrimenti R comunica un errore). Specificando header=TRUE, informiamo R che la prima riga del file da leggere è costituita dal nome delle colonne.
In che formato sono le colonne della tabella? Verifichiamo con il seguente comando:


sapply(dati, class)

chiarezza interesse giudizio frequenza
integer" "integer" "integer" "integer"


Le 4 colonne sono interpretate come numeriche. Eppure solo la colonna frequenza è di tipo numerico; le altre devono essere trasformate in fattori. Per far questo utilizziamo la funcione factor(); possiamo anche specificare il nome di ciascun elemento dei fattori (non è obbligatorio). Procediamo in questo modo e osserviamo la trasformazione, e verifichiamo con la funzione sapply():


dati$giudizio <- factor(dati$giudizio, labels=c("fantastico","buono","discreto","pessimo"))
dati$chiarezza <- factor(dati$chiarezza, labels=c("ottima","buona","discreta","pessima"))
dati$interesse <- factor(dati$interesse, labels=c("moltissimo","molto","poco","pochissimo"))

> dati
chiarezza interesse giudizio frequenza
1 ottima moltissimo fantastico 10
2 ottima moltissimo buono 2
3 ottima moltissimo discreto 3
4 ottima moltissimo pessimo 0
5 ottima molto fantastico 6
6 ottima molto buono 2
7 ottima molto discreto 0
8 ottima molto pessimo 0
9 ottima poco fantastico 4
10 ottima poco buono 10
11 ottima poco discreto 1
12 ottima poco pessimo 0
13 ottima pochissimo fantastico 2
14 ottima pochissimo buono 1
15 ottima pochissimo discreto 3
16 ottima pochissimo pessimo 1
17 buona moltissimo fantastico 10
18 buona moltissimo buono 6
19 buona moltissimo discreto 3
20 buona moltissimo pessimo 2
21 buona molto fantastico 3
22 buona molto buono 6
23 buona molto discreto 4
24 buona molto pessimo 1
25 buona poco fantastico 4
26 buona poco buono 32
27 buona poco discreto 31
28 buona poco pessimo 1
29 buona pochissimo fantastico 6
30 buona pochissimo buono 18
31 buona pochissimo discreto 20
32 buona pochissimo pessimo 6
33 discreta moltissimo fantastico 10
34 discreta moltissimo buono 10
35 discreta moltissimo discreto 21
36 discreta moltissimo pessimo 2
37 discreta molto fantastico 1
38 discreta molto buono 27
39 discreta molto discreto 73
40 discreta molto pessimo 4
41 discreta poco fantastico 1
42 discreta poco buono 41
43 discreta poco discreto 375
44 discreta poco pessimo 47
45 discreta pochissimo fantastico 3
46 discreta pochissimo buono 12
47 discreta pochissimo discreto 195
48 discreta pochissimo pessimo 138
49 pessima moltissimo fantastico 2
50 pessima moltissimo buono 2
51 pessima moltissimo discreto 19
52 pessima moltissimo pessimo 10
53 pessima molto fantastico 1
54 pessima molto buono 12
55 pessima molto discreto 48
56 pessima molto pessimo 26
57 pessima poco fantastico 0
58 pessima poco buono 14
59 pessima poco discreto 288
60 pessima poco pessimo 240
61 pessima pochissimo fantastico 1
62 pessima pochissimo buono 7
63 pessima pochissimo discreto 182
64 pessima pochissimo pessimo 1224


> sapply(dati, class)
chiarezza interesse giudizio frequenza
"factor" "factor" "factor" "integer"


Ora il formato dei dati è quello opportuno per proseguire. In realtà, se il nostro file txt contenesse direttamente il nome, e non il numero, per ogni categoria di ogni variabile, R avrebbe interpretato direttamente come un fattore le varie colonne.

Per eseguire una regressione (logistica) ordinale con il proportional-odds model, installiamo e carichiamo la libreria MASS, e utilizziamo la funzione polr():


library(MASS)
model <- polr(giudizio ~ chiarezza + interesse, data=dati, weights=frequenza, Hess=T)
summary(model)

Call:
polr(formula = giudizio ~ chiarezza + interesse, data = dati,
weights = frequenza, Hess = T)

Coefficients:
Value Std. Error t value
chiarezzabuona 1.1720 0.3509 3.340
chiarezzadiscreta 3.3901 0.3340 10.151
chiarezzapessima 5.3735 0.3404 15.784
interessemolto 0.5239 0.2546 2.058
interessepoco 1.4395 0.2240 6.426
interessepochissimo 3.1973 0.2296 13.924

Intercepts:
Value Std. Error t value
fantastico|buono 0.873 0.345 2.528
buono|discreto 3.041 0.369 8.229
discreto|pessimo 6.862 0.390 17.576

Residual Deviance: 4277.66
AIC: 4295.66


Osserviamo l'input di polr(): viene prima dichiarata la formula da analizzare (il giudizio in funzione di chiarezza e interesse), quindi si specifica la sorgente dei dati (data=dati), le frequenza (weights), e con la dicitura Hess=T abilitiamo la funzione sotto-utilizzata summary() sul modello (se non specifichiamo Hess=TRUE, ci viene restituito errore).

In output invece vengono forniti i valori dei parametri beta [coefficients] (relativi alle due variabili indipendenti) e dei parametri tau [intercepts] (rispettivi delle 4 modalità della variabile dipendente). Per quanto riguarda l'interpretazione dei coefficienti, rimando alla già citata dispensa sulla regressione logistica.

Concentrandoci invece sul programma R, sono disponibili altre funzioni per la ordered logit:

- Nella library Design è disponibile la funzione lrm().
- Nella library VGAM è disponibile la funzione vglm().
- Nella library Zelig è disponibile la funzione zelig().
- Nella library nnet è disponibile la funzione multinom().

Per ognuna di queste funzioni è disponibile l'help relativo (basta caricare la library e quindi digitare help(nome_funzione)). Tutte dovrebbero eseguire delle regressioni logistiche ordinali, ma gli output di tutte le funzioni, sugli stessi dati, sono completamente differenti. L'eccessiva confusione che si è creata (numerosi modelli + diverse modalità di input dati + diversi output), rendono l'argomento alquanto nebuloso.

lunedì 17 agosto 2009

Regressione lineare multivariata

Accolgo volentieri l'invito di Fabio, e mi accingo a cominciare alcuni post sulla statistica multivariata.
Nella regressione lineare semplice, abbiamo immaginato che una certa variabile Y dipendesse dall'andamento di un'altra variabile (X), in maniera lineare con andamento crescente o decrescente. Abbiamo quindi visto come realizzare e disegnare la retta che pone in relazione le due variabili, e come valutare la bontà del modello.
Nella realtà (scientifica, economica, psicometrica, etc.), quasi mai un evento dipende solamente dall'andamento di un certo fattore. Tutti gli eventi (anche i più comuni) sono influenzati da numerosissimi elementi. Risulta pertanto molto più utile formulare un modello che tenga conto di tutte queste influenze (o, come vedremo in seguito, delle influenze maggiori). Ciò si ottiene con lo studio della regressione lineare multipla (o multivariata). In generale si indica con Y la variabile dipendente, e con X seguito da un numero in pedice le variabili indipendenti che si suppone abbiano un effetto. Le X vengono chiamate predittori e la formula generale del modello che cerchiamo è:



è l'intercetta, mentre eccetera, sono i regressori, i quali rappresentano il coefficiente angolare della retta che otterremmo al variare del predittore corrispondente, qualora tutti gli altri predittori fossero costanti.
La rappresentazione grafica dipende dal numero di predittori che si vogliono considerare. Con un solo predittore, si ottiene (come abbiamo visto) una retta; con 2 predittori si ottiene un piano nello spazio a 3 dimensioni; con 3 o più predittori non è possibile la rappresentazione grafica, in quanto occorrerebbe rappresentare uno spazio a più di 3 dimensioni.

Cominciamo a vedere i passaggi da effettuare quando si vuole realizzare un modello di regressione multipla. Utilizzerò qui il dataset disponibile su questo sito. In esso viene raccolto il consumo di gelati pro-capite, considerando il periodo dell'anno, il prezzo del gelato, i guadagni delle famiglie e la temperatura media ambientale. Potremmo ad esempio chiederci: il consumo di gelato aumenta o diminuisce quando la temperatura aumenta? Sul consumo di gelati influisce maggiormente la temperatura o il suo prezzo?
Importiamo i dati. Sul sito è disponibile un file in formato .dat. Possiamo aprirlo con il blocco note per visualizzarlo, e renderci conto di come è strutturato. In R è facile importare questo formato di dati, anche senza la necessità di salvare il file sul nostro computer, utilizzando questo comando:


dati <- read.table("http://www.sci.usq.edu.au/staff/dunn/Datasets/Books/Hand/Hand-R/icecream-R.dat", header = TRUE, row.names = 1)


Adesso in dati è contenuta una tabella, con i valori riportati nel file specificato. Perchè R possa accedere al file .dat, non ci deve essere alcuna limitazione di connessione; se non riesce ad accedere ad internet, conviene controllare il proprio firewall (che spesso blocca la connessione a internet da parte di alcuni programmi). Per comodità trasformiamo la tabella in un dataframe:


datiDF <- data.frame(dati)


Disponiamo quindi di 4 variabili e 30 osservazioni. Per prima cosa cerchiamo di indagare il consumo di gelato in base alla temperatura e al prezzo (ignoriamo le entrate). Le colonne che ci interessano, per ora, sono datiDF$Consumption, datiDF$Price, datiDF$Temp. Cerchiamo un modello di regressione multipla a 3 predittori. La tecnica per stimare i regressori è detta ordinary least square OLS, che ripercorre gli stessi principi della tecnica dei minimi quadrati della regressione semplice.
In R stimiamo il modello in questo modo:


model <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
summary(model)

Call:
lm(formula = datiDF$Consumption ~ datiDF$Price + datiDF$Temp)

Residuals:
Min 1Q Median 3Q Max
-0.08226 -0.02051 0.00184 0.02272 0.10076

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.59655 0.25831 2.309 0.0288 *
datiDF$Price -1.40176 0.92509 -1.515 0.1413
datiDF$Temp 0.00303 0.00047 6.448 6.56e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04132 on 27 degrees of freedom
Multiple R-squared: 0.6328, Adjusted R-squared: 0.6056
F-statistic: 23.27 on 2 and 27 DF, p-value: 1.336e-06


Il model contiene la formula che vogliamo indagare. Dapprima si specifica la variabile dipendente (il consumo), quindi si specificano le variabili indipendenti, tra loro unite da un più.
Con la funzione summary() otteniamo una serie di informazioni. Innanzitutto il valore dell'intercetta e dei regressori; in questo caso:



Pertanto il modello di regressione multipla stimato è:



Osserviamo il segno del primo regressore, che descrive il predittore "prezzo": esso è negativo, il che significa che per ogni aumento unitario del prezzo del gelato, il consumo procapite diminuisce di circa una unità e mezzo (1.40); viceversa la temperatura ha un effetto positivo (seppure molto scarso) sul consumo di gelato (all'aumentare della temperatura, aumenta anche il consumo).

Possiamo rappresentare graficamente il modello creato. Installiamo il pacchetto scatterplot3d, e utilizziamo l'unica funzione in esso compreso, specificando le variabili nell'ordine x, y, z (nel nostro caso meglio specificare (X1, X2, Y)):


library(scatterplot3d)
sc <- scatterplot3d(datiDF$Price, datiDF$Temp, datiDF$Consumption, pch=16)




Possiamo anche intersecare il piano che abbiamo trovato, in questo modo:


model <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
sc$plane(model, col="red")




Analizziamo adesso il modello. Il coefficiente di determinazione (Multiple R-squared) consente di valutare la bontà del modello, ossia la proporzione di variabilità della Y spiegata dalle variabili esplicative considerate. In questo caso R-squared è pari a 0.63 (il 63% dei dati viene spiegato dalle variabili esplicative). Il valore R-squared adjusted è una correzione del coefficiente di determinazione che tiene conto del numero di predittori utilizzati (secondo alcuni statistici è più preciso quest'ultimo valore per definire la bontà del modello o goodness-of-fit).
Inoltre è opportuno effettuare una verifica dell'intercetta e dei regressori, ossia un t-test con ipotesi nulla H0: beta = 0. Ossia si deve verificare se il valore campionario possa essere esteso all'intera popolazione. Nel nostro caso i p-value associati alle variabili sono:
p-value = 0.0288 (significativo)
p-value = 0.1413 (non significativo - dovremmo ignorarla)
p-value = 6.56e-07 (significativo)

La statistica F ci indica se il modello è da scartare nella sua interessa, oppure se può essere ritenuto valido. Essendo F-statistic maggiore dell'F-tabulato, rifiutiamo l'ipotesi nulla che il modello sia da scartare nella sua interezza.

Possiamo calcolare in questo modo gli intervalli di confidenza dell'intercetta e dei 2 regressori:


confint(model)
2.5 % 97.5 %
(Intercept) 0.066549244 1.126551076
datiDF$Price -3.299896373 0.496375437
datiDF$Temp 0.002066035 0.003994569


Se vogliamo invece calcolare l'intervallo di confidenza di un valore predetto (ossia se vogliamo prevedere il valore che assumerà Y, dati i predittori), utilizziamo la funzione predict():


model <- lm(Consumption ~ Price + Temp, data=datiDF)
predict(model, data.frame(Price=0.280, Temp=30), interval="confidence")

fit lwr upr
1 0.2949663 0.2700105 0.3199220





Proviamo a complicare un pò le cose, e supponiamo di considerare nel nostro modello anche la variabile "entrate". All'aumentare dei guadagni, aumenta anche il consumo di gelato?


model <- lm(Consumption ~ Price + Temp + Income, data=datiDF)
summary(model)

Call:
lm(formula = Consumption ~ Price + Temp + Income, data = datiDF)

Residuals:
Min 1Q Median 3Q Max
-0.065302 -0.011873 0.002737 0.015953 0.078986

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1973151 0.2702162 0.730 0.47179
Price -1.0444140 0.8343573 -1.252 0.22180
Temp 0.0034584 0.0004455 7.762 3.1e-08 ***
Income 0.0033078 0.0011714 2.824 0.00899 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03683 on 26 degrees of freedom
Multiple R-squared: 0.719, Adjusted R-squared: 0.6866
F-statistic: 22.17 on 3 and 26 DF, p-value: 2.451e-07


Il regressore beta3 associato al predittore entrate risulta essere significativo, e di segno positivo. Quindi all'aumentare delle entrate, aumenta anche il consumo di gelati.
Il modello a 3 predittori viene complessivamente accettato (la statistica F risulta significativa), e si assiste ad un miglioramento della goodness-of-fit (è aumentato il valore dell'R quadro aggiustato).

Una domanda che potremmo porci è la seguente: tra le tre variabili considerate, qual è quella che maggiormente influisce sulla variabile consumo? Per far questo possiamo confrontare i valori dei regressori beta, ma solo dopo averli standardizzati (altrimenti non è possibile il confronto). Nel pacchetto QuantPsyc è presenta la funzione lm.beta() che fa al caso nostro:


library(QuantPsyc)
lm.beta(model)

Price Temp Income
-0.1324351 0.8632558 0.3140085


Quelli qui riportati sono i coefficienti beta standardizzati. Osserviamo che la variabile esplicativa che maggiormente influisce sul consumo di gelati è la temperatura.
Se non volete scaricare il pacchetto per utilizzare solo questa funzione, è sufficiente sapere che il calcolo di coefficienti beta standardizzati è pari al prodotto del valore beta per il rapporto tra la deviazione standard della X considerata, e la deviazione standard della Y:


coef(model)[2]*(sd(datiDF$Price))/(sd(datiDF$Consumption))
Price
-0.1324351

coef(model)[3]*(sd(datiDF$Temp))/(sd(datiDF$Consumption))
Temp
0.8632558

coef(model)[4]*(sd(datiDF$Income))/(sd(datiDF$Consumption))
Income
0.3140085


In maniera analoga (ma anche un pò più rischiosa), si potrebbe valutare la dipendenza anche osservando la correlazione tra le variabili. Se calcoliamo il coefficiente di correlazione (con la funzione cor())tra il consumo e le variabili prezzo, entrate, temperatura, osserviamo che il maggior valore è relativo alla variabile temperatura (è maggiore la correlazione tra queste due variabili).

Spesso si finisce col cercare una quantità di possibili variabili esplicative, da rendere poco agevole lo sviluppo di un modello di regressione multipla, che teoricamente dovrebbe contenere il minor numero di predittori possibile. La scelta non si basa solo sulla valutazione "ad occhio" del peso di ciascun predittore (come abbiamo fatto prima), ma viene utilizzato il test F parziale, che in pratica va a testare se c'è un incremento dell'R-quadro nel momento in cui si aggiunge un'altra variabile esplicativa al modello.
In R per eseguire un test F-parziale, si confronta la tabella dell'analisi della varianza dei due modelli. Supponiamo di voler confrontare le bontà dei modelli (Price + Temp) e (Price + Temp + Income):


model1 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
model2 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp + datiDF$Income)
anova(model1, model2)

Analysis of Variance Table

Model 1: datiDF$Consumption ~ datiDF$Price + datiDF$Temp
Model 2: datiDF$Consumption ~ datiDF$Price + datiDF$Temp + datiDF$Income
Res.Df RSS Df Sum of Sq F Pr(>F)
1 27 0.046090
2 26 0.035273 1 0.010817 7.9734 0.008989 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Se rifiutiamo l'ipotesi H0, vuol dire che l'inserimento della variabile esplicativa ha apportato un miglioramento al modello. In questo caso F = 7.97, con p-value < 0.05. Rifiutiamo quindi l'ipotesi nulla: l'aggiunta del predittore entrate migliora il modello.

Possiamo infine calcolare i coefficiente di determinazioni parziali. Fissata una delle variabili esplicative, a quanto ammonta la bontà del modello considerando l'altra variabile?
Riprendiamo a considerare solo 2 predittori (Temp e Price). Costruiamo prima i due modelli lineari semplici, quindi il modello a 2 predittori. Tenendo presente le formule per calcolare il coefficiente di determinazione parziale, abbiamo:


lm.x1 <- lm(datiDF$Consumption ~ datiDF$Price)
lm.x2 <- lm(datiDF$Consumption ~ datiDF$Temp)
lm.x1.x2 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)

ssr.x1.x2 <- sum(anova(lm.x1.x2)$"Sum Sq"[1:2])
ssr.x1 <- anova(lm.x1)$"Sum Sq"[1]
ssr.x2 <- anova(lm.x2)$"Sum Sq"[1]
ssr.x1.x2 - ssr.x1
sse.x1 <- anova(lm.x1)$"Sum Sq"[2]
sse.x2 <- anova(lm.x2)$"Sum Sq"[2]

(ssr.x1.x2 - ssr.x2) / sse.x2
[1] 0.07837315

(ssr.x1.x2 - ssr.x1) / sse.x1
[1] 0.6062858


Il coefficiente di determinazione parziale corrispondente alla variabile X1, quando X2 è tenuta costante, ci dice che per una data temperatura, il 7.83% della variabilità del consumo di gelati può essere spiegata dalla variabilità del prezzo nei negozi. Il coefficiente di determinazione parziale corrispondente alla variabile X2, quando X1 è tenuta costante, ci dice che per un dato prezzo, il 60.63% della variabilità del consumo di gelati può essere spiegata dalla variabilità della temperatura (sono valori in accordo con il peso definito con il calcolo dei regressori standardizzati).


Queste viste finora solo alcune delle analisi preliminari che dobbiamo effettuare quando si cerca un modello di regressione lineare multipla. Nei post successivi vedremo altri esempi, e andremo ad effettuare altre analisi, come l'analisi dei residui, la scelta dei predittori e il confronto tra più modelli.

mercoledì 12 agosto 2009

Operazioni con l'odd-ratio in R

Il software R non mette a disposizione molte funzioni per lavorare con le tabelle di contingenza 2x2, e quindi con gli odd-ratio. Tuttavia vedremo in questo post come calcolare l'OR, come stimare l'intervallo di confidenza dell'OR con i metodi di Wolff e Miettinen, come verificare la significatività dell'OR, e come confrontare tra loro due odd-ratios.

Consideriamo la seguente tabelle 2x2:



Immettiamo questi dati in R:


dati <- c(235, 103, 525, 635)
tabella <- matrix(dati, nrow=2)
colnames(tabella) <- c("malSI", "malNO")
rownames(tabella) <- c("InquiSI", "inquiNO")
table <- as.table(tabella)
table

malSI malNO
InquiSI 235 525
inquiNO 103 635


Abbiamo creato così la tabella a doppia entrata. La formula per ottenere l'OR è questa:



In R possiamo calcolarlo, partendo dalla tabella, in questo modo:


OR <- (table[1,1]*table[2,2]) / (table[2,1]*table[1,2])
OR
[1] 2.759593


In R è quindi abbastanza facile richiamare i valori contenuti in una tabella, semplicemente specificando il nome della tabella, e le coordinate (numero riga, numero colonna) della cella da utilizzare.




Verifichiamo adesso l'ipotesi della significatività dell'OR appena calcolato. Si può procedere con uno Z-test, oppure in termini di chi-quadro corretto.

Cominciamo a vedere la verifica di ipotesi dell'OR con lo Z-test.

L'ipotesi nulla è che OR = 1, il che significa, in termini di logaritmo naturale, che H0: ln(OR) = 0.
La statistica test è:



L'errore standard di ln(OR) è:



Quindi la statistica test in R diventa:


esOR <- sqrt(sum( 1/dati ))
z <- log(OR) / esOR
z

[1] 7.685698


Poichè Z è maggiore del valore critico 1.96, rifiuto l'ipotesi nulla: l'odd ratio è significativo.

Adesso verifichiamo l'OR in termini di chi-quadro corretto. L'ipotesi nulla è H0: OR=1, e la statistica test è:



E(a) è il valore atteso della cella a, calcolato con la formula:
var(a) è la varianza della cella a, calcolata con questa formula:
Portiamo le formule in R e calcoliamo il chi-quadro corretto:


E.a <- (table[1,1]+table[1,2])*(table[1,1]+table[2,1]) / sum(dati)
E.a
[1] 171.4820

var.a <- (table[1,1]+table[1,2])*(table[1,1]+table[2,1])*
(table[2,1]+table[2,2])*(table[1,2]+table[2,2]) /
(sum(dati)*sum(dati)*(sum(dati)-1))
var.a
[1] 65.4635

chi.square <- ((abs(table[1,1]-E.a)-0.5)^2) / var.a
chi.square
[1] 60.6639


Anche con il test del chi-quadro corretto, viene rifiutata l'ipotesi nulla: l'OR è significativo.




Calcoliamo ora l'intervallo di confidenza dell'odd-ratio con il metodo di Wolff. Per far questo ci occorre calcolare l'errore standard del logaritmo naturale dell'OR. Semplicemente questo valore è pari alla somma dei reciproci delle celle, sotto radice quadrata. Quindi in R diventa:


esOR <- sqrt(sum( 1/dati ))


L'intervallo di confidenza con il metodo di Wolff (alpha = 0.05), è:


ciMin <- log(OR) - (1.96 * esOR)
> ciMin
[1] 0.7562176

ciMax <- log(OR) + (1.96 * esOR)
> ciMax
[1] 1.273949

ORmin <- exp(ciMin)
ORmax <- exp(ciMax)

> ORmin
[1] 2.130204
> ORmax
[1] 3.574942


Quindi l'odd-ratio risulta essere compreso tra 2.13 e 3.57.

Possiamo calcolare inoltre l'intervallo di confidenza dell'OR con il metodo di Miettinen:



Il chi quadro che qui compare è il chi-quadro corretto già calcolato; Z è il calore critico per alpha stabilito. In R i due limiti dell'intervallo di confidenza sono:


espMin <- 1 - 1.96/(sqrt(chi.square))
espMax <- 1 + 1.96/(sqrt(chi.square))

mORmin <- OR^espMin
mORmax <- OR^espMax

> mORmin
[1] 2.137509

> mORmax
[1] 3.562724


Con il metodo di Miettinen, l'OR risulta compreso tra 2.14 e 3.56, valori molti vicini all'intervallo di confidenza calcolato con il metodo di Wolff.




Possiamo calcolare l'OR anche partendo dalla regressione logistica semplice, OR è pari all'esponenziale del regressore beta1:


table

malSI malNO
InquiSI 235 525
inquiNO 103 635

dft <- as.data.frame(table)
dft

Var1 Var2 Freq
1 InquiSI malSI 235
2 inquiNO malSI 103
3 InquiSI malNO 525
4 inquiNO malNO 635

fit <- glm(Var2 ~ Var1, weights = Freq, data = dft, family = binomial(logit))
OR <- exp(fit$coefficient[2])
OR

Var1inquiNO
2.759593


Per calcolare l'intervallo di confidenza con il metodo di Wolff, si utilizza la funzione confint():


ci <- confint(fit)

2.5 % 97.5 %
(Intercept) 0.6514110 0.9592309
Var1inquiNO 0.7588576 1.2770034


Calcoliamo l'esponenziale dei valori ottenuti in questa tabella, e otteniamo l'intervallo di confidenza:


> exp(ci[2,1])
[1] 2.135835

> exp(ci[2,2])
[1] 3.585878


L'intervallo di confidenza è (2.13, 3.58). Possiamo specificare anche il metodo con cui calcolare l'intervallo di confidenza: confint(fit, method="miettinen"), oppure confint(fit, method="wolff").




Il test Z utilizzato prima per verificare la significatività dell'OR può anche essere utilizzato per confrontare tra loro gli OR di due tabelle. Gli OR vengono calcolati indipendentemente; l'ipotesi nulla è OR1 = OR2, e la statistica test di confronto è la seguente:



La varianza dell'odd ratio è pari alla somma dei reciproci delle celle.
Supponiamo di avere le seguenti due generiche tabelle, di cui si vogliono confrontare gli odd-ratios:



Immettiamo i dati in R, calcoliamo gli OR, quindi le loro varianze, e infine il valore Z:


dati1 <- c(47, 21, 62, 77)
tabella1 <- matrix(dati1, nrow=2)
colnames(tabella1) <- c("SI", "NO")
rownames(tabella1) <- c("A", "B")
table1 <- as.table(tabella1)

dati2 <- c(54, 21, 95, 61)
tabella2 <- matrix(dati2, nrow=2)
colnames(tabella2) <- c("SI", "NO")
rownames(tabella2) <- c("A", "B")
table2 <- as.table(tabella2)

dft1 <- as.data.frame(table1)
dft2 <- as.data.frame(table2)

fit1 <- glm(Var2 ~ Var1, weights = Freq, data = dft1, family = binomial(logit))
OR1 <- exp(fit1$coefficient[2])

fit2 <- glm(Var2 ~ Var1, weights = Freq, data = dft2, family = binomial(logit))
OR2 <- exp(fit2$coefficient[2])

varOR1 <- sum(1/dati1)
varOR2 <- sum(1/dati2)

#ricordo che ln(OR) = coefficiente beta1 della regressione logistica
num <- fit1$coefficient[2] - fit2$coefficient[2]
num <- OR1-OR2
den <- sqrt(varOR1 + varOR2)
Z.OR <- num / den

pnorm(-abs(Z.OR))
[1] 0.004917582


Il p-value è 0.005, inferiore all'alpha prefissato di 0.05. Quindi rifiutiamo l'ipotesi nulla: gli OR delle due tabelle sono significativamente diversi.

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é

venerdì 7 agosto 2009

Test di analisi di omogeneità della varianza

Una delle applicazioni più classiche di questi test, è l'analisi della varianza con il metodo ANOVA. Per poter procedere con una ANOVA omoschedastica, è necessario fare l'assunzione che le varianze dei k gruppi in analisi siano omogenee (altrimenti occorre procedere con una ANOVA eteroschedastica). I test per verificare l'omogeneità delle varianze sono numerosi; alcuni vengono applicati su campioni che seguono una distribuzione normale, altri vengono utilizzati per l'analisi non parametrica (campioni non gaussiani). Nel post sull'ANOVA parametrica omoschedastica, abbiamo già visto 2 test di verifica dell'omogeneità delle varianze. Vediamoli nuovamente in questo post assieme ad altri test (sono tutti test parametrici, che pertanto possono essere applicati se i dati seguono una distribuzione normale; la verifica di normalità è discussa in questo post):

  • Test di Bartlett

  • Test di Fligner-Killeen

  • Brown e Forsyth test

  • Test di Hartley

  • Metodo di Cochran

  • Test di Levene parametrico


Un test per la verifica dell'omogeneità delle varianze, non parametrico (quando i gruppi non seguono una distribuzione normale), è il test di Levene non parametrico, che sfrutta il test di Kruskal-Wallis (ne parlo in fondo al post).

Innanzitutto vediamo come presentare i dati. In generale si deve creare un vettore contenente tutti i dati, e una seconda variabile di tipo factor, per mantenere la suddivisione in gruppi dei dati.
Consideriamo ad esempio i seguenti 4 gruppi, costituiti ciascuno da 5 elementi, e creiamo i vettori dati e gruppi:


a <- c(191, 211, 204, 209, 207)
b <- c(137, 148, 127, 137, 153)
c <- c(172, 165, 186, 178, 189)
d <- c(199, 206, 233, 183, 194)

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


Se tutti i gruppi hanno la stessa numerosità campionaria, è sufficiente modificare il parametro each, dichiarando la lunghezza dei gruppi. Se i gruppi sono invece di numerosità campionarie differenti, la procedura è leggermente diversa: ne parlerò più giù, a proposito del test di Levene.

Può essere utile creare un dataframe (alcuni test accettano anche questa classe di dati):


data <- data.frame(dati, gruppi)
data

dati gruppi
1 191 a
2 211 a
3 204 a
4 209 a
5 207 a
6 137 b
7 148 b
8 127 b
9 137 b
10 153 b
11 172 c
12 165 c
13 186 c
14 178 c
15 189 c
16 199 d
17 206 d
18 233 d
19 183 d
20 194 d





Test di Bartlett

Questo è sicuramente uno dei più utilizzati, specie in biostatica. E' disponibile nel pacchetto standard (stats) di R, e quindi può essere richiamato velocemente:


# potremmo anche digitare: bartlett.test(dati ~ gruppi, data=data)
> bartlett.test(dati, gruppi)

Bartlett test of homogeneity of variances

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


In questo caso, otteniamo p-value > 0.05. Quindi accettiamo (o per meglio dire, non rifiutiamo) l'ipotesi nulla: le varianze sono omogenee.




Test di Fligner-Killeen

Anche questo test è disponibile nel pacchetto stats. Riceve i dati con le stesse modalità della funzione per il test di Bartlett:


fligner.test(dati ~ gruppi, data=data)

Fligner-Killeen test of homogeneity of variances

data: dati by gruppi
Fligner-Killeen:med chi-squared = 1.3866, df = 3, p-value = 0.7087


Anche qui abbiamo ottenuto p-value > 0.05, quindi le varianze risultano essere omogenee.




Test di Brown e Forsyth

Per poter utilizzare questo test occorre scaricare, installare e quindi aprire la libreria HH. Quindi si procede come seguen:


library(HH)
hov(dati ~ gruppi, data=data)

hov: Brown-Forsyth

data: dati
F = 0.6509, df:gruppi = 3, df:Residuals = 16, p-value = 0.5939
alternative hypothesis: variances are not identical


I dati possono essere presentati solamente sottoforma di una formula (come quelle da fittare per una regressione). Il p-value risulta maggiore di 0.05, quindi le varianze sono omogenee.




Test di Hartley

Non ho trovato questo test implementato in nessun pacchetto. Si procede in questo modo: si calcolano le varianze di tutti i gruppi e quindi si calcola il rapporto (varianza massima) / (varianza minima):


varianze <- c(var(a), var(b), var(c), var(d))
varianze
[1] 62.8 104.8 97.5 351.5

max(varianze) / min(varianze)
[1] 5.597134


Il valore calcolato va confrontato con le tavole tabulate della distribuzione F-max (che non è la ben più nota distribuzione F di Fisher-Snedocor), essendo noti il numero di gruppi k e i gradi di libertà di ciascun gruppo con stessa numerosità campionaria (n-1).
Potete trovare una tavola tabulata a questo link.

Andiamo a cercare in orizzontale il numero k di varianze a confronto (4) e in verticale di gradi di libertà (4), e osserviamo che il valore tabulato è 20.6 per alpha = 0.05.
Poichè il valore calcolato è inferiore al valore tabulato, accettiamo l'ipotesi nulla: le varianze sono omogenee.




Metodo di Cochran

Questo test è disponibile nel pacchetto outliers. Dopo averlo scaricato e installato, possiamo applicare la funzione cochran.test:


library(outliers)
cochran.test(dati ~ gruppi)

Cochran test for outlying variance

data: dati ~ gruppi
C = 0.5701, df = 5, k = 4, p-value = 0.1117
alternative hypothesis: Group d has outlying variance
sample estimates:
a b c d
62.8 104.8 97.5 351.5


Il p-value è maggiore di 0.05, quindi anche secondo questo test le varianze sono omogenee. Se non vogliamo scaricare il pacchetto, possiamo facilmente calcolare il valore della statistica test, da confrontare con le tavole tabulate della distribuzione per le varianze di Cochran.
Se sospettiamo che è presente una varianza più grande delle altre, si calcola il rapporto (varianza massima) / (somma di tutte le varianze). In alternativa si calcola il rapporto (varianza minima) / (somma di tutte le varianze). In questo caso è presente una varianza più grande delle altre, quindi calcoliamo il rapporto massimo:


varianze <- c(var(a), var(b), var(c), var(d))
varianze
[1] 62.8 104.8 97.5 351.5

max(varianze) / sum(varianze)
[1] 0.5700616


Questo valore va confrontato con le tavole di Cochran, disponibili in questo pdf.
Il valore tabulato per 4 gruppi, e 4 gradi di libertà, con alpha = 0.05, è 0.6287, che è maggiore del valore calcolato. Quindi accettiamo l'ipotesi nulla: le varianze sono omogenee.




Test di Levene, per campioni normali

Esistono diverse rielaborazioni di questo test, che si basa sull'analisi della varianza degli scarti dei valori campionari dalla media, o dalla mediana o dalla media-trimmed (gli ultimi due sono anche chiamati test di Levene modificati secondo Brown-Forsyth). A differenze degli altri test visti finora, il test di Levene può essere applicato anche quando i k gruppi hanno diverse numerosità campionarie.
Esaminiamo ad esempio i seguenti 4 gruppi (sono gli stessi di prima, cui è stato tolto l'ultimo elemento dai gruppi a e c). Per prima cosa creiamo il vettore contenente tutti i dati, e il fattore gruppi, nel modo seguente:


a <- c(191, 211, 204, 209)
b <- c(137, 148, 127, 137, 153)
c <- c(172, 165, 186, 178)
d <- c(199, 206, 233, 183, 194)

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


A questo punto dobbiamo scaricare la libreria lawstat e applicare la funzione levene.test, scegliendo quale dei 3 test applicare. Personalmente utilizzo il metodo classico basato sugli scarti dalla media, ma possiamo cambiare questa scelta, modificando il parametro location (che può assumere i valori mean, median e trim.mean):


library(lawstat)
levene.test(dati, gruppi, location="mean")

classical Levene's test based on the absolute deviations from the mean ( none not applied
because the location is not set to median )

data: dati
Test Statistic = 0.8531, p-value = 0.4879


Il p-value è ancora una volta maggiore di 0.05, quindi le varianze risultano omogenee.
Se volessimo calcolare manualmente il valore della statistica test (0.8531), dovremmo procedere come segue. Dapprima si calcola la media dei 4 gruppi (oppure la mediana, o la media trimmed). Quindi si calcolano la differenze in valore assoluto (funzione abs) tra le misure e il valore calcolato; infine si avvia una ANOVA a una via, e il risultato della statistica F è il valore che cercavamo:


a <- c(191, 211, 204, 209)
b <- c(137, 148, 127, 137, 153)
c <- c(172, 165, 186, 178)
d <- c(199, 206, 233, 183, 194)

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

a.mean <- mean(a)
b.mean <- mean(b)
c.mean <- mean(c)
d.mean <- mean(d)

a.scarto <- abs(a - a.mean)
b.scarto <- abs(b - b.mean)
c.scarto <- abs(c - c.mean)
d.scarto <- abs(d - d.mean)

dati.scarto <- c(a.scarto, b.scarto, c.scarto, d.scarto)

anova(lm(dati.scarto ~ gruppi))

Analysis of Variance Table

Response: dati.scarto
Df Sum Sq Mean Sq F value Pr(>F)
gruppi 3 139.71 46.57 0.8531 0.4879
Residuals 14 764.26 54.59


La statistica F ha valore 0.8531, esattamente quello calcolato con il test di Levene; anche il p-value è lo stesso.
In alternativa, possiamo calcolare le differenze tra i valori e la mediana (median()), o tra i valori e la media trimmed.




Test di Levene non parametrico per campioni non normali

Fino abbiamo visto una lunga serie di test che può essere applicata quando tutti i gruppi seguono una distribuzione normale. Se però siamo interessati a verificare l'omogeneità delle varianze in gruppi che non seguono una disitribuzione non normale, dobbiamo applicare un test non parametrico di verifica della omogeneità delle varianze: è il test di Levene non parametrico.
Fondamentalmente la teoria è molto semplice: calcolati gli scarti dalla media dei valori rilevati, anzichè applicare una ANOVA parametrica ai dati ottenuti (vedi sopra il test di Levene parametrico), si applica l'analisi della varianza non-parametrica, ossia il test di Kruskal-Wallis.

Osserviamo il seguente esempio: le varianze di 4 gruppi, con numerosità campionarie diverse, vengono confrontate:


a <- c(2, 8, 8, 4)
b <- c(8, 7, 14)
c <- c(33, 59, 48, 56)
d <- c(60, 101, 67)

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

levene.test(dati, gruppi, location="mean", kruskal.test=T)

rank-based (Kruskal-Wallis) classical Levene's test based on the absolute deviations from the
mean ( none not applied because the location is not set to median )

data: dati
Test Statistic = 6.5705, p-value = 0.08692


Essendo il p-value > 0.05, accettiamo l'ipotesi di omogeneità delle varianze.

Alla luce delle considerazioni fatte sul test di Levene parametrico, osservate il codice seguente, che porta ad ottenere lo stesso risultato della statistica test:


a.mean <- mean(a)
b.mean <- mean(b)
c.mean <- mean(c)
d.mean <- mean(d)

a.scarto <- abs(a - a.mean)
b.scarto <- abs(b - b.mean)
c.scarto <- abs(c - c.mean)
d.scarto <- abs(d - d.mean)

dati = list(g1=a.scarto, g2=b.scarto, g3=c.scarto, g4=d.scarto)

kruskal.test(dati)

Kruskal-Wallis rank sum test

data: dati
Kruskal-Wallis chi-squared = 6.5705, df = 3, p-value = 0.08692

mercoledì 5 agosto 2009

I test di verifica di normalità

Numerose volte nei posts precedenti abbiamo proceduto con l'analisi assumendo che i dati seguissero una distribuzione normale (applicando quindi la statistica parametrica). Vediamo adesso come verificare che i propri dati seguano una distribuzione gaussiana. I test a disposizione sono molto numerosi, e la maggior parte di questi sono disponibili nei vari pacchetti di R. Vediamo qui di seguito l'applicazione di alcuni di essi, ed in particolare:


  • D'agostino's K-squared test

  • Jarque-Bera test

  • Anderson-Darling test for normality

  • Cramer-von Mises normality test

  • Lilliefors (Kolmogorov-Smirnov) test for normality

  • Shapiro-Francia normality test

  • Pearson chi-square normality test

  • Shapiro-Wilk normality test

  • Kolmogorov-Smirnov test



Per tutti gli esempi che analizzeremo, userò un set di dati ottenuto applicando la funzione rnorm(), la quale ci fornisce 50 numeri che seguono una distribuzione normale, nota la media (10) e la deviazione standard (2). Ci aspettiamo quindi che tutti i test che vedremo affermino che i dati seguono una distribuzione normale.

x <- rnorm(50, 10, 2)




D'agostino's K-squared test (D'Agostino-Pearson normality test)

Che io sappia questo test non è stato implementato in nessun pacchetto di R. Tuttavia, conoscendo le formule (disponibili su WikiPedia), è facilmente ricostruibile una funzione. Qui di seguito riporto il codice:


dagostino.pearson.test <- function(x) {
# from Zar (1999), implemented by Doug Scofield, scofield at bio.indiana.edu
DNAME <- deparse(substitute(x))
n <- length(x)
x2 <- x * x
x3 <- x * x2
x4 <- x * x3
# compute Z_g1
k3 <- ((n*sum(x3)) - (3*sum(x)*sum(x2)) + (2*(sum(x)^3)/n)) /
((n-1)*(n-2))
g1 <- k3 / sqrt(var(x)^3)
sqrtb1 <- ((n - 2)*g1) / sqrt(n*(n - 1))
A <- sqrtb1 * sqrt(((n + 1)*(n + 3)) / (6*(n - 2)))
B <- (3*(n*n + 27*n - 70)*(n+1)*(n+3)) / ((n-2)*(n+5)*(n+7)*(n+9))
C <- sqrt(2*(B - 1)) - 1
D <- sqrt(C)
E <- 1 / sqrt(log(D))
F <- A / sqrt(2/(C - 1))
Zg1 <- E * log(F + sqrt(F*F + 1))
# compute Z_g2
G <- (24*n*(n-2)*(n-3)) / (((n+1)^2)*(n+3)*(n+5))
k4 <- (((n*n*n + n*n)*sum(x4)) - (4*(n*n + n)*sum(x3)*sum(x)) -
(3*(n*n - n)*sum(x2)^2) + (12*n*sum(x2)*sum(x)^2) -
(6*sum(x)^4)) /(n*(n-1)*(n-2)*(n-3))
g2 <- k4 / var(x)^2
H <- ((n-2)*(n-3)*abs(g2)) / ((n+1)*(n-1)*sqrt(G))
J <- ((6*(n*n - 5*n + 2)) / ((n+7)*(n+9))) * sqrt((6*(n+3)*(n+5)) /(n*(n-2)*(n-3)))
K <- 6 + (8/J)*(2/J + sqrt(1 + 4/(J*J)))
L <- (1 - 2/K) / (1 + H*sqrt(2/(K-4)))
Zg2 <- (1 - 2/(9*K) - (L^(1/3))) / (sqrt(2/(9*K)))
K2 <- Zg1*Zg1 + Zg2*Zg2
pk2 <- pchisq(K2, 2, lower.tail=FALSE)
RVAL <- list(statistic = c(K2 = K2), p.value = pk2, method =
"D'Agostino-Pearson normality test\n\nK2 is distributed as Chi-squared
with df=2", alternative = "distribution is not normal", data.name =
DNAME)
class(RVAL) <- "htest"
return(RVAL)
}


L'ipotesi nulla è che i nostri dati seguono una distribuzione normale. Applichiamo la funzione ai nostri dati:


dagostino.pearson.test(x)

D'Agostino-Pearson normality test

K2 is distributed as Chi-squared with df=2

data: x
K2 = 1.2788, p-value = 0.5276
alternative hypothesis: distribution is not normal


Essendo p-value > 0.05, accettiamo l'ipotesi nulla: siamo di fronte a una distribuzione normale.




Jarque-Bera normality test

E' disponibile nel pacchetto tseries, che deve essere scaricato, ed installato.


library(tseries)
jarque.bera.test(x)

Jarque Bera Test

data: x
X-squared = 1.1181, df = 2, p-value = 0.5717


Anche per questo test, l'ipotesi nulla è la normalità del vettore x. Essendo p-value > 0.05, accettiamo l'ipotesi nulla.




Anderson-Darling test for normality

Test disponibile nella libreria nortest (libreria che raccoglie numerosi test di normalità).


library(nortest)
ad.test(x)

Anderson-Darling normality test

data: x
A = 0.1931, p-value = 0.8898


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0





Cramer-von Mises normality test


library(nortest)
cvm.test(x)

Cramer-von Mises normality test

data: x
W = 0.0289, p-value = 0.8576


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Lilliefors (Kolmogorov-Smirnov) test for normality


library(nortest)
lillie.test(x)

Lilliefors (Kolmogorov-Smirnov) normality test

data: x
D = 0.0808, p-value = 0.5715


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Shapiro-Francia normality test


library(nortest)
sf.test(x)

Shapiro-Francia normality test

data: x
W = 0.9918, p-value = 0.9454


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Pearson chi-square test for normality


library(nortest)
pearson.test(x)

Pearson chi-square normality test

data: x
P = 1.2, p-value = 0.991


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Shapiro-Wilk normality test

Questo è uno dei test più utilizzati, ed è già disponibile nel pacchetto di base di R (stats), quindi non necessita di alcuna installazione.


shapiro.test(x)

Shapiro-Wilk normality test

data: x
W = 0.9865, p-value = 0.834


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Kolmogorov–Smirnov test

Anche questo test è disponibile nel pacchetto di base. Può essere utilizzato per effettuare due diverse analisi. Possiamo ad esempio chiederci "i dati X e Y provengono dalla stessa distribuzione?". In questo caso la verifica si effettua in questo modo:


x <- rnorm(50)
y <- rnorm(30)

ks.test(x, y)

Two-sample Kolmogorov-Smirnov test

data: x and y
D = 0.1533, p-value = 0.7195
alternative hypothesis: two-sided


Essendo p-value > 0.05, accettiamo l'ipotesi nulla: i due vettori provengono dalla stessa distribuzione.

Un'altra verifica che possiamo effettuare è la seguente: "il vettore X proviene da una distribuzione normale, o da una distribuzione chi-quadro a 1 grado di libertà?". In R si procede così:


ks.test(x, "pnorm")

One-sample Kolmogorov-Smirnov test

data: x
D = 0.0931, p-value = 0.7436
alternative hypothesis: two-sided

ks.test(x, "pchisq", 1)

One-sample Kolmogorov-Smirnov test

data: x
D = 0.54, p-value = 3.564e-14
alternative hypothesis: two-sided


In questo caso, oltre a specificare il vettore da analizzare, si deve pure specificare la funzione di confronto, messa tra virgolette. Per testare altre distribuzioni, potete trovare la sintassi precisa digitando ??distribution: così leggete l'elenco delle distribuzioni implementate in R; quindi, richiamando l'help relativo alla distribuzione di vostro interesse (ad esempio help(Normal)), potete vedere la corretta sintassi da utilizzare (generalmente è il nome della funzione, preceduto dalla lettera p: pnorm, pgamma, pchisq).
Come potevamo prevedere, il test di normalità dà p-value > 0.05, mentre il confronto con la distribuzione chi-quadro dà p-value < 0.05. Quindi il vettore x segue una distribuzione normale.