---
title: "Algemeen Lineair Model: niet-additief meervoudig lineair regressie model"
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 (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.


## Laad de libraries

```{r}
library(dplyr)
library(ggplot2)
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))
```

## Vraag 1
Geef een scatterplot voor de log2-overlevingstijd ten opzichte van de dosis, voor elke soort apart. Trek ook een best passende rechte door elke figuur.

```{r}
scatterplot1<-poison%>% ggplot(...)+
  ...+
  stat_smooth(....)+
  facet_wrap(...)
scatterplot1
```

## Vraag 2 (Meerkeuzevraag)
Beschouw de volgende stellingen over bovenstaande figuur

A. De log2-overlevingstijden liggen gemiddeld ongeveer even hoog voor elke soort.

B. Het effect van dosis op de log2-overlevingstijden is voor elke soort ongeveer hetzelfde

Welke stelling(en) is/zijn correct?

1. Zowel A als B
2. Enkel A
3. Enkel B
4. A noch B


## Vraag 3
Geef een scatterplot voor de log2-overlevingstijd ten opzichte van het gewicht, voor elke soort apart. Trek ook een best passende rechte door elke figuur.

```{r}
scatterplot2<-...
scatterplot2
```


## Vraag 4 (Meerkeuzevraag)
Beschouw de volgende stellingen over de figuur:

A. Het bereik van het gewicht is voor alle soorten quasi gelijk.

B. Het effect van gewicht op de log2-overlevingstijden is voor elke soort ongeveer hetzelfde

Welke stelling(en) is/zijn correct?

1. Zowel A als B
2. Enkel A
3. Enkel B
4. A noch B

# Het lineaire model met interacties soort-dosis en soort-gewicht.
## Vraag 5 (Leesopdracht)
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.
We zullen ook toelaten dat elke soort een verschillend effect van gewicht kan hebben.

Het model zal dus het volgende zijn:
$$
y_i=\beta_0+\beta_d x_{id} + \beta_g x_{ig} +\beta_{sg} x_{isg} +\beta_{sz} x_{isz} + \beta_{d:sg} x_{id}x_{isg} + \beta_{d:sz} x_{id}x_{isz} + \beta_{g:sg} x_{ig}x_{isg} + \beta_{g:sz} x_{ig}x_{isz} + \epsilon_i,
$$

met $y_i$ het gewicht, $x_{id}$ de dosis en $x_{ig}$ het gewicht van vis $i$. $x_{isg}$ is een dummy-variabele die aangeeft of vis $i$ een goudvis is en $x_{isz}$ een dummy-variabele die aangeeft of de vis een zebravis is.
De referentieklasse is dus voor de soort dojovissen (als  $x_{isg}=0$ en $x_{isz}=0$).
Verder is $\epsilon_i \text{ i.i.d. } N(0,\sigma^2)$.


Dit model zal eigenlijk drie regressievlakken oplevert: 1 voor dojovissen, 1 voor goudvissen en 1 voor zebravissen.

$$
\begin{array}{ll}
\text{dojovis } (x_{isg}=0\text{ en }x_{isz}=0):&E[y\vert \text{dojovis}]=\beta_0+\beta_d x_{id} + \beta_g x_{ig}\\
\text{goudvis }  (x_{isg}=1\text{ en }x_{isz}=0):&E[y\vert \text{goudvis}]=\beta_0+\beta_{sg}+(\beta_d+ \beta_{d:sg}) x_{id} + (\beta_g+\beta_{g:sg}) x_{ig}\\
\text{zebravis }  (x_{isg}=0\text{ en }x_{isz}=1):&E[y\vert \text{zebravis}]=\beta_0+\beta_{sz}+(\beta_d+ \beta_{d:sz}) x_{id} + (\beta_g+\beta_{g:sg}) x_{ig}\\
\end{array}
$$



Het hoofdeffect voor soort zorgt dus dat elke soort een verschillend intercept heeft ($\beta_0$, $\beta_0 + \beta_{sg}$, $\beta_0+\beta_{sz}$ voor dojo-, goud- en zebravissen, respectievelijk).
De interactie tussen soort en dosis zorgt ervoor dat het dosiseffect (helling) verschillend kan zijn voor elke soort ($\beta_d$, $\beta_d+\beta_{d:sg}$, $\beta_d+\beta_{d:sz}$ voor dojo-, goud- en zebravissen, respectievelijk).
De interactie tussen soort en dosis zorgt ervoor dat het gewichtseffect (helling) verschillend kan zijn voor elke soort ($\beta_g$, $\beta_g+\beta_{g:sg}$, $\beta_g+\beta_{g:sz}$ voor dojo-, goud- en zebravissen, respectievelijk).
Het model heeft verder als voordeel dat het alle data gebruikt om de variantie te schatten van de residuen.

# Voorwaarden van het lineaire model

## Vraag 6 (Meerkeuzevraag)
Ga de voorwaarden van het lineaire model na
```{r}
lmInt <- lm(log2minsurv~ soort + dosis + gewicht + soort:dosis + soort:gewicht ,data = poison)
plot(lmInt)
```

Welk van de voorwaarden van het lineaire regressiemodel is mogelijks niet perfect voldaan, op basis van deze figuren?

1. Onafhankelijkheid 
2. Lineariteit
3. Normaliteit
4. Homoscedasticiteit

## Vraag 7 (Leesopdracht)

De QQ-plot suggereert mogelijks lichte afwijkingen van normaliteit.
Er zijn echter veel observaties waardoor we kunnen aannemen dat de gemiddelden approximatief normaal verdeeld zullen zijn. Bovendien zijn de afwijkingen symmetrisch.
Voor symmetrische distributies convergeert de verdeling van de parameterschatters sneller naar de normale verdeling en kunnen we de centrale limietstelling sneller toepassen.

Verder tonen we via simulatie aan dat de afwijkingen in de QQ-plot vallen binnen hetgene we kunnen verwachten o.b.v. gegevens uit een normale verdeling. (Merk op dat orig slaat op de QQ-plot van de residuen voor het model met interacties)

```{r}
set.seed(1025)
nobs <- nrow(poison)

data.frame(
  y = c(lmInt$res,
        rnorm(nobs*8,
              sd = sigma(lmInt)
             )
      ),
  label = rep(
              c("orig",
                paste0("sim",1:8)),
              each = nobs)) %>%
  ggplot(aes(sample = y)) +
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~ label)
```

# Interpretatie parameterschattingen

## Vraag 8 (Meerkeuzevraag)
Bekijk de summary van ons lineaire model

```{r}
summary(lmInt)
```
Beschouw de volgende interpretatie van de schatting van het intercept:

"De gemiddelde log2-overlevingstijd bij dojovissen met een gewicht van 0 gram die geen gif toegediend kregen is 0.8294."

Wat kan je zeggen over deze stelling?

1. De interpretatie is correct.
2. De interpretatie is fout, aangezien de p-waarde bij intercept niet significant is. 
3. De interpretatie is fout, aangezien het hier niet om dojovissen gaat, maar om alle vissen in het algemeen.
4. De interpretatie is fout, aangezien het berust is op extrapolatie en biologisch geen steek houdt.

## Vraag 9 (Meerkeuzevraag)
Beschouw de volgende interpretaties over het effect van dosis. Welke interpretatie is volledig correct?


1. De overlevingstijd van *vissen* is gemiddeld 0.959 eenheden lager op log2 schaal als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is, *na correctie voor gewicht*.
2. De overlevingstijd van *dojovissen* is gemiddeld 0.959 eenheden lager op log2 schaal als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is, *na correctie voor gewicht*.
3. De overlevingstijd van *vissen* is gemiddeld 0.959 eenheden lager op log2 schaal als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is.
4. De overlevingstijd van *dojovissen* is gemiddeld 0.959 eenheden lager op log2 schaal als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is.
  

## Vraag 10 (Meerkeuzevraag)
Beschouw de volgende interpretaties over het effect van gewicht. Welke interpretatie is volledig correct?

1. De log2-overlevingstijd bij goudvissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is gemiddeld 0.48 eenheden hoger per gram bij de zwaardere vis.
2. De log2-overlevingstijd bij goudvissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is gemiddeld 0.59 eenheden lager per gram bij de zwaardere vis.
3. De log2-overlevingstijd bij goudvissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is gemiddeld 1.08 eenheden hoger per gram bij de zwaardere vis.
4. De log2-overlevingstijd bij goudvissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is gemiddeld 1.56 eenheden hoger per gram bij de zwaardere vis.

## Vraag 11 (Meerkeuzevraag)
Welke interpretatie op **originele schaal** is correct.

1. Het geometrisch gemiddelde van de overlevingstijd bij zebravissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is een factor 0.074 lager per gram bij de zwaardere vis.
2. Het geometrisch gemiddelde van de overlevingstijd bij zebravissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is een factor 1.052 hoger per gram bij de zwaardere vis.
3. Het geometrisch gemiddelde van de overlevingstijd bij zebravissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is een factor 2.007 hoger per gram bij de zwaardere vis.
4. Het geometrisch gemiddelde van de overlevingstijd bij zebravissen die verschillen in gewicht en die aan eenzelfde dosis werden blootgesteld is een factor 2.112 hoger per gram bij de zwaardere vis.

## Vraag 12 (Meerkeuzevraag)

Bekijk volgende interpretaties. Welke interpretatie(s) is/zijn volledig correct?

A.  Het gemiddeld verschil in log2-overlevingstijd tussen goudvissen en dojovissen die beiden een gewicht hebben van 0 gram en zich in een aquarium zonder gif bevinden is -0.59 . 
B.  Het gemiddeld verschil in log2-overlevingstijd tussen goudvissen en dojovissen die eenzelfde gewicht hebben en evenveel gif toegediend kregen is -0.59. 


1. Zowel A als B
2. Enkel A
3. Enkel B
4. A noch B

## Vraag 13 (Meerkeuzevraag)
Wat kan je zeggen over volgende interpretatie.

De log2-overlevingstijd tussen goudvissen neemt gemiddeld met 0.169 eenheden minder snel af per mg/l toename van de dosis waaraan ze worden blootgesteld dan bij dojovissen, na correctie voor gewicht. Het effect van de dosis is dus minder groot bij goudvissen dan bij dojovissen.

1. De interpretatie is correct, maar het verschil in effect is niet significant
2. De stelling is theoretisch gezien correct, maar biologisch niet relevant.
3. De stelling is niet correct.
4. Geen van bovenstaande antwoorden is juist.

# Effect van dosis
## Vraag 14 (Leesopdracht)

De onderzoekers zijn vooral geïnteresseerd in het effect van de dosis.
Aan de hand van de summary output, zouden we testen kunnen doen om het dosis effect te evalueren in dojo-vissen.
En of het effect van de dosis verschilt tussen goud(zebra) vissen en dojo vissen , de interactie term.
Het dosiseffect houdt dus de evaluatie in van meerdere model parameters.

We voeren eerst een test uit om na te gaan of het effect van de dosis significant is.
Onder de nulhypothese veronderstellen we dat er geen effect is van de dosis (bij geen enkele vissoort):

\[H_0: \beta_{d} = \beta_{d} + \beta_{d:sg} = \beta_{d} + \beta_{d:sz}=0\]
we zullen dit testen tegen de alternatieve hypothese dat er een effect is in minstens 1 vissoort:
\[H_1: \beta_{d} \neq 0 \text{ en/of } \beta_{d} + \beta_{d:sg} \neq 0 \text{ en/of } \beta_{d} + \beta_{d:sz}\neq{0}\]

waarbij $\beta_{d}$, $\beta_{d} + \beta_{d:sg}$ en $\beta_{d}+ \beta_{d:sz}$ het dosis effect is voor dojo-, goud- en zebravissen, respectievelijk.

De nullhypothese impliceert dus eveneens dat er geen hoofdeffect en geen soort x dosis interactie:
$H_0: \beta_{d} = \beta_{d:sg} =\beta_{d:sz}=0$

## Vraag 15

We kunnen deze test uitvoeren door middel van een F-test waarin we het model met en het model zonder hoofdeffect voor dosis en zonder soort x dosis interactie met elkaar vergelijken.

```{r}
lmNoDose<- lm(log2minsurv~ soort + gewicht + soort:gewicht ,data=poison)
anova_test <- anova(lmNoDose,lmInt)
anova_test
```

## Vraag 16
Wat kunnen we besluiten?

1. We zien dat het dosis effect niet significant is ($p<<0.001$).
2. We zien dat het dosis effect suggestief is ($p<<0.001$).
3. We zien dat het dosis effect sterk significant is ($p<<0.001$).
4. We zien dat het dosis effect extreem significant is ($p<<0.001$).

# Interactie-effect
## Vraag 17 
We testen nu of er een significante interactie is. Hiervoor moeten we opnieuw meerdere parameters evalueren:
$H_0: \beta_{d:sg} = \beta_{d:sz} =0$ vs $H_1: \beta_{d:sg} \neq 0 \text{ en/of } \beta_{d:sz}\neq{0}$

We kunnen dat opnieuw doen a.d.h.v. een F-test tussen het volledige model en een model waar de interactie tussen soort en dosis is verwijderd maar waarin het hoofdeffect voor dosis wordt behouden.

```{r}
lmNoInt<- ...
anova_tabel1<-anova(...)
anova_tabel1
```
## Vraag 18

Merk op dat we deze analyse ook konden doen a.d.h.v. een anova tabel met type III kwadratensommen.
Deze kwadratensommen kunnen niet verkregen worden via de standaard anova functie.
De `Anova` functie uit het car package laat dat wel toe.

```{r}
anovalmInt<-...
anovalmInt
```

## Vraag 19
Wat kunnen we besluiten omtrent de interactie dosis x soort?

1. We zien dat de dosis x soort interactie niet significant is ($p=0.69$).
2. We zien dat de dosis x soort interactie suggestief is ($p=0.69$).
3. We zien dat de dosis x soort interactie significant is ($p=0.69$).
4. We zien dat de dosis x soort interactie extreem significant is ($p=0.69$).

## Vraag 20
Wat kunnen we besluiten omtrent de interactie gewicht x soort?

1. We zien dat de gewicht x soort interactie niet significant is ($p=0.028$).
2. We zien dat de gewicht x soort interactie suggestief is ($p=0.028$).
3. We zien dat de gewicht x soort interactie significant is ($p=0.028$).
4. We zien dat de gewicht x soort interactie sterk significant is ($p=0.028$).

# Interactie-effect: Conclusie op basis van vereenvoudigd model

## Vraag 21 (Leesopdracht)

Verder includeren we de gewichtseffecten vooral in het model om het dosiseffect zuiver te kunnen schatten en om de precisie van de standard error op het dosis effect te vergroten.

De conventionele aanpak is om nu de niet significante interactie term uit het model te verwijderen en de analyse te doen a.d.h.v. het model met enkel een hoofdeffect voor de dosis.

```{r}
summary(lmNoInt)
```

## Vraag 22 (Meerkeuzevraag)
Welke conclusie is juist op basis van het vereenvoudigd model

1. De log2-overlevingstijd tussen vissen is gemiddeld 0.8416 eenheden lager als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is.
2. De log2-overlevingstijd tussen vissen is gemiddeld 0.8416 eenheden lager als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is, bij gelijk gewicht.
3. De log2-overlevingstijd tussen vissen is gemiddeld 0.8416 eenheden lager als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is, bij gelijk gewicht en soort.
4. De log2-overlevingstijd tussen vissen is gemiddeld 0.8416 eenheden lager als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is, bij dojovissen met gewicht nul.


## Vraag 23
Bereken de schattingen van de effecten op originele schaal en de bijhorende betrouwbaarheidsintervallen. 
```{r}
schattingen_log2<-coef(lmNoInt)
schattingen_orig<-...
schattingen_orig

intervallen_log2_int<-confint(lmNoInt)
intervallen_orig_int<-...
intervallen_orig_int
```

## Vraag 24
Welke van volgende uitspraken is geen conclusie van onze analyses?

1. Er is een extreem significant effect van de dosis op de overleving bij vissen ($p<<0.001$).
Het dosiseffect is niet significant verschillend tussen de soorten ($p=0.69$).
2. Het geometrisch gemiddelde van de overlevingstijd halveert ongeveer (factor $0.558$) bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.472,0.660]), na correctie voor gewicht.
3. Merk op dat we in de conclusie niets over de significantie van het gewichtseffect hebben vermeld gezien de onderzoeksvraag zich focust op het effect van de dosis.
4. Alle bovenstaande uitspraken zijn correcte conclusies van onze analyses.

