3. Covariantieanalyse (ANCOVA)

Het verhaal van dit thema

De bever zat aan de bosrand, een schrift open op zijn knieën, en keek peinzend de heuvel af. Negentig dammen had hij in beeld, negentig verschillende bevers — elk hun eigen werkstuk, elk hun eigen verhaal eronder. “Klaar,” mompelde een bever in de verte tegen niemand in het bijzonder. “Tevreden? Dat is een ándere vraag.”

Sommige bevers (zo wist hij van horen knagen) waren perfectionist: alles moest klopjes-precies, en zelfs dan was er bijna altijd iets dat knaagde. Andere werkten als pragmaticus: goed-genoeg-is-goed, een gerust streepje eronder. En sommigen, zelden maar onmiskenbaar, gingen impulsief te werk — een paar van die dammen waren hun trots, andere hun diepste ergernis. Bij elke bever noteerde hij óók hoeveel familieleden er tijdens het bouwen hadden meegeholpen.

“Want,” zei hij tegen niemand in het bijzonder, “tevredenheid is geen losse score. Ze hangt er ook vanaf hoeveel hulp je kreeg. Als de hele familie meebouwt, is het sowieso prettiger werk dan in je eentje. Onthoud dat.”

De bever wilde dus iets weten dat de mol nog niet had hoeven weten: het effect van een categorische factor (bouwwoede), gecorrigeerd voor een continue voorspeller (meedoende_familie). Niet apart toetsen, maar samen — in één model. Wat hij berekent heet analysis of covariance ANCOVA.

Bij ANCOVA stel je dezelfde drie vragen aan je gegevens als bij ANOVA, maar één vraag krijgt een tweekoppige draai:

1. Voorspelt het samen iets? — de $F$-toets op het hele model, met covariaat én factor.
2. Wie draagt bij — de covariaat, de factor, of beide? — twee aparte $F$-toetsen, een voor elk effect, en bij significantie van de factor: paarsgewijze vergelijking van de adjusted means estimated marginal means / EMM tussen niveaus.
3. Klopt het allemaal wel? — de aannames van ANOVA (normaliteit, homoscedasticiteit, onafhankelijke observaties) plus de drie specifieke aan ANCOVA: een covariaat die zonder meetfout gemeten is, lineariteit van de covariaat-Y-relatie en homogeniteit van regressiehellingen homogeneity of regression slopes — kortweg parallellisme parallelism: de slopes van Y op de covariaat horen ongeveer gelijk te zijn over groepen.

TipWat je na dit thema moet kunnen

Na dit thema kun je:

1. herkennen wanneer ANCOVA de juiste techniek is, en nagaan of een variabele een geschikte covariaat is — ook door te toetsen of de groepen op de covariaat verschillen;
2. uitleggen waarom je een covariaat toevoegt: om de foutvariantie te verkleinen, en om een verschil op de covariaat tussen de groepen weg te corrigeren;
3. de aannames nagaan — de drie van ANOVA en de drie van ANCOVA zelf (een covariaat zonder meetfout, lineariteit, parallellisme) — met parallellisme als de aanname die je toetst via de interactie factor $\times$ covariaat;
4. een ANCOVA draaien en lezen, de helling binnen de groepen ($b_w$) aflezen en interpreteren, en zeggen wat er aan de vrijheidsgraden verandert als er een covariaat bij komt;
5. gecorrigeerde gemiddelden uitrekenen met de formule, ze aflezen in een plot, ze tussen de groepen vergelijken, en zien of de correctie een verschil kleiner, groter of omgekeerd maakt;
6. zeggen waar ANCOVA ophoudt: bij hellingen die niet parallel lopen is het niet het juiste model, en bij groepen die niet door loting zijn gevormd beantwoordt de correctie niet vanzelf een oorzaak-vraag.

In R doe je dat met options(contrasts = …) (en aan het eind weer terug), lm(), Anova(., type = 3), etaSquared(., type = 3), summary() (voor $b_w$), en emmeans_test() met get_emmeans() uit het pakket rstatix voor de gecorrigeerde gemiddelden en de paarsgewijze vergelijking met Bonferroni-correctie.

OpmerkingA. Wat is ANCOVA?

ANCOVA combineert ANOVA en regressie. Je hebt één factor (categorisch, zoals ANOVA), één covariaat covariate (continu, zoals MRA), en één afhankelijke variabele $Y$. Twee doelen tegelijk: (1) foutreductie — de covariaat trekt een deel van de Y-variantie weg, het residu wordt kleiner, en de $F$-toets op de factor krachtiger (zolang de groepen niet ook op de covariaat verschillen; doen ze dat wel, dan verandert ook het effect van de factor zelf, en dat kan twee kanten op — zie T4). (2) Bias-correctie — als groepen op de covariaat van elkaar verschillen, corrigeer je daarvoor zodat het pure factor-effect overblijft.

OpmerkingB. Adjusted means versus raw means

Bij ANCOVA verschuift je interpretatie van “wat is het gemiddelde Y per groep?” naar “wat zou het gemiddelde Y per groep zijn als alle groepen dezelfde covariaat-waarde hadden?”. Die geconditioneerde gemiddelden heten adjusted means of estimated marginal means (EMM):

$$\bar{Y}_j^* = \bar{Y}_j - b_w(\bar{C}_j - \bar{C})$$

waar $b_w$ de pooled-within slope is — de helling van Y op de covariaat binnen de groepen, over de groepen samengenomen. ANCOVA toetst altijd de adjusted means, niet de raw means.

OpmerkingVoor wie verder wil — centreren, een techniek met een naam

In de ANCOVA-formule trek je de grand mean van de covariaat $\bar{C}$ af van de groeps-covariaat-mean $\bar{C}_j$. Dat is centreren centering / mean-centering — een gangbare regressie-techniek waarbij je de covariaat verschuift zodat zijn gemiddelde nul wordt: $C^c_i = C_i - \bar{C}$.

Twee voordelen:

1. Interpretatie van het intercept — zonder centreren is $b_0$ de voorspelde Y-waarde bij covariaat $= 0$, en die nul valt vaak buiten de data (een tak van 0 cm bestaat niet als steiger). Met centreren is $b_0$ de voorspelde Y-waarde bij de gemiddelde covariaat-waarde — meestal midden in de data, dus inhoudelijk uit te leggen.
2. Minder multicollineariteit bij interactie-modellen ($X_1 \cdot X_2$) of polynomen ($X^2$). Een ongecentreerde $X$ correleert sterk met $X^2$ (voor de hele getallen $X = 3, 4, \ldots, 12$ is die correlatie $.99$); dat geeft hoge VIF-waarden en onstabiele coëfficiënten. Een gecentreerde $X^c$ correleert veel zwakker met $(X^c)^2$ (bij een symmetrische verdeling: nul). Dezelfde logica geldt bij interactie-termen: $X_1 \cdot X_2$ correleert sterk met $X_1$ en $X_2$ als die ongecentreerd zijn.

Wat centreren níét verandert. In het gewone ANCOVA-model (het additieve model, zonder interactie) verandert centreren de schatting van $b_0$ (en alleen die); de $F$-toetsen op de andere termen blijven gelijk. In een model mét interactie factor $\times$ covariaat, zoals het model waarmee je verderop parallellisme toetst, ligt dat anders. Daar meet de toets op de factor het groepsverschil in $Y$ bij covariaat $= 0$; centreren verplaatst die nul naar het gemiddelde, en dus verschuift ook die toets. Voor ANCOVA is centreren een didactische lens — het helpt je begrijpen waarom adjusted means zijn wat ze zijn — meer dan een verplichte stap.

OpmerkingHoe dit hoofdstuk leest

Dit hoofdstuk volgt dezelfde drie-laagse structuur als thema 1 en 2:

1. Verhaal-frame — wat de bever doet of denkt.
2. Algemene vorm — abstract, statistiek-Latijn met Y, factor, covariaat, mijn_data.
3. Voor onze dieren — uitvoerend, met # hekjes en concrete namen.

Plus de antwoorden in een dichtgeklapt kader: die open je pas ná je eigen poging. T-vragen vlechten zich tussen de hoofdvragen door om je rekenintuïtie te oefenen.

Werkmaterialen — R-pakketten en functies

Wat de bever uit zijn gereedschapskist pakt

Gereedschap	Waar het voor is	Pakket
lm()	Lineair model — voor zowel ANOVA, ANCOVA als regressie	base
Anova(model, type = 3)	ANCOVA-tabel met type III sums of squares	car
etaSquared(model, type = 3)	Effectgrootte $\eta^2$ en $\eta^2_p$	lsr
summary(model)	Coëfficiënten, waaronder de helling $b_w$ van de covariaat	base
aggregate()	Groepsgemiddelden en standaardafwijkingen per groep	base
emmeans_test(..., covariate = , p.adjust.method = "bonferroni")	Paarsgewijze vergelijking van de adjusted means, met Bonferroni-correctie	rstatix
get_emmeans()	De adjusted means (EMM) per niveau, uit de uitkomst van emmeans_test()	rstatix
emmeans(model, ~ factor)	Adjusted means (EMM) per niveau, rechtstreeks uit het model	emmeans
cooks.distance()	Diagnostische waarde per observatie	base
apa_ancova_table()	Adjusted-means-tabel — hulpfunctie van dit oefenboek; zonder die hulpfunctie haal je dezelfde getallen uit aggregate() en emmeans()	eigen

OpmerkingWaarom Type III in ANCOVA?

Type III sums of squares geven elke term zijn unieke bijdrage na correctie voor de andere termen. Dat is precies wat je wilt: de factor toetsen ná aftrek van de covariaat (en omgekeerd). Bij Type I (sequentieel) zou de volgorde in de formule het resultaat veranderen; bij Type III niet. De regel options(contrasts = c("contr.sum", "contr.poly")) bovenaan zorgt dat Anova(., type = 3) de juiste antwoorden geeft.

OpmerkingTechnisch detail: 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. Wat er gebeurt als je hem vergeet, hangt af van het model:

• In het ANCOVA-model zelf (factor + covariaat, zonder interactie) zijn de toetsen dezelfde: $F$, $p$, $\eta^2_p$ en de EMM veranderen niet. Alleen de coëfficiënten van de factor in summary(lm) zien er anders uit.
• In het model met interactie (factor * covariaat, in 3.2) veranderen de toetsen op de hoofdeffecten wél. Maar in dat model lees je alleen de rij van de interactie, en die blijft gelijk.

Onthoud: zet deze regel bovenaan je script, vóór je eerste lm(). R legt de codering vast op het moment dat het model geschat wordt.

OpmerkingVoor wie verder wil — contr.sum-coëfficiënten lezen in summary(lm)

In thema 1 (MRA, in de sectie Dummy-codering — een categorische voorspeller in je regressie) zag je dummy-codering in zijn basisvorm — R’s default contr.treatment, met één groep als referentie (de $0/0$-rij) en de andere groepen als $0/1$-dummies die het verschil vanaf die referentie meten. In ANCOVA gebruiken we een variant op die codering: contr.sum, ook wel sum-to-zero coding. Geen drama — de adjusted means uit emmeans() blijven identiek aan wat je onder contr.treatment zou krijgen, en wat je inhoudelijk afleest is hetzelfde. Maar de getallen in summary(m_ancova) zien er anders uit, en dat moet je herkennen anders gaat de hand-rekening mis.

De twist. Bij contr.sum krijgt de “extra” groep (die in contr.treatment de $0/0$-referentie was) een $-1/-1$-rij in plaats van $0/0$. Dat zorgt ervoor dat de drie groep-deviaties netjes samen optellen tot nul — handig voor Anova(., type = 3) bij ongebalanceerde designs. Stel we kiezen pragmaticus als $d_1$ en impulsief als $d_2$ — dan wordt perfectionist de derde rij (geen referentie meer in de treatment-zin, wel impliciet via $-(d_1 + d_2)$):

	$d_1$	$d_2$
pragmaticus	1	0
impulsief	0	1
perfectionist	$-1$	$-1$

R noemt de kolommen bouwwoede1 en bouwwoede2. Perfectionist krijgt geen eigen kolom — die wordt afgeleid via $-(d_1 + d_2)$.

De regressievergelijking heeft formeel dezelfde vorm als bij de dummy-codering in thema 1:

$$\hat{Y} = b_0 + b_1 \cdot d_1 + b_2 \cdot d_2 + b_w \cdot C.$$

Maar de interpretatie van $b_0$, $b_1$ en $b_2$ kantelt — en juist daar zit waarom je contr.sum-output niet kunt lezen alsof het contr.treatment is:

• (Intercept) ($b_0$) — niet de adjusted mean van een referentiegroep (zoals bij de dummy-codering in thema 1, met contr.treatment), maar het grand-mean-niveau over de drie groepen bij $C = 0$.
• bouwwoede1 ($b_1$) — deviatie van pragmaticus t.o.v. dat grand-mean-niveau (geen verschil-met-referentie meer).
• bouwwoede2 ($b_2$) — deviatie van impulsief t.o.v. dat grand-mean-niveau.
• De derde groep (perfectionist) heeft geen eigen parameter; haar deviatie is $-(b_1 + b_2)$.

Het advies Codeer lekker zelf uit thema 1 (bij de dummy-codering) geldt hier nog steeds: zet je factor-levels expliciet met factor(..., levels = c(...)) vóór je het model fit, anders ordent R ze alfabetisch en weet je niet zeker welke groep welke kolom krijgt. Dat geldt voor contr.treatment én contr.sum — alleen de betekenis van de coefficienten verschilt, niet de noodzaak om zelf de volgorde te kiezen.

Adjusted means herrekenen uit deze coefficienten — bij grand mean covariaat $\bar{C}$:

$$\bar{Y}_j^* = b_0 + (\text{deviatie}_j) + b_w \cdot \bar{C}.$$

Groep	Deviatie
pragmaticus	$b_{\text{bouwwoede1}}$
impulsief	$b_{\text{bouwwoede2}}$
perfectionist	$-(b_{\text{bouwwoede1}} + b_{\text{bouwwoede2}})$

Minisommetje — fictieve getallen. Stel summary(m_ancova) rapporteert (let op: deze cijfers zijn verzonnen om de mechaniek te zien — niet de echte bever-output):

• (Intercept) = $65$
• bouwwoede1 = $+3$
• bouwwoede2 = $+5$
• meedoende_familie ($b_w$) = $1.5$

Plus: grand mean covariaat $\bar{C} = 4$. Perfectionist-deviatie staat niet in de output — die leid je zelf af:

$$-(b_{\text{bouwwoede1}} + b_{\text{bouwwoede2}}) = -(3 + 5) = -8.$$

Adjusted means bij $C = \bar{C} = 4$:

• $\bar{Y}_{\text{pragmaticus}}^* = 65 + 3 + 1.5 \cdot 4 = 74$
• $\bar{Y}_{\text{impulsief}}^* = 65 + 5 + 1.5 \cdot 4 = 76$
• $\bar{Y}_{\text{perfectionist}}^* = 65 + (-8) + 1.5 \cdot 4 = 63$

Check — som van de drie deviaties: $3 + 5 + (-8) = 0$ (kenmerk van contr.sum). Ongewogen gemiddelde adjusted means: $(74 + 76 + 63)/3 = 71$ = $b_0 + b_w \cdot \bar{C} = 65 + 6$. Klopt — dáár staat het grand-mean-niveau in (Intercept) voor, na correctie voor de covariaat-effect tot grand mean.

Als je dezelfde drie adjusted means uit emmeans(m_ancova, ~ bouwwoede) zou krijgen — en dat doe je — dan is de uitleg compleet. Als ze niet matchen, las je summary() als contr.treatment terwijl contr.sum aanstaat. Klassieke val: emmeans() doet niks raars, jij las de coefficienten verkeerd.

OpmerkingVoor wie verder wil — AN-mono-COVA, AN-bi-COVA, …

Voor wie chemisch wil denken. Geen officiële terminologie — wel een geheugenanker.

Een ANCOVA met één factor + één covariaat zou je chemisch kunnen noemen: AN-mono-COVA (één covariaat = mono).

• AN-mono-COVA — $1$ factor $+ 1$ covariaat (de standaard-ANCOVA in dit thema).
• AN-bi-COVA — $1$ factor $+ 2$ covariaten.
• AN-tri-COVA — $1$ factor $+ 3$ covariaten.
• (M)AN-poly-COVA — meer factoren én meer covariaten; opgeschaald naar (M)ANCOVA in thema 5.

De “mono” / “bi” / “tri” telt het aantal covariaten. Geen literatuur zal je dit ooit zo zien noemen — maar voor het zelf-onthouden van wat ANCOVA structureel doet (factor mét covariaat) is het een aardige denkframe. Wie de structuur van een vraag herkent, weet meteen welke functie-aanroep daarbij hoort.

Notatie — symbolen voor ANCOVA

OpmerkingSleutelsymbolen in dit hoofdstuk

In dit oefenboek zie je telkens dezelfde notatie:

• $\bar{Y}_j$ — raw steekproef-gemiddelde van $Y$ in groep $j$ (gewoon mean() op de Y-waarden in die groep).
• $\bar{C}_j$ — gemiddelde van de covariaat $C$ in groep $j$.
• $\bar{C}$ — grand mean van de covariaat over alle observaties.
• $\bar{Y}_j^*$ — adjusted mean EMM: het gemiddelde van Y in groep $j$ na correctie voor de covariaat, alsof iedereen dezelfde covariaat-waarde $\bar{C}$ had: $$\bar{Y}_j^* = \bar{Y}_j - b_w (\bar{C}_j - \bar{C}).$$
• $b_w$ — pooled-within regressie-slope. Niet de gewone OLS-slope op alle data, maar het gewogen gemiddelde van de slopes binnen elke groep. Past bij de ANCOVA-aanname dat die binnen-groep-slopes gelijk zijn (parallellisme).
• $F$, $df$ — toetsstatistiek en vrijheidsgraden voor zowel covariaat als factor.
• $\eta^2_p$ partial eta-squared — effectgrootte, proportie variantie die het effect verklaart na aftrek van de andere termen in het model.

