Come eseguire l'anova bidirezionale in r

Un’ANOVA a due vie (“analisi della varianza”) viene utilizzata per determinare se esiste o meno una differenza statisticamente significativa tra le medie di tre o più gruppi indipendenti che sono stati suddivisi in due fattori.

Questo tutorial spiega come eseguire un’ANOVA bidirezionale in R.

Esempio: ANOVA a due vie in R

Diciamo che vogliamo determinare se l’intensità dell’esercizio e il sesso influiscono sulla perdita di peso. In questo caso, i due fattori che stiamo esaminando sono l’esercizio fisico e il sesso , e la variabile di risposta è la perdita di peso, misurata in libbre.

Possiamo eseguire un’ANOVA a due vie per determinare se l’esercizio fisico e il sesso influiscono sulla perdita di peso e per determinare se esiste un’interazione tra esercizio fisico e genere sulla perdita di peso.

Stiamo reclutando 30 uomini e 30 donne per partecipare a un esperimento in cui assegniamo casualmente 10 di ciascuno a seguire un programma di non esercizio fisico, esercizio leggero o esercizio intenso per un mese.

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

 #make this example reproducible
set.seed(10)

#create data frame
data <- data.frame(gender = rep(c("Male", "Female"), each = 30),
                   exercise = rep(c("None", "Light", "Intense"), each = 10, times = 2),
                   weight_loss = c(runif(10, -3, 3), runif(10, 0, 5), runif(10, 5, 9),
                                   runif(10, -4, 2), runif(10, 0, 3), runif(10, 3, 8)))

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

# gender exercise weight_loss
#1 Male None 0.04486922
#2 Male None -1.15938896
#3 Male None -0.43855400
#4 Male None 1.15861249
#5 Male None -2.48918419
#6 Male None -1.64738030

#see how many participants are in each group
table(data$gender, data$exercise)

# Intense Light None
# Female 10 10 10
# Male 10 10 10

Esplora i dati

Prima ancora di adattare il modello ANOVA a due vie, possiamo comprendere meglio i dati trovando la media e la deviazione standard della perdita di peso per ciascuno dei sei gruppi di trattamento utilizzando il pacchetto dplyr :

 #load dplyr package
library(dplyr)

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

# A tibble: 6 x 4
# Groups: gender [2]
# gender exercise means sd
#          
#1 Female Intense 5.31 1.02 
#2 Female Light 0.920 0.835
#3 Female None -0.501 1.77 
#4 Male Intense 7.37 0.928
#5 Male Light 2.13 1.22 
#6 Male None -0.698 1.12 

Possiamo anche creare un boxplot per ciascuno dei sei gruppi di trattamento per visualizzare la distribuzione della perdita di peso per ciascun gruppo:

 #set margins so that axis labels on boxplot don't get cut off
by(mar=c(8, 4.1, 4.1, 2.1))

#create boxplots
boxplot(weight_loss ~ gender:exercise,
data = data,
main = "Weight Loss Distribution by Group",
xlab = "Group",
ylab = "Weight Loss",
col = "steelblue",
border = "black", 
las = 2 #make x-axis labels perpendicular
)

Possiamo subito notare che i due gruppi che hanno partecipato ad esercizio intenso sembrano avere valori di perdita di peso più elevati. Possiamo anche vedere che gli uomini tendono ad avere valori di perdita di peso più elevati rispetto alle donne sia nei gruppi di esercizi intensi che in quelli leggeri .

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

Adattamento del modello ANOVA a due vie

La sintassi generale per adattare un modello ANOVA a due vie in R è:

aov(variabile di risposta ~predictor_variable1 *predictor_variable2, data = dataset)

Si noti che l’asterisco * tra le due variabili predittive indica che si desidera testare anche un effetto di interazione tra le due variabili predittive.

Nel nostro esempio, possiamo utilizzare il codice seguente per adattare il modello ANOVA a due vie, utilizzando perdita_peso come variabile di risposta e sesso ed esercizio fisico come due variabili predittive.

Possiamo quindi utilizzare la funzione summary() per visualizzare il risultato del nostro modello:

 #fit the two-way ANOVA model
model <- aov(weight_loss ~ gender * exercise, data = data)

#view the model output
summary(model)

# Df Sum Sq Mean Sq F value Pr(>F)    
#gender 1 15.8 15.80 11.197 0.0015 ** 
#exercise 2 505.6 252.78 179.087 <2e-16 ***
#gender:exercise 2 13.0 6.51 4.615 0.0141 *  
#Residuals 54 76.2 1.41                   
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dai risultati del modello, possiamo vedere che il genere , l’esercizio fisico e l’interazione tra le due variabili sono tutti statisticamente significativi al livello di significatività di 0,05.

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 a due vie presuppone:

1. Indipendenza – le osservazioni di ciascun gruppo devono essere indipendenti l’una dall’altra. Poiché abbiamo utilizzato un disegno randomizzato , questo presupposto dovrebbe essere soddisfatto, quindi non dobbiamo preoccuparcene troppo.

