Analisi della Varianza
Introduzione
Confronto fra variabili categoriali e a intervalli
Abbiamo visto, nel capitolo dedicato al T test, come il t-test ci permetta di confrontare le medie di due gruppi, e di valutare se la loro differenza è significativa.
Il confronto su coppie di elementi è applicabile se abbiamo una variabile indipendente di tipo categoriale con due soli livelli. Vi sono circostanze, però, in cui il confronto deve avvenire fra più di due gruppi. Il caso più semplice è quando la variabile indipendente ha tre o più livelli. Una circostanza più complessa emerge quando le variabili indipendenti sono più di due.
In quale modo possiamo affrontare una simile eventualità?
Confronto a coppie
Una possibile risposta è quella di confrontare ogni possibile combinazione di coppie di gruppi.
Nell'esempio semplice di una variabile indipendente a tre livelli, confrontare il gruppo 1 con il gruppo 2, il gruppo 1 con il 3, il 2 con il 3.
Sebbene questi confronti a coppie siano non solo possibili, ma realmente utilizzati nei confronti post-hoc, l'utilizzo esclusivo di questa metodologia incorre in due limiti.
Problemi del confronto a coppie
Il primo è che, moltiplicando il numero dei confronti, si accresce la probabilità di incorrere in un errore di tipo I.
Immaginiamo, ad esempio, di accettare un valore $\alpha$ pari a 0.05. Questo significa accettare la possibilità che, il 5% delle volte, si rifiuti indebitamente l'ipotesi nulla.
Ma se, invece di un solo confronto, ne facciamo due, qual è la probabilità che si rifiuti indebitamente almeno una volta l'ipotesi nulla?
Ricordiamo che, se l'ipotesi che si sta valutando è se la variabile indipendente influisce su quella dipendente, basta rilevare una differenza significativa fra i gruppi per inferire una influenza.
Ma se i confronti sono numerosi, la probabilità di incorrere in un errore del primo tipo aumenta. Nell'esempio con tre gruppi, poiché i confronti sono tre, con un $\alpha$ pari a 0.05 la probabilità di rifiutare erroneamente l'ipotesi nulla diventa pari a 0.1023.
Più variabili indipendenti
Il secondo inconveniente emerge quando le variabili indipendenti sono più di una. In questo caso, dal confronto a coppie è difficile far emergere quali variabili indipendenti hanno un'influenza significativa sulla variabile dipendente e quali no, così come diventa difficile far emergere l'eventuale interazione fra le variabili indipendenti.
Questi due inconvenienti ci costringono ad identificare una metodologia capace di generalizzare ai confronti con più variabili indipendenti e capace di mantenere sotto controllo l'errore di tipo I.
Varianze
Nel capitolo correlazione sulla regressione lineare sono stati introdotti i concetti di varianza totale, spiegata e residua.
La varianza di una variabile ci offre una misura della distribuzione di quella variabile. Conoscendo la media e la varianza (o la deviazione standard) di una variabile, ed assumendo una distribuzione di tipo normale, possiamo fare delle previsioni su Y. Possiamo assumere che $\bar$ sia il valore atteso di y più probabile, e possiamo stimare la probabilità dell'occorrenza di una osservazione y.
Naturalmente, la varianza di una variabile influisce sulla nostra capacità di fare delle previsioni. Se la varianza è prossima allo zero, noi possiamo prevedere con certezza che il valore atteso di una osservazione y sarà molto vicino alla media $\bar$
Varianza spiegata e previsioni
Riprendiamo l'esempio banale delle precipitazioni atmosferiche: la variazione delle osservazioni può essere molto ampia: da pochi millimetri a decine di centimetri. Tuttavia, tutta la variazione di peso può essere spiegata dal volume dell'acqua (attenzione: spiegata non significa causata).
Passiamo al caso meno banale di dover stimare il consumo di carburante durante l'uso di un'auto in una giornata, conoscendo i km percorsi.
In questo caso, come abbiamo visto, la scommessa sarà meno certa, in quanto il consumo è legato anche al tipo di guida, al tipo di percorso, alle condizioni di traffico e così via. Ciononostante, conoscere il numero di km percorsi mi permette una stima molto più accurata di quanto potrei fare semplicemente tirando ad indovinare.
La varianza residua, ovvero la varianza dei residui, sarà molto più bassa della varianza totale. Il numero di km costituisce dunque un predittore molto utile del consumo di carburante.
Varianza residua e previsioni
Il vantaggio di poter applicare la regressione lineare è che la relazione fra due variabili può essere espressa attraverso due soli parametri: $\beta_0$ e $\beta_1$, ovvero l'intercetta e la pendenza della linea.
Nel caso di variabili indipendenti di tipo categoriale, naturalmente, non è possibile assumere alcuna linearità, e dunque non possono bastarci quei due parametri.
Un esempio: gli affitti in una città
Immaginiamo di voler prendere in affitto un appartamento in una città di medie dimensioni, e vogliamo capire se, in differenti quartieri, i prezzi sono significativamente diversi.
Immaginiamo dunque di fare una ricerca sistematica, usando alcuni siti specializzati, raccogliendo le informazioni relative a 60 appartamenti, distribuiti su 3 quartieri: 20 nel quartiere A, 20 nel B, 20 nel C.
Quello che vogliamo capire è se, nei differenti quartieri, i prezzi sono significativamente diversi.
R: generare l'esempio
Con R, possiamo generare il dataset di 60 valori, invocando rnorm.
Decidiamo di generare 20 valori con media 550 e sd 30 (quartiere A), 20 con media 590 e sd 30 (quartiere B), 20 con media 620 e sd 30 (quartiere C).
Usando le funzioni factor, levels e data.frame creiamo il dataframe con 2 colonne (prezzo, quartiere) e 60 righe.
R: il codice
# creiamo il dataframe degli affitti delle 3 zone
zonaA <- round (rnorm (20,550,30))
zonaB <- round (rnorm (20,560,30))
zonaC <- round (rnorm (20,620,30))
zone <- c (zonaA,zonaB,zonaC)
fZone <- factor (c(rep('A',20),rep('B',20),rep('C',20)))
affitti <- data.frame (prezzo=zone,quartiere=fZone)
Grafico e varianza
Nel grafico fig:esempio1, rappresentiamo in un plot le 60 osservazioni: in rosso le venti del quartiere A, verde il quartiere B, blu il C. La linea rossa marca il valore medio di A, la verde la media di B, blu la media di C. Le due righe rosse tratteggiate, l'intervallo di confidenza al 95% delle osservazioni in A, le verdi in B, le blu in C. La riga continua nera rappresenta la media generale, le due righe tratteggiate nere l'intervallo di confidenza generale, al 95%.
plot(c(1,60),c(min(affitti$prezzo),max(affitti$prezzo)),
type = "n", xlab = "Zone", ylab = "Affitti",
main="Affitti per zona")
points(affitti$prezzo,col=as.integer(affitti$quartiere)+1)
lines (c(1,20),rep(mean(zonaA),2),col=2)
lines (c(21,40),rep(mean(zonaB),2),col=3)
lines (c(41,60),rep(mean(zonaC),2),col=4)
lines (c(1,60),rep(mean(affitti$prezzo),2),col=1)
lines (c(1,20),rep(mean(zonaA)+2*sd(zonaA),2),col=2,lty=2)
lines (c(1,20),rep(mean(zonaA)-2*sd(zonaA),2),col=2,lty=2)
lines (c(21,40),rep(mean(zonaB)+2*sd(zonaB),2),col=3,lty=2)
lines (c(21,40),rep(mean(zonaB)-2*sd(zonaB),2),col=3,lty=2)
lines (c(41,60),rep(mean(zonaC)+2*sd(zonaC),2),col=4,lty=2)
lines (c(41,60),rep(mean(zonaC)-2*sd(zonaC),2),col=4,lty=2)
lines (c(1,60),rep(mean(affitti$prezzo)+2*sd(affitti$prezzo),2),col=1,lty=2)
lines (c(1,60),rep(mean(affitti$prezzo)-2*sd(affitti$prezzo),2),col=1,lty=2)
Nonostante le variazioni dovute al caso, è chiaro che l'intervallo di confidenza dei singoli gruppi è minore (e diverso) dell'intervallo di confidenza totale. Questo significa che conoscere il quartiere dove un appartamento è collocato mi permette di fare delle previsioni migliori in merito al prezzo che mi aspetto di pagare.
Inferenza e previsioni
L'analisi bivariata (descrittiva e inferenziale) ci permette dunque innanzitutto di capire se una variabile influisce su di un'altra.
In secondo luogo, se l'influenza è statisticamente significativa, conoscere il valore della prima variabile ci permette di fare delle previsioni più accurate sulla seconda.
Rapporto fra varianza spiegata e residua
L'analisi della varianza è la statistica inferenziale che valuta se vi è una relazione fra una (o più) variabili indipendenti, di tipo categoriale, e una variabile quantitativa (almeno ad intervalli). Il principio su cui si basa la statistica è proprio la percentuale di varianza spiegata dal modello riespetto alla varianza totale.
La misura che viene presa in considerazione in questa statistica è dunque un rapporto: il rapporto fra varianza spiegata dal modello e varianza residua (ovvero la differenza fra la varianza totale e quella spiegata). Se il rapporto supera un determinato valore critico, si rifiuta l'ipotesi nulla (secondo cui non vi è relazione fra la variabile indipendente e quella dipendente).
L'analisi della Varianza
Quello che l'analisi della varianza ci permette di capire è se le medie della variabile dipendente osservate nei diversi gruppi sono o meno statisticamente diverse. Più precisamente, ci permette di stabilire se esistono almeno due gruppi la cui media sia statisticamente diversa [^1].
Il vantaggio di questo approccio, rispetto al confronto fra coppie di gruppi, è duplice:
- non vi è una proliferazione dell'errore di tipo I, in quanto il confronto è unico
- nel caso di più variabili indipendenti, è possibile stimare l'influenza di ognuna delle variabili indipendenti, nonché della loro interazione.
L'ipotesi nulla
L'ipotesi nulla assume che la media dei gruppi non sia fra loro significativamente diversa, e dunque che le medie dei gruppi siano approssimativamente pari alla media generale, salvo l'errore di campionamento.
Se le medie dei vari gruppi sono tutte perfettamente uguali alla media generale, anche la varianza dei gruppi sarà pari alla varianza generale, e dunque la varianza spiegata sarà pari a zero.
Errore di campionamento
A causa dell'errore di campionamento, però, sappiamo che, anche qualora l'ipotesi nulla sia vera, le medie dei gruppi potranno discostarsi dalla media generale, e dunque la varianza spiegata misurata sarà superiore a zero.
Come nei casi già visti (t test, correlazione, chi quadro), il valore del rapporto fra varianze va dunque confrontato con una distribuzione (teorica o generata empiricamente, ad esempio attraverso una simulazione) in modo da valutare se la proporzione di varianza spiegata è da attribuire al caso (errore di campionamento), e dunque va accettata l'ipotesi nulla, oppure no.
Il calcolo
Somme dei quadrati
Per calcolare il test dell'analisi della varianza, dobbiamo calcolare tre valori.
- la somma dei quadrati dell'errore totale, $SS_T$;
- la somma dei quadrati dell'errore residuo, $SS_R$;
- la somma dei quadrati del modello, $SS_M$;
Per calcolare le varianze totale, residua e spiegata dobbiamo dividere gli SS per i rispettivi gradi di libertà
Somma dei quadrati e varianza totale
La somma dei quadrati dell'errore totale si calcola con la formula $SS_T = \sum_^N (Y_n - \bar{Y_})^2$ dove N è il numero totale di osservazioni e $\bar{Y_}$ è la media totale.
I gradi di libertà della varianza totale sono $df_T = N - 1$.
La varianza totale è pari a $MS_T = SS_T / df_T$.
Somma dei quadrati e varianza residua
La somma dei quadrati dell'errore residuo si calcola con la formula $SS_R = \sum_^I \sum_^ (Y_ - \bar{Y_})^2$ dove I sono i livelli della variabile indipendente, $J_i$ il numero di osservazioni del livello i e $\bar{Y_}$ la media delle osservazioni per il livello i.
I gradi di libertà della varianza residua sono $df_R = N-I$.
La varianza residua è pari a $MS_R = SS_R / df_R$.
Somma dei quadrati e varianza spiegata
La somma dei quadrati del modello si calcola con $SS_M = \sum_^I (\bar{Y_} - \bar{Y_})^2 \cdotp J_i$ Ovvero, si calcola la differenza fra la media del gruppo i e la media totale, la si eleva al quadrato, e la si moltiplica per il numero di osservazioni di quel gruppo.
I gradi di libertà della varianza del modello sono $df_M = I-1$.
La varianza spiegata è pari a $MS_M = SS_M / df_M$.
Identità principale dell'ANOVA
Proprio come per il modello di regressione lineare, $SS_T = SS_M + SS_R$. Questa uguaglianza viene definita identità principale dell'ANOVA.
La significatività del rapporto fra la variabile indipendente e quella dipendente viene misurata mettendo a rapporto la varianza spiegata dal modello con la varianza residua: $F = MS_M / MS_R$.
Distribuzione dell'errore, inferenza
Per introdurre il calcolo dell'analisi della Varianza, usiamo la consueta sequenza logica
- identifichiamo una misura della relazione
- identifichiamo una popolazione virtuale
- estraiamo a caso $k \cdotp I$ campioni, calcoliamo la misura per ogni estrazione, e salviamo nel vettore delle misure
- osserviamo la distribuzione delle misure generate
- calcoliamo il p-value attraverso il confronto con il vettore delle misure
- identifichiamo di una distribuzione teorica che, previo opportuna trasformazione, mappa quella osservata
- calcoliamo il p-value utilizzando la probabilità della distribuzione teorica identificata
- utilizziamo la funzione di R
La simulazione
Per la nostra simulazione, immaginiamo un disegno sperimentale bivariato, dove la variabile indipendente ha tre livelli.
- La misura della relazione è la statistica F identificata sopra;
- generiamo per k volte tre campioni di numerosità m, con stessa media e deviazione standard (media e deviazione standard sono arbitrarie);
- calcoliamo $SS_T, SS_R, SS_M, df_T, df_R, df_M, MS_T, MS_R, MS_M$
- Calcoliamo la statistica $F = MS_M / MS_R$, e la salviamo nel vettore delle misure.
k<- 10000
distribuzione <- vector ("numeric",k)
for (i in 1:k) {
n <- 60
#osservazioni <- rnorm (n,100,6)
# di fatto mescolo i prezzi
osservazioni <- sample(affitti$prezzo)
# e assegno i prezzi mescolati
osservazioniA <- osservazioni [1:20]
osservazioniB <- osservazioni [21:40]
osservazioniC <- osservazioni [41:60]
# calcolo le medie dei 3 gruppi (mescolati)
meanA <- mean(osservazioniA)
meanB <- mean(osservazioniB)
meanC <- mean(osservazioniC)
# calcolo la media totale
meanTot <- mean (osservazioni)
# calcolo la somma dei quadrati dell'errore residuo
SSRA <- sum ((osservazioniA-meanA)^2)
SSRB <- sum ((osservazioniB-meanB)^2)
SSRC <- sum ((osservazioniC-meanC)^2)
SSR <- SSRA + SSRB + SSRC
# calcolo la somma dei quadrati del modello
# (varianza spiegata)
SSMA <- n/3 * (meanA-meanTot)^2
SSMB <- n/3 * (meanB-meanTot)^2
SSMC <- n/3 * (meanC-meanTot)^2
SSM <- SSMA + SSMB + SSMC
# calcolo la somma dei quadrati totale
SST <- sum((osservazioni-meanTot)^2)
# la varianza del modello (spiegata)
MSM <- SSM/(3-1)
# la varianza residua
MSR <- SSR/(60-3)
# F = varianza spiegata / varianza residua
Fvalue <- MSM/MSR
distribuzione [i] <- Fvalue
}
La distribuzione Fisher-Snedecor
hist(distribuzione,probability=TRUE)
xrange <- seq (min(distribuzione),max(distribuzione), by=.2)
lines (xrange,df(xrange, 2, 57))
Nella figura distribuzione, l'istogramma rappresenta la distribuzione dell'errore di campionamento ottenuto con la simulazione.
La linea sovrapposta all'istogramma rappresenta la distribuzione teorica F di Fisher-Snedecor.
Similmente alla distribuzione t di Student e alla distribuzione $\chi^2$, anche la F è una famiglia di distribuzioni, che variano a seconda dei gradi di libertà.
La distribuzione F, però, varia in base a due gradi di libertà. Nell'esempio della simulazione, la linea corrisponde alla distribuzione F (2,57), ovvero ai gradi di libertà della varianza spiegata e della varianza residua.
R: calcolo di F value
Mostriamo il calcolo della statistica F con R, calcolando le tre medie, la media totale, $SS_R$, $SS_M$, $SS_T$, $MS_R$, $MS_M$, ed infine F.
osservazioniA <- zonaA
osservazioniB <- zonaB
osservazioniC <- zonaC
osservazioni <- zone
meanA <- mean(osservazioniA)
meanB <- mean(osservazioniB)
meanC <- mean(osservazioniC)
meanTot <- mean (osservazioni)
SSRA <- sum ((osservazioniA-meanA)^2)
SSRB <- sum ((osservazioniB-meanB)^2)
SSRC <- sum ((osservazioniC-meanC)^2)
SSR <- SSRA + SSRB + SSRC
SSMA <- n/3 * (meanA-meanTot)^2
SSMB <- n/3 * (meanB-meanTot)^2
SSMC <- n/3 * (meanC-meanTot)^2
SSM <- SSMA + SSMB + SSMC
SST <- sum((osservazioni-meanTot)^2)
MSR <- SSR/(n-3)
MSM <- SSM/(3-1)
F_Affitti <- MSM/MSR
c(SSM,SSR,SST)
## [1] 27801.3 51451.1 79252.4
c(MSM,MSR)
## [1] 13900.6500 902.6509
F_Affitti
## [1] 15.39981
F è dunque pari a 15.3998078. Ora, possiamo calcolare il p-value nel modo consueto, confrontando la posizione di questa statistica con la distribuzione calcolata prima.
# utilizzo la distribuzione simulata
p_value_empirica_1 <- 1-rank (c(F_Affitti,distribuzione))[1]/(k+1)
p_value_empirica_1
## [1] 0
Ora, calcoliamo il p-value usando la funzione pf, ovvero in base alla distribuzione F.
# utilizzo la distribuzione statistica
(p_value_F_1 <- 1 - pf (F_Affitti,2,57))
## [1] 4.4967e-06
Il risultato è sostanzialmente simile: 0 vs 4.4966996 × 10-6.
R: uso di aov
R mette a disposizione, per il calcolo dell'analisi della varianza, la funzione $aov(y~x)$, dove y è la variabile dipendente, numerica, e x è il fattore.
aovAffitti <- aov(prezzo~quartiere,data=affitti)
summary(aovAffitti)
## Df Sum Sq Mean Sq F value Pr(>F)
## quartiere 2 27801 13901 15.4 4.5e-06 ***
## Residuals 57 51451 903
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Leggere l'output
utilizzando la funzione summary su aov è possibile avere il dettaglio dei risultati dell'analisi. Nel caso di una analisi ad una via, avremo una tabella con due righe. La seconda riga calcola i gradi di libertà, la somma dei quadrati, e la media dei quadrati dei residui (ovvero $df_R, SS_R, MS_R$). La prima riga calcola i gradi di libertà, la somma dei quadrati, e la media dei quadrati del modello: $df_M, SS_M, MS_M$); Inoltre, calcola $F = MS_M / MS_R$; infine, calcola il p-value. I codici simbolici (gli asterischi) ci suggeriscono la significatività: *** significa che $p-value < 0.01$.
Anova a due vie
Due variabili indipendenti
L'analisi della varianza che abbiamo appena introdotto può essere estesa anche ai casi in cui le variabili indipendenti sono più di una.
Nell'analisi della varianza a due vie, ad esempio, si indaga la relazione fra due variabili indipendenti, entrambe categoriali, ed una variabile dipendente, quantitativa.
In questa sezione analizziamo la circostanza in cui le variabili indipendenti sono due, ma la logica rimane la stessa anche nelle circostanze in cui le variabili indipendenti sono più di due.
Le ipotesi
Nell'Anova a due vie, le domande che il ricercatore si pone sono tre:
- La prima delle variabili indipendenti, influisce significativamente sulla variabile dipendente?
- La seconda delle variabili indipendenti, influisce significativamente sulla variabile dipendente?
- Vi è una interazione fra le due variabili?
Un esempio: antidepressivi e attività aerobica
Introduciamo l'analisi della varianza a due vie con un esempio, anche questo del tutto inventato. Dei ricercatori sono interessati ad analizzare l'influenza dell'attività aerobica e di un tipo di farmaci antidepressivi (ad esempio, gli RSSI, gli inibitori del reuptake della serotonina) sul tono dell'umore di pazienti con diagnosi di depressione maggiore.
Decidono pertanto di selezionare 120 pazienti con diagnosi di depressione, e di assegnarli casualmente a 4 gruppi sperimentali, in un disegno 2x2.
I fattori
Un fattore è dunque l'antidepressivo: a 60 pazienti verrà somministrato, per 30 giorni, una dose di RSSI, mentre agli altri 60 verrà somministrato, per lo stesso periodo, un placebo.
L'altro fattore, l'attività aerobica: a 60 pazienti verrà chiesto di fare 20 minuti di attività aerobica due volte al giorno per 30 giorni. Agli altri 60 pazienti verrà chiesto di fare una attività non aerobica, di controllo, per lo stesso periodo di tempo.
I gruppi sperimentali
Avremo dunque 30 pazienti con placebo e attività non aerobica, 30 con placebo e attività aerobica, 30 con farmaco e attività non aerobica, 30 con farmaco e attività aerobica.
Alla fine dei 30 giorni, verrà somministrato un questionario (ad esempio il BDI, Beck depression inventory), per valutare il loro tono dell'umore alla fine del trattamento.
Il calcolo
Somma dei quadrati totale e residua
La somma dei quadrati totale è identico al caso dell'anova ad una via: si somma il quadrato della differenza di ogni osservazione con la media generale. Per calcolare la varianza totale si divide il tutto per i gradi di libertà della varianza totale, pari a n-1.
Anche la somma dei quadrati residui è identico al caso dell'anova ad una via: per ogni condizione sperimentale, si sommano i quadrati delle differenze fra i valori osservati e la media di quel gruppo, si sommano i valori così ottenuti da ogni gruppo. Di nuovo, la varianza è data dalla somma dei quadrati divisa per i gradi di libertà
Somma dei quadrati dei fattori
Il calcolo di $SS_$: per ogni livello della variabile A, si calcola la media delle osservazioni di quel livello.
Si calcola il quadrato della differenza fra questa media e la media generale.
Si moltiplica questo risultato per il numero di osservazioni del livello.
Alla fine, si sommano i valori ottenuti per ognuno dei livelli. Si dividono per i gradi di libertà (pari al numero di livelli meno uno) per ottenere la varianza $MS_$
In pratica, nel calcolare $SS_$ si fa come se il fattore B non esistesse. Lo stesso procedimento viene usato per calcolare la varianza spiegata dal fattore B.
Somma dei quadrati dell'interazione
Nel caso dell'anova ad una via, la varianza totale era pari alla somma della varianza residua e della varianza spiegata dall'unica variabile indipendente.
Nel caso dell'anova a due vie, però, la somma di varianza residua, varianza spiegata da A e varianza spiegata da B sarà minore della varianza totale.
La differenza è data dalla varianza spiegata dall'interazione fra i due fattori A e B. Più forte è l'interazione fra le due variabili indipendenti, più alta sarà la varianza spiegata dall'interazione (e dunque maggiore sarà la differenza fra la somma delle varianze residue, di A, di B e la varianza totale).
L'interazione fra le variabili indipendenti
Per introdurre il calcolo della varianza spiegata dall'interazione fra A e B, può essere utile riprendere il concetto di frequenza attesa introdotta nel test del $\chi^2$. Anche in quel caso si trattava di valutare l'interazione fra due variabili categoriali. La differenza è che mentre nel $\chi^2$ si misurano le frequenze, in questo caso la misura è data da una variabile quantitativa.
In maniera simile al $\chi^2$, però, è possibile, conoscendo le medie marginali dei livelli di A e B, costruire una tabella delle medie attese, che costituisce il caso perfetto di assenza di interazione fra i due fattori. Questa tabella costituisce il caso ottimale di accettazione dell'ipotesi nulla relativa all'interazione fra le due variabili.
L'esempio: calcolo delle somme dei quadrati
Per esemplificare, torniamo all'esempio di un disegno 2x2. Quello che vogliamo fare è creare una tabella 2x2 delle medie attese.
Il primo passaggio consiste nel calcolare le medie marginali per ogni livello della variabile A, e dunque la media di A1 e la media di A2.
Nel secondo passaggio si calcolano le medie marginali per i livelli di B, e dunque la media di B1 e la media di B2.
Infine, per le quattro celle [1,1] [1,2] [2,1] e [2,2], si calcola la media attesa.
La media attesa della cella [1,1] è pari alla media generale + la differenza fra A1 e la media generale e la differenza fra B1 e la media generale. Dunque, A1 + B1 - media.
Calcolo delle somme dell'interazione
$SS_$ si basa sul quadrato delle differenze fra la tabella delle medie attese e la tabella delle medie osservate, moltiplicata per il numero di osservazioni per gruppo.
Dunque, più la tabella delle medie osservate è simile alla tabelle delle medie attese, minore è l'interazione fra le due variabili indipendenti, e dunque minore è la varianza spiegata dall'interazione.
Viceversa, maggiore è la differenza, maggiore l'interazione, maggiore la varianza spiegata dall'interazione.
Il calcolo, formalizzazione
Somma dei quadrati e varianza totale
La somma dei quadrati dell'errore totale si calcola con la formula $SS_T = \sum_^N (Y_i - \bar{Y_})^2$ dove N è il numero totale di osservazioni e $\bar{Y_}$ è la media totale.
I gradi di libertà della varianza totale sono $df_T = N - 1$.
La varianza totale è pari a $MS_T = SS_T / df_T$.
Somma dei quadrati e varianza residua
La somma dei quadrati dell'errore residuo si calcola con la formula $SS_R = \sum_^I \sum_^J \sum_^K (Y_ - \bar{Y_})^2$ dove I sono i livelli di A, J i livelli di B, K il numero di osservazioni per ogni gruppo e $\bar{Y_}$ la media delle osservazioni per il gruppo ij.
I gradi di libertà della varianza residua sono $df_T = N-I*J$.
La varianza residua è pari a $MS_R = SS_R / df_R$.
Somma dei quadrati e varianza spiegata
La somma dei quadrati del modello si calcola con $SS_A = K \cdotp J \cdotp \sum_^I (\bar{Y_} - \bar{Y_})^2$ $SS_B = K \cdotp I \cdotp \sum_^J (\bar{Y_} - \bar{Y_})^2$ I gradi di libertà della varianza del modello sono $df_A = I-1, df_B = J-1$.
Le varianze spiegate sono $MS_A = SS_A / df_A$, $MS_B = SS_B / df_B$.
Somma dei quadrati e varianza dell'interazione
$SS_ = K \cdotp \sum_^I \sum_^J (\bar{Y_} - \bar{Y_} - \bar{Y_} + \bar{Y_})^2$ ovvero la media osservata meno la media marginale di $A_i$, meno la media marginale di $B_j$, più la media totale.
I gradi di libertà della varianza dell'interazione sono $df_ = (I-1) \cdotp (J-1)$.
La varianza dell'interazione è $MS_ = SS_ / df_$.
Le ipotesi inferenziali
Le ipotesi inferenziali sono tre:
- $H_0^A$: l'influenza della variabile A sulla variabile dipendente non è significativa
- $H_0^B$: l'influenza della variabile B sulla variabile dipendente non è significativa
- $H_0^$: l'interazione fra A e B non è significativa.
Ipotesi inferenziali e F
Per valutare le tre ipotesi inferenziali, vengono calcolati i tre rapporti:
- $F_A = MS_A / MS_R$;
- $F_B = MS_B / MS_R$;
- $F_ = MS_ / MS_R$;
Ognuno dei rapporti viene confrontato con la distribuzione F di Fisher-Snedecor.
Modello lineare
In maniera simile alla regressione lineare, anche l'analisi della varianza (sia a una che a due vie) può essere rappresentata attraverso un modello lineare. Il modello generale per l'anova a due vie è $Y_ = \mu + \alpha_i + \beta_j + \delta_ + \epsilon_, i = 1 ...I, j = 1 ... J, k = 1 ... K$ dove I è il numero di livelli del fattore A, J il numero di livelli del fattore B, K il numero di osservazioni per ogni gruppo.
$\mu$ corrisponde alla media totale di tutte le osservazioni.
$\alpha_i$ corrisponde allo scostamento dalla media totale del livello $A_i$
$\beta_j$ corrisponde allo scostamento dalla media totale del livello $B_j$
$\delta_$ corrisponde alla differenza fra la media del campione osservata e quella attesa in base all'ipotesi di non interazione fra A e B.
$\epsilon_$ è la componente di errore, ovvero la differenza fra il valore atteso dal modello e il valore osservato.
L'esempio dei trattamenti per la depressione
Torniamo all'esempio dei trattamenti per la depressione, e generiamo tre diversi scenari.
Primo scenario
punteggioBase <- 22
effettoFarmaco <- 8
effettoAerobico <- 7
# gl(n,k,length) genera un fattore di n livelli, k volte
# farmaco: due livelli ripetuti 60 volte (=120)
farmaco <- gl (2,60,labels=c("placebo","farmaco"))
# aerobica: due livelli ripetuti 30 volte, ma per un totale di 120
# (e dunque i due da 30 ripetuti due volte)
aerobica <- gl (2,30,120,labels=c("non aerobico","aerobico"))
placeboNA <- rnorm (30,punteggioBase,8)
placeboA <- rnorm (30,punteggioBase+effettoAerobico,8)
farmacoNA <- rnorm (30,punteggioBase+effettoFarmaco,8)
farmacoA <- rnorm (30,punteggioBase+effettoFarmaco+effettoAerobico,8)
punteggi <- c(placeboNA,placeboA,farmacoNA,farmacoA)
expDep1 <- data.frame (punteggi=punteggi,farmaco=farmaco,aerobica=aerobica)
medie.osservate <- tapply (expDep1$punteggi,list(expDep1$farmaco,expDep1$aerobica),mean)
medie.osservate
## non aerobico aerobico
## placebo 21.37817 27.56091
## farmaco 30.52215 36.59854
Rappresentiamo le medie dei quattro gruppi (2x2) usando la funzione interaction.plot
interaction.plot (expDep1$farmaco,expDep1$aerobica,expDep1$punteggi)
La funzione interaction.plot rappresenta la media della variabile dipendente (numerica) nelle diverse condizioni sperimentali in un disegno con due variabili indipendenti categoriche.
Dal grafico si intuisce che entrambe le variabili indipendenti hanno un effetto sulla variabile dipendente, ma che probabilmente non vi sarà interazione fra le due.
Applichiamo la funzione aov
aovScenario1 <- aov(punteggi~farmaco+aerobica+farmaco:aerobica,
data = expDep1)
(summaryAov1 <- summary(aovScenario1))
## Df Sum Sq Mean Sq F value Pr(>F)
## farmaco 1 2479 2479.3 38.479 8.79e-09 ***
## aerobica 1 1127 1127.1 17.494 5.63e-05 ***
## farmaco:aerobica 1 0 0.1 0.001 0.971
## Residuals 116 7474 64.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov()) genera una tabella che riassume i gradi di libertà, la somma dei quadrati, la media, la statistica F ed il p-value del fattore farmaco, del fattore aerobica, dell'interazione farmaco:aerobica e dei residui.
Dal modello dell'analisi della varianza possiamo dedurre che vi è una influenza significativa sia del primo fattore (farmaco) che del secondo (attività aerobica); non vi è, però, interazione significativa fra i due fattori.
L'esempio, secondo scenario
In questo esempio manipoliamo i valori del gruppo farmaco e attività aerobica, aggiungendo 10 a tutti i punteggi. In questo modo stiamo forzando l'interazione fra le due variabili indipendenti.
# aggiungo 10 al gruppo farmacoA(erobica)
#farmacoA <- rnorm (30,punteggioBase+effettoFarmaco+effettoAerobico+10,8)
farmacoAmod <- farmacoA+10
punteggi <- c(placeboNA,placeboA,farmacoNA,farmacoAmod)
expDep2 <- data.frame (punteggi=punteggi,farmaco=farmaco,aerobica=aerobica)
summary (expDep2)
## punteggi farmaco aerobica
## Min. : 4.795 placebo:60 non aerobico:60
## 1st Qu.:23.006 farmaco:60 aerobico :60
## Median :28.744
## Mean :31.515
## 3rd Qu.:41.296
## Max. :58.839
interaction.plot (expDep2$farmaco,expDep2$aerobica,expDep2$punteggi)
medie.osservate <- tapply (expDep2$punteggi,list(expDep2$farmaco,expDep2$aerobica),mean)
medie.osservate
## non aerobico aerobico
## placebo 21.37817 27.56091
## farmaco 30.52215 46.59854
aov_sc2 <- aov(punteggi~farmaco+aerobica+farmaco:aerobica)
risultati <- matrix(unlist(summary(aov_sc2)),ncol = 5)
summary(aov_sc2)
## Df Sum Sq Mean Sq F value Pr(>F)
## farmaco 1 5957 5957 92.45 < 2e-16 ***
## aerobica 1 3716 3716 57.67 8.61e-12 ***
## farmaco:aerobica 1 734 734 11.39 0.001 **
## Residuals 116 7474 64
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In questo caso
- rifiutiamo l'ipotesi nulla $H_0^A$, in quanto l'effetto del farmaco è significativo: $F=92.4471649$, $p= <0.001$;
- rifiutiamo l'ipotesi nulla $H_0^B$, in quanto l'effetto dell'attività aerobica è significativo: $F=57.673775$, $p= <0.001$;
- rifiutiamo l'ipotesi nulla $H_0^$, in quanto è significativa anche l'interazione fra i due fattori: $F=11.3939936$, $p= 0.001$;
Confronti multipli
Confronti multipli ed errore
L'analisi della varianza ci permette di verificare se le differenze fra le medie di tre o più campioni sono da attribuire all'errore campionario, o se sono significative.
Una volta rifiutata l'ipotesi nulla, però, resta da determinare quali differenze sono significative. L'analisi della varianza, infatti, ci dice se vi è almeno una coppia di gruppi la cui differenza è significativa, ma non ci dice quali differenze lo sono.
Per poter determinare quali differenze sono significative, diventa necessario confrontare i gruppi a due a due.
Come abbiamo visto all'inizio del capitolo, si potrebbe decidere di utilizzare, per confrontare a due a due i diversi gruppi, il t-test. Ma, come abbiamo già accennato, applicare ripetutamente il t-test aumenta la probabilità di incorrere in un errore del primo tipo.
Diventa dunque necessario adottare dei test di confronti multipli capaci di mantenere sotto controllo l'errore del I tipo.
La correzione di Bonferroni
Un possibile approccio, finalizzato a controllare la proliferazione dell'errore del I tipo, è quello di adottare la correzione di Bonferroni, che consiste nel dividere il valore $\alpha$ per il numero di confronti effettuati (o, in maniera corrispondente, moltiplicare il p-value per il numero di confronti).
Se, ad esempio, dobbiamo confrontare le medie di 4 gruppi e decidiamo per un valore $\alpha = 0.05$, in base alla correzione di Bonferroni dovremo considerare significativi soltanto quei confronti il cui p-value sia inferiore a $0.05/6=0.0083$ (in quanto i confronti previsti sono 6).
Il problema di questo metodo è che tende ad essere eccessivamente conservativo.
Il test di Tukey
Un metodo di confronto multiplo meno conservativo è il test di Tukey.
Anche attraverso il test di Tukey è possibile mantenere l'errore di tipo I entro un predeterminato valore di $\alpha$ (generalmente pari a 0.05).
Il test di Tukey permette di correggere il p-value in base al numero di confronti che vengono effettuati nel confronto multiplo, senza però penalizzare eccessivamente la statistica.
Per correggere l'errore, il test di Tukey confronta la statistica calcolata con la distribuzione studentized range
Il test di Tukey: calcolo
Per calcolare la significatività della differenza fra due gruppi con il metodo di Tukey si utilizza il seguente algoritmo:
- si calcola l'errore standard, con la formula $SE = \sqrt {(MS_R/n)}$ dove n è il numero di osservazioni per gruppo. Nel caso di gruppi con numerosità diversa, la formula diventa $SE = \sqrt { \frac{2} \cdotp (\frac{1} + \frac{1})}$ dove $n_a$ e $n_b$ sono la numerosità del primo e del secondo gruppo
- si calcola la statistica $Q = \frac{| Y_a - Y_b |}$
- si calcola il p-value, usando la funzione
ptukey (Q,k,Df_R)dove k è il numero di confronti effettuati, e $Df_R$ i gradi di libertà della varianza residua. La funzione ptukey calcola la probabilità sulla distribuzione studentized range.
Affitti: confronti multipli
Torniamo all'esempio degli affitti, e mostriamo il calcolo di uno dei confronti multipli con il test di Tukey. Calcoliamo le tre medie.
medieAffitti <- tapply(affitti$prezzo,affitti$quartiere,mean)
medieAffitti
## A B C
## 555.20 568.55 606.05
Tukey: il calcolo di un confronto
Calcoliamo il conftonto fra le medie dei gruppi A e B.
confronto1 <- abs(medieAffitti[1] - medieAffitti[2])
SE <- sqrt (MSR/20)
Q <- confronto1/SE
p_value <- 1 -ptukey (Q,3,57)
c(confronto1,SE,Q,p_value)
## A A A
## 13.350000 6.718076 1.987176 0.345060
Il p-value è alto, e dunque la differenza fra i gruppi A e B non è significativa.
La funzione R TukeyHSD
La funzione di R per il calcolo del confronto con il metodo Tukey è TukeyHSD. Coerentemente con l'uso dei confronti multipli, la funzione si applica sul risultato della corrispondente analisi della varianza.
(confronti1 <- TukeyHSD(aovAffitti, ordered = TRUE))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
## factor levels have been ordered
##
## Fit: aov(formula = prezzo ~ quartiere, data = affitti)
##
## $quartiere
## diff lwr upr p adj
## B-A 13.35 -9.512883 36.21288 0.3450600
## C-A 50.85 27.987117 73.71288 0.0000048
## C-B 37.50 14.637117 60.36288 0.0006350
La funzione ritorna una tabella, con una riga per ogni confronto, dove vengono mostrate:
- la coppia confrontata (es, il confronto fra il gruppo B ed il gruppo A); l'ordine è tale che il gruppo con media più alta è davanti all'altro;
- la differenza (positiva) fra i due gruppi;
- l'intervallo di confidenza della differenza; ad esempio, nel secondo confronto (C-A), la differenza è di 50.85, l'intervallo di confidenza va da un minimo di 27.9871167 ad un massimo di 73.7128833. p adj è il p-value aggiustato, che nel confronto
C - Aè <0.001; nell'esempio, i confronti C-A e C-B sono significativi, il confronto B-A no (0.35).
Analisi della Varianza: assunti
Come ogni approccio parametrico, anche l'analisi della varianza fa delle assunti:
- indipendenza delle osservazioni
- distribuzione normale dei residui
- omoschedasticità: la varianza dell'errore è costante
- gli errori sono fra loro indipendenti
Distribuzione dei residui
Si assume che i residui abbiano una distribuzione normale, con media pari a 0, e varianza costante fra i gruppi. Per testare l'ipotesi di normalità, è possibile usare il test di Shapiro-Wilk sui residui del modello dell'analisi della varianza:
shapiro.test(aovAffitti$residuals)
##
## Shapiro-Wilk normality test
##
## data: aovAffitti$residuals
## W = 0.9937, p-value = 0.9896
Per testare l'ipotesi di omoschedasticità, si può usare il test di Bartlett:
bartlett.test(prezzo ~ quartiere, data = affitti)
##
## Bartlett test of homogeneity of variances
##
## data: prezzo by quartiere
## Bartlett's K-squared = 0.59846, df = 2, p-value = 0.7414
Test non parametrico
Vi sono circostanze in cui l'analisi della varianza non può essere applicata, in quanto vengono meno alcuni assunti o condizioni:
- non si può assumere la normalità della distribuzione degli errori
- il numero di osservazioni per ogni gruppo è minore di 10
- la variabile dipendente non è ad intervalli, ma ordinale
In questi casi è possibile applicare il test non parametrico di Kruskal-Wallis
Il test di Kruskal-Wallis
il test di Kruskal-Wallis è un'estensione del test di Wilcoxon, che abbiamo visto nel capitolo dedicato al t-test [@ChanWalmsley:1997]. Nel test di Kruskal-Wallis, la prima operazione da compiere è quella di trasformare i punteggi osservati nel loro rango. A questo punto, si applica la formula $K = (n-1)\frac{\sum_^I n_i (\bar{r_} - \bar)^2}{\sum_^I \sum_^ (r_ - \bar)^2} \label$ dove $n_i$ è il numero di osservazioni nel gruppo i, $r_$ è la posizione ordinale dell'osservazione j del gruppo i, N è il numero totale delle osservazioni, $\bar{r_}$ è la media dei rank del gruppo i.
Semplificazioni
L'equazione Kruskal può in realtà essere semplificata, in quanto $\bar = (N+1)/2$ e denominatore è pari a $(N-1) N (N+1) / 12$, e dunque otteniamo
$K = \frac{12}{N(N+1)} \sum_^I n_i (\bar{r_} - \frac{N+1}{2})^2 \label$
La statistica K assume una distribuzione $\chi^2$ con i-1 gradi di libertà.
R: la funzione kruskal.test
Applichiamo il test di Kruskal-Wallis al nostro esempio degli affitti.
(kt1 <- kruskal.test(prezzo ~ quartiere, data = affitti))
##
## Kruskal-Wallis rank sum test
##
## data: prezzo by quartiere
## Kruskal-Wallis chi-squared = 20.062, df = 2, p-value = 4.401e-05
Leggere i risultati
La funzione restituisce la statistica, Kruskal-Wallis chi-squared = 20.0623121; I gradi di libertà: df = 2; il p-value = 4.4007256 × 10-5.
Conclusioni
Da fare.
[^1]: questo tipo di statistica viene definito omnibus