# Dvosmjerna analiza varijance

Nekoliko uvodnih riječi...

Nekakvi tekst pa nabrajanje

- opcija 1 *tekst*
- opcija 2
- opcija 3

Neka fusnota^[ovdje je opis u fusnoti]

## Primjer (bez interakcije)

Primjer u kojemu pretpostavljamo da nema interakcije između dvije faktorske varijable.

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

**Ručni izračun**  

```{r}
tablica_ab <- ToothGrowth %>%
  group_by(supp) %>% mutate(Ya = mean(len)) %>% ungroup() %>%
  group_by(dose) %>% mutate(Yb = mean(len)) %>% ungroup()
Ybar <- mean(ToothGrowth$len)
SSA <- tablica_ab %>% summarize(SSA = sum((Ya - Ybar)^2)) %>% pull(SSA)
SSB <- tablica_ab %>% summarize(SSB = sum((Yb - Ybar)^2)) %>% pull(SSB)
SST <- tablica_ab %>% summarize(SST = sum((len - Ybar)^2)) %>% pull(SST)
SSR <- SST - SSA - SSB
dfA <- nlevels(ToothGrowth$supp) - 1
dfB <- nlevels(ToothGrowth$dose) - 1
dfR <- nrow(ToothGrowth) - 1 - dfA - dfB
MSA <- SSA / dfA
MSB <- SSB / dfB
MSR <- SSR / dfR
FA <- MSA / MSR
FB <- MSB / MSR
pA <- pf(FA, df1=dfA, df2=dfR, lower.tail=FALSE)
pB <- pf(FB, df1=dfB, df2=dfR, lower.tail=FALSE)
list(SSA = SSA, SSB = SSB, SSR = SSR, dfA = dfA, dfB = dfB, dfR = dfR, 
     MSA = MSA, MSB = MSB, MSR = MSR, FA = FA, FB = FB, 
     pA = pA, pB = pB)
```

**Effect size**

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

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

```{r}
list(eta.sq.SUPP = SSA / SST, eta.sq.part.SUPP = SSA / (SSA + SSR),
     ets.sq.DOSE = SSB / SST, eta.sq.part.DOSE = SSB / (SSB + SSR))
```

## Primjer (s interakcijom)

Primjer u kojemu pretpostavljamo da ima interakcije između dvije faktorske varijable.

```{r, warning = F}
ToothGrowth %>% 
  ggplot(aes(x = supp, y=len, color = dose, group = dose)) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line")
```

```{r, warning = F}
ToothGrowth %>% 
  ggplot(aes(x = dose, y=len, color = supp, group = supp)) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line")
```

```{r}
a2_inter <- aov(len ~ supp * dose, data = ToothGrowth)
print(summary(a2_inter), digits = 7)
```

**Ručni izračun**  

```{r}
tablica_ab_inter <- ToothGrowth %>% group_by(supp, dose) %>% 
  mutate(Yab = mean(len)) %>% ungroup %>%
  group_by(supp) %>% mutate(Ya = mean(len)) %>% ungroup() %>%
  group_by(dose) %>% mutate(Yb = mean(len)) %>% ungroup()
Ybar <- mean(ToothGrowth$len)
SSA <- tablica_ab_inter %>% summarize(SSA = sum((Ya - Ybar)^2)) %>% pull(SSA)
SSB <- tablica_ab_inter %>% summarize(SSB = sum((Yb - Ybar)^2)) %>% pull(SSB)
SSAB <- tablica_ab_inter %>% summarize(SSAB = sum((Yab - Ya - Yb + Ybar)^2)) %>% 
  pull(SSAB)
SST <- tablica_ab_inter %>% summarize(SST = sum((len - Ybar)^2)) %>% pull(SST)
SSR <- SST - SSA - SSB - SSAB
dfA <- nlevels(ToothGrowth$supp) - 1
dfB <- nlevels(ToothGrowth$dose) - 1
dfAB <- dfA * dfB
dfR <- nrow(ToothGrowth) - (dfA + 1) * (dfB + 1)
MSA <- SSA / dfA
MSB <- SSB / dfB
MSAB <- SSAB / dfAB
MSR <- SSR / dfR
FA <- MSA / MSR
FB <- MSB / MSR
FAB <- MSAB / MSR
pA <- pf(FA, df1=dfA, df2=dfR, lower.tail=FALSE)
pB <- pf(FB, df1=dfB, df2=dfR, lower.tail=FALSE)
pAB <- pf(FAB, df1=dfB, df2=dfR, lower.tail=FALSE)
list(SSA = SSA, SSB = SSB, SSAB = SSAB, SSR = SSR, 
     dfA = dfA, dfB = dfB, dfAB = dfAB, dfR = dfR, 
     MSA = MSA, MSB = MSB, MSAB = MSAB, MSR = MSR, 
     FA = FA, FB = FB, FAB = FAB, 
     pA = pA, pB = pB, pAB = pAB)
```

**Effect size**

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

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

```{r}
list(eta.sq.SUPP = SSA / SST, eta.sq.part.SUPP = SSA / (SSA + SSR),
     ets.sq.DOSE = SSB / SST, eta.sq.part.DOSE = SSB / (SSB + SSR),
     eta.sq.SUPP.DOSE = SSAB / SST, eta.sq.part.SUPP.DOSE = SSAB / (SSAB + SSR))
```

## Levene test

This function expects that you have a saturated model (i.e., included all of the 
relevant terms), because the test is primarily concerned with the within-group 
variance, and it doesn’t really make a lot of sense to calculate this any way 
other than with respect to the full model.

```{r, error=TRUE}
leveneTest(a2_model, center = mean)
```

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

