Eenrichtings-anova uitvoeren in r

Een eenrichtings-ANOVA wordt gebruikt om te bepalen of er al dan niet een statistisch significant verschil bestaat tussen de gemiddelden van drie of meer onafhankelijke groepen.

Dit type test wordt eenrichtings- ANOVA genoemd omdat we de impact van een voorspellende variabele op een responsvariabele analyseren.

Opmerking : als we in plaats daarvan geïnteresseerd zouden zijn in de impact van twee voorspellende variabelen op een responsvariabele, zouden we een tweerichtings-ANOVA kunnen uitvoeren.

Eenrichtings-ANOVA uitvoeren in R

Het volgende voorbeeld illustreert hoe u een eenrichtings-ANOVA in R uitvoert.

Achtergrond

Stel dat we willen bepalen of drie verschillende oefenprogramma’s een verschillende impact hebben op gewichtsverlies. De voorspellende variabele die we bestuderen is het oefenprogramma en deresponsvariabele is gewichtsverlies, gemeten in kilo’s.

We kunnen een one-way ANOVA uitvoeren om te bepalen of er een statistisch significant verschil is tussen het gewichtsverlies als gevolg van de drie programma’s.

We werven 90 mensen om deel te nemen aan een experiment waarbij we willekeurig 30 mensen toewijzen om een maand lang Programma A, Programma B of Programma C te volgen.

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

 #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

De eerste kolom van het dataframe toont het programma waaraan de persoon een maand lang heeft deelgenomen en de tweede kolom toont het totale gewichtsverlies dat de persoon aan het einde van het programma heeft ervaren, gemeten in kilo’s.

Verken de gegevens