Voor inferentie praat je over de populatie-waarden ($\mu_j$, en de helling binnen de groepen in de populatie); voor rapportage over de steekproef ($\bar{Y}_j$, $\bar{Y}_j^*$, $b_w$, $M$, $SD$).

Let op β. In dit oefenboek staat β voor het gestandaardiseerde regressiegewicht, zoals in thema 1: de helling als $X$ en $Y$ allebei in standaardafwijkingen gemeten zijn. Het is hier dus geen teken voor een helling in de populatie. In ANCOVA werk je met de gewone, ongestandaardiseerde helling $b_w$, in de eenheden van $Y$ per eenheid van de covariaat.

TipVuistregels zijn afspraken, geen wetten

De drempels in dit oefenboek zijn breed gangbare conventies — niet universeel:

• robuustheid van $F$: dezelfde drie vuistregels als bij ANOVA — $n \geq 15$ per groep, $n_{\max}/n_{\min} \leq 1.5$, en $s_{\max}/s_{\min} < 2$ voor ongeveer gelijke varianties (zie Robuustheid van $F$ in 3.1.b);
• parallellisme-aanname: de interactie tussen factor en covariaat moet niet significant zijn ($p > .05$). Bij kleine samples is “niet significant” zwak bewijs voor gelijkheid, want dan mist de toets veel; check ook visueel of de slopes per groep ongeveer parallel lopen;
• lineariteit covariaat-Y: per groep zou de relatie ongeveer recht moeten zijn — check met een scatter Y vs $C$ per groep;
• geen meetfout in de covariaat: ANCOVA neemt aan dat $C$ zonder ruis gemeten is; bij sterke meetfout wordt $b_w$ onderschat (meetfout maakt de geschatte helling platter) en de bias-correctie onbetrouwbaar (zie T7);
• Cook’s distance (uit thema 1): boven $1$ is een observatie invloedrijk genoeg om je zorgen over te maken.

Handboeken kiezen hier iets andere getallen. Dit oefenboek volgt consequent één set, maar dat is een keuze, geen natuurwet: werk je met ander materiaal, kijk dan welke drempels dat aanhoudt.

OpmerkingCode in dit hoofdstuk — wat moet je kunnen typen?

Code die je zelf moet kunnen typen staat open: options(contrasts = c("contr.sum", "contr.poly")) en aan het eind weer terug naar "contr.treatment", lm(), de parallellisme-check lm(Y ~ factor * covariaat), Anova(., type = 3), etaSquared(., type = 3), summary() voor $b_w$, aggregate() en table(), emmeans_test() en get_emmeans(). Code die alleen een plaatje of een tabel mooi maakt — kleuren en labels, Venn-diagrammen, de tabel-hulpfunctie van dit oefenboek — zit achter een knopje “Toon code”. Klap hem open als je nieuwsgierig bent.

Base R of dplyr. Waar je gegevens samenvat (aantallen, gemiddelden en standaardafwijkingen per groep), staan twee tabbladen naast elkaar: base R en dplyr. Ze geven dezelfde getallen, al staan de rijen soms in een andere volgorde en drukt dplyr minder cijfers af (een tibble laat er standaard drie zien; het getal erachter is hetzelfde); kies de schrijfwijze die je al gewend bent. De modelaanroepen zelf (lm(), Anova()) hebben maar één vorm.

3.0 Project- en datavoorbereiding

De bever houdt zijn aantekeningen netjes

Download de bron van dit oefenboek (de zip op de voorpagina), pak hem uit en dubbelklik op mvda_werkboek.Rproj — één project voor alle thema’s. Open daarna 03_covariantieanalyse/03_covariantieanalyse.qmd en draai de code vanuit dat bestand: RStudio voert een chunk uit in de map van het hoofdstuk, dus de paden kloppen. De data van dit thema ligt in 03_covariantieanalyse/data/, de helper-functies in common/R/. Niets te installeren als de vijf pakketten (car, lsr, emmeans, rstatix, tidyverse) op je computer staan; anders één keer (tidyverse brengt dplyr mee, voor de tabbladen):

install.packages(c("car", "lsr", "emmeans", "rstatix", "tidyverse"))

1. Open de hoofd-dataset.

load("data/bever_tevredenheid.RData")
str(bever_tevredenheid)

'data.frame':   90 obs. of  3 variables:
 $ damtevredenheid  : int  67 64 57 64 50 75 62 50 52 60 ...
 $ meedoende_familie: int  5 6 5 8 5 7 8 0 3 0 ...
 $ bouwwoede        : Factor w/ 3 levels "perfectionist",..: 1 1 1 1 1 1 1 1 1 1 ...

Per bever drie waarden: bouwwoede (perfectionist/pragmaticus/impulsief), meedoende_familie (aantal familieleden, $0$–$15$) en damtevredenheid (score $0$–$100$). bouwwoede is een factor; meedoende_familie en damtevredenheid zijn getallen (int in de uitvoer: hele getallen).

3.1 Het ANCOVA-model fitten en interpreteren

De bever vergelijkt drie werkstijlen, gegeven de hoeveelheid familie-hulp

“Als ik alleen de drie werkstijlen vergelijk,” zei de bever, “dan vergeet ik dat de pragmatici toevallig vaker hulp van familie kregen. Niet eerlijk. Ik moet het gecorrigeerd voor meedoende familie zien.”

OpmerkingVoor je gaat rekenen — vijf vragen aan jezelf

Het stappenplan, opnieuw, met de techniek-keuze als laatste stap.

1. Wie of wat wordt er gemeten? $90$ verschillende bevers, elk met één eigen dam.

2. Wat wordt er gemeten? Drie variabelen per bever: bouwwoede, meedoende_familie, damtevredenheid.

3. Onafhankelijk of afhankelijk? Onderzoeksvraag: voorspelt bouwwoede de damtevredenheid, gecorrigeerd voor meedoende_familie? Dan is damtevredenheid afhankelijk ($Y$), bouwwoede de hoofd-onafhankelijke (factor), en meedoende_familie de covariaat ($C$) waarvoor we corrigeren.

4. Meetniveau van elke variabele? damtevredenheid interval, meedoende_familie interval, bouwwoede nominaal (3 niveaus).

5. 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>thema 2</i>]
    C -->|gemengd| G[ANCOVA<br/><i>dit thema</i>]
    style G fill:#faf3e2,stroke:#c9a05a,stroke-width:2px

OpmerkingPijltjes-overzicht — alle MVDA-technieken in één tabel

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 (damtevredenheid), $X$-en zijn gemengd (covariaat meedoende_familie is interval, factor bouwwoede is nominaal). Antwoord: ANCOVA.

T1 — Techniek-keuze: ANCOVA, ANOVA, MRA of factorial ANOVA?

OpmerkingVraag T1 — Pen-en-papier

Voor elk vignet: bepaal afhankelijke variabele, meetniveau van elke voorspeller, en kies dan de techniek.

a) Een vleermuis vergelijkt drie jachtgebieden (open / dicht / boven water) op vangstsucces (aantal insecten/uur). Bij elke vlucht meet hij óók de luchttemperatuur (in graden Celsius). Hij wil weten of het jachtgebied uitmaakt na correctie voor temperatuur.

b) Een mus vergelijkt drie nest-types (boom / dakgoot / klimop) op kuikenoverleving (aantal jongen dat uitvliegt). Geen continue covariaten in zijn dataset.

c) Een arbeidspsycholoog meet bij $90$ medewerkers de werkprestatie (KPI-score, schaal $0$–$100$). Voorspellers: aantal jaren ervaring, aantal trainingen gevolgd en zelf-gerapporteerde motivatie — allemaal continu.

d) Een onderwijskundige vergelijkt drie didactische methoden (frontaal / probleemgestuurd / blended) en twee niveaus (HBO / WO) op tentamen-cijfer. Bij $200$ studenten, geen continue covariaat.

OpgeletAntwoord T1 — open na je eigen poging

a) $Y$ = vangstsucces (INT), één factor (jachtgebied, NOM, 3 niveaus) + één continue covariaat (temperatuur, INT). $\Rightarrow$ ANCOVA.

b) $Y$ = kuikenoverleving (INT), één factor (nest-type, NOM, 3 niveaus). $\Rightarrow$ eenweg-ANOVA (thema 2).

c) $Y$ = werkprestatie (INT), drie continue voorspellers (ervaring, trainingen, motivatie). Geen factor. $\Rightarrow$ multiple regressie-analyse (thema 1).

d) $Y$ = tentamen-cijfer (INT), twee factoren (didactische methode 3 niveaus, niveau 2 niveaus), geen continue covariaat. $\Rightarrow$ factorial / tweeweg-ANOVA (thema 2).

Vuistregel. Eén $Y$ + één factor + één continue covariaat = ANCOVA. Verwissel niet met factorial ANOVA: daar zijn beide voorspellers nominaal.

3.1.a Verkenning — scatter Y vs covariaat per groep

Voordat de bever toetst, kijkt hij eerst.

Algemene vorm.

table(mijn_data$factor)
aggregate(Y ~ factor, data = mijn_data, FUN = mean)
aggregate(C ~ factor, data = mijn_data, FUN = mean)

# Scatter Y vs C, gekleurd per niveau van de factor.
ggplot(mijn_data, aes(x = C, y = Y, colour = factor)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

Voor onze dieren.

Base R

# Aantallen per bouwwoede.
table(bever_tevredenheid$bouwwoede)

perfectionist   pragmaticus     impulsief 
           30            32            28 

# Raw groepsgemiddelden van de afhankelijke en de covariaat.
aggregate(damtevredenheid   ~ bouwwoede, data = bever_tevredenheid, FUN = mean)

      bouwwoede damtevredenheid
1 perfectionist        59.66667
2   pragmaticus        80.84375
3     impulsief        75.71429

aggregate(meedoende_familie ~ bouwwoede, data = bever_tevredenheid, FUN = mean)

      bouwwoede meedoende_familie
1 perfectionist             4.900
2   pragmaticus             8.375
3     impulsief             7.000

dplyr

# Alle drie in één tabel: per werkstijl het aantal, het gemiddelde van Y en
# het gemiddelde van de covariaat.
bever_tevredenheid |>
  group_by(bouwwoede) |>
  summarise(n                 = n(),
            damtevredenheid   = mean(damtevredenheid),
            meedoende_familie = mean(meedoende_familie))

# A tibble: 3 × 4
  bouwwoede         n damtevredenheid meedoende_familie
  <fct>         <int>           <dbl>             <dbl>
1 perfectionist    30            59.7              4.9 
2 pragmaticus      32            80.8              8.38
3 impulsief        28            75.7              7   

# Tijdelijke factor met hoofdletter-niveaus voor leesbare legend.
.stijl_lbl <- factor(tools::toTitleCase(as.character(bever_tevredenheid$bouwwoede)),
                     levels = tools::toTitleCase(levels(bever_tevredenheid$bouwwoede)))

ggplot(bever_tevredenheid,
       aes(x = meedoende_familie, y = damtevredenheid, colour = .stijl_lbl)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_colour_manual(values = c("#0077BB", "#EE7733", "#009988")) +
  labs(x = "Meedoende familieleden (aantal)",
       y = "Damtevredenheid (score 0-100)",
       colour = "Bouwwoede") +
  theme_minimal(base_size = 11)

OpmerkingVragen 3.1

a) Hoeveel dammen per werkstijl? Gebalanceerd?

b) Wat zijn de raw groepsgemiddelden voor damtevredenheid? In welke richting wijst het ruwe verschil?

c) Wat zijn de groeps-covariaat-means? Verschillen ze opvallend?

d) Lopen de regressielijnen in de scatter ongeveer parallel? Wat zegt dat?

OpgeletAntwoord 3.1 — open na je eigen poging

In gewone woorden.

1. Per werkstijl $30$, $32$ en $28$ dammen — lichtelijk ongebalanceerd ($n_{\max}/n_{\min} = 32/28 \approx 1.14$, ruim binnen de $1.5$-robuustheidsgrens), in de praktijk fijn.

2. Raw means: perfectionist $M = 59.67$, pragmaticus $M = 80.84$, impulsief $M = 75.71$ — pragmatici zijn het tevredenst, perfectionisten het minst.

3. Covariaat-means: $\bar{C}_{\text{perfectionist}} = 4.90$, $\bar{C}_{\text{pragmaticus}} = 8.38$, $\bar{C}_{\text{impulsief}} = 7.00$ — perfectionisten kregen systematisch minder familie-hulp dan de andere twee (familie haakt af, of perfectionist wil het zelf doen). Een deel van het ruwe tevredenheids-verschil zou dus ook door familie-hulp kunnen komen, niet door werkstijl. Daar gaat ANCOVA over.

4. De drie regressielijnen lopen ongeveer parallel — dat is een visuele eerste check op de parallellisme-aanname. We toetsen hem in 3.2.

APA-stijl.

Per werkstijl werden ongeveer $30$ dammen gemeten ($n_{\text{perfectionist}} = 30$, $n_{\text{pragmaticus}} = 32$, $n_{\text{impulsief}} = 28$). De raw tevredenheid was het hoogst bij pragmatici ($M = 80.84$), gevolgd door impulsieven ($M = 75.71$) en perfectionisten ($M = 59.67$). De groepen verschilden in gemiddelde familie-hulp (perfectionist $M = 4.90$; pragmaticus $M = 8.38$; impulsief $M = 7.00$), wat illustreert waarom ANCOVA-correctie hier zinvol is.

3.1.a-bis Systematische bias op de covariaat — toetsen vóór ANCOVA

Visueel zagen we al dat de perfectionisten gemiddeld minder familie-hulp kregen dan de andere werkstijlen. Voor ANCOVA is dat een observatie die de hele analyse draagt: alleen als er systematisch verschil zit in de covariaat over de factor-niveaus, dóét de ANCOVA-correctie iets noemenswaardigs. Toets het expliciet — niet via oogwerk alleen.

OpmerkingWaarom toetsen vóór ANCOVA?

Drie scenario’s:

• Geen bias ($\bar{C}$ gelijk over groepen): ANCOVA-resultaat ≈ ANOVA-resultaat. De covariaat doet weinig groeps-relevant werk; je gebruikt ’m dan vooral voor foutreductie.
• Wel bias (groepen verschillen in $\bar{C}$): ANCOVA en ANOVA kunnen tot andere conclusies komen. Dat is de situatie waarvoor ANCOVA ontworpen is.
• Te grote bias (groepen overlappen nauwelijks in $C$): adjusted means liggen ver buiten het waarneembare bereik per groep — extrapolatie. Wantrouw de uitkomst.

Helder maken in welk scenario je zit, voordat je de ANCOVA draait — dan weet je wat je verwacht en kun je achteraf herkennen of de uitkomst klopt.

Algemene vorm.

# Eenweg-ANOVA van de covariaat op de factor.
lm_cov <- lm(C ~ factor, data = mijn_data)
Anova(lm_cov, type = 3)

# Pairwise vergelijkingen tussen factor-niveaus op de covariaat.
emmeans(lm_cov, pairwise ~ factor, adjust = "tukey")

Voor onze dieren.

# Stap 1: eenweg-ANOVA van meedoende_familie op bouwwoede.
lm_cov_bias <- lm(meedoende_familie ~ bouwwoede, data = bever_tevredenheid)
Anova(lm_cov_bias, type = 3)

Anova Table (Type III tests)

Response: meedoende_familie
            Sum Sq Df F value      Pr(>F)    
(Intercept) 4098.6  1 614.572   < 2.2e-16 ***
bouwwoede    188.8  2  14.154 0.000004765 ***
Residuals    580.2 87                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Stap 2: pairwise via emmeans om te zien tussen welke werkstijlen het verschil zit.
emmeans(lm_cov_bias, pairwise ~ bouwwoede, adjust = "tukey")$contrasts

 contrast                    estimate    SE df t.ratio p.value
 perfectionist - pragmaticus    -3.48 0.656 87  -5.295 <0.0001
 perfectionist - impulsief      -2.10 0.679 87  -3.095  0.0074
 pragmaticus - impulsief         1.38 0.668 87   2.058  0.1048

P value adjustment: tukey method for comparing a family of 3 estimates 

OpmerkingVragen 3.1.a-bis

a) Is de eenweg-ANOVA op meedoende_familie ~ bouwwoede significant? Welk $F$, welk $p$?

b) Welke werkstijlen verschillen pairwise in de hoeveelheid familie-hulp die ze krijgen?

c) Welk scenario uit het kader hierboven past hier (geen bias / wel bias / te grote bias)?

d) Wat verwacht je dat ANCOVA met de tevredenheids-rangorde gaat doen — versterken, afzwakken, of omkeren ten opzichte van de raw means? Beargumenteer.

OpgeletAntwoord 3.1.a-bis — open na je eigen poging

a) De ANOVA op meedoende_familie ~ bouwwoede is significant, $F(2, 87) = 14.15$, $p < .001$ — er is wel systematische bias op de covariaat.

b) Pragmatici krijgen significant meer familie-hulp dan perfectionisten; impulsieven zitten ertussenin (significant verschil met perfectionist, niet met pragmaticus). Het verschil zit dus vooral bij perfectionisten — die werken meer alleen.

c) Wel bias — een serieus mean-verschil in $C$ over de groepen, maar de groepen overlappen nog steeds in waarneembaar bereik. Dit is de situatie waarvoor ANCOVA bedoeld is.

