---
title: "Algemeen lineair model: 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 (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)
#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))

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 is 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 %>% as.factor, y = gewicht)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter() +
  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 dosisgroepen.

- 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. Interacties modelleren komt later aan bod.

Om de data toch correct te modelleren, zullen we het dosiseffect in dit practicum apart modelleren per soort.

We starten hierbij met de dojovissen.  

# Schat het effect van de dosis na correctie voor het gewicht bij dojovissen.

## Aanpak van vorig practicum met enkelvoudige regressie

Model met enkel dosis effect:

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

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

## Aanpak waarbij we corrigeren voor het gewicht

We voegen nu een extra predictor (onafhankelijke variabele of covariaat) toe 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)$.

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

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

### Aannames

- We zien geen problemen m.b.t. lineariteit. De residuen zijn mooi verspreid rond 0 voor alle gefitte waarden. 

- We zien geen problemen mbt de aanname van gelijkheid van variantie. De residuen hebben een gelijke spreiding voor alle gefitte waarden. 

- De data lijken wat afwijkingen van normaliteit te vertonen in de linkerstaart.

- We gaan eerst na welke afwijkingen we kunnen verwachten in een dataset met evenveel residuen (`r sum(poison$soort == 0)`) als in onze dataset wanneer de residuen worden getrokken uit een normale verdeling met dezelfde standaard deviatie als wat we observeerden in onze dataset.

```{r}
set.seed(1406)
nobs <- sum(poison$soort == 0)

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

- We observeren dat de afwijkingen die we zien in de QQ-plot voor de residuen van het lineaire model voor dojovissen in de lijn ligt van de afwijkingen die we kunnen verwachten wanneer de residuen normaal verdeeld zouden zijn. De afwijkingen zijn dus niet extreem en bovendien beschikken we reeds over `r nobs` observaties zodat we er zeker van uit mogen gaan dat de parameterschatters approximatief normaal verdeeld zijn.

### Inferentie

We observeren dat het effect van de dosis veel significanter is na correctie voor gewicht.

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!

Het dosiseffect krijgt in dit model de betekenis van de gemiddelde verandering in log2-overlevingstijd tussen twee groepen dojovissen met hetzelfde gewicht die worden blootgesteld aan een dosis die 1 mg/l verschilt:

\begin{eqnarray}
\hat{\mu_1}&=& \beta_0 + \beta_d x_{1d} + \beta_g x_g \text{ (gemiddelde log2-overlevingstijd voor dosis 1)}\\
\hat{\mu_2}&=& \beta_0 + \beta_d x_{2d} + \beta_g x_g \text{ (gemiddelde log2-overlevingstijd voor dosis 2)}\\
\hat \mu_2- \hat \mu_1&=&\beta_0 + \beta_d x_{2d} + \beta_g x_g - (\beta_0 + \beta_d x_{1d} + \beta_g x_g) \text{ (verschil in gemiddelde log2-overlevingstijd tussen dosis 2 en dosis 1)}\\
\hat \mu_2-\hat \mu_1&=&\beta_d (x_{2d}-x_{1d})
\end{eqnarray}

Als het concentratieverschil $(x_{2d}-x_{1d})=1$ mg/l dan bekomen we inderdaad dat $\hat \mu_2-\hat \mu_1=\beta_d$. $\beta_d$ kwantificeert dus inderdaad het verschil in gemiddelde log2-overlevingstijd bij een dosisverschil van 1 mg/l voor dojovissen met hetzelfde gewicht.

Dus de log2-overlevingstijd bij dojovissen is gemiddeld `r format(abs(lmDojo$coef["dosis"]),digits=3)` lager bij vissen die blootgesteld worden aan een dosis van het gif 1mg/l hoger is, na correctie voor gewicht.

Als we data modelleren op log2-schaal bekomen we geometrische gemiddeldes, $\widehat{minsurv}$ na terugtransformatie:
\begin{eqnarray}
\log2 \widehat{minsurv}&=&\beta_0+  \beta_d x_{d} + \beta_g x_g\\
\widehat{minsurv}&=&2^{(\beta_0+  \beta_d x_{d} + \beta_g x_g)}\\
\widehat{minsurv}&=&2^{\beta_0}2^{\beta_d x_{d}}2^{\beta_g x_g}\\
\end{eqnarray}

Voor twee groepen dojovissen met hetzelfde gewicht die aan een verschillende dosis worden blootgesteld bekomen we dan:
\begin{eqnarray}
\frac{\widehat{minsurv}_2}{\widehat{minsurv}_1}&=& \frac{2^{\beta_dx_{2d}}}{2{\beta_dx_{1d}}}\\
&=& 2^{(\beta_dx_{2d} - \beta_dx_{1d})} \\
&=&2^{\left[\beta_d (x_{2d}-x_{1d})\right]}\\
\end{eqnarray}


```{r}
2^(lmDojo$coef)
2^(confint(lmDojo))
```


## 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)`)  als de dosis van het gif met 1 mg/l toeneemt (95% BI [`r paste(format(2^(confint(lmDojo)["dosis",]),digits=3),collapse=",")`]), na correctie voor het gewicht.

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



# 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 approximatief 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


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)`)  als de dosis van het gif met 1 mg/l toeneemt (95% BI [`r paste(format(2^(confint(lmDojo)["dosis",]),digits=3),collapse=",")`]), na correctie voor het gewicht.

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_d} = \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_g+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$).



