Niet-parametrische testen: Plantengroei in bladgroenten

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. Éé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.

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.

Onderzoeksvraag

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

Waarom een niet-parametrische test?

Voor het project lossen 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.

#Lees de data in
dataset=read.csv("~/Dropbox/PhD/Onderwijs/Statistiek Biochemie/201819/dropboxStats1819/class6_nonPar/half/versgewicht_sla.txt",header=TRUE)
m <- lm(versgewicht ~ behandeling, data=dataset)
par(mfrow=c(2,2))
plot(m)
par(mfrow=c(1,1))

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. Anderzijds, kan men argumenteren dat (i) de QQ-plot licht afwijkt, en (ii) we bovendien een Normale distributie moeten veronderstellen binnen elke groep, en de bovenstaande QQ-plot geen sluitend bewijs geeft. Zoals je ziet, is dit een dataset waarbij zowel de argumentatie voor een parametrische ANOVA test als de argumentatie voor een niet-parametrische Kruskal-Wallis test steek kan houden. De veilige optie is duidelijk de Kruskal-Wallis test, aangezien die geen parametrsiche assumpties veronderstelt. De enige veronderstelling is dat de data onafhankelijk zijn. Bovendien kan het, om een algemeen inzicht te verkrijgen, ook nuttig zijn om de resultaten te interpreteren in functie van de probabilistische index.

Stel dat de versgewichten van de planten geen Normale maar bvb een gamma distributie volgen. Een gamma distributie heeft een aantal eigenschappen die een betere keuze zou maken on gewichten te modelleren. Bijvoorbeeld, gewichten kunnen niet negatief zijn, een normalie distributie kan tot min oneindig lopen waar een gamma distributie een ondergrens heeft bij nul.

x = seq(-1,10,.01)
ynorm = dnorm(x,2)
ygam = dgamma(x,5,2)
plot(x,ynorm,type = 'l')
lines(x,ygam,col = 'red')

abline(v = 0)
legend('topright',c('Normale', 'Gamma'),text.col = c('black','red'))

We zullen 9 samples van elk 7 datapunten genereren uit een normale distributie.

par(mfrow = c(3,3))
set.seed(11) ## zorgt ervoor dat bij herlopen van de code exact dezelfde 'random' waarden genereerd worden.
for(i in 1:9){
  x = rnorm(7)
  qqnorm(x, main = i)
  qqline(x)
}

We zien in de QQ-plots dat de meeste samples normaal verdeeld lijken, terwijl sommigen toch niet-lineaire patronen vertonen (zoals QQ-plot 5) en dus misschien (incorrect) als niet normaal beoordeeld zouden worden.

We zullen nu 9 samples van elk 7 datapunten genereren uit een gamma distributie.

par(mfrow = c(3,3))
set.seed(11)
for(i in 1:9){
  x = rgamma(7,5,2)
  qqnorm(x, main = i)
  qqline(x)
}

We zien in de QQ-plots (bijvoorbeeld QQ-plot 4) dat sommige stalen samples afwijken van een rechte wat correct tot de conlusie van niet-normaal verdeelde data zou leiden. Aan de andere kant hebben de datapunten in bvb QQ-plot 7 wel een linear verband en zou deze data incorrect wel als normaal verdeeld beschouwd kunnen worden.

Nu is de assumptie van normaliteit belangrijk bij ANOVA (en ook t-test) omdat men hier ervan uitgaat dat je steekproefgemiddelde normaal verdeeld is rond het populatiegemidelde. Als je data normaal verdeeld is, zal dit altijd het geval zijn. Als je data niet normaal verdeeld is, dan gaat dit alleen maar op als je dataset groot genoeg is (zie de Centrale Limietstelling in de cursus).

We illustreren dit met een simulatie. Een groot aantal datasets van elk 4 datapunten worden gegenereerd uit een normale verdeling. We bereken per dataset het gemiddelde en we construeren een QQ-plot van deze gemiddelden om te controleren of ze normaal verdeeld zijn. We herhalen dit, maar nu worden de 4 datapunten gegenereerd uit een Gamma distributie. Om het effect van dataset grootte te bestuderen zullen we ook datasets uit een Gamma distributie genereren van grootte 7 en 10.

set.seed(11)
samples_normaal = replicate(30000, mean(rnorm(4)))
samples_gamma_4 = replicate(30000, mean(rgamma(3,5,2)))
samples_gamma_7 = replicate(30000, mean(rgamma(7,5,2)))
samples_gamma_10 = replicate(30000, mean(rgamma(30,5,2)))
par(mfrow = c(2,2))
qqnorm(samples_normaal,pch = 20, main = 'normaal; n = 4')
qqline(samples_normaal, col = 'red')

