---
title: "PC-practicum 4: Project"
output: html_document
---

#Filmpje

https://www.youtube.com/watch?v=IGb5Ra48xXI&feature=youtu.be

# 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 de 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. Eé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)
#install.packages("tidyr")
library(tidyr)
#install.packages("multcomp")
library(multcomp)
library(readr)
```


# Vraag 1: Opdracht (leesopdracht)

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. Veldrond waaraan compost gemengd met biochar werd toegevoegd (cobc)

De dataset 'versgewicht_sla.txt' bevat het versgewicht in gram voor 28 slaplanten en welke behandeling ze ondergingen. Gebruik een ANOVA om na te gaan of er een effect is van de verschillende behandelingen op de plantengroei van bladgroenten. Gebruik post hoc-testen om na te gaan welke behandelingen effectief zijn. Voer deze analyse uit.

# Vraag 2: Inlezen van de data (leesopdracht)
```{r}
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()))
```

# Vraag 3: Dataverkenning (frequentietabel)

Tel het aantal metingen per behandeling aan de hand van de `count()` functie. Sla de frequentietabel op in `tabel`.
```{r}
tabel <- ...
tabel
```

# Vraag 4: Dataverkenning (boxplot)

Maak een figuur `boxplot` met daarin de boxplots met de versgewichten voor alle behandelingen, met elk van de boxplots in een andere kleur en waarin de individuele datapunten getoond worden.


```{r}
boxplot <- ... %>% ggplot(aes(...))+
  ... + # maak boxplot. Gebruik het argument outlier.shape = NA
  ... + # toon individuele punten
  theme_bw() 
boxplot
```

# Vraag 5: Conclusie dataverkenning (leesopdracht)

Uit de frequenctietabel uit vraag 3 kunnen we afleiden dat iedere behandelingsgroep evenveel metingen bevat, namelijk 7.
Uit de boxplot in vraag 4 kunnen we een aantal zaken afleiden.

- De gemeten versgewichten binnen de refoak behandelingsgroep zijn zeer vergelijkbaar met die van de controlegroep.
- De gemeten versgewichten binnen de cobc en de compost behandelingsgroepen zijn merkelijk hoger dan die binnen de controlegroep.
- De gemeten versgewichten binnen de compost behandelingsgroep liggen gemiddeld gezien iets hoger dan die binnen  de cobc groep, maar er grote overlap tussen gemeten waardes.
- Geen enkele behandelingsgroep bevat uitschieters.
- De metingen binnen iedere behandelingsgroep lijkt symmetrisch verdeeld.
- De variantie van de metingen lijkt niet sterk af te wijken tussen de behandelingsgroepen.
- De range en spreiding van de waarden tussen de verschillende groepen is zeer vergelijkbaar.

# Vraag 6: Statistische test uitvoeren: nulhypothese (meerkeuzevraag)
We kunnen een ANOVA uitvoeren om te testen of er een verschil is in gemiddeld versgewicht van de slaplanten tussen de verschillende behandelingsgroepen.

De nulhypothese kan dan geformuleerd worden als

1 Het versgewicht in gram is gelijk over de verschillende behandelingsgroepen heen.
2 Het gemiddelde versgewicht in gram is gelijk tussen ten minste 2 van de verschillende behandelingsgroepen.
3 Het versgewicht in gram tussen de verschillende behandelingsgroepen is gelijk aan 0.
4 Het gemiddelde versgewicht in gram is gelijk over de verschillende behandelingsgroepen heen.
5 Het gemiddelde versgewicht in gram tussen de verschillende behandelingsgroepen is gelijk aan 0.

# Vraag 7: Statistische test uitvoeren: alternatieve hypothese (meerkeuzevraag)

De alternatieve hypothese kan geformuleerd worden als

1 Het gemiddeld versgewicht in gram voor minstens 1 behandelingsgroep is verschillend van het gemiddelde versgewicht in gram van minstens 1 andere behandelingsgroep.
2 Het versgewicht in gram voor minstens 1 behandelingsgroep is verschillend van het gemiddelde versgewicht in gram binnen de andere behandelingsgroepen.
3 Het gemiddeld versgewicht in gram is verschillend tussen alle behandelingsgroepen.
4 Het gemiddeld versgewicht in gram is gelijk tussen alle behandelingsgroepen.

# Vraag 8: Lineair model voor ANOVA
We zullen voor de ANOVA een lineair model fitten waar het gemiddeld versgewicht in gram geschat wordt naargelang de type grond. Fit het correct lineair model. Geef tot slot een summary van dit gefitte model.

```{r}
fit <- lm(...) # lineair model
sum <- ... # summary van het gefitte model
```

# Vraag 9: Lineair model (leesopdracht)

De geschatte waarden voor de regressie parameters lijken reeds te suggereren dat er een verschil is in gemiddeld versgewicht naargelang het type grond. 
Let wel op, hier is nog geen correctie voor meervoudig testen op de p-waarden gebeurd!

In verband met de hypothesen kan men dus ook zeggen dat de nulhypothese stelt dat alle regressie parameters behalve het intercept nul zijn. De alternatieve hypothese stelt dan dat minstens 1 regressieparameter behalve het intercept verschillend is van nul.

# Vraag 10: Assumpties -- Onafhankelijkheid (leesopdracht)

Voor we de ANOVA uitvoeren zullen we nagaan of alle assumpties van de test voldaan zijn.

Uit het opgegeven studiedesign kunnen we afleiden dat de gemeten waarden onafhankelijk zijn van elkaar.

# Vraag 11: Assumpties -- Normaliteit

Ga na dat de metingen van versgewicht binnen elke behandelingsgroep normaal verdeeld zijn. Sla de plot op in `plot`.

```{r}
plot <- ... + 
  facet_wrap(~behandeling) # plot voor elke behandelingsgroep
