---
title: "Algemeen lineair model: half uitgewerkte oefening"
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, libraries laden en data-exploratie

We herhalen aan het begin van deze oefening nogmaals de situatie en bijhorende data-exploratie. Net zoals in de vorige oefening zullen jullie hier een additief model opstellen, waarbij het effect van de dosis na correctie voor het gewicht bij **goudvissen** (soort 1) geschat zal worden.

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 (de uitkomst, `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.


```{r}
#gebruik de install.packages() functie voor packages die nog niet geïnstalleerd zijn.
library(dplyr)
library(ggplot2)
library(GGally)
library(car)
library(multcomp)
```


```{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))
```

In de vorige oefening onderzochten we de invloed van dosis en gewicht op de overlevingstijd voor goudvissen We zullen nu een gelijkaardige analyse uitvoeren, maar we zullen deze toepassen op de goudvissen (soort==1).

# Vraag 1 Aanpak van vorig practicum met enkelvoudige regressie (leesopdracht)

Herinner het model met enkel dosis-effect uit het practicum enkelvoudige lineaire regressie, hier enkel voor de goudvissen:

$$
y_i=\beta_0+\beta_d x_d + \epsilon_i,
$$

```{r}
lm1 <- lm(log2minsurv~dosis,
          data = poison %>%
            filter(soort == 1)) # Filter op goudvissen
summary(lm1)
```

# Vraag 2: Aanpassen model door correctie voor gewicht bij goudvissen

We willen nu een extra predictor (onafhankelijke variabele of covariaat) toevoegen aan het model, namelijk 'gewicht'.

$$
y_i=\beta_0+\beta_d x_d + \beta_g x_g + \epsilon_i,
$$

met $\epsilon_i \text{ i.i.d. } N(0,\sigma^2)$.

Corrigeer het enkelvoudige lineair model model `lm1`uit vraag 1 voor het gewicht bij de goudvissen.

```{r}
lmGoud <- lm(... ~ ... ,
             data = poison %>%
              filter(...))