qqnorm(samples_gamma_4, pch = 20, main = 'gamma; n = 4')
qqline(samples_gamma_4, col = 'red')

qqnorm(samples_gamma_7,pch = 20, main = 'gamma; n = 7')
qqline(samples_gamma_7, col = 'red')

qqnorm(samples_gamma_10,pch = 20, main = 'gamma; n = 30')
qqline(samples_gamma_10, col = 'red')

We zien in de QQ-plot linksboven dat de gemiddelden van de normaal verdeelde datasets inderdaad op een rechte liggen en dus normaal verdeeld lijken. Bij de gemiddelden van de gamma datasets met grootte 4 zien we duidelijke afwijken van de normaliteit (de distributie is scheef verdeeld naar links). Teststatistieken en p-waarden berekend met ANOVA en t-testen zijn in dat geval ongeldig.

In de onderste rij QQ-plots kan men zien dat hoe groter de steekproeven worden, hoe meer de gemiddelden een normale distributie benaderen. Bij grotere datasets zal dus men kunnen argumenteren dat het steekproefgemiddelde een normale verdeling volgt. Dit effect, waarbij de gemiddelden van een steekproef uit om het even welke distributie een normale distributie gaan benaderen als de steekproef maar groot genoeg is, noemt men de centrale limietstelling. Als je steekproef groot genoeg is, kan je dus de centrale limietstelling inroepen om de normaliteitsassumptie te aanvaarden.

Als er geen sterke argumenten om de normaliteitsassumptie te aanvaarden (bijvoorbeeld gelijkaardige grotere datasets in het verleden lijken normaal verdeeld), is het doorgaans veiliger om dit niet te doen en een niet-parametrische test te gebruiken.

Merk op: ook als de assumptie van homoscedasticiteit geschonden is of moeilijk kan worden nagegaan door een te kleine steekproef, zullen de teststatistieken en p-waarden vertekend zijn. Voor een t-test kan je nog gebruik maken van de Welch t-test, voor een ANOVA zal je bij zeer sterke afwijkingen van homoscedasticiteit ook een niet-parametrische test moeten gebruiken en je resulaten in termen van de probabilistische index moeten interpreteren.

Data analyse

We zullen nu de slaplanten dataset herevalueren en niet-parametrische testen overwegen.

#Lees de data in
dataset=read.csv("~/Dropbox/PhD/Onderwijs/Statistiek Biochemie/201819/dropboxStats1819/class6_nonPar/half/versgewicht_sla.txt",header=TRUE)

Dataverkenning

## Tel het aantal metingen per behandeling
table(dataset$behandeling)

    cobc  compost controle   refoak 
       7        7        7        7 

Tabel 1: frequentietabel van het aantal observaties per behandelingsgroep.

boxplot(versgewicht~behandeling,dataset)
set.seed(10)
stripchart(versgewicht~behandeling, data = dataset,
            vertical = TRUE, method = "jitter", 
            pch = 19, col =c("coral","steelblue","gold","aquamarine3"), 
            add = TRUE)

Figuur 1: Boxplot van het versgewicht in gram per behandelingsgroep. De punten zijn de individuele metingen.

In de frequentietabel (Tabel 1) zien we dat er slechts 7 observaties per groep zijn. 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 in de cobc groep, maar er is een grote overlap tussen gemeten waardes. Geen enkele behandelingsgroep bevat uitschieters. De metingen binnen iedere behandelijksgroep lijkt symmetrisch verdeeld. De variantie van de metingen lijkt niet sterk af te wijken tussen de behandelingsgroepen. Het bereik en de spreiding van de waarden tussen de verschillende groepen is zeer vergelijkbaar.

Kruskal-Wallis Rank Test

Door de weinig datapunten is het moeilijk om de distributionele assumpties voor de parametrische ANOVA te controleren (normaliteit en gelijke variantie tussen de groepen). De Kruskal-Wallis Rank Test (KW-test) is een niet-parametrische variant van de ANOVA. Deze test laat ons toe om de distributionele assumpties te relaxeren. Let erop dat de assumptie van onafhankelijkheid tussen de observaties nog steeds moet opgaan!

Als we willen testen voor een verschil tussen de mediaan van de verschillende behandelingsgroepen, dan moeten we uitgaan van de locatie shift assumptie. Onder de locatie shift assumptie veronderstellen we 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).

Wanneer je kijkt naar de boxplots in Figuur 1, denk je dan dat de locatie shift assumptie aannemelijk is?

De range en spreiding van de waarden tussen de verschillende groepen is op het eerste zicht zeer vergelijkbaar, maar je ziet in de boxplots wel grote verschillen in interquantiel afstanden (bvb tussen COBC en refoak boxplot). Eigenlijk hebben we te weinig data om deze assumptie te verifieren.

