---
title: "Analiza wariancji"
output:
  pdf_document: default
  html_document: default
---

Na początek kilka ważnych faktów o analizie wariancji.

* Jest powszechnie stosowana metoda statystyczna pozwalająca na ocenę istotności różnic wielu średnich.
* Opracował ją i upowszechnił Ronald Fisher w latach 20 XX wieku.
* Jest to technika badania wyników (doświadczeń, obserwacji), które zależą od jednego (jednoczynnikowa analiza wariancji) lub kilku czynników działających równocześnie (dwuczynnikowa lub wieloczynnikowa analiza wariancji).
    * przykłady takich czynników: metody leczenia, płeć, sposoby nawożenia (Fisher!).
    * czynniki takie nazywamy *czynnikami grupującymi, klasyfikacyjnymi lub zmiennymi manipulacyjnymi*.
    * Każdy czynnik ma kilka poziomów lub wariantów (częstym wzorcem jest np. mało-średnio-dużo czegoś)
* Założenia analizy wariancji są podobne, jak w przypadku regresji liniowej - homogeniczność wariancji w grupach i normalność w grupach.

## Model teoretyczny

Model strukturalny leżący u podstaw analizy wariancji ma następującą postać:

$$X_{ij} = \mu + \tau_j + \epsilon_{ij}$$

* $X_{ij}$ - indywidualna obserwacja X numer $i$ należąca do grupy $j$
* $\mu$ - ogólna średnia
* $\tau_j$ - efekt przynależenia do grupy $j$
* $\epsilon_{ij}$ - błąd

## Układ hipotez przy analizie wariancji

Przy analizie wariancji mamy do czynienia z następującym układem hipotez:

$$
H_0: \mu_1 = \mu_2 = ... = \mu_k
$$

W tym wypadku $\mu_1$ do $\mu_k$ odnoszą się do średnich w populacji dla poszczególnych poziomów czynnika grupującego. Hipoteza zerowa stwierdza, że średnie te są sobie równe.

$$
H_A: \neg H_0
$$

Hipoteza alternatywna (ta, którą chcemy uzasadnić obalając hipotezę zerową!) głosi, że $H_0$ jest fałszem. Warto zwrócić uwagę, że $H_0$ jest sformułowana w bardzo specyficzny sposób i może okazać się fałszywa na wiele sposobów. Wszystkie średnie mogą być od siebie różne albo tylko niektóre. Sam test jednak powie nam jedynie tyle, że powinniśmy odrzucić hipotezę zerową. Czy możemy mimo wszystko dowiedzieć się czegoś na temat tego, w jaki sposób $H_0$ okazała się fałszywa? Tak! Żeby to zrobić, musimy przeprowadzić analizę *post-hoc*. O tym jednak powiemy sobie więcej na koniec tej części kursu (Uwaga! Możemy to zrobić wyłącznie wówczas, gdy ANOVA dała istotny statystycznie wynik).

## Wariancja międzygrupowa i wewnątrzgrupowa

Logika analizy wariancji polega na porównywaniu wariancji *międzygrupowej* i *wewnątrzgrupowej*. Korzystamy z faktu, że -- przy założeniu hipotezy zerowej -- na podstawie obu możemy skonstruować estymator wariancji populacji. Jeśli hipoteza zerowa jest zaś fałszywa, to estymator ,,działa"' tylko w jednym przypadku. Mechanizm ten będzie jaśniejszy, jeśli prześledzimy analizę wariancji na przykładzie.

Stwórzmy ramkę danych, w której mamy poziomów czynnika (jakiegokolwiek, tutaj użyhemy niewiele znaczacych poziomów ,,bardzo mało'', ,,mało'', ,,średnio'', ,,dużo'' oraz ,,bardzo dużo''), ale wszystkie wartości pochodzą z populacji o tych samych parametrach (takiej samej średniej oraz odchyleniu standardowym). Łatwo zauważyć, że jest to sytuacja, w której hipoteza zerowa jest prawdziwa.

```{r}
# Nasz czynnik grupujący ma 5 poziomów, po 200 obserwacji każdy
faktor <- rep(c('bardzo mało', 'mało', 'średnio', 'dużo', 'bardzo dużo'), 200)
# Losujemy 1000 obserwacji z rozkładu o mu = 10 i sd = 5; po 200 na każdy poziom czynnika
wartosci <- rnorm(1000, mean = 10, sd = 5)
# Tworzymy z naszych wylosowanych obserwacji ramkę danych
df <- data.frame(condition = faktor, value = wartosci)
# Wyświetlamy pierwsze 6 wierszy
head(df)
```

