---
title: 'Oefening Correlatiecoëfficiënt'
output:
  html_document:
    code_download: yes
    highlight: tango
    number_sections: yes
    theme: cosmo
    toc: yes
    toc_float: yes
  word_document:
    toc: yes
---

# Colocalisatie aan de hand van correlatiecoefficient

De dataset in deze oefening is afgeleid uit 12 microscopie opnames die de colocalisatie van transferrine (gekleurd in texasrood) en Rab10 proteïnen (gekleurd in groen m.b.v. GFP) in kaart wilden brengen volgens het werk van McDonald & Dunn (2013).
Om de colocalisatie van beide eiwitten te kwantificeren wordt een biologisch staal aan de hand van probes (vaak antilichamen) gekleurd met beide stoffen.
De bepaling van hun colocalisatie werd vroeger typisch visueel uitgevoerd: indien de meeste pixels geel (=groen+rood) zijn op een donkere achtergrond, dan spreken we van een sterke colocalisatie.
Dit is echter subjectief en laat geen mate van sterkte in de colocalisatie toe.
Een objectievere manier is om de intensiteit van het rode en het groene signaal te meten in iedere pixel, en diens Pearson correlatiecoëfficiënt te berekenen voor het beeld.
Voor elk biologisch staal bekomt men aldus een microscopiebeeld, dat men zal samenvatten in diens correlatiecoëfficiënt.
Om colocalisatie te testen, kan men aldus een statistische test uitvoeren of de correlatiecoëfficiënt al dan niet groter is dan nul.
Naast een objectieve analyse, laat deze methode ook toe om de mate van colocalisatie te kwantificeren aan de hand van de gemiddelde correlatiecoëfficiënt.

# R libraries inlezen

```{r}
library(ggplot2)
library(dplyr)
#install.packages("tidyr")
library(tidyr)
library(readr)
```

# Dataset correlation.dat inlezen

Het inlezen van de dataset kan via onderstaand commando met behulp van de weblink waar de data is opgeslagen:

```{r}
correlatie <- read_table("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum2/correlation.dat")
correlatie
```

# Vraag 1: Stel de nulhypothese en alternatieve hypothese op.

Vertaal de vraagstelling naar een nulhypothese $H_0$ en een alternatieve hypothese $H_A$. 
 
$H_0$: De correlatiecoëfficiënt $\rho$ is ...

$H_A$: De correlatiecoëfficiënt $\rho$ is ...


# Vraag 2: Maak een histogram van de correlaties.

Gebruik een ggplot visualisatie om een histogram te maken (met bins = 10) van de correlaties. 

```{r}
gg_hist <- ... %>% ggplot() + # specifieer de dataset en de aesthetics
...   # om een histogram te maken, gebruik bins = 10
gg_hist
```

# Vraag 3: Bereken de  gemiddelde correlatiecoëfficiënt en diens standard error.

Sla je resultaten respectivelijk op in `gemiddelde` en `se`.


```{r}
gemiddelde <- ...
gemiddelde
```


```{r}
stdev <- ... # standaarddeviatie correlatiecoëfficiënt
se <- ... 
se
```

# Vraag 4: Ga na of de data afkomstig kunnen zijn uit een normale verdeling.

Gebruik hiervoor een ggplot visualisatie en sla deze op in `plot`.

```{r}
plot <- correlatie %>% ... 
plot
```


**Merk op** dat we slechts twaalf datapunten hebben. Dit zijn weinig datapunten om de assumptie van normaliteit na te gaan, maar indien we geen sterke afwijkingen observeren, kunnen we wel de assumptie aannemen dat de data bij benadering een normale verdeling volgen. 

# Vraag 5: Afwijkingen van normaliteit of niet?

Het is niet altijd eenvoudig om op basis van een klein aantal datapunten de normaliteit te beoordelen. Op de QQ-plot in de vorige opgave zouden we een lichte afwijking in de linkse staart kunnen opmerken. De vraag is dan of de geobserveerde afwijking hier valt onder de te verwachten afwijking ven een normale distributie. We kunnen dit nagaan door QQ-plots te simuleren van 12 datapunten onder een normale verdeling.

```{r}
set.seed(1)
n_simulations <- 9
n_datapoints <- 12

SimulationDataFrame <- sapply(1:n_simulations, 
                              function(x){
                                rnorm(n=n_datapoints, mean=gemiddelde, sd=stdev)
                                }) %>%
                              as_tibble() %>% 
                              pivot_longer(cols = 1:n_simulations, 
                                           names_to = "SimulationIndex",
                                           values_to = "DataPoint") 

SimulationDataFrame %>% 
  ggplot(aes(sample = DataPoint)) + 
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~SimulationIndex,
             scales = "fixed") + 
  ylab("Sample Quantiles") +
  xlab("Theoretical Quantiles") +
  ggtitle("Simulation QQ Plots") + 
  theme_bw()
```