d) Verwachting: ANCOVA gaat het ruwe tevredenheids-verschil tussen perfectionisten en de andere twee afzwakken. Een deel van het ruwe verschil komt namelijk niet door werkstijl maar door minder familie-hulp — en familie-hulp helpt de tevredenheid. De adjusted means corrigeren daarvoor: ze laten zien wat de tevredenheid per werkstijl zou zijn als alle drie de groepen evenveel familie-hulp hadden gehad.

Inzicht. Deze toets-vooraf is wat onderscheid maakt tussen “we draaien ANCOVA omdat het hoort” en “we draaien ANCOVA omdat we weten waarom”. Bij geen bias is ANCOVA-correctie cosmetisch. Bij wel bias verandert het verhaal. Bij te grote bias kun je de uitkomst niet meer vertrouwen.

3.1.b Aannamechecks — kort

ANCOVA leunt op de drie aannames van ANOVA en drie eigen:

• Normaal verdeelde residuen binnen elke groep (bij grote $n$ minder kritisch);
• Gelijke varianties homogeneity of variance — met de vuistregel: zet de standaardafwijkingen per groep naast elkaar ($s_{\max}/s_{\min} < 2$). Verschillen ze toch, dan beslissen de groepsgroottes of dat erg is (zie Robuustheid van $F$ hieronder);
• Onafhankelijke observaties — een kwestie van het design, geen toets;
• Een covariaat zonder meetfout — ook geen toets, wel een vraag aan je meetinstrument (zie T7);
• Lineariteit van de covariaat-Y-relatie binnen elke groep — visueel via de scatter;
• Parallellisme — de slopes binnen elke groep zijn gelijk. Apart getoetst in 3.2.

Algemene vorm.

aggregate(Y ~ factor, data = mijn_data, FUN = sd)   # SD per groep
plot(lm(Y ~ C + factor, data = mijn_data), which = 1)

Voor onze dieren. De aantallen per werkstijl heb je uit 3.1.a; nu de standaardafwijkingen.

Base R

# Standaardafwijking van damtevredenheid per werkstijl.
aggregate(damtevredenheid ~ bouwwoede, data = bever_tevredenheid, FUN = sd)

      bouwwoede damtevredenheid
1 perfectionist        8.441823
2   pragmaticus        8.861331
3     impulsief       17.602309

dplyr

bever_tevredenheid |>
  group_by(bouwwoede) |>
  summarise(damtevredenheid = sd(damtevredenheid))

# A tibble: 3 × 2
  bouwwoede     damtevredenheid
  <fct>                   <dbl>
1 perfectionist            8.44
2 pragmaticus              8.86
3 impulsief               17.6 

m_ancova <- lm(damtevredenheid ~ bouwwoede + meedoende_familie, data = bever_tevredenheid)
plot(m_ancova, which = 1, sub.caption = "",
     caption = "Residuen tegen voorspelde waarden")

plot(m_ancova, which = 4, sub.caption = "",
     caption = "Cook's distance per observatie")

OpmerkingRobuustheid van $F$ — wanneer mag je doorgaan

Voor ANCOVA gelden dezelfde robuustheidsregels als voor ANOVA — plus de extra parallellisme-aanname. De korte versie: zijn je groepen groot genoeg en ongeveer even groot, maak je dan niet druk. De lange versie:

• Normaliteit — $F$ is robuust zolang $n \geq 15$ in elke groep. Onder die grens helpt de centrale limietstelling onvoldoende.
• Homoscedasticiteit — een echt probleem is het pas als twee dingen tegelijk gelden: de varianties verschillen ($s_{\max}/s_{\min} \geq 2$) én de groepen zijn ongelijk groot ($n_{\max}/n_{\min} \geq 1.5$). Zijn de groepen ongeveer even groot ($n_{\max}/n_{\min} \leq 1.5$), dan is $F$ robuust, ook bij ongelijke varianties. Waarom: dan telt elke groep even zwaar mee in de geschatte foutvariantie, en niets trekt de toets één kant op. Het gaat pas mis als de grootste variantie in de kleinste groep zit ($F$ te liberaal: te snel significant, vaker een Type-I-fout) of in de grootste ($F$ te conservatief) — de richtingsregel uit thema 2. Een gecorrigeerde toets zoals Brown-Forsythe uit thema 2 is gemaakt voor een eenweg-ANOVA zonder covariaat; bij ANCOVA houd je het, net als bij een factorial design, bij de richtingsregel.
• Onafhankelijkheid — niet te redden met een correctie; designkeuze + repeated-measures (thema 6).
• Parallellisme — extra aan ANCOVA. Apart getoetst in 3.2 via de interactie covariaat $\times$ factor. Is die geschonden, dan helpt geen correctie meer, want dan is ANCOVA niet het juiste model (zie 3.2.c).

OpmerkingVragen 3.1

e) Mag je de $F$-toets vertrouwen? Loop de vuistregels langs: $n$ per groep, $n_{\max}/n_{\min}$ en $s_{\max}/s_{\min}$.

f) Wat zie je in de residual plot — een rustige horizontale band, of iets verdachts?

g) Eén observatie heeft een opvallend hoge Cook’s distance. Hoe ga je daarmee om in een rapportage? (Geen R-actie nodig: bedenk de keuze.)

OpgeletAntwoord 3.1 — open na je eigen poging

In gewone woorden.

5. Hier zit een haakje. De standaardafwijkingen lopen flink uiteen: impulsief $17.60$, perfectionist $8.44$, pragmaticus $8.86$. Dan is $s_{\max}/s_{\min} = 17.60/8.44 = 2.09$ — boven de $2$, dus de varianties zijn níét ongeveer gelijk. Dat zag je al in de scatter van 3.1.a: de impulsieven liggen veel verder uiteen. Een paar impulsieven zijn extreem tevreden, een paar balen. Maar de groepen zijn ongeveer even groot ($n_{\max}/n_{\min} = 32/28 \approx 1.14$), en elke groep telt minstens $15$ bevers. Dan blijft de $F$-toets volgens de robuustheidsvuistregel verdedigbaar — maar in een rapportage meld je de ongelijke spreiding eerlijk en bespreek je de implicatie.

6. De residual plot toont een lichte trompet (grotere spreiding rechts) door de impulsieven; geen kromming, wel ongelijke variantie.

7. Eén observatie steekt boven de andere uit met een Cook’s distance van $0.31$ — niet alarmerend (vuistregel: $> 0.5$ is reden tot inspectie, $> 1$ tot zorg). Toch netjes: rapporteer dat je hem hebt opgemerkt, draai het model met-en-zonder die observatie, en vermeld of de conclusies veranderen. Niet doen: stilletjes weglaten zonder dat te melden.

APA-stijl.

De varianties verschilden tussen de werkstijlen ($s_{\max}/s_{\min} = 2.09$), voornamelijk doordat de impulsief-groep een veel grotere spreiding had ($SD = 17.60$) dan perfectionisten ($SD = 8.44$) en pragmatici ($SD = 8.86$). Omdat de groepen ongeveer even groot waren ($n_{\max}/n_{\min} = 1.14$), is $F$ in deze configuratie nog robuust; we rapporteren de heteroscedasticiteit expliciet bij de bevindingen. De residual plot toonde een lichte trompetvorm in lijn hiermee. De maximale Cook’s distance was $D = 0.31$.

OpmerkingVoor wie verder wil — Levene’s toets voor gelijke varianties

Er bestaat ook een toets voor gelijke varianties: Levene’s toets Levene’s test. Dit hoofdstuk werkt met de vuistregels en niet met die toets. Levene zegt óf de varianties verschillen; de vuistregels zetten $s_{\max}/s_{\min}$ naast $n_{\max}/n_{\min}$ en zeggen of dat verschil $F$ in gevaar brengt. In SPSS-uitvoer staat Levene er vaak wel bij, dus je komt hem tegen.

# Levene's toets zit in het pakket car (dat laadde je al voor Anova()).
leveneTest(damtevredenheid ~ bouwwoede, data = bever_tevredenheid)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value   Pr(>F)   
group  2  5.9785 0.003691 **
      87                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De nulhypothese is dat de varianties in de populatie gelijk zijn. Bij de bevers is Levene significant, $F(2, 87) = 5.98$, $p = .004$: hetzelfde oordeel als de vuistregel $s_{\max}/s_{\min} \geq 2$. En ook hier geldt: zijn de groepen ongeveer even groot, dan blijft de $F$-toets robuust, ook bij een significante Levene.

3.1.c Het ANCOVA-model en de Type-III ANOVA-tabel

OpmerkingVolgorde in de formule — factor eerst, of covariaat eerst?

Bij Type-III SS via Anova(., type = 3) maakt de volgorde voor de uiteindelijke $F$- en $p$-waarden geen verschil — Type III is invariant onder volgorde. Of je nu lm(Y ~ factor + covariaat) of lm(Y ~ covariaat + factor) schrijft, je krijgt dezelfde toets-statistieken. Wat afwijkt is alleen de Type-I (sequentiële) SS-uitsplitsing, en daar werken we hier sowieso niet mee.

Dit oefenboek schrijft consequent lm(Y ~ factor + covariaat): factor eerst, covariaat daarna. Eén vaste volgorde, dan lees je elke tabel in dit thema op dezelfde manier. Bij Type-I SS zou de volgorde wél uitmaken; daarom blijft Anova(., type = 3) de juiste keuze.

Algemene vorm.

# Eerst de codering, dan het model: R legt hem vast bij lm().
options(contrasts = c("contr.sum", "contr.poly"))
mijn_model <- lm(Y ~ factor + covariaat, data = mijn_data)
Anova(mijn_model, type = 3)            # ANCOVA-tabel, type III
etaSquared(mijn_model, type = 3)       # eta-squared en partial eta-squared
summary(mijn_model)                    # de rij 'covariaat' geeft b_w

Voor onze dieren.

# Het ANCOVA-model: factor + covariaat (zonder interactie), factor eerst,
# covariaat daarna. De codering staat al bovenaan op contr.sum.
m_ancova <- lm(damtevredenheid ~ bouwwoede + meedoende_familie, data = bever_tevredenheid)
Anova(m_ancova, type = 3)

Anova Table (Type III tests)

Response: damtevredenheid
                  Sum Sq Df F value     Pr(>F)    
(Intercept)        38882  1 311.577  < 2.2e-16 ***
bouwwoede           2808  2  11.252 0.00004562 ***
meedoende_familie   2135  1  17.105 0.00008200 ***
Residuals          10732 86                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Effectgroottes.
etaSquared(m_ancova, type = 3)

                     eta.sq eta.sq.part
bouwwoede         0.1382361   0.2074074
meedoende_familie 0.1050706   0.1659019

# De helling binnen de groepen, b_w: de rij van de covariaat in summary().
summary(m_ancova)$coefficients["meedoende_familie", ]

     Estimate    Std. Error       t value      Pr(>|t|) 
1.91808859014 0.46376971285 4.13586428133 0.00008199635 

Die laatste regel geeft je $b_w$, de helling binnen de groepen: $b_w = 1.92$. Lees hem zo: binnen dezelfde werkstijl is een bever met één meedoend familielid extra gemiddeld 1.92 punten tevredener. Met dat getal rekent ANCOVA de gemiddelden om — dat doe je zelf in T2.

TipVrijheidsgraden bij ANCOVA — hoofdrekenbaar

Bij één covariaat en één factor met $I$ niveaus, $N$ observaties:

Bron	df
Covariaat	$1$
Factor	$I - 1$
Residu	$N - I - 1$
Totaal	$N - 1$

Som-check: $1 + (I - 1) + (N - I - 1) = N - 1$.

Hoofdrekenbaar voorbeeld. $I = 3$ groepen, $N = 60$: $df_C = 1$, $df_F = 2$, $df_w = 60 - 3 - 1 = 56$, $df_{\text{tot}} = 59$. Som: $1 + 2 + 56 = 59$.

Voor onze bevers. $I = 3$, $N = 90$: $df_C = 1$, $df_F = 2$, $df_w = 90 - 3 - 1 = 86$, $df_{\text{tot}} = 89$. Check in de Anova-output hierboven of de tweede df van elke $F$-toets inderdaad $86$ is.

OpmerkingVragen 3.1

h) Is de covariaat meedoende_familie significant? Rapporteer $F$, $df$, $p$ en $\eta^2_p$.

i) Is bouwwoede significant na correctie voor meedoende_familie? Rapporteer $F$, $df$, $p$ en $\eta^2_p$.

j) Vergelijk de effectgroottes: wie verklaart meer, de covariaat of de factor?

k) Klopt $df_w$ met de hoofdreken-formule?

OpgeletAntwoord 3.1 — open na je eigen poging

In gewone woorden.

8. Meedoende_familie heeft een sterk significant effect: $F(1, 86) = 17.11$, $p < .001$, $\eta^2_p = .17$.

9. Bouwwoede is ook significant na correctie voor familie-hulp: $F(2, 86) = 11.25$, $p < .001$, $\eta^2_p = .21$.

10. De factor verklaart $\sim 21\%$ van de variantie die na de covariaat overblijft, de covariaat $\sim 17\%$ van wat na de factor overblijft — werkstijl weegt iets zwaarder, en familie-hulp doet er duidelijk óók na correctie toe.

11. $df_w = 86 = N - I - 1 = 90 - 3 - 1$, klopt.

APA-stijl.

Een ANCOVA met bouwwoede als factor en meedoende_familie als covariaat toonde een significant effect van bouwwoede, $F(2, 86) = 11.25$, $p < .001$, $\eta^2_p = .21$. Het effect van familie-hulp was eveneens significant, $F(1, 86) = 17.11$, $p < .001$, $\eta^2_p = .17$.

OpmerkingVoor wie verder wil — wanneer heet een effect klein, middelgroot of groot?

Je leest vaak de woorden klein, middelgroot en groot zonder dat er grenzen bij staan. Wie grenzen wil, vindt ze bij Cohen (1988): $\eta^2_p \approx .01$ klein, $\approx .06$ middelgroot, $\approx .14$ groot. Houd rekening met je vakgebied: in experimentele psychologie zijn dat realistische maten, in toegepaste contexten ligt de lat soms lager. En lees een $\eta^2_p$ altijd als wat hij is: het deel van de variantie dat na de andere termen overblijft.

T2 — Adjusted means hand-bereken

OpmerkingVraag T2 — Pen-en-papier

Adjusted means worden geschat met de formule: $$\bar{Y}_j^* = \bar{Y}_j - b_w (\bar{C}_j - \bar{C})$$

Voor de bever-data ben je gegeven:

• $\bar{Y}_{\text{perfectionist}} = 59.67$, $\bar{Y}_{\text{pragmaticus}} = 80.84$, $\bar{Y}_{\text{impulsief}} = 75.71$
• $\bar{C}_{\text{perfectionist}} = 4.90$, $\bar{C}_{\text{pragmaticus}} = 8.38$, $\bar{C}_{\text{impulsief}} = 7.00$
• Grand mean covariaat $\bar{C} = 6.79$ — het gemiddelde over alle 90 bevers (mean(bever_tevredenheid$meedoende_familie)), niet het gemiddelde van de drie groepsgemiddelden; bij ongelijke groepen scheelt dat iets
• Pooled-within slope $b_w = 1.92$ — de rij meedoende_familie in summary(m_ancova), zie 3.1.c

a) Bereken $\bar{Y}_{\text{perfectionist}}^*$, $\bar{Y}_{\text{pragmaticus}}^*$ en $\bar{Y}_{\text{impulsief}}^*$.

b) Welke groep krijgt de grootste opwaartse correctie? Welke de grootste neerwaartse? Klopt dat met wat je zou verwachten gegeven de covariaat-means?

c) Vergelijk de drie raw means met de drie adjusted means: welke twee groepen liggen na correctie dichter bij elkaar dan ervoor?

OpgeletAntwoord T2 — open na je eigen poging

a)

$\bar{Y}_{\text{perfectionist}}^* = 59.67 - 1.92 \cdot (4.90 - 6.79) = 59.67 - 1.92 \cdot (-1.89) = 59.67 + 3.63 = 63.30$.

$\bar{Y}_{\text{pragmaticus}}^* = 80.84 - 1.92 \cdot (8.38 - 6.79) = 80.84 - 1.92 \cdot (1.59) = 80.84 - 3.05 = 77.79$.

$\bar{Y}_{\text{impulsief}}^* = 75.71 - 1.92 \cdot (7.00 - 6.79) = 75.71 - 1.92 \cdot (0.21) = 75.71 - 0.40 = 75.31$.

De tabel met gecorrigeerde gemiddelden in 3.1.d geeft perfectionist 63.29, pragmaticus 77.80, impulsief 75.31 — in de tweede decimaal soms net anders. Dat is afronding: R rekent met de onafgeronde $b_w = 1.9181$ en onafgeronde gemiddelden. Kom jij met de hand op $63.30$ uit, dan heb je het goed.

b) De perfectionist-groep krijgt de grootste opwaartse correctie ($+3.63$), omdat hun covariaat-mean het laagst was (minder familie-hulp — een handicap die je wegcorrigeert). Pragmaticus krijgt de grootste neerwaartse correctie ($-3.05$): hun covariaat-mean lag duidelijk boven gemiddeld (meer familie-hulp), dus je trekt af. Impulsief zat dicht bij het overall covariaat-gemiddelde; vrijwel geen correctie. Dat past bij wat je van ANCOVA verwacht: groepen met betere covariaat-startpositie krijgen een correctie omlaag.

c) De raw means lopen $59.67 \to 80.84$ (range $\approx 21.2$). De adjusted means $63.30 \to 77.79$ (range $\approx 14.5$). De groepen liggen na correctie dichter bij elkaar — een deel van het ruwe verschil zat in de covariaat (familie-hulp), niet in werkstijl.