### Wariancja wewnętrz grup

Obliczmy wariancje w każdej z grup za pomocą funkcji `tapply`. 

Krótkie przypomnienie jak działa `tapply`. Funkcja ta bierze wartości przekazane jako pierwszy argument (tutaj: `df$value`) i dzieli je według czynnika przekazanego jako drugi argument (tutaj: `df$condition`). Następnie do każdego elementu takiego podziału stosuje funkcję przekazaną jako trzeci argument -- tutaj jest to `var` obliczające wariancję. Funkcja ta zwróci więc tyle wartości wariancji z próby (`var`) ile mieliśmy poziomów czynnika grupującego (`df$condition`) -- w naszym wypadku 5.

```{r}
tapply(df$value, df$condition, var)
```

W następnym kroku obliczymy średnią z obliczonych właśnie wariancji. W ten sposób skonstruowaliśmy pierwszy z dwóch estymator wariancji populacji -- obliczyliśmy wariancję w każdej z grup z osobna i wyciągnęliśmy z nich średnią.

```{r}
mse <- mean(tapply(df$value, df$condition, var))
mse
```

Częściej spotkają się Państwo ze sposobem obliczania, w którym najpierw obliczamy sumę kwadratów a następnie średni kwadrat według następujących wzorów:


$$
SSE  = \sum_{i=1}^{k} \sum_{j=1}^{n_i}(X_{ij} - \bar{X_i})^2
$$

$$
MSE = \frac{SSE}{n-k}
$$
Sposób ten jest równoważny temu, co zrobiliśmy wcześniej, co można prześledzić wykonując umieszczony poniżej kod. W moim przekonaniu w tym sposobie nieco gorzej widać,,o co chodzi''.

```{r}
# Dzielimy wyjściową ramkę danych według czynnika z kolumny `condition` 
df_s <- split(df, df$condition)
# Obliczamy sumę kwadratów odchyleń
sse <- sum(sapply(df_s, function(x){sum((x$value - mean(x$value))^2)}))
sse

# Obliczamy średni kwadrat
mse <- sse/(1000-5)
mse
```

### Wariancja między grupami

Nasz drugi estymator wariancji w populacji skonstruujemy, korzystajac z nabytej wcześniej wiedzy dotyczącej odchylenia standardowego średnich z próby. Innymi słowy - kolejny raz skorzystamy z Centralnego Twierdzenia Granicznego!

Wiemy z CTG, że:

$$
\frac{\sigma^2}{n} = s^2_{\bar{X}}
$$

Ten wzór możemy przekształcić w taki sposób, że uzyskamy:

$$
\sigma^2 = ns^2_{\bar{X}}
$$
Widzimy więc, że wariancja w populacji ($\sigma^2$) to po prostu wariancja średnich z próby ($s^2_{\bar{X}}$) pomnożona przez liczebność próby.

W takim razie, żeby wyestymować wariancję w populacji, powinniśmy obliczyć wariancję średnich z próby:

```{r}
# Za pomocą `tapply` obliczamy średnią dla każdego poziomu `df$condition`
# Takich średnich w naszym wypadku otrzymujemy 5
# `var` oblicza wariancje dla tych 5 średnich
var(tapply(df$value, df$condition, mean))
# Żeby otrzymać estymator wariancji w populacji, musimy przemnożyć nasz wynik przez n = 200
# Bo tyle obserwacji mamy na każdy poziom czynnika grupującego
mstreat <- var(tapply(df$value, df$condition, mean))*200
mstreat
```

Częściej spotkają się Państwo ze sposobem obliczania, w którym najpierw obliczamy sumę kwadratów a następnie średni kwadrat według następujących wzorów:

$$
SSA = \sum_{i=1}^{k} (\bar{X_i} - \bar{X})^2 n_i
$$

$$
MSA = \frac{SSA}{k-1}
$$

```{r}
ssa <- (sum((tapply(df$value, df$condition, mean) - mean(df$value))^2))*200
ssa

msa <- ssa/(5-1)
msa
```

Zwróćmy uwagę, że obliczona przez nas wartość będzie dobrym estymatorem wariancji w populacji wyłącznie wtedy, gdy próby pochodzą z tej samej populacji, ponieważ tylko takiej sytuacji dotyczy przywołane sformułowanie Centralnego Twierdzenia Granicznego! W innym wypadku tak obliczony ,,estymator'' nie będzie działać i ten fakt wykorzystamy do skonstruowania testu statystycznego! 