2. Normalità – la variabile dipendente dovrebbe avere una distribuzione approssimativamente normale per ciascuna combinazione di gruppi dei due fattori.

Un modo per verificare questa ipotesi è creare un istogramma dei residui del modello. Se i residui sono distribuiti approssimativamente normalmente, questa ipotesi dovrebbe essere soddisfatta.

 #define model residuals
reside <- model$residuals

#create histogram of residuals
hist(resid, main = "Histogram of Residuals", xlab = "Residuals", col = "steelblue")

I residui sono distribuiti approssimativamente normalmente, quindi possiamo supporre che il presupposto di normalità sia soddisfatto.

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

Un modo per verificare questa ipotesi è eseguire un test di Levene per l’uguaglianza delle varianze utilizzando il pacchetto auto :

 #load car package
library(car)

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

#Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
#group 5 1.8547 0.1177
#54  

Poiché il valore p del test è maggiore del nostro livello di significatività pari a 0,05, possiamo presumere che la nostra assunzione di uguaglianza delle varianze tra i gruppi sia soddisfatta.

Analizzare le differenze di trattamento

Una volta verificato che le ipotesi del modello sono 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 ~ gender * exercise, data = data)
#
#$gender
# diff lwr upr p adj
#Male-Female 1.026456 0.4114451 1.641467 0.0014967
#
#$exercise
# diff lwr upr p adj
#Light-Intense -4.813064 -5.718493 -3.907635 0.0e+00
#None-Intense -6.938966 -7.844395 -6.033537 0.0e+00
#None-Light -2.125902 -3.031331 -1.220473 1.8e-06
#
#$`gender:exercise`
# diff lwr upr p adj
#Male:Intense-Female:Intense 2.0628297 0.4930588 3.63260067 0.0036746
#Female:Light-Female:Intense -4.3883563 -5.9581272 -2.81858535 0.0000000
#Male:Light-Female:Intense -3.1749419 -4.7447128 -1.60517092 0.0000027
#Female:None-Female:Intense -5.8091131 -7.3788841 -4.23934219 0.0000000
#Male:None-Female:Intense -6.0059891 -7.5757600 -4.43621813 0.0000000
#Female:Light-Male:Intense -6.4511860 -8.0209570 -4.88141508 0.0000000
#Male:Light-Male:Intense -5.2377716 -6.8075425 -3.66800066 0.0000000
#Female:None-Male:Intense -7.8719429 -9.4417138 -6.30217192 0.0000000
#Male:None-Male:Intense -8.0688188 -9.6385897 -6.49904786 0.0000000
#Male:Light-Female:Light 1.2134144 -0.3563565 2.78318536 0.2185439
#Female:None-Female:Light -1.4207568 -2.9905278 0.14901410 0.0974193
#Male:None-Female:Light -1.6176328 -3.1874037 -0.04786184 0.0398106
#Female:None-Male:Light -2.6341713 -4.2039422 -1.06440032 0.0001050
#Male:None-Male:Light -2.8310472 -4.4008181 -1.26127627 0.0000284
#Male:None-Female:None -0.1968759 -1.7666469 1.37289500 0.9990364

Il valore p indica se esiste o meno una differenza statisticamente significativa tra ciascun gruppo.

Ad esempio, nell’ultima riga sopra, vediamo che il gruppo di uomini senza esercizio non ha riscontrato una differenza statisticamente significativa nella perdita di peso rispetto al gruppo di donne senza esercizio (valore p: 0,990364).

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

 #set axis margins so labels don't get cut off
by(mar=c(4.1, 13, 4.1, 2.1))

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

Segnalazione dei risultati ANOVA bidirezionali

Infine, possiamo riportare i risultati dell’ANOVA a due vie in un modo che riassume i risultati:

È stata eseguita un’ANOVA a due vie per esaminare gli effetti del genere ( maschio, femmina) e del programma di esercizi (nessuno, leggero, intenso) sulla perdita di peso (misurata in libbre). È stata riscontrata un’interazione statisticamente significativa tra gli effetti del genere e dell’esercizio fisico sulla perdita di peso (F(2, 54) = 4,615, p = 0,0141). Sono stati eseguiti i test HSD post-hoc di Tukey.

Per gli uomini, un programma di esercizi intensi ha comportato una perdita di peso significativamente maggiore rispetto a un programma leggero (p < 0,0001) o a nessun programma di esercizi (p < 0,0001). Inoltre, negli uomini, una dieta leggera ha comportato una perdita di peso significativamente maggiore rispetto all’assenza di un regime di esercizio fisico (p < 0,0001).

Per le donne, un programma di esercizi intensi ha comportato una perdita di peso significativamente maggiore rispetto a un programma leggero (p < 0,0001) o a nessun programma di esercizi (p < 0,0001).

Sono stati eseguiti controlli di normalità e test di Levene per verificare che le ipotesi dell’ANOVA fossero soddisfatte.