Wprowadzenie
ANOVA z powtarzanymi pomiarami przyjmuje następujące założenia dotyczące danych:
Brak znaczących wartości odstających w jakiejkolwiek komórce projektu. Można to sprawdzić wizualizując dane za pomocą metody box plot oraz za pomocą funkcji
identify_outliers()[pakiet rstatix].Normalność: zmienna wynikowa (lub zależna) powinna być w przybliżeniu normalnie rozłożona w każdej komórce projektu. Może to być sprawdzone przy użyciu testu normalności Shapiro-Wilka (
shapiro_test()[rstatix]) lub przez wizualną inspekcję przy użyciu QQ plot (gggqqplot()[pakiet ggpubr]).Założenie o sferyczności: wariancja różnic między grupami powinna być równa. Można to sprawdzić za pomocą testu sferyczności Mauchly’ego, który jest automatycznie podawany przy użyciu funkcji R
anova_test()[pakiet rstatix].
Zauważ, że jeśli powyższe założenia nie są spełnione istnieje nieparametryczna alternatywa (test Friedmana) dla jednokierunkowej ANOVY powtarzanych pomiarów!
Niestety, nie istnieją nieparametryczne alternatywy dla dwukierunkowej i trójkierunkowej ANOVY powtarzanych pomiarów! Tak więc, w sytuacji, gdy założenia nie są spełnione, można rozważyć przeprowadzenie dwukierunkowego/trzykierunkowego testu ANOVA dla danych przekształconych i nieprzekształconych, aby sprawdzić, czy istnieją jakieś znaczące różnice (np. za pomocą logarytmu).
Jeśli oba testy prowadzą do tych samych wniosków, można nie decydować się na przekształcenie zmiennej wynikowej i kontynuować dwu-/trzy-kierunkową ANOVA dla danych oryginalnych.
Możliwe jest również wykonanie odpornego testu ANOVA przy użyciu pakietu WRS2 w R.
Anova dla prób zależnych w R
Kluczowe funkcje:
anova_test()
[pakiet rstatix], wrapper wokółcar::Anova()` dla ułatwienia obliczeń ANOVA z powtarzanymi miarami. Kluczowe argumenty do wykonania ANOVA z powtarzanymi miarami:data: ramka danychdv: (numeryczna) badana zmienna zależnawid: identyfikator rekorduwithin: within-subjects czynnik, zmienna grupująca
get_anova_table()[rstatix package]. Funkcja ta wyodrębnia tabelę ANOVA z funkcjianova_test(). Zwraca tabelę ANOVA, która jest automatycznie skorygowana o ewentualne odchylenia od założenia o sferyczności. Domyślnie automatycznie stosowana jest korekta sferyczności Greenhouse’a-Geissera tylko do czynników wewnątrzprzedmiotowych naruszających założenie sferyczności (tzn. wartość p testu Mauchly’ego jest istotna, p <= 0.05).
Anova 1-czynnikowa dla prób zależnych
W zbiorze “selfesteem” zawarto wyniki samooceny 10 osób w trzech punktach czasowych podczas stosowania określonej diety w celu określenia, czy ich samoocena uległa poprawie.
data("selfesteem", package = "datarium")
head(selfesteem, 3)## # A tibble: 3 x 4
## id t1 t2 t3
## <int> <dbl> <dbl> <dbl>
## 1 1 4.01 5.18 7.11
## 2 2 2.56 6.91 6.31
## 3 3 3.24 4.44 9.78Jednoczynnikowa ANOVA może być użyta do określenia, czy średnie wyniki samooceny różnią się znacząco pomiędzy trzema punktami czasowymi. Przekształćmy więc tę ramkę danych na długi format:
selfesteem <- selfesteem %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
convert_as_factor(id, time)
head(selfesteem, 3)## # A tibble: 3 x 3
## id time score
## <fct> <fct> <dbl>
## 1 1 t1 4.01
## 2 2 t1 2.56
## 3 3 t1 3.24Statystyki opisowe
Tabela statystyk opisowych wg czasu:
depression <- depression %>%
gather(key="time", value="score", t0, t1, t2, t3) %>%
convert_as_factor(id, time, treatment)
head(depression, 3)## # A tibble: 3 x 4
## id treatment time score
## <fct> <fct> <fct> <dbl>
## 1 1 ctr t0 296
## 2 2 ctr t0 376
## 3 3 ctr t0 309Wykresy ramkowe:
ggboxplot(depression, x = "time", y = "score", color = "treatment")
Założenia
Zbadajmy, czy w zbiorze mamy obserwacje odstające:
depression %>%
group_by(time, treatment) %>%
identify_outliers(score)## # A tibble: 4 x 6
## treatment time id score is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <lgl> <lgl>
## 1 ctr t0 5 150 TRUE FALSE
## 2 ctr t0 8 447 TRUE FALSE
## 3 treated t0 18 138 TRUE FALSE
## 4 treated t0 22 150 TRUE FALSESprawdźmy założenie normalności zmiennej badanej testem Shapiro-Wilka wg czasu:
depression %>%
group_by(time, treatment) %>%
shapiro_test(score)## # A tibble: 8 x 5
## treatment time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 ctr t0 score 0.945 0.572
## 2 treated t0 score 0.836 0.0248
## 3 ctr t1 score 0.918 0.268
## 4 treated t1 score 0.953 0.680
## 5 ctr t2 score 0.951 0.652
## 6 treated t2 score 0.957 0.734
## 7 ctr t3 score 0.873 0.0706
## 8 treated t3 score 0.932 0.398Anova
test <- anova_test(data = depression, formula = score~treatment*time, wid = id, within = time)## Coefficient covariances computed by hccm()test## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 treatment 1 88 43.423 3.13e-09 * 0.330
## 2 time 3 88 19.713 7.31e-10 * 0.402
## 3 treatment:time 3 88 3.755 1.40e-02 * 0.113Testy post-hoc
depression %>%
group_by(treatment) %>%
emmeans_test(score~time)## # A tibble: 12 x 10
## treatment term .y. group1 group2 df statistic p p.adj
## * <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 ctr time score t0 t1 88 1.60 0.113 0.676
## 2 ctr time score t0 t2 88 2.98 0.00375 0.0225
## 3 ctr time score t0 t3 88 3.60 0.000520 0.00312
## 4 ctr time score t1 t2 88 1.38 0.172 1
## 5 ctr time score t1 t3 88 2.00 0.0485 0.291
## 6 ctr time score t2 t3 88 0.625 0.534 1
## 7 treated time score t0 t1 88 6.25 0.0000000142 0.0000000855
## 8 treated time score t0 t2 88 6.05 0.0000000338 0.000000203
## 9 treated time score t0 t3 88 5.84 0.0000000854 0.000000513
## 10 treated time score t1 t2 88 -0.196 0.845 1
## 11 treated time score t1 t3 88 -0.409 0.683 1
## 12 treated time score t2 t3 88 -0.213 0.832 1
## # ... with 1 more variable: p.adj.signif <chr>Wnioski i podsumowanie
ggwithinstats(data = depression, x = time, y = score, grouping.var = treatment, type = "P")