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 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 de assumpties voldaan wordt. Echter, hoe kleiner de dataset, hoe moeilijker het is om deze assumptie te controleren. Een mogelijke oplossing kan zijn om de Normaliteitsassumptie te checken op de volledige dataset aan de hand van een lineair model, zoals we hieronder doen.

#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 ze een betere keuze zouden maken on gewichten te modelleren. Bijvoorbeeld, gewichten kunnen niet negatief zijn. Een normalie distributie kan negatief worden, maar 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/versgewicht_sla.txt",header=TRUE)

cannot open file '/Users/koenvandenberge/Dropbox/PhD/Onderwijs/Statistiek Biochemie/201819/dropboxStats1819/class6_nonPar/versgewicht_sla.txt': No such file or directoryError in file(file, "rt") : cannot open the connection

Dataverkenning

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

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?

…

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 is dan een willekeurige meetwaarde uit de tweede groep.

Hypothese

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

…

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)

…

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

Post-hoc analyse

Zoals bij de parametrische ANOVA zal je ook bij na 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 om 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?

…

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

library(coin)

…

Kun je nu ook aan de hand van de Mann-Whitney U1 statistiek een puntschatting geven van de Probabilistische index?

…

Hoe interpreteer je de resultaten van de WMW testen?

…

Conclusie

…