---
title: "Algemeen lineair model met interactie: uitgewerkt voorbeeld Dojovissen"
output:
  html_document:
    code_download: yes
    highlight: tango
    number_sections: yes
    theme: cosmo
    toc: yes
    toc_float: yes
  pdf_document:
    toc: yes
  word_document:
    toc: yes
---

# Achtergrond

Resistentie tegen het gif EI-43,064 wordt getest bij 96 vissen (dojovissen (0), goudvissen (1) en zebravissen (2)).
Elke vis wordt apart in een aquarium gestopt die een bepaalde `dosis` (in mg) van het gif bevat.
Naast de overlevingstijd in minuten (`minsurv`) werd ook het `gewicht` van de vis gemeten (in gram).
De onderzoekers weten uit vorige experimenten dat de overlevingstijd vaak sterk afhangt van het gewicht en dat de resistentie dikwijls soortafhankelijk is.
De onderzoekers wensen inzicht te krijgen in het effect van de dosis op de overlevingstijd en of resistentie tegen het gif verschillend is bij de verschillende soorten.


## Laad de libraries

```{r}
library(dplyr)
library(ggplot2)
#install.packages("GGally")
library(GGally)
library(car)
library(multcomp)
```

# Data-exploratie

Lees de dataset poison.dat in via `read.table`.
We zagen in de voorgaande practica al dat de overlevingstijd beter op log2-schaal wordt gemodelleerd.

```{r}
poison <- read.table("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum8/poison.dat", sep="", header = TRUE)

# We vormen de vissoort om in een factor en log2 transformeren de overlevingstijd
poison <- poison %>%
  mutate(soort = as.factor(soort),
         log2minsurv = log2(minsurv))

nobs <- sum(poison$soort == 0)

poison %>%
  ggpairs
```

- De overlevingstijd lijkt geassocieerd met het gewicht, soort en de dosis.

- We observeren een sterke positieve associatie tussen de log2-overlevingstijd en het gewicht.

- Bij lage gewichten lijkt de log2-overlevingstijd wat af te vlakken.

- Daarnaast zien we ook dat het gewicht niet gelijk verdeeld is binnen elke dosis.

- Er lijkt ook een associatie tussen de gewicht en soort.

We bestuderen de effecten verder:

```{r}
poison %>%
  ggplot(aes(x = dosis, y = log2minsurv, col = soort)) +
  geom_point() +
  stat_smooth(method = "loess",col="red") + # fit een kromme door de puntenwolk (rode lijn)
  ylab("Overlevingstijd (log2)") +
  xlab("Dosis (mg)") +
  theme_bw()

poison %>%
  ggplot(aes(x = gewicht, y = log2minsurv, col = soort)) +
  geom_point() +
  stat_smooth(method = "loess",col="red") + # fit een kromme door de puntenwolk (rode lijn)
  ylab("Overlevingstijd (log2)") +
  xlab("Gewicht") +
  theme_bw()

scatterplot1 <- poison %>%
  ggplot(aes(x = dosis, y = gewicht)) +
  geom_point() +
  theme_bw()
scatterplot1

boxplot2 <- poison %>% ggplot(
    aes(x = soort,
        y = log2minsurv,
        col = soort)
    ) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter() +
  theme_bw()
boxplot2

boxplot3 <- poison %>% ggplot(
    aes(x = soort,
        y = gewicht,
        col = soort)
    ) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter() +
  theme_bw()
boxplot3
```

- We zien een heel sterke relatie tussen het gewicht en de log2-overlevingstijd. Het gewicht is eveneens ongelijk tussen de soorten.

- We observeren dat gewichten niet uniform verdeeld zijn over de verschillende dosissen.

- Het patroon wordt eveneens weerspiegeld in de plot waarbij de log2survival wordt uitgezet t.o.v. dosis.

- Deze plot geeft ook weer dat er een effect is van de soort op de overleving.