## Statystyka testowa

Statystyka testowa $F$ będzie proporcją między dwiema obliczonymi własnie wartościami. W przypadku, kiedy grupa nie ma znaczenia (czyli, biorac pod uwagę nasz model teoretyczny, $\theta = 0$), spodziewamy się, że statystyka ta będzie wynosić 1 (taka jest jej wartość oczekiwana). Dlaczego? Ponieważ obie wartości są estymatorami tego samego parametru!

$$
F = \frac{MSA}{MSE}
$$

```{r}
f_stat <- mstreat/mse
f_stat
```

Czy wielkośc statystyki testowej, którą otrzymaliśmy, jest duża czy mała? Możemy zobaczyć, jaki rozkład ma interesująca nas statystyka testowa, przeprowadzając odpowiednią symulację. 

Możemy również posłużyc się rozkładem teoretycznym $F$. Ma on dwa parametry - liczbę stopni swobody licznika (czyli $k-1$, gdzie $k$ to liczba czynników) oraz liczbę stopni swobody licznika ($n-k$, gdzie $k$ to liczba czynników a $n$ to całkowita liczebność próby).

Na początek symulacja. Sprawdzamy, jaki rozkład ma statystyka testowa $F$ przy założeniu prawdziwości hipotezy zerowej. W tym celu powtórzymy nasz eksperyment 10 000 razy i obliczamy dla każdego powtórzenia statystykę $F$. Następnie nakładamy na histogram statystyk $F$ uzyskanych w symulacji odpowiednią krzywą z teoretycznego rozkładu $F$. Finalnie możemy zobaczyć, że histogram z symulowanymi wartościami $F$ oraz krzywa rozkładu teoretycznego świetnie do siebie pasują.

```{r}
f_val <- replicate(10000,
                   {
                     # Najpierw symulujemy naszą próbę
                     faktor = rep(c('bardzo mało', 'mało', 'średnio', 'dużo', 'bardzo dużo'),
                                  200)
                     wartosci <- rnorm(1000, mean = 10, sd = 5)
                     df <- data.frame(condition = faktor, value = wartosci)
                     # Obliczamy MSE i MS_Treatment 
                     mse <- mean(tapply(df$value, df$condition, var))
                     mstreat <- var(tapply(df$value, df$condition, mean)) * 200
                     # Obliczamy statystykę testową F czyli iloraz tych dwóch wartości
                     mstreat / mse
                   })
# Wyniki symulacji możemy przedstawić na histogramie
hist(f_val, freq = FALSE)
curve(df(x, df1 = 5-1, df2 = 995), add = TRUE)
```

Skoro obliczyliśmy już statystykę testową $F$, to musimy jeszcze sprawdzić, czy uzyskanie takiej lub bardziej skrajnej wartości jest mniej prawdopodobne niż $5\%$ (przyjmijmy taki konwencjonalny poziom $\alpha$). Żeby to zrobić, musimy posłużyć się dystrybuantą rozkładu $F$, czyli funkcją `pf`. Rozkład $F$ ma dwa parametry (dwie liczby stopni swobody). `df1` będzie tutaj wynosiła $k-1$ gdzie $k$ jest liczbą poziomów czynnika grupującego, a `df2` będzie wynosiło $n - k$ gdzie $n$ to (całkowita) liczebność próby.

```{r}
pf(f_stat, df1 = 5-1, df2 = 995, lower.tail = FALSE)
```

Obliczyliśmy właśnie p-wartość i tym samym przeprowadziliśmy naszą pierwszą analizę wariancji. Oczywiście przeprowadzając analizę wariancji nie musimy wykonywać ręcznie obliczeń. Jak zwykle możemy się posłużyć wbudowaną funkcją R, którą w tym wypadku jest `aov`.

```{r}
summary(aov(value ~ condition, df))
```

## Analiza wariancji na przykładzie

Prześledźmy analizę wariancji na przykładzie. W bibliotece standardowej R znajduje się zbiór danych `PlantGrowth`, dotyczący wpływu pewnych czynników na obfitość plonów. W jednej kolumnie mamy wagę roślin, w drugiej to, do jakiej grupy przynależy dana obserwercja. Grupy były trzy - dwie eksperymentalne (`trt1` oraz `trt2`) i oraz grupa kontrolna (`ctrl`). Za pomocą analizy wariancji możemy zbadać, czy warunki eksperymentalne miały wpływ na wagę roślin.

