Hoe tweeweg-anova uit te voeren in r

Een tweewegs-ANOVA (“variantieanalyse”) wordt gebruikt om te bepalen of er al dan niet een statistisch significant verschil bestaat tussen de gemiddelden van drie of meer onafhankelijke groepen die over twee factoren zijn verdeeld.

In deze tutorial wordt uitgelegd hoe u een tweerichtings-ANOVA uitvoert in R.

Voorbeeld: Tweeweg-ANOVA in R

Laten we zeggen dat we willen bepalen of de trainingsintensiteit en het geslacht invloed hebben op het gewichtsverlies. In dit geval zijn de twee factoren waar we naar kijken lichaamsbeweging en geslacht , en is de responsvariabele gewichtsverlies, gemeten in kilo’s.

We kunnen een tweerichtings-ANOVA uitvoeren om te bepalen of lichaamsbeweging en geslacht invloed hebben op gewichtsverlies en om te bepalen of er een interactie bestaat tussen lichaamsbeweging en geslacht op gewichtsverlies.

We rekruteren 30 mannen en 30 vrouwen om deel te nemen aan een experiment waarbij we willekeurig 10 van elk een maand lang een programma zonder lichaamsbeweging, lichte lichaamsbeweging of intensief bewegen toewijzen.

De volgende code creëert het dataframe waarmee we zullen werken:

 #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

Verken de gegevens

Voordat we zelfs maar het tweerichtings-ANOVA-model passen, kunnen we de gegevens beter begrijpen door de gemiddelde en standaardafwijking van het gewichtsverlies te vinden voor elk van de zes behandelingsgroepen met behulp van het dplyr- pakket:

 #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 

We kunnen ook een boxplot maken voor elk van de zes behandelingsgroepen om de verdeling van het gewichtsverlies voor elke groep te visualiseren:

 #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
)

We kunnen meteen zien dat de twee groepen die aan intensieve training deelnamen hogere gewichtsverlieswaarden lijken te hebben. We kunnen ook zien dat mannen doorgaans hogere gewichtsverlieswaarden hebben dan vrouwen in zowel de intensieve als de lichte trainingsgroepen.

Vervolgens passen we het tweerichtings-ANOVA-model aan onze gegevens toe om te zien of deze visuele verschillen daadwerkelijk statistisch significant zijn.

Passend bij het tweeweg ANOVA-model

De algemene syntaxis voor het aanpassen van een tweerichtings-ANOVA-model in R is:

aov(responsvariabele ~predictor_variable1 *predictor_variable2, data = dataset)

Merk op dat de * tussen de twee voorspellende variabelen aangeeft dat we ook willen testen op een interactie-effect tussen de twee voorspellende variabelen.

In ons voorbeeld kunnen we de volgende code gebruiken om het tweerichtings-ANOVA-model aan te passen, waarbij we gewichtsverlies gebruiken als de responsvariabele en geslacht en beweging als de twee voorspellende variabelen.

We kunnen vervolgens de functie summary() gebruiken om het resultaat van ons model weer te geven:

 #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

Uit de modelresultaten kunnen we zien dat geslacht , lichaamsbeweging en de interactie tussen de twee variabelen allemaal statistisch significant zijn op het significantieniveau van 0,05.

Het controleren van modelaannames

Voordat we verder gaan, moeten we verifiëren dat aan de aannames van ons model wordt voldaan, zodat onze modelresultaten betrouwbaar zijn. Een tweerichtings-ANOVA gaat met name uit van:

1. Onafhankelijkheid – de observaties van elke groep moeten onafhankelijk van elkaar zijn. Omdat we een gerandomiseerd ontwerp hebben gebruikt , zou aan deze veronderstelling moeten worden voldaan, dus we hoeven ons er niet al te veel zorgen over te maken.

2. Normaliteit – de afhankelijke variabele moet een ongeveer normale verdeling hebben voor elke combinatie van groepen van de twee factoren.

Eén manier om deze aanname te testen is door een histogram van de modelresiduen te maken. Als de residuen bij benadering normaal verdeeld zijn, moet aan deze aanname worden voldaan.

 #define model residuals
reside <- model$residuals

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

De residuen zijn bij benadering normaal verdeeld, dus we kunnen aannemen dat aan de normaliteitsaanname is voldaan.

3. Gelijke variantie – de varianties voor elke groep zijn gelijk of ongeveer gelijk.

Eén manier om deze aanname te controleren is door een Levene-test uit te voeren voor gelijkheid van varianties met behulp van het autopakket :

 #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  

Omdat de p-waarde van de test groter is dan ons significantieniveau van 0,05, kunnen we aannemen dat aan onze aanname van gelijkheid van varianties tussen groepen wordt voldaan.

Analyseer behandelingsverschillen

Zodra we hebben geverifieerd dat aan de aannames van het model is voldaan, kunnen we een post-hoctest uitvoeren om precies te bepalen welke behandelingsgroepen van elkaar verschillen.

Voor onze post-hoc-test zullen we de functie TukeyHSD() gebruiken om de Tukey-test uit te voeren voor meerdere vergelijkingen:

 #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

De p-waarde geeft aan of er al dan niet een statistisch significant verschil bestaat tussen elke groep.

In de laatste rij hierboven zien we bijvoorbeeld dat de groep mannen die niet trainde geen statistisch significant verschil in gewichtsverlies ervoer vergeleken met de groep vrouwen die niet trainde (p-waarde: 0,990364).

We kunnen ook de 95% betrouwbaarheidsintervallen visualiseren die het resultaat zijn van de Tukey-test met behulp van de plot()- functie 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)

Tweerichtings-ANOVA-resultaten rapporteren

Ten slotte kunnen we de resultaten van de tweeweg-ANOVA rapporteren op een manier die de resultaten samenvat:

Er werd een tweeweg-ANOVA uitgevoerd om de effecten van geslacht ( man, vrouw) en trainingsprogramma (geen, licht, intens) op gewichtsverlies te onderzoeken (gemeten in kilo’s). Er was een statistisch significante interactie tussen de effecten van geslacht en lichaamsbeweging op gewichtsverlies (F(2, 54) = 4,615, p = 0,0141). Post-hoc werden de HSD-tests van Tukey uitgevoerd.

Bij mannen resulteerde een intensief oefenprogramma in significant meer gewichtsverlies dan een licht programma (p < 0,0001) of geen oefenprogramma (p < 0,0001). Bovendien resulteerde een licht dieet bij mannen in een significant groter gewichtsverlies dan geen trainingsregime (p < 0,0001).

Bij vrouwen resulteerde een intensief oefenprogramma in significant meer gewichtsverlies dan een licht programma (p < 0,0001) of geen oefenprogramma (p < 0,0001).

Normaliteitscontroles en de Levene-test werden uitgevoerd om te verifiëren dat aan de aannames van de ANOVA werd voldaan.