---
title: "Analiza wariancji - ANOVA"
author: "Bartosz Kozak"
date: "6 maja 2019"
output:
  html_document:
    df_print: paged
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Wprowadzenie

Na dzisiejszych zajęciach zapoznamy się z jednoczynnikową analizą wariancji - ANOVA. Nasze rozważania rozpoczniemy od szybkiego zapoznania z teorią stanowiącą podstawę tej analizy.

Niech $x_{ij}$ oznacza *j-otą* obserwację w *i-tej* grupie, wtedy $x_{35}$, oznacza piątą obserwację w trzeciej grupie; $\overline{x}_i$ oznacza średnią dla *i-tej* grupy, a $\overline{x}_.$ oznacza średnią całkowitą (średnią z wszystkich obserwacji).

Możemy rozbić obserwację na czynniki jako:
\begin{equation}
\label{eq:1}
x_{ji}=\overline{x}_. + (\overline{x}_{i}-\overline{x}_.) + (x_{ij}-\overline{x}_i)
\end{equation}
gdzie:

$\overline{x}_{i}-\overline{x}_.$ - odpowiada różnicy średniej grupowej od średniej całkowitej (wpływ grupy)

$x_{ij}-\overline{x}_i$ - odpowiada różnicy obserwacji od średniej grupowej (wpływ błędu obserwacji)

Koresponduje to z modelem
\begin{equation}
\label{eq:2}
X_{ij} = \mu + \alpha_i + \epsilon_{ij}, \:\:\:\: \epsilon_{ij} ~ N(\mu,\sigma)
\end{equation}
gdzie:

$\alpha_i$ - oznacza wpływ czynnika (grupy)

$\epsilon_{ij}$ - oznacza błąd pomiaru 

W modelu tym zakładamy hipotezę zerową, że wpływ czynnika ($\alpha_i$) jest równy zero dla wszystkich badanych czynników. Zauważmy, że zakładamy niezależność czynnika błędu, oraz że ma on tą samą wariancję dla wszystkich czynników.

Rozważmy teraz sumy kwadratów powyższych czynników znanych jako *wariancję wewnątrzgrupową* - $SSD_{W}$
\begin{equation}
\label{eq:3}
SSD_{W} = \sum_i\sum_j(x_{ij}-\overline{x}_i)^2
\end{equation}
oraz *wariancję międzygrupową* $SSD_{B}$
\begin{equation}
\label{eq:4}
SSD_{B} = \sum_i\sum_j(x_{i}-\overline{x}_.)^2=\sum_in_i(x_i-\overline{x}_.)^2
\end{equation}
Możemy też wykazać, że
\begin{equation}
\label{eq:5}
SSD_{W}+SSD_{B}=SSD_{całkowita}= \sum_i\sum_j(x_{ij}-\overline{x}_.)^2
\end{equation}

Należy jednak zaznaczyć, że suma kwadratów (SS) może być tylko pozytywna, więc nawet kompletnie nieinformatywne (przypadkowe) grupowanie wyjaśni jakąś część obserwowanej zmienności.  Nasuwa się więc pytanie jak małe wyjaśnienie całkowitej wariancji możemy uznać za przypadkowe? Okazuje się, że w przypadku braku systematycznych różnic między grupami należy oczekiwać, że suma kwadratów podzielona przez liczbę stopni swobody dla każdego czynniki, $k-1$ dla $SSD_{B}$ oraz $N-k$ dla $SSD_{W}$, gdzie $k$ oznacza liczbę grup, a $N$ oznacza całkowitą liczbę obserwacji (sumę liczby obserwacji z wszystkich grup). 

Możemy więc obliczyć *średnie kwadraty* jako:
\begin{equation}
\label{eq:6}
MS_W=SSD_W/(N-k)
\end{equation}

\begin{equation}
\label{eq:7}
MS_B=SSD_B/(k-1)
\end{equation}

