---
title: "Niet-parametrische testen: Plantengroei in bladgroenten"
output:
    html_document:
      code_download: true    
      theme: cosmo
      toc: true
      toc_float: true
      highlight: tango
      number_sections: true
---


# Situering

In de landbouw is het belangrijk om een goede productie van gewassen te bekomen. Voor het kweken van bladgroenten zoals sla en spinazie houdt dit in om zo groot mogelijke kroppen of bladeren te bekomen. De consument zal namelijk eerder kiezen voor grote, volle kroppen met veel blad.

Om dit te bekomen, gebruiken boeren vaak verschillende types van mest. Er is meer en meer een tendens om over te schakelen van kunstmeststoffen naar meer duurzame organische meststoffen. De meest gekende organische meststof is compost. In het ILVO worden er echter ook nieuwe types van organische meststoffen getest. Één van deze nieuwe stoffen is biochar. Biochar wordt gevormd tijdens een pyrolyse proces van biomassa (zoals bv. houtafval) waarbij energie wordt opgewekt. Het restmateriaal van de pyrolyse wordt de biochar genoemd, een stof die sterk gelijkt op houtskool maar nuttige eigenschappen heeft zoals vasthouden van het water in de bodem en het beïnvloeden van de nuttige bacteriën in de bodem.

```{r}
library(ggplot2)
library(dplyr)
library(coin)
library(tidyr)
library(readr)
```

## Data

De onderzoekers willen nagaan of biochar, compost en compost gemengd met biochar invloed heeft op de groei van bladgroenten. Daarvoor groeiden ze sla op in potten met veldgrond in een groeikamer. Het versgewicht van de plant werd na acht weken gemeten. Versgewicht wordt gebruikt als een maat voor de plantengroei. De planten zijn in vier verschillende types grond opgegroeid:
 
1. Veldgrond (controle) 
2. Veldgrond waaraan biochar werd toegevoegd (refoak)
3. Veldgrond waaraan compost werd toegevoegd (compost)
4. Veldgrond waaraan compost gemengd met biochar werd toegevoegd (cobc)

Het bestandje 'versgewicht_sla.txt' bevat het versgewicht in gram voor 28 slaplanten en welke behandeling ze ondergingen.


```{r}
#Lees de data in
biochar <- read_csv("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum5/versgewicht_sla.txt", col_types = list(col_integer(), col_factor(levels = c("cobc","compost","controle","refoak")), col_integer()))
```

## Onderzoeksvraag
Ga na of er een effect is van de verschillende behandelingen op de plantengroei van bladgroenten.
Indien dit zo is, welke behandelingen zijn effectief?


# Vraag 1: Data-exploratie
Tel het aantal metingen per behandeling, gebruik de count functie.

```{r}
count1<-...
count1
```

# Vraag 2: Data-exploratie
Maak een vergelijkende boxplot voor de versgewichten voor de verschillende behandelingen. Geef elke boxplot een andere kleur en plot de individuele datapunten. Gebruik voor deze **en alle andere plots in dit practicum** de theme_bw()-functie voor de layout van de plots. Gebruik ook het argument `outlier.shape = NA`.


```{r}
boxplot1<-...

boxplot1
```

# Vraag 3: Data-exploratie

Maak een QQ-plot met lijn voor elke groep, gebruik makend van de facet_wrap()-functie.

```{r}
#QQplots voor elke groep apart
qqplot2<-...
qqplot2
```

# Vraag 4: Data-exploratie
Maak een lineair model om de normaliteit van de residuen na te gaan.

```{r}
#QQplot residuen
model1<-lm(..)
plot(model1, which = 2)
```

# Vraag 5: Data-exploratie
Wat kan je **niet** besluiten uit de data-exploratie?