Gezien de onderzoekers op basis van voorgaande studies vermoeden dat het effect van de dosis kan variëren van soort tot soort zouden we in de modellen interacties moeten voorzien tussen dosis en soort, zodat elke soort een verschillende dosisrespons kan tonen. Daarom modelleren we momenteel enkel voor dojovissen.

***Vorige oefening hebben we enkel de hoofdeffecten in het model opgenomen, zonder interactietermen. Dit was louter voor didactische redenen. Normaal gezien moet je altijd eerst de interactie meenemen in het model, en bekijken of deze significant is. Indien je deze interactie niet weerhoudt, dan pas kan je een model opstellen met enkel de hoofdeffecten. We verwachten dus, wanneer je een analyse doet, dat je de methode van deze file gebruikt, in plaats van het lineair additief model.***

# Schat het effect van de dosis na correctie voor het gewicht bij dojovissen met een mogelijkse interactie tussen dosis en gewicht.

## Model

$$
y_i=\beta_0+\beta_d x_d + \beta_g x_g +\beta_{d:g} x_d x_g+ \epsilon_i,
$$
met $\epsilon_i \text{ i.i.d. } N(0,\sigma^2)$

We kunnen de volgende interpretaties geven aan de geschatte coëffiënten:

- Voor een vis met gewicht $x_g$, zal de verwachte log-overlevingstijd $\beta_d+\beta_{d:g}*x_g$ hoger liggen, wanneer de dosis met eenheid verhoogt wordt, ten opzichte van een vis met hetzelfde gewicht. 

- Voor een vis die dosis $x_d$ kreeg, zal de verwachte log-overlevingstijd $\beta_g+\beta_{d:g}*x_d$ hoger liggen, wanneer het gewicht van die vis eenheid hoger is, vergeleken met een vis die dezelfde dosis kreeg.

De parameter $\beta_{d:g}$ kan je dus op twee manieren interpreteren. De meest voor de hand liggende is de volgende: "Het geeft aan hoe de invloed van dosis op de log-overlevingstijd afhankelijk is van het gewicht van een vis". Maar je kan ook zeggen: "Het geeft aan hoe de invloed van het gewicht op de log-overlevingstijd afhankelijk is van de dosis".

We schatten de parameters van het model, evalueren de aannames en doen besluitvorming in R:

```{r}
lmDojoInt <- lm(log2minsurv ~ dosis * gewicht, data = poison %>%filter(soort == 0))
summary(lmDojoInt)
plot(lmDojoInt)
```

## Aannames

- We zien geen problemen m.b.t. lineariteit.

- We zien geen problemen mbt de aanname van gelijkheid van variantie.

- De data lijken wat afwijkingen van normaliteit te vertonen in de linkerstaart. De afwijkingen zijn echter niet extreem en bovendien beschikken we reeds over `r nobs` observaties zodat we er zeker van uit mogen gaan dat de parameterschatters benaderend normaal verdeeld zijn.


## Inferentie

Het dosiseffect wordt geparameteriseerd door meerdere modelparameters.
We evalueren eerst de omnibushypothese dat er geen effect is van de dosis: dus geen hoofdeffect en geen interactie.
Dat kan d.m.v. een F-test waarbij we het model vergelijken met hoofdeffect voor dosis en interactie tussen dosis en gewicht, met een model waarin enkel gewicht is opgenomen.

```{r}
lmDojoH0 <- lm(log2minsurv ~ gewicht,
               data = poison %>%
                filter(soort == 0))
anova(lmDojoH0,lmDojoInt)
```

We zien dat er een extreem effect is van de dosis op de overleving bij vissen.

## Conventionele aanpak

We hebben reeds gezien dat er een significant dosis effect is.

We gaan nu eerst na of er een interactie is.
Omdat er maar 1 interactie parameter is kan dat d.m.v.

1. Summary functie
2. F-test waarbij we model met en zonder interactie vergelijken
3. anova tabel met type III kwadratensommen  