Nasza hipoteza zerowa będzie stwierdzać, że średnia waga jest taka sama we wszystkich grupach. Innymi słowy, że obserwacje dla wszystkich trzech poziomów czynnika grupującego pochodzą z populacji o tym samym parametrze średniej. Hipoteza alterantywna będzie zaś temu zaprzeczać. Oto zbiór danych (proszę wybaczyc, że wydrukuje go w całości, ale jest dość krótki):

```{r}
PlantGrowth
```

Na początek zawsze warto stworzyć wizualizację danych - między innymi po to, aby skontrolować to, czy nasze dane spełniają założenia analizy wariancji. Zasadniczo założenia są, podobnie jak przy regresji, dwa - równa wariancja oraz normalność w grupach. Na podstawie wykresu pudełkowego możemy "na oko" ocenić, czy mamy prawo zakładać, że założenia są spełnione. Te dane nie wyglądają bardzo źle -- nie mamy powodów by powiedzieć, że w sposób oczywisty dane nie spełniają założeń.

```{r}
library(ggplot2)
qplot(x = group, y = weight, geom = 'boxplot', data = PlantGrowth)
```

Normalność w grupach możemy "na oko" ocenić za pomocą wykresu kwantyl-kwantyl dla każdej grupy z osobna. Przy małych próbach może to okazać się to bardziej pomocne niż formalne testy na normalność rozkładu. Nie chcemy, żeby w sposób wyraźny i podpadający pod jakiś wzór odbiegały od przekątnej.

```{r fig.width=5, fig.height=12}
p <- ggplot(PlantGrowth, aes(sample = weight))
# scales = "free" sprawia, że skale na każdym wykresie są osobne.
p + 
  stat_qq() + 
  facet_grid(group~., scales = 'free') 
```

Następnie możemy "ręcznie" przeprowadzić analizę wariancji w R. Najpierw musimy policzyć sumę kwadratów między grupami i wewnatrz grup, następnie obliczyć średni kwadrat w obu przypadkach. Po obliczeniu tych dwóch rzeczy możemy obliczyć statystykę testową F, dzieląc jedną przez drugą, i sprawdzić, korzystając z dystrybuanty rozkładu F, czy odrzucamy hipotezę zerową.

```{r collapse= T }
df_s <- split(PlantGrowth, PlantGrowth$group)

sse <- sum(sapply(df_s, function(x){sum((x$weight - mean(x$weight))^2)}))
print('SSE')
sse

mse <- sse/(nrow(PlantGrowth) - (length(levels(PlantGrowth$group))))
print('MSE')
mse


ssa <-(sum((tapply(PlantGrowth$weight, PlantGrowth$group, mean) - 
           mean(PlantGrowth$weight))^2))*(nrow(PlantGrowth)/3)
print('SSA')
ssa

msa <- ssa/(length(levels(PlantGrowth$group))-1)
print('MSA')
msa

f_stat <- msa/mse

print('Statystyka testowa F')
f_stat

pval <- pf(f_stat, (length(levels(PlantGrowth$group))-1), (nrow(PlantGrowth) - (length(levels(PlantGrowth$group)))), 
          lower.tail = FALSE)
print('Wartość p')
pval
```

To samo możemy wykonać za pomoca wbudowanej funkcji `aov`. Jak widzimy, uzyskane w eksperymencie rezultaty uprawniają nas do odrzucenia hipotezy zerowej. Możemy więc stiwerdzić, że proby nie pochodzą z populacji o takich samych średnich.

```{r}
fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)
```

## Testy *post-hoc*

Często nie tylko interesuje nas fakt, że odrzuciliśmy hipotezę, zgodnie z którą wszystkie próby pochodzą z populacji o  tej samej średniej, ale chcielibyśmy móc ustalić, między którymi grupami występują statytycznie istotne różnice. Tutaj przydatne okazują się tzw. testy *post-hoc*. Możemy ich użyć tylko jeśli nasza główna metoda (tutaj -- ANOVA) dała statystycznie istotny wynik. Zaletą testów *post-hoc* jest to, że uwzględniają fakt, że są *post-hoc* i podają bardziej wiarygodne wyniki niż np. seria testów *t*. Testów *post-hoc*, które możemy zastosować przeprowadziwszy analizę wariancji jest kilka, najbardziej popularnym jest chyba test HSD (*honest significant difference*) Tukeya. Możemy taki test łatwo przeprowadzić za pomocą wbudowanej w R funkcji `TukeyHSD`. Na poniższym wydruku widzimy, że statystycznie istotna jest różnica pomiędzy dwoma grupami eksperymentalnymi (ale już nie między każdą z nich a grupą kontrolną!). 