1. Er zijn slechts 7 observaties per groep.
2. De gemeten versgewichten binnen de cobc en de compost behandelingsgroepen zijn merkelijk hoger dan die binnen de controlegroep.
3. Geen enkele behandelingsgroep bevat uitschieters.
4. De variantie van de metingen lijkt niet sterk af te wijken tussen de behandelingsgroepen.
5. De gegevens in elke groep komen uit een normale verdeling


# Vraag 6: Keuze test (Leesopdracht)
Vorige les losten jullie deze onderzoeksvraag op met behulp van parametrische ANOVA (de F-test).
Dit is een geldige test als aan een aantal assumpties voldaan worden die geverifieerd moeten worden in de dataset. Echter, hoe kleiner de dataset hoe moeilijker het is om deze assumptie te controleren. Laten we als voorbeeld kijken naar de asumptie van normaliteit.

```{r}
plot(model1, which = 2)
qqplot2
```

De QQ-plot is nogal twijfelachtig, en men zou beide kanten kunnen beargumenteren.
Inderdaad, enerzijds kan men argumenteren dat de QQ-plot niet bijzonder sterk afwijkt van wat men zou verwachten op basis van een Normale distributie. Dit geldt zowel voor de QQ-plot van de residuen als voor de QQ-plots in de aparte groepen.
Anderzijds, kan men argumenteren dat (i) de QQ-plot bij de residuen wel degelijk licht afwijkt, en (ii) de QQ-plots binnen elke groep weinig bewijs leveren voor normaliteit, aangezien er zo weinig datapunten in elke groep zitten.

Ter vergelijking zullen we onze steekproef vergelijken met simulaties uit een normale verdeling en een gamma verdeling. Een gamma distributie heeft een aantal eigenschappen die een betere zijn om gewichten te modelleren. Bijvoorbeeld, gewichten kunnen niet negatief zijn, een normale distributie kan tot min oneindig lopen waar een gamma distributie een ondergrens heeft bij nul. 

We tonen eerst visueel het verschil tussen de normale verdeling en de gamma verdeling

```{r}
x= seq(-1,10,.01)
data1<-data.frame(x ,normaal = dnorm(x,2),gamma = dgamma(x,5,2))
data2<-data1 %>% pivot_longer(cols = c("normaal","gamma"),
                              names_to = "verdeling", 
                              values_to = "dichtheid")
density_plot<-data2%>%ggplot(aes(x=x, y=dichtheid, col=verdeling))+
  geom_line()+
  geom_line(aes(x=0), color="black")
density_plot
```

We zullen nu zowel steekproeven genereren uit een normale verdeling en een gamma verdeling met dezelfde grootte als onze oorspronkelijke steekproef.
```{r}
set.seed(1)
#Normale verdeling
data1<-data.frame(versgewicht_sim=rnorm(28),behandeling_sim=as.factor(rep(1:4,7)))
model1<-lm(versgewicht_sim~behandeling_sim, data=data1)
qqplot1_sim<-ggplot(mapping=aes(sample=model1$residuals))+
  geom_qq()+geom_qq_line()+theme_bw()
qqplot1_sim
qqplot2_sim<-data1%>%ggplot(aes(sample=versgewicht_sim))+
  geom_qq()+geom_qq_line()+facet_wrap(~behandeling_sim)+theme_bw()
qqplot2_sim

#Gamma verdeling
data1<-data.frame(versgewicht_sim=rgamma(28,5,2),behandeling_sim=as.factor(rep(1:4,7)))
model1<-lm(versgewicht_sim~behandeling_sim, data=data1)
qqplot1_sim<-ggplot(mapping=aes(sample=model1$residuals))+
  geom_qq()+geom_qq_line()+theme_bw()
qqplot1_sim
qqplot2_sim<-data1%>%ggplot(aes(sample=versgewicht_sim))+
  geom_qq()+geom_qq_line()+facet_wrap(~behandeling_sim)+theme_bw()
qqplot2_sim
```