```{r}
summary(lmDojoInt)
lmDojo <- lm(log2minsurv ~ dosis + gewicht,
             data = poison %>%
              filter(soort == 0))
anova(lmDojo,lmDojoInt)
Anova(lmDojoInt,type="III")
```

Er is geen significante interactie tussen de dosis en het gewicht. Het effect van de dosis op de overleving van vissen verandert dus niet significant in functie van het gewicht van de vis.

De conventionele analyse is om de interactie uit het model te weren.
We hebben dit reeds gedaan a.d.h.v model lmDojo.

```{r}
summary(lmDojo)
```

### Conclusie

We krijgen, doordat we de interactie weerden uit het model, hetzelfde model als bij het additief lineair model. Daarom komt de conclusie overeen met de vorige conclusie:

Het gif heeft een extreem significant effect op de overlevingstijd bij dojovissen ($p<<0.001$).
Het geometrisch gemiddelde van de overlevingstijd van dojovissen wordt gehalveerd (factor $2^{\beta_d}=$ `r format(2^(lmDojo$coef["dosis"]),digits=3)`) voor vissen met hetzelfde gewicht als de dosis van het gif met 1 mg/l toeneemt (95% BI [`r paste(format(2^(confint(lmDojo)["dosis",]),digits=3),collapse=",")`]).

De verandering van het effect van de dosis op de overleving in functie van het gewicht van de vis is niet significant (p = `r format(anova(lmDojo,lmDojoInt)[2,"Pr(>F)"], digits=2)`).

## Alternatieve aanpak zonder het vereenvoudigen van het model

We hebben reeds gezien dat er een significant dosis effect is.
We testen nu de interactie.


```{r}
summary(lmDojoInt)
```

Er is geen significante interactie. Het aanvaarden van de nulhypothese is een zwakke conclusie.
We weten ook dat de power op het vinden van een interactie eerder laag is.
Daarom kunnen we opteren om de interactie in het model te laten.


We kunnen echter wel een uitspraak doen over het gemiddelde dosis effect bij de Dojo vissen in het experiment d.m.v. te marginaliseren over alle Dojo vissen in het experiment.

$$
\frac{\sum\limits_{i=1}^{n_\text{dojo}} (\beta_d + \beta_{d:w} x_{iw})}{n_\text{dojo}} = \beta_d + \beta_{d:w} \bar{x}_w
$$

We berekenen m.a.w. het effect bij de gemiddelde dojovis in het experiment.

```{r}
wBar <- poison %>%
  filter(soort == 0) %>%
  pull(gewicht) %>%
  mean
wBar
```


Het gemiddeld gewicht voor een dojo-vis in het experiment is
`r round(wBar, 2)` g.


We evalueren nu het contrast voor het marginale effect. We testen dus of $\beta_d+2.1313*\beta_{d:g}$ significant verschilt van nul.


```{r}
marginaleDosisEffect <- glht(lmDojoInt, linfct = c("dosis+ 2.1313*dosis:gewicht = 0"))
summary(marginaleDosisEffect)
confint(marginaleDosisEffect)
2^confint(marginaleDosisEffect)$confint
```

Merk op dat het marginale effect van de dosis bijna exact gelijk is aan het effect die wordt geschat in het model zonder interactie.

### Conclusie

Er is een extreem significant effect van de dosis op de overleving bij vissen ($p<<0.001$).

De verandering van het effect van de dosis op de overleving  in functie van het gewicht van de vis is niet significant (p = `r format(anova(lmDojo,lmDojoInt)[2,"Pr(>F)"], digits=2)`).

Het marginaal geometrisch gemiddelde van de overlevingstijd halveert ongeveer (factor `r format( 2^confint(marginaleDosisEffect)$confint[1,1], digits=3)`) bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [`r paste(format(2^confint(marginaleDosisEffect)$confint[1,-1],digits=3),collapse=", ")`]), na correctie voor gewicht ($p<<0.001$).