Inzicht. ANCOVA “egaliseert” niet altijd; soms maakt hij verschillen groter (als groepen in tegengestelde richting verschillen op de covariaat). Hier kleiner, maar nog steeds substantieel — perfectionisten zijn ook na correctie minder tevreden dan pragmatici.

T3 — $b_w$ pooled-within: wat is gepooled?

OpmerkingVraag T3 — Pen-en-papier

In de ANCOVA-formule staat $b_w$, de pooled-within slope.

a) Wat betekent “pooled-within” letterlijk? Welke slopes worden gepooled?

b) Hoe verschilt $b_w$ van de gewone OLS-slope $b_{\text{totaal}}$ die je krijgt als je lm(Y ~ C) doet zonder factor? (Hint: bij verschillende groeps-Y-means kunnen ze flink verschillen.)

c) Onder welke aanname is $b_w$ een zinvolle parameter? Wat als die aanname niet houdt?

OpgeletAntwoord T3 — open na je eigen poging

a) “Pooled-within” = gewogen gemiddelde van de regressie-slopes binnen elke groep. Per groep $j$ heb je een slope $b_j$ (de slope van Y op C als je alleen die groep bekijkt). $b_w$ is daar het pooled (gewogen) gemiddelde van.

b) $b_{\text{totaal}}$ uit lm(Y ~ C) zonder factor mengt twee bronnen: (1) de echte binnen-groep-relatie van C op Y, en (2) het feit dat groepen zowel op C als op Y verschillen (cross-group-variatie). Als groepen met hogere $\bar{C}$ ook hogere $\bar{Y}$ hebben omdat de groep beter is (niet alleen omdat C ze omhoog trekt), wordt $b_{\text{totaal}}$ opgepompt. $b_w$ kijkt door de groepen heen en negeert die between-group-trend — vandaar dat het de juiste parameter is voor adjusted means.

c) Onder de parallellisme-aanname: alle binnen-groep-slopes zijn (in de populatie) gelijk. Dan bestaat $b_w$ als één getal en heeft hij betekenis. Als de slopes echt verschillen tussen groepen (significante interactie covariaat $\times$ factor), bestaat er geen één pooled $b_w$ — het effect van de covariaat hangt dan af van de groep, en dan bestaan er ook geen gecorrigeerde gemiddelden meer die voor elke covariaat-waarde gelden. ANCOVA is dan niet de juiste analyse.

3.1.d Adjusted means en paarsgewijze vergelijking

Na een significante factor wil je weten welke groepen van elkaar verschillen, nu op de gecorrigeerde gemiddelden. Dat doe je hier met emmeans_test() uit het pakket rstatix, met Bonferroni-correctie: elke $p$-waarde maal het aantal paren (met $1$ als plafond). Dat is dezelfde Bonferroni als in thema 2, waar je $\alpha$ deelde door het aantal toetsen. Of je de $\alpha$ drie keer zo klein maakt en daarmee vergelijkt, of de $p$-waarde drie keer zo groot: hetzelfde. Bij p.adj blijf je dus gewoon op $.05$ denken. Na een ANOVA gebruikte je Tukey, en Tukey zou hier ook kunnen; dat het Bonferroni wordt, is een keuze, geen statistische wet. Dit oefenboek kiest Bonferroni omdat je hem met de hand kunt narekenen — $p$ maal het aantal paren, bij drie groepen dus drie — terwijl je voor Tukey een tabel of de computer nodig hebt. Hieronder zie je dat narekenen voor één rij voorgedaan.

Algemene vorm.

library(rstatix)
em <- emmeans_test(data = mijn_data, Y ~ factor, covariate = C,
                   p.adjust.method = "bonferroni")
em                # paarsgewijze vergelijking: lees de kolom p.adj
get_emmeans(em)   # de adjusted means per niveau

emmeans_test() schat zelf het ANCOVA-model, met de covariaat erin; je hoeft het model niet eerst met lm() te maken.

Voor onze dieren.

library(rstatix)

# Paarsgewijze vergelijking van de adjusted means, met Bonferroni-correctie.
em_bever <- emmeans_test(data = bever_tevredenheid, damtevredenheid ~ bouwwoede,
                         covariate = meedoende_familie, p.adjust.method = "bonferroni")
em_bever

# A tibble: 3 × 9
  term      .y.             group1        group2         df statistic         p     p.adj p.adj.signif
* <chr>     <chr>           <chr>         <chr>       <dbl>     <dbl>     <dbl>     <dbl> <chr>       
1 bouwwoede damtevredenheid perfectionist pragmaticus    86    -4.45  0.0000260 0.0000781 ****        
2 bouwwoede damtevredenheid perfectionist impulsief      86    -3.89  0.000200  0.000599  ***         
3 bouwwoede damtevredenheid pragmaticus   impulsief      86     0.842 0.402     1         ns          

# De adjusted means (EMM) per bouwwoede.
get_emmeans(em_bever)

# A tibble: 3 × 7
  bouwwoede     emmean    se    df conf.low conf.high method      
  <fct>          <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>       
1 perfectionist   63.3  2.22    86     58.9      67.7 Emmeans test
2 pragmaticus     77.8  2.11    86     73.6      82.0 Emmeans test
3 impulsief       75.3  2.11    86     71.1      79.5 Emmeans test

Zo lees je de uitvoer. Elke rij is een paar (group1 tegen group2). De kolom p is de gewone $p$-waarde, p.adj de $p$-waarde na Bonferroni-correctie: die rapporteer je. Reken het na in de derde rij: $.402 \cdot 3 = 1.21$, en boven de $1$ wordt dat $1$. Een p.adj van $1$ schrijf je als $p > .99$. De kolom statistic is de $t$-waarde; een minteken betekent dat group1 lager scoort dan group2. get_emmeans() geeft de adjusted means — dezelfde als emmeans(m_ancova, ~ bouwwoede) uit het model van 3.1.c, maar die laatste hoef je hier niet te typen. Zie je bij jezelf afgekapte namen (perfe…), maak dan je console breder of typ as.data.frame(em_bever).

Waarschuwingemmeans_test() knipt namen op het streepje

In dit hoofdstuk heeft geen niveau een streepje; dit gaat over je eigen data. Heeft een niveau een streepje in zijn naam (zoals hands-free), dan knipt emmeans_test() elk paar op dat streepje, en kloppen de kolommen group1 en group2 niet meer. Je ziet het aan een halve naam: dan staat er bijvoorbeeld free bij group2, of heten twee rijen hetzelfde paar. De $p$-waarden kloppen wel. Geef zo’n niveau vooraf een naam zonder streepje. Welk nummer het heeft, zie je met levels(mijn_data$factor); staat hands-free daar als tweede, dan wordt het levels(mijn_data$factor)[2] <- "handsfree".

Adjusted-means-tabel — raw versus gecorrigeerd in één blik

apa_ancova_table(
  m_ancova,
  factor          = "bouwwoede",
  covariate       = "meedoende_familie",
  data            = bever_tevredenheid,
  dependent_label = "damtevredenheid",
  factor_label    = "Bouwwoede",
  covariate_label = "meedoende_familie",
  caption         = "Bever — damtevredenheid per werkstijl, raw versus adjusted"
)

Bever — damtevredenheid per werkstijl, raw versus adjusted
Bouwwoede	n	Mmeedoende_familie	Mdamtevredenheid, ruw	Mdamtevredenheid, adjusted
Perfectionist	30	4.90	59.67	63.29
Pragmaticus	32	8.38	80.84	77.80
Impulsief	28	7.00	75.71	75.31
Note. Mdamtevredenheid, adjusted = estimated marginal mean voor damtevredenheid bij gemiddelde meedoende_familie (6.79). Pooled-within slope bw = 1.92; $\eta^2_p$ = partiële eta-kwadraat. meedoende_familie: F(1, 86) = 17.11, p < .001, $\eta^2_p$ = .17; Bouwwoede: F(2, 86) = 11.25, p < .001, $\eta^2_p$ = .21.

OpmerkingWat de tabel laat zien — drie kolommen, één verschil

De tabel zet de raw-mean van Y, de adjusted-mean van Y en de groeps-covariaat-mean naast elkaar. Het verschil tussen kolom $M_{Y, \text{ruw}}$ en $M_{Y, \text{adjusted}}$ is precies $b_w \cdot (\bar{C}_j - \bar{C})$ — de bias-correctie. Als de groep boven het overall covariaat-gemiddelde zat, gaat de adjusted-mean omlaag; eronder, omhoog.

OpmerkingVragen 3.1

l) Welke twee werkstijlen verschillen significant na correctie? Welke twee niet?

m) Schrijf één APA-zin met de adjusted means en de paarsgewijze conclusies.

OpgeletAntwoord 3.1 — open na je eigen poging

In gewone woorden. Adjusted means: perfectionist $\bar{Y}^* = 63.29$, pragmaticus $77.80$, impulsief $75.31$. Paarsgewijze vergelijkingen met Bonferroni-correctie: perfectionist versus pragmaticus significant ($p < .001$), perfectionist versus impulsief significant ($p < .001$), pragmaticus versus impulsief niet significant ($p > .99$). De perfectionist staat na correctie nog steeds duidelijk apart; pragmatici en impulsieven verschillen onderling niet betrouwbaar (l).

APA-stijl (m).

Na correctie voor meedoende familieleden waren de geschatte gemiddelde tevredenheids-scores $\bar{Y}^*_{\text{perfectionist}} = 63.29$, $\bar{Y}^*_{\text{pragmaticus}} = 77.80$ en $\bar{Y}^*_{\text{impulsief}} = 75.31$. Paarsgewijze vergelijkingen met Bonferroni-correctie lieten zien dat perfectionisten significant minder tevreden waren dan pragmatici ($p < .001$) en dan impulsieven ($p < .001$); het verschil tussen pragmatici en impulsieven was niet significant ($p > .99$).

T4 — ANOVA tegenover ANCOVA: wat doet de covariaat met $F$?

OpmerkingVraag T4 — Pen-en-papier (en daarna één R-chunk)

Hier zie je wat de covariaat doet met de twee helften van $F$: de noemer (de foutvariantie) en de teller (het effect van de factor). We vergelijken twee modellen op dezelfde data:

• ANOVA: $Y \sim \text{bouwwoede}$ (factor alleen).
• ANCOVA: $Y \sim \text{bouwwoede} + \text{meedoende\_familie}$ (factor eerst, covariaat daarna).

Stel je leest af:

• ANOVA-tabel: $SS_{\text{factor}} = 7449$, $SS_w = 12867$, $df_{\text{factor}} = 2$, $df_w = 87$.
• ANCOVA-tabel: $SS_{\text{factor}} = 2808$, $SS_{\text{covariaat}} = 2135$, $SS_w = 10732$, $df_{\text{factor}} = 2$, $df_{\text{covariaat}} = 1$, $df_w = 86$.

a) Bereken $F_{\text{factor}}$ in beide modellen.

b) Wat is er gebeurd met $SS_w$? Met $MS_w$?

c) Het $SS_{\text{factor}}$ daalde van $7449$ naar $2808$ — minder unieke variantie voor de factor. En toch is $F$ nog steeds duidelijk significant. Hoe kan dat?

d) Als meedoende_familie geen variantie in damtevredenheid zou hebben verklaard, wat zou er met $F_{\text{factor}}$ in ANCOVA gebeuren ten opzichte van ANOVA?

OpgeletAntwoord T4 — open na je eigen poging

a) ANOVA: $F = (7449/2) / (12867/87) = 3724.50 / 147.90 = 25.18$. ANCOVA: $F = (2808/2) / (10732/86) = 1404.00 / 124.79 = 11.25$.

b) $SS_w$ daalde van $12867$ naar $10732$ — een serieuze reductie. Meedoende_familie verklaarde een stuk van de Y-variantie dat in ANOVA “binnen-groep-ruis” was. $MS_w$ daalde van $147.90$ naar $124.79$ — de noemer van de $F$-toets is kleiner geworden.

c) Het $SS_{\text{factor}}$ kromp óók — en flink (van $7449$ naar $2808$). Dat is bias-correctie: een groot deel van wat in de ANOVA aan de factor leek toegekend, ging eigenlijk via familie-hulp (perfectionisten kregen minder hulp, dus hun raw-mean was extra laag). Na correctie blijft een kleiner maar zuiverder factor-effect over. De teller $MS_{\text{factor}}$ daalt dus ook, maar het effect blijft groot genoeg om sterk significant te zijn. In data waar de groepen niet op de covariaat verschillen, blijft de teller stabiel en daalt alleen de noemer — dan stijgt $F$ duidelijk. Dat is het zuivere foutreductie-effect.

d) Als meedoende_familie niets verklaarde, zou $SS_w^{\text{ancova}} \approx SS_w^{\text{anova}}$ — geen winst. Erger: je verliest één df ($df_w = 86$ in plaats van $87$), wat $F$ ietsje conservatiever maakt. Vandaar: een covariaat alleen toevoegen als je verwacht (of weet) dat hij Y-variantie verklaart.

Inzicht. Twee mechanismen lopen tegelijk in een ANCOVA: (1) foutreductie (de noemer daalt) en (2) bias-correctie (de teller wordt aangepast als groepen op C verschillen). Beide tegelijk maken het effect van de factor cleaner — niet automatisch sterker, maar wel eerlijker.

# Live: vergelijk ANOVA en ANCOVA op dezelfde data.
m_anova_only <- lm(damtevredenheid ~ bouwwoede, data = bever_tevredenheid)
Anova(m_anova_only, type = 3)

Anova Table (Type III tests)

Response: damtevredenheid
            Sum Sq Df  F value         Pr(>F)    
(Intercept) 466144  1 3151.922      < 2.2e-16 ***
bouwwoede     7449  2   25.185 0.000000002349 ***
Residuals    12867 87                            
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(m_ancova, type = 3)

Anova Table (Type III tests)

Response: damtevredenheid
                  Sum Sq Df F value     Pr(>F)    
(Intercept)        38882  1 311.577  < 2.2e-16 ***
bouwwoede           2808  2  11.252 0.00004562 ***
meedoende_familie   2135  1  17.105 0.00008200 ***
Residuals          10732 86                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Venn voor ANCOVA — drie cirkels, drie verhalen

Hoofdrekenbare illustratie van foutreductie

Stel een hoofdrekenbare ANCOVA met $SS_{\text{totaal}} = 100$. Daarvan:

• $SS_{\text{covariaat}} = 30$ — de covariaat trekt eerst $30\%$ variantie weg;
• $SS_{\text{factor}} = 20$ — de factor verklaart, na correctie, nog eens $20\%$;
• $SS_w = 50$ — residu, $50\%$.

In een gewone ANOVA zonder covariaat zou de factor het hele “covariaat + residu”-stuk in de noemer hebben: $SS_w^{\text{anova}} \approx 80$. Door de covariaat eruit te halen, zakt de noemer naar $50$ — dat is het foutreductie-mechanisme.

library(eulerr)

fit_ancova <- euler(c(
  "Covariaat" = 30,
  "Factor"    = 20,
  "Residu"    = 50
))
plot(fit_ancova,
     fills      = list(fill = c("#0077BB", "#CCBB44", "#BBBBBB"),
                       alpha = 0.55),
     labels     = list(col = "black", font = 2),
     quantities = list(type = "counts", cex = 0.9),
     edges      = list(col = "white", lwd = 2))

ANCOVA-Venn: $Y$-totaal verdeeld in covariaat ($30$), factor ($20$, na correctie), en residu ($50$).

TipWat het plaatje zegt

Drie disjuncte segmenten — analoog aan het $\eta^2$-plaatje uit thema 2. ANCOVA onderscheidt zich van een gewone tweeweg-ANOVA doordat de covariaat een continue variabele is, niet een tweede factor. De $F$-toets op de factor gebruikt $SS_{\text{factor}} / df_{\text{factor}}$ in de teller en $SS_w / df_w$ in de noemer — niet $SS_{\text{factor}} + SS_{\text{covariaat}}$ in de noemer. Het residu is “schoner” geworden.

Dit plaatje is het geval zonder bias: factor en covariaat raken elkaar niet, omdat de groepen niet op de covariaat verschillen. Bij de bevers verschillen ze wél (3.1.a-bis), en dan overlappen de cirkels. Tel maar na: de totale kwadratensom van damtevredenheid is $20316$, de drie stukken uit de ANCOVA-tabel samen zijn $15675$. Het verschil, $4641$, is variantie die factor en covariaat delen. Type III geeft dat gedeelde stuk aan geen van beide — en daarom zakte $SS_{\text{factor}}$ in T4 van $7449$ naar $2808$.

3.2 De parallellisme-aanname expliciet toetsen

“Lopen de slopes binnen elke groep wel ongeveer gelijk?”

“Voordat ik mijn ANCOVA-conclusies serieus neem,” zei de bever, “moet ik checken of het effect van familie-hulp binnen elke werkstijl ongeveer gelijk is. Anders mengt mijn $b_w$ ongelijke dingen door elkaar.”

De parallellisme-aanname zegt: binnen elke factor-groep is de regressie van Y op de covariaat dezelfde. Visueel: de drie regressielijnen (een per groep) lopen parallel. Statistisch: de interactie tussen factor en covariaat is niet significant.

3.2.a Het model met interactie

Algemene vorm.

options(contrasts = c("contr.sum", "contr.poly"))   # vóór lm()
m_par <- lm(Y ~ factor * covariaat, data = mijn_data)
Anova(m_par, type = 3)

De * in factor * covariaat betekent: factor-hoofdeffect + covariaat-hoofdeffect + interactie factor $\times$ covariaat. Het is de interactie waar het ons om gaat — en alléén die. De andere regels van deze tabel lees je niet: in een model mét interactie meet de toets op de factor het groepsverschil in $Y$ bij covariaat $= 0$, en dat is een andere vraag dan de ANCOVA stelt. Is de interactie niet significant, dan leg je dit model weg en ga je verder met het model zonder interactie.

