---
title: "Post hoc test"
author: "Bartosz Kozak"
date: "13 maja 2019"
output:
  html_document:
    df_print: paged
header-includes:
- \usepackage{multirow}
- \usepackage[table,xcdraw]{xcolor}
---

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

## Wprowadzenie

W czasie testowania hipotezy możemy popełnić błąd polegający na odrzuceniu prawdziwej hipotezy zerowej (**błąd pierwszego typu**), lub błąd polegający na przyjęciu hipotezy zerowej kiedy jest ona nieprawdziwa (**błąd drugiego rodzaju**). Sytuację tą ilustruje tabela \ref{tab:1}. W tabeli tej wiersze korespondują do nieznanej badaczowi prawdy, natomiast kolumny reprezentują wyniki przeprowadzonych testów statystycznych. 

\begin{table}[h]
\caption{Błąd I i II rodzaju w testowaniu hipotez}
\label{tab:1}
\begin{tabular}{|
>{\columncolor[HTML]{FFFFFF}}l 
>{\columncolor[HTML]{FFFFFF}}l |
>{\columncolor[HTML]{FFFFFF}}l 
>{\columncolor[HTML]{FFFFFF}}l |}
\hline
\multicolumn{2}{|l|}{\cellcolor[HTML]{FFFFFF}{\color[HTML]{343434} }}                                 & \multicolumn{2}{c|}{\cellcolor[HTML]{FFFFFF}{\color[HTML]{343434} Test}}                                                                                     \\ \cline{3-4} 
\multicolumn{2}{|l|}{\multirow{-2}{*}{\cellcolor[HTML]{FFFFFF}{\color[HTML]{343434} }}}               & \multicolumn{1}{c}{\cellcolor[HTML]{FFFFFF}{\color[HTML]{343434} Nie istotny}} & \multicolumn{1}{c|}{\cellcolor[HTML]{FFFFFF}{\color[HTML]{343434} Istotny}} \\ \hline
\cellcolor[HTML]{FFFFFF}{\color[HTML]{343434} }                         & {\color[HTML]{343434} $H_0$} & {\color[HTML]{343434} Naprawdę nieistotny}                                     & {\color[HTML]{343434} \textbf{Błąd typu I}}                                          \\
\multirow{-2}{*}{\cellcolor[HTML]{FFFFFF}{\color[HTML]{343434} Prawda}} & {\color[HTML]{343434} $H_A$} & {\color[HTML]{343434} \textbf{Błąd typu II}}                                            & {\color[HTML]{343434} Naprawdę istotny}                                     \\ \hline
\end{tabular}
\end{table}

W analizach statystycznych nie możemy kontrolować równocześnie błędu typu I i typu II. Tradycyjnie w testach statystycznych staramy się kontrolować błąd typu I. Wartość $\alpha$ w testach statystycznych oznacza poziom (maksymalny poziom) dla błędu typu I, jaki dopuszczamy przy testowaniu hipotezy. Wartość $\alpha$ jest zatem prawdopodobieństwem odrzucenia poprawnej hipotezy zerowej ($H_0$).  Zwyczajowo wartość ta ustalana jest ana poziomie 0.05. Oznacza to, że dopuszczamy odrzucenie prawidłowej hipotezy zerowej z prawdopodobieństwem wynoszącym 5%.

W przypadku porównań wielokrotnych, np. przy porównywaniu średnich konieczne jest zastosowanie odpowiednich metod, aby kontrolować poziom **błędu typu I** na odpowiednim poziomie, (poziomie zadanym przez $\alpha$).

### Test Tukeya

Odrzucenie hipotezy zerowej w analizie ANOVA pozwala nam stwierdzić, że istnieją różnice między średnimi w badanych grupach. Wynik ten nie mówi nam jednak, które z badanych średnich są różne. W celu określenia różnic pomiędzy średnimi musimy dokonać porównania (wielokrotnego) tychże średnich. W celu kontroli **błędu typu I** opracowana została metoda, zwana testem Tukeya. 

Przypuśćmy, że opracowaliśmy model w analizie ANOVA dany wzorem:

\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 analizie wariancji testujemy hipotezę postaci $H_0: \mu_1 = \mu_2 = ... = \mu_m$. Po obliczeniu statystyki testowej F i odrzuceniu hipotezy zerowej na poziomie istotnosci $\alpha$ chcemy stwierdzić średnie których z badanych grup są różne. Chcemy więc przeprowadzić serię "M wyborów 2"" hipotez danych jako: $H_0: \mu_i = \mu_j$.

Test HSD Tukeya (ang. *Tukey's honest significant difference test*) bazuje na rozkładzie "studentized range ($q$) distribution" - rozkładzie zbliżonym do rozkładu t studenta. Dla każdej hipotezy zerowej: $H_0: \mu_i = \mu_j$, statystyka testowa dana jest wzorem:
\begin{equation}
\label{eq:3}
t_s = \frac{\sqrt{2}\widehat{D}}{\sqrt{Var[\widehat{D}]}} ~ q_{m, (m \times n)-m}
\end{equation}

gdzie $\widehat{D}$ jest różnicą średnich grup $i$ oraz $j$. Wzór \ref{eq:3} można przekształcić do postaci:
\begin{equation}
\label{eq:4}
t_s = \frac{\sqrt{n}\widehat{D}}{\sqrt{MSE}} ~ q_{m, (m \times n)-m}
\end{equation}

Wzór ten jest bardzo podobny do wzoru na statystykę $t$ dla porównania dwóch średnich. Można wykazać, że $t_s = \sqrt{2}t$, gdzie $t$ jest statystyką testu dla dwóch grup o równej liczbie obserwacji. Poprawka dla wielu porównań wiąże się z liczbą stopni swobody dla rozkładu $q$ i statystyki $t_s$. Przy wzroście liczby porównań ($m$) roślinie wartość krytyczna $q_{m,(m \times n)-m}$ także rośnie, przez co trudniej jest odrzucić hipotezę zerową.

**Uwaga**

Test Tukeya należy do grupy testów post-hoc. Oznacza to, że test ten może zostać wykonany **tylko** po odrzuceniu hipotezy zerowej w analizie wariancji!

## Test Tukeya w R

Kontynujacja przykładu 1 z ćwiczeń 5.

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? Która pożywka jest najlepsza?

```{r, echo=FALSE}
dane <- read.csv('~/konas13@gmail.com/Dokumenty/Dydaktyka/Statystyka_Ochrona_Srodowiska/Cw_5_ANOVA/Cw5_1.txt', sep = "\t")
dane1 <- stack(dane)
colnames(dane1) <- c("liczba.kolonii","pozywka")
model <- lm(liczba.kolonii~pozywka,data=dane1)
anova(model)
```

Do wykonania testu Tukeya w R wykorzystujemy funkcję `TukeyHSD`.

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

Uzyskane wyniki (kolumna p adj) pozwalają stwierdzić, że statystycznie istotnie różnice występują jedynie pomiędzy pożywką C i pożywką A (p adj < $\alpha$). Nie ma istotnych różnic między pożywką A i B oraz B i C (p adj > $\alpha$).

**Ćwiczenie 1**

Proszę za pomocą testu t porównać średnie:
- z pożywki A i B
- z pożywki A i C
- z pożywki B i C

Czy uzyskane wyniki porównań testu t są zgodne z wynikami testu Tukeya?

*Podpowiedź*

W celu selekcji wyników dla pożywki A oraz B można wykorzystać następujący kod R:

```{r}
dane2 <- dane1[dane1$pozywka %in% c("M_A", "M_B"),]
dane2
```

**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. Które podłoże jest najlepsze?
