Come eseguire l'anova unidirezionale in r

Un’ANOVA unidirezionale viene utilizzata per determinare se esiste o meno una differenza statisticamente significativa tra le medie di tre o più gruppi indipendenti.

Questo tipo di test è chiamato ANOVA unidirezionale perché analizziamo l’impatto di una variabile predittrice su una variabile di risposta.

Nota : se fossimo invece interessati all’impatto di due variabili predittive su una variabile di risposta, potremmo eseguire un’ANOVA a due vie .

Come eseguire l’ANOVA unidirezionale in R

L’esempio seguente illustra come eseguire un’ANOVA unidirezionale in R.

Sfondo

Supponiamo di voler determinare se tre diversi programmi di esercizi hanno un impatto diverso sulla perdita di peso. La variabile predittiva che studiamo è il programma di esercizi e la variabile di risposta è la perdita di peso, misurata in libbre.

Possiamo eseguire un’ANOVA unidirezionale per determinare se esiste una differenza statisticamente significativa tra la perdita di peso risultante dai tre programmi.

Reclutiamo 90 persone per partecipare a un esperimento in cui assegniamo casualmente 30 persone a seguire il Programma A, il Programma B o il Programma C per un mese.

Il codice seguente crea il frame di dati con cui lavoreremo:

 #make this example reproducible
set.seed(0)

#create data frame
data <- data.frame(program = rep(c("A", "B", "C"), each = 30),
                   weight_loss = c(runif(30, 0, 3),
                                   runif(30, 0, 5),
                                   runif(30, 1, 7)))

#view first six rows of data frame
head(data)

# program weight_loss
#1 A 2.6900916
#2 A 0.7965260
#3 A 1.1163717
#4 A 1.7185601
#5 A 2.7246234
#6 A 0.6050458

La prima colonna del riquadro dati mostra il programma a cui la persona ha partecipato per un mese e la seconda colonna mostra la perdita di peso totale sperimentata dalla persona alla fine del programma, misurata in libbre.

Esplora i dati

Prima ancora di adattare il modello ANOVA unidirezionale, possiamo comprendere meglio i dati trovando la media e la deviazione standard della perdita di peso per ciascuno dei tre programmi utilizzando il pacchetto dplyr :

 #load dplyr package
library (dplyr)

#find mean and standard deviation of weight loss for each treatment group
data %>%
  group_by (program) %>%
  summarize (mean = mean(weight_loss),
            sd = sd(weight_loss))

# A tibble: 3 x 3
# program mean sd
#      
#1 A 1.58 0.905
#2 B 2.56 1.24 
#3 C 4.13 1.57  

Possiamo anche creare un boxplot per ciascuno dei tre programmi per visualizzare la distribuzione della perdita di peso per ciascun programma:

 #create boxplots
boxplot(weight_loss ~ program,
data = data,
main = "Weight Loss Distribution by Program",
xlab = "Program",
ylab = "Weight Loss",
col = "steelblue",
border = "black")

Da questi box plot, possiamo vedere che la perdita di peso media è più alta per i partecipanti al Programma C e che la perdita di peso media è più bassa per i partecipanti al Programma A.

Possiamo anche vedere che la deviazione standard (la “lunghezza” del boxplot) per la perdita di peso è leggermente più alta nel programma C rispetto agli altri due programmi.

Successivamente, adatteremo il modello ANOVA unidirezionale ai nostri dati per vedere se queste differenze visive sono effettivamente statisticamente significative.

Raccordo del modello ANOVA unidirezionale

La sintassi generale per adattare un modello ANOVA unidirezionale in R è:

aov(variabile di risposta ~ variabile_predittore, dati = set di dati)

Nel nostro esempio, possiamo utilizzare il codice seguente per adattare il modello ANOVA unidirezionale, utilizzando perdita_peso come variabile di risposta e programma come variabile predittrice. Possiamo quindi utilizzare la funzione summary() per visualizzare il risultato del nostro modello:

 #fit the one-way ANOVA model
model <- aov(weight_loss ~ program, data = data)

#view the model output
summary(model)

# Df Sum Sq Mean Sq F value Pr(>F)    
#program 2 98.93 49.46 30.83 7.55e-11 ***
#Residuals 87 139.57 1.60                     
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dai risultati del modello, possiamo vedere che il programma delle variabili predittive è statisticamente significativo al livello di significatività 0,05.

In altre parole, esiste una differenza statisticamente significativa tra la perdita di peso media risultante dai tre programmi.

Verifica delle ipotesi del modello

Prima di andare oltre, dobbiamo verificare che le ipotesi del nostro modello siano soddisfatte in modo che i risultati del nostro modello siano affidabili. In particolare, un’ANOVA unidirezionale presuppone:

1. Indipendenza – le osservazioni di ciascun gruppo devono essere indipendenti l’una dall’altra. Poiché abbiamo utilizzato un disegno randomizzato (ovvero, abbiamo assegnato i partecipanti ai programmi di esercizi in modo casuale), questo presupposto dovrebbe essere soddisfatto in modo da non doverci preoccupare troppo.

2. Normalità : la variabile dipendente dovrebbe avere una distribuzione approssimativamente normale per ciascun livello della variabile predittrice.

3. Varianza uguale : le varianze per ciascun gruppo sono uguali o approssimativamente uguali.

