# Jednosmjerna analiza varijance

Neka uvodna priča... ^[samo primjer za fusnotu :)]

**PDF od F-distribucije**

```{r, warning = F, fig.width = 9, fig.height = 6, fig.cap="F-distribucije s različitim stupnjevima slobode"}
Fdist <- expand_grid(d1 = c(1,2,5,10,15,50), d2 = c(1,2,5,10,15,50)) %>% 
  group_by(d1,d2) %>%
  summarize(x = seq(0,5,0.01), fun = df(x, df1 = d1, df2 = d2), .groups = "drop")
Fdist %>%
  ggplot(aes(x = x, y = fun)) + 
  geom_line(color = "blue") + ylim(0,2.5) + 
  ylab("df(x)") + facet_grid(vars(d1), vars(d2), labeller = label_both)
```

## Malo matematike

Neka je $x_{ij}$ $j$-to opažanje u $i$-toj grupi, $\bar{x}_i$ aritmetička sredina $i$-te grupe i $\bar{x}$ aritmetička sredina svih opažanja. Tada je
$$x_{ij}=\bar{x}+(\bar{x}_i-\bar{x})+(x_{ij}-\bar{x}_i)$$
pri čemu je $\bar{x}_i-\bar{x}$ odstupanje aritmetičke sredine $i$-te grupe od aritmetičke sredine svih opažanja, a $x_{ij}-\bar{x}_i$ odstupanje $j$-tog opažanja iz $i$-te grupe od aritmetičke sredine te grupe.

Neka je $m$ ukupni broj grupa, a $n_i$ broj opažanja u $i$-toj grupi. **Varijacije unutar grupa** računaju se po formuli $$SSD_W=\sum_{i=1}^{m}{\sum_{j=1}^{n_i}{(x_{ij}-\bar{x}_i)^2}},$$
a **varijacije između grupa** po formuli
$$SSD_B=\sum_{i=1}^{m}{n_i(\bar{x}_{i}-\bar{x})^2}.$$

### Dodatna napomena {-}

Nenumerirana cjelina se označava s `{.unnumbered}` ili kraće s `{-}` na kraju imena u kodu.

## Primjer

```{r}
ToothGrowth$dose <- factor(ToothGrowth$dose)
a1_model <- aov(len ~ dose, data = ToothGrowth)
print(summary(a1_model), digits = 7)
```

**Ručni izračun**  

```{r}
Ybar <- mean(ToothGrowth$len)
tablica <- ToothGrowth %>% group_by(dose) %>% 
  mutate(Yk = mean(len), Nk = n()) %>% ungroup()
SSB <- tablica %>% summarize(SSB = sum((Yk - Ybar)^2)) %>% pull(SSB)
SSW <- tablica %>% summarize(SSW = sum((len - Yk)^2)) %>% pull(SSW)
dfB <- nlevels(factor(ToothGrowth$dose)) - 1
dfW <- nrow(tablica) - nlevels(factor(ToothGrowth$dose))
MSB <- SSB / dfB
MSW <- SSW / dfW
Fvalue <- MSB / MSW
pvalue <- pf(Fvalue, dfB, dfW, lower.tail = FALSE)
list(SSB = SSB, SSW = SSW, dfB = dfB, dfW = dfW, MSB = MSB, 
     MSW = MSW, Fvalue = Fvalue, pvalue = pvalue)
```

```{r, warning = F, fig.cap = "Prihvaćamo alternativnu hipotezu da postoji značajna razlika srednjih vrijednosti po grupama"}
tibble(x = seq(1,40,0.01)) %>%
  ggplot(aes(x)) + 
  geom_text(aes(x=Fvalue, label="F", y=0.1), nudge_x = 0.9, colour="red",
            text=element_text(size=11)) +
  geom_area(stat = "function", fun = df, args = list(df1 = dfB, df2 = dfW), fill = "pink", 
            xlim = c(qf(0.95, df1 = dfB, df2 = dfW), 40), alpha=0.75) +
  geom_vline(xintercept = Fvalue, linetype = "dotted") +
  stat_function(fun=df, args = list(df1 = dfB, df2 = dfW), color="blue", xlim=c(1,40)) +
  labs(title = "one way ANOVA")
```

**Effect size**

Trebamo učitati paket `lsr`.

```{r}
library(lsr)
```

Nakon toga možemo koristiti funkciju `etaSquared`.

```{r}
etaSquared(x = a1_model)
```

Isti rezultat možemo dobiti i ručnim izračunom.

```{r}
SSB / (SSB + SSW) 
```

## Post hoc tests

```{r}
pairwise.t.test( x = ToothGrowth$len, g = ToothGrowth$dose, p.adjust.method = "none")
```

```{r}
pairwise.t.test( x = ToothGrowth$len, g = ToothGrowth$dose, 
                 p.adjust.method = "bonferroni")
```

```{r}
pairwise.t.test( x = ToothGrowth$len, g = ToothGrowth$dose, p.adjust.method = "holm")
```

```{r}
TukeyHSD(a1_model)
```

## Levene test

Trebamo najprije učitati paket `car`.

```{r, message = F}
library(car)
```

Nakon toga možemo koristiti funkciju `leveneTest` za ispitivanje homogenosti varijanci po grupama.

```{r}
leveneTest(a1_model, center = mean)
```

