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

De snelheid van een enzymatische reactie werd gemeten door Treloar (1974). 
Hierbij werd gebruik gemaakt van een radioactief product en kon de reactie snelheid worden gemeten door gebruik te maken van het een geigerteller. De reactiesnelheid werd gemeten in aantal tellingen per minuut voor verschillende substraatconcentratie die werden geregistreerd in deeltjes per miljoen (parts per million, ppm). 
De data bestaat uit de initiële reactiesnelheid en de initiële substraatconcentratie voor elke reactie die werd uitgevoerd. 
De reacties werden uitgevoerd onder twee condities: met het onbehandeld enzym en met het enzym na behandeling met puromycine.

1. Ga na of er een verband is tussen de substraatconcentratie en de snelheid
**voor zowel de behandelde als de onbehandelde enzymen.** 
2. Ga na of het verband tussen de substraatconcentratie en de reactiesnelheid wijzigt na behandeling met pyromycine. 

# Importeer data

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

# Data verwerking

Voor een duidelijkere interpretatie van de modelparameters, definiëren we onbehandelde enzymes als referentieklasse.

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

# Vraag 1: Data Exploratie

Eerst visualizeren we de associatie tussen concentratie en de reactiesnelheid, voor beide enzyme statussen.

Maak een plot voor iedere enzyme status (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.
Wat neem je waar? 

```{r}
assocation_plot1 <- ...  %>%
    ggplot(...) +
    ... + #Voeg de punten toe
    geom_smooth(method = "loess", col = "red") + # fit een kromme door de punten (rode lijn)
    geom_smooth(method = "lm", col = "black") + # fit een rechte door de punten aan de hand van de kleinstekwadratenmethode
    ylab("Reaction rate (counts/min)") +
    xlab("Substrate concentratie (ppm)") +
    facet_wrap(...)
```

# Vraag 2: transformatie

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

Omdat het verband niet lineair is transformeren we de predictor substraatconcentratie. 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 omdat dit de interpretatie vergemakkelijkt.

```{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(...) + # fit een kromme door de punten (rode lijn)
    geom_smooth(...) + # fit een rechte door de punten aan de hand van de kleinstekwadratenmethode
    ylab("Reaction rate (counts/min)") +
    xlab("Substrate concentratie (ppm)") +
    facet_wrap(...)
```

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

# Vraag 3: Het lineair model

Het is mogelijk dat het verband tussen de substraatconcentratie en de reactiesnelheid verandert naar gelang de staat van het enzyme (onbehandeld of met puromycine behandeld). 
Er is dus mogelijks een interactie tussen de behandeling en de substraatconcentratie. 

We zullen daarom een model met hoofdeffecten voor behandeling en log10 -substraatconcentratie en een behandeling x log10-substraatconcentratie interactie. 

$$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 Puromyci

- $\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 regressiemodellen

- 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 hoofdeffect voor behandeling heeft de interpretatie als de verandering in intercept tussen behandelde en niet behandelde samples.
De interactie term heeft de interpretatie van 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 4: Assumpties

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

1. Onafhankelijkheid  (We nemen aan dat het experiment goed ontworpen was en dat de verschillende reacties die gebruikt werden in het experiment onafhankelijk zijn.)
2. lineariteit
3. Normale verdeling van de residuals
4. Homoscedasticiteit


```{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 5: Interpretatie parameters 
```{r}
summary(mod1)
```

Welke interpretatie van de schattingen is juist

1. De reactiesnelheid van een enzymatische reactie na behandeling van Puromycin is gemiddeld 62.1 counts/min hoger wanneer de substraat concentratie tien keer hoger is.

2. De reactiesnelheid van een enzymatische reactie na behandeling van Puromycin is gemiddeld 164.6  counts/min hoger wanneer de substraat concentratie tien keer hoger is.

3. De reactiesnelheid van een enzymatische reactie na behandeling van Puromycin is gemiddeld 44.6  counts/min hoger wanneer de substraat concentratie tien keer hoger is.


4. De reactiesnelheid van een enzymatische reactie na behandeling van Puromycin is gemiddeld 85.5  counts/min hoger wanneer de substraat concentratie tien keer hoger is.

# Vraag 6: Interpretatie parameters

Welke interpretatie van de schattingen is juist

1. De reactiesnelheid van een enzymatische reactie met het onbehandelde enzyme is gemiddeld 62.1 counts/min hoger wanneer de substraat concentratie tien keer hoger is.

2. De reactiesnelheid van een enzymatische reactie met het onbehandelde enzyme is gemiddeld 164.6  counts/min hoger wanneer de substraat concentratie tien keer hoger is.

3. De reactiesnelheid van een enzymatische reactie met het onbehandelde enzyme is gemiddeld 44.6  counts/min hoger wanneer de substraat concentratie tien keer hoger is.


4. De reactiesnelheid van een enzymatische reactie met het onbehandelde enzyme is gemiddeld 85.5  counts/min hoger wanneer de substraat concentratie tien keer hoger is.

# Vraag 7: Test of er een effect is van de log10 substraatconcentratie op de reactiesnelheid

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 8: Interactie test

Test of het effect van de substraatconcentratie op reactiesnelheid verschillend is bij reacties met behandeld en onbehandelde enzyme met de Anova functie. Sla het resultaat op in anova_interaction.

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

# Vraag 9: Hypotheses (leesopdracht)

We kunnen de 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. Is er een lineaire associatie tussen reactiesnelheid en de substraatconcentratie wanneer de reactie doorgaat met het onbehandelde enzyme?

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

2. Is er een lineaire associatie tussen reactiesnelheid en de substraatconcentratie wanner de reactie doorgaat met het enzyme dat werd behandeld met puromycine?

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

3. Is de lineaire associatie tussen reactiesnelheid en de substraatconcentratie verschillend tussen de reacties met het behandeld en onbehandeld enzyme? 

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


# Vraag 10: Test elk van 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 (leesopdracht)

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

Voor reacties met het onbehandelde enzyme is het effect van de initiële substraatconcentratie op de initiële reactiesnelheid extreem significant  (p << 0.001). Een reactie met het onbehandelde enzyme die doorgaat bij een initiële substraatconcentratie die 10 keer hoger ligt, heeft een  initiële reactiesnelheid die 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). 

Voor reacties met het enzyme dat met puromycine is behandeld, is het effect van de initiële substraatconcentratie op de initiële reactiesnelheid extreem significant  (p << 0.001). Een reactie met het behandelde enzyme die doorgaat bij een initiële substraatconcentratie die 10 keer hoger ligt, heeft een  initiële reactiesnelheid die 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 substraatconcentratie op reactiesnelheid is zeer significant hoger voor reacties met het behandelde enzyme dan voor reacties met het onbehandelde enzyme  (p = `r round(summary(mcp1)$test$pvalue[3],3)`). Een reactie bij een initiële substraatconcentratie die tien keer hoger is, zal gemiddeld een reactiesnelheid hebben die `r round(confint(mcp1)$confint[3,1],1)` counts/min hoger is bij reacties met het behandelde enzyme dan bij reacties met het onbehandelde enzyme (95% CI [`r round(confint(mcp1)$confint[3,-1],1)`] counts/min).