--- title: 'PC-practicum 5: one-way anova' output: html_document --- Het is bekend dat de koekoek niet zelf een nest bouwt maar zijn eieren legt in de nesten van andere vogels. Sinds 1892 weet men reeds dat het soort koekoekseieren eigen is aan de locatie waar ze gevonden worden. Een studie in 1940 toonde aan dat de koekoeken elk jaar terugkeren naar hetzelfde grondgebied en eieren leggen in de nesten van welbepaalde "pleegouder"-vogels. Bovendien paren koekoeken enkel binnen hun grondgebied. Op die manier zijn geografische subsoorten ontwikkeld, elk met een dominante pleegouder-soort. Hierdoor kon een specialisatie optreden van de koekoek aan de pleegouder-soort via natuurlijke selectie, zodat de koekoekseieren een hogere kans kregen om geadopteerd te worden door de pleegouder-soort. De dataset koekoek.txt bevat de lengte (variabele `lengte`) van de koekoekseieren (in mm) van willekeurig gekozen geparasiteerde nesten. In totaal bevat de dataset 120 observaties en voor elk ei is aangegeven van welke vogelsoort (variabele `soort`) het nest is. De codering voor soort is als volgt: - `soort=1`: graspieper - `soort=2`: boompieper - `soort=3`: heggemus - `soort=4`: roodborstje - `soort=5`: witte kwikstaart - `soort=6`: winterkoning In deze analyse zullen we nagaan of de pleegouder-soort een invloed heeft op de gemiddelde lengte van de koekoekseieren. Dataset koekoek.txt inlezen ```{r} koekoek<-read.table("koekoek.txt",header=TRUE) soort = koekoek$soort lengte = koekoek$lengte ``` ### Hoeveel observaties zijn er voor elke soort? Frequentietabel en barplot maken voor variabele soort ```{r} table(soort) barplot(table(soort)) ``` ### Data-exploratie ```{r} boxplot(lengte~soort) set.seed(10) stripchart(lengte~soort, vertical = TRUE, method = "jitter", pch = 19, col =c("bisque","coral","darkcyan","steelblue","gold","aquamarine3"), add = TRUE) # eendimensionaal scatterplot ``` ### Welke test kan men uitvoeren om de gemiddelde lengte simultaan te vergelijken tussen alle soorten? In vorige lessen zagen we enkel de two-sample t-test om twee gemiddelden met elkaar vergelijken. We hebben echter ook reeds gezien dat de two-sample t-test een specifieke versie is van een lineair model, namelijk van een lineair model waarbij de covariaat een categorische variabele $X$ is met $2$ levels, i.e. \[ E(Y_i) = \beta_0 + \beta_1 X_i \] Bijvoorbeeld, indien $Y_i$ de lengte van persoon $i$ voorstelt en $X_i$ het geslacht van die persoon waarbij $X_i=0$ indien persoon $i$ een vrouw is, en $X_i=1$ indien niet. In dat geval, stelt $\beta_0$ de gemiddelde lengte voor vrouwen voor, en $\beta_1$ staat voor het verschil in gemiddelde lengte tussen vrouwen en mannen. De gemiddelde lengte voor een man kan men dan bekomen door $E(Y_{male}) = \beta_0 + \beta_1$. Men kan dit ook schrijven als \[ E(Y | female) = \beta_0\] \[ E(Y | male) = \beta_0 + \beta_1\] Dit lineair model kan echter ook makkelijk veralgemeend worden naar factoren met meerdere levels (bvb. het voorbeeld over het verschil in de overlevingstijd bij dojo-/goud-/zebravissen in acquaria met toxische stoffen). We konden dan op basis van het model testen voor verschillen in gemiddelden van een bepaalde level ten opzichte van het intercept, maar het was met de standaard parametrisatie niet steeds mogelijk om alle groepen te vergelijken met elkaar. Voor dit laatste moesten we soms de levels van onze categorische variabele aanpassen, zodat het intercept een andere level voorstelde. Dit is een ruwe manier en kan veel eleganter gedaan worden, door bvb. contrasten van parameters te testen (zie lessen regressie voor Biochemie aan het eind van de cursus). Er bestaat echter ook nog een manier waarbij we **alle levels simultaan kunnen testen**, men zal namelijk testen of de gehele factor variabele een invloed heeft op de respons. In de context van ons voorbeeld, zal men kunnen testen of de pleegouder-soort ??berhaupt een effect heeft op de gemiddelde lengte van koekoekseieren. Zulk een test heet een one-way ANOVA. Men noemt de test 'one-way' omdat het enkel 'main effects' test, met andere woorden het model bevat geen interacties. ### Geef de nul- en alternatieve hypothese voor de test Stel dat $\mu_1$ de gemiddelde lengte van koekoekseieren voor graspiepers (`soort=1`) voorstelt, en idem voor $\mu_2 , \ldots , \mu_6$. De nul- en alternatieve hypothese voor een ANOVA kan men dan voorstellen als $H_0$: $\mu_1=\mu_2=\mu_3=\mu_4=\mu_5=\mu_6$ $H_A$: Voor minstens één $i \ne j$ is $\mu_i \neq \mu_j$ In woorden, zegt de nulhypothese dat de gemiddelde lengte van koekoekseiren onafhankelijk is van de pleegouder-soort: er is geen systematisch verschil in gemiddelde lengte van het ei. De alternatieve hypothese zegt dat de gemiddelde lengte verschilt tussen **minstens twee pleegouder-soorten**. Merk op dat men bij het verwerpen van de nulhypothese **niet weet tussen welke soorten** er een verschil is! ### Fit het model voor de analyse ```{r model} m <- lm(lengte~soort) summary(m) #OPGELET: geen dummy's! Komt omdat soort als continue variabele gezien wordt in R. soort #continu soort <- factor(soort) soort #categorisch m <- lm(lengte~soort) summary(m) ``` De output van het model suggereert dat er inderdaad verschillen lijken te zijn in gemiddelde lengte tussen de pleegoudersoorten. Merk op dat in de standaard output op basis van dit model de p-waarden echter niet aangepast worden voor meervoudig toetsen. ### Ga de assumpties voor een ANOVA na. Merk op dat de assumpties checken niet steeds eenduidig zijn, zeker niet bij een ANOVA test. Zoals beschreven in de cursus, veronderstelt ANOVA een locatie-shift model. Dit wil zeggen dat elke groep een gelijke distributie heeft en er enkel shifts in gemiddelde kunnen optreden. In het bijzonder nemen we de aanname dat elke groep een Normale verdeling zou volgen. Dit impliceert dat - elke groep Normaal verdeeld moet zijn. - elke groep een gelijke variantie moet hebben. Bovendien neemt de test nog aan dat - de groepen onafhankelijk zijn van elkaar. - de gegevens binnen een groep onafhankelijk zijn van elkaar. Deze laatste twee assumpties zijn voldaan; de nesten werden willekeurig gekozen. De eerste twee assumpties kunnen we nagaan indien er niet te veel groepen zijn. Hier hebben we zes groepen en is het checken van assumpties binnen elke groep nog doenbaar. Voor het nagaan van homoscedasticiteit werken we met boxplots: ```{r} boxplot(lengte~soort, xlab="soort", ylab="lengte van het ei (mm)") ``` De data lijken gelijke varianties te hebben, en deze assumptie lijkt alvast niet geschonden. We zullen nu de assumptie van Normale verdeling binnen elke groep nagaan: ```{r} par(mfrow=c(3,2), mar=c(3,2,2,1)) #grafische parameters: verdeel blad in 3 rijen en 2 kolommen en maak margins kleiner. qqnorm(lengte[soort==1]) qqline(lengte[soort==1]) qqnorm(lengte[soort==2]) qqline(lengte[soort==2]) qqnorm(lengte[soort==3]) qqline(lengte[soort==3]) qqnorm(lengte[soort==4]) qqline(lengte[soort==4]) qqnorm(lengte[soort==5]) qqline(lengte[soort==5]) qqnorm(lengte[soort==6]) qqline(lengte[soort==6]) par(mfrow=c(1,1)) #zet terug naar 1 plot per pagina. ``` Elke groep lijkt een Normale verdeling te volgen en deze assumptie is ook voldaan. Met behulp van wat copy-pasting was het nagaan van de assumpties nog doenbaar. Indien men bvb. 20 groepen hoeft te vergelijken is deze manier van werken hinderlijk en mogelijks zelf ambigue. In dat geval kan men ook de residuen van het lineair model checken. Merk op dat men dan checkt voor een normale distributie van alle residuen van de respons variabele ten opzichte van hun groepsgemiddelde, en dus niet voor een normale distributie binnen elke groep... ```{r} par(mfrow=c(2,2)) plot(m) par(mfrow=c(1,1)) ``` De QQ-plot vertoont geen systematische afwijkingen van een Normale distributie. Het nagaan van de assumpties op deze manier vormt zeker geen sluitend bewijs (dat hebben we namelijk zelden) dat de data Normaal verdeeld is binnen elke groep, maar het is een benadering die we kunnen gebruiken in het geval dat - er te veel groepen zijn om de assumpties te checken binnen elke groep; - er te weinig observaties zijn per groep om binnen elke groep de assumpties na te gaan. Merk op dat je in principe de assumptie van homoscedasticiteit ook op basis van de plot linksboven zou kunnen checken: Elke 'kolom' van punten stelt een soort voor (1 soort heeft 1 geschat gemiddelde) en de punten stellen de residuen voor ten opzichte van hun groepsgemiddelde. Men kan deze plot dus ook gebruiken om te kijken of er groepen (soorten) zijn die een verschillende variantie hebben ten opzichte van andere groepen. ### Voer de ANOVA uit We voeren de ANOVA test uit aan de hand van het lineair regressiemodel. In principe testen we dan volgende nulhypothese \[ H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = \beta_5 = 0 \] met de alternatieve hypothese dat minstens twee regressieparameters verschillen van elkaar. Merk op dat deze nulhypothese evenwaardig is aan de nulhypothese hierboven! ```{r} anova(m) ``` Onze p-waarde is bijzonder klein. We besluiten dat we de nulhypothese kunnen verwerpen ($p<<0.001$) en dat de gemiddelde lengte van koekoekseieren verschilt tussen minstens twee van de bestudeerde pleegoudersoorten op het 5% significantieniveau. Aan de hand van dit resultaat weten we echter niet tussen 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. ### Post-hoc analyse De post-hoc analyse bestaat eruit om paarsgewijze testen uit te voeren. Indien men over $k$ groepen beschikt is het totaal aantal paarsgewijze vergelijkingen gelijk aan $k(k-1)/2$. Bij ons is $k=6$ waardoor we $15$ paarsgewijze vergelijkingen zullen uitvoeren. We kunnen echter niet elke test op het 5% significantieniveau uitvoeren vanwege het meervoudig toetsen probleem. Inderdaad, indien men 15 vergelijkingen zou doen, elk op het 5% significantieniveau, dan is de kans dat we minstens één nulhypothese zouden verwerpen terwijl die eigenlijk waar is niet langer gelijk aan ons significantieniveau (5%). In ons geval, zou deze kans gelijk zijn aan ```{R} alpha <- 0.05 nComparisons <- 15 1-(1-alpha)^nComparisons ``` Dus indien we elke test op het 5% significantieniveau zouden uitvoeren hebben we, als alle nulhypotheses waar zouden zijn, een kans van $\sim 54\%$ dat we minstens één nulhypothese verkeerd zouden verwerpen! Om deze kans globaal gezien (dit is, over alle paarsgewijze vergelijkingen) op 5% te houden, kunnen we bijvoorbeeld de Bonferroni correctie uitvoeren. In `R` kunnen we de post-hoc analyse uitvoeren met behulp van het `multcomp` package aan de hand van de `glht` functie. We specifiëren hier in het `linfct` argument dat we *multiple comparisons* (`mcp`) willen uitvoeren waarbij we alle paarsgewijze vergelijkingen voor de `soort` variabele willen testen aan de hand van de `"Tukey"` methode. Je kan de Tukey test aanschouwen als een vorm van de two-sample t-test waarbij de variantie over alle groepen geschat wordt (hierdoor vereist het ook de assumptie van gelijke variantes). Het resultaat van deze test slaan we vervolgens op in het object `mcp`, waarop we een `summary` opvragen van dat object. Het `multcomp` package zorgt ervoor dat deze p-waarden automatisch gecorrigeerd worden voor meervoudig toetsen aan de hand van een efficiënte manier (efficiënter dan de Bonferroni methode). ```{r} library(multcomp) mcp <- glht(m,linfct=mcp(soort="Tukey")) summary(mcp) ``` In de output hiervan zien we de verschillende paarsgewijze vergelijkingen die werden uitvoerd, bijvoorbeeld `2 -1 == 0` duidt erop dat voor dit contrast wordt getest of het gemiddelde voor soort `2` min het gemiddelde voor soort `1` gelijk is aan nul of niet. In de tweede kolom wordt het verschil in gemiddelden weergegeven, met hun standaard error en teststatistiek in de respectievelijk derde en vierde kolom. De laatste kolom geeft aangepaste p-waarden weer op een globaal significantieniveau van 5%. Aan de hand van de aangepaste p-waarden zien we dat de gemiddelde lengte van soort 6 (winterkoning) verschilt van alle andere soorten. De effectgrootte is voor alle soorten negatief, hetgeen erop duidt dat de gemiddelde lengte van koekoekseieren lager is in nesten van winterkoning in vergelijking met andere soorten. Voor de rapportering zullen we ook betrouwbaarheidsintervallen voor elke paarsgewijze vergelijking opvragen. We kunnen deze ook makkelijk grafisch voorstellen aan de hand van de nuttige `plot` functie die zo op een `glht` object kan toegepast worden. ```{r} confint(mcp) plot(mcp) ``` Men kan de resultaten zelf ook nog eens bekijken op basis van de ruwe data: ```{r} boxplot(lengte~soort) set.seed(10) stripchart(lengte~soort, vertical = TRUE, method = "jitter", pch = 19, col =c("bisque","coral","darkcyan","steelblue","gold","aquamarine3"), add = TRUE) ``` ### Besluit We vinden een extreem significante afhankelijkheid tussen de gemiddelde lengte van koekoekseieren en de pleegoudersoort (one-way ANOVA test, $p<<0.001$). Op een globaal 5% significantieniveau vinden we verschillen in gemiddelde lengte van de koekoekseieren tussen verschillende soorten. De gemiddelde lengte van koekoekseieren in nesten van winterkoning is kleiner in vergelijking met deze in nesten van alle andere bestudeerde soorten: graspieper (Tukey test, verschil=-1.16, aangepaste p-waarde = 0.005, 95% BI: [-2.07 ; -0.25]), boompieper (Tukey test, verschil=-1.96, aangepaste p-waarde < 0.001 , 95% BI: [-3.07 ; -0.85]), heggenmus (Tukey test, verschil=-1.99, aangepaste p-waarde < 0.001 , 95% BI: [-3.12 ; -0.86]), roodborstje (Tukey test, verschil=-1.45, aangepaste p-waarde = 0.003 , 95% BI: [-2.53 ; -0.35]) en witte kwikstaart (Tukey test, verschil=-1.77, aangepaste p-waarde < 0.001 , 95% BI: [-2.88 ; -0.66]). We vinden onvoldoende bewijs voor een verschil in gemiddelde lengte van de koekoekseieren tussen de overige soorten.