Uit de simulaties blijkt dat we onze gegevens niet kunnen onderscheiden van normaal verdeelde gegevens, maar bijvoorbeeld ook niet van gamma verdeelde gegevens. We hebben dus geen bewijs dat de data niet normaal verdeeld zijn, maar zeker ook geen bewijs dat ze wél normaal verdeeld zijn. Aangezien we hier ook te weinig gegevens hebben om gebruik te maken van de Centrale limietstelling, zijn we niet zeker of de resultaten van ANOVA betrouwbaar zullen zijn.

Zoals je ziet, is dit een dataset waarbij zowel de argumentatie voor een parametrische ANOVA test als de argumentatie voor een niet-parametrische test steek kan houden. 

# Vraag 7: Keuze test (Meerkeuzevraag)

We hebben gezien dat men kan argumenteren dat normaliteit mogelijks niet opgaat en we hebben ook te weinig data om de assumptie van homoscedasticiteit na te gaan. Welke test zullen we kiezen als we willen vermijden om assumpties te doen die we niet goed kunnen nagaan?

1. ANOVA
2. Een permutatietest met de ANOVA F-teststatistiek
3. Een Kruskal-Wallis test met een interpretatie op basis van de probabilistische index
4. Een Kruskal-Wallis test met een interpretatie op basis van de mediaan


# Vraag 8: Voorwaarden (Meerkeuzevraag)

Welke voorwaarden moeten voldaan zijn voor de Kruskal-Wallis test 

1. Er zijn geen voorwaarden voor deze test.
2. De gegevens moeten onafhankelijk zijn.
3. Er mag enkel sprake zijn van een locatie-shift.
4. De gegevens moeten onafhankelijk zijn en er mag enkel sprake zijn van een locatie-shift.

# Vraag 9: Voorwaarden (leesopdracht)

Onder de locatie shift assumptie dat de verschillende groepen onderling nog steeds dezelfde distributie volgen (dit hoeft uiteraard niet noodzakelijk een normale verdeling te zijn, elke verdeling is goed, zolang het maar in beide groepen dezelfde verdeling is),zouden we de resultaten van de rank sum test kunnen interpreteren in termen van een verschil in mediaan. Merk op, dat dat opnieuw moeilijk kan worden geëvalueerd op basis een kleine dataset.

Bovendien is het bij rangtesten niet nodig om deze aanname te maken. Daarom zullen we de resultaten van de test steeds interpreteren in termen van probabilistische indices, de kans dat een willekeurige meetwaarde uit de ene groep groter is dan of gelijk aan ("$\geq$") een willekeurige meetwaarde uit een andere groep.

# Vraag 10: Hypotheses (Meerkeuzevraag)

Welke nulhypothese is correct geformuleerd?

1. De nulhypothese stelt dat de gemiddeldes van de versgewichten van slaplanten gelijk zijn in de vier types grond.
2. De nulhypothese stelt dat de medianen van de versgewichten van slaplanten gelijk zijn in de vier types grond.
3. De nulhypothese stelt dat de verdelingen van de versgewichten van slaplanten gelijk zijn in de vier types grond.


# Vraag 11: Hypotheses (Meerkeuzevraag)

1. De alternatieve hypothese stelt de gemiddeldes van de versgewichten van slaplanten verschillen  tussen de vier behandelingen.
2. De alternatieve hypothese stelt de medianen van de versgewichten van slaplanten verschillen  tussen de vier behandelingen.
3. De alternatieve hypothese stelt dat de kans dat het versgewicht van een willekeurige slaplant uit minstens één van de vier soorten grond groter is dan of gelijk is aan  ("$\geq$") het versgewicht van een willekeurige slaplant uit minstens één van de andere soorten grond verschilt van 50%.


# Vraag 12: Test uitvoeren

We zullen de test zowel eens uitvoeren met de kruskal.test() functie als de kruskal_test functie. Je hoeft, behalve de formule en de data, nog geen andere argumenten in te geven. 
```{r}
test1<-kruskal.test(..., data=...)
test1
test2<-kruskal_test(..., data=...)
test2
```