Wanneer we niet bereid zijn om de locatie shift assumptie te aanvaarden, dan kunnen we de distributionele assumpties nog verder relaxeren, en testen we in termen van probabilistische indexen. Hier testen we in het geval van twee groepen de kans dat een willekeurige meetwaarde uit de eerste groep komt groter of gelijk (“\(\geq\)”) is dan een willekeurige meetwaarde uit de tweede groep.

Hypothese

Probeer nu eens de nulhypothese op te stellen in termen van de onderzoeksvraag?

De nulhypothese stelt dat de verdelingen van de versgewichten van slaplanten gelijk zijn in de vier types grond.

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

Statistische test

Sinds we te maken hebben met een klein aantal observaties per groep is het 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 neigt (dus voor een groot aantal observaties).

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.725183510^{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 een groot aantal permutaties (bv. 10000). (Tip: zie cursus http://users.ugent.be/~lclement/statistiek/ voor code voorbeelden) Wat kan je besluiten uit dit resultaat?

library(coin)
kwPerm <- kruskal_test(versgewicht~behandeling,dataset, distribution=approximate(B=100000))
kwPerm

    Approximative Kruskal-Wallis Test

data:  versgewicht by behandeling (cobc, compost, controle, refoak)
chi-squared = 20.715, p-value < 2.2e-16

We zien een extreem significant effect (p << .001) van het type bemesting op het versgewicht van slaplanten. Op het 5% significantieniveau stellen we datde kans dat een willekeurige slaplant uit minstens één van de vier types grond een groter dan of gelijk (“\(\geq\)”) versgewicht heeft vergeleken met een willekeurige slaplant uit minstens één van de andere types grond verschilt van 50%.

Merk op: doordat we de p-waarde benaderd hebben aan de hand van simulaties, kan de p-waarde 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.

Als we willen weten in welke soorten grond er nu juist een significant verschil in versgewicht optreedt, zal men een post-hoc analyse moeten doen.

Post-hoc analyse

Zoals bij de parametrische ANOVA zal je ook bij de niet-parametrische KW test paarsgewijze testen moeten uitvoeren om onderlinge verschillen in versgewicht tussen de slaplanten te detecteren. Dit zullen uiteraard ook niet-parametrische testen zijn zoals de Wilcoxon-Mann-Whitney Test (WMW test). Als je de keuze gemaakt hebt om bij de KW test de locatie shift assumptie te verwerpen dan zal je dit hier terug doen en de resultaten in termen van de probabilistische index interpreteren.

Wat zal nu de nul en alternatieve hypothese zijn voor iedere paarsgewijze test?

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 uit 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%.

Voer een WMW post-hoc analyse uit op de data.

We doen een paarsgewijze WMW test zoals beschreven in de curus.

pairwise.wilcox.test(dataset$versgewicht,dataset$behandeling)

cannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with ties

    Pairwise comparisons using Wilcoxon rank sum test 

data:  dataset$versgewicht and dataset$behandeling 

         cobc  compost controle
compost  0.400 -       -       
controle 0.013 0.013   -       
refoak   0.013 0.013   0.400   

P value adjustment method: holm 

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 de help 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 en vervolgens te corrigeren voor multiple testing met de p.adjust() functie.

library(coin)
## tel het aantal observaties per groep
nGroep <- table(dataset$behandeling)
## bereken met de wilcoxon test een p-waarde voor iedere combinatie van grondsoorten
behandelingen = (levels(dataset$behandeling))
pvalues <- combn(behandelingen,2,function(x){
  ## Doe een paarsgewijze WMW test.
  test = wilcox_test(versgewicht~behandeling,subset(dataset,behandeling%in%x), distribution = 'exact')
  ## Neem de pwaarde uit output van de test
  pvalue(test)
})
## Doe een correctie voor multiple testing met Holm (je kan eventueel ook kiezen voor Bonferroni etc ...)
pvalues_holm = p.adjust(pvalues,method = 'holm') 
 ## link de pwaarde met de juiste paarsgewijze test
names(pvalues_holm) <- combn(levels(dataset$behandeling),2,paste,collapse="_VS_")
pvalues_holm

    cobc_VS_compost    cobc_VS_controle      cobc_VS_refoak compost_VS_controle   compost_VS_refoak  controle_VS_refoak 
        0.393939394         0.005244755         0.003496503         0.003496503         0.003496503         0.405594406 

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 tussen de slaplanten die gegroeid zijn in compost ten opzichte van de planten die zonder compost opgegroeid zijn.

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

## tel het aantal observaties per groep
nGroep <- table(dataset$behandeling)
## bereken een probabilistische index voor iedere mogelijke combinatie van grondsoorten
behandelingen = (levels(dataset$behandeling))
probInd <- combn(behandelingen,2,function(x){
  ## Berkeken de U1 statistiek
  U1 = wilcox.test(versgewicht~behandeling,subset(dataset,behandeling%in%x))$statistic
  ## Bereken probabilistische index
  U1/prod(nGroep[x])
})

cannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with ties

names(probInd) <- combn(levels(dataset$behandeling),2,paste,collapse="_VS_")
probInd

    cobc_VS_compost    cobc_VS_controle      cobc_VS_refoak compost_VS_controle   compost_VS_refoak  controle_VS_refoak 
          0.2857143           0.9795918           1.0000000           1.0000000           1.0000000           0.6428571 

We zien in de output opnieuw dezelfde foutmelding: cannot compute exact p-value with ties 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 gedefinieerd bij ties (zie ook cursus).

Hoe interpreteer je de resultaten van de WMW testen?

We vinden significante verschillen in de verdelingen van versgewicht in vergelijking met de controle behandeling. De kans dat het versgewicht voor slaplanten opgegroeid in grond met compost hoger of gelijk is dan (“\(\geq\)”) het versgewicht voor slaplanten opgegroeid in de controle behandeling is gelijk aan 100% (WMW test, aangepaste p-waarde < 0.01). Gelijkaardig, de kans dat het versgewicht voor slaplanten opgegroeid in grond met compost en biochar hoger of gelijk (“\(\geq\)”) is dan het versgewicht voor slaplanten opgegroeid in de controle behandeling is gelijk aan 98.0% (WMW test, aangepaste p-waarde < 0.01).

Bovendien vinden we ook significante verschillen in de verdelingen van versgewicht in vergelijking met de biochar behandeling. De kans dat het versgewicht voor slaplanten opgegroeid in grond met biochar hoger of gelijk (“\(\geq\)”) is dan het versgewicht voor slaplanten opgegroeid in grond met compost is gelijk aan 100% (WMW test, aangepaste p-waarde < 0.01). Ook, de kans dat het versgewicht voor slaplanten opgegroeid in grond met biochar hoger of gelijk (“\(\geq\)”) is dan het versgewicht voor slaplanten opgegroeid in grond met compost en biochar is gelijk aan 100% (WMW test, aangepaste p-waarde < 0.01).

We vinden geen significant verschillende verdelingen in versgewicht van slaplanten tussen grond met biochar en de controle grond en tussen de grond met compost en de grond met compost en biochar.

Conclusie

We zien dat er een extreem significant effect (p << .001) is van het type grond op het versgewicht van slaplanten. Op het 5% significantieniveau stellen we dat de kans dat het versgewicht van een willekeurige slaplant opgegroeid in de ene soort grond hoger is dan of gelijk is (“\(\geq\)”) aan het versgewicht van een willekeurige slaplant uit de andere soorten grond extreem significant verschilt van 50%.

We vinden significante verschillen in de verdelingen van versgewicht in vergelijking met de controle behandeling. De kans dat het versgewicht voor slaplanten opgegroeid in grond met compost hoger of gelijk (“\(\geq\)”) is dan het versgewicht voor slaplanten opgegroeid in de controle behandeling is gelijk aan 100% (WMW test, aangepaste p-waarde < 0.01). Gelijkaardig, de kans dat het versgewicht voor slaplanten opgegroeid in grond met compost en biochar hoger of gelijk (“\(\geq\)”) is dan het versgewicht voor slaplanten opgegroeid in de controle behandeling is gelijk aan 98.0% (WMW test, aangepaste p-waarde < 0.01).

Bovendien vinden we ook significante verschillen in de verdelingen van versgewicht in vergelijking met de biochar behandeling. De kans dat het versgewicht voor slaplanten opgegroeid in grond met biochar hoger of gelijk (“\(\geq\)”) is dan het versgewicht voor slaplanten opgegroeid in grond met compost is gelijk aan 100% (WMW test, aangepaste p-waarde < 0.01). Ook, de kans dat het versgewicht voor slaplanten opgegroeid in grond met biochar hoger of gelijk (“\(\geq\)”) is dan het versgewicht voor slaplanten opgegroeid in grond met compost en biochar is gelijk aan 100% (WMW test, aangepaste p-waarde < 0.01).

We vinden geen significant verschillende kans op een verschillend versgewicht van slaplanten tussen grond met biochar en de controle grond en tussen de grond met compost en de grond met compost en biochar.

We kunnen besluiten dat grond met compost en grond 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 niet besluiten voor bladgroenten in het algemeen.