# Interactie-effect: Alternatieve aanpak zonder reductie van het model

## Vraag 25 (Leesopdracht)

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 effect bij de vissen in het experiment d.m.v. te marginaliseren over alle vissen in het experiment.


$$\frac{\sum\limits_{i=1}^n (\beta_d + \beta_{d:sg} X_{isg} +  \beta_{d:sz} X_{isz})}{n} = \beta_d +\beta_{d:sg} \frac{\sum\limits_{i=1}^n X_{isg}}{n} + \beta_{d:sz} \frac{\sum\limits_{i=1}^n X_{isz}}{n} = \beta_d +\beta_{d:sg} n_{sg}/n + \beta_{d:sg} n_{sz}/n$$

## Vraag 26
Tel het aantal vissen van elke soort met behulp van de count functie

```{r}
nS <- ...
nS
```
## Vraag 27 (Meerkeuzevraag)
Interpreteer de schatting op log2-schaal.
We evalueren nu het contrast voor het marginale effect:

```{r}
marginaleDosisEffect <- glht(lmInt, linfct = c("dosis+38/96*soort1:dosis + 19/96*soort2:dosis = 0"))
summary(marginaleDosisEffect)
```
We concluderen dat het marginaal gemiddelde van de log2-overlevingstijd van vissen 0.8419 eenheden lager is als de dosis van het gif waaraan ze werden blootgesteld 1 mg/l hoger is, bij gelijk gewicht en soort.