# Vraag 13: Test uitvoeren (Leesopdracht)

We zien dat we twee keer quasi dezelfde output krijgen. De output bij de tweede functie geeft echter aan dat het hier (bij allebei!) gaat om een asymptotische test. Omdat we te maken hebben met een klein aantal observaties per groep is het echter niet aan te raden om de asymptotische KW-test te gebruiken om een p-waarde te berekenen. Deze stelt namelijk dat je teststatistiek onder de nulhypothese een $\chi^2$ distributie volgt (de nuldistributie) wanneer het aantal oberservaties per groep naar oneindig gaat (dus voor een groot aantal observaties). Dit is hier duidelijk niet het geval.

# Vraag 14: Test uitvoeren 

Wanneer je een beperkt aantal observaties hebt, is het beter de nuldistributie exact te berekenen via permutaties. Dit is uiteraard computationeel onmogelijk, want met 4 groepen van elk 7 observaties bestaan er 4.725e+14 mogelijke permutaties. Gelukkig is het meestal voldoende om de exacte nuldistributie te benaderen met een kleiner aantal permutaties.

Voer een KW-test uit op de data met behulp van 10000 permutaties. Wat kan je besluiten uit het resultaat?

```{r}
set.seed(0)
test3 <- kruskal_test(...,data=...,distribution=...)
test3
```

# Vraag 15: Conclusie (leesopdracht)

Het type bemesting heeft een extreem significant effect op de distributie van het versgewicht van slaplanten (p << 0.001). In minstens één van de vier types grond is er dus een hogere kans op slaplanten met een hoger versgewicht dan in één van de overige types grond.
 
 **Merk op:** doordat we de p-waarde benaderd hebben aan de hand van simulaties, kan de p-waarde verschillen indien we de code opnieuw uitvoeren (en dus ook met je collega naast jou). Merk op dat deze variatie in p-waarde zou moeten verminderen als we het aantal permutaties zouden opdrijven, omdat we dan dichter in de buurt komen bij de exacte nuldistributie. 


# Vraag 16: Post-hoc analyse (Meerkeuzevraag)

Als we willen weten in welke soorten grond er nu verschillen in de distributie van het versgewicht voorkomen, zal men een post-hoc analyse moeten doen. Zoals bij de parametrische ANOVA zal je ook bij de niet-parametrische Kruskal-Wallis test paarsgewijze testen moeten uitvoeren om onderlinge verschillen in versgewicht tussen de slaplanten te detecteren.

Op welke manier zou jij een post-hoc analyse uitvoeren?

1. Paarsgewijze t-testen met Welch-benadering en gecorrigeerde p-waarden
2. Paarsgewijze wilcoxon-testen met asymptotische en gecorrigeerde p-waarden
3. Paarsgewijze wilcoxon-testen met exacte en gecorrigeerde p-waarden
4. Paarsgewijze wilcoxon-testen met approximatieve en gecorrigeerde p-waarden


# Vraag 17: Post-hoc analyse (Leesopdracht)

Dit zullen uiteraard ook niet-parametrische testen zijn zoals de Wilcoxon-Mann-Whitney Test (WMW test). Hoewel het bij de Kruskal-Wallis test niet mogelijk was om een exacte p-waarde te berekenen, zal dit wel mogelijk zijn bij de Wilcoxon testen.

Voor elke paarsgewijze wilcoxon test gelden de volgende hypothesen:

* De nulhypothese stelt dat de verdeling van het versgewicht van slaplanten gelijk is voor beide types grond.
 
* De alternatieve nulhypothese stelt dat de kans, dat het versgewicht van een willekeurige slaplant opgegroeid in de ene soort grond hoger is dan of gelijk is aan ("$\geq$") het versgewicht van een willekeurige slaplant uit de andere soort grond, verschillend is van 50%.