```{r}
TukeyHSD(fit)
```

Możemy wyniki naszej analizy przedstawić w formie graficznej za pomocą funkcji `plot`. Na osi Y mamy zestawienia ze sobą poziomów czynnika (tutaj - `group`), na osi X zaś 95% przedziały ufności różnicy średnich między nimi. Istotne statystycznie są te różnice, które nie obejmują zera (proszę przypomnieć sobie testowanie hipotez za pomocą przedziałów ufności).

```{r}
plot(TukeyHSD(fit))
```

## Wielkość efektu

W przypadku analizy wariancji standardowo używa się miary $\eta^2$. Oblicza się ją dzięląc międzygrupową sumę kwadratów przez całkowitą sumę kwadratów.

```{r}
# Całkowita suma kwadratów
sst <- sum((PlantGrowth$weight - mean(PlantGrowth$weight)) ** 2)
eta_kw = ssa/sst
eta_kw
```

Interpretacja tej miary jest analogiczna jak w przypadku znanego już Państu $R^2$ i jest to stosunek zmienności wyjaśnianej przez nasz predyktor do całkowitej zmienności zmiennej zależnej. Jeżeli mamy taką ochotę, to mozemy skorzystać z funkcji obliczającą tę miarę w jednej z wielu bibliotek, które taką funkcję oferują. Ja skorzystałem tutaj z `DescTools` i funkcji `EtaSq` (ale jest ich wiele!).

```{r}
DescTools::EtaSq(fit)
```

## Analiza wariancji a model liniowy

Istnieją bardzo podobne podobieństwa między analizą wariancji i modelem liniowym. Można nawet stwierdzić, że analiza wariancji to po prostu szczególny przypadek modelu liniowego. Fakt, że omawiane są oddzielnie podyktowane jest w znacznej mierze zaszłościami historycznymi. W artykule dostępnym w sekcji dla zalogowanych "Multiple Regression As A General Data-Analytic System" autorstwa Jacoba Cohena znajduje się dobre wyjaśnienie wielu zawiłości z tym związanych. 

Zobaczmy sobie podsumowanie jednoczynnikowej analizy wariancji oraz modelu liniowego z jednym kategorialnym predyktorem.

**Analiza wariancji:**

```{r}
fit_aov <- aov(weight ~ group, data = PlantGrowth)
summary(fit_aov)
```

```{r}
DescTools::EtaSq(fit_aov)
```

**Model liniowy: **

```{r}
fit_lm <- lm(weight ~ group, data = PlantGrowth)
summary(fit_lm)
```

Zwróćmy uwagę na kilka rzeczy. Dla analizy wariancji hipoteza zerowa wygląda tak:

$$
H_0: \mu_1 = \mu_2 = \mu_3
$$

Dla modelu liniowego hipotezą zerową było:

$$
H_0: R^2 = 0
$$

Hipotezy te, chociaż pozornie brzmiące różnie, w rzeczywistości są sobie matematycznie równoważne i można je przetestować tym samym testem (na co wskazuje nam w wydruku R taka sama wartość statystyki testowej $F$). Możemy to zilustrować takim oto rozumowaniem. 

Jeżeli hipoteza zerowa analizy wariancji jest prawdziwa (czyli średnie w grupach są w populacji identyczne), to znajomość tego, do której grupy dana obserwacja należy, nie pozwala nam przewidywać wartości tej obserwacji lepiej niż po prostu znajomość ogólnej średniej. W takiej sytuacji nasz model nie wyjaśnia żadnej części wariancji zmiennej zależnej czyli $R^2$ jest równe zeru. Warto również zwrócić uwagę na to, że $R^2$ jest równy obliczonemu przez nas $\eta^2$, co nie powinno dziwić biorać pod uwagę sposób w jaki oblicza się obie te wartości.

Na kolejnych zajęciach, gdy będziemy zajmować się nieco bardziej skomplikowanymi wariantami analizy wariancji (w których uwzględniamy więcej czynników), zobaczymy w jaki sposób tę paralelę między modelem liniowym a analizą wariancji można pociągnąć dalej.