$MS_W$ jest połączoną wartością wariancji poszczególnych grup i jest estymatorem $\sigma^2$. W przypadku braku wpływu czynnika grupy $MS_B$ także jest estymatorem $\sigma^2$, jednakże jeżeli wpływ czynnika grupy będzie istotny wtedy różnice między grupami, a co za tym idzie $MS_B$ będzie "duże". Dlatego test na różnice między średnimi grup może zostać wykonany jako porównanie dwóch estymatorów wariancji. Dlatego opisywana procedura pomimo tego że służy do porównania średnich między grupami nosi nazwę *analizy wariancji*. Z poprzednich zajęć wiemy, że do porównania wariancji wykorzystujemy $test F$ oraz dystrybucję $F$. Wartość statystyki testowej obliczamy ze wzoru 
\begin{equation}
\label{eq:8}
F=MS_B/MS_W
\end{equation}
i porównujemy z wartością graniczną $F_{\alpha}$ z rozkładu F o liczbie stopni swobody: $(k-1)$ oraz $(N-k)$.

Testujemy hipotezę zerową postaci: $$H_0:\frac{MS_B}{MS_W}=1$$ wobec hipotezy alternatywnej: $$H_A:\frac{MS_B}{MS_W}\neq1$$

## Analiza wariancji - algorytm postępowania

**Przykład 1**

Przeprowadzono eksperyment w celu porównania trzech pożywek (A, B, C). Obliczono liczbę koloni bakterii (w trzech powtórzeniach) dla każdej pożywki. Otrzymane dane zapisano w pliku Cw5_1.txt. Czy pożywki wpływają na liczbę kolonii bakterii?

Do odpowiedzi na postawione pytanie wykorzystamy jednoczynnikową analizę wariancji "One-way ANOVA".

**Algorytm postępowania**

1. Formujemy $H_0$ i $H_A$.

   $H_0:$ pożywka **nie wpływa** na liczbę kolonii. Jest to równoważne stwierdzeniu, że średnie liczby kolonii dla poszczególnych pożywek **są równe**.

   $H_A$: Pożywka **wpływa** na liczbę kolonii bakterii. Jest to równoważne stwierdzeniu, że średnie liczby kolonii dla poszczególnych pożywek **nie są równe**.

2. Obliczamy statystykę $F$

3. Znajdujemy  wartość krytyczną $F\alpha$ (dla zadanej $\alpha$ i odpowiednich df)

4. Porównujemy $F$ z $F\alpha$

5. Wnioskujemy

```{r, echo=FALSE}
dane <- read.csv('~/konas13@gmail.com/Dokumenty/Dydaktyka/Statystyka_Ochrona_Srodowiska/Cw_5_ANOVA/Cw5_1.txt', sep = "\t")
```

```{r}
# wczytanie danych do R
dane
```

**ad 2.** Obliczamy Statystykę $F$.

Korzystamy z wzoru \ref{eq:3} i obliczamy $SSD_W$.
```{r, echo=FALSE}
sr_a <- mean(dane$M_A)
sr_b <- mean(dane$M_B)
sr_c <- mean(dane$M_C)

ss_a <- sum((dane$M_A-sr_a)^2)
ss_b <- sum((dane$M_B-sr_b)^2)
ss_c <- sum((dane$M_C-sr_c)^2)

ssd_w <- sum(c(ss_a,ss_b,ss_c))
```
```{r}
# obliczamy średnie dla grup A,B oraz C
c(sr_a,sr_b,sr_c)
# obliczamy sumy kwadratów dla grup A, B oraz C
c(ss_a,ss_b,ss_c)
# obliczamy ssd_w
ssd_w
```
Korzystamy z wzoru \ref{eq:5} i obliczamy $SSD_{całkowity}$
```{r, echo=FALSE}
GM <- sum(dane)/(dim(dane)[1]*dim(dane)[2])
ssd_c <- sum((dane-GM)^2)
```

```{r}
# obliczamy średnią całkowitą
GM
# obliczamy ssd_c
ssd_c
```
Następnie ponownie korzystamy z wzoru \ref{eq:5} i obliczamy $SSD_B$
```{r, echo=FALSE}
ssd_b <- ssd_c - ssd_w
```

```{r}
# obliczamy ssd_b
ssd_b
```