Wat is het voornaamste verschil tussen deze conclusie en de de conclusie die we trokken bij het vereenvoudigde model (zonder interactie tussen soort en dosis)?

1. De conclusies zijn identiek
2. De schatting van het effect is sterk verschillend
3. De standaard error van het effect is sterk verschillend
4. We spreken nu over het marginaal gemiddelde, omdat we marginaliseren over alle vissen.*


## Vraag 28
Geef het 95% betrouwbaarheidsinterval voor het marginale effect van dosis op log2-schaal.
```{r}
confint1<-...
confint1
```

## Vraag 29
Geef het marginale dosiseffect op originele schaal met bijhorend betrouwbaarheidsinterval
```{r}
interval_log2<-confint1$confint
interval_log2
interval_orig<-...
interval_orig
``` 

## Vraag 30 (Leesopdracht)

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

Het dosiseffect is niet significant verschillend tussen de soorten (p=0.69).

Het marginaal geometrisch gemiddelde van de overlevingstijd halveert ongeveer (het wordt vermenigvuldigd met een factor 0.558) bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.471, 0.661]), na correctie voor gewicht en soort. ($p<<0.001$)

# Meer uitgebreide rapportering van gemiddelden 
## Vraag 31 (Leesopdracht)

Veronderstel dat de onderzoekers over het gemiddelde dosiseffect voor elke vissoort alsook over het marginale dosiseffect wensten te rapporteren met 95% betrouwbaarheidsintervallen.

De effecten van interesse zijn opnieuw lineaire combinaties van modelparameters, ook lineaire contrasten genoemd. Het dosis effect voor de verschillende vissoorten is:
$$
\begin{array}{ll}
\text{dojovis:}&\beta_d\\
\text{goudvis:}&\beta_d+ \beta_{d:sg}\\
\text{zebravis:}&\beta_d+ \beta_{d:sz}\\
\end{array}
$$

We kunnen dan ook alle contrasten van interesse evalueren en hierbij corrigeren voor multiple testing.

```{r}
lmIntGlht <- glht(lmInt,
	linfct = c("dosis = 0",
		"dosis+soort1:dosis = 0",
		"dosis+soort2:dosis = 0",
		"dosis+38/96*soort1:dosis + 19/96*soort2:dosis = 0")
		)
confint(lmIntGlht)
```


Door middel van contrasten kunnen we dus over meerdere gemiddelden rapporteren en kunnen we hierbij ook corrigeren voor multiple testing.
