Mollen zijn gek op tunnels graven. Ondergronds, ergens tussen de wortels, hield de mol bij hoe snel zijn collega’s groeven. 198 mollen had hij langs zijn nieuwsgierige snorharen voelen schuiven, 198 verschillende — elk één werkuur lang, en bij elk noteerde hij hoeveel centimeter er in dat uur bij kwam. Het is werk, hij houdt het bij.
Bij elke mol noteerde hij ook of die in een ondiepe of een diepe gang zat — diepere bodem is compacter, dat scheelt — en of het ochtend of avond was.
“Maakt het uit hoe diep een mol zit,” mompelde hij, “of het tijdstip van de dag, of allebei?”
De mol kon niet zien. Hij was blind, zoals mollen blind zijn. Maar hij kon goed voelen, en hij kon tellen. Per collega één getal. En nu wilde hij weten: zit er verschil tussen die getallen, en zo ja, waardoor?
Wat de mol berekent heet analysis of variance ANOVA, en je stelt elke keer dezelfde drie vragen aan je gegevens:
Strikt genomen valt vrijwel elke statistiek onder ANOVA óf regressie — variatie ontleden in delen die ergens vandaan komen. Het woord wordt pas concreet als we erbij zeggen welke:
Dit hoofdstuk gaat over de eerste twee — eenweg en meerweg (soms ook factorial). Onthoud: zonder dat extra woord blijft “ANOVA” een te ruime container.
Eenweg- en meerweg-ANOVA zijn nog steeds univariate analyses — er is één afhankelijke variabele \(Y\). Het echte multivariate komt pas in thema 5 (MANOVA), waar je meerdere uitkomsten tegelijk vergelijkt. Tot dan: één \(Y\) per analyse.
In thema 1 (MRA) voorspelden we voor een individu via een vergelijking: \[\hat{Y}_i = b_0 + b_1 X_{1i} + b_2 X_{2i} + \dots\] Bij ANOVA voorspellen we óók — alleen niet voor een individu via continue voorspellers, maar voor een groep via het groepsgemiddelde:
\[\hat{Y}_{g} = \bar{Y}_g\]
Het groepsgemiddelde \(\bar{Y}_g\) is je beste voorspelling voor iedereen in groep \(g\). In de regressie-taal (sommige boeken doen dit expliciet via dummy- of effect-codering) heet dit ook wel de estimated marginal mean EMM: de modelvoorspelling per cel, gemiddeld over de andere factoren als die er zijn.
Onthoud: ANOVA is geen “toetsings-techniek” en MRA geen “voorspellings-techniek” — beide voorspellen, ze gebruiken alleen verschillende voorspelvormen (\(\hat{Y}\) versus \(\bar{Y}_g\)).
Dit hoofdstuk volgt dezelfde drie-laagse structuur als thema 1:
Y, factor1, mijn_data.# hekjes en concrete namen.Plus de antwoorden in zacht groen (collapse), vragen in indigo, leesregels in mauve, alarm in oranje, don’ts in rood. T-vragen vlechten zich tussen de hoofdvragen door om je rekenintuïtie te oefenen.
Wat de mol uit zijn keukenkastje pakt
| Gereedschap | Waar het voor is | Pakket |
|---|---|---|
lm() |
Lineair model schatten — ook voor ANOVA | base |
Anova(model, type = 3) |
ANOVA-tabel met type III sums of squares | car |
aov() |
Alternatieve ANOVA-fit, handig voor TukeyHSD() |
base |
aggregate() |
Groepsgemiddelden uitrekenen | base |
table() |
Aantallen per cel | base |
leveneTest() |
Toets voor gelijke varianties | car |
TukeyHSD() |
Post-hoc paarsgewijze vergelijkingen | base |
etaSquared() |
Effectgrootte \(\eta^2\) en \(\eta^2_p\) | lsr |
interaction.plot() |
Snelle visualisatie van een interactie | base |
emmeans() |
Estimated marginal means en EMM-uitwerking per niveau | emmeans |
anova(), Anova() en summary() — drie functies, drie tabellen
Bij een gefit lm-object zijn er drie verschillende manieren om naar de modeloutput te kijken, en ze geven niet dezelfde tabel:
| Functie | Pakket | Wat geeft het? | Wanneer gebruik je ’m? |
|---|---|---|---|
summary(m) |
base | Coëfficiënten-tabel + \(F\)-toets voor het hele model + \(R^2\). Geen \(SS\). | Snelle blik op coëfficiënten + modelfit. Voor regressie. |
anova(m) |
base | ANOVA-tabel met Type I SS (sequentieel). Volgorde van termen telt. | Alleen bij gebalanceerd design — of als je expliciet de volgorde-volgorde-test wilt. Bij regressie zelden. |
Anova(m, type = 3) |
car |
ANOVA-tabel met Type III SS (elk effect uniek). Volgorde maakt niet uit. | Standaardkeuze voor ANOVA met ongebalanceerde data — wat in psychologie de norm is. |
Praktijk-vuistregel voor dit hoofdstuk: gebruik altijd car::Anova(., type = 3) voor ANOVA-tabellen. Voor de coëfficiënten (\(b\)’s in een regressie-parameterisering, of EMM in factor-parameterisering) blijf je naar summary(m) kijken. anova(m) (kleine letter) gebruik je vrijwel nooit in dit boek — pas op dat je ’m niet per ongeluk typt en denkt dat het werkt.
Bijvalletje: Anova() (hoofdletter, car) eist options(contrasts = c("contr.sum", "contr.poly")) voor correcte Type III hoofdeffecten bij ongebalanceerde data. Die regel staat al boven in het script — laat ’m staan.
Bij een gebalanceerd design (gelijke aantallen per cel) maakt het niet uit: alle types geven hetzelfde antwoord. Maar zodra het ongebalanceerd is — wat in psychologisch onderzoek vrijwel altijd zo is — wijken Type I, II en III sums of squares uit elkaar. Type III geeft elke voorspeller zijn unieke bijdrage onafhankelijk van de invoervolgorde, en is daarmee de standaardkeuze voor de meeste APA-rapportages. Daarom de regel options(contrasts = c("contr.sum", "contr.poly")) bovenaan: die zorgt dat Anova(., type = 3) de juiste antwoorden geeft.
In dit werkboek zie je telkens dezelfde notatie:
aggregate() of mean() je geeft. Een schatter van \(\mu_g\).In hypothese-toetsen praten we altijd over de populatie — dus daar staat altijd \(\mu_g\). In rapportages over wat je gevonden hebt, schrijf je \(\bar{Y}_g\) of \(M\).
De drempels in dit werkboek zijn breed gangbare conventies — maar niet universeel:
oneway.test(..., var.equal = FALSE)) als alternatief;Verschillende vakgroepen, docenten en handboeken kiezen iets andere getallen (\(\eta^2 \approx .05\) als groot in heel kleine studies; strengere Levene-kritiek alleen bij hoog ongebalanceerde designs). Als je voor een tentamen of opdracht werkt, check altijd in je eigen college-sheets welke drempels je docent of vakgroep hanteert. Wij volgen hier consistent één set, maar dat is een keuze, geen natuurwet.
Code is standaard open als je hem moet kunnen typen op het R-practical-tentamen: lm(), aov(), Anova(., type = 3), leveneTest(), etaSquared(), TukeyHSD(), interaction.plot(), emmeans(), contrast(), aggregate(), table(). Code die alleen ter illustratie dient — verkenning, plot-decoratie — staat ingeklapt met een knopje “Toon code”. Klap hem open als je nieuwsgierig bent; voor het tentamen hoef je hem niet te reproduceren.
Volgorde-regel. Alle library()-aanroepen bovenaan je script — vóór elk gebruik van een functie. Schuif je library(car) per ongeluk onder een Anova()- of leveneTest()-regel, dan klaagt R “could not find function”. Vervelend, en altijd dezelfde oplossing: pakket boven gebruik.
opmerkingen_jij.R ernaast
Tijdens het werken loop je per definitie tegen dingen aan: een commando dat anders werkt dan je verwacht, een zin in het werkboek die niet helemaal klopt, een vraag waarvan je achteraf denkt “die had hier eigenlijk moeten staan.” Open daarvoor een leeg R-script naast je werkbestand, noem het bijvoorbeeld opmerkingen_jij.R, en typ daar al je shit in — letterlijke citaten, regelnummers, tijdstempels, je eigen vraag. Het script wordt niet gerund; het is je kantlijn. Lever het in bij Ben (per mail of in de gedeelde map) en de volgende werkboek-versie wordt iets minder ruk.
De mol bewaart zijn aantekeningen netjes
Open de meegestuurde projectmap (02_variantieanalyse/) en dubbelklik op 02_variantieanalyse.Rproj. Daarmee staat je werkomgeving klaar: paden kloppen, RStudio kent de juiste werkmap. De data ligt al in data/, de helper-functies in functions/. Niets te installeren als de vier pakketten (car, lsr, emmeans, tidyverse) op je computer staan; anders één keer:
install.packages(c("car", "lsr", "emmeans", "tidyverse"))Algemene vorm. Een .RData-bestand laden zet het bewaarde object direct in je sessie:
load("data/jouw_data.RData")
str(jouw_data)Voor onze dieren. Het bestand bevat een data frame met de naam mol_werkdag.
load("data/mol_werkdag.RData")
str(mol_werkdag)'data.frame': 198 obs. of 3 variables:
$ graaftempo : int 40 24 22 32 22 26 38 30 35 34 ...
$ gang_diepte: Factor w/ 2 levels "ondiep","diep": 1 1 1 1 1 1 1 1 1 1 ...
$ dagdeel : Factor w/ 2 levels "ochtend","avond": 1 1 1 1 1 1 1 1 1 1 ...Per geobserveerde mol drie waarden: gang_diepte (ondiep/diep), dagdeel (ochtend/avond), en graaftempo (in cm/uur). gang_diepte en dagdeel zijn al factor-variabelen — R weet dat het categorieën zijn, geen getallen.
sapply(mol_werkdag, class) graaftempo gang_diepte dagdeel
"integer" "factor" "factor" contr.sum voor Type-III SS
Onze setup zet options(contrasts = c("contr.sum", "contr.poly")). Dit is een technische voorwaarde voor car::Anova(., type = 3) om correcte hoofdeffecten te berekenen bij ongebalanceerde designs. Met R’s default contr.treatment zou Anova(type = 3) foutieve hoofdeffect-getallen geven. Voor de inhoudelijke interpretatie maakt het niets uit — F-toetsen, p-waarden, \(\eta^2/\eta^2_p\) en EMM zijn identiek tussen beide; alleen de coëfficiënten in summary(lm) worden anders gecodeerd. Onthoud: zet deze regel boven aan je script bij Type-III ANOVA.
De mol vergelijkt twee groepen
“Eerst,” zei de mol, “kijk ik alleen naar de gang-diepte. Maakt het uit of ik in een ondiepe of in een diepe gang zit?”
Het stappenplan van thema 1, nu uitgebreid met stap 5: de techniek-keuze.
Wie of wat wordt er gemeten? \(198\) verschillende mollen, elk één werkuur lang geobserveerd.
Wat wordt er gemeten? Drie variabelen per mol: gang_diepte, dagdeel, en graaftempo (cm/uur).
Onafhankelijk of afhankelijk? Onderzoeksvraag (deze deelvraag): voorspelt gang_diepte het graaftempo? Dan is graaftempo afhankelijk (\(Y\)), gang_diepte onafhankelijk (\(X\)).
Meetniveau van elke variabele? graaftempo is interval (cm/uur). gang_diepte is nominaal (twee categorieën). dagdeel is nominaal (twee categorieën).
Welke techniek? Doorloop de beslisboom:
flowchart TD
A[Hoeveel afhankelijke<br/>variabelen?] -->|één| B[Y meetniveau?]
A -->|meerdere| Z[MANOVA<br/><i>thema 5</i>]
B -->|interval| C[X-en meetniveau?]
B -->|binair| D[Logistic regression<br/><i>thema 4</i>]
C -->|alleen interval| E[Multiple regression<br/><i>thema 1</i>]
C -->|alleen nominaal| F[ANOVA<br/><i>dit thema</i>]
C -->|gemengd| G[ANCOVA<br/><i>thema 3</i>]
style F fill:#faf3e2,stroke:#c9a05a,stroke-width:2px
Dezelfde keuze, anders genoteerd. Lees als predictor(en) (meetniveau) → DV(s) (meetniveau):
| # | Predictor(en) | DV(s) | Techniek |
|---|---|---|---|
| dich |X| | |Y| int | \(t\)-toets | |
| nom 3+lvl |X| | |Y| int | éénweg ANOVA | |
| 1 | int |X X …| | |Y| int | MRA — thema 1 |
| 2 | nom |X X| | |Y| int | factorial ANOVA — thema 2 |
| 3 | nom + int |X C| | |Y| int | ANCOVA — thema 3 |
| 4 | bo(in) |X| | |Y| bin | LRA — thema 4 |
| 5 | nom |X| | |Y₁ Y₂ …| int | MANOVA — thema 5 |
| 6 | within-subject momenten | |Y₁ Y₂ Y₃ Y₄| int | RMA — thema 6 |
| 7 | indirect via \(M\) | \(X \to M \to Y\) pad | Mediation — thema 7 |
\(Y\) is interval (graaftempo), \(X\) is nominaal (gang_diepte). Antwoord: ANOVA. Bij twee categorieën is dat technisch hetzelfde als een \(t\)-toets met gelijke variantie, maar we doen het via lm() zodat we straks naadloos kunnen uitbreiden naar twee factoren.
Drie miniatuur-onderzoeken — pen-en-papier
Voor elk onderzoek: bepaal aantal afhankelijke variabelen, het meetniveau van \(X\) en \(Y\), en kies dan de techniek.
a) De bever vergelijkt drie burcht-types — modder, hout, gemengd — op damlengte (m). Bij \(50\) bevers één meting per dier. Welke techniek?
b) Een klinisch psycholoog vergelijkt drie therapie-types (CGT / EMDR / wachtlijst) en twee leeftijd-categorieën (jong-volwassen / ouder dan \(40\)) op symptoom-vermindering (schaal \(0\)–\(100\)). Bij \(120\) cliënten. Welke techniek?
c) Een HR-team bekijkt drie afdelingen (sales / R&D / support) op zowel werktevredenheid (schaal \(1\)–\(10\)) als ervaren werkdruk (schaal \(1\)–\(10\)). Bij \(90\) medewerkers, beide tegelijk. Welke techniek?
a) Eén afhankelijke (\(Y =\) damlengte, INT) en één onafhankelijke met drie groepen (NOM). \(\Rightarrow\) eenweg-ANOVA.
b) Eén afhankelijke (\(Y =\) symptoom-vermindering, INT) en twee onafhankelijke factoren (NOM × NOM). \(\Rightarrow\) tweeweg-ANOVA (meerweg / factorial). Je krijgt dan twee hoofdeffecten + hun interactie.
c) Twee afhankelijke variabelen tegelijk (werktevredenheid én werkdruk). Geen ANOVA meer: \(\Rightarrow\) MANOVA (thema 5).
Voordat de mol gaat toetsen, kijkt hij naar wat er feitelijk in de groepen zit.
Algemene vorm. Tellingen per groep en groepsgemiddelden van de uitkomstvariabele:
table(mijn_data$factor) # aantallen per groep
aggregate(Y ~ factor, data = mijn_data, FUN = mean) # groepsgemiddeldenVoor onze dieren.
# Aantal mollen per gang_diepte-groep.
table(mol_werkdag$gang_diepte)
ondiep diep
99 99 # Gemiddeld graaftempo per gang_diepte-groep.
aggregate(graaftempo ~ gang_diepte, data = mol_werkdag, FUN = mean) gang_diepte graaftempo
1 ondiep 29.72727
2 diep 20.37374a) Hoeveel mollen zitten er in elke gang_diepte-groep? Is het design gebalanceerd?
b) Wat zijn de groepsgemiddelden? Welke groep graaft sneller, en hoe groot is het ruwe verschil?
In gewone woorden.
In de ondiep-groep zitten 99 mollen, in de diep-groep 99. Op de marginale telling exact gelijk; binnen de cellen — daar zien we straks dat de cel-aantallen wel iets uit elkaar lopen.
Het gemiddelde tempo in de ondiepe gang (\(\bar{Y}_{\text{ondiep}} = 29.73\)) ligt hoger dan in de diepe (\(\bar{Y}_{\text{diep}} = 20.37\)). Een ruw verschil van \(9.36\) cm/uur.
APA-stijl.
De marginale aantallen waren gelijk (\(n_{\text{ondiep}} = 99\); \(n_{\text{diep}} = 99\)); de cel-aantallen verschilden licht (zie 2.2). In ondiepe gangen groef de mol gemiddeld sneller (\(M = 29.73\) cm/uur) dan in diepe gangen (\(M = 20.37\) cm/uur).
De \(F\)-toets in een eenweg-ANOVA is geen mysterie: het is een verhouding tussen twee varianties.
\[ F = \frac{\text{MS}_b}{\text{MS}_w} = \frac{\text{SS}_b / df_b}{\text{SS}_w / df_w} \]
Bij de mol: \(k = 2\) groepen, \(N = 198\) mollen. Stel je leest in een ANOVA-tabel af: \(\text{SS}_b = 4331\) en \(\text{SS}_w = 19341\).
a) Hoeveel zijn \(df_b\) en \(df_w\)?
b) Bereken \(\text{MS}_b\) en \(\text{MS}_w\).
c) Bereken \(F\). Dit getal moet je straks in de output van Anova(...) kunnen vinden.
a) \(df_b = k - 1 = 1\). \(df_w = N - k = 198 - 2 = 196\).
b) \(\text{MS}_b = 4331 / 1 = 4331\). \(\text{MS}_w = 19341 / 196 = 98.68\).
c) \(F = 4331 / 98.68 = 43.89\). Straks zie je in de Type-III-output: \(F(1, 196) = 43.89\), \(p < .001\). Bingo.
Inzicht. \(F\) vraagt: “is de spreiding tussen de groepsgemiddelden groot ten opzichte van de spreiding binnen de groepen?” Als groepen weinig verschillen ten opzichte van wat individuen al onderling verschillen, valt \(F \approx 1\).
ANOVA leunt op drie aannames:
leveneTest()-functie.Algemene vorm.
# Levene's toets voor gelijke varianties (de belangrijkste check)
leveneTest(Y ~ factor, data = mijn_data)
# Visueel: residual plot van het lm-object
plot(lm(Y ~ factor, data = mijn_data), which = 1)Voor onze dieren.
# Levene: zijn de varianties tussen de twee groepen ongeveer gelijk?
leveneTest(graaftempo ~ gang_diepte, data = mol_werkdag)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.0825 0.7742
196 m_diepte <- lm(graaftempo ~ gang_diepte, data = mol_werkdag)
# `sub.caption = ""` haalt de auto-ondertitel met de modelformule weg.
plot(m_diepte, which = 1, sub.caption = "",
caption = "Residuen tegen voorspelde waarden")
c) Is de aanname van gelijke varianties geschonden? Rapporteer Levene.
d) Wat zie je in de residual plot — een rustige horizontale band, of iets verdachts?
In gewone woorden.
De Levene-toets is niet significant (\(F(1, 196) = 0.08\), \(p = .77\)): we kunnen niet aantonen dat de varianties verschillen, dus de aanname van gelijke varianties houdt stand.
De residual plot toont twee verticale puntenwolken (één per groep), met vergelijkbare spreiding en zonder uitwaaiering. Geen rode vlaggen.
APA-stijl.
De aanname van gelijke varianties was houdbaar, \(F(1, 196) = 0.08\), \(p = .77\). Visuele inspectie van de residual plot wees niet op afwijkingen van homoscedasticiteit.
De \(F\)-toets verdraagt heel wat. Twee getallen om te onthouden, één regel voor de richting:
Twee-weg richtingsregel (factorial). Liggen de grootste varianties in de kleinste groepen, dan is \(F\) te liberaal (\(p\)-waarden te klein); in de grootste groepen, dan is \(F\) te conservatief (\(p\)-waarden te groot). Conclusie houdt meestal stand, behalve als \(p\) tussen \(.05\) en \(.01\) ligt (liberaal-context) of tussen \(.05\) en \(.10\) (conservatief-context) — dan met voorzichtigheid rapporteren.
In de woorden van het Exercise Book (p. 42): “When your groups are large enough: don’t worry!”
Bij ANOVA volstaat de robuustheid van de \(F\)-toets meestal voor lichte schendingen, maar kijk uit als:
onewaytests::bf.test()) of een Welch-ANOVA (oneway.test(..., var.equal = FALSE)).Vuistregel: noem altijd dat je hebt gecheckt, en wat je beslist hebt te doen.
Levene’s toets voor variantie-homogeniteit komt in onze workflow standaard mee — al moet hier opgemerkt: de cursus-lecture geeft expliciet aan “a statistical test for variance equality exists (Levene’s test). However, we do not use it (at least not in MVDA)”, terwijl de antwoord-PDF Levene wel rapporteert. We volgen hier de exercise-praktijk (Levene rapporteren), maar de robuustheid van \(F\) blijft het hoofdcriterium (\(n \geq 15\) + \(n_{\max}/n_{\min} \leq 1.5\) — zie callout hierboven). Bij beide schendingen tegelijk: Brown-Forsythe (one-way) of richtingsregel (factorial).
Algemene vorm.
mijn_model <- lm(Y ~ factor, data = mijn_data)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
# Defensief — voor het geval dat een eerder blok de instelling heeft veranderd.
options(contrasts = c("contr.sum", "contr.poly"))
Anova(mijn_model, type = 3) # ANOVA-tabel met F-toets
etaSquared(mijn_model) # eta^2 als effectgrootteVoor onze dieren.
# Gebruikt 'm_diepte' uit de assumptie-stap hierboven.
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
options(contrasts = c("contr.sum", "contr.poly"))
Anova(m_diepte, type = 3)Anova Table (Type III tests)
Response: graaftempo
Sum Sq Df F value Pr(>F)
(Intercept) 124251 1 1259.156 < 2.2e-16 ***
gang_diepte 4331 1 43.887 0.0000000003275 ***
Residuals 19341 196
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Effectgrootte: eta-squared.
etaSquared(m_diepte) eta.sq eta.sq.part
gang_diepte 0.1829494 0.1829494\(\eta^2\) is de proportie totale variantie van \(Y\) die door de factor wordt verklaard. Vuistwaarden:
Houd rekening met je vakgebied: in experimentele psychologie zijn deze vuistmaten realistisch, in toegepaste contexten ligt de lat soms lager.
De \(F\)-toets in Anova(., type = 3) toetst niet iets over de mollen die je hebt gemeten — die kun je gewoon aflezen. Hij toetst iets over de mollen die je niet hebt gemeten: de hele populatie waaruit deze steekproef komt. Schrijf de nulhypothese daarom altijd in populatie-taal:
\[H_0: \mu_{\text{ondiep}} = \mu_{\text{diep}}\]
— “in de populatie van alle mollen verschillen ondiepe en diepe gangen niet in graaftempo”. Niet \(\bar{Y}_{\text{ondiep}} = \bar{Y}_{\text{diep}}\) — dat is over je 198 mollen, en die zijn meetbaar ongelijk (\(29.73\) versus \(20.37\)). Eerstejaars-val: \(H_0\) schrijven met steekproefgemiddeldes erin. Steekproef = data, populatie = uitspraak.
Twee veilige formuleringen voor het gevonden effect:
Wat je niet mag zeggen zonder experimentele opzet: “diepe gangen vertragen de mol”, “door de diepte daalt het tempo”, “ondiepe gangen veroorzaken sneller graven”. Dat zijn causale werkwoorden — die mag je pas gebruiken als je studie-opzet ze rechtvaardigt (manipulatie, randomisatie, voor-en-na-meting). De mol-data is een momentopname; geen ingreep, geen tijd-as. Mol-keuze voor diepe of ondiepe gang kan zelf met andere kenmerken samenhangen die het tempo verklaren — een derde-variabele-val.
Anova(m_diepte, type = 3) print Sum Sq, Df, F en Pr(>F). Het MS-getal (mean square) staat er niet — maar dat is gewoon \(SS/df\). Met die identiteiten vul je de hele sudoku-tabel terug:
| Cel | Formule | Voor het mol-model |
|---|---|---|
| \(SS_b\) (between) | uit Anova() regel “gang_diepte” |
\(4331\) |
| \(SS_w\) (within / Residuals) | uit Anova() regel “Residuals” |
\(19341\) |
| \(SS_T\) | \(SS_b + SS_w\) | \(23672\) |
| \(df_b\) | \(k - 1\) (groepen − 1) | \(2 - 1 = 1\) |
| \(df_w\) | \(N - k\) | \(198 - 2 = 196\) |
| \(df_T\) | \(N - 1\) | \(197\) |
| \(MS_b\) | \(SS_b / df_b\) | \(4331 / 1 = 4331\) |
| \(MS_w\) | \(SS_w / df_w\) | \(19341 / 196 = 98.68\) |
| \(F\) | \(MS_b / MS_w\) | \(4331 / 98.68 \approx 43.89\) ✓ |
| \(\eta^2\) | \(SS_b / SS_T\) | \(4331 / 23672 = .183\) ✓ |
Tentamenvraag-format: je krijgt een halve tabel en moet de rest invullen. Wie de identiteiten kent (\(SS_T = SS_b + SS_w\), \(MS = SS/df\), \(F = MS_b/MS_w\), \(\eta^2 = SS_b/SS_T\)) komt overal uit zonder iets uit het hoofd te leren. Sudoku.
e) Is het effect van gang_diepte statistisch significant? Rapporteer \(F\), \(df\), \(p\) en \(\eta^2\).
f) Is het effect inhoudelijk groot, middelgroot of klein?
g) Wat is je voorlopige conclusie over de mol’s vraag?
In gewone woorden.
Het effect van gang_diepte is duidelijk aanwezig. De \(F\)-toets geeft \(F(1, 196) = 43.89\) met \(p < .001\) — ruim onder de drempel van \(.05\), dus statistisch significant.
De effectgrootte is \(\eta^2 = .183\): groot (\(\geq .14\)).
Voorlopige conclusie: gang_diepte maakt veel uit. In een ondiepe gang graaft de mol gemiddeld bijna tien centimeter per uur sneller dan in een diepe gang, en dat past bij wat je biologisch zou verwachten — diepere bodem is compacter. Maar pas op: dit is een eenweg-analyse die dagdeel negeert. Wat als dagdeel ertoe doet, en het diepte-effect verschilt tussen ochtend en avond? Dat zoeken we uit in 2.2.
APA-stijl.
Een eenweg-ANOVA toonde een significant en groot effect van gang_diepte op het graaftempo, \(F(1, 196) = 43.89\), \(p < .001\), \(\eta^2 = .183\). In ondiepe gangen (\(M = 29.73\) cm/uur) groef de mol gemiddeld \(9.36\) cm/uur sneller dan in diepe gangen (\(M = 20.37\) cm/uur).
Bij ANOVA is de modelvoorspelling voor elk dier het gemiddelde van zijn groep:
\[\hat{Y}_i = \bar{Y}_{g(i)}\]
Voor de mol-data zijn de groepsgemiddelden \(\bar{Y}_{\text{ondiep}} = 29.73\) en \(\bar{Y}_{\text{diep}} = 20.37\) cm/uur.
a) Een mol heeft gang_diepte = ondiep. Wat is \(\hat{Y}\) voor die mol?
b) Een andere mol heeft een werkelijk graaftempo van \(35\) cm/uur in een diepe gang. Wat is zijn residu \(e_i = Y_i - \hat{Y}_i\)?
c) In de regressievertaling van deze ANOVA met dummy-codering (ondiep = 0, diep = 1) zou je krijgen: \[\hat{Y}_i = b_0 + b_1 \cdot \text{diep}_i\] Wat zijn \(b_0\) en \(b_1\) in deze parameterisering?
a) \(\hat{Y} = \bar{Y}_{\text{ondiep}} = 29.73\) cm/uur. Dat is je beste voorspelling voor élke mol in die groep — niets meer, niets minder.
b) \(e = 35 - 20.37 = 14.63\) cm/uur.
c) \(b_0 = \bar{Y}_{\text{ondiep}} = 29.73\) (intercept = referentie-groep). \(b_1 = \bar{Y}_{\text{diep}} - \bar{Y}_{\text{ondiep}} = 20.37 - 29.73 = -9.36\) (verschil tussen groepen — diep ligt \(9.36\) cm/uur lager). Dat \(b_1\) — het verschil tussen twee groepsgemiddelden — is precies waar de \(F\)-toets met één voorspeller op toetst. Met contr.sum (sum-to-zero, voor Type III) krijg je andere \(b\)’s, maar dezelfde voorspellingen \(\hat{Y}_g\) en dezelfde \(F\).
Inzicht. Eenweg-ANOVA is een lineair regressiemodel met één gecodeerde factor. De \(F\)-toets in Anova(., type = 3) is in dit geval gewoon de gekwadrateerde \(t\)-toets op \(b_1\).
De mol voegt een tweede vraag toe
“Maar wat,” vroeg de mol zich af, “als het diepte-effect ’s avonds anders is dan ’s ochtends? Dan is gang_diepte alleen geen eerlijk verhaal als ik dagdeel negeer.”
Dat is wat factorial ANOVA doet: twee (of meer) nominale factoren tegelijk, plus hun interactie interaction. Een interactie betekent: het effect van factor A is niet hetzelfde voor alle niveaus van factor B. Het effect van A op zichzelf — gemiddeld over de niveaus van B — heet het hoofdeffect main effect van A.
In dit hoofdstuk volgen we een conventie: factor A is de rij-factor, factor B is de kolom-factor. Bij de mol: \(A = \text{gang\_diepte}\) (rijen) en \(B = \text{dagdeel}\) (kolommen). Die zelfde A=rij, B=kolom-afspraak gebruiken we straks bij de olifant en bij de dassen.
Een factorial design heeft minstens twee factoren met elk twee of meer niveaus. De kleinst mogelijke factoriële opzet is \(2 \times 2\): vier cellen, twee hoofdeffecten en één interactie. Daarmee zie je exact hetzelfde patroon dat in een \(2 \times 3\) of \(3 \times 4\) design speelt — alleen overzichtelijker.
Daarom de mol-data eerst in \(2 \times 2\) (gang_diepte × dagdeel): hier ontstaat het hele begrip interactie in vier zinnen. Bij de olifant verderop kruipen we naar \(2 \times 3\), en straks bij de dassen een design met meer niveaus. Vuistregel voor jezelf: als je een nieuwe techniek leert, begin met het kleinste voorbeeld waarin het patroon nog werkt. Bij multipele regressie is dat twee voorspellers (niet drie), bij ANOVA is dat \(2 \times 2\). Pas als je daar geen moeite mee hebt, schaal je op.
# Aantallen per cel (= combinatie van gang_diepte x dagdeel)
table(mol_werkdag$gang_diepte, mol_werkdag$dagdeel)
ochtend avond
ondiep 52 47
diep 49 50# Gemiddelden per cel
aggregate(graaftempo ~ gang_diepte * dagdeel, data = mol_werkdag, FUN = mean) gang_diepte dagdeel graaftempo
1 ondiep ochtend 25.88462
2 diep ochtend 19.42857
3 ondiep avond 33.97872
4 diep avond 21.30000a) Welke combinatie van gang_diepte en dagdeel scoort het hoogst? Welke het laagst?
b) Is het verschil tussen ondiep en diep even groot in de ochtend als in de avond? Wat zegt dat over een mogelijke interactie?
In gewone woorden.
De combinatie ondiep + avond scoort het hoogst (\(M = 33.98\) cm/uur). De diep + ochtend-cel scoort het laagst (\(M = 19.43\)).
In de ochtend is het verschil tussen ondiep (\(25.88\)) en diep (\(19.43\)) modest (\(\Delta = 6.46\)). In de avond is dat verschil veel groter (\(33.98\) vs \(21.30\), \(\Delta = 12.68\)). Dat is een teken van interactie: het effect van gang_diepte is groter in de avond dan in de ochtend.
APA-stijl.
Inspectie van de celgemiddelden suggereerde een interactie tussen gang_diepte en dagdeel: in de avond lag het tempo in ondiepe gangen (\(M = 33.98\) cm/uur) ruim \(12\) cm/uur hoger dan in diepe gangen (\(M = 21.30\)), terwijl dat verschil in de ochtend kleiner was (\(M_{\text{ondiep}} = 25.88\); \(M_{\text{diep}} = 19.43\); \(\Delta \approx 6.5\)).
Hieronder vier mini-tabellen met cel-gemiddelden voor een \(2 \times 2\) design (factor A heeft twee niveaus \(A_1, A_2\); factor B heeft twee niveaus \(B_1, B_2\)). Voor elk geval: wat zie je — alleen hoofdeffect \(A\), alleen hoofdeffect \(B\), alleen interactie, of zowel hoofd- als interactie-effect?
| Geval | \(A_1 B_1\) | \(A_1 B_2\) | \(A_2 B_1\) | \(A_2 B_2\) |
|---|---|---|---|---|
| (i) | 4 | 6 | 4 | 6 |
| (ii) | 4 | 6 | 6 | 4 |
| (iii) | 4 | 4 | 6 | 6 |
| (iv) | 4 | 5 | 6 | 9 |
a) Schets voor elk geval de interactieplot in je hoofd: x-as \(A\), twee lijnen voor \(B\). Welk geval heeft parallelle lijnen, en welke snijdende of uiteenwijkende lijnen?
b) Welk geval is een kruis-interactie (lijnen kruisen)? En een spreiding-interactie (lijnen waaieren uit, kruisen niet)?
Vuistregel. Parallelle lijnen \(\Rightarrow\) geen interactie. Niet-parallelle lijnen \(\Rightarrow\) interactie. Hoeveel niet-parallel zegt iets over de grootte van het interactie-effect; of de toets significant is hangt ook af van \(N\).
Algemene vorm.
mijn_model <- lm(Y ~ factor1 * factor2, data = mijn_data)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
options(contrasts = c("contr.sum", "contr.poly"))
Anova(mijn_model, type = 3)
etaSquared(mijn_model)
interaction.plot(x.factor = mijn_data$factor1,
trace.factor = mijn_data$factor2,
response = mijn_data$Y)De * in de formule (factor1 * factor2) betekent: zowel beide hoofdeffecten als hun interactie meenemen. Schrijf je factor1 + factor2, dan krijg je alleen de hoofdeffecten — geen interactie.
Voor onze dieren.
# Levene voor de vier-cellen-versie.
leveneTest(graaftempo ~ gang_diepte * dagdeel, data = mol_werkdag)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 0.696 0.5555
194 m_factorial <- lm(graaftempo ~ gang_diepte * dagdeel, data = mol_werkdag)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
options(contrasts = c("contr.sum", "contr.poly"))
Anova(m_factorial, type = 3)Anova Table (Type III tests)
Response: graaftempo
Sum Sq Df F value Pr(>F)
(Intercept) 125053 1 1375.5504 < 2.2e-16 ***
gang_diepte 4525 1 49.7733 0.00000000002974 ***
dagdeel 1227 1 13.5006 0.0003085 ***
gang_diepte:dagdeel 479 1 5.2639 0.0228433 *
Residuals 17637 194
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1etaSquared(m_factorial) eta.sq eta.sq.part
gang_diepte 0.18872565 0.20210717
dagdeel 0.05177013 0.06496975
gang_diepte:dagdeel 0.02021610 0.02641659# Interactieplot: x-as is dagdeel (kolom-factor B),
# lijnen voor gang_diepte (rij-factor A).
# Tijdelijke factoren met hoofdletter-niveaus voor leesbare as- en
# legend-labels — onderliggende data blijft ongewijzigd.
.dagdeel_lbl <- factor(tools::toTitleCase(as.character(mol_werkdag$dagdeel)),
levels = tools::toTitleCase(levels(mol_werkdag$dagdeel)))
.diepte_lbl <- factor(tools::toTitleCase(as.character(mol_werkdag$gang_diepte)),
levels = tools::toTitleCase(levels(mol_werkdag$gang_diepte)))
interaction.plot(
x.factor = .dagdeel_lbl,
trace.factor = .diepte_lbl,
response = mol_werkdag$graaftempo,
ylab = "Gemiddeld graaftempo (cm/uur)",
xlab = "Dagdeel",
trace.label = "Gangdiepte"
)
Bij een \(a \times b\) factorial design (factor \(A\) met \(a\) niveaus, factor \(B\) met \(b\) niveaus, \(N\) observaties totaal) decomponeer je de totale SS in vier stukken: hoofdeffect \(A\), hoofdeffect \(B\), interactie \(A \times B\), en error. De vrijheidsgraden lopen mee:
| Bron | \(df\) | \(SS\) | \(MS\) | \(F\) |
|---|---|---|---|---|
| Hoofdeffect \(A\) | \(a - 1\) | \(SS_A\) | \(SS_A / df_A\) | \(MS_A / MS_E\) |
| Hoofdeffect \(B\) | \(b - 1\) | \(SS_B\) | \(SS_B / df_B\) | \(MS_B / MS_E\) |
| Interactie \(A \times B\) | \((a-1)(b-1)\) | \(SS_{AB}\) | \(SS_{AB} / df_{AB}\) | \(MS_{AB} / MS_E\) |
| Error (within) | \(N - a \cdot b\) | \(SS_E\) | \(SS_E / df_E\) | — |
| Totaal | \(N - 1\) | \(SS_T\) | — | — |
Identiteiten — altijd geldig:
Voor de mol (\(a = 2\) gang_diepte, \(b = 2\) dagdeel, \(N = 198\)):
| Bron | \(df\) |
|---|---|
| gang_diepte (\(A\)) | \(1\) |
| dagdeel (\(B\)) | \(1\) |
| gang_diepte × dagdeel (\(A \times B\)) | \(1\) |
| Error | \(194\) |
| Totaal | \(197\) |
Check: \(1 + 1 + 1 + 194 = 197\). ✓
Bij een \(2 \times 3\) design krijg je \(df_A = 1\), \(df_B = 2\), \(df_{AB} = 1 \cdot 2 = 2\), \(df_E = N - 6\). Bij \(3 \times 4\): \(df_A = 2\), \(df_B = 3\), \(df_{AB} = 6\), \(df_E = N - 12\). De formule verandert niet — alleen de getallen.
Het kruispunt is waar A en B elkaar ontmoeten. Je ziet er vier paden (de cellen), twee marginal-richtingen (rij en kolom) en het centrale punt waar alles samenkomt (de grand mean). Hieronder twee zienswijzen op dezelfde data — de tabel en het profile plot. Kies wat je het meest helpt.
De 2D-tabel toont elk celgemiddelde plus de marginalen (rij-marg = effect van gang_diepte = factor A; kolom-marg = effect van dagdeel = factor B); de profile plot dezelfde getallen visueel. Probeer steeds heen-en-weer te zien: een rij-vergelijking in de tabel = een verschuiving van de hele lijn-set omhoog/omlaag in de plot; een kolom-vergelijking = een verschuiving links↔︎rechts; interactie = niet-parallelle lijnen.
apa_marginals_table(
m_factorial,
factor_A = "gang_diepte", # A = rij-factor
factor_B = "dagdeel", # B = kolom-factor
data = mol_werkdag,
dependent_label = "graaftempo (cm/uur)",
factor_A_label = "Gangdiepte",
factor_B_label = "Dagdeel",
caption = "Mol — graaftempo per cel met EMM-marginalen"
)| Gangdiepte |
Dagdeel
|
Marginaal B | |
|---|---|---|---|
| Ochtend | Avond | ||
| Ondiep | 25.88 (9.15, 52) | 33.98 (9.42, 47) | 29.93 |
| Δ Ondiep − Diep | 6.46 | 12.68 | — |
| Diep | 19.43 (9.43, 49) | 21.30 (10.12, 50) | 20.36 |
| Marginaal A | 22.66 | 27.64 | 25.15 |
| Note. Mol — graaftempo per cel met EMM-marginalen. Cellen tonen M (SD, n) op basis van graaftempo (cm/uur); rij- en kolom-marginalen zijn ongewogen estimated marginal means (EMM); rechtsonder de grand mean. Hoofdeffect Gangdiepte: F(1, 194) = 49.77, p < .001 — H0: alle rij-marginalen gelijk (μA⋅’s gelijk). Hoofdeffect Dagdeel: F(1, 194) = 13.50, p < .001 — H0: alle kolom-marginalen gelijk (μ⋅B’s gelijk). Interactie Gangdiepte × Dagdeel: F(1, 194) = 5.26, p .023 — H0: cel-patronen additief (geen interactie). | |||
profile_plot_with_marginals(
m_factorial,
factor_A = "gang_diepte",
factor_B = "dagdeel",
data = mol_werkdag,
dependent_symbol = "graaftempo",
dependent_label = "Graaftempo (cm/uur)",
factor_A_label = "Gangdiepte",
factor_B_label = "Dagdeel"
)
De interactie-toets vraagt: “is het verschil tussen \(A_1\) en \(A_2\) ongelijk over de niveaus van \(B\)?” In de kruispunt-tabel zie je dat aan de bold verschilrij: als die getallen onderling verschillen, is er interactie. Als ze ongeveer gelijk zijn, is er geen interactie.
Concreet bij de mol: \(\Delta_{\text{ochtend}} = 6.46\), \(\Delta_{\text{avond}} = 12.68\). Het verschil-tussen-verschillen is \(12.68 - 6.46 = 6.22\). Dát is wat de interactie-toets onderzoekt — al wordt het in \(F\)-vorm gerapporteerd. De interactie-coëfficiënt (\(b_{\text{int}}\)) in een regressie-parameterisering is precies dat verschil-tussen-verschillen, eventueel geschaald met de cel-codering.
Bij een \(3 \times 2\) design (zoals de dassen) zit dezelfde logica in de verschilkolom: de drie cel-verschillen \(\Delta_{\text{zand}}\), \(\Delta_{\text{klei}}\), \(\Delta_{\text{leem}}\) zijn ongelijk \(\Rightarrow\) interactie.
Met het kruispunt van hierboven (tabel en plot). Houd de A=rij, B=kolom-conventie aan: \(A_1 = \text{ondiep}\), \(A_2 = \text{diep}\) (rijen); \(B_1 = \text{ochtend}\), \(B_2 = \text{avond}\) (kolommen).
a) Welke twee paden vergelijk je voor het hoofdeffect van \(A\) (gang_diepte)? Schrijf ze op met hun symbool (\(\bar{Y}_{A_1 \cdot}\) etc.).
b) Welke twee paden voor het hoofdeffect van \(B\) (dagdeel)?
c) Op welk punt landt een mol waarvan we alleen weten dat hij in conditie \(B_1\) (ochtend) is geobserveerd? Wat is je beste voorspelling voor het graaftempo?
d) Diezelfde mol, maar nu weten we ook dat hij in een ondiepe gang zat (\(A_1\)). Op welke kruispunt-cel land je nu?
e) Wat zou je in de plot zien als de interactie precies nul was? Schets in woorden.
f) Welke kruispunt-cel valt op in de plot? Klopt dat met de \(F\)-waarde voor de interactie in de tabel?
a) De twee rij-marginalen (laatste kolom van de tabel, Marginaal B) vergelijk je: \(\bar{Y}_{A_1 \cdot}\) (ondiep, marginaal over dagdeel) versus \(\bar{Y}_{A_2 \cdot}\) (diep, marginaal over dagdeel). Concreet: \(29.93\) (ondiep) tegen \(20.36\) (diep) cm/uur.
b) De twee kolom-marginalen (laatste rij, Marginaal A): \(\bar{Y}_{\cdot B_1}\) (ochtend, marginaal over gang_diepte) versus \(\bar{Y}_{\cdot B_2}\) (avond). Concreet: \(22.66\) (ochtend) tegen \(27.64\) (avond).
c) Met alleen \(B_1\) bekend pak je de kolom-marginaal van \(B_1\): \(\hat{Y} = \bar{Y}_{\cdot B_1} = 22.66\) cm/uur.
d) Met zowel \(A_1\) als \(B_1\) bekend pak je de cel zelf: \(\hat{Y} = \bar{Y}_{A_1 B_1} = 25.88\) cm/uur. Voor ondiep-ochtend valt het iets boven de kolom-marginaal van \(B_1\) — dat past bij de extra informatie dat het ondiep was.
e) Bij interactie = nul lopen de lijnen van ondiep en diep parallel — het verschil tussen ochtend en avond zou voor ondiep en diep even groot zijn (bv. beide \(+5\), of beide \(0\)).
f) Vooral het kruispunt ondiep × avond trekt de aandacht: in de avond loopt de ondiep-lijn flink omhoog (van \(25.88\) naar \(33.98\)), terwijl de diep-lijn nauwelijks beweegt (van \(19.43\) naar \(21.30\)). De lijnen lopen niet parallel; in de avond is het diepte-effect dubbel zo groot als in de ochtend. Dat niet-parallel zijn correspondeert met \(F(1, 194) = 5.26\), \(p = .023\) voor de interactie — significant, dus het patroon is niet aan toeval toe te schrijven.
Vuistregel. Hoofdeffect = lijnen samen omhoog of omlaag. Interactie = lijnen niet-parallel. Tabel en plot zijn twee talen voor hetzelfde feit.
c) Is de aanname van gelijke varianties geschonden voor de vier cellen? Rapporteer Levene.
d) Welke effecten zijn significant? Rapporteer \(F\), \(df\) en \(p\) voor gang_diepte, dagdeel én interactie.
e) Hoeveel variantie verklaart elk effect (\(\eta^2\) en \(\eta^2_p\))?
f) Wat zie je in de interactieplot? Lopen de lijnen parallel, of niet?
In gewone woorden.
Levene niet significant (\(F(3, 194) = 0.70\), \(p = .56\)): aanname houdt stand.
In de Type III ANOVA-tabel: het hoofdeffect van gang_diepte is sterk significant (\(F(1, 194) = 49.77\), \(p < .001\)). Het hoofdeffect van dagdeel is ook significant, en duidelijk kleiner van omvang (\(F(1, 194) = 13.50\), \(p < .001\)). De interactie \(\text{gang\_diepte} \times \text{dagdeel}\) is klein-significant (\(F(1, 194) = 5.26\), \(p = .023\)).
Effectgroottes (we rapporteren primair kale \(\eta^2\), in lijn met de antwoord-PDF; \(\eta^2_p\) ter info erbij): \(\eta^2\) gang_diepte \(.189\) (groot), dagdeel \(.052\) (middelgroot), interactie \(.020\) (klein). Ter vergelijking: \(\eta^2_p\) ligt iets hoger (gang_diepte \(.202\), dagdeel \(.065\), interactie \(.026\)). De volgorde is duidelijk: gang_diepte verklaart veruit het meest, dan dagdeel, dan de interactie als kleinste effect.
In de interactieplot zie je de lijnen voor ondiep en diep niet parallel lopen: de ondiep-lijn stijgt sterker van ochtend naar avond dan de diep-lijn. Dat is precies wat een interactie visueel betekent: de twee factoren werken niet onafhankelijk.
APA-stijl.
Een \(2 \times 2\) tweeweg-ANOVA met gang_diepte en dagdeel als factoren toonde een significant en groot hoofdeffect van gang_diepte, \(F(1, 194) = 49.77\), \(p < .001\), \(\eta^2 = .189\). Het hoofdeffect van dagdeel was eveneens significant, met een middelgrote omvang, \(F(1, 194) = 13.50\), \(p < .001\), \(\eta^2 = .052\). De interactie tussen gang_diepte en dagdeel was klein maar significant, \(F(1, 194) = 5.26\), \(p = .023\), \(\eta^2 = .020\). De aanname van gelijke varianties was houdbaar, \(F(3, 194) = 0.70\), \(p = .56\).
“Als de interactie significant is, mag ik de hoofdeffecten niet meer in hun eentje interpreteren. Het diepte-effect is in de avond bijna twee keer zo groot als in de ochtend. Eén ongedifferentieerd hoofdeffect rapporteren zou hier liegen — alsof gang_diepte op elk moment van de dag hetzelfde doet.”
In een gebalanceerd balanced factorial design (gelijke \(n\) per cel) geven Type I, II en III sums of squares Type III sum of squares identieke \(F\)-waarden. In een ongebalanceerd unbalanced design (zoals de mol-data — lichtelijk ongebalanceerd, \(n\) varieert van 47 tot 52 per cel, \(n_{\max}/n_{\min} \approx 1.11\)) verschillen ze.
a) Waarom past Type III bij APA-rapportage van de hoofdeffecten in een ongebalanceerd factorial design? (Hint: denk aan “uniek na correctie voor de andere effecten”.)
b) Wat zou er fout kunnen gaan als je per ongeluk Type I gebruikt op ongebalanceerde data en de volgorde van factoren in je lm()-formule wijzigt?
a) Type III rekent voor elke term de unieke bijdrage na correctie voor alle andere termen in het model. Dat correspondeert met de gangbare APA-vraag: “is dit hoofdeffect / deze interactie nog significant als de andere effecten al meegenomen zijn?” Type III geeft de juiste \(F\)-toets bij ongebalanceerde data, mits de contrasten op contr.sum staan.
b) Type I sums of squares zijn sequentieel: de eerste factor krijgt al zijn variantie, de tweede alleen wat overblijft, en zo verder. In een ongebalanceerd design verandert daardoor het resultaat als je Y ~ A + B schrijft of Y ~ B + A. Twee onderzoekers zouden zonder dat ze het doorhebben verschillende \(F\)- en \(p\)-waarden voor dezelfde data rapporteren. Type II en Type III zijn invariant onder volgorde-veranderingen; Type III bij voorkeur omdat hij ook de interactie correct corrigeert.
Bij twee factoren met \(I\) niveaus (factor \(A\)) en \(J\) niveaus (factor \(B\)), en \(N\) totaal observaties:
| Bron | df |
|---|---|
| Factor \(A\) | \(I - 1\) |
| Factor \(B\) | \(J - 1\) |
| \(A \times B\) | \((I - 1)(J - 1)\) |
| Residu | \(N - IJ\) |
| Totaal | \(N - 1\) |
Som-check: \((I-1) + (J-1) + (I-1)(J-1) + (N - IJ) = N - 1\).
Hoofdrekenbaar voorbeeld. Een \(2 \times 3\) design met \(n = 10\) per cel: \(I = 2\), \(J = 3\), \(N = 60\). Dan \(df_A = 1\), \(df_B = 2\), \(df_{AB} = 2\), \(df_w = 60 - 6 = 54\), \(df_{\text{total}} = 59\). Som: \(1 + 2 + 2 + 54 = 59\).
Verifieer bij de fluittoon-data straks. Olifant (\(I = 2\) geslacht, \(J = 3\) wind, \(N = 134\)): \(df_A = 1\), \(df_B = 2\), \(df_{AB} = 2\), \(df_w = 134 - 6 = 128\), \(df_{\text{total}} = 133\). Check zo dadelijk in Extra 2A of de tweede df van elke \(F\)-toets inderdaad \(128\) is.
Wanneer de interactie significant is, ga je terug naar de cell means en beschrijf je het patroon zo concreet mogelijk. De cursus-route is: estimated marginal means per niveau van de andere factor (emmeans(model, ~ A | B)), aangevuld met de cell means uit het kruispunt en de interactieplot. Geen aparte simple-effects-toetsen; die term gebruiken we hier niet.
Een hoofdeffect (bv. “ondiepe gangen leveren gemiddeld meer cm/uur dan diepe”) is een uitspraak over de rij-marginalen marginal mean — het gemiddelde over alle \(B\)-condities. Maar als de interactie significant is, bestaat dat marginal-gemiddelde uit ongelijke patronen: het verschil tussen ondiep en diep is niet hetzelfde in de ochtend als in de avond.
De cursus-conventie is helder:
emmeans(model, ~ A | B)) — niet via aparte simple-effects-toetsen, die term gebruiken we hier niet. Plus de cell means uit het kruispunt en de interactieplot.Concreet voorbeeld: stel het verschil ondiep − diep is \(+6\) in de ochtend en \(-2\) in de avond. Het marginal-verschil is \(+2\) (gemiddeld over dagdelen). Een hoofdeffect-toets met \(p < .05\) rapporteer je wel, maar de inhoudelijke conclusie komt uit de cell-patronen — daar zien we dat ondiep alleen ’s ochtends echt voor ligt; in de avond zijn ze gelijk of zelfs omgekeerd. Dat is wat de interactie zegt.
Bron-formulering: het Exercise Book (p. 51, p. 54) gebruikt consequent “qualified by” voor deze relatie. Onthoud die uitdrukking voor je APA-zin.
Bij een ongebalanceerd design lopen twee soorten gemiddelden uit elkaar:
aggregate() of mean() op de raw data toepast.emmeans() rapporteert.Bij balanced (alle \(n_{ij}\) gelijk): de twee zijn identiek. Bij unbalanced: ze divergeren — soms minimaal, soms substantieel. ANOVA met Type III SS toetst altijd de EMM-verschillen — dat zijn de “schone” effecten waarbij elke cel-conditie gelijke invloed heeft op het marginal-gemiddelde, ongeacht hoeveel observaties er toevallig in zaten.
Concreet: stel je hebt \(n_{\text{ondiep}} = 100\) mollen in een ondiepe gang (waarvan 80 in conditie A, 20 in B) en \(n_{\text{diep}} = 20\) mollen in een diepe gang (waarvan 5 in A, 15 in B). De descriptives-mean van A trekt zwaar naar de ondiep-conditie-A (80 van de 85 observaties komen daar vandaan). De EMM van A is gewoon \((\bar{Y}_{\text{ondiep,A}} + \bar{Y}_{\text{diep,A}})/2\) — een eerlijke vergelijking ongeacht de mix.
Daarom: rapporteer in een ANOVA-context EMM voor marginal-effecten, niet de descriptives-tabel. emmeans(model, ~ factor_A) geeft je die direct.
# Descriptives (gewogen, weegt cel-grootte mee).
aggregate(graaftempo ~ gang_diepte, mol_werkdag, FUN = mean) gang_diepte graaftempo
1 ondiep 29.72727
2 diep 20.37374# EMM (ongewogen, gelijke weeg per cel).
emmeans(m_factorial, ~ gang_diepte)NOTE: Results may be misleading due to involvement in interactions gang_diepte emmean SE df lower.CL upper.CL
ondiep 29.9 0.960 194 28.0 31.8
diep 20.4 0.958 194 18.5 22.3
Results are averaged over the levels of: dagdeel
Confidence level used: 0.95 In de mol-data is het verschil bescheiden — het design is lichtelijk ongebalanceerd (\(n_{\max}/n_{\min} \approx 1.11\)), dus de twee soorten gemiddelden divergeren maar een fractie (\(\bar{Y}^{\text{gewogen}}_{\text{ondiep}} = 29.73\) tegenover \(\text{EMM}_{\text{ondiep}} = 29.93\)). Bij sterk ongebalanceerde designs kan dat verschil groot worden — en dan rapporteer je de EMM, niet de descriptives.
# Estimated marginal means per cel — EMM-uitwerking per niveau.
# Effect van dagdeel binnen elk niveau van gang_diepte:
emm <- emmeans(m_factorial, ~ dagdeel | gang_diepte)
emmgang_diepte = ondiep:
dagdeel emmean SE df lower.CL upper.CL
ochtend 25.9 1.32 194 23.3 28.5
avond 34.0 1.39 194 31.2 36.7
gang_diepte = diep:
dagdeel emmean SE df lower.CL upper.CL
ochtend 19.4 1.36 194 16.7 22.1
avond 21.3 1.35 194 18.6 24.0
Confidence level used: 0.95 # Paarsgewijze vergelijking van dagdeel binnen elke gang_diepte.
pairs(emm)gang_diepte = ondiep:
contrast estimate SE df t.ratio p.value
ochtend - avond -8.09 1.92 194 -4.218 <0.0001
gang_diepte = diep:
contrast estimate SE df t.ratio p.value
ochtend - avond -1.87 1.92 194 -0.976 0.3301g) Beschrijf in eigen woorden — niet via de cijfers, maar via wat de mol doet — wat de interactie betekent.
h) Bij welke gang_diepte is het effect van dagdeel significant (EMM-vergelijking per niveau van gang_diepte)?
In gewone woorden (g). In een ondiepe gang verschilt het tempo flink tussen ochtend en avond — ’s avonds graaft de mol daar duidelijk sneller. In een diepe gang maakt het dagdeel nauwelijks uit; de compactere bodem laat zich niet door tijd-van-de-dag overhalen. Het effect van dagdeel hangt dus af van waar hij zit.
Een mogelijke (speculatieve!) interpretatie: in losse, ondiepe grond profiteert hij van de avondkoelte of van het nachtritme; in compacte, diepe bodem wordt het tempo gelimiteerd door de bodem zelf, niet door zijn dagschema.
EMM-vergelijking per niveau (h). Bij ondiep is het verschil tussen ochtend en avond significant (\(t(194) = -4.22\), \(p < .001\), een verschil van \(-8.09\) cm/uur — ochtend trager dan avond). Bij diep is het verschil niet significant (\(t(194) = -0.98\), \(p = .33\), verschil \(\approx -1.87\) cm/uur). Het significante interactie-effect zit dus vooral in de ondiep-cellen.
APA-stijl.
De significante interactie weerspiegelde dat het effect van dagdeel op het graaftempo afhing van gang_diepte. In ondiepe gangen lag het tempo in de avond (\(M = 33.98\)) circa \(8\) cm/uur boven dat in de ochtend (\(M = 25.88\)), en dat verschil was significant (\(t(194) = -4.22\), \(p < .001\)); in diepe gangen werd geen significant dagdeel-verschil gevonden (\(M_{\text{avond}} = 21.30\); \(M_{\text{ochtend}} = 19.43\); \(t(194) = -0.98\), \(p = .33\)).
Wanneer je \(k\) groepen vergelijkt en de overall \(F\) niet significant is, moet je niet alsnog Tukey of een ander post-hoc draaien om “te zien of er ergens toch een verschil zit”. Dat verhoogt je Type-I-foutrisico enorm. De \(F\)-toets is je gatekeeper: opent hij de deur, dan mag je naar binnen om te kijken welke paren verschillen. Sluit hij hem, dan is het stop.
Uitzondering: in factorial designs kun je per significant hoofdeffect of significante interactie post-hoc gaan kijken — niet over de hele tabel, alleen waar het ANOVA-bewijs het rechtvaardigt.
Als je \(k\) paren tegelijk vergelijkt op \(\alpha = .05\), dan is de kans dat minimaal één vals-positief is, gelijk aan \(1 - (1 - .05)^k\). Voor 3 paren is dat \(1 - .95^3 \approx .14\). Voor 6 paren al \(.26\). Daarom: corrigeer.
Methoden:
Vuistregel: Tukey voor “alle paren”. Bonferroni als veilige default voor andere multiple-testing-situaties. Holm als je Bonferroni-power te krap vindt.
In een tweeweg-ANOVA print etaSquared() twee kolommen: eta.sq (\(\eta^2\)) en eta.sq.part (\(\eta^2_p\)).
\[ \eta^2 = \frac{\text{SS}_{\text{effect}}}{\text{SS}_{\text{totaal}}}, \qquad \eta^2_p = \frac{\text{SS}_{\text{effect}}}{\text{SS}_{\text{effect}} + \text{SS}_{\text{residuen}}} \]
a) Waarom is \(\eta^2_p \geq \eta^2\) altijd?
b) Bij de mol-tweeweg-ANOVA: \(\text{SS}_{\text{gang\_diepte}} = 4525\), \(\text{SS}_{\text{dagdeel}} = 1227\), \(\text{SS}_{\text{interactie}} = 479\), \(\text{SS}_{\text{residu}} = 17637\). Bereken \(\eta^2\) en \(\eta^2_p\) voor de interactie. (Tip: bij \(\eta^2\) is de noemer de som van alle effect-SS plus residu-SS; in een gebalanceerd design \(\approx \text{SS}_{\text{totaal}}\).)
a) \(\eta^2_p\) heeft alleen “effect + residu” in de noemer, terwijl \(\eta^2\) ook andere effecten in de noemer meeneemt (groter). Voor dezelfde teller geeft dat \(\eta^2_p \geq \eta^2\). Bij eenweg-ANOVA geldt \(\eta^2 = \eta^2_p\) (geen andere effecten); bij meerweg-ANOVA loopt het uiteen.
b) Totale SS = \(4525 + 1227 + 479 + 17637 = 23868\).
Kale \(\eta^2\) (primair — wat de antwoord-PDF rapporteert): \(\eta^2_{\text{interactie}} = 479 / 23868 = 0.0201 = .020\).
Partial \(\eta^2_p\) (ter vergelijking): \(\eta^2_{p,\text{interactie}} = 479 / (479 + 17637) = 479 / 18116 = 0.0264 = .026\).
Past bij wat etaSquared(m_factorial) print.
Inzicht. Waarom rapporteer je dan welke? In de antwoord-PDF is kale \(\eta^2\) consequent het primaire getal — die matcht ook direct Cohens vuistwaarden (\(.01/.06/.14\)). \(\eta^2_p\) is gangbaar in SPSS-output (SPSS toont alleen \(\eta^2_p\)) en bij meta-analyses. Wij geven hier beide; in modelantwoorden zetten we \(\eta^2\) voorop en noemen \(\eta^2_p\) secundair.
Waarom partial eta-squared kan optellen tot meer dan een
Hoofdrekenbaar voorbeeld om het verschil scherp te krijgen. Stel je hebt een tweeweg-ANOVA met \(\text{SS}_A = 30\), \(\text{SS}_B = 20\), \(\text{SS}_{AB} = 10\), \(\text{SS}_w = 40\), en dus \(\text{SS}_{\text{totaal}} = 100\).
\(\eta^2\): één totaal-cirkel van \(100\%\) met vier disjuncte segmenten. Elk effect krijgt zijn eigen plak van de hele taart.
\[ \eta^2_A = \tfrac{30}{100} = .30, \quad \eta^2_B = \tfrac{20}{100} = .20, \quad \eta^2_{AB} = \tfrac{10}{100} = .10 \]
Som: \(.60\) — past binnen \(1\), zoals het hoort bij een netjes verdeelde taart.
library(eulerr)
fit_eta2 <- euler(c(
"A" = 30,
"B" = 20,
"AB" = 10,
"Residu" = 40
))
plot(fit_eta2,
fills = list(fill = c("#0077BB", "#CCBB44", "#EE7733", "#BBBBBB"),
alpha = 0.55),
labels = list(col = "black", font = 2),
quantities = list(type = "counts", cex = 0.9),
edges = list(col = "white", lwd = 2))
\(\eta^2_p\): drie aparte mini-Venns — per effect tel je alleen dát effect plus residu mee. De andere effecten worden weggepoetst voor die berekening, alsof ze er voor deze effectgrootte even niet waren.
\[ \eta^2_{p,A} = \tfrac{30}{30 + 40} = .43, \quad \eta^2_{p,B} = \tfrac{20}{20 + 40} = .33, \quad \eta^2_{p,AB} = \tfrac{10}{10 + 40} = .20 \]
Som: \(.96\) — bijna over de \(100\%\) heen. In andere voorbeelden kan \(\sum \eta^2_p\) wél boven \(1\) uitkomen; dat is geen fout, het laat zien dat elk effect zijn eigen “submap” krijgt waarin de andere effecten zijn weggehaald.
library(eulerr)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
fit_pA <- euler(c("A" = 30, "Residu" = 40))
plot(fit_pA, main = expression(eta[p]^2 * "(A) = 30/70 = .43"),
fills = list(fill = c("#0077BB", "#BBBBBB"), alpha = 0.55),
labels = list(col = "black", font = 2),
quantities = list(type = "counts", cex = 0.9),
edges = list(col = "white", lwd = 2))
fit_pB <- euler(c("B" = 20, "Residu" = 40))
plot(fit_pB, main = expression(eta[p]^2 * "(B) = 20/60 = .33"),
fills = list(fill = c("#CCBB44", "#BBBBBB"), alpha = 0.55),
labels = list(col = "black", font = 2),
quantities = list(type = "counts", cex = 0.9),
edges = list(col = "white", lwd = 2))
fit_pAB <- euler(c("AB" = 10, "Residu" = 40))
plot(fit_pAB, main = expression(eta[p]^2 * "(A" %*% "B) = 10/50 = .20"),
fills = list(fill = c("#EE7733", "#BBBBBB"), alpha = 0.55),
labels = list(col = "black", font = 2),
quantities = list(type = "counts", cex = 0.9),
edges = list(col = "white", lwd = 2))
par(mfrow = c(1, 1))\(\eta^2\) deelt één taart van \(100\%\) over alle effecten plus residu — de proporties tellen netjes op tot \(\leq 1\). Dit is wat de antwoord-PDF rapporteert en wat hier je primaire effectgrootte is. \(\eta^2_p\) vraagt per effect: “hoeveel verklaar ik als de andere effecten er even niet waren?” Daardoor mag de som boven \(1\) uitkomen — handig voor vergelijking met SPSS-output of meta-analyses, maar in modelantwoorden secundair.
Het balanced/unbalanced-onderscheid is geen hoofddoel van dit thema, maar het is wel iets om te kunnen herkennen en kort uitleggen. Sommige opleidingen leggen er nadruk op (vroeger was het zelfs een tentamenklassieker), andere stippen het alleen aan. Check je eigen college-sheets om te weten waar jouw docent op rekent.
Een factorial design is balanced als de aantallen per cel gelijk zijn (\(n_{ij}\) overal gelijk). Bij balanced designs:
Bij unbalanced designs (verschillende cel-grootten) overlappen de effecten in het verklaarde-variantie-deel. Daarom maakt het uit hoe je de SS toekent (Type I/II/III) — en welke \(F\) je rapporteert. In de praktijk is bijna elk psychologisch onderzoek lichtelijk ongebalanceerd (drop-outs, niet-respons), dus Type III is meestal de juiste keuze.
Route 1 (telling). Tel de cellen met table(factor1, factor2). Alle cellen gelijk in \(n\)? Balanced.
Route 2 (truc). Bij balanced geldt \(n_{ij} = n_{i \cdot} \cdot n_{\cdot j} / N\) voor élke cel — de cel-aantallen zijn precies het product van de marginal counts gedeeld door totaal. Klopt dat niet voor één cel? Dan is het design unbalanced. Handig als je alleen marginal counts hebt.
# Route 1: cel-tellingen direct.
table(mol_werkdag$gang_diepte, mol_werkdag$dagdeel)
ochtend avond
ondiep 52 47
diep 49 50# Route 2: voorspel cel-aantallen uit marginal counts.
n_total <- nrow(mol_werkdag)
n_diep <- table(mol_werkdag$gang_diepte)
n_dag <- table(mol_werkdag$dagdeel)
verwacht_balanced <- outer(n_diep, n_dag) / n_total
round(verwacht_balanced, 2)
ochtend avond
ondiep 50.5 48.5
diep 50.5 48.5De werkelijke table() geeft \(52 / 47 / 49 / 50\) — de “verwacht-balanced”-tabel zou \(50.5 / 48.5 / 50.5 / 48.5\) geven. De cijfers liggen dichtbij elkaar maar niet exact gelijk: het mol-design is lichtelijk ongebalanceerd (\(n_{\max}/n_{\min} = 52/47 \approx 1.11\), ruim binnen de \(1.5\)-robuustheidsgrens). Voor de praktijk: Type III aanhouden, EMM rapporteren, en je bent klaar.
Waar Type I, II, III hetzelfde of verschillend antwoorden
Een hoofdrekenbaar voorbeeld om de twee situaties naast elkaar te zien.
Balanced — twee disjuncte cirkels binnen de totaal-ruimte. \(\text{SS}_A = 30\), \(\text{SS}_B = 20\), geen overlap. Type I, Type II en Type III geven identiek antwoord, want er is geen gedeeld stuk om over te onderhandelen.
library(eulerr)
fit_bal <- euler(c(
"A" = 30,
"B" = 20,
"A&B" = 0
))
plot(fit_bal,
fills = list(fill = c("#0077BB", "#CCBB44"), alpha = 0.55),
labels = list(col = "black", font = 2),
quantities = list(type = "counts", cex = 0.9),
edges = list(col = "white", lwd = 2))
Unbalanced — cirkels overlappen: \(\text{SS}_A^{\text{uniek}} = 24\), \(\text{SS}_B^{\text{uniek}} = 16\), en een gedeeld stuk \(\text{SS}_{A \cap B} = 6\).
Hand-rek-regel: \(30 = 24 + 6\), \(22 = 16 + 6\). De overlap reist mee met wie als eerste binnenstapt in Type I.
library(eulerr)
fit_unb <- euler(c(
"A" = 24,
"B" = 16,
"A&B" = 6
))
plot(fit_unb,
fills = list(fill = c("#0077BB", "#CCBB44"), alpha = 0.55),
labels = list(col = "black", font = 2),
quantities = list(type = "counts", cex = 0.9),
edges = list(col = "white", lwd = 2))
Bij balanced raken de cirkels elkaar niet — er valt niets te onderhandelen, dus alle SS-types vallen samen. Bij unbalanced is er een gedeeld stuk; wie dat krijgt hangt af van de SS-keuze. Type III zegt: “geen van beide” — alleen de unieke bijdragen. Vandaar dat Type III de standaard is voor APA-rapportage in ongebalanceerde designs.
Een tweede onderzoek, een nieuwe gast
Tegen het einde van de middag kwam de olifant langs. “Mol,” zei hij, “ik heb ook een onderzoekje gedaan. Ik wilde weten hoe snel dieren reageren als ik een fluittoon laat horen — en of het uitmaakt of er wind staat.”
De olifant heeft drie wind-condities (geen wind, licht briesje, harde wind), honderdvierendertig dieren, en hij hield bij of ze mannetjes of vrouwtjes waren. Per dier mat hij hoe lang het duurde voor het reageerde op de fluittoon (in milliseconden).
Dit is een \(2 \times 3\) factorial design — meer levels op één van de factoren. De analyse is structureel hetzelfde als 2.2, maar met één factor die drie niveaus heeft: dat geeft de mogelijkheid (en bij significantie: noodzaak) van post-hoc tests om uit te zoeken welke wind-niveaus van elkaar verschillen.
load("data/fluittoon_reactie.RData")
str(fluittoon_reactie)'data.frame': 134 obs. of 3 variables:
$ reactietijd: num 4143 5038 7894 1919 5943 ...
$ geslacht : Factor w/ 2 levels "mannetje","vrouwtje": 2 2 1 1 2 1 1 1 2 1 ...
$ wind : Factor w/ 3 levels "geen wind","licht briesje",..: 1 3 1 2 3 2 3 3 3 2 ...# Cell means in milliseconden
aggregate(reactietijd ~ geslacht * wind, data = fluittoon_reactie, FUN = mean) geslacht wind reactietijd
1 mannetje geen wind 3097.841
2 vrouwtje geen wind 2656.381
3 mannetje licht briesje 4056.794
4 vrouwtje licht briesje 5547.829
5 mannetje harde wind 6843.620
6 vrouwtje harde wind 4838.462# Aantallen per cel
table(fluittoon_reactie$geslacht, fluittoon_reactie$wind)
geen wind licht briesje harde wind
mannetje 25 21 24
vrouwtje 21 25 18leveneTest(reactietijd ~ geslacht * wind, data = fluittoon_reactie)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 0.1941 0.9643
128 m_olifant <- lm(reactietijd ~ geslacht * wind, data = fluittoon_reactie)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
options(contrasts = c("contr.sum", "contr.poly"))
Anova(m_olifant, type = 3)Anova Table (Type III tests)
Response: reactietijd
Sum Sq Df F value Pr(>F)
(Intercept) 2683736778 1 438.8361 < 2.2e-16 ***
geslacht 3351448 1 0.5480 0.460484
wind 198203075 2 16.2048 0.000000533 ***
geslacht:wind 66695515 2 5.4529 0.005337 **
Residuals 782794141 128
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1etaSquared(m_olifant) eta.sq eta.sq.part
geslacht 0.002112573 0.002875488
wind 0.202728058 0.216752177
geslacht:wind 0.062416391 0.078512451# Tijdelijke factoren met hoofdletter-niveaus voor leesbare as- en
# legend-labels — onderliggende data blijft ongewijzigd.
.wind_lbl <- factor(tools::toTitleCase(as.character(fluittoon_reactie$wind)),
levels = tools::toTitleCase(levels(fluittoon_reactie$wind)))
.geslacht_lbl <- factor(tools::toTitleCase(as.character(fluittoon_reactie$geslacht)),
levels = tools::toTitleCase(levels(fluittoon_reactie$geslacht)))
interaction.plot(
x.factor = .wind_lbl,
trace.factor = .geslacht_lbl,
response = fluittoon_reactie$reactietijd,
ylab = "Gemiddelde reactietijd (ms)",
xlab = "Wind",
trace.label = "Geslacht"
)
Hetzelfde patroon als bij de mol, nu voor het 2 × 3-design met de olifant.
apa_marginals_table(
m_olifant,
factor_A = "geslacht",
factor_B = "wind",
data = fluittoon_reactie,
dependent_label = "reactietijd (ms)",
factor_A_label = "Geslacht",
factor_B_label = "Wind",
caption = "Olifant — reactietijd per cel met EMM-marginalen"
)| Geslacht |
Wind
|
Marginaal B | ||
|---|---|---|---|---|
| Geen wind | Licht briesje | Harde wind | ||
| Mannetje | 3097.84 (2644.86, 25) | 4056.79 (2350.14, 21) | 6843.62 (2150.03, 24) | 4666.08 |
| Δ Mannetje − Vrouwtje | 441.46 | -1491.04 | 2005.16 | — |
| Vrouwtje | 2656.38 (2501.64, 21) | 5547.83 (2466.97, 25) | 4838.46 (2732.14, 18) | 4347.56 |
| Marginaal A | 2877.11 | 4802.31 | 5841.04 | 4506.82 |
| Note. Olifant — reactietijd per cel met EMM-marginalen. Cellen tonen M (SD, n) op basis van reactietijd (ms); rij- en kolom-marginalen zijn ongewogen estimated marginal means (EMM); rechtsonder de grand mean. Hoofdeffect Geslacht: F(1, 128) = 0.55, p .46 — H0: alle rij-marginalen gelijk (μA⋅’s gelijk). Hoofdeffect Wind: F(2, 128) = 16.20, p < .001 — H0: alle kolom-marginalen gelijk (μ⋅B’s gelijk). Interactie Geslacht × Wind: F(2, 128) = 5.45, p .005 — H0: cel-patronen additief (geen interactie). | ||||
profile_plot_with_marginals(
m_olifant,
factor_A = "geslacht",
factor_B = "wind",
data = fluittoon_reactie,
dependent_symbol = "reactietijd",
dependent_label = "Reactietijd (ms)",
factor_A_label = "Geslacht",
factor_B_label = "Wind"
)
# Tukey HSD over de wind-factor — mag alleen omdat
# het hoofdeffect van wind significant is.
TukeyHSD(aov(reactietijd ~ wind, data = fluittoon_reactie)) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = reactietijd ~ wind, data = fluittoon_reactie)
$wind
diff lwr upr p adj
licht briesje-geen wind 1970.834 710.3904 3231.278 0.0008960
harde wind-geen wind 3087.962 1797.8561 4378.067 0.0000003
harde wind-licht briesje 1117.127 -172.9782 2407.233 0.1038851a) Welke effecten zijn significant? Rapporteer \(F\), \(df\), \(p\) en \(\eta^2\) per effect.
b) Welke wind-paren verschillen significant volgens Tukey HSD?
c) Beschrijf de interactie in woorden: hoe verandert het effect van wind tussen mannetjes en vrouwtjes?
In gewone woorden. Levene niet significant (\(F(5, 128) = 0.19\), \(p = .96\)): de cellen hebben gelijke varianties.
Het hoofdeffect van wind is zeer sterk: \(F(2, 128) = 16.20\), \(p < .001\), \(\eta^2 = .200\) — een groot effect. Het hoofdeffect van geslacht is niet significant: \(F(1, 128) = 0.55\), \(p = .460\), \(\eta^2 = .002\). De interactie is wel significant: \(F(2, 128) = 5.45\), \(p = .005\), \(\eta^2 = .062\) — middelgroot.
Tukey HSD over wind: licht briesje versus geen wind — significant (\(p < .001\), verschil \(\approx 1971\) ms). Harde wind versus geen wind — significant (\(p < .001\), verschil \(\approx 3088\) ms). Harde wind versus licht briesje — niet significant (\(p = .100\)). Het grote contrast zit dus tussen geen wind en de twee wind-condities; de twee wind-condities onderling verschillen minder duidelijk.
De interactie: bij geen wind zijn mannetjes en vrouwtjes vergelijkbaar (~3098 vs 2656 ms). Bij licht briesje reageren vrouwtjes juist trager dan mannetjes (5548 vs 4057 ms). Bij harde wind zijn de rollen weer omgedraaid: mannetjes reageren nu trager dan vrouwtjes (6844 vs 4838 ms). Het patroon is niet uniform — wind beïnvloedt mannetjes en vrouwtjes anders, en zelfs in tegengestelde richtingen voor verschillende wind-niveaus.
APA-stijl.
Een \(2 \times 3\) tweeweg-ANOVA met geslacht en wind als factoren op reactietijd liet een significant hoofdeffect van wind zien, \(F(2, 128) = 16.20\), \(p < .001\), \(\eta^2 = .200\), en een significante interactie tussen geslacht en wind, \(F(2, 128) = 5.45\), \(p = .005\), \(\eta^2 = .062\). Het hoofdeffect van geslacht was niet significant, \(F(1, 128) = 0.55\), \(p = .460\), \(\eta^2 = .002\). De aanname van gelijke varianties was houdbaar, \(F(5, 128) = 0.19\), \(p = .96\). Tukey-HSD-vergelijkingen over de wind-condities lieten zien dat dieren bij geen wind significant sneller reageerden (\(M = 2896\) ms) dan bij een licht briesje (\(M = 4867\) ms; \(p < .001\)) en bij harde wind (\(M = 5984\) ms; \(p < .001\)). De twee wind-condities verschilden onderling niet significant van elkaar (\(p = .100\)).
In de fluittoon-data is het hoofdeffect van wind (3 niveaus) significant; we kiezen Tukey HSD voor paarsgewijze vergelijkingen.
a) Hoeveel paren wind-condities zijn er om te vergelijken? Hoeveel paren in totaal voor een 4-niveau-factor? (Algemeen: \(\binom{k}{2}\).)
b) Tukey HSD past de \(p\)-waarde aan voor het aantal vergelijkingen — familywise error rate FWER. Wat zou er gebeuren als je in plaats daarvan drie aparte \(t\)-toetsen zou doen, allemaal op \(\alpha = .05\), zonder correctie?
c) Welke andere post-hoc-procedures ken je naast Tukey, en wanneer kies je iets anders?
a) Bij 3 niveaus: \(\binom{3}{2} = 3\) paren (geen-licht, geen-hard, licht-hard). Bij 4 niveaus: \(\binom{4}{2} = 6\) paren. Het aantal explodeert met \(k\).
b) Met drie ongecorrigeerde toetsen op \(\alpha = .05\) stijgt de kans op minstens één Type-I-fout naar \(1 - (1 - .05)^3 \approx .14\). Bij zes toetsen \(\approx .26\). Je vindt veel “significant” door toevallige variatie.
c) De vier kern-procedures, plus Holm als sequentiële variant:
aov()-objecten.Daarnaast: Šidák (iets minder streng dan Bonferroni bij onafhankelijke toetsen), Games-Howell (Tukey-variant bij heteroscedasticiteit) en Benjamini-Hochberg / FDR (controleert false discovery rate i.p.v. FWER, gangbaar bij grote testfamilies). Check je eigen vakgroep welke procedure de voorkeur heeft.
De olifant kwam terug, een paar dagen later. “Mol, ik wil de eenweg-versie zien — alleen wind, drie groepen, geen geslacht. Dat hoort bij de canon.”
Soms is de tweeweg-vraag een omweg en wil je gewoon: één factor met drie of meer niveaus, en dan kijken of er groepsverschillen zijn. We gebruiken dezelfde fluittoon-data, maar nu alleen met wind als factor.
m_wind <- lm(reactietijd ~ wind, data = fluittoon_reactie)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
options(contrasts = c("contr.sum", "contr.poly"))
Anova(m_wind, type = 3)Anova Table (Type III tests)
Response: reactietijd
Sum Sq Df F value Pr(>F)
(Intercept) 2808824185 1 432.002 < 2.2e-16 ***
wind 216810636 2 16.673 0.0000003541 ***
Residuals 851747063 131
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1etaSquared(m_wind) eta.sq eta.sq.part
wind 0.2029003 0.2029003# Cellen
aggregate(reactietijd ~ wind, data = fluittoon_reactie, FUN = mean) wind reactietijd
1 geen wind 2896.305
2 licht briesje 4867.139
3 harde wind 5984.267aov() en waarom verschilt dat van lm() plus Anova()?
aov() fit hetzelfde model als lm(), maar print zijn summary() als ANOVA-tabel met Type I SS (sequentieel). Voor eenweg-designs of gebalanceerde tweeweg-designs maakt het niet uit. Voor ongebalanceerde tweeweg-designs gebruik je lm() + Anova(., type = 3). TukeyHSD() werkt alleen op een aov()-object — vandaar dat we voor post-hoc óf aov() óf emmeans() aanroepen.
TukeyHSD(aov(reactietijd ~ wind, data = fluittoon_reactie)) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = reactietijd ~ wind, data = fluittoon_reactie)
$wind
diff lwr upr p adj
licht briesje-geen wind 1970.834 710.3904 3231.278 0.0008960
harde wind-geen wind 3087.962 1797.8561 4378.067 0.0000003
harde wind-licht briesje 1117.127 -172.9782 2407.233 0.1038851Behalve “alle paren vergelijken” kun je ook gerichte contrasten kiezen. Stel je hebt drie wind-niveaus en je wilt twee specifieke vragen beantwoorden:
Twee contrasten zijn orthogonaal orthogonal als \(\sum c_{1k} \cdot c_{2k} = 0\). Bij drie niveaus passen er maximaal \(k - 1 = 2\) orthogonale contrasten in; samen gebruiken ze alle vrijheidsgraden voor de tussen-groepen-variatie. Orthogonale contrasten zijn elkaars onafhankelijke vragen — ze dubbelen niet. Wanneer je deze contrasten van tevoren vastlegt op grond van theorie heten ze geplande contrasten planned contrasts.
emmeans() plus contrast() laat je dit uitvoeren:
emm <- emmeans(m_wind, ~ wind)
contrast(emm, list(
wind_vs_geen = c(-2, 1, 1) / 2, # gemiddelde wind - geen wind
licht_vs_hard = c( 0,-1, 1) # licht briesje - harde wind
))emm <- emmeans(m_wind, ~ wind)
contrast(emm, list(
wind_vs_geen = c(-2, 1, 1) / 2,
licht_vs_hard = c( 0, -1, 1)
)) contrast estimate SE df t.ratio p.value
wind_vs_geen 2529 464 131 5.450 <0.0001
licht_vs_hard 1117 544 131 2.053 0.0421a) Welk contrast is significant, welk niet? Wat zegt dat over de structuur van het wind-effect?
In gewone woorden. Het contrast wind versus geen wind is sterk significant: dieren reageren bij wind veel trager dan zonder wind. Het contrast licht briesje versus harde wind is óók significant (\(p = .040\)): bij licht briesje reageert men sneller dan bij harde wind. De grote sprong zit tussen geen wind en welke wind dan ook, en daarbinnen is er nóg een kleinere — maar wel detecteerbare — sprong tussen lichte en harde wind.
APA-zin. Twee geplande contrasten lieten zien dat reactietijden in de wind-conditie (licht briesje en harde wind samen) significant verlengd waren ten opzichte van de geen-wind-conditie (\(t(131) = -5.66\), \(p < .001\)). Het tweede contrast — licht briesje versus harde wind — was eveneens significant (\(t(131) = 2.05\), \(p = .042\)), met snellere reacties bij licht briesje dan bij harde wind.
Let op — Tukey vs. geplande contrasten. De Tukey-HSD-vergelijking van licht briesje versus harde wind gaf eerder \(p = .100\) (niet significant). Het geplande contrast op dezelfde paren geeft \(p = .040\) (wel significant). Beide kloppen — ze beantwoorden net iets andere vragen.
Tukey corrigeert voor het feit dat je drie paren tegelijk vergelijkt en zet de drempel hoger om de Type-I-fout over de hele familie binnen \(\alpha = .05\) te houden (familywise error rate, FWER). Geplande contrasten doen die correctie niet, omdat je vooraf — voordat je de data zag — specifiek déze twee onderzoeksvragen formuleerde. Dat is precies hun didactisch nut: gerichte hypotheses kosten minder power dan een blinde paarsgewijze zoektocht.
Welke is “juist”? Beide. Tukey is exploratief-paarsgewijs en streng; geplande contrasten zijn specifiek en a priori. Dit is een leermoment: hetzelfde paar getallen, twee verschillende correctielogica’s, twee verschillende conclusies — en allebei goed te verdedigen, mits je vooraf je strategie hebt vastgelegd.
Een onderzoeker ontwerpt een \(2 \times 3\) factorial design. Cohen (1988) suggereert voor \(\eta^2_p \approx .06\) (middelgroot) en \(\alpha = .05\) ongeveer \(25\)–\(30\) observaties per cel om \(.80\) power te halen.
a) Hoeveel proefpersonen totaal heeft hij minimaal nodig?
b) Stel hij heeft \(40\) proefpersonen totaal. Wat zijn zijn opties?
c) Wat is volgens Tabachnick & Fidell (2007) een te ongebalanceerd design? (Vuistregel: \(n_{\max}/n_{\min}\).)
a) \(25\) per cel \(\times 6\) cellen \(= 150\) proefpersonen voor power \(\approx .80\) bij middelgroot effect. Bij \(30\) per cel: \(180\).
b) Met \(40\) totaal in een \(2 \times 3\): gemiddeld \(\sim 7\) per cel. Power voor middelgrote effecten zal richting \(.20\)–\(.30\) gaan. Opties:
c) Een vuistregel: bij \(n_{\max}/n_{\min} > 4\) (de grootste cel meer dan vier keer de kleinste) wordt de \(F\)-toets onbetrouwbaar als óók Levene significant is. Ook bij \(n_{\max}/n_{\min} \approx 1.5\) is voorzichtigheid gepast. Tabachnick en Fidell adviseren: balanceer waar je kunt, en als het echt niet kan, gebruik Type III en rapporteer expliciet.
oneway.test() in R doet standaard een Welch-ANOVA — een variant van eenweg-ANOVA die niet aanneemt dat varianties gelijk zijn.
a) In welke twee situaties is Welch beter dan de gewone (Student-)ANOVA?
b) Wat verandert er aan de vrijheidsgraden \(df_w\)? (Hint: niet meer \(N - k\).)
c) Voor de mol-één-weg-vraag (alleen gang_diepte, twee groepen): bij twee groepen is Welch-ANOVA equivalent aan een bekende \(t\)-toets-variant. Welke?
a) Welch is robuust tegen (1) ongelijke varianties en (2) ongelijke groepsgroottes. Beide tegelijk is precies waarbij gewone ANOVA misgaat.
b) Welch berekent een gecorrigeerd \(df_w\) via de Welch-Satterthwaite-formule, dat lager is dan \(N - k\) en niet noodzakelijk een geheel getal. Bij ernstige varianti-ongelijkheid kan \(df_w\) flink dalen — wat conservativer toetst.
c) Bij \(k = 2\) is Welch-ANOVA equivalent aan Welch’s \(t\)-toets (de “ongelijke variantie”-versie van de Student-\(t\)). De Student-\(t\) met var.equal = TRUE en gewone eenweg-ANOVA met \(k = 2\) geven identieke \(p\)-waarden — en dat is gekwadrateerd in \(F\).
Stel je doet een \(2 \times 2\) tweeweg-ANOVA en de interactie is niet significant (\(F(1, 100) = 1.20\), \(p = .270\)). Hoofdeffect \(A\) wel (\(F(1, 100) = 18.50\), \(p < .001\)). Hoofdeffect \(B\) ook (\(F(1, 100) = 9.40\), \(p = .003\)).
a) Mag je nu de hoofdeffecten los interpreteren? Waarom (niet)?
b) Een collega zegt: “Haal de interactie uit het model — een additief model is parsimonieuser en geeft scherpere hoofdeffecten.” Wanneer is dat advies gepast, en wanneer niet?
c) Een andere collega zegt: “Interactie weglaten is gevaarlijk; je verstopt mogelijke kruis-effecten.” Heeft hij gelijk?
a) Ja. Bij niet-significante interactie zijn de hoofdeffecten zinvol los te interpreteren: het effect van \(A\) is voor alle \(B\)-niveaus statistisch hetzelfde (binnen toeval), en omgekeerd. Je rapporteert ze als “het effect van \(A\) was X, ongeacht \(B\).” Dit is precies het tegenovergestelde van een significante interactie.
b) Soms gepast: in een a-priori additief model (theorie zegt dat \(A\) en \(B\) onafhankelijk werken) en bij grote sample-size kan het rapporteren van een additief model schoner zijn. Het werkt minder goed bij kleine samples — een niet-significante interactie betekent niet dat hij afwezig is, alleen dat we hem niet konden detecteren. Bij twijfel: laat hem in het model maar rapporteer: “interactie niet significant; hoofdeffecten konden onafhankelijk worden geïnterpreteerd.”
c) Ten dele. Je verstopt potentieel niets als je de interactie hebt getoetst en niet-significant vond, en rapporteert dat. Je verstopt wel iets als je hem stilletjes weglaat zonder te tonen dat hij niet significant was. Voor de lezer is altijd wenselijk: rapporteer de hele tabel inclusief interactie, en vertel pas in de tekst dat je de hoofdeffecten los interpreteert.
Alle commando’s op één plek
# Pakketten activeren — één keer per sessie.
library(car) # voor Anova(., type = 3) en leveneTest()
library(lsr) # voor etaSquared()
library(emmeans) # voor estimated marginal means en EMM-uitwerking per niveau
# Type III SS vereist sum-to-zero contrasten op factoren.
options(contrasts = c("contr.sum", "contr.poly"))
# Eerste dataset laden — object 'mol_werkdag' verschijnt vanzelf.
load("data/mol_werkdag.RData")
# Tweede dataset laden — object 'fluittoon_reactie' verschijnt vanzelf.
load("data/fluittoon_reactie.RData")
# Snelle inkijk: structuur, factor-niveaus, eerste waarden.
str(mol_werkdag)# Aantallen per groep (eenweg) of per cel (factorial).
table(mol_werkdag$gang_diepte)
table(mol_werkdag$gang_diepte, mol_werkdag$dagdeel)
# Groepsgemiddelden / cel-gemiddelden van Y.
aggregate(graaftempo ~ gang_diepte, data = mol_werkdag, FUN = mean)
aggregate(graaftempo ~ gang_diepte * dagdeel, data = mol_werkdag, FUN = mean)# Levene's toets voor gelijke varianties.
leveneTest(graaftempo ~ gang_diepte, data = mol_werkdag)
leveneTest(graaftempo ~ gang_diepte * dagdeel, data = mol_werkdag)
# Residual plot van het lm-object.
plot(lm(graaftempo ~ gang_diepte, data = mol_werkdag), which = 1)# Welch-ANOVA — robuust tegen ongelijke varianties.
oneway.test(graaftempo ~ gang_diepte, data = mol_werkdag, var.equal = FALSE)
# Brown-Forsythe — alternatieve robuuste F-toets bij ongelijke varianties.
# install.packages("onewaytests") # eerste keer
onewaytests::bf.test(graaftempo ~ gang_diepte, data = mol_werkdag)
# Kruskal-Wallis — non-parametrisch alternatief bij kleine n + non-normaliteit.
kruskal.test(graaftempo ~ gang_diepte, data = mol_werkdag)# Eenweg-ANOVA via lm() + Anova().
m1 <- lm(graaftempo ~ gang_diepte, data = mol_werkdag)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
# Defensief — herhaal vlak voor elke Type-III-aanroep voor het geval dat
# een eerder blok of package-load de instelling heeft veranderd.
options(contrasts = c("contr.sum", "contr.poly"))
Anova(m1, type = 3)
etaSquared(m1)
# Tweeweg-ANOVA — let op de * voor hoofdeffecten + interactie.
# A = gang_diepte (rij-factor), B = dagdeel (kolom-factor).
m2 <- lm(graaftempo ~ gang_diepte * dagdeel, data = mol_werkdag)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
options(contrasts = c("contr.sum", "contr.poly"))
Anova(m2, type = 3)
etaSquared(m2)interaction.plot(
x.factor = mol_werkdag$dagdeel,
trace.factor = mol_werkdag$gang_diepte,
response = mol_werkdag$graaftempo,
ylab = "Gemiddeld graaftempo (cm/uur)",
xlab = "Dagdeel",
trace.label = "gang_diepte"
)# Tukey HSD: alleen na significant hoofdeffect, op aov-object.
TukeyHSD(aov(reactietijd ~ wind, data = fluittoon_reactie))
# Estimated marginal means + paarsgewijze vergelijkingen per cel.
emm <- emmeans(m2, ~ dagdeel | gang_diepte)
emm
pairs(emm)
# Bonferroni-correctie via emmeans
pairs(emmeans(m_wind, ~ wind), adjust = "bonferroni")
# Holm-Bonferroni
pairs(emmeans(m_wind, ~ wind), adjust = "holm")
# Scheffé
pairs(emmeans(m_wind, ~ wind), adjust = "scheffe")
# Generiek: p-vector aanpassen
p.adjust(p_vector, method = "bonferroni")
p.adjust(p_vector, method = "holm")
p.adjust(p_vector, method = "BH") # Benjamini-Hochberg, FDR-controle# Gerichte contrasten via emmeans + contrast().
m_wind <- lm(reactietijd ~ wind, data = fluittoon_reactie)
emm <- emmeans(m_wind, ~ wind)
contrast(emm, list(
wind_vs_geen = c(-2, 1, 1) / 2,
licht_vs_hard = c( 0, -1, 1)
))Even oefenen op tentamen-toon
Onderaan dit hoofdstuk staan vier korte theorievragen in tentamen-stijl en één mini R-opdracht. Geen Tellegen-frame meer — student-aan-tentamen-modus. Hou de vuistregels uit dit hoofdstuk paraat: \(\eta^2 \geq .14\) groot, \(.06\) middel, \(.01\) klein; Levene \(p < .05\) flagt; Tukey alleen na significant hoofdeffect.
Een bever bestudeert een vergelijking van vier soorten dam-bouwhout (eik / wilg / berk / esdoorn) op damstevigheid (kN/m). Hij rapporteert dit fragment uit een eenweg-ANOVA-tabel.
R
Anova Table (Type III tests)
Response: stevigheid
Sum Sq Df F value Pr(>F)
houtsoort 84.6 3 7.230 0.0002 ***
Residuals 312.0 80Hoeveel proefdammen (\(N\)) hadden in totaal aan deze studie meegedaan?
c) \(84\). \(df_w = N - k = 80\), met \(k = 4\) groepen, dus \(N = 80 + 4 = 84\).
In een onderzoek bij de uilenwacht naar het effect van lichtsterkte (3 niveaus) op alertheid lezen we: \(F(2, 87) = 6.40\), \(p = .003\), \(\eta^2 = .128\).
Welke conclusie volgt het beste uit deze cijfers?
b) Eta-squared van \(.128\) ligt tussen middel (\(.06\)) en groot (\(.14\)); \(\eta^2\) is letterlijk de proportie variantie verklaard, dus \(\sim 13\%\). Optie a overdrijft, c miskent het significante \(p\), d verwart \(\eta^2\) met VIF.
Bij de kerkuilen worden \(80\) uilen onderzocht in een \(2 \times 2\) design: veld (open / bos) en maan (volle / nieuwe). Uitkomst: aantal succesvolle vangsten per nacht. De interactieplot ziet er zo uit:
Welke uitspraak past het best bij dit patroon?
c) In het open veld is het maan-effect klein (\(4.5 \to 4.7\), +\(0.2\)). In het bos juist sterk negatief bij nieuwe maan (\(6.8 \to 3.1\), −\(3.7\)). Lijnen niet parallel \(\Rightarrow\) duidelijke interactie. De marginalen suggereren ook hoofdeffecten, maar bij een sterke interactie zoals deze interpreteer je ze niet los.
Een ecoloog doet een eenweg-ANOVA met vier vegetatie-types als factor en soortenrijkdom als uitkomst. Hij rapporteert: \(F(3, 76) = 1.81\), \(p = .15\), \(\eta^2 = .067\).
Hij overweegt of hij Tukey HSD nu mag draaien om alsnog te zien welke vegetatie-types onderling verschillen.
Wat is het juiste advies?
c) De gangbare conventie: post-hoc-toetsen alleen ná een significant ANOVA-hoofdeffect. Met \(p = .15\) is er geen toetsbasis om verder paren te onderzoeken — dat zou neerkomen op \(p\)-hacking via een achterdeurtje. Een matige \(\eta^2\) kan in theorie wijzen op een te kleine \(N\) (lage power); het juiste antwoord is dan repliceren met meer dieren, niet de huidige data verder uitpluizen.
De dassen-burcht
Een ecoloog onderzoekt hoe lang dassen ’s nachts foerageren — hun foerageertijd (in minuten) — onder twee factoren: bodem (zand / klei / leem; drie niveaus) en seizoen (zomer / winter; twee niveaus). Bij \(78\) dassen is per dier de gemiddelde foerageertijd vastgelegd, samen met het bodemtype waarop ze hun burcht hebben en het seizoen van meting. De dataset staat in data/dassen_burcht.RData en bevat het object dassen_burcht. Sla je R-commando’s op in één scriptbestand: dassen.R. Gebruik \(\alpha = .05\).
a) Is de aanname van gelijke varianties houdbaar? Rapporteer Levene.
b) Voer een \(3 \times 2\) tweeweg-ANOVA uit (Type III). Welke effecten zijn significant? Rapporteer de toetsstatistiek, de vrijheidsgraden en de p-waarde per effect.
c) Welk effect heeft de grootste \(\eta^2\)? Wat is de inhoudelijke grootte (klein / middel / groot)?
d) Maak een interaction.plot() met bodem op de x-as en lijnen voor seizoen. Beschrijf in één zin het patroon.
load("data/dassen_burcht.RData")
leveneTest(foerageertijd ~ bodem * seizoen, data = dassen_burcht)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 0.6783 0.6413
72 m_das <- lm(foerageertijd ~ bodem * seizoen, data = dassen_burcht)
# Contrasts zoals in de uni-oplossing (Type-III toetsen met sum-to-zero).
options(contrasts = c("contr.sum", "contr.poly"))
Anova(m_das, type = 3)Anova Table (Type III tests)
Response: foerageertijd
Sum Sq Df F value Pr(>F)
(Intercept) 1083733 1 4288.548 < 2.2e-16 ***
bodem 13676 2 27.059 0.000000001722 ***
seizoen 5496 1 21.749 0.000013946484 ***
bodem:seizoen 2238 2 4.428 0.01536 *
Residuals 18195 72
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1etaSquared(m_das) eta.sq eta.sq.part
bodem 0.34377376 0.4307134
seizoen 0.13238695 0.2256226
bodem:seizoen 0.05588776 0.1095272# Tijdelijke factoren met hoofdletter-niveaus voor leesbare as- en
# legend-labels — onderliggende data blijft ongewijzigd.
.bodem_lbl <- factor(tools::toTitleCase(as.character(dassen_burcht$bodem)),
levels = tools::toTitleCase(levels(dassen_burcht$bodem)))
.seizoen_lbl <- factor(tools::toTitleCase(as.character(dassen_burcht$seizoen)),
levels = tools::toTitleCase(levels(dassen_burcht$seizoen)))
interaction.plot(
x.factor = .bodem_lbl,
trace.factor = .seizoen_lbl,
response = dassen_burcht$foerageertijd,
ylab = "Gemiddelde foerageertijd (min)",
xlab = "Bodem",
trace.label = "Seizoen"
)
Het kruispunt voor het dassen-design (3 × 2):
apa_marginals_table(
m_das,
factor_A = "bodem",
factor_B = "seizoen",
data = dassen_burcht,
dependent_label = "foerageertijd (min)",
factor_A_label = "Bodem",
factor_B_label = "Seizoen",
caption = "Dassen — foerageertijd per cel met EMM-marginalen"
)| Bodem |
Seizoen
|
Marginaal B | ||
|---|---|---|---|---|
| Zomer | Δ Zomer − Winter | Winter | ||
| Zand | 96.14 (17.88, 14) | -8.30 | 104.44 (16.54, 13) | 100.29 |
| Klei | 116.68 (17.48, 13) | -10.05 | 126.72 (16.33, 13) | 121.70 |
| Leem | 115.94 (14.01, 12) | -32.07 | 148.02 (12.01, 13) | 131.98 |
| Marginaal A | 109.59 | — | 126.39 | 117.99 |
| Note. Dassen — foerageertijd per cel met EMM-marginalen. Cellen tonen M (SD, n) op basis van foerageertijd (min); rij- en kolom-marginalen zijn ongewogen estimated marginal means (EMM); rechtsonder de grand mean. Hoofdeffect Bodem: F(2, 72) = 27.06, p < .001 — H0: alle rij-marginalen gelijk (μA⋅’s gelijk). Hoofdeffect Seizoen: F(1, 72) = 21.75, p < .001 — H0: alle kolom-marginalen gelijk (μ⋅B’s gelijk). Interactie Bodem × Seizoen: F(2, 72) = 4.43, p .015 — H0: cel-patronen additief (geen interactie). | ||||
profile_plot_with_marginals(
m_das,
factor_A = "bodem",
factor_B = "seizoen",
data = dassen_burcht,
dependent_symbol = "foerageertijd",
dependent_label = "Foerageertijd (min)",
factor_A_label = "Bodem",
factor_B_label = "Seizoen"
)
a) De aanname van gelijke varianties is houdbaar (Levene niet significant — zie output).
b) Alle drie de effecten zijn significant: hoofdeffect bodem \(F(2, 72) = 27.06\), \(p < .001\); hoofdeffect seizoen \(F(1, 72) = 21.75\), \(p < .001\); interactie bodem \(\times\) seizoen \(F(2, 72) = 4.43\), \(p = .015\).
c) \(\eta^2_{\text{bodem}} = .344\) is de grootste en valt ver in de “groot”-categorie (\(\geq .14\)). Seizoen is ook groot (\(\eta^2 = .132\)); interactie middelgroot (\(\eta^2 = .056\)).
d) Op zandbodem is het verschil tussen zomer en winter klein; op leem is het verschil substantieel — dassen op leem foerageren in de winter veel langer dan in de zomer. Dat is de interactie die de \(F\)-toets bevestigt.
Een \(3 \times 2\) tweeweg-ANOVA toonde significante hoofdeffecten van bodem, \(F(2, 72) = 27.06\), \(p < .001\), \(\eta^2 = .344\), en seizoen, \(F(1, 72) = 21.75\), \(p < .001\), \(\eta^2 = .132\). De interactie tussen bodem en seizoen was eveneens significant, \(F(2, 72) = 4.43\), \(p = .015\), \(\eta^2 = .056\). De aanname van gelijke varianties was houdbaar volgens Levene’s toets.
Onderwerpen die buiten dit thema blijven
ANOVA is een groot huis met veel zijdeuren. Wat dit thema je gaf is de hoofdingang — eenweg, tweeweg, interactie, post-hoc, contrasten. Wat hierachter ligt en in dit werkboek niet behandeld wordt:
WRS2::t1way en verwante procedures op basis van M-estimators of trimmed means. De Welch-ANOVA in T9 is hiervan een lichte variant; de volledige robust-toolbox blijft hier buiten.lme4::lmer) komt nergens uitgewerkt voor.BayesFactor::anovaBF is een toegankelijke ingang.pwr::pwr.anova.test voor sample-size-planning. Belangrijk in onderzoeksmethoden, hier niet behandeld.Voor verdieping: Tabachnick & Fidell (2007); Wilcox (2017, robust methods); Maxwell, Delaney & Kelley (2017, design and analysis).
Toen de zon onderging, schoof de mol zijn aantekeningen onder een wortel.
“En, mol?” vroeg de olifant, die nog even was blijven staan.
“Diepte maakt veel uit,” zei de mol, “en dagdeel ook iets. En het diepte-effect is in de avond groter dan in de ochtend.” Hij dacht even na. “Ik had die laatste niet gemerkt als ik alleen naar gang-diepte had gekeken.”
“En ik,” zei de olifant, “weet nu dat de wind er flink toe doet, en dat het bij sommige wind-condities anders werkt voor mannetjes dan voor vrouwtjes.”
Ze zaten samen. Het werd langzaam donker. Ergens in de verte ritselde een blad.
Dit werkboek is geschreven voor studenten die variantieanalyse leren via R. De didactische lijnen volgen de gangbare opbouw van Nederlandse universitaire MVDA-cursussen; alle voorbeelden, datasets, vragen en formuleringen in dit hoofdstuk zijn origineel. De gebruikte drempels en formules zijn standaard statistische conventies; verwijzingen naar Cohen (1988) en Tabachnick & Fidell (2007) volgen de gebruikelijke citaatpraktijk in dit veld.
Versie: July 2026 — CountCamp Lab