plot
```

# Vraag 12: Assumpties -- Normaliteit, conclusie (leesopdracht)
De QQplots in figuur 2 lijken te suggereren dat de metingen binnen iedere behandelingsgroep normaal verdeeld zijn.
Echter, QQplots zullen veelal normaliteit suggereren bij weinig datapunten. 

Een tweede controle zou de distributie van de residuen van het lineair model kunnen zijn.
Deze moeten namelijk ook normaal verdeeld zijn als de data per behandelingsgroep normaal verdeeld is.

```{r}
plot(fit, which = 2,col = fit$model$behandeling)
legend('bottomright', levels(fit$model$behandeling), text.col = 1:4)
```

De qq-plot suggereert dat de residuen ietswat scheef verdeeld zijn naar links.
28 datapunten is echter nog altijd niet bijzonder veel. Deze qqplot zou men dan ook kunnen verwachten wanneer de data uit een normale distributie zou komen.
We zullen voor deze oefening aannemen dat de data normaal verdeeld is.

# Vraag 13: Assumpties -- Homoscedasticiteit (leesopdracht)
Tijdens de dataverkenning hebben we reeds opgemerkt dat de spreiding en range van de meetwaarden zeer vergelijkbaar is tussen de verschillende behandelingsgroepen (zie boxplots in vraag 4).
Daarom is de assumptie dat de varianties gelijk zijn tussen de groepen zeer aannemelijk.

Wanneer we echter de varianties per groep uitrekenen zien we dat er grote verschillen tussen de groepen (zie tabel 2).

```{r}
## bereken de variantie per behandelingsgroep. Met sapply itereer je over iedere behandelings groep. Zie google of de help als je wil weten hoe dit juist werkt. Je zou dit ook kunnen doen via een for loop.
variance_versgewicht <- sapply(levels(biochar$behandeling), 
                                function(b){var(biochar$versgewicht[b == biochar$behandeling])})