Voor onze dieren.

# De codering staat al bovenaan op contr.sum.
m_par <- lm(damtevredenheid ~ bouwwoede * meedoende_familie, data = bever_tevredenheid)
Anova(m_par, type = 3)

Anova Table (Type III tests)

Response: damtevredenheid
                            Sum Sq Df  F value    Pr(>F)    
(Intercept)                  38284  1 299.8281 < 2.2e-16 ***
bouwwoede                      495  2   1.9395 0.1501661    
meedoende_familie             2113  1  16.5519 0.0001064 ***
bouwwoede:meedoende_familie      6  2   0.0243 0.9759587    
Residuals                    10726 84                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

OpmerkingVragen 3.2

a) Is de interactie bouwwoede $\times$ meedoende_familie significant? Rapporteer $F$, $df$ en $p$.

b) Wat betekent dat voor de parallellisme-aanname?

c) Mag je dan terug naar het additieve ANCOVA-model (bouwwoede + meedoende_familie)?

OpgeletAntwoord 3.2 — open na je eigen poging

In gewone woorden.

1. Interactie niet significant: $F(2, 84) = 0.02$, $p = .98$.

2. De parallellisme-aanname is dus houdbaar — er is geen overtuigend bewijs dat de slopes binnen de drie werkstijlen verschillen.

3. Conclusie: ANCOVA zoals in 3.1 mag.

APA-stijl.

De aanname van homogene regressielijnen (parallellisme) was houdbaar: de interactie tussen bouwwoede en meedoende_familie was niet significant, $F(2, 84) = 0.02$, $p = .98$. ANCOVA werd dus gerechtvaardigd toegepast.

WaarschuwingEen niet-significante interactie is geen bewijs van parallellisme

“Niet significant” betekent: te weinig bewijs om de nul (gelijke slopes) te verwerpen. Bij grote $N$ wordt de toets gevoeliger; bij kleine $N$ minder. Dus: kijk altijd ook naar de visualisatie en de effectgrootte. Als de slopes met het blote oog duidelijk verschillen ($\eta^2_p$ van de interactie groot) maar de toets niet significant is door kleine $N$ — wees voorzichtig met de ANCOVA-conclusie.

3.2.b Visualisatie — slopes per groep

# Per-groep regressielijnen om parallellism visueel te checken.
.stijl_lbl <- factor(tools::toTitleCase(as.character(bever_tevredenheid$bouwwoede)),
                     levels = tools::toTitleCase(levels(bever_tevredenheid$bouwwoede)))

ggplot(bever_tevredenheid,
       aes(x = meedoende_familie, y = damtevredenheid, colour = .stijl_lbl)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_colour_manual(values = c("#0077BB", "#EE7733", "#009988")) +
  labs(x = "Meedoende familieleden (aantal)",
       y = "Damtevredenheid (score 0-100)",
       colour = "Bouwwoede",
       title  = "Per-groep slopes — parallellisme visueel") +
  theme_minimal(base_size = 11)

OpmerkingDe drie lijnen lopen ongeveer parallel

In de plot zie je drie regressielijnen die ongeveer dezelfde helling hebben. Ze liggen niet op elkaar (dat zou betekenen: geen factor-effect), en ze lopen ongeveer evenwijdig. Dat laatste is wat parallellisme vraagt — niet of de lijnen elkaar kruisen. Twee lijnen met een duidelijk verschillende helling schenden de aanname ook als ze elkaar binnen de data nergens raken. Verticaal verschoven, horizontaal evenwijdig — dat is precies wat ANCOVA aanneemt. Het $b_w \approx 1.92$ van het additieve model is hun gewogen gemiddelde slope.

T5 — Bias-correctie: wat als groepen verschillen op de covariaat?

OpmerkingVraag T5 — Pen-en-papier

Stel twee groepen, A en B, op een uitkomst $Y$, met covariaat $C$.

• Groep A: $\bar{Y}_A = 50$, $\bar{C}_A = 4$.
• Groep B: $\bar{Y}_B = 60$, $\bar{C}_B = 8$.
• Pooled-within slope: $b_w = 2$.
• Grand mean covariaat: $\bar{C} = 6$.

a) Wat is het raw verschil tussen groepen A en B op Y?

b) Bereken de adjusted means $\bar{Y}_A^*$ en $\bar{Y}_B^*$.

c) Wat is het adjusted verschil? Hoe verhoudt het zich tot het raw verschil?

d) Wat zou het adjusted verschil zijn als $b_w = 0$ (covariaat verklaart niets)? En als $b_w = 2.5$?

e) Hoe noemen we deze tweede functie van ANCOVA, naast foutreductie?

OpgeletAntwoord T5 — open na je eigen poging

a) Raw verschil: $\bar{Y}_B - \bar{Y}_A = 60 - 50 = 10$.

b) $\bar{Y}_A^* = 50 - 2 \cdot (4 - 6) = 50 - 2 \cdot (-2) = 50 + 4 = 54$. $\bar{Y}_B^* = 60 - 2 \cdot (8 - 6) = 60 - 4 = 56$.

c) Adjusted verschil: $56 - 54 = 2$. Veel kleiner dan het raw verschil van $10$. Een groot deel van het ruwe verschil zat dus niet in groep-zijn maar in covariaat-zijn: groep B had toevallig een hogere covariaat-waarde, en omdat de covariaat positief samenhangt met Y, lag hun Y vanzelf hoger.

d) $b_w = 0$: geen correctie, adjusted = raw, dus verschil blijft $10$. $b_w = 2.5$: $\bar{Y}_A^* = 50 + 5 = 55$; $\bar{Y}_B^* = 60 - 5 = 55$. Adjusted verschil = $0$ — het hele “groepseffect” verdwijnt na correctie. Dat is de kern: het verschil tussen A en B was eigenlijk een C-verschil.

e) Bias-correctie (of: confounding-correctie). De covariaat confoundeerde het groepseffect; ANCOVA haalt dat eruit.

Inzicht. ANCOVA werkt vóór jou (foutreductie: kleinere noemer, meer power) en tegen misleidende interpretaties (bias-correctie: andere conclusie als groepen op C verschillen). Beide zitten in dezelfde formule.

T6 — $\eta^2_p$ in ANCOVA: wat zit in de noemer?

OpmerkingVraag T6 — Pen-en-papier

etaSquared() print twee kolommen: eta.sq ($\eta^2$, met totale SS in de noemer) en eta.sq.part ($\eta^2_p$, met effect-SS + residu-SS in de noemer).

In een ANCOVA met covariaat $C$ en factor $F$ heb je $SS_C$, $SS_F$ en $SS_w$.

a) Wat is de noemer voor $\eta^2_p$ van de factor — alleen $SS_F + SS_w$, of $SS_C + SS_F + SS_w$?

b) Kan $\sum \eta^2_p$ over alle effecten boven $1$ uitkomen? Waarom (niet)?

c) Voor de bever-data: $SS_F = 2808$, $SS_C = 2135$, $SS_w = 10732$. Bereken $\eta^2_p$ voor zowel de factor als de covariaat.

OpgeletAntwoord T6 — open na je eigen poging

a) Voor $\eta^2_p$ van de factor: alleen $SS_F + SS_w$. Dezelfde logica als bij meerweg-ANOVA: het partial eta-squared sluit de andere effecten uit de noemer.

b) Ja, $\sum \eta^2_p$ kan boven $1$ uitkomen — net als in meerweg-ANOVA. Elk effect “doet alsof” de andere er niet zijn voor zijn eigen noemer; daardoor mag de som de $1$ overschrijden. Dat is geen fout, het past bij de definitie. Bij $\eta^2$ (totale SS in de noemer) blijft de som onder de $1$.

c) Factor: $\eta^2_p = 2808 / (2808 + 10732) = 2808 / 13540 = .21$. Covariaat: $\eta^2_p = 2135 / (2135 + 10732) = 2135 / 12867 = .17$. Past bij etaSquared(m_ancova, type = 3).

Inzicht. $\eta^2_p$ van de factor plus dat van de covariaat is $.21 + .17 = .38$ — ruim boven het gewone $\eta^2$-totaal, $\mathit{\eta}^2_F + \mathit{\eta}^2_C = .14 + .11 = .25$. Logisch: ze rekenen elkaar uit hun eigen noemer.

3.2.c Wat als parallellisme is geschonden?

WaarschuwingSlopes ongelijk? Dan geen ANCOVA

Als de interactie covariaat $\times$ factor wél significant is, betekent dat: het effect van de covariaat verschilt tussen groepen. Of equivalent: het effect van de factor hangt af van de waarde van de covariaat. Dan bestaat er geen één $b_w$, en dus ook geen gecorrigeerd gemiddelde per groep: het verschil tussen de groepen is bij elke waarde van de covariaat anders. Rapporteer dan niet de ANCOVA met $b_w$ en adjusted means — de variabele werkt in dit model niet als covariaat.

OpmerkingVoor wie verder wil — als de hellingen verschillen: een moderatie-effect

Een significante interactie covariaat $\times$ factor is geen “schending van een ANCOVA-aanname die je negeert” maar een substantieve bevinding — een moderatie-effect. De factor modereert de covariaat-Y-relatie. Je kunt dat effect zelf onderzoeken, en dan rapporteer je:

• de interactie-toets ($F$, $df$, $p$, $\eta^2_p$) als hoofdresultaat;
• per-groep slopes (uit het volledige interactie-model, of via emtrends());
• simple slopes of EMM bij specifieke covariaat-waarden (emmeans(model, ~ factor | covariaat, at = list(covariaat = c(...)))).

Het wordt dan eerder een regressie-met-interactie-verhaal dan een ANCOVA-verhaal.

BelangrijkDon’t: ANCOVA op observationele non-equivalente groepen, en daarna causaal praten

Twee groepen die op de covariaat niet vergelijkbaar zijn — en jij zet er ANCOVA op om “te corrigeren” en concludeert dat groep A en B verschillen. Lord (1967) liet zien dat dezelfde data tegengestelde conclusies kan opleveren, afhankelijk van of je naar het verschil-tussen-tijdstippen of de adjusted-post-score kijkt. Beide rekenkundig correct. Welke de “echte” is, hangt af van wat je vooraf veronderstelt — niet van de getallen.

Wie ANCOVA op zulke groepen draait en daarmee een oorzaak claimt, moet eerst uitleggen welke andere verklaringen hij heeft uitgesloten — en hoe. Doe je dat niet — dan rapporteer je geen analyse, dan rapporteer je een wens. T8 hieronder laat zien waarom.

T7 — Aannamechecks specifiek voor ANCOVA

OpmerkingVraag T7 — Pen-en-papier

Naast de “klassieke” ANOVA-aannames (normaliteit residuen, gelijke varianties, onafhankelijkheid) heeft ANCOVA er drie extra:

a) Welke drie?

b) Voor één van de drie bestaat een toets (significantie van de interactie); voor een andere voornamelijk een visuele check. Welke is welke?

c) Wat is de aanname over meetfout in de covariaat, en waarom doet die ertoe? (Hint: bij sterke meetfout in $C$ wordt de schatting van $b_w$ vertekend en de bias-correctie onbetrouwbaar.)

OpgeletAntwoord T7 — open na je eigen poging

a) (1) Een covariaat zonder meetfout. (2) Lineariteit van de covariaat-Y-relatie binnen elke groep. (3) Parallellisme.

b) Parallellisme toets je via de interactie-significantie. Lineariteit check je vooral visueel (scatter Y vs C per groep, residual plot). Sommige methoden gebruiken een toets op de kwadratische term (lm(Y ~ C + I(C^2) + factor)); we doen dat hier niet routinematig.

c) ANCOVA neemt aan dat de covariaat zonder meetfout is gemeten. Bij sterke meetfout in $C$ wordt $b_w$ onderschat (regression-attenuation: meetfout dempt de geschatte slope), en de bias-correctie in de adjusted means is dan onvolledig. Dit is een serieus probleem in psychologisch onderzoek waar covariaten vaak vragenlijst-scores zijn met substantiële meetfout. Oplossingen (latent-variable modellen) vallen buiten dit thema.

T8 — Waar ANCOVA ophoudt: Lord’s paradox en non-equivalent groups

In een observationele studie naar gewichtsverandering bij studenten over een jaar (Frederick Lord, 1967) bekeken twee onderzoekers dezelfde data. De scenario-cijfers:

• Mannen wogen aan begin van het jaar gemiddeld $\bar{C}_{\text{man}} = 80$ kg, vrouwen $\bar{C}_{\text{vrouw}} = 60$ kg.
• Een jaar later: mannen $\bar{Y}_{\text{man}} = 80$ kg, vrouwen $\bar{Y}_{\text{vrouw}} = 60$ kg. Niemand veranderde gemiddeld iets in gewicht.
• Binnen elke groep wel regressie naar het geslachts-gemiddelde: studenten die zwaar boven hun groepsmean begonnen, dreven terug naar dat gemiddelde over het jaar (binnen-groep slope $b_w \approx 0.5$). Die $0.5$ zie je in de plot hieronder terug als de helling van de paarse en groene lijnen — flauwer dan de identiteits-diagonaal ($Y = X$, slope $1$), en juist omdat hij flauwer is kruipen de gevulde driehoeken (adjusted means bij $\bar{C} = 70$) tien kilo uit elkaar.

library(ggplot2)
library(dplyr)

set.seed(42)
n <- 30
maak_lord_data <- function(mean_pre, n, bw) {
  pre  <- rnorm(n, mean_pre, 4)
  post <- mean_pre + bw * (pre - mean_pre) + rnorm(n, 0, 2)
  tibble(pre = pre, post = post)
}

lord <- bind_rows(
  maak_lord_data(80, n, 0.5) %>% mutate(geslacht = "man"),
  maak_lord_data(60, n, 0.5) %>% mutate(geslacht = "vrouw")
)

raw_means <- tibble(
  geslacht = c("man", "vrouw"),
  pre  = c(80, 60),
  post = c(80, 60)
)

# Adjusted post-gewicht bij grand mean pre = 70
adj_means <- tibble(
  geslacht = c("man", "vrouw"),
  pre  = c(70, 70),
  post = c(80 - 0.5 * (80 - 70),
           60 - 0.5 * (60 - 70))
)

ggplot(lord, aes(pre, post, colour = geslacht)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey50") +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_vline(xintercept = 70, linetype = "dotted", colour = "grey50") +
  geom_point(data = raw_means, shape = 5, size = 4, stroke = 1.1) +
  geom_point(data = adj_means, shape = 17, size = 4) +
  scale_colour_manual(values = c(man = "#332288", vrouw = "#117733")) +
  labs(x = "Pre-gewicht (kg, begin jaar)",
       y = "Post-gewicht (kg, einde jaar)") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

Lord’s paradox in beeld. Paars = mannen, groen = vrouwen. Stippellijn = identiteitslijn $Y = X$ (geen verandering). Doorgetrokken lijnen = binnen-groep regressielijnen (slope $b_w \approx 0.5$). Open ruitjes = raw groepsgemiddelden (pre, post). Gevulde driehoeken = adjusted post-gewicht bij grand mean pre ($\bar{C} = 70$). De ruitjes liggen óp de identiteitslijn (geen verandering); de driehoeken liggen tien kilo uit elkaar.

Twee analyses op dezelfde data:

• Onderzoeker 1 vergeleek het verschil in gewichtsverandering ($\Delta = \text{post} - \text{pre}$) tussen mannen en vrouwen. Voor mannen $\bar{\Delta}_{\text{man}} = 0$ kg, voor vrouwen $\bar{\Delta}_{\text{vrouw}} = 0$ kg. Conclusie: geen verschil in gewichtsverandering tussen geslachten.
• Onderzoeker 2 deed ANCOVA met post-gewicht als $Y$ en pre-gewicht als covariaat. Bij grand mean pre ($\bar{C} = 70$ kg): adjusted post-gewicht mannen $= 80 - 0.5 (80 - 70) = 75$ kg; adjusted post-gewicht vrouwen $= 60 - 0.5 (60 - 70) = 65$ kg. Conclusie: wél een verschil in adjusted post-gewicht tussen geslachten — tien kilo.

TipValstrik — “wie had gelijk?”

De voor de hand liggende reflex is om er één goed antwoord uit te kiezen. Dat is de valstrik. Lord’s punt is dat beide kloppen — ze beantwoorden gewoon andere vragen. Welke “echt” is, hangt niet af van de getallen maar van wát je vooraf veronderstelt over de groepen en de covariaat.

OpmerkingVraag T8 — Pen-en-papier

a) Zet de twee conclusies onder elkaar in volledige zinnen (object-altijd-erbij — verschil in wat?). Welke vraag stelt onderzoeker 1, welke onderzoeker 2? Probeer ze als twee verschillende onderzoeksvragen op te schrijven.

b) Onderzoeker 2’s ANCOVA-redenering luidt impliciet: “als mannen en vrouwen even zwaar waren begonnen, zouden ze…”. Bij geslacht is dat een contrafeitelijke aanname — je kunt iemand niet random ander geslacht toewijzen. Wat zegt dat over de causale interpretatie van Onderzoeker 2’s conclusie?

c) In de bever-data zijn de drie werkstijlen niet gerandomiseerd — een bever pakt zélf zijn dammen op een bepaalde manier aan, en familieleden besluiten zelf of ze wel of niet meehelpen. Wat betekent dat voor de causale interpretatie van “effect van werkstijl op damtevredenheid, na correctie voor familie-hulp”?

d) Wanneer is ANCOVA op een pre-test-covariaat wel een legitieme analyse, zonder Lord-paradox-zorgen?

OpgeletAntwoord T8 — open na je eigen poging

a) Twee volledige zinnen:

