Jak wykonać jednokierunkową anova w r

Jednoczynnikową ANOVA stosuje się do określenia, czy istnieje statystycznie istotna różnica między średnimi z trzech lub więcej niezależnych grup.

Ten typ testu nazywa się jednokierunkową ANOVA, ponieważ analizujemy wpływ zmiennej predykcyjnej na zmienną odpowiedzi.

Uwaga : gdybyśmy zamiast tego byli zainteresowani wpływem dwóch zmiennych predykcyjnych na zmienną odpowiedzi, moglibyśmy przeprowadzić dwuczynnikową analizę ANOVA .

Jak wykonać jednokierunkową ANOVA w R

Poniższy przykład ilustruje sposób przeprowadzenia jednokierunkowej analizy ANOVA w R.

Tło

Załóżmy, że chcemy ustalić, czy trzy różne programy ćwiczeń mają różny wpływ na utratę wagi. Zmienną predykcyjną, którą badamy, jest program ćwiczeń , azmienną odpowiedzi jest utrata masy ciała mierzona w funtach.

Możemy wykonać jednoczynnikową analizę ANOVA, aby określić, czy istnieje statystycznie istotna różnica pomiędzy utratą masy ciała wynikającą z trzech programów.

Rekrutujemy 90 osób do udziału w eksperymencie, w ramach którego losowo przydzielamy 30 osób do udziału w Programie A, Programie B lub Programie C przez miesiąc.

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

 #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

Pierwsza kolumna ramki danych pokazuje program, w którym dana osoba uczestniczyła przez miesiąc, a druga kolumna pokazuje całkowitą utratę wagi, której doświadczyła osoba na koniec programu, mierzoną w funtach.

Przeglądaj dane

