Jak wykonać dwukierunkową anova w r

Dwukierunkową ANOVA („analizę wariancji”) stosuje się w celu ustalenia, czy istnieje statystycznie istotna różnica między średnimi z trzech lub więcej niezależnych grup, które zostały podzielone na dwa czynniki.

W tym samouczku wyjaśniono, jak wykonać dwukierunkową analizę ANOVA w języku R.

Przykład: dwukierunkowa ANOVA w R

Załóżmy, że chcemy ustalić, czy intensywność ćwiczeń i płeć wpływają na utratę wagi. W tym przypadku bierzemy pod uwagę dwa czynniki: aktywność fizyczną i płeć , a zmienną odpowiedzi jest utrata masy ciała mierzona w funtach.

Możemy przeprowadzić dwuczynnikową analizę ANOVA, aby określić, czy ćwiczenia fizyczne i płeć wpływają na utratę wagi oraz czy istnieje interakcja między ćwiczeniami fizycznymi a płcią na utratę wagi.

Rekrutujemy 30 mężczyzn i 30 kobiet do udziału w eksperymencie, w którym losowo przydzielamy po 10 osób do wykonywania przez miesiąc programu bez ćwiczeń, lekkich lub intensywnych ćwiczeń.

Poniższy kod tworzy ramkę danych, z którą będziemy pracować:

 #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

Przeglądaj dane

Jeszcze przed dopasowaniem dwukierunkowego modelu ANOVA możemy lepiej zrozumieć dane, znajdując średnią i odchylenie standardowe utraty masy ciała dla każdej z sześciu grup terapeutycznych przy użyciu pakietu 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 

Możemy również utworzyć wykres pudełkowy dla każdej z sześciu grup terapeutycznych, aby zwizualizować rozkład utraty wagi dla każdej grupy:

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

Od razu widać, że dwie grupy, które uczestniczyły w intensywnych ćwiczeniach, wydają się mieć wyższe wartości utraty wagi. Widzimy również, że mężczyźni mają zwykle wyższe wartości utraty wagi niż kobiety, zarówno w grupach intensywnych , jak i lekkich ćwiczeń.

Następnie dopasujemy dwukierunkowy model ANOVA do naszych danych, aby sprawdzić, czy te różnice wizualne są rzeczywiście istotne statystycznie.

Dopasowanie dwukierunkowego modelu ANOVA

Ogólna składnia dopasowywania dwukierunkowego modelu ANOVA w R jest następująca:

aov(zmienna odpowiedzi ~zmienna_predyktora1 *zmienna_predyktora2, dane = zbiór danych)

Należy zauważyć, że * między dwiema zmiennymi predykcyjnymi wskazuje, że chcemy również przetestować efekt interakcji między dwiema zmiennymi predykcyjnymi.

W naszym przykładzie możemy użyć poniższego kodu, aby dopasować dwukierunkowy model ANOVA, używając Weight_loss jako zmiennej odpowiedzi oraz płci i ćwiczeń jako dwóch zmiennych predykcyjnych.

Możemy następnie użyć funkcji podsumowania(), aby wyświetlić wynik naszego modelu:

 #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

Z wyników modelu widzimy, że płeć , aktywność fizyczna i interakcja między tymi dwiema zmiennymi są statystycznie istotne na poziomie istotności 0,05.

Sprawdzanie założeń modelu

Zanim przejdziemy dalej, musimy sprawdzić, czy założenia naszego modelu są spełnione, aby wyniki naszego modelu były wiarygodne. W szczególności dwukierunkowa ANOVA zakłada:

1. Niezależność – obserwacje każdej grupy muszą być od siebie niezależne. Ponieważ zastosowaliśmy projekt losowy , założenie to powinno być spełnione, więc nie musimy się tym zbytnio przejmować.

2. Normalność – zmienna zależna powinna mieć rozkład w przybliżeniu normalny dla każdej kombinacji grup obu czynników.

Jednym ze sposobów sprawdzenia tego założenia jest utworzenie histogramu reszt modelu. Jeżeli reszty mają rozkład w przybliżeniu normalny, założenie to powinno być spełnione.

 #define model residuals
reside <- model$residuals

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

Reszty mają w przybliżeniu rozkład normalny, więc możemy założyć, że założenie normalności jest spełnione.

3. Równa wariancja – wariancje dla każdej grupy są równe lub w przybliżeniu równe.

Jednym ze sposobów sprawdzenia tego założenia jest wykonanie testu Levene’a na równość wariancji przy użyciu pakietu samochodu :

 #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  

Ponieważ wartość p testu jest większa niż nasz poziom istotności wynoszący 0,05, możemy założyć, że spełnione jest nasze założenie o równości wariancji między grupami.

Przeanalizuj różnice w leczeniu

Po sprawdzeniu, czy założenia modelu zostały spełnione, możemy przeprowadzić test post hoc , aby dokładnie określić, które grupy leczenia różnią się od siebie.

W naszym teście post hoc użyjemy funkcji TukeyHSD() do wykonania testu Tukeya dla wielokrotnych porównań:

 #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

Wartość p wskazuje, czy istnieje statystycznie istotna różnica pomiędzy każdą grupą.

Na przykład w ostatnim wierszu powyżej widzimy, że grupa mężczyzn niećwiczących nie doświadczyła statystycznie istotnej różnicy w utracie wagi w porównaniu z grupą kobiet niećwiczących (wartość p: 0,990364).

Możemy również wizualizować 95% przedziały ufności wynikające z testu Tukeya za pomocą funkcji plot() w 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)

Raportowanie wyników dwuczynnikowej analizy ANOVA

Na koniec możemy zgłosić wyniki dwuczynnikowej analizy ANOVA w sposób podsumowujący wyniki:

Przeprowadzono dwuczynnikową analizę ANOVA w celu zbadania wpływu płci ( mężczyzna, kobieta) i programu ćwiczeń (brak, lekkie, intensywne) na utratę wagi (mierzoną w funtach). Stwierdzono statystycznie istotną interakcję pomiędzy wpływem płci i ćwiczeń na utratę masy ciała (F(2, 54) = 4,615, p = 0,0141). Wykonano post-hoc testy HSD Tukeya.

W przypadku mężczyzn intensywny program ćwiczeń spowodował znacznie większą utratę wagi niż program lekki (p < 0,0001) lub brak programu ćwiczeń (p < 0,0001). Dodatkowo u mężczyzn dieta lekkostrawna powodowała istotnie większą utratę masy ciała niż brak ćwiczeń fizycznych (p < 0,0001).

W przypadku kobiet intensywny program ćwiczeń spowodował znacznie większą utratę wagi niż program lekki (p < 0,0001) lub brak programu ćwiczeń (p < 0,0001).

Przeprowadzono kontrole normalności i test Levene’a, aby sprawdzić, czy spełniono założenia analizy ANOVA.