# Vraag 18: Post-hoc analyse
Voer een WMW post-hoc analyse uit op de data.

We doen eerst een paarsgewijze WMW test zoals beschreven in de cursus. (Het is normaal dat je een warning krijgt bij het uitvoeren van de test.)
```{r}
WMW_test <- pairwise.wilcox.test(...,...)
WMW_test
```

# Vraag 19: Post-hoc analyse
We zien in de output een foutmelding: `cannot compute exact p-value with ties`. Dit is omdat de `pairwise.wilcox.test()` functie de standaard `wilcox.test()` functie gebruikt in R. In het helpbestand van deze functie zien we dat deze functie standaard een exacte test doet maar in de aanwezigheid van ties dit niet kan doen en terugrijpt naar een asymptotische test.

We kunnen echter wel exacte p-waarden bekomen met de `wilcox_test()` uit de `coin` package voor iedere combinatie van behandeling. Vul de argumenten van de test aan in onderstaande code en probeer ook de overige code goed te begrijpen.

**Hiervoor zullen we dus een exacte test gebruiken aangezien het aantal combinaties niet te groot is!**

**Lees de onderstaande code aandachtig, welke data moeten we gebruiken bij de wilcox_test? We filteren de data zodat we telkens maar twee behandelingen vergelijken, dus niet simpelweg biochar ingeven!**

```{r}
behandelingen<-levels(biochar$behandeling)  #de verschillende behandelingen
pairwise_comparisons<-combn(behandelingen, 2)  #we kiezen telkens 2 behandelingen uit
pairwise_comparisons   #Elke kolom staat voor één paarsgewijze vergelijking.
n<-ncol(pairwise_comparisons)   #aantal paarsgewijze vergelijkingen

pvalues<-c() #lege vector creëren voor de p-waarden
names_comparisons<-c() #lege vector creëren voor de namen van de vergelijkingen
for(i in 1:n){
  behandelingen_pc<-pairwise_comparisons[,i] # twee behandelingen nemen
  data_pc<-biochar%>%filter(behandeling %in% behandelingen_pc) #de bijbehorende data selecteren
  test<-wilcox_test(..., data=..., distribution=...) #de test uitvoeren
  pvalue1<-pvalue(test) #de p-waarde opslaan ..
  pvalues<-c(pvalues, pvalue1) #.. en toevoegen aan onze vector
  name1<-paste0(behandelingen_pc, collapse="_vs_") # ook de namen van de vergelijking opslaan..
  names_comparisons<-c(names_comparisons, name1) #.. en toevoegen aan de vector met namen
}
names(pvalues)<- names_comparisons
pvalues
```

# Vraag 20: Post-hoc analyse (Leesopdracht)
De gevonden p-waarden zijn nog niet gecorrigeerd voor multiple testing. Dit kunnen we doen met de `p.adjust()` functie.
```{r}
pvalues_holm = p.adjust(pvalues,method = 'holm')
pvalues_holm
```
We zien inderdaad dat de exacte p-waarden afwijken van de p-waarden berekend door de `pairwise.wilcox.test()` functie. We zullen verder werken met deze laatste p-waarden.
We zien dat er alleen significante verschillen zijn in de distributie van het versgewicht tussen de slaplanten die gegroeid zijn in de behandelingen met compost in vergelijking met planten die werden opgegroeid in de behandelingen zonder compost.


# Vraag 21: Post-hoc analyse

Kan je nu ook aan de hand van de Mann-Whitney teststatistiek een puntschatting geven voor de probabilistische index?**

Je zal in de output opnieuw de foutmelding: `cannot compute exact p-value with ties` zien, omdat we terug de `wilcox.test()` functie gebruiken. Maar dit is hier geen probleem omdat we niet de p-waarde nodig hebben maar de U1 statistiek. U1 is het aantal keer dat een observatie uit de eerste groep groter of gelijk is aan een observatie uit de tweede groep en is dus ook gedefinieerd bij ties (zie ook cursus).