• Verschil in gewichtsverandering ($\Delta = \text{post} - \text{pre}$) tussen mannen en vrouwen is nul.
• Verschil in adjusted post-gewicht (bij gelijke pre-waarde van $70$ kg) tussen mannen en vrouwen is tien kilo.

De onderzoeksvragen erachter:

• Onderzoeker 1: “Verandert het gewicht binnen mannen anders dan binnen vrouwen over dit jaar?” — empirisch te beantwoorden zonder veronderstellingen.
• Onderzoeker 2: “Als mannen en vrouwen even zwaar waren begonnen, zouden ze nu even zwaar zijn?” — vereist een contrafeitelijke aanname.

De twee zinnen lijken op elkaar, maar verschillen op het cruciale object van de vergelijking. Verschil in $\Delta$ versus verschil in adjusted $Y$ bij gelijke $C$. De eerste is een feit; de tweede een gedachte-experiment.

b) Onderzoeker 2’s vraag is bij geslacht causaal onzinnig — je kunt iemand niet random ander geslacht geven en kijken wat er gebeurt. ANCOVA “corrigeert” hier rekenkundig wel, maar interpreteert iets dat in de werkelijkheid niet bestaat: een hypothetische mens-die-evenzwaar-begint-maar-een-ander-geslacht-heeft. Dat is Lord’s paradox: de paradox is geen rekenfout, het is een herinnering dat een covariaat-keuze altijd een causale veronderstelling impliceert, en als die veronderstelling niet houdt, beantwoordt ANCOVA misschien een vraag zonder zinnig antwoord.

c) In de bever-data is bouwwoede niet gerandomiseerd. Een bever kiest (of wordt door zijn natuur gestuurd naar) een werkstijl, en familieleden besluiten zelf of ze meehelpen. De ANCOVA-conclusie “er is een effect van werkstijl op damtevredenheid, na correctie voor familie-hulp” is daarmee een uitspraak over adjusted means in deze observationele steekproef — niet automatisch een uitspraak over wat een willekeurig in een andere werkstijl geplaatste bever zou ervaren. Voor causale claims helpt randomisatie, of een gerichte causale analyse (DAG’s, propensity score matching, instrumentele variabelen — buiten dit oefenboek).

d) Bij gerandomiseerde experimenten zijn groepen op verwachting gelijk op alle pre-treatment-variabelen — inclusief de pre-test. ANCOVA met die pre-test als covariaat is dan een legitieme power-boost: je vermindert restvariantie zonder dat de bias-correctie iets onzinnigs vooronderstelt (de groepen zijn random toegewezen, de adjusted-means-vergelijking is causaal zinvol). Het beste anti-Lord-paradox-recept is dus: randomiseer waar je kunt. Waar dat niet kan, blijf scherp op wat je covariaat-keuze impliciet veronderstelt — en wees terughoudend met causale taal in de conclusie.

Inzicht. Lord’s paradox is geen rekenfout; het is een herinnering dat statistiek niet zonder onderzoeksontwerp werkt. Welke covariaat je kiest — en welke niet — is een theoretische beslissing met empirische consequenties. Bij elke ANCOVA op niet-gerandomiseerde groepen geldt: rapporteer wat je modelleert (adjusted means) en wees voorzichtig met wat je concludeert (causale claims).

T9 — Meerdere covariaten

OpmerkingVraag T9 — Pen-en-papier

ANCOVA kan natuurlijk uitgebreid worden met meerdere covariaten: $Y \sim C_1 + C_2 + \text{factor}$.

a) Hoe verandert de df-tabel? Bij twee covariaten in plaats van één?

b) Bij meerdere covariaten is multicollineariteit tussen $C_1$ en $C_2$ een nieuw probleem. Welke maat ken je daarvoor uit thema 1 (MRA)?

c) Bedenk een bever-uitbreiding: de bever heeft naast meedoende_familie óók bouwtijd (uren) bijgehouden. Hoe zou je dat in het model verwerken? Welk effect verwacht je (foutreductie, bias-correctie, of beide)?

OpgeletAntwoord T9 — open na je eigen poging

a) Met twee covariaten: $df_{C_1} = 1$, $df_{C_2} = 1$, $df_F = I - 1$, $df_w = N - I - 2$, $df_{\text{tot}} = N - 1$. Algemeen voor $K$ covariaten: $df_w = N - I - K$.

b) $\text{VIF}_j$ — de variance inflation factor voor elke covariaat ($\text{VIF}_j > 5$ of $> 10$ als veel-gebruikte vuistregels). Of equivalent de tolerance $T_j = 1/\text{VIF}_j$.

c) Model: lm(damtevredenheid ~ bouwwoede + meedoende_familie + bouwtijd). Als bouwtijd nog wat unieke Y-variantie verklaart, krijg je extra foutreductie (kleinere $SS_w$, hogere power op de factor-toets). Als bouwtijd ook tussen werkstijlen verschilt (perfectionisten doen er waarschijnlijk veel langer over), krijg je extra bias-correctie. Beide tegelijk is normaal.

3.3 Cartografie — zes scenes uit bever-land

Zes hypothetische bever-onderzoeken. Telkens twee kolonies, telkens damtevredenheid en familie-hulp — maar telkens iets anders eronder.

“Vorig jaar,” zei de bever, “ben ik bij collega’s wezen kijken in andere dammen. Andere kolonies, andere familie-regels, andere damrijkheid. Soms hielp familie-hulp om damtevredenheid te begrijpen. Soms helemaal niet. En in één geval — let op — keerde de conclusie zich helemaal om.”

Tot nu toe heb je ANCOVA gezien op één kolonie. Maar één resultaat is geen ANCOVA-vermogen. Wie de mechaniek wil leren, moet de techniek zien werken in zes verschillende werelden: van een kolonie waar de covariaat niks doet, tot een kolonie waar de covariaat het hele groepseffect verslindt.

Hieronder zes mini-onderzoeken. Telkens twee bever-kolonies — A en B. Telkens damtevredenheid ($Y$) als uitkomst, familie-hulp ($C$) als covariaat. De groepsgemiddelden $\bar{Y}$ en $\bar{C}$ verschillen per scenario, evenals de pooled-within slope $b_w$.

Zo lees je een gecorrigeerd gemiddelde af in zo’n plot. Begin bij het open ruitje van een groep: dat staat op $(\bar{C}_j, \bar{Y}_j)$. Volg vandaar de gestreepte lijn van die groep tot de verticale stippellijn bij $\bar{C}$. De hoogte waar je uitkomt, is $\bar{Y}_j^*$ — daar staat de gevulde driehoek. Dat is de formule $\bar{Y}_j^* = \bar{Y}_j - b_w(\bar{C}_j - \bar{C})$ in beeld: de lijn heeft helling $b_w$, en je schuift over de afstand $\bar{C}_j - \bar{C}$. In het antwoord bij 3.3 staat het voor één paneel helemaal uitgeschreven; daar heet het de vinger-volgt-lijn.

# Per scenario: vier vaste groepsgemiddelden + slope, simulatie rond die centra.
scenarios <- tribble(
  ~scene, ~titel,             ~y_a, ~y_b, ~c_a, ~c_b, ~bw,
  1,      "1. Identiek",        55,   55,    6,    6,  2.0,
  2,      "2. Puur effect",     50,   70,    6,    6,  2.0,
  3,      "3. Suppressie",      65,   55,    4,    8,  2.0,
  4,      "4. Demping",         50,   65,    4,    8,  2.0,
  5,      "5. Simpson",         70,   65,   10,    5,  2.0,
  6,      "6. Confound",        50,   60,    4,    8,  2.5
)

set.seed(303)
maak_data <- function(rij, n = 25) {
  ca <- rnorm(n, rij$c_a, 1.2)
  cb <- rnorm(n, rij$c_b, 1.2)
  ya <- rij$y_a + rij$bw * (ca - rij$c_a) + rnorm(n, 0, 4)
  yb <- rij$y_b + rij$bw * (cb - rij$c_b) + rnorm(n, 0, 4)
  tibble(
    titel = rij$titel,
    Y     = c(ya, yb),
    C     = c(ca, cb),
    Groep = rep(c("A", "B"), each = n)
  )
}

scenes <- bind_rows(lapply(seq_len(nrow(scenarios)),
                            function(i) maak_data(scenarios[i, ])))

# Per scene: grand mean C, raw groepsmiddens, adjusted means.
markers <- scenarios %>%
  mutate(
    cm      = (c_a + c_b) / 2,
    ya_star = y_a - bw * (c_a - cm),
    yb_star = y_b - bw * (c_b - cm)
  )

raw_means <- bind_rows(
  markers %>% transmute(titel, Groep = "A", C = c_a, Y = y_a),
  markers %>% transmute(titel, Groep = "B", C = c_b, Y = y_b)
)

adj_means <- bind_rows(
  markers %>% transmute(titel, Groep = "A", C = cm, Y = ya_star),
  markers %>% transmute(titel, Groep = "B", C = cm, Y = yb_star)
)

vlines <- markers %>% select(titel, cm)

# De gestreepte lijnen: per groep de helling b_w uit het ontwerp, door het ruwe
# groepsgemiddelde, over het stuk C waar die groep zit. Dan staat de gevulde
# driehoek precies op zijn lijn, zoals de vinger-volgt-lijn hieronder belooft.
# Een lm()-lijn per groep zou b_w schatten uit 25 punten met ruis, en komt dan
# overal tussen ongeveer 0.3 en 3.4 uit.
binnen <- scenes %>%
  group_by(titel, Groep) %>%
  summarise(c_min = min(C), c_max = max(C), .groups = "drop") %>%
  left_join(bind_rows(
    markers %>% transmute(titel, Groep = "A", cj = c_a, yj = y_a, bw),
    markers %>% transmute(titel, Groep = "B", cj = c_b, yj = y_b, bw)
  ), by = c("titel", "Groep")) %>%
  mutate(y_min = yj + bw * (c_min - cj), y_max = yj + bw * (c_max - cj))

ggplot(scenes, aes(C, Y, colour = Groep)) +
  geom_point(alpha = 0.4, size = 1.4) +
  geom_segment(data = binnen, aes(x = c_min, xend = c_max, y = y_min, yend = y_max),
               linetype = "dashed", linewidth = 0.6) +
  geom_smooth(aes(group = 1), method = "lm", se = FALSE,
              colour = "black", linewidth = 0.6) +
  geom_vline(data = vlines, aes(xintercept = cm),
             linetype = "dotted", colour = "grey50") +
  geom_point(data = raw_means, shape = 5, size = 3.0, stroke = 0.9) +
  geom_point(data = adj_means, shape = 17, size = 3.0) +
  scale_colour_manual(values = c(A = "#332288", B = "#117733")) +
  # Alle zes op dezelfde assen: dan ziet dezelfde helling b_w = 2.0 er in elke
  # scene even steil uit, en de b_w = 2.5 van scene 6 ook echt steiler.
  facet_wrap(~ titel, ncol = 3) +
  labs(x = "C (familie-hulp)", y = "Y (damtevredenheid)") +
  theme_minimal(base_size = 10) +
  theme(legend.position = "bottom",
        strip.text = element_text(face = "bold"))

Zes scenes uit bever-land. Gekleurde punten = observaties (paars = kolonie A, groen = kolonie B); gestreepte lijnen = binnen-groep regressielijnen (slope $b_w$); doorgetrokken zwarte lijn = totale regressielijn (slope $b_{\text{tot}}$). Open ruitjes = raw groepsgemiddelden bij eigen $\bar{C}_j$; gevulde driehoeken = adjusted means bij grand-mean $\bar{C}$ (verticale stippellijn).

OpmerkingVraag 3.3 — Sudoku van bever-land

Bestudeer de zes scatters. Onderstaande tabel toont per scene de raw groepsgemiddelden ($\bar{Y}_A$, $\bar{Y}_B$, $\bar{C}_A$, $\bar{C}_B$) en de pooled-within slope $b_w$. Vul de overige zes kolommen in. De grand mean $\bar{C}$ is steeds het ongewogen gemiddelde van $\bar{C}_A$ en $\bar{C}_B$ (gelijke groepsgrootte).

Formules ter herinnering:

$$\bar{Y}_j^* = \bar{Y}_j - b_w (\bar{C}_j - \bar{C}), \qquad \Delta_{\text{raw}} = \bar{Y}_B - \bar{Y}_A, \qquad \Delta_{\text{adj}} = \bar{Y}_B^* - \bar{Y}_A^*.$$

Zes scenes in bever-land: wat je krijgt, en wat je zelf uitrekent
Scene	$\bar{Y}_A$	$\bar{Y}_B$	$\bar{C}_A$	$\bar{C}_B$	$b_w$	$\bar{C}$	$\bar{Y}_A^*$	$\bar{Y}_B^*$	$\Delta_{\text{raw}}$	$\Delta_{\text{adj}}$
1	55	55	6	6	2.0	?	?	?	?	?
2	50	70	6	6	2.0	?	?	?	?	?
3	65	55	4	8	2.0	?	?	?	?	?
4	50	65	4	8	2.0	?	?	?	?	?
5	70	65	10	5	2.0	?	?	?	?	?
6	50	60	4	8	2.5	?	?	?	?	?
Note. $Y$ = damtevredenheid, $C$ = familie-hulp, A en B = de twee kolonies; $b_w$ = de helling binnen de groepen, een * = aangepast gemiddelde, $\bar{C}$ = het gemiddelde van $\bar{C}_A$ en $\bar{C}_B$. Elke $\Delta$ is B min A: ruw (raw) en na de correctie (adj). Elk ? reken je zelf uit.

Vervolgens — uit de tabel die je net hebt ingevuld:

a) In welke scenes is er systematische bias op de covariaat (groepen verschillen op $\bar{C}$)? In welke richting?

b) In welke scene heeft ANCOVA geen meerwaarde ten opzichte van een gewone $t$-toets, en waarom?

c) In welke scene dempt de bias het ruwe effect (adjusted verschil kleiner dan raw)? In welke scene versterkt de bias het effect (adjusted verschil groter)?

d) In welke scene keert het teken van $\Delta$ om na correctie? (Dat heet ook wel Simpson’s paradox.)

e) In welke scene verklaart de covariaat het hele groep-effect (adjusted verschil $\approx 0$ terwijl raw $\neq 0$)?

f) Vergelijk de positie van de open ruitjes (raw means) met de gevulde driehoeken (adjusted means bij $\bar{C}$) in elke plot. Welke scenes laten zien dat de adjusted means dichter bij elkaar liggen dan de raw means? Welke verder uit elkaar?

OpgeletAntwoord 3.3 — open na je eigen poging

Tabel ingevuld.

Zes scenes in bever-land, ingevuld
Scene	$\bar{Y}_A$	$\bar{Y}_B$	$\bar{C}_A$	$\bar{C}_B$	$b_w$	$\bar{C}$	$\bar{Y}_A^*$	$\bar{Y}_B^*$	$\Delta_{\text{raw}}$	$\Delta_{\text{adj}}$
1	55	55	6	6	2.0	6.0	55	55	0	0
2	50	70	6	6	2.0	6.0	50	70	20	20
3	65	55	4	8	2.0	6.0	69	51	−10	−18
4	50	65	4	8	2.0	6.0	54	61	15	7
5	70	65	10	5	2.0	7.5	65	70	−5	5
6	50	60	4	8	2.5	6.0	55	55	10	0
Note. $Y$ = damtevredenheid, $C$ = familie-hulp, A en B = de twee kolonies; $b_w$ = de helling binnen de groepen, een * = aangepast gemiddelde, $\bar{C}$ = het gemiddelde van $\bar{C}_A$ en $\bar{C}_B$. Elke $\Delta$ is B min A: ruw (raw) en na de correctie (adj).

a) Scenes 3, 4, 5 en 6 hebben bias. Richtingen: scene 3: $\bar{C}_A < \bar{C}_B$ (kolonie B had hogere covariaat — terwijl A juist hogere Y had); scene 4: $\bar{C}_A < \bar{C}_B$ (B hogere C, in dezelfde richting als Y); scene 5: $\bar{C}_A \gg \bar{C}_B$ (A had veel hogere C); scene 6: $\bar{C}_A < \bar{C}_B$.

b) Scene 1 — geen groep-effect, geen bias. Een $t$-toets zou ook geen verschil vinden. Scene 2 — wel groep-effect maar geen bias, dus de correctie laat het verschil ongemoeid; $t$-toets en ANCOVA komen praktisch op hetzelfde uit (ANCOVA wint alleen op power als $b_w > 0$ binnen elke groep — foutreductie via een kleinere $MS_w$).

c) Demping: scene 4 (raw 15 → adjusted 7) en scene 6 (raw 10 → adjusted 0). Versterking: scene 3 (raw $-10$ → adjusted $-18$, absoluut groter) — bias liep tégen het groep-effect in. Dit is suppressie: ANCOVA legt het echte groep-effect bloot dat door de covariaat verstopt werd.

d) Scene 5. Raw: $\bar{Y}_A > \bar{Y}_B$ ($\Delta_{\text{raw}} = -5$ ten opzichte van B-min-A, dus A > B). Adjusted: $\bar{Y}_A^* < \bar{Y}_B^*$ ($\Delta_{\text{adj}} = +5$, dus B > A). Het teken keert om — Simpson’s paradox in ANCOVA-vorm. Reden: groep A had veel hogere C, en omdat $b_w$ positief is, krijgt A een grote neerwaartse correctie en B een grote opwaartse correctie.

e) Scene 6. Raw verschil 10, adjusted verschil 0 — de hele $Y$-discrepantie tussen kolonie A en B zat in hun verschil op familie-hulp, niet in de kolonie-zelf. Dit is het klassieke confounding-voorbeeld; in echte data zouden we hier moeten concluderen dat “kolonie” geen toegevoegde verklaarwaarde heeft gegeven familie-hulp.