Un modo per verificare le ipotesi di normalità e di uguale varianza è utilizzare la funzione plot() , che produce quattro grafici di controllo del modello. In particolare, ci interessano particolarmente i seguenti due lotti:

• Residui vs. adattato : questo grafico mostra la relazione tra residui e valori adattati. Possiamo utilizzare questo grafico per valutare approssimativamente se la varianza tra i gruppi è approssimativamente uguale o meno.
• Grafico QQ : questo grafico mostra i residui standardizzati rispetto ai quantili teorici. Possiamo utilizzare questo grafico per valutare approssimativamente se il presupposto di normalità è soddisfatto o meno.

Il seguente codice può essere utilizzato per produrre questi grafici di controllo del modello:

 plot(model)

Il grafico QQ sopra ci consente di verificare l’ipotesi di normalità. Idealmente, i residui standardizzati dovrebbero trovarsi lungo la linea diagonale del grafico. Tuttavia, nel grafico sopra possiamo vedere che i residui deviano leggermente dalla linea verso l’inizio e la fine. Ciò indica che il nostro presupposto di normalità potrebbe essere violato.

I residui vs. Il grafico corretto qui sopra ci consente di verificare la nostra ipotesi di varianze uguali. Idealmente, vorremmo che i residui fossero distribuiti equamente per ciascun livello dei valori adattati.

Possiamo vedere che i residui sono molto più dispersi per i valori adattati più alti, indicando che la nostra assunzione di uguaglianza delle varianze potrebbe essere violata.

Per testare formalmente le varianze uguali, potremmo eseguire il test di Levene utilizzando il pacchetto auto :

 #load car package
library (car)

#conduct Levene's Test for equality of variances
leveneTest(weight_loss ~ program, data = data)

#Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)  
#group 2 4.1716 0.01862 *
#87                  
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il valore p del test è 0,01862 . Se utilizzassimo un livello di significatività pari a 0,05, rifiuteremmo l’ipotesi nulla secondo cui le varianze sono uguali nei tre programmi. Tuttavia, se utilizziamo un livello di significatività pari a 0,01, non rifiuteremo l’ipotesi nulla.

Sebbene potremmo tentare di trasformare i dati per garantire che le nostre ipotesi di normalità e uguaglianza delle varianze siano soddisfatte, per ora non ce ne preoccuperemo troppo.

Analizzare le differenze di trattamento

Una volta verificato che le ipotesi del modello sono soddisfatte (o ragionevolmente soddisfatte), possiamo eseguire un test post hoc per determinare esattamente quali gruppi di trattamento differiscono l’uno dall’altro.

Per il nostro test post hoc, utilizzeremo la funzione TukeyHSD() per eseguire il test Tukey per confronti multipli:

 #perform Tukey's Test for multiple comparisons
TukeyHSD(model, conf.level=.95) 

#Tukey multiple comparisons of means
# 95% family-wise confidence level
#
#Fit: aov(formula = weight_loss ~ program, data = data)
#
#$program
# diff lwr upr p adj
#BA 0.9777414 0.1979466 1.757536 0.0100545
#CA 2.5454024 1.7656076 3.325197 0.0000000
#CB 1.5676610 0.7878662 2.347456 0.0000199

Il valore p indica se esiste o meno una differenza statisticamente significativa tra ciascun programma. I risultati mostrano che esiste una differenza statisticamente significativa tra la perdita di peso media di ciascun programma al livello di significatività di 0,05.

Possiamo anche visualizzare gli intervalli di confidenza al 95% risultanti dal test di Tukey utilizzando la funzione plot(TukeyHSD()) in R:

 #create confidence interval for each comparison
plot(TukeyHSD(model, conf.level=.95), las = 2)

I risultati degli intervalli di confidenza sono coerenti con i risultati dei test di ipotesi.

In particolare, possiamo vedere che nessuno degli intervalli di confidenza per la perdita di peso media tra i programmi contiene il valore zero , indicando che esiste una differenza statisticamente significativa nella perdita media tra i tre programmi.

Ciò è coerente con tutti i valori p per i nostri test di ipotesi inferiori a 0,05.

Segnalazione dei risultati ANOVA unidirezionali

Infine, possiamo riportare i risultati dell’ANOVA unidirezionale in un modo che riassume i risultati:

È stata eseguita un’ANOVA unidirezionale per esaminare gli effetti del programma di esercizi   sulla perdita di peso (misurata in libbre). È stata riscontrata una differenza statisticamente significativa tra gli effetti dei tre programmi sulla perdita di peso (F(2, 87) = 30,83, p = 7,55e-11). Sono stati eseguiti i test HSD post-hoc di Tukey.

La perdita di peso media dei partecipanti al programma C è significativamente maggiore della perdita di peso media dei partecipanti al programma B (p < 0,0001).

La perdita di peso media dei partecipanti al programma C è significativamente maggiore della perdita di peso media dei partecipanti al programma A (p < 0,0001).

Inoltre, la perdita di peso media dei partecipanti al programma B è stata significativamente maggiore della perdita di peso media dei partecipanti al programma A (p = 0,01).

Risorse addizionali

Le esercitazioni seguenti forniscono informazioni aggiuntive sugli ANOVA unidirezionali:

Un’introduzione all’ANOVA unidirezionale
Una guida all’utilizzo dei test post-hoc con ANOVA
La guida completa: come riportare i risultati ANOVA