Zapravo, Levene test radi anovu na varijabli $Z$ koja je definirana s $z_{ik}=|x_{ik}-\bar{x}_k|$.

```{r}
tablica1 <- tablica %>% mutate(Zk = abs(len - Yk))
tablica1
```

```{r}
hom <- aov(Zk ~ dose, data = tablica1)
print(summary(hom), digits = 7)
```

## Normalnost reziduala

```{r}
a1.residuals <- tibble(res = residuals(object = a1_model))
a1.residuals
```

**Histogram**

```{r}
ggplot(a1.residuals, aes(x=res)) + geom_histogram(binwidth=1, fill="pink", color="red")
```

**QQ-plot**

```{r}
ggplot(a1.residuals, aes(sample = res)) + stat_qq() + stat_qq_line()
```

**Shapiro-Wilk test za ispitivanje normalnosti**

Početna hipoteza: *"Populacija ima normalnu distribuciju."*

```{r}
shapiro.test(a1.residuals$res)
```

## Welchova ANOVA

Koristi se kada nije zadovoljen uvjet homogenosti varijanci grupa.

```{r}
oneway.test(len ~ dose, data = ToothGrowth)
```

Ako stavimo `var.equal = TRUE`, dobivamo standardnu 1-ANOVU.

```{r}
oneway.test(len ~ dose, data = ToothGrowth, var.equal = TRUE)
```

**Ručni izračun**

```{r}
shap <- ToothGrowth %>% group_by(dose) %>% 
  summarize(Xk = mean(len), v = var(len), n = n(), w = n / v)
J <- nlevels(ToothGrowth$dose)
U <- shap %>% summarize(U = sum(w)) %>% pull(U)
Xtilda <- shap %>% summarize(Xtilda = 1 / U * sum(w * Xk)) %>% pull(Xtilda)
A <- shap %>% 
  summarize(A = 1 / (J - 1) * sum(w * (Xk - Xtilda)^2)) %>% 
  pull(A)
B <- shap %>% 
  summarize(B = 2 * (J - 2) / (J^2 - 1) * sum((1 - w / U)^2 / (n - 1))) %>% 
  pull(B)
Fw <- A / (1+ B)
v1 <- J - 1
v2 <- shap %>% 
  summarize(v2 = (J^2 - 1) / (3 * sum((1 - w / U)^2 / (n - 1)))) %>% 
  pull(v2)
pvalue <- pf(Fw, df1=v1, df2=v2, lower.tail = FALSE)
list(tablica = shap, U = U, Xtilda = Xtilda, A = A, B = B, Fw = Fw, 
     df1 = v1, df2 = v2, pvalue = pvalue)
```


## Kruskal-Wallis test

Koristi se kada nije zadovoljen uvjet normalnosti.

```{r}
kruskal.test(len ~ dose, data = ToothGrowth)
```

**Ručni izračun**  

```{r}
tooth.rank <- ToothGrowth %>% mutate(Rk = rank(len))
kw.tab <- tooth.rank %>% group_by(dose) %>% 
  summarize(Rkbar = mean(Rk), Nk = n()) %>% ungroup()
Rbar <- mean(tooth.rank$Rk)
RSS.tot <- tooth.rank %>% summarize(RSS.tot = sum((Rk - Rbar)^2)) %>% pull(RSS.tot)
RSS.b <- kw.tab %>% summarize(RSS.b = sum(Nk * (Rkbar - Rbar)^2)) %>% pull(RSS.b)
N <- sum(kw.tab$Nk)
K <- (N - 1) * RSS.b / RSS.tot
dfB <- nlevels(tooth.rank$dose) - 1
pvalue <- pchisq(K, dfB, lower.tail = FALSE)
list(RSS.tot = RSS.tot, RSS.b = RSS.b, dfB = dfB,
     K = K, pvalue = pvalue)
```

```{r, warning = F, fig.cap = "Prihvaćamo alternativnu hipotezu da postoji značajna razlika srednjih vrijednosti po grupama"}
tibble(x = seq(0,50,0.01)) %>%
  ggplot(aes(x)) + 
  geom_text(aes(x=K, label="K", y=0.1), nudge_x = 0.9, colour="red",
            text=element_text(size=11)) +
  geom_area(stat = "function", fun = dchisq, args = list(df = dfB), fill = "pink", 
            xlim = c(qchisq(0.95, df = dfB), 40), alpha=0.75) +
  geom_vline(xintercept = K, linetype = "dotted") +
  stat_function(fun=dchisq, args = list(df = dfB), color="blue", xlim=c(0,50)) +
  labs(title = "Kruskal-Wallis test")
```

**Dunn test (post hoc test)**

**Link**  
[Dunn test](https://stats.stackexchange.com/tags/dunn-test/info){target="_blank"}  

```{r warning = F, message=F}
library(FSA)
```

```{r warning = F}
dunnTest(len ~ dose, data=ToothGrowth, method="bonferroni")
```

```{r warning = F}
library(dunn.test)
dunn.test(ToothGrowth$len, ToothGrowth$dose, altp = TRUE, table = FALSE, list = TRUE)
```

```{r warning = F}
library(dunn.test)
dunn.test(ToothGrowth$len, ToothGrowth$dose, table = FALSE, list = TRUE)
```