f) Adjusted dichter bij elkaar dan raw: scenes 4 en 6. Adjusted verder uit elkaar dan raw: scenes 3 en 5. Scenes 1 en 2: even ver uit elkaar (de correctie verschuift nul of niets).

Inzicht — de zes scenes als kaart.

Scene	Bias?	Effect na correctie	Naam in de literatuur
1	nee	niets verandert; niets te zien	basislijn
2	nee	gelijk aan raw effect	pure foutreductie
3	ja	versterkt	suppressor
4	ja	verzwakt	partiële confounding
5	ja	keert om	Simpson’s paradox
6	ja	verdwijnt	volledige confounding

Wie ANCOVA als instrument heeft, kan in elke kolonie aankomen en eerst de kaart maken — welk scene tref ik aan? — vóór hij conclusies trekt. Dat is het verschil tussen techniek toepassen en techniek begrijpen.

Vinger-volgt-lijn — paneel 4 (Demping) als voorbeeld. Zoek paneel 4 op in de figuur. Het paarse open ruitje (kolonie A’s raw mean) staat op $(\bar{C}_A, \bar{Y}_A) = (4, 50)$. Schuif nu met je vinger horizontaal naar de stippellijn op $\bar{C} = 6$ (de grand mean). Volg dan de gestreepte paarse lijn omhoog volgens slope $b_w = 2$: $50 + 2 \cdot (6 - 4) = 54$. Daar staat de paarse gevulde driehoek — de adjusted mean. Hetzelfde voor groen: ruitje op $(8, 65)$, schuif naar $C = 6$, volg de lijn omlaag: $65 - 2 \cdot (8 - 6) = 61$. Driehoek op $(6, 61)$. Het raw verschil ($65 - 50 = 15$) wordt na correctie het adjusted verschil ($61 - 54 = 7$) — demping zichtbaar. Eén scene zo doorgewerkt, de andere vijf werken op dezelfde manier — alleen de richting en de grootte van de schuif verschillen.

Wat de doorgetrokken zwarte lijn vertelt. Dat is de totale regressielijn (alle punten samen, groep genegeerd, $b_{\text{tot}}$). In de meeste scenes loopt hij ongeveer parallel aan de gestreepte binnen-groep-lijnen — saai. Maar in paneel 3 (Suppressie) loopt hij juist de andere kant op: gestreepte lijnen omhoog, zwarte lijn omlaag. Dáár zie je visueel hoe de totale lijn iets anders kan vertellen dan de lijnen binnen de groepen — alsof je twee verhalen tegelijk vertelt en ze elkaar tegenspreken. In paneel 6 (Confound) liggen zwart en gestreept boven op elkaar — covariaat verklaart alles, geen unieke groep-effect over.

OpmerkingVraag 3.3 g — een zevende scene waar ANCOVA niet meer past

Stel je krijgt een 7e scene waar de twee binnen-groep regressielijnen duidelijk verschillende hellingen hebben — de gestreepte lijnen van kolonie A en B kruisen elkaar over het zichtbare C-bereik in plaats van parallel te lopen.

Waarom is ANCOVA hier niet de juiste keuze? Denk aan hoe je een gecorrigeerd gemiddelde afleest (vanaf het ruitje langs de gestreepte lijn naar de stippellijn): wat gebeurt er met het verschil tussen de twee gevulde driehoeken als je de stippellijn een stukje naar links of naar rechts schuift?

OpgeletAntwoord 3.3 g — open na je eigen poging

Verschillende binnen-groep slopes = parallellisme-schending = significante interactie covariaat $\times$ factor. De ANCOVA-aanname dat één gepoolde $b_w$ bestaat, geldt dan niet — het effect van de covariaat hangt af van de groep.

Schuif je de stippellijn op, dan verandert het verschil tussen de driehoeken mee: links van het kruispunt ligt A boven B, rechts ervan B boven A. Er is dus niet één gecorrigeerd verschil tussen de kolonies, maar bij elke waarde van $C$ een ander. Daarmee valt de vraag weg die ANCOVA beantwoordt. Conclusie: geen ANCOVA — de variabele werkt hier niet als covariaat.

OpmerkingVoor wie verder wil — een tweede covariaat in R

“Vorig jaar,” zei de bever, “noteerde ik er ook bij hoeveel uur ik aan elke dam werkte. Misschien helpt dat ook.”

Een tweede covariaat past in dezelfde mal — in deze data hebben we geen bouwtijd-kolom, maar het patroon werkt voor elke continue covariaat:

options(contrasts = c("contr.sum", "contr.poly"))   # vóór lm()
m_uitbreiding <- lm(damtevredenheid ~ bouwwoede + meedoende_familie + bouwtijd,
                    data = bever_tevredenheid_uitgebreid)
Anova(m_uitbreiding, type = 3)
emmeans(m_uitbreiding, ~ bouwwoede)

Wat verandert er, en wat blijft hetzelfde?

• Adjusted means worden nu berekend bij de gemiddelde $\bar{C}_1$ en $\bar{C}_2$ tegelijk.
• df dalen één extra: $df_w = N - I - K$ met $K$ covariaten.
• Multicollineariteit tussen covariaten kan optreden: check VIF.
• Parallellisme-aanname: nu over beide covariaat-factor-interacties.
• De rest (Type-III SS, EMM, paarsgewijze vergelijking, $\eta^2_p$) blijft hetzelfde principe.

R-spiekblad bij ANCOVA

Alle commando’s op één plek

Pakketten en data laden

# Pakketten activeren — één keer per sessie.
library(car)         # voor Anova(., type = 3)
library(lsr)         # voor etaSquared()
library(emmeans)     # voor adjusted means rechtstreeks uit het model
library(rstatix)     # voor emmeans_test(): paarsgewijs, met Bonferroni
library(dplyr)       # alleen als je de dplyr-schrijfwijze gebruikt

# Type III SS vereist sum-to-zero contrasten op factoren.
# Bovenaan, vóór je eerste lm(): R legt de codering vast bij lm().
options(contrasts = c("contr.sum", "contr.poly"))

# Hoofd-dataset laden — object 'bever_tevredenheid' verschijnt vanzelf.
load("data/bever_tevredenheid.RData")
str(bever_tevredenheid)

# Mini-dataset voor E5.
load("data/otter_huishouden.RData")

Verkenning en vuistregels — aantallen, gemiddelden en standaardafwijkingen

Base R

# Aantallen per factor-niveau.
table(bever_tevredenheid$bouwwoede)

# Raw groepsgemiddelden (Y en covariaat).
aggregate(damtevredenheid   ~ bouwwoede, data = bever_tevredenheid, FUN = mean)
aggregate(meedoende_familie ~ bouwwoede, data = bever_tevredenheid, FUN = mean)

# Standaardafwijkingen per groep, voor de vuistregel s_max / s_min < 2.
aggregate(damtevredenheid   ~ bouwwoede, data = bever_tevredenheid, FUN = sd)

dplyr

# Per werkstijl: aantal, gemiddelde Y, gemiddelde covariaat en SD van Y.
# Eigen namen links van de =: heet de eerste weer 'damtevredenheid', dan
# rekent sd() daarna met dat ene gemiddelde en krijg je NA.
bever_tevredenheid |>
  group_by(bouwwoede) |>
  summarise(n               = n(),
            M_tevredenheid  = mean(damtevredenheid),
            M_familie       = mean(meedoende_familie),
            SD_tevredenheid = sd(damtevredenheid))

Verkenning in beeld — boxplot en scatter

# Boxplot van Y per factor-niveau.
boxplot(damtevredenheid ~ bouwwoede, data = bever_tevredenheid)

# Scatter Y vs covariaat, lijn per groep — visuele parallellism-check.
library(ggplot2)
ggplot(bever_tevredenheid,
       aes(x = meedoende_familie, y = damtevredenheid, colour = bouwwoede)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

Parallellisme-check (interactie factor × covariaat)

# Model met interactie — toetst homogeniteit van regressielijnen.
# De contrasten staan al bovenaan op contr.sum.
m_par <- lm(damtevredenheid ~ bouwwoede * meedoende_familie, data = bever_tevredenheid)
Anova(m_par, type = 3)
# Lees alleen de regel van de interactie.
# Interactie niet significant => parallellism houdbaar => ANCOVA mag.

ANCOVA-model fitten — factor eerst, covariaat daarna

# Het additieve ANCOVA-model — factor eerst, covariaat daarna.
m_ancova <- lm(damtevredenheid ~ bouwwoede + meedoende_familie,
               data = bever_tevredenheid)

# b_w, de helling binnen de groepen: de rij van de covariaat.
summary(m_ancova)

Type-III SS met car::Anova en sum-to-zero contrasten

# ANCOVA-tabel met F-toets per term. Werkt alleen goed als
# options(contrasts = c("contr.sum", "contr.poly")) vóór lm() stond.
Anova(m_ancova, type = 3)

Adjusted means en paarsgewijze vergelijking

# Paarsgewijze vergelijking van de adjusted means, met Bonferroni-correctie.
# Lees de kolom p.adj.
em <- emmeans_test(data = bever_tevredenheid, damtevredenheid ~ bouwwoede,
                   covariate = meedoende_familie, p.adjust.method = "bonferroni")
em

# De adjusted means (EMM) per niveau van de factor.
get_emmeans(em)

Effectgrootte: partial eta-squared

# eta-squared en partial eta-squared per term.
etaSquared(m_ancova, type = 3)

Diagnostiek — vuistregels, residual en Cook’s distance

# De vuistregels lees je af uit de aantallen en SD's bij Verkenning:
# n >= 15 per groep, n_max / n_min <= 1.5, s_max / s_min < 2.

# Residual en Cook's distance plots.
plot(m_ancova, which = 1)   # Residuen vs voorspeld
plot(m_ancova, which = 4)   # Cook's distance per observatie

Plot — regressielijnen per groep en adjusted means

library(ggplot2)

# Per-groep regressielijnen — visualiseer parallellism.
ggplot(bever_tevredenheid,
       aes(x = meedoende_familie, y = damtevredenheid, colour = bouwwoede)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

# Adjusted means met 95%-CI.
plot(emmeans(m_ancova, ~ bouwwoede))

Klaar — codering terug

# Klaar met ANCOVA? Zet de codering terug op R's standaard, anders geldt
# contr.sum ook in je volgende regressie.
options(contrasts = c("contr.treatment", "contr.poly"))

Voorbeeld-tentamenvragen

Even oefenen op tentamen-toon

Onderaan dit hoofdstuk staan vier korte theorievragen in tentamen-stijl en één mini R-opdracht. Geen verhaal-frame meer: hier oefen je zoals op het tentamen. Hou de vuistregels uit dit hoofdstuk paraat: $\eta^2_p$ lees je als een proportie van de variantie die na de andere term overblijft; parallellisme-aanname via $p$-waarde van interactie covariaat $\times$ factor; bias-correctie als groepen op de covariaat verschillen.

Theorie en handreken

OpmerkingVraag E1 — Adjusted versus raw means uit output

Een eekhoorn vergelijkt drie nest-types op nestlocatie-veiligheid (score $0$–$10$), met de continue covariaat boomhoogte (m). Hij rapporteert:

R
              raw_M   adj_M     n
boomhol       6.20    6.85     22
takkennest    7.10    7.05     19
holboom       8.40    7.80     21

Welke uitspraak past het best?

1. De raw-mean en de adjusted-mean zijn identiek; ANCOVA heeft geen verschil gemaakt.
2. De holboom-groep had een hogere boomhoogte dan gemiddeld; daarom is hun adjusted-mean lager dan hun raw-mean.
3. De boomhol-groep had een hogere boomhoogte dan gemiddeld; daarom ligt hun adjusted-mean hoger.
4. Adjusted means hebben niets met de covariaat te maken; ze corrigeren alleen voor onbalans in $n$.

OpgeletAntwoord E1 — open na je eigen poging

b) $\bar{Y}_j^* = \bar{Y}_j - b_w(\bar{C}_j - \bar{C})$. De holboom raw $>$ adjusted ($8.40 > 7.80$): de groep zat boven het overall covariaat-gemiddelde, dus correctie omlaag. De boomhol raw $<$ adjusted ($6.20 < 6.85$): die zat onder gemiddeld, dus correctie omhoog. Optie a is feitelijk fout, c heeft de richting verkeerd, d verwart adjusted means met onbalans-correctie.

OpmerkingVraag E2 — $b_w$ interpreteren uit output

Bij een ANCOVA op vissoortenrijkdom (aantal soorten/poel) met poelopp (m²) als covariaat en poel-type (3 niveaus) als factor leest een ecoloog dit fragment.

R
Coefficients:
              Estimate  Std. Error  t value  Pr(>|t|)
(Intercept)    11.40       1.20      9.50    < .001
poelopp         0.045      0.011     4.09    < .001
poel_type1     -1.30       0.78     -1.67    .098
poel_type2      0.40       0.75      0.53    .597

Wat is de meest verdedigbare interpretatie van de coëfficiënt $b_w = 0.045$?

1. Per extra m² oppervlakte stijgt de vissoortenrijkdom met $0.045$ soorten, gemiddeld over de drie poel-types (pooled-within slope).
2. De gemiddelde poelopp is $0.045$ m².
3. $0.045$ is de standaardafwijking van de covariaat.
4. De covariaat verklaart $4.5\%$ van de variantie.

OpgeletAntwoord E2 — open na je eigen poging

a) $b_w$ in een ANCOVA-output is de pooled-within slope: het gemiddelde effect van de covariaat op Y, binnen de groepen. Optie b verwart coëfficiënt met gemiddelde, c verwart met SD, d verwart met $\eta^2$.

OpmerkingVraag E3 — Parallellisme-check beoordelen

Een vleermuis vergelijkt drie jachtgebieden op vangstsucces, met luchttemperatuur als covariaat. Voor de parallellisme-check rapporteert hij:

R
Anova Table (Type III tests)

                          Sum Sq  Df  F value  Pr(>F)
(Intercept)               145.30   1  98.12    < .001
temperatuur                42.10   1  28.43    < .001
jachtgebied                 5.20   2   1.76    .180
temperatuur:jachtgebied     2.90   2   0.98    .378
Residuals                 152.40 103

Wat is de juiste conclusie voor de parallellisme-aanname?

1. De aanname is geschonden ($F$ van temperatuur is significant).
2. De aanname is houdbaar; de interactie temperatuur $\times$ jachtgebied is niet significant ($p = .38$).
3. De aanname is niet te toetsen zonder Levene’s toets.
4. De aanname is per definitie alleen geldig als ook het hoofdeffect van de factor significant is.

OpgeletAntwoord E3 — open na je eigen poging

b) De parallellisme-aanname wordt getoetst via de interactie covariaat $\times$ factor. $F(2, 103) = 0.98$, $p = .38$ — niet significant, dus geen bewijs tegen gelijke slopes. Conclusie: de aanname is houdbaar, ANCOVA mag worden uitgevoerd. Levene gaat over gelijke varianties van de residuen, niet over slopes.

OpmerkingVraag E4 — $\eta^2_p$ aflezen en interpreteren

Een uil rapporteert deze ANCOVA-tabel:

R
Anova Table (Type III tests)

Response: prooi_aantal
              Sum Sq Df  F value   Pr(>F)
(Intercept)   1820.5  1  214.7    < .001
maan_helderheid 412.3  1   48.6    < .001
gebied          78.6  2    4.6     .013
Residuals      678.4 80

Welke partial eta-squared waarden horen bij deze tabel?

1. Maan-helderheid: $\eta^2_p = .38$; gebied: $\eta^2_p = .10$.
2. Maan-helderheid: $\eta^2_p = .35$; gebied: $\eta^2_p = .07$.
3. Maan-helderheid: $\eta^2_p = .61$; gebied: $\eta^2_p = .12$.
4. Maan-helderheid: $\eta^2_p = .10$; gebied: $\eta^2_p = .38$.

OpgeletAntwoord E4 — open na je eigen poging

a) Voor maan-helderheid: $\eta^2_p = 412.3 / (412.3 + 678.4) = 412.3 / 1090.7 = .38$. Voor gebied: $\eta^2_p = 78.6 / (78.6 + 678.4) = 78.6 / 757.0 = .10$. Optie b deelt door alle kwadratensommen samen: dat is de gewone $\eta^2$, niet de partiële. Optie c vergeet het effect zelf in de noemer en deelt alleen door $SS_{\text{residu}}$. Optie d verwisselt de twee effecten.

R-practical opdrachtje

Het huishouden van de otter

OpmerkingVraag E5 — Mini-ANCOVA bij de otter

Een onderzoeker bekijkt de huishoudtevredenheid (score $0$–$100$) van een otter in een relatie met drie typen partner-stijlen — dagslaper, nachtbraker en wisselend — met de gedeelde tijd (uren samen per week) als continue covariaat. Bij $78$ koppels werd zowel de tevredenheid als de gedeelde tijd vastgelegd. De dataset staat in data/otter_huishouden.RData en bevat het object otter_huishouden. Sla je R-commando’s op in één scriptbestand: otter.R. Gebruik $\alpha = .05$.

a) Toets de parallellisme-aanname met lm(... * ...) en Anova(., type = 3). Is de aanname houdbaar? Rapporteer de toetsstatistiek, de vrijheidsgraden en de p-waarde van de interactie.

b) Voer de ANCOVA uit met factor eerst en covariaat daarna (lm(huishoudtevredenheid ~ partner_stijl + gedeelde_tijd)). Welke effecten zijn significant? Rapporteer de toetsstatistiek, de vrijheidsgraden en de p-waarde per effect.

c) Welk effect heeft de grootste $\eta^2_p$? Hoeveel procent van de variantie die na de andere term overblijft, verklaart het?