Voordat we zelfs maar het eenrichtings-ANOVA-model passen, kunnen we de gegevens beter begrijpen door de gemiddelde en standaardafwijking van het gewichtsverlies te vinden voor elk van de drie programma’s 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 (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  

We kunnen ook een boxplot maken voor elk van de drie programma’s om de verdeling van het gewichtsverlies voor elk programma te visualiseren:

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

Uit deze boxplots kunnen we zien dat het gemiddelde gewichtsverlies het hoogst is voor deelnemers aan Programma C en het gemiddelde gewichtsverlies het laagst is voor deelnemers aan Programma A.

We kunnen ook zien dat de standaardafwijking (de „lengte“ van de boxplot) voor gewichtsverlies iets hoger is in programma C vergeleken met de andere twee programma’s.

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

Eenrichtings-ANOVA-modelfitting

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

aov(responsvariabele ~ voorspellervariabele, data = dataset)

In ons voorbeeld kunnen we de volgende code gebruiken om het eenrichtings-ANOVA-model aan te passen, waarbij we gewicht_verlies gebruiken als de responsvariabele en programma als de voorspellende variabele. We kunnen vervolgens de functie summary() gebruiken om het resultaat van ons model weer te geven:

 #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

Uit de modelresultaten kunnen we zien dat het programma van voorspellende variabelen statistisch significant is op het significantieniveau van 0,05.

Met andere woorden: er is een statistisch significant verschil tussen het gemiddelde gewichtsverlies als gevolg van de drie programma’s.

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 eenrichtings-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 (dat wil zeggen, we hebben deelnemers willekeurig aan de oefenprogramma’s toegewezen), moet aan deze veronderstelling worden voldaan, zodat we ons er niet al te veel zorgen over hoeven te maken.

2. Normaliteit – de afhankelijke variabele moet een ongeveer normale verdeling hebben voor elk niveau van de voorspellende variabele.

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

Eén manier om de aannames van normaliteit en gelijke variantie te controleren, is door de functie plot() te gebruiken, die vier modelcontroleplots produceert. In het bijzonder zijn wij bijzonder geïnteresseerd in de volgende twee percelen:

• Residuen vs. aangepast – deze grafiek toont de relatie tussen residuen en aangepaste waarden. We kunnen deze grafiek gebruiken om grofweg te beoordelen of de variantie tussen groepen ongeveer gelijk is of niet.
• QQ-plot – deze plot toont de gestandaardiseerde residuen ten opzichte van de theoretische kwantielen. We kunnen deze grafiek gebruiken om grofweg te beoordelen of aan de normaliteitsaanname wordt voldaan.

De volgende code kan worden gebruikt om deze modelcontroleplots te maken:

 plot(model)

Met de bovenstaande QQ-grafiek kunnen we de normaliteitsaanname verifiëren. Idealiter zouden de gestandaardiseerde residuen langs de rechte diagonale lijn van de grafiek liggen. In de bovenstaande grafiek kunnen we echter zien dat de residuen naar het begin en het einde enigszins afwijken van de lijn. Dit geeft aan dat onze normaliteitsaanname mogelijk wordt geschonden.

De residuen vs. Met de aangepaste grafiek hierboven kunnen we onze aanname van gelijke varianties verifiëren. Idealiter zouden we willen dat de residuen gelijk verdeeld worden voor elk niveau van de gefitte waarden.

We kunnen zien dat de residuen veel meer verspreid zijn voor de hoger aangepaste waarden, wat aangeeft dat onze aanname van gelijkheid van varianties mogelijk geschonden is.

Om formeel op gelijke varianties te testen, kunnen we de Levene-test uitvoeren met behulp van het autopakket :

 #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

De p-waarde van de test is 0,01862 . Als we een significantieniveau van 0,05 gebruiken, zouden we de nulhypothese verwerpen dat de varianties voor de drie programma’s gelijk zijn. Als we echter een significantieniveau van 0,01 hanteren, zullen we de nulhypothese niet verwerpen.

Hoewel we zouden kunnen proberen de gegevens te transformeren om ervoor te zorgen dat aan onze aannames van normaliteit en gelijkheid van varianties wordt voldaan, zullen we ons daar voorlopig niet al te veel zorgen over maken.

Analyseer behandelingsverschillen

Zodra we hebben geverifieerd dat aan de modelaannames wordt voldaan (of redelijkerwijs wordt 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 ~ 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

De p-waarde geeft aan of er al dan niet een statistisch significant verschil bestaat tussen elk programma. De resultaten laten zien dat er een statistisch significant verschil is tussen het gemiddelde gewichtsverlies van elk programma op het significantieniveau van 0,05.

We kunnen ook de 95% betrouwbaarheidsintervallen visualiseren die het resultaat zijn van de Tukey-test met behulp van de plot(TukeyHSD())- functie in R:

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

De resultaten van de betrouwbaarheidsintervallen komen overeen met de resultaten van de hypothesetoetsen.

We kunnen met name zien dat geen van de betrouwbaarheidsintervallen voor het gemiddelde gewichtsverlies tussen de programma’s de waarde nul bevat, wat aangeeft dat er een statistisch significant verschil is in het gemiddelde gewichtsverlies tussen de drie programma’s.

Dit komt overeen met het feit dat alle p-waarden voor onze hypothesetoetsen kleiner zijn dan 0,05.

Rapporteren van One-Way ANOVA-resultaten

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

Er werd een one-way ANOVA uitgevoerd om de effecten van het oefenprogramma te onderzoeken   op gewichtsverlies (gemeten in kilo’s). Er was een statistisch significant verschil tussen de effecten van de drie programma’s op gewichtsverlies (F(2, 87) = 30,83, p = 7,55e-11). Post-hoc werden de HSD-tests van Tukey uitgevoerd.

Het gemiddelde gewichtsverlies van deelnemers aan programma C is significant groter dan het gemiddelde gewichtsverlies van deelnemers aan programma B (p < 0,0001).

Het gemiddelde gewichtsverlies van deelnemers aan programma C is significant groter dan het gemiddelde gewichtsverlies van deelnemers aan programma A (p < 0,0001).

Bovendien was het gemiddelde gewichtsverlies van deelnemers aan programma B significant groter dan het gemiddelde gewichtsverlies van deelnemers aan programma A (p = 0,01).

Aanvullende bronnen

De volgende zelfstudies bieden aanvullende informatie over eenrichtings-ANOVA’s:

Een inleiding tot One-Way ANOVA
Een handleiding voor het gebruik van post-hoctesten met ANOVA
De complete gids: ANOVA-resultaten rapporteren