Anova dla prób zależnych

2022-01-17

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 danych

    • dv: (numeryczna) badana zmienna zależna

    • wid: identyfikator rekordu

    • within: within-subjects czynnik, zmienna grupująca

  • get_anova_table() [rstatix package]. Funkcja ta wyodrębnia tabelę ANOVA z funkcji anova_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 2-czynnikowa dla prób zależnych

W zbiorze “depression” zawarto wyniki poziom stanu depresjnego u 24 osób w czterech punktach czasowych. połowa badanych otrzymała lek, a pozostali są poddani kontroli. Badanie ma na celu orzec, czy na wynik depresji wpływa czas i leczenie.

data("depression", package = "datarium")
head(depression, 3)
## # A tibble: 3 x 6
##   id    treatment    t0    t1    t2    t3
##   <fct> <fct>     <dbl> <dbl> <dbl> <dbl>
## 1 1     ctr         296   175   187   242
## 2 2     ctr         376   329   236   126
## 3 3     ctr         309   238   150   173

Przekształcenie ramki w potzrebną do wykonania obliczeń ramkę danych na długi format:

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      309

Statystyki opisowe

Wykresy ramkowe:

ggboxplot(depression, x = "time", y = "score", color = "treatment", palette = "jco")

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       FALSE

Dane depression posiadają 4 zmienne odstają, które nie odstają na szczęście ekstremalnie. Ten sam wniosek, można wydedukować dzięki zamieszczonemu wyżej wykresowi ramka-wąs.

Sprawdź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.398

Tylko dla jednej zmiennej (osoba leczona w czasie 0) hipoteza zerowa zostaje odrzucona. Pozostałe podlegają normalności.

Wykres normalności:

ggqqplot(depression, "score") + facet_grid(time~treatment)

Anova

anova <- anova_test(data=depression, formula = score~treatment*time,wid = id, within = time)
## Coefficient covariances computed by hccm()
get_anova_table(anova)
## 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.113

Zgodnie z testem ANOVA na wynik depresji istotnie wpływa czas oraz leczenie.

Testy post-hoc

Co sprawiło, że nastapiła widoczna różnica?

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>

Dla pięciu z dwunastu powiązanych zmiennych różniece w wynikach poziomu stanu depresyjnego okazują sie istotne (p_value<0.5). Są to: 1. osoba pod obserwacjąT0T2 2. osoba pod obserwacjąT0T3 3. osoba leczonaT0T1 4. osoba leczonaT0T2 5. osoba leczonaTOT3

Wykres podsumuwujący

ggwithinstats(data = depression, x = time, y = score, grouping.var = treatment, type = “P”)