Vul in onderstaande code de formule in om de probabilistische index te berekenen. Merk op dat de groepsgroottes overal gelijk zijn (in elke groep zitten 7 observaties).

```{r}
probInds<-c()
ng<-7 # de grootte van elke behandelingsgroep is 7
for(i in 1:n){
  behandelingen_pc<-pairwise_comparisons[,i] # twee behandelingen nemen
  data_pc<-biochar%>%filter(behandeling %in% behandelingen_pc) #de bijbehorende data selecteren
  U1 <- wilcox.test(versgewicht~behandeling,data_pc)$statistic #U-teststatistiek opslaan
  probInd1<-...   #de probabilistische index berekenen
  probInds<-c(probInds, probInd1)
}
names(probInds) <- names_comparisons
probInds
```


# Vraag 22: Post-hoc analyse (Leesopdracht)

**Hoe interpreteer je de resultaten van de WMW testen?**

We vinden significante verschillen in de verdelingen van versgewicht van slaplanten opgegroeid in beide behandelingen met compost in vergelijking met de controlebehandeling.

Het is meer waarschijnlijk om slaplanten met een hoger versgewicht te observeren in de compostbehandeling dan in controlebehandeling. De puntschatting op deze kans is 100%  (WMW test, aangepaste p-waarde < 0.01).

Het is meer waarschijnlijk om slaplanten met een hoger versgewicht te observeren in de compost met biochar behandeling dan in controlebehandeling. De puntschatting op deze kans is 98.0%  (WMW test, aangepaste p-waarde < 0.01).

Bovendien vinden we ook significante verschillen in de verdelingen van versgewicht voor beide behandelingen met compost in vergelijking met de behandeling met enkel biochar.

Het is meer waarschijnlijk om slaplanten met een hoger versgewicht te observeren in de compostbehandeling dan in de biochar behandeling. De puntschatting op deze kans is 100%  (WMW test, aangepaste p-waarde < 0.01).
Het is meer waarschijnlijk om slaplanten met een hoger versgewicht te observeren in de compost met biochar behandeling dan in de behandeling met enkel biochar. De puntschatting op deze kans is 100%  (WMW test, aangepaste p-waarde < 0.01).

De verdelingen van versgewicht van slaplanten in de compost behandeling en in de compost met biochar behandeling zijn niet significant verschillend op het 5% significantie-niveau (WMW test, aangepaste p-waarde 0.39).

De verdelingen van versgewicht van slaplanten in de controlebehandeling en in de behandeling met enkel biochar zijn niet significant verschillend op het 5% significantie-niveau (WMW test, aangepaste p-waarde 0.41).


# Vraag 23: Conclusie (Leesopdracht)

Er is een extreem significant effect van het type grond op de verdeling van het versgewicht van slaplanten (p << .001).

Het is meer waarschijnlijk om slaplanten met een hoger versgewicht te observeren in de compostbehandeling en de compost met biochar behandeling dan in de controlebehandeling (beide p-waarden p < 0.01). De puntschattingen op deze kansen zijn respectievelijk 100% en 98%.  

Het is meer waarschijnlijk om slaplanten met een hoger versgewicht te observeren in de compostbehandeling en de compost met biochar behandeling dan in behandeling met enkel biochar (beide p-waarden p < 0.01). De puntschattingen op deze kansen zijn beiden 100%.  

Op het 5% significantie niveau is de verdeling van versgewicht van slaplanten niet verschillend tussen de compost behandeling vs de compost met biochar behandeling en tussen de controle behandeling vs de behandeling met enkel biochar (respectievelijke p-waarden 0.39 en 0.41).

We kunnen besluiten dat behandeling van grond met compost of met compost en biochar een positieve invloed heeft op de groei van slaplanten. Aangezien in deze studie alleen gekeken werd naar slaplanten, kunnen we dit uiteraard geen conclusies trekken voor bladgroenten in het algemeen.

