---
title: "8.3. Multiple regression with interaction: puromycin example"
author: "Lieven Clement & Jeroen Gilis"
date: "statOmics, Ghent University (https://statomics.github.io)"
output:
  html_document:
    code_download: yes
    theme: cosmo
    toc: yes
    toc_float: yes
    highlight: tango
    number_sections: yes
---

<a rel="license" href="https://creativecommons.org/licenses/by-nc-sa/4.0"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png" /></a>

```{r, message=FALSE, warning=FALSE}
library(tidyverse)
library(multcomp)
```

# Puromycin data

Gegevens over de snelheid van een enzymatische reactie werden verkregen door Treloar (1974).
Het aantal tellingen per minuut radioactief product van de reactie was
gemeten als functie van substraatconcentratie in delen per miljoen (ppm) en
uit deze tellingen werd de initiële snelheid (of snelheid) van de reactie berekend (tellingen/min/min). 
Het experiment werd één keer uitgevoerd met het behandelde enzym met Puromycin, en eenmaal met het enzym onbehandeld.

Ga na of er een verband is tussen de substraatconcentratie en de snelheid
**voor zowel de behandelde als de onbehandelde enzymen.**

# Importeer data

```{r}
data(Puromycin)
```

# Data verwerking

Voor een duidelijkere interpretatie van de modelparameters, maken we de onbehandelde enzymes de referentie categorie.

```{r}
Puromycin <- Puromycin %>%
    mutate(state =fct_relevel(state,c("untreated","treated")))
```

## Vraag1: Data Exploratie

Eerst visualizeren we de associatie tussen concentration en enzyme rate, voor beide enzyme states.

Maak een plot voor iedere enzyme state (met facet_wrap) waarin we alle datapunten kunnen zien. Fit een rode kromme lijn en zwarte rechte lijn door de punten. Sla deze figuur op in assocation_plot1

```{r}
assocation_plot1 <- ...  %>%
    ggplot(...) +
    ... + #Voeg de punten toe
    geom_smooth(...) + #Voeg een rode kromme lijn toe
    geom_smooth(...) + #Voeg een zwarte rechte lijn toe
    ylab("Reaction rate (counts/min)") +
    xlab("Substrate concentratie (ppm)") +
    facet_wrap(...)
```


In de plot kunnen we zien dat er een relatie is tussen rate en concentration, deze relatie ziet er wel niet lineair uit.

Om de impact van een log transformatie van concentratie na te gaan maken we een nieuwe kolom log10conc aan in Puromycin. Omdat de concentratie gemeten is in ppm zullen we een log10 transformatie gebruiken.

```{r}
Puromycin <- ...
```

Plot de vorige figuur nogmaals maar met de log10conc variabele. Sla deze figuur op in assocation_plot2. 

```{r}
assocation_plot2. <- ...  %>%
    ggplot(...) +
    ... + #Voeg de punten toe
    geom_smooth(...) + #Voeg een rode kromme lijn toe
    geom_smooth(...) + #Voeg een zwarte rechte lijn toe
  ylab("Reaction Rate (counts/min)") +
  xlab("log10(Substrate concentration) (log10 ppm)") +
  facet_wrap(...)
```

De relatie tussen rate en log10 getransformeerde concentratie ziet er lienair uit.

# Vraag2: Het lineair model

We zullen het volgende model fitten.

$$Y_i = \beta_0 + \beta_c x_c+ \beta_s x_s +\beta_{c:s}x_{c}x_{s} + \epsilon_i$$

met

- $$Y_i$$ de reactie snelheid (rate),

- $$\beta_0$$ het intercept,

- $$\beta_{c}$$ het hoofd effect voor de log10 concentratie,

- $$x_c$$ de log10 concentratie

- $$\beta_{s}$$ het hoofd effect voor de behandeling (treatment)

- $$x_s$$ een dummy variable voor "state" dat 0 is wanneer de enzymes niet behandeld zijn en 1 als ze behandeld zijn met Puromycin,

- $$\beta_{c:s}$$ Het interactie effect tussen concentratie en behandeling,

- $$\epsilon_i$$ i.i.d. normaal verdeeld met gemiddelde 0 en variantie $$\sigma^2$$.

Let op, we schrijven de substraat concentratie met een kleine letter omdat de predictor niet random is. 
De onderzoekers hebben de substraat concentraties gekozen tijdens het ontwerpen van het experiment en is dus geen random variable.


Het model impliceert twee verschillende regressie lijnen

- geen behandeling ($$x_s = 0$$) 
$$
Y_i = \beta_0 + \beta_c x_c + \epsilon  
$$
- behandeling ($$x_s = 1$$)

$$
Y_i = (\beta_0 + \beta_s)  + (\beta_c+\beta_{c:s}) x_c + \epsilon  
$$

Dus het hoofd effect voor behandeling heeft de interpretatie als de verandering in intercept tussen behandelde en niet behandelde samples.
De interactie term heeft de interpratie als de verandering in helling tussen behandelde en niet behandelde samples

Fit het correcte lineair model en sla dit op in mod1.
```{r}
mod1 <- ...
summary(mod1)
```

# Vraag 3: Assumpties

Voor we verder gaan, zullen we eerste de assumpties na gaan.