We kunnen vaststellen dat de bekomen QQ-plot niet afwijkt van de 9 plots die we hier bekomen voor data die uit een normale verdeling gesimuleerd werden. We hebben dus geen aanwijzing dat de 12 correlatiecoëfficiënten niet afkomstig zouden zijn uit een normale verdeling.

# Vraag 6: Voer een t-test uit om na te gaan of je de nulhypothese kan verwerpen op het 5% significantieniveau.

We werken hier met een eenzijdige t-test. Aangezien de wetenschappers enkel geïnteresseerd zijn of de correlatiecoëfficiënt groter is dan nul, zullen we een eenzijde t-test moeten uitvoeren. Bij `?t.test` zien wat dat standaard een tweezijdige test zal uitgevoerd worden. Door dit argument aan te passen naar "greater" of "less" zal een eenzijdige test uitgevoerd worden. Bepaal zelf welke van de twee opties hier van toepassing is. Sla de t-test op in `t_test`.

```{r}
t_test <- t.test(..., alternative = "...")
t_test
```

# Vraag 7: De p-waarde van de test opslaan

De p-waarde van de test uit de vorige oefening kunnen we extraheren uit de output van `t.test`. We slaan deze waarde op in `p_waarde` op de volgende manier: 

```{r}
p_waarde <- t_test$p.value
p_waarde
```


# Vraag 8: Zal je de nulhypothese kunnen verwerpen op het 5% significantieniveau?

Bij vraag 7 stelden we vast dat de p-waarde van deze test gelijk is aan 5.63e-11. Zullen we de nulhyposthese kunnen verwerpen op het 5% significantieniveau?

- Ja, want de p-waarde is kleiner dan 0.05.
- Nee, want de p-waarde is kleiner dan 0.05.


# Vraag 9

Bereken het 95% betrouwbaarheidsinterval met behulp van de`t.test` functie. Het 95% betrouwbaarheidsinterval wordt meegegeven met de output van deze functie op de volgende manier:

```{r}
t_test$conf.int[1:2]
```


**Opmerking** Het betrouwbaarheidsniveau van het interval dat wordt meegegeven met de `t.test` functie kan aangepast worden met het argument `conf.level`. Zie hiervoor ook `?t.test`.

# Vraag 10

Bereken de teststatistiek, het 95% betrouwbaarheidsinterval en de corresponderende p-waarde zonder gebruik te maken van de functie `t.test`. Sla de waarden respectievelijk op in `teststat`, `conf_95` en `pval`. Voor extra informatie rond het bepalen van betrouwbaarheidsintervallen verwijzen we naar de theorie uit de Dodona cursus statistiek "5.4 Intervalschatters (Vervolg 1)" en de bijhorende kennisclip. Let wel; Let wel; in 5.4 worden enkel tweezijdige betrouwbaarheidsintervallen besproken.


```{r}
alpha <- 0.05 #significantieniveau
n <- 12 #steekproefgrootte

#teststatistiek T
teststat <- gemiddelde/se #gemiddelde en se werden reeds berekend hierboven.
teststat

#95% betrouwbaarheidsinterval
t <- qt(1-alpha, n-1) # 95% percentiel uit t-verdeling met 11

ondergrens <- gemiddelde-t*se # ondergrens BI
bovengrens <- Inf # bovengrens BI
conf_95 <- c(ondergrens, Inf)
conf_95

#p-waarde
pval <- 1-pt(teststat, df=n-1)
pval
```


# Vraag 11 Interpreteer het 95\% betrouwbaarheidsinterval en formuleer een conclusie gebaseerd op bovenstaande analyses. (leesopdracht)

**Interpretatie van het 95\% betrouwbaarheidsinterval**

In onze analyse verkrijgen we een 95% betrouwbaarheidsinterval van [0.53, 1]. Onder het herhaaldelijk uitvoeren van het experiment, zal de werkelijke gemiddelde correlatiecoëfficiënt in 95% van de berekende 95% betrouwbaarheidsintervallen zitten.

Merk op dat we het 95\% betrouwbaarheidsinterval geven van 0.53 tot 1 en niet tot oneindig. Een correlatiecoëfficiënt kan immers nooit groter zijn dan 1. 
Merk ook op dat er steeds een één-op-één correspondentie is tussen het betrouwbaarheidsinterval en de p-waarde. Wanneer p < 0.05, dan is 0 niet ingesloten is in het 95\% betrouwbaarheidsinterval en omgekeerd.

Conclusie:

De gemiddelde correlatiecoëfficiënt tussen transferrine en Rab10 is gelijk aan 0.58 (95\% BI: 0.53 tot 1). 
Deze gemiddelde correlatiecoëfficiënt is sterk significant groter dan nul op het 5\% significantieniveau (p < 0.0001). 
Dit suggereert dat er in zekere mate colocalisatie tussen de proteïnen plaats vindt.
