Anova 3-czynnikowa
Wprowadzenie
Trójczynnikowa ANOVA jest rozszerzeniem dwukierunkowej ANOVA do oceny, czy istnieje efekt interakcji pomiędzy trzema niezależnymi zmiennymi kategorycznymi (jakościowymi) na ciągłą zmienną wynikową (ilościową).
Wykorzystamy zbiór danych headache [z pakietu datarium], który zawiera ankietę dot. miary bólu w epizodzie migrenowego bólu głowy u 72 uczestników leczonych trzema różnymi metodami. Wśród uczestników jest 36 mężczyzn i 36 kobiet. Mężczyźni i kobiety zostali dalej podzieleni na grupy o niskim lub wysokim ryzyku migreny.
Chcemy zrozumieć, jak każda niezależna zmienna (rodzaj leczenia, ryzyko migreny i płeć) współdziałają w przewidywaniu wyniku bólu.
Statystyki opisowe
Zobaczmy jak wygląda rozkład wyników - poziomu bólu głowy osób o różnej płci, ryzyku migreny i rodzaju leczenia.
headache %>%
group_by(gender, risk, treatment) %>%
get_summary_stats(pain_score, type = "mean_sd")## # A tibble: 12 x 7
## gender risk treatment variable n mean sd
## <fct> <fct> <fct> <chr> <dbl> <dbl> <dbl>
## 1 male high X pain_score 6 92.7 5.12
## 2 male high Y pain_score 6 82.3 5.00
## 3 male high Z pain_score 6 79.7 4.05
## 4 male low X pain_score 6 76.1 3.86
## 5 male low Y pain_score 6 73.1 4.76
## 6 male low Z pain_score 6 74.5 4.89
## 7 female high X pain_score 6 78.9 5.32
## 8 female high Y pain_score 6 81.2 4.62
## 9 female high Z pain_score 6 81.0 3.98
## 10 female low X pain_score 6 74.2 3.69
## 11 female low Y pain_score 6 68.4 4.08
## 12 female low Z pain_score 6 69.8 2.72data(headache)
ggplot(headache, aes(x=pain_score)) +
geom_histogram(bins=10) +
facet_grid(gender~risk~treatment)
Rozkład wyników przedstawiony wkresem pudełkowym.
pudelkowy<-ggboxplot(headache, x="treatment",
y="pain_score",color="risk",palette="jco",facet.by="gender")
pudelkowy
Założenia
Sprawdzamy, czy mamy prawo użyć Anovy parametrycznej.
Obserwacje odstające
Sprawdzamy, czy wśród zmiennych występują obserwacje odstające.
headache %>%
group_by(gender,risk,treatment) %>%
identify_outliers(pain_score)## # A tibble: 4 x 7
## gender risk treatment id pain_score is.outlier is.extreme
## <fct> <fct> <fct> <int> <dbl> <lgl> <lgl>
## 1 female high X 57 68.4 TRUE TRUE
## 2 female high Y 62 73.1 TRUE FALSE
## 3 female high Z 67 75.0 TRUE FALSE
## 4 female high Z 71 87.1 TRUE FALSEWystępują 4 obserwacje odstające, a jedna z nich jest obserwacją odstajacą ekstremalnie.
Normalność
W tym kroku sprawdzamy normalność rozkładu zmiennych
headache %>%
group_by(gender,risk,treatment) %>%
shapiro_test(pain_score)## # A tibble: 12 x 6
## gender risk treatment variable statistic p
## <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 male high X pain_score 0.958 0.808
## 2 male high Y pain_score 0.902 0.384
## 3 male high Z pain_score 0.955 0.784
## 4 male low X pain_score 0.982 0.962
## 5 male low Y pain_score 0.920 0.507
## 6 male low Z pain_score 0.924 0.535
## 7 female high X pain_score 0.714 0.00869
## 8 female high Y pain_score 0.939 0.654
## 9 female high Z pain_score 0.971 0.901
## 10 female low X pain_score 0.933 0.600
## 11 female low Y pain_score 0.927 0.555
## 12 female low Z pain_score 0.958 0.801Z powyższego testu wynika, że odrzucamy hipotezę zerową o normalności rozkładu tylko dla jednej z badanych grup (female,high,X). Dla całej reszty nie ma podstaw, aby odrzucić hipotezę o normalności rozkładów poziomu bólu głowy.
Zobaczmy, jak wyglądają wykresy normalności wg podgrup:
ggqqplot(headache, "pain_score") +
facet_grid(gender~risk~treatment)
Jednorodność wariancji
Założenie jednorodności wariancji
headache %>% levene_test(pain_score ~ gender*risk*treatment)## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 11 60 0.179 0.998P-value > 0.05, a więc stwierdzamy brak istotnej różnicy wariancji (nie możemy odrzucić hipotezy zerowej). Zatem wariancje są jednorodne.
Anova
Przeprowadzamy test Anova.
wyniki <- headache %>% anova_test(pain_score~gender*risk*treatment)## Coefficient covariances computed by hccm()wyniki## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 gender 1 60 16.196 0.000163000000000 * 0.213
## 2 risk 1 60 92.699 0.000000000000088 * 0.607
## 3 treatment 2 60 7.318 0.001000000000000 * 0.196
## 4 gender:risk 1 60 0.141 0.708000000000000 0.002
## 5 gender:treatment 2 60 3.338 0.042000000000000 * 0.100
## 6 risk:treatment 2 60 0.713 0.494000000000000 0.023
## 7 gender:risk:treatment 2 60 7.406 0.001000000000000 * 0.198Z powyższeego testu dowiadujemy się, że występuje statystycznie istotna trójczynnikowa interakcja między płcią, ryzykiem i leczeniem (ponieważ dla #7 F>p).
Testy post-hoc
Sprawdźmy, jak dokładnie wyniki poziomu bólu głowu różnią się wg sposobu leczenia, ryzyka występowania migreny i płci:
wyniki2 <- lm(pain_score ~ gender*risk*treatment, data = headache)
headache %>%
group_by(gender) %>%
anova_test(pain_score ~ risk*treatment, error = wyniki2)## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()## # A tibble: 6 x 8
## gender Effect DFn DFd F p `p<.05` ges
## * <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 male risk 1 60 50.0 0.00000000187 "*" 0.455
## 2 male treatment 2 60 10.2 0.000157 "*" 0.253
## 3 male risk:treatment 2 60 5.25 0.008 "*" 0.149
## 4 female risk 1 60 42.8 0.000000015 "*" 0.416
## 5 female treatment 2 60 0.482 0.62 "" 0.016
## 6 female risk:treatment 2 60 2.87 0.065 "" 0.087Jak widzimy powyżej, dla czterech par średnich różnice w ocenach poziomu bólu głowy okazały się istotne (p<0.05).Dla dwóch par różnice nie są istotne (female:treatment, female:risk:treatment).
##Podsumowanie
Podsumujmy testy Anova za pomocą wykresu ujmującego razem wszystkie poziomy badanych czynników z etykietami i podpisami.
wyniki3 <- lm(pain_score ~ gender*risk*treatment, data = headache)
emmip(wyniki3,gender~risk~treatment,cov.reduce=range,xlab="treatment",
ylab="pain_score")