## Normalnost reziduala

Pogledajmo reziduale za model s interakcijom. 

```{r}
a2.residuals <- tibble(res = residuals(object = a2_inter))
a2.residuals
```

**Histogram**

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

**QQ-plot**

```{r}
ggplot(a2.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(a2.residuals$res)
```

## Biranje modela

Bolji model u ovom slučaju je model s uključenom interakcijom.

```{r}
anova(a2_model, a2_inter)
```

**Ručni izračun**

* Stupnjevi slobode definirani su u smislu broja parametara koji se moraju procijeniti u modelu.
* F-test je uspoređivanje dva linearna modela.
* Prvi stupanj slobobde $\mathrm{d}f_1$ u F-testu jednak je razlici broja parametara koje treba procijeniti
u dva modela koji se uspoređuju.
* Drugi stupanj slobode $\mathrm{d}f_2$ u F-testu odnosi se na stupnjeve slobode povezane s rezidualima.
To je zapravo jednako razlici ukupnog broja opažanja (broj ispitanika) i broja parametara koje treba procijeniti.

U slučaju 2-anove, ukoliko faktorske varijable imaju redom $a$ i $b$ faktora, tada se svaka od njih može
redom reprezentirati s $a-1$ i $b-1$ binarnih varijabli. Za model bez interakcije treba procijenti ukupno 
$(a-1)+(b-1)+1=a+b-1$ parametara (parametar uz svaku binarnu varijablu i slobodni član).  
Za puni model sa svim interakcijama treba procijeniti ukupno $ab$ parametara.
$$ab=((a-1)+1)((b-1)+1)=(a-1)+(b-1)+(a-1)(b-1)+1$$
Naime, uz svaku binarnu varijablu imamo $(a-1)+(b-1)$ parametara, uz sve moguće interakcije imamo
$(a-1)(b-1)$ parametara i konačno još jedan parametar za slobodni član.

```{r}
a2_model_tidy <- tidy(a2_model)
a2_inter_tidy <- tidy(a2_inter)
a2_model_tidy
a2_inter_tidy
```

```{r}
a <- nlevels(ToothGrowth$supp)
b <- nlevels(ToothGrowth$dose)
SSR.null <- a2_model_tidy %>% pull(sumsq) %>% last()
SSR.full <- a2_inter_tidy %>% pull(sumsq) %>% last()
SSR.diff <- SSR.null - SSR.full
df1 <- a * b - (a + b - 1)
df2 <- nrow(ToothGrowth) - a * b
MS.diff <- SSR.diff / df1
MS.res <- SSR.full / df2
Fv <- MS.diff / MS.res
pvalue <- pf(Fv, df1, df2, lower.tail = FALSE)
list(SSR.null = SSR.null, SSR.full = SSR.full, SSR.diff = SSR.diff,
     df1 = df1, df2 = df2, Fv = Fv, pvalue = pvalue)
```

Zapravo sve se radi preko lineranog modela i `lm` funkcije.

```{r}
lm.1 <- lm(len ~ supp + dose, data = ToothGrowth)
lm.2 <- lm(len ~ supp * dose, data = ToothGrowth)
sum(residuals(lm.1)^2)
sum(residuals(lm.2)^2)
anova(lm.1,lm.2)
```

## Tukey test

```{r}
TukeyHSD(a2_model)
```

```{r}
TukeyHSD(a2_inter)
```

## Nebalansirani dizajn

[Strategije link](https://md.psych.bio.uni-goettingen.de/mv/unit/lm_cat/lm_cat_unbal_ss_explained.html){target="_blank"} 

```{r}
load("coffee.Rdata")
coffee %>% kable()
coffee %>% count(milk, sugar)
```

### Strategija tipa I

* Strategija tipa I gradi model postupno, počevši od najjednostavnijeg mogućeg modela i postupnog dodavanja novih varijabli.
* Nedostatak: ovisi o redoslijedu dodavanja novih varijabli u model.

```{r}
model.1 <- lm(babble ~ sugar + milk + sugar:milk, coffee)
model.2 <- lm(babble ~ milk + sugar + sugar:milk, coffee)
anova(model.1)
anova(model.2)
```

### Strategija tipa III

* Uvijek pokreće F-test s alternativnom hipotezom koja je jednaka punom modelu kojeg ispitujemo.
* Nulta hipoteza je jednaka modelu kojeg dobijemo tako da iz punog modela maknemo varijablu koju ispitujemo.
* Ovdje nemamo problem s redoslijedom varijabli koji postoji u strategiji tipa I.
* Međutim, ova strategija ovisi o načinu na koji se reprezentiraju faktorske varijable.

Za strategije tipa III i tipa II treba koristiti funkciju `Anova` iz `car` biblioteke.

```{r}
Anova(model.1, type = 3)
Anova(model.2, type = 3)
```

Ukoliko promijenimo način reprezentacije faktorskih varijabli, rezultati se mogu značajno promijeniti.

```{r}
model.H <- lm( babble ~ sugar * milk, coffee, 
               contrasts = list(milk = "contr.Helmert", sugar = "contr.Helmert"))
Anova(model.H, type = 3)
```

### Strategija tipa II

* Uklanja probleme koji se javljau u prethodne dvije strategije pa je lakše interpretirati rezultate.
* Strategija je u potpunosti slična strategiji tipa III uz dodatni uvjet da se poštuje **načelo marginalnosti**.
* Načelo marginalnosti zabranjuje izostavljanje varijabli nižeg reda iz modela ako postoje varijable višeg reda
koje ovise o njima.

```{r}
Anova(model.1, type = 2)
Anova(model.2, type = 2)
Anova(model.H, type = 2)
```
