Wprowadzenie

ANOVA z powtarzanymi środkami 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.

RM Anova in 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\]. Wyodrębnia tabelę ANOVA z funkcji anova_test(). Zwraca ona 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.78

Jednoczynnikowa 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.24

Statystyki opisowe

W poniższej tabeli przedstawiono statystyki opisowe badanych próbek.

to <- selfesteem %>% 
                    select(time,score)

raport <- list(
  "Samoocena" =
    list(
      "Min" = ~ min(score),
      "Max" = ~ max(score),
      "Q1"  = ~ quantile(score, 0.25),
      "Mediana" = ~ round(median(score,0.5),2),
      "Średnia" = ~ round(mean(score),2),
      "Odchylenie std." = ~ round(sd(score),2),
      "Skośność" = ~ round(skew(score),2),
      "Kurtoza" = ~ round(kurtosi(score),2)
    )
)


tabela<-summary_table(to, summaries=raport, by=c("time"))
kbl(tabela,align = "cc", digits=2, caption="Statystyki opisowe ocen samopoczucia w poszczególnych punktach czasowych diety", ) %>%
  kable_classic(full_width=F) %>%
    kable_styling(bootstrap_options=c("striped","hover"))

Na podstawie powyższej tabeli można zauważyć, iż wraz z upływem czasu stosowania diety następował wzrost wartości odnotowywanych statystyk opisowych. Zatem można domniemać, iż samoocena badanych ulegała poprawie.

Powyższe statystyki opisowe zwizualizowano na poniższych wykresach pudełkowych.

wykres1 <- ggboxplot(selfesteem , x = "time", y = "score",
          color = "time", palette = c("darkred", "darksalmon", "deeppink4"),
          order = c("t1", "t2", "t3"), xlab="Punkty czasowe diety", 
          ylab="Ocena samopoczucia", main= "Wizualizacja statystyk opisowych")

wykres2 <- ggplotly(wykres1)
wykres2

Analizując powyższy rysunek można zauważyć, że wraz z upływem czasu następował wzrost poprawy samopoczucia. Jednakże, aby stwierdzić czy był to statystycznie istotny wzrost poprawy samopoczucia należy przeprowadzić test jednoczynnikowej ANOVY dla prób zależnych. Ponadto dla punktu czasowego t2 możemy zidentyfikować outlier.

Założenia

Przed wykonaniem testu ANOVA należy zweryfikować normalność rozkładu analizowanej zmiennej zależnej. W tym celu należy pogrupować dane wg zmiennej “time”, a następnie wykonać test Shapiro-Wilka na normalność rozkładu badanej zmiennej zależnej. W niniejszym teście:

H0: Próba pochodzi z populacji o rozkładzie normalnym.

H1: Próba nie pochodzi z populacji o rozkładzie normalnym..

selfesteem  %>%
    group_by(time) %>%
    shapiro_test(score)

## # A tibble: 3 x 4
##   time  variable statistic     p
##   <fct> <chr>        <dbl> <dbl>
## 1 t1    score        0.967 0.859
## 2 t2    score        0.876 0.117
## 3 t3    score        0.923 0.380

Na podstawie powyższych wyników możemy zauważyć, że w każdej podgrupie poziom p >alfa zatem nie ma podstaw do odrzucenia H0. W związku z powyższym wszystkie analizowane powyżej próbki pochodzą z populacji o rozkładzie normalnym.

ANOVA

Ze względu na to, że założenie o normalności jest spełnione możemy wykonać jednoczynnikową parametryczną ANOVĘ dla prób zależnych. W niniejszym teście:

H0: Z czasem samopoczucie badanych nie uległo istotnie statystycznej poprawie. H1: Z czasem samopoczucie badanych uległo istotnie statystycznej poprawie.

wyniki <- anova_test(data=selfesteem, dv=score, wid=id, within=time)

get_anova_table(wyniki)

## ANOVA Table (type III tests)
## 
##   Effect DFn DFd      F        p p<.05   ges
## 1   time   2  18 55.469 2.01e-08     * 0.829

Na podstawie powyższych wyników można zauważyć, że poziom p < alfa zatem należy odrzucić H0. W związku z powyższym stwierdzamy, że z czasem samopoczucie badanych uległo istotnie statystycznej poprawie.

Powyższy test można również wizualizować:

ggwithinstats(data=selfesteem, 
              y=score,
              x=time,
              type="p",
              xlab="Punkty czasowe diety",
              ylab="Ocena samopoczucia",
              legend.title="Momenty czasowe czasowe",
               caption="ANOVA jednoczynnikowa dla prób zależnych")

W oparciu o powyższy rysunek można ponownie odrzucić H0 o braku istotnie statystycznej poprawy samopoczucia następującej wraz z upływem czasu.

Test post-hoc

Ze względu na to, iż odrzucono H0 w jednoczynnikowej ANOVIE należy wykonać test post-hoc, aby zidentyfikować między którymi punktami czasowymi podczas stosowania diety wystąpiły istotne statystycznie różnice w otrzymanych ocenach samopoczucia. W tym przypadku będzie to test prawdziwie istotnej różnicy HSD Tukeya. W prezentowanym teście:

H0: Między średnimi ocen samopoczucia w poszczególnych punktach czasowych diety nie występują statystycznie istotne różnice.

H1: Między średnimi ocen samopoczucia w poszczególnych punktach czasowych diety występują statystycznie istotne różnice.

model <- aov(score~time,data=selfesteem)

hsd<-TukeyHSD(model,which="time", ordered=TRUE)
print(hsd)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##     factor levels have been ordered
## 
## Fit: aov(formula = score ~ time, data = selfesteem)
## 
## $time
##          diff       lwr      upr     p adj
## t2-t1 1.79382 0.8114139 2.776226 0.0003104
## t3-t1 4.49622 3.5138136 5.478626 0.0000000
## t3-t2 2.70240 1.7199935 3.684806 0.0000007

Na podstawie powyższych wyników można zauważyć, że poziom p < alfa dla każdej z badanych par punktów czasowych diety. W związku z powyższym należy odrzucić H0. Zatem możemy uznać, że występują statystycznie istotne różnice w ocenach samopoczucia, które dokonano w poszczególnych punktach czasowych diety.

Wnioski i podsumowanie

W efekcie przeprowadzonego badania sformułowano następujące konkluzje:

1. Z czasem stosowania diety samopoczucie badanych uległo istotnie statystycznej poprawie.

2. Wzrost poprawy samopoczucia badanych był istotny statystycznie między wszystkimi punktami czasowymi stosowania diety.

Podsumowując, statystyka matematyczna to potężne i zarazem fascynujące narzędzie, które pozwala stwierdzić czy obserwowane w otaczającym świecie różnice rzeczywiście mają miejsce czy są one jedynie kwestią przypadku. Ponadto do sformułowania niniejszych konkluzji niepotrzebne jest zbadanie całej populacji - co w większości przypadków byłoby niemożliwe - lecz wystarczą jedynie próbki!

Autorka prezentowanego projektu wyraża wdzięczność za możliwość zrealizowania tak fascynującego przedmiotu z zastosowaniem profesjonalnego oprogramowania analitycznego jakim jest język R. Składa serdeczne podziękowania dr. inż. Karolowi Flisikowskiemu za przekazaną wiedzę i umiejętności, a także wsparcie merytoryczne na konsultacjach oraz życzy dalszych sukcesów w rozwijaniu przedmiotów analitycznych i nie tylko!