Jeszcze przed dopasowaniem jednokierunkowego modelu ANOVA możemy lepiej zrozumieć dane, znajdując średnią i odchylenie standardowe utraty wagi dla każdego z trzech programów 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 (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  

Możemy również utworzyć wykres pudełkowy dla każdego z trzech programów, aby zwizualizować rozkład utraty wagi dla każdego programu:

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

Z tych wykresów pudełkowych widzimy, że średnia utrata masy ciała jest najwyższa w przypadku uczestników Programu C, a średnia utrata masy ciała jest najniższa w przypadku uczestników Programu A.

Widzimy również, że odchylenie standardowe („długość” wykresu pudełkowego) utraty wagi jest nieco wyższe w programie C w porównaniu z dwoma pozostałymi programami.

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

Jednokierunkowe dopasowanie modelu ANOVA

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

aov(zmienna odpowiedzi ~ zmienna_predykcyjna, dane = zbiór danych)

W naszym przykładzie możemy użyć poniższego kodu, aby dopasować jednokierunkowy model ANOVA, używając Weight_loss jako zmiennej odpowiedzi i programu jako zmiennej predykcyjnej. Możemy następnie użyć funkcji podsumowania(), aby wyświetlić wynik naszego modelu:

 #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

Z wyników modelu wynika, że program zmiennych predykcyjnych jest istotny statystycznie na poziomie istotności 0,05.

Innymi słowy, istnieje statystycznie istotna różnica pomiędzy średnim ubytkiem masy ciała wynikającym z trzech programów.

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 jednokierunkowa ANOVA zakłada:

1. Niezależność – obserwacje każdej grupy muszą być od siebie niezależne. Ponieważ zastosowaliśmy projekt losowy (czyli losowo przypisaliśmy uczestników do programów ćwiczeń), to założenie powinno być spełnione, abyśmy nie musieli się tym zbytnio przejmować.

2. Normalność – zmienna zależna powinna mieć w przybliżeniu rozkład normalny dla każdego poziomu zmiennej predykcyjnej.

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

Jednym ze sposobów sprawdzenia założeń o normalności i równej wariancji jest użycie funkcji plot() , która generuje cztery wykresy sprawdzające model. W szczególności interesują nas dwie następujące działki:

• Pozostałości vs. dopasowane – ten wykres przedstawia zależność pomiędzy wartościami resztowymi i dopasowanymi. Możemy użyć tego wykresu, aby z grubsza ocenić, czy wariancja między grupami jest w przybliżeniu równa, czy nie.
• Wykres QQ – ten wykres przedstawia reszty standaryzowane względem teoretycznych kwantyli. Możemy użyć tego wykresu, aby z grubsza ocenić, czy założenie normalności jest spełnione.

Do utworzenia wykresów sprawdzających model można użyć poniższego kodu:

 plot(model)

Powyższy wykres QQ pozwala nam zweryfikować założenie normalności. W idealnym przypadku reszty standaryzowane leżałyby wzdłuż prostej ukośnej linii wykresu. Jednakże na powyższym wykresie widzimy, że reszty odbiegają nieco od linii w kierunku początku i końca. Oznacza to, że nasze założenie o normalności może zostać naruszone.

Pozostałości vs. Skorygowany powyższy wykres pozwala nam zweryfikować nasze założenie o równych wariancjach. Idealnie byłoby, gdyby reszty były równomiernie rozłożone na każdym poziomie dopasowanych wartości.

Widzimy, że reszty są znacznie bardziej rozłożone w przypadku wyższych dopasowanych wartości, co wskazuje, że nasze założenie o równości wariancji może zostać naruszone.

Aby formalnie przetestować równe wariancje, moglibyśmy wykonać test Levene’a, korzystając z pakietu samochodu :

 #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

Wartość p testu wynosi 0,01862 . Jeśli zastosujemy poziom istotności 0,05, odrzucimy hipotezę zerową, że wariancje są równe we wszystkich trzech programach. Jeśli jednak zastosujemy poziom istotności 0,01, nie odrzucimy hipotezy zerowej.

Chociaż moglibyśmy podjąć próbę przekształcenia danych, aby upewnić się, że nasze założenia dotyczące normalności i równości wariancji zostaną spełnione, na razie nie będziemy się tym zbytnio przejmować.

Przeanalizuj różnice w leczeniu

Po sprawdzeniu, czy założenia modelu są spełnione (lub w rozsądnym stopniu 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 ~ 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

Wartość p wskazuje, czy istnieje statystycznie istotna różnica pomiędzy każdym programem. Wyniki pokazują, że istnieje statystycznie istotna różnica pomiędzy średnią utratą masy ciała w każdym programie na poziomie istotności 0,05.

Możemy również wizualizować 95% przedziały ufności wynikające z testu Tukeya za pomocą funkcji plot(TukeyHSD()) w R:

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

Wyniki przedziałów ufności są zgodne z wynikami testów hipotez.

W szczególności widać, że żaden z przedziałów ufności dla średniej utraty wagi pomiędzy programami nie zawiera wartości zero , co wskazuje, że istnieje statystycznie istotna różnica w średniej utracie wagi pomiędzy trzema programami.

Jest to zgodne z tym, że wszystkie wartości p dla naszych testów hipotez są mniejsze niż 0,05.

Raportowanie wyników jednokierunkowej analizy ANOVA

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

W celu zbadania efektów programu ćwiczeń przeprowadzono jednoczynnikową analizę ANOVA   utraty wagi (mierzonej w funtach). Stwierdzono statystycznie istotną różnicę pomiędzy wpływem trzech programów na utratę wagi (F(2, 87) = 30,83, p = 7,55e-11). Przeprowadzono post-hoc testy HSD Tukeya.

Średni ubytek masy ciała uczestników programu C jest istotnie większy niż średni ubytek masy ciała uczestników programu B (p < 0,0001).

Średni ubytek masy ciała uczestników programu C jest istotnie większy niż średni ubytek masy ciała uczestników programu A (p < 0,0001).

Dodatkowo średnia utrata masy ciała uczestników programu B była istotnie większa niż średnia utrata masy ciała uczestników programu A (p = 0,01).

Dodatkowe zasoby

Poniższe samouczki zawierają dodatkowe informacje na temat jednokierunkowej analizy ANOVA:

Wprowadzenie do jednokierunkowej ANOVA
Przewodnik po korzystaniu z testów post-hoc z ANOVA
Kompletny przewodnik: Jak zgłaszać wyniki ANOVA