---
title: "Uitgewerkt voorbeeld niet-parametrische statistiek: Garnalen"
output:
    html_document:
      code_download: true    
      theme: cosmo
      toc: true
      toc_float: true
      highlight: tango
      number_sections: true
---


# Introductie

PCBs (Polychlorinated biphenyls, o.a. gebruikt in koelelementen en smeermiddel) stapelen zich makkelijk op in het vetweefsel van garnalen. In dit experiment werden twee groepen van elk 18 samples van garnalen gekweekt in verschillende condities. De garnalen werden gerandomiseerd over de twee groepen, en de PCB-concentratie werd gemeten door het vetweefsel te extraheren van 100g garnalen in een sample.
De onderzoeksvraag is om na te gaan of de PCB-concentratie in garnalen afhankelijk is van de groeicondities. De PCB-concentraties werden gemeten in pg/g vet.

## Libraries laden 
```{r}
library(ggplot2)
library(dplyr)
library(coin)
library(readr)
```


## Data inlezen

```{r}
shrimps <- read_delim("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum6/shrimps.txt",
                      delim=" ",
                      skip=1,
                      col_names = c("row","PCB.conc","group"),
                      col_types = list(col_integer(), col_double(), col_factor()),
                      col_select =c("PCB.conc","group")
                      )
```

## Data-exploratie

```{r}
boxplot1<-shrimps %>% ggplot(aes(x=group, y=PCB.conc, col=group))+
  geom_boxplot(outlier.shape=NA)+
  geom_jitter()+
  theme_bw()
boxplot1
```

De "outliers" in de data werden dubbel gecheckt door de onderzoekers en zijn correcte metingen.

```{r}
qqplot1<-shrimps%>%ggplot(aes(sample=PCB.conc))+
  geom_qq()+
  geom_qq_line()+
  facet_wrap(~group)
qqplot1
```

```{r}
hist1<-shrimps%>%ggplot(aes(x=PCB.conc, fill=group))+
  geom_histogram(alpha=0.2, position="identity", bins=6)
hist1
```

Conclusies data-exploratie:

 - De data lijken niet afkomstig te zijn uit een Normale verdeling;
 - Er zijn onvoldoende observaties om te berusten op de centrale limietstelling;
 - De varianties in beide groepen zijn mogelijks niet gelijk;
 - Er zijn grote outliers.
 
Geen enkel van de parametrische toetsen die we reeds gezien hebben kunnen gebruikt worden voor deze data analyse. We zullen dus beroep doen op niet-parametrische methoden.

# Wilcoxon rank-sum test/ Mann–Whitney U test


In deze sectie zullen we dus een alternatieve niet-parametrische toets uitvoeren, die de Wilcoxon rank-sum test of de Mann–Whitney U test wordt genoemd.
Beide teams van onderzoekers ontwikkelden een equivalente test in parallel. 

Indien we een locatie-shift model veronderstellen, of met andere woorden, veronderstellen dat de vorm van beide distributies gelijk is en dat deze enkel verschillen in locatie, dan kunnen we het resultaat van de rank test interpreteren in functie van een verschil in mediaan. Maar, de Wilcoxon-Mann-Withney test is ook geldig indien het locatie-shift model niet opgaat. In dat geval, kunnen we het resultaat interpreteren in functie van de probabilistische index $P(Y_{1} \ge Y_2)$.