d) Vergelijk de partner-stijlen paarsgewijs met emmeans_test() (Bonferroni-correctie) en haal de adjusted means op met get_emmeans(). Welke partner-stijlen verschillen significant na correctie? Schrijf één APA-zin.

OpgeletAntwoord E5 — open na je eigen poging

load("data/otter_huishouden.RData")
str(otter_huishouden)

'data.frame':   78 obs. of  3 variables:
 $ huishoudtevredenheid: int  48 41 62 59 63 59 72 61 73 40 ...
 $ gedeelde_tijd       : num  19.3 13 20.9 18.8 19.3 17.2 18.9 24.7 17.9 9.9 ...
 $ partner_stijl       : Factor w/ 3 levels "dagslaper","nachtbraker",..: 1 1 1 1 1 1 1 1 1 1 ...

# Parallellism-check. Codering eerst, vóór lm().
options(contrasts = c("contr.sum", "contr.poly"))
m_par_otter <- lm(huishoudtevredenheid ~ partner_stijl * gedeelde_tijd, data = otter_huishouden)
Anova(m_par_otter, type = 3)

Anova Table (Type III tests)

Response: huishoudtevredenheid
                            Sum Sq Df  F value          Pr(>F)    
(Intercept)                 9126.3  1 114.5784       < 2.2e-16 ***
partner_stijl                  1.8  2   0.0110          0.9890    
gedeelde_tijd               4231.6  1  53.1268 0.0000000003232 ***
partner_stijl:gedeelde_tijd   59.9  2   0.3761          0.6879    
Residuals                   5734.9 72                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ANCOVA-model — factor eerst, covariaat daarna.
m_otter <- lm(huishoudtevredenheid ~ partner_stijl + gedeelde_tijd, data = otter_huishouden)
Anova(m_otter, type = 3)

Anova Table (Type III tests)

Response: huishoudtevredenheid
              Sum Sq Df  F value           Pr(>F)    
(Intercept)   9134.6  1 116.6497        < 2.2e-16 ***
partner_stijl  963.7  2   6.1529         0.003373 ** 
gedeelde_tijd 4608.7  1  58.8528 0.00000000005403 ***
Residuals     5794.8 74                              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

etaSquared(m_otter, type = 3)

                  eta.sq eta.sq.part
partner_stijl 0.07985615   0.1425845
gedeelde_tijd 0.38191162   0.4429925

# Paarsgewijze vergelijking met Bonferroni-correctie, en de adjusted means.
em_otter <- emmeans_test(data = otter_huishouden, huishoudtevredenheid ~ partner_stijl,
                         covariate = gedeelde_tijd, p.adjust.method = "bonferroni")
em_otter

# A tibble: 3 × 9
  term          .y.                  group1      group2         df statistic       p   p.adj p.adj.signif
* <chr>         <chr>                <chr>       <chr>       <dbl>     <dbl>   <dbl>   <dbl> <chr>       
1 partner_stijl huishoudtevredenheid dagslaper   nachtbraker    74    -3.10  0.00271 0.00814 **          
2 partner_stijl huishoudtevredenheid dagslaper   wisselend      74    -0.230 0.818   1       ns          
3 partner_stijl huishoudtevredenheid nachtbraker wisselend      74     2.86  0.00553 0.0166  *           

get_emmeans(em_otter)

# A tibble: 3 × 7
  partner_stijl emmean    se    df conf.low conf.high method      
  <fct>          <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>       
1 dagslaper       66.5  1.78    74     62.9      70.0 Emmeans test
2 nachtbraker     74.2  1.70    74     70.8      77.6 Emmeans test
3 wisselend       67.1  1.81    74     63.5      70.7 Emmeans test

Dezelfde getallen naast elkaar, met de tabel-hulpfunctie van dit oefenboek:

apa_ancova_table(
  m_otter,
  factor          = "partner_stijl",
  covariate       = "gedeelde_tijd",
  data            = otter_huishouden,
  dependent_label = "huishoudtevredenheid",
  factor_label    = "Partner-stijl",
  covariate_label = "gedeelde_tijd",
  caption         = "Otter — huishoudtevredenheid per partner-stijl, raw versus adjusted"
)

Otter — huishoudtevredenheid per partner-stijl, raw versus adjusted
Partner-stijl	n	Mgedeelde_tijd	Mhuishoudtevredenheid, ruw	Mhuishoudtevredenheid, adjusted
Dagslaper	26	18.45	63.38	66.48
Nachtbraker	27	21.03	74.56	74.16
Wisselend	25	22.80	69.88	67.08
Note. Mhuishoudtevredenheid, adjusted = estimated marginal mean voor huishoudtevredenheid bij gemiddelde gedeelde_tijd (20.73). Pooled-within slope bw = 1.36; $\eta^2_p$ = partiële eta-kwadraat. gedeelde_tijd: F(1, 74) = 58.85, p < .001, $\eta^2_p$ = .44; Partner-stijl: F(2, 74) = 6.15, p = .003, $\eta^2_p$ = .14.

a) De parallellisme-aanname is houdbaar: de interactie partner_stijl $\times$ gedeelde_tijd is niet significant, $F(2, 72) = 0.38$, $p = .69$.

b) Beide effecten zijn significant na correctie: factor partner_stijl $F(2, 74) = 6.15$, $p = .003$; covariaat gedeelde_tijd $F(1, 74) = 58.85$, $p < .001$.

c) Gedeelde_tijd heeft de grootste: $\eta^2_p = .44$, dus 44% van de variantie die na partner_stijl overblijft. Partner_stijl: $\eta^2_p = .14$.

d) De paarsgewijze vergelijking met Bonferroni-correctie laat zien dat de nachtbraker significant hoger scoort dan zowel de dagslaper ($p = .008$) als wisselend ($p = .017$); dagslaper en wisselend verschillen onderling niet significant ($p > .99$).

Een ANCOVA met partner_stijl als factor en gedeelde_tijd als covariaat toonde een significant effect van gedeelde tijd op huishoudtevredenheid, $F(1, 74) = 58.85$, $p < .001$, $\eta^2_p = .44$, en een significant effect van partner_stijl na correctie, $F(2, 74) = 6.15$, $p = .003$, $\eta^2_p = .14$. De parallellisme-aanname was houdbaar, $F(2, 72) = 0.38$, $p = .69$. Adjusted means: dagslaper $\bar{Y}^* = 66.48$, nachtbraker $\bar{Y}^* = 74.16$, wisselend $\bar{Y}^* = 67.08$. Paarsgewijze vergelijkingen met Bonferroni-correctie lieten zien dat nachtbrakers significant hoger scoorden dan dagslapers ($p = .008$) en dan wisselende partners ($p = .017$).

Klaar? Zet de codering terug

Je hebt dit hele thema gerekend met options(contrasts = c("contr.sum", "contr.poly")). Die instelling blijft staan zolang je R-sessie openstaat — ook als je straks aan een ander thema begint. Zet hem na een ANOVA of ANCOVA daarom terug op R’s standaard, want in een regressie wil je weer de gewone dummycodering:

# Terug naar R's standaardcodering (dummy-codering).
options(contrasts = c("contr.treatment", "contr.poly"))

Vergeet je dit, dan rekent thema 4 (logistische regressie) in dezelfde sessie met de verkeerde codering, en krijg je bij een categorische voorspeller andere coëfficiënten dan in de uitwerkingen.

3.4 Techniekkeuze — van bevers naar mensen

Slot van dit thema. Bevers en otters zijn een anker — geen einddoel. Wat je hier hebt geleerd, moet ook werken als de proefpersonen mensen zijn.

Een ANCOVA op een beverkolonie is niet anders dan een ANCOVA op een psychotherapie-trial. De techniek is identiek; alleen de context wisselt. Wie alleen oefent op dieren wordt goed in dier-statistiek — niet in statistiek-die-toevallig-bij-dieren-werd-uitgelegd. Daarom dit slot: zeven vignetten, waarvan twee uit het dierenbos en vijf uit een spreekkamer, klaslokaal, bedrijf of ziekenhuis. De vraag is steeds dezelfde: welke techniek past hier?

Loop bij elk vignet je interne stappenplan af:

1. Wie of wat wordt er gemeten?
2. Wat wordt er gemeten? Welke variabelen doen mee?
3. Wat is de onafhankelijke variabele (IV), wat is de afhankelijke (DV)?
4. Meetniveaus per variabele.
5. Dan pas kiezen.

De keuzes komen uit deze cumulatieve set:

• $t$-toets (twee groepen, één continue DV).
• Eenweg-ANOVA ($\geq 3$ groepen, één continue DV, geen covariaat).
• Factorial ANOVA (twee of meer factoren, één continue DV, eventueel interactie).
• MRA (geen factor; continue DV met meerdere continue voorspellers).
• ANCOVA (factor + continue covariaat, één continue DV).

OpmerkingVraag 3.4 — Welke techniek past hier?

Lees onderstaande zeven vignetten en kies per stuk de meest passende techniek.

Vignet 1 — De spreekkamer. Een klinisch psycholoog onderzoekt het effect van drie soorten therapie (CGT, schema-therapie, controle-wachtlijst) op depressie-score (BDI-II) na zes weken. Omdat patiënten verschillend zwaar binnenkomen, wordt ook de baseline-BDI gemeten ($N = 120$, $n \approx 40$ per groep). De onderzoeksvraag: zit er verschil tussen de drie condities, na correctie voor verschil in startniveau?

Vignet 2 — De studieadviseur. Bij $86$ eerstejaars-studenten worden tentamencijfer (op een $10$-schaal), aantal studie-uren in de week vóór het tentamen, en gemiddelde slaapduur per nacht in diezelfde week vastgelegd. De onderzoeksvraag: kunnen studie-uren en slaapduur samen het tentamencijfer voorspellen, en wat is hun relatieve gewicht?

Vignet 3 — De otter-onderzoeker. Bij $78$ otter-koppels wordt huishoudtevredenheid gemeten in drie partner-stijlen (dagslaper / nachtbraker / wisselend). Geen covariaten verzameld. De onderzoeksvraag: verschilt huishoudtevredenheid tussen de drie partner-stijlen?

Vignet 4 — De HR-afdeling. Een organisatiepsycholoog volgt $200$ werknemers in een $2 \times 3$-design: wel/geen leiderschapstraining gevolgd, gekruist met drie afdelingen (sales, IT, support). DV: jaarlijkse werktevredenheids-score. De onderzoeksvraag: is er een hoofd-effect van training, een hoofd-effect van afdeling, en/of een interactie?

Vignet 5 — De IC-arts. Tijdens een corona-piek constateert een onderzoeker dat rokende IC-patiënten minder vaak overlijden dan niet-rokers. Maar rokers zijn gemiddeld jonger dan niet-rokers ($M_{\text{rokers}} = 52$, $M_{\text{niet-rokers}} = 71$). DV: overleving-score (continue klinische ernst-schaal, niet binair). De onderzoeksvraag: blijft het ‘roken-voordeel’ staan na correctie voor leeftijd?

Vignet 6 — Het bever-bos. De bever-onderzoeker (je kent hem inmiddels) vergelijkt damtevredenheid tussen drie werkstijlen (perfectionist / pragmaticus / impulsief), terwijl hij meedoende-familie meeneemt als achtergrondvariabele. De onderzoeksvraag: verschillen werkstijlen in damtevredenheid, ook nadat we voor familie-hulp corrigeren?

Vignet 7 — De schoolpsycholoog. Een onderwijsonderzoeker vergelijkt drie didactische methoden (frontale les / werkvormen / online) op leerprestatie (score $0$–$100$). Hij wil “rekening houden met” de sociaal-economische status (SES) van de leerlingen, gemeten als drie categorieën: laag / midden / hoog. Steekproef $N = 180$, ongeveer $20$ leerlingen per cel. De onderzoeksvraag: verschilt leerprestatie tussen de drie methodes, gegeven verschillen in SES?

Kies per vignet: t-toets / eenweg-ANOVA / factorial ANOVA / MRA / ANCOVA.

Bonus-vraag. Bij vignet 5 — als de onderzoeker ook co-morbiditeit (continue ernst-score) had meegenomen naast leeftijd, welke techniek wordt het dan?

OpgeletAntwoord 3.4 — open na je eigen poging

Vignet	Context	Techniek	Waarom
1	klinisch	ANCOVA	factor (3 condities) + continue covariaat (baseline) + continue DV
2	onderwijs	MRA	geen factor; continue DV met twee continue voorspellers
3	dier	Eenweg-ANOVA	factor (3 stijlen) + continue DV, geen covariaat
4	organisatie	Factorial ANOVA ($2 \times 3$)	twee factoren, één continue DV, vraag naar hoofd-effecten + interactie
5	medisch	ANCOVA	factor (rookstatus) + continue covariaat (leeftijd) + continue DV. Klassieke confounding — typevoorbeeld waar bias-correctie het verhaal kantelt
6	dier	ANCOVA	factor (3 werkstijlen) + continue covariaat (familie-hulp) + continue DV
7	onderwijs	Factorial ANOVA ($3 \times 3$)	Instinker. SES is hier categoriaal (laag/midden/hoog), niet continu. “Rekening houden met” suggereert ANCOVA — maar ANCOVA vereist een continue covariaat. Met een categoriale tweede variabele zit je in factorial-ANOVA-land: methode $\times$ SES, met hoofd-effecten + eventueel interactie. Was SES als IQ-score (continu) gemeten, dan was het wel ANCOVA

Bonus. Twee covariaten ($K = 2$): nog steeds ANCOVA-structuur, alleen met meerdere covariaten — in dit oefenboek aangestipt in T9. Geen aparte techniek, wel een paar extra punten: multicollineariteit tussen covariaten (VIF), $df_w = N - I - K$, en de adjusted means worden geconditioneerd op de grand means van alle covariaten tegelijk.

Inzicht — wat de zeven vignetten samen je geven. Zeven verschillende werelden, vier verschillende technieken. De context stuurt niet de keuze; de structuur doet dat. Een $2 \times 3$-design in een HR-afdeling vraagt dezelfde factorial ANOVA als een $2 \times 3$-design bij dassen-burchten (thema 2). Een baseline-correctie in een depressie-trial vraagt dezelfde ANCOVA als familie-hulp-correctie bij bevers. Wie de structuur leest, kiest goed — onafhankelijk van of de proefpersoon vier poten of twee heeft.

De val: aannemen dat “mensen” een eigen statistiek hebben. Dat is verwarrend; daarom oefenen we met dieren — schone vraag, geen ruis van persoonlijke voorkeuren. Maar de transfer-stap (vignet 1, 2, 4, 5 en 7 hier) moet je net zo soepel kunnen maken. Dat is wat dit slot je vraagt.

Wat blijft liggen

Onderwerpen die buiten dit thema blijven

ANCOVA is een flexibele techniek met veel variaties. Wat dit thema je gaf is de hoofdingang — één factor, één covariaat, additief model, parallellisme-check. Wat hierachter ligt en in dit oefenboek niet behandeld wordt:

• Meerdere covariaten in R uitwerken ($K \geq 2$) — wat er aan de vrijheidsgraden verandert, zie je in T9; hoe het model in R eruitziet, staat in het kader Voor wie verder wil — een tweede covariaat in R. Verder dan dat gaat dit oefenboek niet.
• Polynomiale of niet-lineaire covariaat — lm(Y ~ poly(covariaat, 2) + factor) voor kwadratische trends, of GAM (mgcv::gam) voor flexibele smoothers. Niet behandeld; centreren wordt dan extra belangrijk om VIF onder controle te houden.
• MANCOVA — meerdere afhankelijke variabelen tegelijk, met covariaat-correctie. Thema 5 (MANOVA) behandelt hem niet; daar staat hij, net als hier, onder Wat blijft liggen.
• Causale alternatieven voor ANCOVA bij groepen die niet door loting zijn gevormd — DAG’s, propensity score matching, instrumentele variabelen. Waarom je ze nodig kunt hebben, staat bij T8 (Lord’s paradox); hoe ze werken, valt buiten dit oefenboek.
• Heterogene slopes als ontwerpkeuze — geen “schending” maar een moderation-analyse: het effect van de factor is dan niet één getal maar hangt af van de covariaat-waarde. Te schatten via emtrends() of via simple-slope-uitwerking bij specifieke covariaat-niveaus.
• Robuuste ANCOVA — bij sterke meetfout in de covariaat of niet-normale residuen; zie het WRS2-pakket (Wilcox).
• Bayesiaanse ANCOVA — een alternatief inferentie-frame; BayesFactor of brms als toegankelijke ingangen.

Voor verdieping: Tabachnick & Fidell (2007); Maxwell, Delaney & Kelley (2017, design and analysis); Field (2018, hoofdstuk ANCOVA); Wilcox (2017) voor robuuste varianten.

Aan het eind van de dag

Toen het laat werd, klapte de bever zijn schrift dicht.

“De perfectionisten zijn echt minder tevreden,” zei hij, “ook na correctie voor familie-hulp. Niet omdat ze toevallig minder hulp kregen — wel omdat de werkstijl zelf knaagt.”

“En de impulsievelingen?” vroeg de otter, die langs zwom om iets in zijn eigen schriftje op te zoeken.

“Niet betrouwbaar anders dan de pragmatici,” zei de bever. “Allebei boven het perfectionisme. En meer familie meebouwen helpt sowieso — twee punten tevredener per familielid extra, ongeveer.”

De otter knikte traag. “Bij ons koppels was het de partnerstijl. De nachtbrakers waren het tevredenst, ook na correctie voor de tijd die de partners samen waren. Met de wisselenden? Die zijn moeilijker te lezen.”

Ze zaten bij het water. Iemand had ergens een tak laten vallen. Hij dreef voorbij.

Verantwoording

Dit oefenboek is geschreven voor studenten die covariantieanalyse 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: October 2026 — CountCamp Lab