W oparciu o wzory \ref{eq:6} i \ref{eq:7} obliczamy $MS_B$ oraz $MS_W$
```{r, echo=FALSE}
N <- dim(dane)[1]*dim(dane)[2]
k <- ncol(dane)
ms_b <- ssd_b/(k-1)
ms_w <- ssd_w/(N-k)
```

```{r}
# liczba obserwacji N
N
# liczba grup k
k
# obliczamy ssd_w
ms_b
# obliczamy ms_w
ms_w
```

W oparciu o wzór \ref{eq:8} obliczamy statystykę testową $F$
```{r, echo=FALSE}
F <- ms_b/ms_w
```

```{r}
# obliczamy F
F
```

**ad3** Znajdujemy wartość krytyczną $F_{\alpha}$, dla $\alpha$ = 0.05 i liczby stopni swobody $k-1$ oraz $N-k$
```{r}
# znajdujemy Fa
Fa <- qf(0.95,(k-1),(N-k))
Fa
```

Po porównaniu $F$ z $F_{\alpha}$. Wyniki analizy anova zwyczajowo prezentuje się w tabeli analizy wariancji.
```{r}
# wektor z wartościami sum kwadratów
ss <- c(ssd_b,ssd_w,ssd_c)
# wektor z wartościami średnich sum kwadratów
ms <- c(ms_b,ms_w,"")
# wektor z liczbami stopni swobody
df <- c((k-1),(N-k),(N-1))
# wektor z wartościami statystyki F
f <- c(F,"","")
# wektor z wartościami F krytyczne
fa <- c(round(Fa,2), "","")
# tworzymy tabelę
anova_tabela <- as.data.frame(matrix(c(df,ss,ms,f,fa),nrow = 3))
# zmieniamy nazwy kolumn i wierszy
colnames(anova_tabela) <- c("df", "SS", "MS", "F", "Fa")
row.names(anova_tabela) <- c("pożywka", "błąd", "całkowita")
# wyświetlamy tabele analizy wariancji
anova_tabela
```
Jaki wniosek możemy wyciągnąć względem testowanej hipotezy ($H_0$)?

## Anova w R

ANOVA może być wykonana w R wyłącznie na danych typu `stack`. Musimy więc przekształcić dane.
```{r}
dane1 <- stack(dane)
colnames(dane1) <- c("liczba.kolonii","pozywka")
dane1
```

Do wykonania analizy ANOVA w R użyjemy funkcji `lm` oraz `anova`

```{r}
# budujemy model
model <- lm(liczba.kolonii~pozywka,data=dane1)
# wykonujemy ANOVA
anova(model)
```
Wartość `Pr` w tebeli wariancji odpowiada wartości pvalue znanej z testu $t$.
Ponieważ wartość `Pr` jest mniejsza od $\alpha$ odrzucamy $H_0$ na rzecz $H_A$ z prawdopodobieństwem 0.95%.

## Graficzna prezentacja wyników - wykres pudełkowy

Często wykorzystywaną metodą prezentacji wyników w pakietach statystycznych są wykresy pudełkowe "box-plot". W programie R do stowrzenia tego typu wykresów używamy funkcji `boxplot`.

```{r}
boxplot(liczba.kolonii~pozywka,data=dane1)
```

**Zadanie do samodzielnego rozwiązania**

W celu zbadania plonowania pomidora na różnego rodzaju podłożach zostało w szklarni założone doświadczenie metodą kompetnej randomizacji w 6 powtórzeniach (na powtórzenie składała się skrzynka lub mata z trzema roślinami). Badano następujące rodzaje podłoża: 

- A słoma żytnia
- B słoma żytnia + torf (3:1)
- C słoma pszenna
- D słoma pszenna + kora (3:1)
- E wełna mineralna

Badano plon z rośliny (kg/roślina) w 22 tygodniu, dane zamieszczono w pliku Cw5_2.txt. Wykonaj analizę wariancji wg właściwego modelu, aby stwierdzić czy zastosowane podłoże ma pływ na uzyskany plon pomidorów ($\alpha$=0.05)? Przedstaw uzyskane średnie plony dla badanych podłoży graficznie (wykres pudełkowy, lub wykres słupkowy) .