summary(...)
```



# Vraag 3: Ga de aannames voor lineaire regressie na aan de hand van de `plot()`functie.

```{r}
plot(...)
```


# Vraag 4: Lineariteit (meerkeuzevraag)

Op basis van de plot uit vraag 3 kunnen we stellen dat 

1. De lineariteit niet voldaan is omdat er veel variatie zit op de gefitte waarden.
2. De lineariteit niet voldaan is omdat de rode smoother een duidelijke niet-lineaire trend vertoont.
3. De lineariteit wel voldaan is.
4. De lineariteit moeilijk te beoordelen is op basis van de plots.

# Vraag 5: Gelijke varianties (meerkeuzevraag)

Op basis van de plot uit vraag 3 kunnen we stellen dat 

1. De aanname van gelijke varianties niet voldaan is omdat er veel variatie zit op de gefitte waarden.
2. De aanname van gelijke varianties niet voldaan is omdat de puntenwolk van residuen een duidelijk patroon vertoont waarbij er een lage variantie zit bij lage en hoge gefitte waarden, maar een hoge variantie bij tussenliggende gefitte waarden. 
3. De aanname van gelijke varianties wel voldaan is.
4. De aanname van gelijke varianties moeilijk te beoordelen is op basis van de plots.

# Vraag 6: Normaliteit (meerkeuzevraag)

Op basis van de plot uit vraag 3 kunnen we stellen dat 

1. De aanname van normaliteit niet voldaan is.
2. De aanname van normaliteit voldaan is.
3. De aanname van normaliteit moeilijk te beoordelen is op basis van slechts 38 datapunten.

# Vraag 7: Hoe verandert het effect van dosis na correctie voor gewicht bij de goudvissen?

1. Het effect van dosis was al significant wanneer niet werd gecorrigeerd voor gewicht, maar is nog significanter na correctie voor gewicht.
2. Het effect van dosis na correctie voor gewicht blijft gelijk.
3. Het effect van dosis na correctie voor gewicht wordt minder sterk (p-waarde van 0.721 i.p.v. 2.43e-10)

# Vraag 8: Effect van dosis na correctie voor gewicht (leesopdracht)
We observeren dat het effect van de dosis ook bij de goudvissen significanter is na correctie voor gewicht. Hierbij geldt dezelfde redenering als bij het model met enkel goudvissen (herhaling):

1. Dat is logisch want uit de data-exploratie bleek dat gewicht geassocieerd was met de overleving.
Wanneer we een extra variabele opnemen in het model die veel van de variabiliteit in de uitkomstvariabele kan verklaren, zal de residuele variabiliteit afnemen, waardoor de standaard errors op de parameterschatters sterk af kunnen nemen.

2. Bovendien observeerden we in de patronen van gewicht en dosis, gewicht en overlevingstijd en dosis en overlevingstijd ook dat de gewichten niet goed gebalanceerd zijn tussen de verschillende dosissen waardoor het effect van gewicht de schatting van het dosiseffect wat verstoort. We zullen het effect van de dosis dus beter kunnen schatten als we corrigeren voor het gewicht.

Uit (1) en (2) volgt dus dat als we corrigeren voor gewicht dat we het dosiseffect dus beter en met een grotere precisie kunnen schatten!

# Vraag 9: Schatting van de log2-overlevingstijd bij goudvissen

Scaht de log2-overlevingstijd bij goudvissen op basis van het opgestelde model `lmGoud` in vraag 2.

1. De log2-overlevingstijd bij goudvissen is gemiddeld 1.556 lager bij vissen die blootgesteld worden aan een dosis van het gif dat 1mg/l hoger is, na correctie voor gewicht.
2. De log2-overlevingstijd bij goudvissen is gemiddeld 1.556 hoger bij vissen die blootgesteld worden aan een dosis van het gif dat 1mg/l hoger is, na correctie voor gewicht.
3. De log2-overlevingstijd bij goudvissen is gemiddeld 0.790 lager bij vissen die blootgesteld worden aan een dosis van het gif dat 1mg/l hoger is, na correctie voor gewicht.
4. De log2-overlevingstijd bij goudvissen is gemiddeld 0.790 hoger bij vissen die blootgesteld worden aan een dosis van het gif dat 1mg/l hoger is, na correctie voor gewicht.



# Vraag 10: Terugtransformeren

Transformeer de data terug zodat een interpretatie in termen van geometrische gemiddelde mogelijk is.
```{r}
coef_transf <- 2^(...) # Terugtransformeren van de regressiecoëfficiënten
conf_transf <- 2^(confint(...)) # Terugtransformeren van de betrouwbaarheidsintervallen
```

# Vraag 11: Conclusie effect dosis

Geef aan welke conclusie correct is met betrekking van het effect van dosis op overlevingstijd op basis van de voorgaande analyse.

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


# Vraag 12: Conclusie effect gewicht

Geef aan welke conclusie correct is met betrekking van het effect van gewicht op overlevingstijd op basis van de voorgaande analyse.

1. Het effect van het gewicht op de overlevingstijd is extreem significant ($p<<0.001$).
Het geometrisch gemiddelde van de overlevingstijd van een goudvis die een 1 gram meer weegt dan een andere goudvis en die aan dezelfde dosis wordt blootgesteld, is dubbel zo groot (factor $2^{ \beta_g}$= `r format(2^(lmGoud$coef["gewicht"]),digits=3)`, 95% BI [`r paste(format(2^(confint(lmGoud)["gewicht",]),digits=3),collapse=",")`]).

2. Het effect van het gewicht op de overlevingstijd is extreem significant ($p<<0.001$).
Het geometrisch gemiddelde van de overlevingstijd van een goudvis die een 1 gram meer weegt dan een andere goudvis en die aan dezelfde dosis wordt blootgesteld, is drie keer zo groot (factor $2^{ \beta_g}$= `r format(2^(lmGoud$coef["gewicht"]),digits=3)`, 95% BI [`r paste(format(2^(confint(lmGoud)["gewicht",]),digits=3),collapse=",")`]).

3. Het effect van het gewicht op de overlevingstijd is extreem significant ($p<<0.001$).
Het geometrisch gemiddelde van de overlevingstijd van een goudvis die een 1 gram meer weegt dan een andere goudvis is dubbel zo groot (factor $2^{ \beta_g}$= `r format(2^(lmGoud$coef["gewicht"]),digits=3)`, 95% BI [`r paste(format(2^(confint(lmGoud)["gewicht",]),digits=3),collapse=",")`]).

4.  Het effect van het gewicht op de overlevingstijd is extreem significant ($p<<0.001$).
Het geometrisch gemiddelde van de overlevingstijd van een goudvis die een 1 gram meer weegt dan een andere goudvis is drie keer zo groot (factor $2^{ \beta_g}$= `r format(2^(lmGoud$coef["gewicht"]),digits=3)`, 95% BI [`r paste(format(2^(confint(lmGoud)["gewicht",]),digits=3),collapse=",")`]).


# Vraag 13: Model met interactietermen

Stel het model op met interactietermen.
```{r}
lmGoudInt <- lm(..., data = ... %>%filter(...))
```

# Vraag 14: Ga de voorwaarden van het model na (Meerkeuzevraag)

```{r}
plot(lmGoudInt)
```

In het vorige model was aan alle voorwaarden voldaan. Welke van de voorwaarden is niet meer voldaan voor dit model?

1. Onafhankelijke gegevens
2. Lineariteit
3. Normaal verdeelde residuen
4. Geen van bovenstaande

# Vraag 15: Interpretatie (Meerkeuzevraag)

Welke interpretatie van de schattingen is juist
```{r}
summary(lmGoudInt)
```
1. De verwachte log2-overlevingstijd van een vis van 2 kg zal gemiddeld 1.66 lager liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt.

2. De verwachte log2-overlevingstijd van een vis van 2 kg zal gemiddeld 0.37 hoger liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt.

3. De verwachte log2-overlevingstijd van een vis van 2 kg zal gemiddeld 0.51 lager liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt.

4. De verwachte log2-overlevingstijd van een vis van 2 kg zal gemiddeld 0.66 lager liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt.

# Vraag 16: Test of er een effect is van dosis.

Test of er een effect is van dosis (hoofdeffect en interactie-effect samen). Stel daarvoor eerst een lineair model op zonder dosis en vergelijk dit model met het lmGoudInt model via de anova-functie

```{r}
lmGoudH0 <- lm(..., data = ... %>%filter(...))
tabel_anova<-anova(...)
tabel_anova
```
# Vraag 17: Interactie testen

Test of er een interactie-effect aanwezig is. Geef ook de tabel Type III kwadratensommen.

```{r}
tabel_kwadratensommen<-...
tabel_kwadratensommen
```

# Vraag 18: Bereken het gemiddelde gewicht van de goudvissen

```{r}
wBar <- ...
wBar
```

# Vraag 19: Marginale dosis-effect (Leesopdracht)

Hoewel de interactieterm niet significant is, zullen we deze toch in ons model laten. We zullen het marginale dosis-effect schatten. Dit is het gemiddelde effect over alle vissen heen en is gelijk aan het effect voor een vis met gemiddeld gewicht, m.a.w. $\beta_d +\beta_{d:g}*\bar{x}_g$  

Eerst kunnen we dit doen op basis van de coëfficiënten die we vinden in de summary:
```{r}
summary(lmGoudInt)
0.3691+wBar*(-0.5146)
```

# Vraag 20: Schat nu het marginaal effect met behulp van de glht-functie en geef de summary met de p-waarde.
```{r}
marginaleDosisEffect <- glht(..., linfct = c("..."))
summary1<-...
summary1
```


# Vraag 21: Geef de geschatte coëfficienten met betrouwbaarheidsinterval op originele schaal (niet log-getransformeerd) en interpreteer. (Meerkeuzevraag)
```{r}
2^confint(marginaleDosisEffect)$confint
```

1. De gemiddelde overlevingstijd zal 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$).

2. Het geometrisch gemiddelde van de overlevingstijd zal een factor 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$).

3. Het marginaal gemiddelde van de overlevingstijd zal 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$).

4. Het marginaal geometrisch gemiddelde van de overlevingstijd zal een factor 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$).