1. Onafhankelijkheid 
2. lineariteit
3. Normale verdeling van de residuals
4. Homoscedasticiteit

## Onafhankelijkheid

We nemen aan dat het experiment goed ontworpen was en dat de verschillende reacties die gebruikt werden in het experiment onafhankelijk zijn.


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

Welke stelling is correct

1. We zien in de residual plot dat er niet voldaan wordt aan de aanname van lineariteit
2. De QQ-plot toont grote afwijken van de normaliteit
3. We zien dat er een grote spreiding is in de residuals, we mogen dus homoscedasticiteit aannemen
4. Aan alle aannames is voldaan

# Vraag 4: Interpretatie
```{r}
summary(mod1)
```

Welke interpretatie van de schattingen is juist

1. De verwachte rate van een enzymatische reactie, onder behandeling van Puromycin, zal een gemiddelde rate hebben die 62.129 counts/min/min hoger is vergeleken met de rate van een enzymatische reactie onder behandeling van Puromycin, indien de substraat concentratie tien keer hoger is.

2. De verwachte rate van een enzymatische reactie, onder behandeling van Puromycin, zal een gemiddelde rate hebben die 164.588 counts/min/min hoger is vergeleken met de rate van een enzymatische reactie onder behandeling van Puromycin, indien de substraat concentratie tien keer hoger is.

3. De verwachte rate van een enzymatische reactie, onder behandeling van Puromycin, zal een gemiddelde rate hebben die 44.606 counts/min/min hoger is vergeleken met de rate van een enzymatische reactie onder behandeling van Puromycin, indien de substraat concentratie tien keer hoger is.

4. De verwachte rate van een enzymatische reactie, onder behandeling van Puromycin, zal een gemiddelde rate hebben die 85.45 counts/min/min hoger is vergeleken met de rate van een enzymatische reactie onder behandeling van Puromycin, indien de substraat concentratie tien keer hoger is.

## Inferentie 

# Vraag 5: Test of er een effect is van de log10 concentratie op de rate

Test of er een effect is van de log10concentratie. Maak hiervoor een gereduceerd model, sla dit model op in mod0. Vergelijk de modellen met elkaar via de anova-functie, sla dit op in anova1.

```{r}
mod0 <- ...
anova1 <- ...
```

# Vraag 6: Interactie test

Test of er een interactie-effect aanwezig is met de Anova functie. Sla het resultaat op in anova_interaction.

```{r}
anova_interaction <- car::Anova(...)
```

# Vraag 7: Hypotheses

We kunnende interactie niet verwijderen van het model. Hierdoor is het niet mogelijk om het effect van concentratie na te gaan zonder rekening te houden met de behandeling. We gaan de volgende onderzoeksvragen na.

1. De associatie tussen rate en concentratie is significant in de niet behandelde groep

$$
H_0: \beta_c = 0 \text{ vs } H_1: \beta_c \neq 0
$$

2. De associatie tussen rate en de concentratie is significant in de behandelde groep

$$

H_0: \beta_c + \beta_{c:s}= 0 \text { vs } H_1: \beta_c + \beta_{c:s}\neq 0
$$

3. De associatie tussen rate en concentratie is verschillend tussen de behandelde en niet behandelde groep

$$
H_0: \beta_{c:s}= 0\text{ vs }H_1: \beta_{c:s}\neq 0
$$


# Vraag 8: Onderzoek de hypotheses

We kunnen de hypotheses onderzoeken door gebruik te maken van de `glht` functie van de `multcomp` library. Sla het resultaat op in mcp1.

Bereken ook de betrouwbaarheidsintervallen, sla dit op in mcp1_confint.

Let op: schrijf de hypotheses uit in de volgorde van leesopdracht 7

```{r}
set.seed(0)
mcp1 <- glht(...)
mcp1_confint <- ...
```


### Conclusie

Er is een extreem signifcant effect van de substraat concentratie op de rate (p<<0.001)

Het effect van de substraat concentratie op de rate is extreem significant voor reacties van niet behandelde enzymes. Een reactie in een concentratie dat tien maal hoger is, zal een rate hebben dat gemiddeld `r round(confint(mcp1)$confint[1,1],1)` counts/min hoger is (95% CI [`r round(confint(mcp1)$confint[1,-1],1)`] counts/min) (p << 0.001). 

Het effect van de substraat concentratie op de rate is extreem significant voor reacties van behandelde enzymes (p << 0.001). Een reactie in een concentratie dat tien maal hoger is, zal een rate hebben dat gemiddeld `r round(confint(mcp1)$confint[2,1],1)` counts/min hoger is (95% CI [`r round(confint(mcp1)$confint[2,-1],1)`] counts/min).

Het effect van substraat concentratie op de rate is zeer significant hoger voor reacties van behandelde enzymes dan niet behandelde enzymes  (p = `r round(summary(mcp1)$test$pvalue[3],3)`). Een reactie in een concentratie dat tien maal hoger is zal gemiddeld een rate hebben dat `r round(confint(mcp1)$confint[3,1],1)` counts/min hoger is bij reacties van behandelde enzymes dan niet behandelde enzymes (95% CI [`r round(confint(mcp1)$confint[3,-1],1)`] counts/min).