variance_versgewicht
```

Met slechts 7 meetwaarden is er weinig power om de echte variantie correct te schatten.
Dergelijk verschillen in varianties zouden we ook kunnen verwachten onder een normale verdeling wanneer we slechts 7 datapunten hebben.

Voor deze oefening zullen we veronderstellen dat alle assumpties voor een ANOVA voldaan zijn.

**Merk op** Wanneer men niet bereid is om deze assumpties te aanvaarden (bijvoorbeeld wanneer eerdere studies reeds aangetoond hebben dat de data geen normale distributie volgt), kunnen we ook een niet-parametrische test uitvoeren.
In deze testen nemen we geen distributionele assumpties (maar hebben over het algemeen minder power als de distributionele assumpties wel zouden opgaan). Hiervoor verwijzen we naar de lessen en practica over niet-parmetrische testen.

# Vraag 14: Voer de ANOVA test uit

In vraag 8 werd het lineair regressiemodel gefit.

Voer nu de ANOVA test uit. Sla het resultaat op in `anova`.
```{r}
# ANOVA test uitvoeren
anova <- ...
```


# Vraag 15: Conclusie ANOVA
De p-waarde (p<<0.001) van de ANOVA is bijzonder klein en we kunnen besluiten dat we de nulhypothese kunnen verwerpen.
Merk op dat we dit ook al konden zien in de output van `summary(fit)` (zie onderste lijn).
Het gemiddelde versgewicht van slaplanten verschilt tussen minstens twee van de bestudeerde soorten grond op het 5% significantieniveau.

Aan de hand van dit resultaat weten we echter niet bij welke soorten er een verschil optreedt, en hiervoor zal men een post-hoc analyse moeten uitvoeren. Een post-hoc analyse voert men enkel uit indien de ANOVA test significant was, en bestaat erin om paarsgewijze vergelijkingen uit te voeren tussen de groepen.

# Vraag 16: Post-hoc analyse (nulhypothese)
De nulhypothese bij iedere paarsgewijze test is dat
- er geen verschil is in gemiddelde versgewicht van de slaplanten op de twee gronden
- er een verschil is in gemiddelde versgewicht van de slaplanten op de twee gronden
- het gemiddelde versgewicht in de ene groep groter is dan het gemiddelde versgewicht in de andere groep

# Vraag 17: Post-hoc analyse (alternatieve hypothese)

De alternatieve hypothese bij iedere paarsgewijze test is dat 
- er een verschil is in gemiddelde versgewicht van de slaplanten op de twee gronden
- er geen verschil is in gemiddelde versgewicht van de slaplanten op de twee gronden
- het gemiddelde versgewicht in de ene groep groter is dan het gemiddelde versgewicht in de andere groep

# Vraag 18: Post-hoc analyse (uitvoeren test)
Voer een Post-hoc analyse uit waar je paarsgewijs test op een verschil in gemiddelde versgewicht aan de hand van een Tukey test, sla dit op als mcp. Bereken ook een betrouwbaarheidsinterval op het gemiddelde verschil voor iedere paarsgewijze test, sla dit op als CI_mcp.

```{r}
library(multcomp, quietly = TRUE)
mcp <- glht(..., linfct = mcp(behandeling = "Tukey"))
summary(mcp)
CI_mcp <- ...  # betrouwbaarheidsinterval voor iedere paarsgewijze test
```

# Vraag 19: Conclusie (leesopdracht)
Er is een extreem significant effect van het soort grond op de gemiddelde versgewicht van slaplanten en het soort grond waarin die gegroeid is (one-way ANOVA test, p«0.001). Op een globaal 5% significantieniveau vinden we significante verschillen in gemiddelde versgewicht van slaplanten tussen de gronden waaraan compost werd toegevoegd (compost, compost en biochar) en de overige gronden (controle en refoak). Het gemiddeld versgewicht van slaplanten opgegroeid in de controle grond is extreem significant verschillend van dat van planten opgegroeid in grond met compost en planten opgegroeid in grond met compost en biochar (beide aangepaste p-waarden < 0.001).  Gemiddeld is het versgewicht van slaplanten respectievelijk  20.9 gram (95% BI: [13.2 ; 28.5]) en 16.1 gram (95% BI: [8.5 ; 23.8]) lager in de controle grond dan in de grond met compost en in de grond met compost en biochar. Het gemiddeld versgewicht van slaplanten opgegroeid in de grond met biochar is extreem significant verschillend van dat van planten opgegroeid in grond met compost en planten opgegroeid in grond met compost en biochar (Tukey test, beide aangepaste p-waarden < 0.001).  Gemiddeld is het versgewicht van slaplanten respectievelijk 22.7 gram (95% BI: [15.0 ; 30.3]) en 18.0 gram (95% BI: [10.4; 25.6]) lager in de grond met biochar dan in de grond met compost en in de grond met compost en biochar. Het verschil in gemiddelde versgewicht van slaplanten in grond met biochar versus in controle grond, en,  in grond met compost versus grond met compost en biochar is niet significant op het 5% significantie-niveau (p-waarden respectievelijk p = 0.91 en p = 0.34).

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