Bij rank tests wordt de data $Y_i$ getransformeerd naar ranks
$$ R_i=R(Y_i) = \#\{Y_j: Y_j\leq Y_i; j=1,\ldots, n\}. $$
Ranks zijn erg robuust voor outliers: of de hoogste observatie nu een waarde van 100 of 10 heeft, het blijft dezelfde rank behouden.
In het geval van *ties*, dit zijn gelijke waarden, worden midranks gebruikt:
$$ R(Y_i) = \frac{\sum\limits_{\forall j : Y_j=Y_i}R(Y_j)}{\#{\forall j:Y_j=Y_i}}  $$
De midrank is het gemiddelde van de ranks van de gelijke observaties.

## Voorwaarden
De Wilcoxon-Mann-Whitney test heeft bijna geen voorwaarden. De enige veronderstelling die we nog steeds maken is dat de gegevens onafhankelijk zijn.

## Hypothesen
Aangezien we geen gemiddelden van de data gebruiken in de test, kunnen we de test uiteraard ook niet interpreteren in termen van deze gemiddelden. Via de ranks kunnen we het resultaat echter wel interpreteren aan de hand van de *probabilistische index*, i.e. de kans dat een willekeurige observatie uit de eerste groep groter is dan of gelijk aan een willekeurige observatie in de tweede groep.
De nulhypothese stelt dat de de verdelingen in beide groepen gelijk zijn.

$$ H_0: f_1 = f_2 $$

In woorden: de verdeling van PCB-concentraties in garnalen zijn in beide condities gelijk.

De alternatieve hypothese in termen van de probabilistische index kan men schrijven als volgt:

$$ H_1: P(Y_{i1} \geq Y_{i2}) \ne 1/2 $$

In woorden: de kans dat een willekeurige observatie van de PCB concentratie van garnalen die opgegroeid werden in conditie 1 groter is dan of gelijk is aan een willekeurige observatie van de PCB concentratie van garnalen die opgegroeid werden in conditie 2 is niet gelijk aan 50%. In een van de twee condities heb je dus een hogere kans om hogere PCB waarden te observeren.

## Uitvoeren van de test en interpretatie

Vervolgens kan men de gemiddelde rank vergelijken tussen de twee groepen aan de hand van de Wilcoxon-Mann-Whitney (WMW) test.

```{r}
wilcox.test(PCB.conc~group, data = shrimps)
```

De Wilcoxon test hierboven is standaard een exacte test (zie vierde paragraaf in sectie Details in `?wilcox.test`) en geeft dus de exacte p-waarde. De test is gebaseerd op alle mogelijke permutaties.

Voor grotere steekproeven is het niet mogelijk om alle permutaties door te rekenen. Dan kunnen we de test doen door gebruik te maken van asymptotische resultaten of door gebruik te maken van een test waarin een random aantal permutaties uitvoeren. 

Dat kan met de functie wilcox_test uit het coin package die toelaat om zowel de exacte test als een p-waarde op basis van gekozen aantal random permutaties te berekenen.

```{r}
set.seed(2021)
wilcox_test(PCB.conc~group, data = shrimps, distribution=approximate(nresample=100000))
wilcox_test(PCB.conc~group, data = shrimps, distribution="exact")
```

We vinden een exacte p-waarde van `r round(pvalue(wilcox_test(PCB.conc~group, data = shrimps, distribution="exact")),3)*100`% en kunnen op het 5% significantie-niveau dus de nulhypothese verwerpen ten voordele van de alternatieve hypothese. We kunnen dus besluiten dat de kans dat de PCB-concentratie van een willekeurige garnaal die opgegroeid werd in conditie 1 groter dan of gelijk is aan de PCB-concentratie van een willekeurige garnaal die opgegroeid werd in conditie 2 verschilt van 50%. Het is m.a.w. dus meer waarschijnlijk om hogere PCB-concentraties te observeren in garnalen in één van beide condities.

Deze kans kunnen we schatten op basis van de steekproef aan de hand van de geobserveerde teststatistiek van de wilcox.test-functie:

```{r}
n1 <- n2 <- 18 #18 observaties in elke groep
WObs <- wilcox.test(PCB.conc~group, data=shrimps)$statistic #geobserveerde teststatistiek
WObs/(n1*n2)
```

De puntschatting voor deze kans is dus `r round(WObs/(n1*n2),3)*100`%. Merk op dat we deze U-teststatistiek enkel terugvinden bij de wilcox.test functie en niet bij de wilcox_test-functie.

We interpreteren dit als volgt:
De kans dat de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 1 groter is dan of gelijk is aan de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 2 is gelijk aan `r round(WObs/(n1*n2),3)*100`%.

Of nog:

De kans dat de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 1 kleiner is dan de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 1 is gelijk aan `r round(1-WObs/(n1*n2),3)*100`%.

Deze kans verschilt significant van 50% op het 5% significantieniveau (p=0.019).


## Asymptotische testen/asymptotische p-waarde

Omdat men vroeger niet zoveel rekenkracht had om veel permutaties te berekenen, maakte men vaak gebruik van asymptotische verdelingen. Deze zullen uiteraard enkel een goede benadering opleveren wanneer het aantal gegevens voldoende groot is. Standaard geeft de wilcoxon_test functie de asymptotische p-waarde mee. Wij zullen echter meestal werken met de exacte p-waarde of een approximatieve p-waarde op basis van permutaties.

```{r}
wilcox_test(PCB.conc~group, data=shrimps)
wilcox.test(PCB.conc~group, data=shrimps, exact=F, correct=F)
```

# Conclusie

Er is een significant verschil in de distributie van de PCB concentraties van garnalen die werden opgegroeid onder twee verschillende condities (p = `r round(pvalue(wilcox_test(PCB.conc~group, data = shrimps, distribution="exact")),3)`). Het is meer waarschijnlijk om hogere PCB concentraties te observeren bij garnalen die opgegroeid zijn in conditie 2 dan in conditie 1. De puntschatting voor deze kans bedraagt `r round(1-WObs/(n1*n2),3)*100`%.