6. Repeated Measures ANOVA (RMA)

Het verhaal van dit thema

De kameleon zat op een tak en wisselde nauwelijks merkbaar van kleur.

“Wie mij op één moment vastlegt, mist het voornaamste,” zei hij. “Op een ander moment ben ik weer iemand anders, en allebei ben ik.”

Negentig kameleons, vier momenten van de dag — ochtend, middag, avond, nacht. Per kameleon hetzelfde meetinstrument: mate_van_zichzelf ($0$-$100$). Niet vier groepen kameleons; dezelfde kameleon vier keer. Vraag: verschilt de mate-van-zichzelf systematisch tussen momenten — niet tussen dieren, maar binnen elk dier?

De techniek heet repeated measures ANOVA RMA.

Bij RMA stel je dezelfde drie vragen als bij ANOVA, maar nu over herhaalde metingen op dezelfde proefpersoon:

Voorspelt het samen iets? — de multivariate $F$-toets (Pillai’s trace en variants) op het binnen-persoon-effect van moment.
Tussen welke momenten zit het verschil? — polynomiale of geplande contrasten, met Bonferroni-correctie voor paarsgewijze vergelijkingen.
Klopt het allemaal wel? — multivariate normaliteit van de within-subject-vector, geen extreme outliers in trajecten, geen missing data per kameleon.

A. Wat is RMA?

Repeated measures ANOVA test of metingen binnen dezelfde persoon systematisch verschillen tussen condities of momenten. Bij gewone factorial ANOVA zaten in elke groep verschillende dieren; bij RMA zit elke kameleon zelf in elke conditie — ze zijn elkaars vergelijkingsbasis.

Voordeel: tussen-personen-variantie wordt uit de fout gehaald. Wat overblijft is binnen-persoon-variatie — typisch veel kleiner dan tussen-personen, en daarmee een krachtiger toets bij dezelfde $N$.

B. Multivariate vs univariate — onze keuze

RMA kent twee toetsroutes:

Univariate: behandelt de herhaalde metingen als één lange vector, vereist sphericity (gelijke varianties van verschilscores tussen alle paren). Bij schending: Greenhouse-Geisser- of Huynh-Feldt-correctie.
Multivariate: behandelt de herhaalde metingen als $p$ aparte DV’s (zoals MANOVA), met Pillai/Wilks/Hotelling-Lawley/Roy. Geen sphericity-aanname nodig.

	Multivariate	Univariate
Sphericity nodig?	nee	ja (anders GG/HF-correctie)
Power bij kleine $n$	minder krachtig	krachtiger als sphericity klopt
Power bij grotere $n$	krachtig genoeg	krachtig
Aansluiting op MANOVA	direct ($V$, $\Lambda$, $T^2$, $\theta$)	indirect
Didactisch	bouwt voort op thema 5	nieuwe denkstap

Wij volgen de cursus: de multivariate route is de hoofdroute. De Heus verantwoordt de keuze in het cursusboek nuchter:

Dit is niet omdat de multivariate aanpak superieur is aan de univariate aanpak — beide leiden meestal (maar niet altijd!) tot ongeveer dezelfde resultaten, en allebei hebben ze sterke en zwakke kanten. De multivariate aanpak is hier vooral om praktische redenen gekozen: omdat het op deze manier het snelst en eenvoudigst uit te leggen is, voortbouwend op de eerder behandelde MANOVA.

— De Heus, Exercise Book, p. 104 (vertaald uit het Engels)

Pillai is de standaard-keuze, conform thema 5.

Opmerking voor de oude garde — vroeger zat sphericity-correctie (Greenhouse-Geisser, Huynh-Feldt) wél in het curriculum. Met de multivariate route is die correctie niet langer nodig: de multivariate $F$ heeft de aanname helemaal niet. Wie oude collegehandouts leest of software gebruikt die alleen de univariate route biedt, weet nu waarom dat hier ontbreekt.

C. Kameleon-frame — per moment iemand anders

Een kameleon laat zich slecht beschrijven door één meting. Op één moment is hij rustig en bedacht, op een ander moment fel en aanwezig — en in beide gevallen is hij zichzelf, maar dan wel een ander stuk van zichzelf. Wie hem op één moment vastlegt, mist de kleur eronder.

RMA pakt deze veelheid serieus: voor elke kameleon meten we op meerdere momenten, en kijken niet of kameleons onderling verschillen, maar of de momenten binnen kameleons verschillen. Niet “wat is een kameleon?”, maar “hoe verandert hij?”

Eén ding om mee te nemen: omdat dezelfde kameleon vier keer wordt gemeten, kan de volgorde van metingen het resultaat beïnvloeden — het carry-over-effect. Een kameleon die net opgeschrikt is in de ochtend zit nog in zijn schrik in de middag. In de praktijk wordt dat ondervangen door de volgorde te randomiseren of te counterbalancen (bij meerdere experimentele condities elke volgorde even vaak toepassen). Bij vaste tijdmomenten zoals onze ochtend-middag-avond-nacht is randomiseren niet mogelijk; dan analyseer je het in de tijd zoals hij is, en houd je in je interpretatie rekening met dat eerdere metingen latere kunnen kleuren.

Hoe dit hoofdstuk leest

Dit hoofdstuk volgt dezelfde drie-laagse structuur als thema 1, 2, 3, 4 en 5:

Verhaal-frame — wat de kameleon doet of denkt.
Algemene vorm — abstract, statistiek-Latijn met Y1, Y2, Y3, Y4, mijn_data.
Voor onze dieren — uitvoerend, met # hekjes en concrete namen.

Plus de antwoorden in zacht groen (collapse), vragen in indigo, leesregels in mauve, alarm in oranje, don’ts in rood. T-vragen vlechten zich tussen de hoofdvragen door om je rekenintuïtie te oefenen.

Werkmaterialen — R-pakketten en functies

De kameleon pakt zijn gereedschap

Gereedschap	Waar het voor is	Pakket
`lm(cbind(Y1, Y2, Y3, Y4) ~ 1)`	Multivariaat lineair model — basis voor multivariate RMA	base
`car::Anova(lm_obj, idata, idesign, type = 3)`	RMA via multivariate route — Pillai/Wilks/Hotelling/Roy	`car`
`summary(res, multivariate = TRUE)`	Multivariate output met alle vier teststatistics	`car`
`contr.poly(k)`	Polynomiale contrasten (lineair, kwadratisch, kubisch, …)	base
`t.test(y1, y2, paired = TRUE)`	Paarsgewijze vergelijking tussen twee momenten	base
`p.adjust(p, method = "bonferroni")`	Multiple-testing-correctie op $p$-waardes	base
`reshape()`	Wide ↔︎ long format omzetten voor visualisatie	base
`ggplot2::geom_line(aes(group = id))`	Spaghetti-plot van individuele trajecten	`ggplot2`

Waarom lm(cbind(...) ~ 1) + Anova() voor RMA?

R kent geen aparte rmanova()-functie; de gangbare route is dezelfde als bij MANOVA in thema 5: cbind() plakt de momenten tot een matrix-uitkomst, en car::Anova() met idata (de within-subject design-tabel) en idesign (de within-subject formule) doet de RMA. Het intercept-only-model ~ 1 zegt: er zijn geen tussen-persoon-voorspellers — alleen het herhaalde-metingen-deel telt. Bij mixed-design (6.A) komt er wel een between-subject factor in de formule.

Notatie — symbolen voor RMA

Sleutelsymbolen in dit hoofdstuk

In dit werkboek zie je telkens dezelfde notatie:

$n$ — aantal proefpersonen (hier: kameleons). $N$ in absolute zin is hetzelfde, maar bij RMA spreken we vaak over $n$ voor de proefpersonen en $k$ voor het aantal momenten.
$k$ — aantal niveaus van de within-subject factor. Voor de kameleon: $k = 4$ (ochtend, middag, avond, nacht).
$\boldsymbol{\mu}$ — vector van populatie-gemiddelden over de $k$ momenten: $\boldsymbol{\mu} = (\mu_1, \mu_2, \ldots, \mu_k)^\top$.
$\Lambda$ — Wilks’ Lambda binnen de RMA. Klein = sterk effect.
$V$ — Pillai’s trace.
$T^2$ — Hotelling-Lawley trace.
$\theta$ — Roy’s largest root.
$F$ — APA-italic, hier als $F$-approximatie van elk van de vier teststatistics.
$df_H$ — hypothesis $df$. Voor eenwegs-RMA met $k$ momenten: $df_H = k - 1$.
$df_E$ — error $df$ in de $F$-approximatie. Variant per teststatistic.
$df_W$ — within-subject totaal $df = n(k - 1)$.
$p$ — $p$-waarde van de $F$-approximatie.
Polynomiale contrasten — orthogonale gewichten voor lineaire ($\text{L}$), kwadratische ($\text{Q}$), kubische ($\text{C}$) trends over geordende niveaus.

Voor inferentie spreek je over de populatie-vector $\boldsymbol{\mu}^*$ (asterisk-notatie); voor rapportage over de steekproef-vector $\bar{\boldsymbol{Y}}$. APA-italic geldt voor enkele Latijnse symbolen ($F$, $t$, $p$, $r$, $N$, $df$); multi-letter symbolen rechtop ($\text{RMA}$, $\text{MS}$, $\text{Pillai}$).

Vuistregels zijn afspraken, geen wetten

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

Aantal proefpersonen: $n_{\text{kleinste}} \geq 20$ per groep maakt de multivariate $F$ robuust tegen schending van multivariate normaliteit (zoals bij MANOVA in thema 5). Onder $n < k + 5$ wordt de Pillai-toets onbetrouwbaar (te weinig $df_E$). Check je college-sheets of jouw vakgroep deze drempel hanteert.
Welke teststatistic: bij RMA met intercept-only-model zijn Pillai, Wilks en Hotelling-Lawley wiskundig identiek (zelfde $F$-waarde). Roy idem als er één hypothese-dimensie is. Standaard rapporteer je Pillai als primary, conform thema 5.
Bonferroni voor paren: bij $k$ momenten zijn er $k(k-1)/2$ paarsgewijze vergelijkingen. Voor $k = 4$: zes paren, $\alpha_{\text{per paar}} = .05/6 \approx .008$. Strenger dan unprotected paired $t$-tests.
Polynomiale contrasten: voor geordende niveaus (tijd, dosis) geven contr.poly orthogonale lineaire / kwadratische / kubische trends. Voor ongeordende niveaus (kleuren, condities) horen ze niet thuis.
Effectgrootte multivariaat: Pillai $V$ zelf is een effectgrootte (analoog aan $\eta^2$); $V \geq .14$ is een groot effect (Cohen-conventie).
Sphericity: in dit werkboek niet als hoofdroute behandeld, maar de univariate-route bestaat. Greenhouse-Geisser- en Huynh-Feldt-correcties horen daar bij — niet bij de multivariate route.

Verschillende vakgroepen, docenten en handboeken kiezen iets andere drempels. Als je voor een tentamen of opdracht werkt, check altijd in je eigen college-sheets welke drempels je docent of vakgroep hanteert. Wij volgen hier consistent één set, maar dat is een keuze, geen natuurwet.

Code in dit hoofdstuk — wat moet je kunnen typen?

Code is standaard open als je hem moet kunnen typen op het R-practical-tentamen: lm(cbind(Y1, Y2, Y3, Y4) ~ 1), Anova(., idata, idesign, type = 3), summary(., multivariate = TRUE), contr.poly(), t.test(., paired = TRUE), p.adjust(). Code die alleen ter illustratie dient — verkenning, plot-decoratie, spaghetti-plot, gt-tabel — staat ingeklapt met een knopje “Toon code”. Klap hem open als je nieuwsgierig bent; voor het tentamen hoef je hem niet te reproduceren.

6.0 Project- en datavoorbereiding

Aantekeningen netjes, paden kloppen

Open de meegestuurde projectmap (06_repeated_measures/) en dubbelklik op 06_repeated_measures.Rproj. Daarmee staat je werkomgeving klaar: paden kloppen, RStudio kent de juiste werkmap. De data ligt al in data/, de helper-functies in functions/. Niets te installeren als de pakketten (car, tidyverse) op je computer staan; anders één keer:

install.packages(c("car", "tidyverse"))

Open de hoofd-dataset.

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

'data.frame':   90 obs. of  5 variables:
 $ id     : Factor w/ 90 levels "k01","k02","k03",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ ochtend: num  54.1 77.4 70.1 44.8 56.6 61.5 72.8 55.1 45.3 61.4 ...
 $ middag : num  27.5 52.4 29.9 33.4 25.4 37.8 41.9 16.4 15.3 10.3 ...
 $ avond  : num  51.7 78.4 50.8 34.6 41.7 47.2 35.9 34.8 35.5 57.8 ...
 $ nacht  : num  79 84 100 88.3 70.4 85.1 71.3 68.8 51.9 73 ...

Per kameleon vijf waarden: een id plus de mate-van-zichzelf op vier momenten — ochtend, middag, avond, nacht (elk continu, schaal $0$-$100$). De vier momenten vormen samen het binnen-persoon-traject.

Wide format versus long format

De data staat in wide format: één rij per kameleon, vier kolommen voor de momenten. Voor RMA via lm(cbind(...) ~ 1) gebruiken we precies dit format — cbind() plakt de momenten tot een matrix-uitkomst.

Voor visualisatie (spaghetti-plot, boxplot per moment) is long format handiger: één rij per meting, met een moment-kolom en een mate-kolom. We zetten dat één keer om in 6.1.a en bewaren beide naast elkaar.

Belangrijke conceptuele knoop: in wide format zien de vier momenten eruit als vier kolommen — alsof het vier verschillende DV’s zijn (zoals bij MANOVA in thema 5). Dat is een illusie. RMA heeft conceptueel één DV (mate-van-zichzelf) die op vier momenten is gemeten. De kolom-vorm is alleen een handige opslag — niet de inhoudelijke structuur. Bij MANOVA had elke leeuw vier aparte gevoelens; hier heeft elke kameleon vier metingen van hetzelfde gevoel. Die structuur-omkering zorgt dat de procedure en de aannames anders zijn: bij RMA komt sphericity erbij, en het multivariate-vs-univariate-pad krijgt een andere betekenis (zie T3 en de Don’t-callout verderop).

6.1 Het RMA-model fitten en interpreteren

Vier momenten op dezelfde dieren — één hypothese

“Vier keer dezelfde meting bij hetzelfde dier,” mompelde de kameleon. “Niet of de ene kameleon meer zichzelf is dan de andere, maar of de momenten binnen mij van elkaar verschillen. Mijn ochtend tegen mijn middag — daar zit de toets.”

Voor je gaat rekenen — vijf vragen aan jezelf

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

Wie of wat wordt er gemeten? De kameleon heeft $90$ kameleons genoteerd, elk vier keer.
Wat wordt er gemeten? Eén variabele per moment: mate_van_zichzelf. Plus de moment-aanduiding zelf (ochtend/middag/avond/nacht).
Onafhankelijk of afhankelijk? Onderzoeksvraag: voorspelt moment de mate-van-zichzelf? Dan is mate_van_zichzelf afhankelijk ($Y$), en moment onafhankelijk (within-subject factor, $k = 4$).
Meetniveau van elke variabele? Mate_van_zichzelf interval ($0$-$100$); moment ordinaal (geordend langs de dag) of nominaal (vier categorieën).
Welke techniek? Doorloop de beslisboom:

flowchart TD
    A[Hetzelfde dier<br/>meerdere keren gemeten?] -->|ja| Z[RMA<br/><i>dit thema</i>]
    A -->|nee| B[Hoeveel afhankelijke<br/>variabelen?]
    B -->|één| C[Y meetniveau?]
    B -->|meerdere<br/>continue| D[MANOVA<br/><i>thema 5</i>]
    C -->|interval| E[X-en meetniveau?]
    C -->|binair| F[Logistic regression<br/><i>thema 4</i>]
    E -->|alleen interval| G[Multiple regression<br/><i>thema 1</i>]
    E -->|alleen nominaal| H[ANOVA<br/><i>thema 2</i>]
    E -->|gemengd| I[ANCOVA<br/><i>thema 3</i>]
    style Z fill:#faf3e2,stroke:#c9a05a,stroke-width:2px

Bens pijltjes-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

Eén $Y$ (mate_van_zichzelf), meerdere metingen op dezelfde proefpersoon. Het herhaalde-metingen-design is wat de techniek bepaalt. Antwoord: RMA.

T1 — Techniek-keuze: RMA, eenwegs-ANOVA, paired t-toets of mixed-effects?

Vraag T1 — Pen-en-papier

Voor elk vignet: bepaal het design (binnen-persoon, tussen-persoon, of gemengd) en kies dan de techniek.

a) Een wolf meet bij $50$ wolfjes hun jachtsucces in week 1 en week 8 van een training. Vraag: verandert het succes?

b) Een hert meet bij $80$ herten op vier seizoenen het hormoonniveau (lente, zomer, herfst, winter). Elk hert wordt op alle vier momenten gemeten.

c) Een onderwijskundige meet bij $60$ studenten de leesmotivatie, met als groep digitaal-onderwijs versus klassikaal-onderwijs (between).

d) Een gezinspsycholoog meet bij $45$ koppels relatie-tevredenheid op drie momenten ($t_0$, $t_1$, $t_2$) plus per koppel of ze stedelijk of rurale woonomgeving hebben.

Antwoord T1 — open na je eigen poging

a) Twee herhaalde metingen op dezelfde wolfjes. $\Rightarrow$ paired $t$-toets (twee niveaus is een speciaal geval van RMA — RMA met $k = 2$ is wiskundig equivalent).

b) Vier herhaalde metingen op dezelfde herten. $\Rightarrow$ RMA ($k = 4$, één within-subject factor seizoen). De multivariate route uit dit thema.

c) Twee groepen, één meting per student. $\Rightarrow$ eenwegs-ANOVA (thema 2). Geen herhaalde metingen.

d) Drie momenten binnen elk koppel (within) plus woonomgeving stedelijk/ruraal (between). $\Rightarrow$ mixed-design RMA (between $\times$ within), zie 6.A. Bij geneste of gekruiste random effects: mixed-effects model met lme4::lmer.

Vuistregel. Hetzelfde dier meerdere keren = RMA-familie. Twee metingen → paired $t$. Meer dan twee → RMA. Plus een tussen-persoon factor → mixed-design RMA. Random effects op meerdere niveaus → mixed-effects.

T2 — Waarom RMA krachtiger dan one-way ANOVA?

Vraag T2 — Pen-en-papier

Een onderzoeker overweegt: hij heeft $90$ kameleons elk vier keer gemeten. Stel hij gooit het design “weg” en behandelt het als een eenwegs-ANOVA met $360$ onafhankelijke observaties (vier groepen van $90$).

a) Welke informatie gooit hij weg door dat te doen?

b) In de RMA-decompositie wordt tussen-persoon-variantie uit het residu gehaald. Wat blijft er over als residu?

c) Bij dezelfde groepsgemiddelden: wordt $F$ groter of kleiner als de residu-MS daalt?

d) In welke situatie levert RMA geen power-winst op boven eenwegs-ANOVA?

Antwoord T2 — open na je eigen poging

a) De koppeling tussen de vier metingen — dat het dezelfde kameleon is op vier momenten. Een kameleon die op ochtend hoog scoort, heeft kans om ook op nacht hoog te scoren (zijn algemene niveau). Door die koppeling te negeren, behandelt hij vier waarnemingen van dezelfde kameleon als vier onafhankelijke. Dat is statistisch incorrect (dependent observations) en weggegooide informatie.

b) RMA splitst de totale $SS$ in: tussen-personen ($SS_{BS} = $ koppen die hoog of laag scoren over alle momenten heen), tussen momenten ($SS_{moment}$, het effect dat we willen toetsen), en within-subject residu ($SS_{resW}$, alleen ruis na verwijdering van persoon en moment). Dat residu is veel kleiner dan het ANOVA-residu.

c) Groter. $F = MS_{moment} / MS_{resW}$. Als de noemer daalt en de teller gelijk blijft, schiet $F$ omhoog. Vandaar de power-winst.

d) Als de tussen-persoon-correlaties $\rho \approx 0$ — dan is er geen koppeling, alle metingen lijken al op onafhankelijke trekkingen. RMA biedt dan geen winst (en kost zelfs $df$). Vuistregel: RMA verdient zijn keuze bij gemiddelde within-subject correlaties $> .3$.

Inzicht. RMA is geen “andere” toets — het is dezelfde $F$-logica met een slimme decompositie van het residu. De winst zit in de noemer.

6.1.a Verkenning — gemiddelden, correlaties en spaghetti

Voordat de kameleon fit, kijkt hij eerst.

Algemene vorm.

# Gemiddelden per moment.
colMeans(mijn_data[, c("Y1", "Y2", "Y3", "Y4")])

# Correlatiematrix tussen momenten — within-subject samenhang.
cor(mijn_data[, c("Y1", "Y2", "Y3", "Y4")])

# Boxplot per moment.
boxplot(mijn_data[, c("Y1", "Y2", "Y3", "Y4")])

Voor onze dieren.

# Gemiddelden per moment.
round(colMeans(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")]), 2)

ochtend  middag   avond   nacht 
  60.09   33.24   46.55   73.49

# SD's per moment.
round(apply(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")], 2, sd), 2)

ochtend  middag   avond   nacht 
  15.37   14.02   13.44   15.15

# Correlatiematrix tussen de vier momenten.
round(cor(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")]), 2)

        ochtend middag avond nacht
ochtend    1.00   0.47  0.50  0.37
middag     0.47   1.00  0.33  0.36
avond      0.50   0.33  1.00  0.40
nacht      0.37   0.36  0.40  1.00

Toon code (boxplot per moment + spaghetti-plot)

par(mfrow = c(1, 2), mar = c(4, 4, 2.5, 1))

# Boxplot per moment.
boxplot(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")],
        col  = c("#0077BB", "#33BBEE", "#EE7733", "#CCBB44"),
        ylab = "Mate van zichzelf",
        main = "Per moment",
        names = c("Ochtend", "Middag", "Avond", "Nacht"))

# Spaghetti-plot — individuele trajecten + groepsgemiddelde.
matY <- as.matrix(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])
matplot(t(matY), type = "l", lty = 1, col = adjustcolor("#888888", alpha.f = 0.35),
        xaxt = "n", xlab = "Moment", ylab = "Mate van zichzelf",
        main = "Trajecten per kameleon")
axis(1, at = 1:4, labels = c("Ochtend", "Middag", "Avond", "Nacht"))
lines(1:4, colMeans(matY), col = "#CC3311", lwd = 3)
points(1:4, colMeans(matY), col = "#CC3311", pch = 16, cex = 1.4)

Toon code (boxplot per moment + spaghetti-plot)

par(mfrow = c(1, 1))

Means-tabel — moment-overzicht met `gt`

Toon code (gt-tabel met means + SD per moment)

means_v <- colMeans(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])
sds_v   <- apply(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")], 2, sd)

tab_df <- data.frame(
  moment = c("Ochtend", "Middag", "Avond", "Nacht"),
  n      = nrow(kameleon_momenten),
  M      = means_v,
  SD     = sds_v
)

tab_df |>
  gt() |>
  tab_header(title = "Kameleon — mate-van-zichzelf per moment") |>
  fmt_number(columns = c(M, SD), decimals = 2) |>
  cols_label(
    moment = "Moment",
    n      = md("*n*"),
    M      = md("*M*"),
    SD     = md("*SD*")
  ) |>
  tab_options(table.font.size = px(12),
              table.border.top.style    = "none",
              table.border.bottom.style = "none")

Moment	n	M	SD
Kameleon — mate-van-zichzelf per moment
Ochtend	90	60.09	15.37
Middag	90	33.24	14.02
Avond	90	46.55	13.44
Nacht	90	73.49	15.15

Vragen 6.1

a) Hoeveel kameleons heeft de kameleon genoteerd? Hoeveel metingen heeft hij in totaal?

b) Op welk moment is de mate-van-zichzelf het hoogst? Het laagst?

c) Past het patroon bij de verwachting (ochtend redelijk → middag dip → avond herstel → nacht piek)?

d) Wat zegt de correlatiematrix over de within-subject samenhang? Zijn de momenten allemaal even sterk gekoppeld?

e) Wat zie je in de spaghetti-plot dat je in de boxplots niet ziet?

Antwoord 6.1 — open na je eigen poging

In gewone woorden.

$n = 90$ kameleons, vier momenten elk = $360$ metingen totaal.
De mate-van-zichzelf is het hoogst in de nacht ($M \approx 73.49$) en het laagst in de middag ($M \approx 33.24$).
Het patroon klopt: ochtend redelijk ($M \approx 60$), middag dip ($M \approx 33$), avond herstel ($M \approx 47$), nacht piek ($M \approx 73$) — een opvallende U-vorm met een dip in de middag, en daarbovenop een algemene stijging over de dag.
De correlatiematrix laat $r$’s tussen $.33$ en $.50$ zien — alle momenten matig positief gekoppeld, maar geen extreme samenhang. De nacht-correlaties met de overige momenten zijn iets lager dan de ochtend-correlaties, wat suggereert dat het nacht-zelf wat losser staat van de rest van de dag.
De spaghetti-plot toont individuele trajecten: niet elke kameleon volgt netjes het gemiddelde patroon. Sommigen zijn in de middag bijna tot $0$ gezakt; anderen blijven hoog. De rode lijn is het gemiddelde-traject; daaronder zit de variatie tussen kameleons.

APA-stijl.

Van $n = 90$ kameleons werden vier metingen verzameld op één dag ($N = 360$ metingen). De gemiddelde mate-van-zichzelf was het hoogst in de nacht ($M = 73.49$, $SD = 15.15$) en het laagst in de middag ($M = 33.24$, $SD = 14.02$). Within-subject correlaties tussen de momenten waren matig positief (range $r = .33$-$.50$).

T3 — Multivariate vs univariate route, en de sphericity-aanname

Vraag T3 — Pen-en-papier

a) Wat houdt sphericity precies in voor RMA?

b) Waarom is multivariate RMA niet afhankelijk van sphericity?

c) In welke situatie kan de univariate route krachtiger zijn dan de multivariate?

d) Wat doet de Greenhouse-Geisser-correctie als sphericity geschonden is?

Antwoord T3 — open na je eigen poging

a) Sphericity = de varianties van alle paarsgewijze verschilscores tussen momenten zijn gelijk. Bij $k = 4$ momenten zijn er $\binom{4}{2} = 6$ verschilparen; sphericity zegt dat $\text{Var}(Y_i - Y_j)$ voor alle $i, j$ ongeveer hetzelfde is. Dit is een sterkere aanname dan “gelijke varianties per moment”. Mauchly’s $W$-toets test dit, met $H_0$: sphericity geldt.

b) Multivariate RMA behandelt de momenten als $k$ aparte DV’s (zoals MANOVA): de toets gebruikt de volledige covariantie-matrix tussen momenten in plaats van te aannemen dat verschilscore-varianties gelijk zijn. Geen aanname over de structuur van die matrix dus geen sphericity-eis.

c) Bij klein $n$ ($n < k + 5$): de multivariate route heeft te weinig $df_E$ en wordt onbetrouwbaar. Univariate route met GG- of HF-correctie behoudt dan power. Bij grotere $n$ wint multivariate vaak.

d) Greenhouse-Geisser past de $df$ van de univariate $F$-toets aan: $df_H \cdot \hat\epsilon$ en $df_E \cdot \hat\epsilon$, waar $\hat\epsilon$ tussen $1/(k-1)$ en $1$ ligt. Bij sterke sphericity-schending is $\hat\epsilon$ klein en worden $df$’s kleiner — strengere $F$-drempel, behoud van type-I-fout. Huynh-Feldt is een ietsje minder conservatieve variant van dezelfde correctie.

Inzicht. In dit werkboek volgen we de multivariate route als hoofdkeuze: bij $n = 90$ ruim genoeg, en geen sphericity-discussie nodig. Voor wie de univariate route wil: zie Wat blijft liggen en je eigen college-aantekeningen.

Mauchly’s $W$, Greenhouse-Geisser, Huynh-Feldt — wat lees je in de output?

In de RMA-output van Leiden’s Anova() zie je vier dingen die met sphericity te maken hebben:

Mauchly’s $W$ — toets op sphericity ($H_0$: sphericity geldt). Significant $\rightarrow$ sphericity geschonden.
Greenhouse-Geisser (GG) $\hat\epsilon$ — correctie-factor tussen $1/(k-1)$ en $1$. Hoe verder van $1$, hoe sterker de schending.
Huynh-Feldt (HF) $\hat\epsilon$ — iets minder strenge variant van GG.
Bijbehorende GG- en HF-$p$-waardes — de univariate $F$-toets met aangepaste $df$.

Vuistregel voor het lezen: kijk eerst naar Mauchly. Bij $p > .05$ kun je in principe sphericity aannemen en de unaangepaste univariate $F$ lezen. Bij $p < .05$ rapporteer je liever de GG- of HF-correctie. In ons werkboek volgen we de multivariate route als hoofdkeuze (geen sphericity-aanname nodig), maar je moet de Mauchly-/GG-/HF-output wel kunnen lezen in de Leiden-uitvoer.

6.1.b RMA fitten — Leiden’s route via `factorial_design()`

Leiden’s R-practical schrijft een wrapper voor: factorial_design() uit rstatix. Die definieert het RMA-design (welke variabelen, welke factor is within, welke is between), en geeft je een object dat je in Anova() van car stopt voor de feitelijke toets.

Algemene vorm (long-format data).

# Stap 0: data in long format brengen — één rij per persoon × moment.
library(tidyr)
long_data <- mijn_wide_data |>
  pivot_longer(cols      = c(Y1, Y2, Y3, Y4),
               names_to  = "moment",
               values_to = "score")

# Stap 1: factorial_design() — vertelt rstatix wat WS-factor is.
library(rstatix)
fd <- factorial_design(
  data   = long_data,
  dv     = score,
  wid    = id,        # persoons-identificatie (within-subject id)
  within = moment     # WS-factor
)

# Stap 2: Anova() op het lm-object dat factorial_design() teruggeeft.
# icontrasts schrijf je hier voor, NIET via options(contrasts=...).
library(car)
res_rma <- Anova(
  fd$model,
  idata      = fd$idata,
  idesign    = fd$idesign,
  icontrasts = c("contr.sum", "contr.poly"),
  type       = "III"
)

# Stap 3: samenvatting — multivariate + univariate samen.
summary(res_rma, multivariate = TRUE)

factorial_design() doet de wide↔︎long-omzetting onder de motorkap, bouwt de juiste idata en idesign, en retourneert een lijst met $model, $idata, $idesign. Je geeft die door aan Anova() met icontrasts als argument — niet via options(contrasts=). Conform Leiden’s R-practical hoofdstuk 6.

Voor onze dieren.

library(tidyr)
library(rstatix)
library(car)

# Stap 0: kameleon_momenten staat in wide format — eerst long maken.
kameleon_long <- kameleon_momenten |>
  pivot_longer(
    cols      = c(ochtend, middag, avond, nacht),
    names_to  = "moment",
    values_to = "mate_van_zichzelf"
  ) |>
  mutate(moment = factor(moment,
                         levels = c("ochtend", "middag", "avond", "nacht")))

# Stap 1: factorial_design — kameleon als wid, moment als within-factor.
fd_kam <- factorial_design(
  data   = kameleon_long,
  dv     = mate_van_zichzelf,
  wid    = id,
  within = moment
)

# Stap 2: Anova met icontrasts en type "III".
res_kam <- Anova(
  fd_kam$model,
  idata      = fd_kam$idata,
  idesign    = fd_kam$idesign,
  icontrasts = c("contr.sum", "contr.poly"),
  type       = "III"
)

# Stap 3: samenvatting — multivariate en univariate gelijk weergegeven.
summary(res_kam, multivariate = TRUE)


Type III Repeated Measures MANOVA Tests:

------------------------------------------
 
Term: (Intercept) 

 Response transformation matrix:
        (Intercept)
avond             1
middag            1
nacht             1
ochtend           1

Sum of squares and products for the hypothesis:
            (Intercept)
(Intercept)     4097579

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1  0.961061 2196.619      1     89 < 2.22e-16 ***
Wilks             1  0.038939 2196.619      1     89 < 2.22e-16 ***
Hotelling-Lawley  1 24.681113 2196.619      1     89 < 2.22e-16 ***
Roy               1 24.681113 2196.619      1     89 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: moment 

 Response transformation matrix:
        moment1 moment2 moment3
avond         1       0       0
middag        0       1       0
nacht         0       0       1
ochtend      -1      -1      -1

Sum of squares and products for the hypothesis:
          moment1   moment2   moment3
moment1  16494.43  32709.98 -16319.79
moment2  32709.98  64866.92 -32363.66
moment3 -16319.79 -32363.66  16147.00

Multivariate Tests: moment
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1  0.883892 220.7672      3     87 < 2.22e-16 ***
Wilks             1  0.116108 220.7672      3     87 < 2.22e-16 ***
Hotelling-Lawley  1  7.612663 220.7672      3     87 < 2.22e-16 ***
Roy               1  7.612663 220.7672      3     87 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

             Sum Sq num Df Error SS den Df F value    Pr(>F)    
(Intercept) 1024395      1    41505     89  2196.6 < 2.2e-16 ***
moment        81118      3    33532    267   215.3 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Mauchly Tests for Sphericity

       Test statistic p-value
moment        0.95598 0.55674


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

        GG eps Pr(>F[GG])    
moment 0.97178  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

         HF eps   Pr(>F[HF])
moment 1.008238 5.813392e-71

Klassieke car-route — alternatief

Voor wie geen rstatix::factorial_design() wil gebruiken: dezelfde toets kan via een direct lm(cbind(...) ~ 1)-object plus Anova() met losse idata/idesign:

options(contrasts = c("contr.sum", "contr.poly"))
lm_kam <- lm(cbind(ochtend, middag, avond, nacht) ~ 1, data = kameleon_momenten)
idata  <- data.frame(moment = factor(c("ochtend","middag","avond","nacht"),
                                     levels = c("ochtend","middag","avond","nacht")))
res    <- Anova(lm_kam, idata = idata, idesign = ~moment, type = 3)
summary(res, multivariate = TRUE)

Wiskundig gelijk; alleen is dit niet wat Leiden’s R-practical voorschrijft. Voor het tentamen ken je de factorial_design()-route.

Vrijheidsgraden bij RMA — hoofdrekenbaar voorbeeld

Voor de eenwegs-RMA met factor moment ($k = 4$) en $n = 90$ kameleons:

Bron	Formule	Hier
Within-subjects (binnen-persoon) totaal	$df_W = (n-1)(k-1)$	$89 \cdot 3 = 267$
Conditie (moment)	$df_{cond} = k - 1$	$3$
Residu within	$df_{resW} = (n-1)(k-1)$	$267$
Multivariate $df_H$	$k - 1$	$3$
Multivariate $df_E$	$n - k + 1$	$87$

Hoofdreken-controle. Bij $n = 90$, $k = 4$: univariate $F(3, 267) = 215.30$ (verschijnt in het stuk over alternatieve route); multivariate $F(3, 87)$ in de Pillai-output.

Bij $n = 20$, $k = 3$ (kleiner voorbeeld): $df_{cond} = 2$, $df_{resW} = 38$, multivariate $F(2, 18)$.

6.1.c De vier teststatistics — Pillai als primary

Vragen 6.1

f) Voor de term moment: welke vier teststatistics rapporteert de output? Welke waardes hebben ze?

g) Welke $df$ en $p$-waardes horen erbij? Wijzen ze allemaal in dezelfde richting?

h) Welke teststatistic rapporteer je als primary? Waarom?

i) Waarom zijn Pillai, Wilks en Hotelling-Lawley voor de term moment allemaal hetzelfde $F$-getal (en Roy ook)?

Antwoord 6.1 — open na je eigen poging

In gewone woorden.

Voor moment: Pillai’s $V = 0.88$, Wilks’ $\Lambda = 0.12$, Hotelling-Lawley $T^2 = 7.61$, Roy’s $\theta = 7.61$.
$df_H = 3$, $df_E = 87$ in de $F$-approximatie; alle vier $F = 220.77$, $p < .001$ — allemaal hetzelfde antwoord.
Standaard rapporteer je Pillai’s trace als primary: het is robuust tegen lichte schendingen van multivariate normaliteit en covariantie-matrix-gelijkheid.
De vier teststatistics geven hetzelfde $F$-getal omdat dit een intercept-only-RMA is — er is maar één hypothese-dimensie (de within-subject factor moment); alle vier reduceren tot dezelfde wiskundige formulering. Bij mixed-design (between $\times$ within) of bij MANOVA met meerdere factoren lopen ze wél uiteen, zoals in thema 5.

APA-stijl.

Een eenwegs-RMA met de multivariate route toonde een sterk effect van moment op de mate-van-zichzelf: Pillai’s $V = .88$, $F(3, 87) = 220.77$, $p < .001$. De andere drie teststatistics gaven hetzelfde $F$-getal (Wilks’ $\Lambda = .12$, Hotelling-Lawley $T^2 = 7.61$, Roy’s $\theta = 7.61$), zoals te verwachten bij een intercept-only-RMA.

T4 — Polynomiale contrasten in RMA

Vraag T4 — Pen-en-papier

Bij geordende within-subject niveaus (tijd, dosis, moment) kun je in plaats van paarsgewijze vergelijkingen ook polynomiale contrasten toetsen.

a) Wat toetst een lineair contrast bij vier momenten?

b) Wat toetst een kwadratisch contrast?

c) Wat toetst een kubisch contrast?

d) Polynomiale contrasten zijn orthogonaal: wat betekent dat hier?

e) Bij vier momenten zijn er drie contrasten ($k - 1$). Waarom precies drie?

Antwoord T4 — open na je eigen poging

a) Lineair contrast (gewichten $-3, -1, 1, 3$ voor vier geordende niveaus): test of de mate-van-zichzelf een monotone stijging of daling over de momenten vertoont. Significant lineair = de mate-van-zichzelf neemt systematisch toe (of af) over de dag.

b) Kwadratisch contrast (gewichten $1, -1, -1, 1$): test of de mate-van-zichzelf een U-vorm of omgekeerde U-vorm over de momenten heeft — dip in het midden of piek in het midden. Significant kwadratisch bij onze kameleon = “mate-van-zichzelf dipt in middag, hoog aan de uiteinden” (of andersom).

c) Kubisch contrast (gewichten $-1, 3, -3, 1$): test of de mate-van-zichzelf een op-neer-op of neer-op-neer-patroon over de momenten heeft — een S-vormige beweging. Vaak kleinste van de drie als de hoofdpatronen in de mate-van-zichzelf al in lineair en kwadratisch zitten.

d) Orthogonaal = de contrasten meten onafhankelijke stukken van het patroon. Het lineaire contrast deelt geen variantie met het kwadratische; samen tellen ze alle variatie tussen momenten op tot $SS_{cond}$. Voordelen: elke trend-toets is interpretabel los van de andere.

e) Bij $k$ niveaus zijn er $k - 1$ orthogonale contrasten — analoog aan $df_{cond} = k - 1$. Vier momenten geven drie onafhankelijke vragen: lineair, kwadratisch, kubisch. Alles wat daarna nog komt zou afhankelijk zijn van eerdere contrasten.

Inzicht. Polynomiale contrasten zijn de voorkeur bij geordende niveaus; bij ongeordende condities (bv. drie ongelijksoortige interventies) horen ze niet thuis — kies dan paarsgewijze of geplande contrasten.

Rangschikking bij meerdere significante contrasten. Meerdere polynomiale contrasten kunnen tegelijk significant zijn — lineair én kwadratisch én kubisch. Statistisch is elk verdedigbaar als losse uitspraak. Voor interpretatie rangschik je op absolute contrast-waarde: het sterkste contrast is de dominante trend, de andere zijn aanvullende patroon-componenten. Rapporteer in volgorde van sterkte, niet in volgorde van rang. Zie de surprise-opgave in 6.4 voor een geval waar de visuele intuïtie een ander contrast aanwijst dan het rekenen.

6.1.d Multivariate effectgrootte

Pillai als effectgrootte voor RMA

Pillai’s $V$ tussen $0$ en $\min(k - 1, df_H)$ is zelf al een effectgrootte: het analoog van $\eta^2$ in de univariate ANOVA. Cohens vuistregels:

klein: $V \approx .01$
matig: $V \approx .06$
groot: $V \geq .14$

Voor onze kameleon: $V = .88$, ver boven de grote-effect-grens — een zeer sterk binnen-persoon-effect.

Een alternatief is de partial $\eta^2$ van het univariate equivalent: $\eta^2_p = SS_{cond} / (SS_{cond} + SS_{resW})$. Voor onze data: $\eta^2_p = 81118 / (81118 + 33532) = .708$.

# Univariate equivalent voor effectgrootte-vergelijking.
# Long format omzetting voor de univariate aov-route.
long_kam <- reshape(
  kameleon_momenten,
  varying   = list(c("ochtend", "middag", "avond", "nacht")),
  v.names   = "mate",
  timevar   = "moment",
  times     = c("ochtend", "middag", "avond", "nacht"),
  direction = "long"
)
long_kam$moment <- factor(long_kam$moment,
                          levels = c("ochtend", "middag", "avond", "nacht"))

# Univariate RMA via aov() — geeft SS_cond en SS_resW direct af.
fit_uni <- aov(mate ~ moment + Error(id/moment), data = long_kam)
summary(fit_uni)


Error: id
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 89  41505   466.4               

Error: id:moment
           Df Sum Sq Mean Sq F value Pr(>F)    
moment      3  81118   27039   215.3 <2e-16 ***
Residuals 267  33532     126                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Partial eta^2 = SS_moment / (SS_moment + SS_residu_within)
ss_moment <- 81118
ss_resW   <- 33532
eta2_p    <- ss_moment / (ss_moment + ss_resW)
round(eta2_p, 3)

[1] 0.708

Vraag 6.1

j) Wat is de univariate $F$-waarde voor moment? Met welke $df$?

k) Welk percentage van de within-subject variantie wordt verklaard door moment (partial $\eta^2$)?

l) Hoe verhoudt $V = .88$ zich tot $\eta^2_p = .71$? Beide zeer groot — waarom verschillen ze?

Antwoord 6.1 — open na je eigen poging

In gewone woorden.

Univariaat: $F(3, 267) = 215.30$, $p < .001$ — sterk significant en met grote $df_E$ door alle binnen-persoon-residuen samen.
Partial $\eta^2 = .708$ — moment verklaart ongeveer $71\%$ van de within-subject-variantie, zeer groot effect.
$V = .88$ en $\eta^2_p = .71$ zijn niet gelijk omdat ze andere maten zijn: $V$ telt eigenvalues van de hypothese-matrix, $\eta^2_p$ is een SS-ratio. Ze convergeren in dezelfde richting (groot effect) maar volgen verschillende schalen. Beide rapporteer je gerust naast elkaar.

APA-stijl.

De multivariate effectgrootte was zeer groot, Pillai’s $V = .88$. De univariate-equivalent partial $\eta^2 = .71$ bevestigt dit beeld.

6.2 Contrasten en follow-up

“Welke momenten zijn de echte motoren?”

“Ik weet nu dat de momenten verschillen,” zei de kameleon. “Maar in welke vorm? Stijgt het over de dag, of dipt het in de middag, of springt het op-en-neer? En welke specifieke paren zijn de drijvende krachten?”

De multivariate $F$ zegt alleen dát moment ertoe doet. Voor welke trends en welke paren zit het effect, heb je follow-up nodig:

Polynomiale contrasten — lineair / kwadratisch / kubisch via contr.poly.
Geplande contrasten — theoretisch gemotiveerde vergelijkingen (bv. nacht versus ochtend, middag versus rest).
Paarsgewijze paired $t$-toetsen met Bonferroni-correctie voor alle $k(k-1)/2$ paren.

Wat is een contrast eigenlijk?

Een contrast is gewoon een ingewikkeld verschil van meer dan twee dingen — met een eentje hier, een halfje daar, een minnetje verderop. Bij twee gemiddelden is een contrast simpel: $\mu_1 - \mu_2$, met gewichten $(+1, -1)$. Voor het toetsen heet dat een gepaarde $t$-toets of een ongepaarde, afhankelijk van design.

Bij drie of meer gemiddelden kun je oneindig veel contrasten verzinnen — sommige zijn theoretisch zinnig (bv. “is de eerste meting anders dan het gemiddelde van de andere drie?”), andere zijn willekeurig. Een gewogen optelling $L = c_1\mu_1 + c_2\mu_2 + c_3\mu_3 + c_4\mu_4$ heet een contrast zodra de gewichten optellen tot $0$: $\sum_j c_j = 0$. Dan meet het een verschil tussen groepen, niet een gemiddelde van groepen.

Voorbeelden: - $\mu_1 - \mu_2$: gewichten $(+1, -1, 0, 0)$ — paarsgewijs, som = 0 ✓ - $\mu_1 - \tfrac{1}{3}(\mu_2 + \mu_3 + \mu_4)$: gewichten $(1, -\tfrac{1}{3}, -\tfrac{1}{3}, -\tfrac{1}{3})$ — eerste vs de rest, som = 0 ✓ - $\mu_1 + \mu_2$: gewichten $(+1, +1, 0, 0)$ — som = 2, geen contrast (dit is een gemiddelde van twee)

In RMA-context heten contrasten ook wel geplande vergelijkingen als ze theoretisch gemotiveerd zijn. Polynomiale contrasten (lineair / kwadratisch / kubisch) zijn een speciale familie geordende contrasten — zie 6.2.a. De cross-link naar callout C van thema 5 (“Sjonge sjonge — al die woorden voor één ding”): een contrast is hetzelfde gereedschap als een gewogen optelling, een discriminant-functie, een lineaire combinatie. Alleen met de extra eis dat de gewichten optellen tot $0$.

6.2.0 Contrasten met de hand opbouwen — rustige bijles-volgorde

Voordat we polynomiale contrasten op de data loslaten via R, bouwen we de contrast-gewichten zelf op. Net zoals in de bijles bij Leiden’s drug-experiment, waar Jeanne stap voor stap de vier contrast-gewichten formuleerde voordat ze er statistiek op zette.

Stappenplan voor een geplande set contrasten:

Formuleer de onderzoeksvraag in normale taal. Bv: “Verschilt de ochtend van de drie latere momenten samen?”
Vertaal de vraag naar een vergelijking tussen groepen of groepen-van-groepen. Bv: ochtend versus gemiddelde van middag/avond/nacht.
Geef gewichten zodanig dat:
- de vergeleken kanten optellen tot $0$ (één kant +, één kant −),
- de absolute som van de positieve gewichten gelijk is aan de absolute som van de negatieve gewichten.
Check: $\sum_j c_j = 0$. Als de som $\neq 0$: geen contrast, geen vergelijking.
Herhaal voor elke volgende vraag. Bij $k$ niveaus kun je maximaal $k-1$ orthogonale contrasten neerzetten — drie vragen voor vier momenten.

Voorbeeld — drie geplande contrasten op vier momenten

Stel je vraagt drie dingen, in volgorde:

Vraag	In gewone taal	Gewichten ochtend / middag / avond / nacht
L₁	Ochtend vs. de rest — wijkt de ochtend af van het gemiddelde van middag, avond en nacht?	$(+3, -1, -1, -1)$
L₂	Middag vs. avond+nacht — binnen de niet-ochtend-momenten: middag vs gemiddeld-laat?	$(0, +2, -1, -1)$
L₃	Avond vs. nacht — binnen de twee laat-op-de-dag-momenten?	$(0, 0, +1, -1)$

Som van elk contrast (controle):

$L_1$: $3 + (-1) + (-1) + (-1) = 0$ ✓
$L_2$: $0 + 2 + (-1) + (-1) = 0$ ✓
$L_3$: $0 + 0 + 1 + (-1) = 0$ ✓

Orthogonaliteit — inproduct van elk paar gewichten moet $0$ zijn:

$L_1 \cdot L_2$: $3 \cdot 0 + (-1) \cdot 2 + (-1) \cdot (-1) + (-1) \cdot (-1) = 0 - 2 + 1 + 1 = 0$ ✓
$L_1 \cdot L_3$: $3 \cdot 0 + (-1) \cdot 0 + (-1) \cdot 1 + (-1) \cdot (-1) = 0 + 0 - 1 + 1 = 0$ ✓
$L_2 \cdot L_3$: $0 \cdot 0 + 2 \cdot 0 + (-1) \cdot 1 + (-1) \cdot (-1) = 0 + 0 - 1 + 1 = 0$ ✓

Alle drie paren = $0$ → de set is orthogonaal. Drie onafhankelijke vragen, elk een eigen stukje van de tussen-momenten-variantie.

Hoe je een contrast-waarde berekent uit groeps-gemiddelden

Voor elk contrast: vermenigvuldig elk gewicht met het bijbehorende groeps-gemiddelde en tel op:

\[L = c_1 \bar{Y}_1 + c_2 \bar{Y}_2 + c_3 \bar{Y}_3 + c_4 \bar{Y}_4\]

Stel onze kameleon-gemiddelden zijn (verzonnen voor uitleg): $\bar{Y}_{\text{ochtend}} = 60$, $\bar{Y}_{\text{middag}} = 30$, $\bar{Y}_{\text{avond}} = 50$, $\bar{Y}_{\text{nacht}} = 80$.

Dan:

$L_1 = 3 \cdot 60 + (-1) \cdot 30 + (-1) \cdot 50 + (-1) \cdot 80 = 180 - 30 - 50 - 80 = 20$. Positief — de ochtend zit gemiddeld boven de drie latere momenten samen (al is dat 20 punten op een schaal die loopt van — kijk in je data).
$L_2 = 0 \cdot 60 + 2 \cdot 30 + (-1) \cdot 50 + (-1) \cdot 80 = 60 - 50 - 80 = -70$. Sterk negatief — middag scoort gemiddeld véél lager dan avond+nacht.
$L_3 = 0 \cdot 60 + 0 \cdot 30 + 1 \cdot 50 + (-1) \cdot 80 = -30$. Negatief — avond zit onder nacht.

Interpretatie: $|L_2| = 70$ is veruit het grootst, dus de middag-versus-laat-tegenstelling is het dominantste patroon. Ochtend-versus-rest en avond-versus-nacht zijn ook present maar zwakker. Welke significant is hangt af van $SE_L$ en de bijbehorende $t$-toets — dat doet R voor je in 6.2.a.

Bij paarsgewijze vergelijkingen ($\mu_1 - \mu_2$ etc.) doe je in feite hetzelfde, maar met gewichten $(+1, -1, 0, \ldots, 0)$. Een paarsgewijze vergelijking is gewoon een speciaal-geval-contrast.

6.2.a Polynomiale contrasten — eerst hand, dan R

Polynomiale contrasten zijn een speciale familie

Bij geordende niveaus (tijd, dosis, leeftijd) kun je drie soorten patronen apart toetsen via een vaste set contrast-gewichten:

Trend	Wat het meet	Gewichten ochtend / middag / avond / nacht
Lineair $L$	monotone stijging of daling	$(-3, -1, +1, +3)$
Kwadratisch $Q$	U- of omgekeerde-U-vorm (piek in het midden)	$(+1, -1, -1, +1)$
Kubisch $C$	S-vorm (op-neer-op of neer-op-neer)	$(-1, +3, -3, +1)$

Deze drie zijn onderling orthogonaal (check zelf: $L \cdot Q = (-3)(1) + (-1)(-1) + (1)(-1) + (3)(1) = -3 + 1 - 1 + 3 = 0$ ✓). Samen vangen ze alle drie de $k - 1 = 3$ vrijheidsgraden van het moment-effect.

Bereken een polynomial-contrast met de hand:

Stel onze kameleon-gemiddelden: $\bar{Y}_{\text{ochtend}} = 60$, $\bar{Y}_{\text{middag}} = 30$, $\bar{Y}_{\text{avond}} = 50$, $\bar{Y}_{\text{nacht}} = 80$.

Lineair: $L = -3(60) - 1(30) + 1(50) + 3(80) = -180 - 30 + 50 + 240 = 80$. Sterk positief → over de dag stijgt het gemiddeld.
Kwadratisch: $Q = 1(60) - 1(30) - 1(50) + 1(80) = 60 - 30 - 50 + 80 = 60$. Positief → uiteinden hoger dan midden (U-vorm).
Kubisch: $C = -1(60) + 3(30) - 3(50) + 1(80) = -60 + 90 - 150 + 80 = -40$. Negatief → S-vorm met dip-in-het-midden.

Welke is dominant? Vergelijk absolute waardes: $|L| = 80 > |Q| = 60 > |C| = 40$. Lineair wint, kwadratisch tweede, kubisch derde. Pas in R krijg je daar $F$- en $p$-waardes bij (zie hieronder).

Voor de volledige polynomial-trends-opgave met vier kameleon-profielen — inclusief Dirks tricky $(3, 7, 5, 9)$-set waar het oog en het rekenen het oneens zijn — zie § 6.3.

In R — met `icontrasts = "contr.poly"`

In Leiden’s factorial_design()-route geef je polynomial-contrasten mee via icontrasts = c("contr.sum", "contr.poly") (al gedaan in 6.1.b). De polynomiale split komt vanzelf in de univariate output — summary() toont aparte rijen voor moment.L, moment.Q, moment.C.

Voor onze dieren.

# Met icontrasts = "contr.poly" splitst Anova het moment-effect automatisch
# in linear / quadratic / cubic. Bij hetzelfde factorial_design-object van 6.1.b:
summary(res_kam, multivariate = FALSE)


Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

             Sum Sq num Df Error SS den Df F value    Pr(>F)    
(Intercept) 1024395      1    41505     89  2196.6 < 2.2e-16 ***
moment        81118      3    33532    267   215.3 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Mauchly Tests for Sphericity

       Test statistic p-value
moment        0.95598 0.55674


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

        GG eps Pr(>F[GG])    
moment 0.97178  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

         HF eps   Pr(>F[HF])
moment 1.008238 5.813392e-71

De univariate tabel toont nu drie rijen onder moment — moment.L (lineair), moment.Q (kwadratisch), moment.C (kubisch). Per contrast krijg je $SS$, $df = 1$, $F$ en $p$.

De drie polynomiale gewichten op een rij

Voor vier geordende niveaus (gewichten orthogonaal-genormaliseerd):

Trend	Ochtend	Middag	Avond	Nacht
Lineair (L)	$-3/\sqrt{20}$	$-1/\sqrt{20}$	$1/\sqrt{20}$	$3/\sqrt{20}$
Kwadratisch (Q)	$1/2$	$-1/2$	$-1/2$	$1/2$
Kubisch (C)	$-1/\sqrt{20}$	$3/\sqrt{20}$	$-3/\sqrt{20}$	$1/\sqrt{20}$

Lineair: stijgend van laag naar hoog. Kwadratisch: hoog-laag-laag-hoog (U-vorm). Kubisch: laag-hoog-laag-hoog (S-vorm).

Toon code (lineair / kwadratisch / kubisch handmatig via t-test op contrast-scores)

# Met de hand om te zien hoe contrasten werken.
ochtend <- kameleon_momenten$ochtend
middag  <- kameleon_momenten$middag
avond   <- kameleon_momenten$avond
nacht   <- kameleon_momenten$nacht

c_lin <- c(-3, -1,  1, 3) / sqrt(20)
c_kw  <- c( 1, -1, -1, 1) / 2
c_kub <- c(-1,  3, -3, 1) / sqrt(20)

s_lin <- c_lin[1]*ochtend + c_lin[2]*middag + c_lin[3]*avond + c_lin[4]*nacht
s_kw  <- c_kw[1] *ochtend + c_kw[2] *middag + c_kw[3] *avond + c_kw[4] *nacht
s_kub <- c_kub[1]*ochtend + c_kub[2]*middag + c_kub[3]*avond + c_kub[4]*nacht

list(
  lineair       = t.test(s_lin),
  kwadratisch   = t.test(s_kw),
  kubisch       = t.test(s_kub)
)

$lineair

    One Sample t-test

data:  s_lin
t = 9.4275, df = 89, p-value = 4.892e-15
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
  9.440226 14.482223
sample estimates:
mean of x 
 11.96122 


$kwadratisch

    One Sample t-test

data:  s_kw
t = 24.865, df = 89, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 24.74069 29.03819
sample estimates:
mean of x 
 26.88944 


$kubisch

    One Sample t-test

data:  s_kub
t = -5.0015, df = 89, p-value = 0.000002837
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -8.289734 -3.575836
sample estimates:
mean of x 
-5.932785

Vragen 6.2

a) Welke polynomiale trends zijn significant? Met welke $t$-waardes?

b) Wat zegt het significante lineaire contrast over de dag-loop?

c) Wat zegt het significante kwadratische contrast?

d) Hoe interpreteer je het kubische contrast in deze data?

Antwoord 6.2 — open na je eigen poging

In gewone woorden.

Alle drie trends zijn significant: lineair $t(89) = 9.43$, kwadratisch $t(89) = 24.87$, kubisch $t(89) = -5.00$ — alle $p < .001$.
Het lineaire contrast wijst op een algemene stijging over de dag — je gemiddeld meer-zichzelf-zijn neemt toe van ochtend naar nacht.
Het kwadratische is het sterkst: er is een U-vormig patroon — hoog aan de uiteinden, dip in het midden. Hier vooral een dip in de middag, met ochtend en nacht relatief hoog.
Het kubische is kleiner maar significant: bovenop de stijging-met-dip zit een licht op-neer-op-neer-patroon (ochtend hoog, middag laag, avond medium, nacht weer hoog). De kameleon hapert in de middag, herpakt zich licht in de avond, en pas in de nacht is hij weer helemaal zichzelf.

APA-stijl.

Polynomiale contrasten over moment lieten alle drie trends significant zien. De kwadratische trend was het sterkst, $t(89) = 24.87$, $p < .001$ — een dip in de middag met hoge waarden aan beide uiteinden. De lineaire trend was eveneens significant, $t(89) = 9.43$, $p < .001$, met een algemene stijging over de dag. De kubische trend was kleiner maar wel significant, $t(89) = -5.00$, $p < .001$.

T5 — Geplande contrasten en orthogonaliteit

Vraag T5 — Pen-en-papier

Naast polynomiale contrasten kun je theorie-gedreven geplande contrasten opstellen. Stel: je vermoedt dat nacht uniek hoog is en de andere drie momenten samen lager.

a) Hoe ziet dat contrast eruit (gewichten over ochtend, middag, avond, nacht)?

b) Wat is het tweede orthogonale contrast als je daarna middag versus ochtend+avond wilt toetsen?

c) Wat betekent orthogonaliteit voor twee contrasten ($c_1 \cdot c_2 = 0$)?

d) Bij vier momenten kun je tot $k - 1 = 3$ orthogonale contrasten neerzetten. Verzin een derde contrast dat orthogonaal is op (a) en (b).

Antwoord T5 — open na je eigen poging

a) Nacht versus rest — gewichten $(-1, -1, -1, 3)$. Som = $0$ (vereiste voor een geldig contrast).

b) Middag versus ochtend+avond — gewichten $(-1, 2, -1, 0)$, met nacht uitgesloten (gewicht $0$). Som = $0$.

c) Inproduct van de gewichten = $0$. Voor (a) en (b): $(-1)(-1) + (-1)(2) + (-1)(-1) + (3)(0) = 1 - 2 + 1 + 0 = 0$. Orthogonaal — beide contrasten meten onafhankelijke stukken van het patroon.

d) Een mogelijkheid: ochtend versus avond — gewichten $(1, 0, -1, 0)$. Toets met (a): $(1)(-1) + (0)(-1) + (-1)(-1) + (0)(3) = -1 + 0 + 1 + 0 = 0$ ✓. Met (b): $(1)(-1) + (0)(2) + (-1)(-1) + (0)(0) = -1 + 0 + 1 + 0 = 0$ ✓. Drie orthogonale contrasten samen — de hele tussen-momenten-variantie zit verdeeld.

Inzicht. Geplande contrasten vóór je de data ziet zijn krachtiger dan post-hoc analyses; ze controleren type-I-fout zonder Bonferroni, mits je ze van tevoren specificeert. Bij vier momenten heb je drie kogels in je geweer.

Tentamen-tip — wegstreep-strategie. Op een tentamen krijg je vaak een set contrasten en moet je checken of ze orthogonaal zijn. Ortho-check is wegstrepen, geen rekenen:

Vermenigvuldig niet als één van de twee gewichten $0$ is — die term valt sowieso weg, sla over.
Vermenigvuldig niet uit als je in een oogopslag ziet dat termen tegen elkaar wegvallen (bv. $+\tfrac{1}{2} \cdot 1 + (-\tfrac{1}{2}) \cdot 1$ direct = $0$).
Begin onderaan in de tabel — examinatoren verstoppen niet-orthogonale paren vaak in latere rijen, en als de eerste paren orthogonaal zijn ben je nog steeds niet klaar.
Als alle paren naar $0$ wegstrepen, is de set orthogonaal — soms hoef je geen enkele vermenigvuldiging uit te rekenen.

6.2.b Paarsgewijze vergelijking met Bonferroni

Hoeveel paren bij $k$ niveaus?

Bij $k$ momenten zijn er $\binom{k}{2} = \dfrac{k(k-1)}{2}$ paarsgewijze vergelijkingen:

$k$	$\binom{k}{2}$	Bonferroni-drempel ($\alpha = .05$)
2	1	$.050$
3	3	$.0167$
4	6	$.0083$
5	10	$.0050$
6	15	$.0033$

Voor onze kameleon met vier momenten zijn dat zes paren: $(O,M), (O,A), (O,N), (M,A), (M,N), (A,N)$. Bonferroni-corrected $\alpha_{\text{per paar}} = .05 / 6 \approx .008$.

Algemene vorm — Leiden’s emmeans_test()-route.

emmeans_test() (uit rstatix) levert pairwise EMM-vergelijkingen met Bonferroni in één call. Maar onder de motorkap fit emmeans() het model — en voor RMA-pairwise heeft emmeans() een model nodig waarin moment als variabele in de modelformule staat met een Error-term per persoon. Dat is precies wat aov(y ~ moment + Error(id/moment)) levert. (Het mlm-object uit factorial_design() werkt hier niet — daar zit moment impliciet in de kolomvolgorde.)

We bouwen daarom twee modellen parallel: factorial_design() voor de multivariate-en-univariate RMA-output (sectie 6.1.b), en aov() met Error-term voor de pairwise vergelijkingen hier:

# Pairwise EMM-vergelijkingen met Bonferroni — model = aov(. + Error(id/within))
library(rstatix)
m_aov <- aov(score ~ moment + Error(id/moment), data = long_data)
em_pairs <- long_data |>
  emmeans_test(
    score ~ moment,
    p.adjust.method = "bonferroni",
    model           = m_aov
  )
em_pairs

Output: één rij per paar met estimate, SE, df, statistic ($t$-ratio), p.adj. De df is hier de gepoolde within-error-df ($(n-1)(k-1)$), niet de paarsgewijze paired-$t$-df ($n-1$) — RMA gebruikt de gepoolde error-schatting voor meer statistische power.

Voor onze dieren.

# Tweede model: aov() met Error(id/moment) — voorwaarde voor emmeans_test().
m_kam_aov <- aov(mate_van_zichzelf ~ moment + Error(id/moment),
                 data = kameleon_long)

em_pairs_kam <- kameleon_long |>
  emmeans_test(
    mate_van_zichzelf ~ moment,
    p.adjust.method = "bonferroni",
    model           = m_kam_aov
  )

em_pairs_kam |>
  select(group1, group2, df, statistic, p, p.adj, p.adj.signif) |>
  print()

# A tibble: 6 × 7
  group1  group2    df statistic        p    p.adj p.adj.signif
  <chr>   <chr>  <dbl>     <dbl>    <dbl>    <dbl> <chr>       
1 ochtend middag   267     16.1  4.08e-41 2.45e-40 ****        
2 ochtend avond    267      8.10 1.93e-14 1.16e-13 ****        
3 ochtend nacht    267     -8.02 3.39e-14 2.03e-13 ****        
4 middag  avond    267     -7.97 4.74e-14 2.84e-13 ****        
5 middag  nacht    267    -24.1  6.61e-69 3.97e-68 ****        
6 avond   nacht    267    -16.1  2.68e-41 1.61e-40 ****

Twee modellen, één design — waarom?

factorial_design() is gebouwd om de RMA-output te produceren: multivariate Pillai/Wilks-tabellen, univariate Type-III met sphericity-correctie. Daarvoor schrijft het achter de schermen een lm(cbind(t1, t2, t3, t4) ~ 1) (een mlm, multivariate lineair model) — moment zit impliciet in de kolomvolgorde van de DV-matrix.

emmeans() daarentegen wil moment als variabele in de modelformule zien, niet als kolom-index. Daarom bouw je voor de pairwise stap een tweede model met aov(. + Error(id/moment)) op de long-vorm. Beide modellen schatten dezelfde gepoolde within-error; de Type-III-Sum-of-Squares (sectie 6.1.b) en de pairwise contrasts (hier) komen overeen.

Geen dubbel werk in de analyse-stap — wel twee modellen in de R-code, omdat de twee output-vormen elk hun eigen model-structuur nodig hebben. Op het tentamen kan je naar de syntax van beide worden gevraagd.

Klassieke paired-$t$-route — alternatief

Voor wie geen emmeans_test() wil gebruiken: dezelfde zes paren via gewone t.test()-aanroepen plus handmatige Bonferroni-correctie:

ps <- c(
  o_m = t.test(kameleon_momenten$ochtend, kameleon_momenten$middag, paired = TRUE)$p.value,
  o_a = t.test(kameleon_momenten$ochtend, kameleon_momenten$avond,  paired = TRUE)$p.value,
  o_n = t.test(kameleon_momenten$ochtend, kameleon_momenten$nacht,  paired = TRUE)$p.value,
  m_a = t.test(kameleon_momenten$middag,  kameleon_momenten$avond,  paired = TRUE)$p.value,
  m_n = t.test(kameleon_momenten$middag,  kameleon_momenten$nacht,  paired = TRUE)$p.value,
  a_n = t.test(kameleon_momenten$avond,   kameleon_momenten$nacht,  paired = TRUE)$p.value
)
round(p.adjust(ps, method = "bonferroni"), 4)

Wiskundig gelijk; alleen niet wat Leiden’s R-practical voorschrijft.

Vragen 6.2

e) Hoeveel paren toets je bij $k = 4$ momenten? Wat is $\alpha_{\text{per paar}}$ na Bonferroni?

f) Welke paren zijn significant na Bonferroni-correctie?

g) Welk paar heeft de grootste $t$-waarde? Wat zegt dat?

Antwoord 6.2 — open na je eigen poging

In gewone woorden.

Bij $k = 4$ momenten zijn er $\binom{4}{2} = 6$ paren. Bonferroni: $\alpha_{\text{per paar}} = .05/6 \approx .008$.
Alle zes paren overleven Bonferroni — alle $p$-Bonf-corrected $< .001$. Alle vier momenten verschillen onderling significant, in alle combinaties.
Het sterkste paar is middag versus nacht: $t(89) = -23.02$ — het verschil is bijna $40$ punten ($M_{\text{middag}} \approx 33$ versus $M_{\text{nacht}} \approx 73$). Dat past bij ons inzicht uit de polynomiale contrasten: de dip-in-middag en de piek-in-nacht zijn de twee uiteinden van het patroon.

APA-stijl.

Paarsgewijze paired $t$-toetsen met Bonferroni-correctie ($\alpha_{\text{per paar}} = .008$) lieten zien dat alle vier momenten onderling significant verschilden (alle $p < .001$). Het sterkste verschil zat tussen middag en nacht, $t(89) = -23.02$, $p < .001$, met de kameleons gemiddeld $40$ schaalpunten meer-zichzelf in de nacht.

T6 — Bonferroni voor RMA-paren handrekenen

Vraag T6 — Pen-en-papier

Bij $k = 5$ momenten doe je paarsgewijze vergelijkingen met Bonferroni-correctie.

a) Hoeveel paren zijn er?

b) Wat is $\alpha_{\text{per paar}}$?

c) Een paar heeft $p = .009$. Significant na Bonferroni?

d) Een paar heeft $p = .003$. Significant na Bonferroni?

Antwoord T6 — open na je eigen poging

a) $\binom{5}{2} = 5 \cdot 4 / 2 = 10$ paren.

b) $\alpha_{\text{per paar}} = .05 / 10 = .005$.

c) $p = .009 > .005$ — niet significant na Bonferroni. (Wel significant op de gewone $.05$.)

d) $p = .003 < .005$ — wel significant na Bonferroni.

Inzicht. Bonferroni voor RMA-paren: hetzelfde mechanisme als bij MANOVA in thema 5, alleen tellen we nu paren-binnen-een-factor in plaats van DV’s. Bij veel momenten wordt de drempel snel streng — overweeg dan polynomiale contrasten als alternatief.

6.2.c Visualisatie — gemiddeld profiel + spaghetti

Toon code (interaction-plot stijl: groepsgemiddelde + 95% CI + spaghetti)

means_v <- colMeans(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])
sds_v   <- apply(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")], 2, sd)
n_k     <- nrow(kameleon_momenten)
ses_v   <- sds_v / sqrt(n_k)
ci_lo   <- means_v - 1.96 * ses_v
ci_hi   <- means_v + 1.96 * ses_v

matY <- as.matrix(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])

par(mar = c(4, 4.5, 2.5, 1))
matplot(t(matY), type = "l", lty = 1,
        col = adjustcolor("#888888", alpha.f = 0.30),
        xaxt = "n", xlab = "Moment", ylab = "Mate van zichzelf",
        ylim = c(0, 100),
        main = "Trajecten + gemiddeld profiel")
axis(1, at = 1:4, labels = c("Ochtend", "Middag", "Avond", "Nacht"))

# Gemiddeld profiel met error bars.
arrows(1:4, ci_lo, 1:4, ci_hi, length = 0.05, angle = 90, code = 3,
       col = "#CC3311", lwd = 2)
lines(1:4, means_v, col = "#CC3311", lwd = 3)
points(1:4, means_v, col = "#CC3311", pch = 16, cex = 1.6)

legend("bottomleft",
       legend = c("Individuele kameleons", "Gemiddelde + 95% CI"),
       col    = c(adjustcolor("#888888", alpha.f = 0.7), "#CC3311"),
       lty    = 1,
       lwd    = c(1, 3),
       bty    = "n",
       cex    = 0.85)

6.3 Polynomial-trends-surprise — vier dagdeel-profielen

Eerst raden, dan rekenen — wat zegt het oog, wat zegt de optelling?

“Vier kameleons, ieder met een eigen dagritme,” vertelde de oude kameleon. “Sommigen volgen een rechte lijn. Sommigen zakken in het midden in. Sommigen springen op en neer. En een paar lijken iets te doen wat niet helemaal is wat het lijkt.”

In deze opgave krijg je vier dagdeel-profielen, één per fictieve kameleon. Per profiel: het gemiddelde mate-van-zichzelf op de vier momenten ochtend / middag / avond / nacht. Voor je gaat rekenen, kijk eerst goed naar elk profiel en gok welke polynomiale trend dominant is: lineair, kwadratisch, kubisch, of geen.

Daarna reken je voor elk profiel de drie contrasten uit en kijk je of je intuïtie klopt. Bij minstens één profiel zit een verrassing.

Wat ga je rekenen?

Drie polynomiale contrasten op vier geordende niveaus, met de niet-genormaliseerde gewichten (handrekenbaar):

Lineair $L_1$: gewichten $(-3, -1, +1, +3)$ — meet stijging of daling
Kwadratisch $L_2$: gewichten $(+1, -1, -1, +1)$ — meet U-vorm of omgekeerde U
Kubisch $L_3$: gewichten $(-1, +3, -3, +1)$ — meet S-vorm (op-neer-op of neer-op-neer)

Per profiel: $L_k = \sum_j c_{kj} \cdot \bar{Y}_j$. Rangschik op absolute waarde — de grootste $|L_k|$ is de dominante trend.

De vier profielen

Toon code (gt-tabel met vier profielen)

profielen <- data.frame(
  Kameleon = c("Anita",        "Berend",     "Clara",       "Dirk"),
  ochtend  = c( 5,              0,             9,             3),
  middag   = c( 5,              3,             7,             7),
  avond    = c( 5,              4,             4,             5),
  nacht    = c( 5,              1,             0,             9)
)

profielen |>
  gt() |>
  tab_header(title = "Vier kameleon-profielen — gemiddelde mate-van-zichzelf per moment") |>
  cols_label(Kameleon = "Kameleon",
             ochtend  = "Ochtend",
             middag   = "Middag",
             avond    = "Avond",
             nacht    = "Nacht") |>
  tab_options(table.font.size = px(13),
              table.border.top.style    = "none",
              table.border.bottom.style = "none")

Kameleon	Ochtend	Middag	Avond	Nacht
Vier kameleon-profielen — gemiddelde mate-van-zichzelf per moment
Anita	5	5	5	5
Berend	0	3	4	1
Clara	9	7	4	0
Dirk	3	7	5	9

Visualisatie — vier profielen in één plot

Toon code (lineplot vier profielen)

profielen_long <- data.frame(
  Kameleon = rep(profielen$Kameleon, each = 4),
  moment   = factor(rep(c("ochtend", "middag", "avond", "nacht"), times = 4),
                    levels = c("ochtend", "middag", "avond", "nacht")),
  waarde   = c(t(profielen[, 2:5]))
)

par(mar = c(4, 4, 2.5, 5), xpd = TRUE)
plot(NA, xlim = c(1, 4), ylim = c(-1, 11),
     xaxt = "n", xlab = "Moment", ylab = "Mate-van-zichzelf",
     main = "Vier kameleon-profielen")
axis(1, at = 1:4, labels = c("ochtend", "middag", "avond", "nacht"))

kleuren <- c(Anita = "#0077BB", Berend = "#33BBEE", Clara = "#EE7733", Dirk = "#CC3311")

for (k in profielen$Kameleon) {
  waardes <- as.numeric(profielen[profielen$Kameleon == k, 2:5])
  lines(1:4, waardes, col = kleuren[k], lwd = 2.5)
  points(1:4, waardes, col = kleuren[k], pch = 19, cex = 1.2)
}

legend("right", inset = c(-0.18, 0), legend = names(kleuren),
       col = kleuren, lwd = 2.5, bty = "n", cex = 0.9)

Vraag 6.3.a — Eerst kijken en raden (pen-en-papier, geen rekenen)

Voor elk van de vier kameleons: welke polynomiale trend lijkt visueel dominant?

Anita: _________________________
Berend: _________________________
Clara: _________________________
Dirk: _________________________

Geef per kameleon je beste gok: lineair, kwadratisch, kubisch, of geen trend. Schrijf ook kort op waarom je dat ziet.

6.3.a Reken de drie contrasten uit

L1 <- c(-3, -1,  1,  3)  # lineair
L2 <- c( 1, -1, -1,  1)  # kwadratisch
L3 <- c(-1,  3, -3,  1)  # kubisch

bereken_contrasten <- function(profile) {
  c(L1 = sum(L1 * profile),
    L2 = sum(L2 * profile),
    L3 = sum(L3 * profile))
}

resultaten <- t(apply(profielen[, 2:5], 1, bereken_contrasten))
rownames(resultaten) <- profielen$Kameleon

cat("Contrast-waardes per kameleon:\n")

Contrast-waardes per kameleon:

print(resultaten)

        L1 L2 L3
Anita    0  0  0
Berend   4 -6 -2
Clara  -30 -2  0
Dirk    16  0 12

cat("\nAbsolute waardes (sterkste eerst per rij):\n")


Absolute waardes (sterkste eerst per rij):

rangorde <- t(apply(abs(resultaten), 1, function(x) names(x)[order(-x)]))
rownames(rangorde) <- profielen$Kameleon
print(rangorde)

       [,1] [,2] [,3]
Anita  "L1" "L2" "L3"
Berend "L2" "L1" "L3"
Clara  "L1" "L2" "L3"
Dirk   "L1" "L3" "L2"

Vragen 6.3.b — Rekenen en rangschikken

a) Wat zijn de drie contrast-waardes voor Berend $(0, 3, 4, 1)$? Welke trend is dominant?

b) Wat zijn de drie contrast-waardes voor Clara $(9, 7, 4, 0)$? Welke trend is dominant?

c) Wat zijn de drie contrast-waardes voor Dirk $(3, 7, 5, 9)$? Welke trend is dominant? Komt dit overeen met je gok in vraag 6.3.a?

d) Wat zijn de drie contrast-waardes voor Anita $(5, 5, 5, 5)$? Wat zegt dat?

e) Rangschik per kameleon op absolute waarde. Welke trend domineert, welke is tweede, welke derde?

Antwoord 6.3.b — open na je eigen poging

Anita $(5, 5, 5, 5)$:

$L_1 = -3(5) - 1(5) + 1(5) + 3(5) = -15 - 5 + 5 + 15 = 0$
$L_2 = 1(5) - 1(5) - 1(5) + 1(5) = 5 - 5 - 5 + 5 = 0$
$L_3 = -1(5) + 3(5) - 3(5) + 1(5) = -5 + 15 - 15 + 5 = 0$

Alle drie nul — er is geen trend. Constante waarde laat geen polynomial-trend zien. Anita is een rust-kameleon.

Berend $(0, 3, 4, 1)$:

$L_1 = -3(0) - 1(3) + 1(4) + 3(1) = 0 - 3 + 4 + 3 = 4$
$L_2 = 1(0) - 1(3) - 1(4) + 1(1) = 0 - 3 - 4 + 1 = -6$
$L_3 = -1(0) + 3(3) - 3(4) + 1(1) = 0 + 9 - 12 + 1 = -2$

Absoluut: $|L_2| = 6 > |L_1| = 4 > |L_3| = 2$. Kwadratisch dominant — Berend heeft een piek in het midden (middag-avond hoog, ochtend en nacht laag). Negatief kwadratisch contrast = omgekeerde U-vorm.

Clara $(9, 7, 4, 0)$:

$L_1 = -3(9) - 1(7) + 1(4) + 3(0) = -27 - 7 + 4 + 0 = -30$
$L_2 = 1(9) - 1(7) - 1(4) + 1(0) = 9 - 7 - 4 + 0 = -2$
$L_3 = -1(9) + 3(7) - 3(4) + 1(0) = -9 + 21 - 12 + 0 = 0$

Absoluut: $|L_1| = 30 \gg |L_2| = 2 > |L_3| = 0$. Lineair dominant, sterk negatief — Clara daalt monotoon over de dag.

Dirk $(3, 7, 5, 9)$:

$L_1 = -3(3) - 1(7) + 1(5) + 3(9) = -9 - 7 + 5 + 27 = 16$
$L_2 = 1(3) - 1(7) - 1(5) + 1(9) = 3 - 7 - 5 + 9 = 0$
$L_3 = -1(3) + 3(7) - 3(5) + 1(9) = -3 + 21 - 15 + 9 = 12$

Absoluut: $|L_1| = 16 > |L_3| = 12 \gg |L_2| = 0$. Lineair dominant, daarna kubisch — Dirk stijgt over de hele dag, met een S-vormige rimpel erin.

Surprise — vraag e) over Dirk. Visueel ziet Dirks profiel er wisselend uit: hij gaat omhoog (3 → 7), zakt in (7 → 5), schiet weer omhoog (5 → 9). Het oog leest dat als kubisch-dominante S-vorm. Maar de rekenkundige sterkste trend is lineair — over de hele dag stijgt hij. De kubische component is er wel, maar minder sterk dan de lineaire.

Inzicht: het oog ziet bewegingen op de korte termijn (op-neer-op), het rekenen registreert de grote lijn (begin laag, eind hoog). Allebei waar, allebei in de drie contrasten verdeeld. Rapporteer in volgorde van absolute sterkte, niet in volgorde van rang. Voor Dirk: “Lineaire stijging dominant, met aanvullende kubische rimpel.”

In Leiden’s exercise book hoofdstuk 7 staat een vergelijkbare opgave met set $(3, 7, 5, 9)$ — daar is dezelfde verrassing besproken. Onze kameleon-versie Dirk is een 1-op-1-vertaling met een naam erbij.

Reken-veilig — let op tekens

Bij polynomiale-gewichten staat in $L_3$ vooraan een min: $-1 \cdot \bar{Y}_1$ is niet $+\bar{Y}_1$. Klein teken-foutje zorgt regelmatig voor de verkeerde rangschikking. Controle: som van gewichten = 0 voor elk contrast.

$L_1$: $-3 - 1 + 1 + 3 = 0$ ✓
$L_2$: $1 - 1 - 1 + 1 = 0$ ✓
$L_3$: $-1 + 3 - 3 + 1 = 0$ ✓

En altijd: contrasten zijn gewogen verschillen, geen gemiddelden. De gewichten optellen tot 0 is het bewijs dat het inderdaad om een vergelijking gaat.

T7 — RMA-aannames

Vraag T7 — Pen-en-papier

a) Welke aanname maakt RMA via de multivariate route over de within-subject vector?

b) Hoe check je dat visueel bij vier momenten?

c) Wat is een serieuze multivariate outlier in een trajectset?

d) Bij geconstateerde schending van multivariate normaliteit: welke twee opties heb je?

Antwoord T7 — open na je eigen poging

a) De gezamenlijke verdeling van de within-subject vector $(Y_1, Y_2, Y_3, Y_4)$ is multivariaat normaal. Sterker dan “elk moment apart normaal” — alle lineaire combinaties moeten ook normaal verdeeld zijn.

b) Per moment een histogram of QQ-plot; en in totaal een chi-square QQ-plot van de Mahalanobis-afstanden tussen elk dier en de gemiddelde vector. Bij sterke afwijkingen van de diagonaal: schending. Of inspecteer de spaghetti-plot — extreem afwijkende trajecten vallen op.

c) Een traject dat in vorm én niveau sterk afwijkt — bv. één kameleon die op alle momenten extreem laag scoort terwijl anderen variëren. Een Mahalanobis $D^2 > \chi^2_{k, .999}$ is een vuistregel. Bij $k = 4$: $\chi^2_{4, .999} \approx 18.47$.

d) (1) Transformatie als de data scheef is (log, sqrt) — past beter bij multivariate normaliteit. (2) Robuuste varianten of bootstrap. Bij gebalanceerde $n$ en $n_i \geq 20$ is RMA-$F$ overigens robuust tegen lichte schendingen, dus paniek is zelden nodig.

Inzicht. RMA’s aannames stapelen MANOVA’s aannames op binnen-persoon-niveau. Grote $n$ + visuele check + één robuustheids-melding in de discussie is in de praktijk genoeg.

T8 — Missing data in RMA

Vraag T8 — Pen-en-papier

Een proefpersoon mist één van de vier metingen. Wat doe je?

a) Wat doet lm(cbind(...) ~ 1) standaard met een rij waarin één moment ontbreekt?

b) Wat is de simpelste oplossing?

c) Wat zijn twee betere oplossingen?

d) Onder welke aanname werkt multiple imputation?

Antwoord T8 — open na je eigen poging

a) Standaard: listwise deletion — de hele rij valt eruit. Bij vier momenten en één missend gegeven verlies je dus alle vier de metingen voor die kameleon.

b) Behoud van design en complete-case-analyse: gewoon doorgaan met $n - 1$. Power-verlies, geen bias als de missing-data-mechanisme MCAR (missing completely at random) is.

c) (1) Mixed-effects model met lme4::lmer(mate ~ moment + (1|id)) — hanteert per-meting missingness automatisch en gebruikt alle beschikbare data. (2) Multiple imputation met mice of Amelia — vul de missing waardes meerdere keren plausibel in, run RMA op elke gevulde dataset, pool de resultaten.

d) Multiple imputation veronderstelt MAR (missing at random): de missing-data-kans hangt af van waargenomen variabelen, niet van de niet-waargenomen waarde zelf. Sterker dan MCAR (kans op missing onafhankelijk van alles), zwakker dan MNAR (kans hangt af van de missende waarde).

Inzicht. Bij RMA’s standaard-route ben je gevoelig voor missing data; bij mixed-effects-modellen wordt het automatisch beter aangepakt. Zie Wat blijft liggen.

T9 — RMA versus mixed-effects models

Vraag T9 — Pen-en-papier

Mixed-effects models (random effects op proefpersoon-niveau) zijn een uitbreiding van RMA. Wanneer kies je welke?

a) Wanneer is RMA voldoende?

b) Wanneer kies je voor mixed-effects (lme4::lmer)?

c) Wat is een random intercept in deze context?

d) Wat is een random slope?

Antwoord T9 — open na je eigen poging

a) Bij gebalanceerd design (alle proefpersonen op alle momenten gemeten), geen missing data, en als je vooral geïnteresseerd bent in gemiddelde binnen-persoon-effecten (niet in tussen-persoon-variabiliteit van die effecten). RMA is dan precies, simpel en standaard.

b) Bij onbalans (verschillende aantallen metingen per persoon), missing data, continue tijd-variabele (in plaats van discrete momenten), of als je tussen-persoon-variatie in trajecten wilt modelleren (“verschillen kameleons in hun ochtend-middag-helling?”). Mixed-effects geeft per-proefpersoon-schattingen, RMA niet.

c) Random intercept = elke kameleon heeft zijn eigen baseline-niveau. Sommige zijn algemeen hoog (“zelfgevoel-rijkere” kameleons), anderen algemeen laag. Het random intercept legt vast hoeveel die baselines variëren.

d) Random slope = elke kameleon heeft zijn eigen ochtend-middag-helling (of moment-trend). Niet alleen variërend in hoogte, maar in vorm van het traject. Sommige kameleons dippen extreem in de middag, anderen nauwelijks.

Inzicht. Mixed-effects is RMA + “per persoon variatie in effecten”. Voor nu: RMA als hoofdroute. Bij design-uitbreiding (geneste metingen, onbalans, missing) ga je verder; zie Wat blijft liggen en lme4.

6.A Uitbreiding — mixed-design (between $\times$ within)

“Wat als ik ook habitat meeneem?”

“Misschien,” zei de kameleon, “verschilt de mate-van-zichzelf niet alleen tussen momenten, maar ook tussen habitats — bos versus woestijn — en is het effect van moment anders in een bos dan in een woestijn. Een interactie tussen between en within.”

In een mixed-design RMA komt één between-subject factor (verschillende kameleons in elke groep) bij de within-subject factor (alle kameleons op alle momenten). Drie effecten te toetsen: het hoofdeffect habitat, het hoofdeffect moment, en de interactie habitat $\times$ moment.

Algemene vorm.

# In factorial_design(): tussen-persoon-factor toegevoegd via `between =`.
fd_mix <- factorial_design(
  data    = long_data,
  dv      = score,
  wid     = id,
  between = habitat,
  within  = moment
)

# Anova zoals voorheen — nu krijg je drie multivariate F's: habitat, moment, habitat:moment.
res_mix <- Anova(
  fd_mix$model,
  idata      = fd_mix$idata,
  idesign    = fd_mix$idesign,
  icontrasts = c("contr.sum", "contr.poly"),
  type       = "III"
)
summary(res_mix, multivariate = TRUE)

Wat verandert er bij een mixed-design?

Drie multivariate $F$-toetsen: voor habitat (between), moment (within), en habitat:moment (interactie).
Bij significante interactie: het hoofd-effect heeft minder betekenis en mag niet los geïnterpreteerd worden — qualified by interaction. Bekijk dan binnen elk niveau van habitat de moment-trend afzonderlijk.
Type-III SS wordt belangrijk bij onbalans + interactie. Vandaar Anova(., type = 3) met contr.sum-contrasten.
Tussen-persoon factor: gewone Type-III $F$ via summary, alleen op de habitat-rij.
Univariate route bij mixed-design: vraagt extra sphericity-aannames (binnen elk habitat-niveau). Multivariate route is en blijft sphericity-vrij.

In dit werkboek gaan we hier niet diep op in — voor de meeste tentamen-opgaven volstaat de eenwegs-RMA. Voor bredere ervaring: zie de aanbevolen literatuur in Wat blijft liggen.

T10 — Between $\times$ within: wanneer interactie interessant?

Vraag T10 — Pen-en-papier

In een mixed-design RMA met habitat (bos/woestijn) en moment (vier momenten) op de mate-van-zichzelf:

a) Wat betekent een significant hoofdeffect van moment op de mate-van-zichzelf zonder interactie met habitat?

b) Wat betekent een significante interactie van habitat $\times$ moment op de mate-van-zichzelf?

c) Een onderzoeker rapporteert: “hoofdeffect van moment significant, interactie habitat $\times$ moment significant” (DV impliciet: mate-van-zichzelf). Hoe interpreteer je dit?

d) Wat doe je als follow-up bij een significante interactie habitat $\times$ moment op de mate-van-zichzelf?

Antwoord T10 — open na je eigen poging

a) Het effect van moment op de mate-van-zichzelf bestaat en is voor beide habitats gelijk — bos-kameleons en woestijn-kameleons volgen parallelle trajecten over de dag in de mate-van-zichzelf.

b) Het effect van moment op de mate-van-zichzelf verandert in grootte en/of richting met het niveau van habitat — bos-kameleons hebben een ander dagpatroon in de mate-van-zichzelf dan woestijn-kameleons. Geen parallelle trajecten meer.

c) Het hoofdeffect van moment op de mate-van-zichzelf mag niet los geïnterpreteerd worden — qualified by interaction met habitat. Het algemene moment-patroon is een gemiddelde van twee verschillende habitat-patronen. Volgende stap: split-by-habitat-analyse — binnen elk habitat-niveau het effect van moment op de mate-van-zichzelf afzonderlijk schatten.

d) Per niveau van de between-factor habitat afzonderlijk het effect van moment op de mate-van-zichzelf toetsen (“simple effects within habitat” — wij noemen dat hier polynomiale contrasten per habitat). Of polynomiale contrasten op de interactie-term zelf: heeft de habitat-verschil-curve in de mate-van-zichzelf zelf een significante lineaire / kwadratische trend?

Inzicht. Mixed-design RMA’s grootste rijkdom zit in de interactie-term. Hoofdeffecten op de DV zijn vaak voorspelbaar; de interactie laat zien of de tijdsdynamiek in de DV habitat-afhankelijk is.

R-spiekblad bij RMA

Alle commando’s op één plek

Pakketten en data laden — wide en long format

# Pakketten activeren — één keer per sessie.
library(tidyr)       # pivot_longer, pivot_wider
library(dplyr)       # |> pipe + mutate
library(rstatix)     # factorial_design, emmeans_test
library(car)         # Anova met idata/idesign/icontrasts
library(emmeans)     # emmeans() voor interactie-EMM (mixed design)

# Hoofd-dataset laden (wide format).
load("data/kameleon_momenten.RData")
str(kameleon_momenten)

# Wide → long via pivot_longer.
kameleon_long <- kameleon_momenten |>
  pivot_longer(
    cols      = c(ochtend, middag, avond, nacht),
    names_to  = "moment",
    values_to = "mate_van_zichzelf"
  ) |>
  mutate(moment = factor(moment,
                         levels = c("ochtend", "middag", "avond", "nacht")))

# Long → wide via pivot_wider (terugweg).
kameleon_wide_again <- kameleon_long |>
  pivot_wider(names_from = moment, values_from = mate_van_zichzelf)

Verkenning — spaghetti, correlaties, boxplot

# Means + SD per moment.
colMeans(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])
apply(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")], 2, sd)

# Correlatiematrix tussen momenten — within-subject samenhang.
cor(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])

# Boxplot per moment.
boxplot(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])

# Spaghetti-plot — individuele trajecten.
matY <- as.matrix(kameleon_momenten[, c("ochtend", "middag", "avond", "nacht")])
matplot(t(matY), type = "l")

RMA fitten — Leiden’s `factorial_design()` + `Anova()`

# Stap 1: factorial_design wrapper definieert het RMA-design.
fd_kam <- factorial_design(
  data   = kameleon_long,
  dv     = mate_van_zichzelf,
  wid    = id,
  within = moment
)

# Stap 2: Anova() met icontrasts ALS argument (NIET via options()).
res_kam <- Anova(
  fd_kam$model,
  idata      = fd_kam$idata,
  idesign    = fd_kam$idesign,
  icontrasts = c("contr.sum", "contr.poly"),
  type       = "III"
)

# Stap 3: samenvatting — multivariate én univariate samen.
summary(res_kam, multivariate = TRUE)

De vier teststatistics — Pillai, Wilks, Hotelling-Lawley, Roy

# Alle vier komen tegelijk uit summary(., multivariate = TRUE).
# Bij intercept-only RMA (geen between-factor) zijn ze meestal identiek.
# Standaard rapporteer: Pillai's V als primary.
#
# RICHTING (mnemonic):
#   Pillai V      — hoger = sterker (range 0 tot min(I-1, p))
#   Wilks Λ       — LAGER = sterker (range 0 tot 1; tegenovergesteld aan F!)
#   Hotelling T²  — hoger = sterker (geen bovengrens)
#   Roy θ         — hoger = sterker (alleen grootste eigenvalue)
# Bij intercept-only-RMA komen alle vier op dezelfde F-waarde uit.

Polynomiale contrasten — `contr.poly()`

# Polynomiale contrasten op een geordende within-subject factor.
# In de Leiden-route gebeurt dit automatisch via icontrasts =:
res_poly <- Anova(fd_kam$model,
                  idata     = idata,
                  idesign   = ~moment,
                  icontrasts = c("contr.sum", "contr.poly"),
                  type       = "III")
summary(res_poly, multivariate = TRUE)

# Polynomial-gewichten voor k = 4 (Leiden-scope, hand-rekenen):
#
#   lineair     (-3, -1, +1, +3)
#   kwadratisch (+1, -1, -1, +1)
#   kubisch     (-1, +3, -3, +1)
#
# In R: contr.poly(k) geeft de genormaliseerde versie als matrix.

Polynomial-gewichten voor $k = 3$ en $k = 5$ — voor wie verder wil

Buiten Leiden-scope, maar handig bij ander cursus-materiaal:

# k = 3:  lineair    (-1,  0, +1)
#         kwadratisch (+1, -2, +1)
#
# k = 5:  lineair    (-2, -1,  0, +1, +2)
#         kwadratisch (+2, -1, -2, -1, +2)
#         kubisch    (-1, +2,  0, -2, +1)
#         quartisch  (+1, -4, +6, -4, +1)

Bij MVDA Leiden komt alleen $k = 4$ expliciet terug in de opdrachten — contr.poly(k) doet de rest voor je in R.

Geplande contrasten + Bonferroni voor paren

# Leiden-route: emmeans_test() met aov(. + Error(id/within))-model.
m_aov <- aov(mate_van_zichzelf ~ moment + Error(id/moment),
             data = kameleon_long)
em_pairs <- kameleon_long |>
  emmeans_test(
    mate_van_zichzelf ~ moment,
    p.adjust.method = "bonferroni",
    model           = m_aov
  )

# Alternatief: paarsgewijze paired t-toetsen + handmatige Bonferroni.
# Wiskundig gelijk; Leiden vraagt naar de emmeans_test-syntax.
ps <- c(
  o_m = t.test(kameleon_momenten$ochtend, kameleon_momenten$middag, paired = TRUE)$p.value,
  o_a = t.test(kameleon_momenten$ochtend, kameleon_momenten$avond,  paired = TRUE)$p.value,
  o_n = t.test(kameleon_momenten$ochtend, kameleon_momenten$nacht,  paired = TRUE)$p.value,
  m_a = t.test(kameleon_momenten$middag,  kameleon_momenten$avond,  paired = TRUE)$p.value,
  m_n = t.test(kameleon_momenten$middag,  kameleon_momenten$nacht,  paired = TRUE)$p.value,
  a_n = t.test(kameleon_momenten$avond,   kameleon_momenten$nacht,  paired = TRUE)$p.value
)
p.adjust(ps, method = "bonferroni")

Mixed-effects alternatief (kort, voor wie verder wil)

# Voor onbalans, missing data, of random slopes: lme4.
# library(lme4); library(lmerTest)
# fit_me <- lmer(mate ~ moment + (1 | id), data = long_kam)
# summary(fit_me)
# anova(fit_me)

Visualisatie — spaghetti + gemiddelde lijn

# ggplot variant — long format vereist.
# library(ggplot2)
# ggplot(long_kam, aes(x = moment, y = mate)) +
#   geom_line(aes(group = id), color = "grey", alpha = 0.3) +
#   stat_summary(fun = mean, geom = "line", color = "red", size = 1.2,
#                aes(group = 1)) +
#   stat_summary(fun.data = mean_cl_normal, geom = "errorbar",
#                color = "red", width = 0.1) +
#   labs(x = "Moment", y = "Mate van zichzelf")

Voorbeeld-tentamenvragen

Even oefenen op tentamen-toon

Onderaan dit hoofdstuk staan vier korte theorievragen in tentamen-stijl en één mini R-opdracht. Geen verhalend frame meer — student-aan-tentamen-modus. Hou de vuistregels paraat: vier teststatistics (Pillai primary), polynomiale contrasten bij geordende niveaus, Bonferroni $\alpha/k(k-1)/2$ voor paren, multivariate route is sphericity-vrij.

Theorie en handreken

Vraag E1 — Pillai uit RMA-output interpreteren

Een onderzoeker meet bij $60$ eenden het rust-niveau op drie momenten van het broedseizoen (vroeg / midden / laat). RMA-output:

Multivariate Tests: moment op rust-niveau
                 Df test stat approx F num Df den Df    Pr(>F)
Pillai            1  0.624     48.20      2     58    < .001
Wilks             1  0.376     48.20      2     58    < .001
Hotelling-Lawley  1  1.662     48.20      2     58    < .001
Roy               1  1.662     48.20      2     58    < .001

Welke uitspraak is correct?

De vier teststatistics geven verschillende $p$-waardes; we moeten Hotelling-Lawley als primary kiezen.
Pillai’s $V = 0.624$, $F(2, 58) = 48.20$, $p < .001$ — een groot multivariaat effect (Cohens drempel $V \geq .14$); het rust-niveau verschilt sterk tussen broedmomenten.
Wilks’ $\Lambda = 0.376$ betekent dat slechts $38\%$ van de variantie in het rust-niveau verklaard wordt; klein effect.
De gelijke $F$-waardes verraden een schending van multivariate normaliteit.

Antwoord E1 — open na je eigen poging

b) Allemaal in dezelfde richting, Pillai is de standaard primary, $V = .62 \gg .14$ dus groot effect van moment op het rust-niveau. Optie a is fout (alle vier identiek door intercept-only-RMA — niet anders kiezen). Optie c verwart Wilks (lager = sterker effect) — $\Lambda = .38$ betekent juist groot effect van moment op het rust-niveau ($1 - \Lambda = .62$ als ruwe $\eta^2$-analoog). Optie d is onzinnig: bij RMA met één within-factor en geen between komen Pillai/Wilks/Hotelling/Roy wiskundig op hetzelfde $F$ uit, geen aanname-schending.

Vraag E2 — Polynomiaal contrast significant — wat betekent het kubische?

Een onderzoeker meet bij dezelfde proefpersonen de reactietijd op vier dosissen ($t_1$ tot $t_4$) van een prikkel. RMA met polynomiale contrasten geeft:

Polynomiale contrasten over dosis op reactietijd (multivariaat):
  L (lineair):    F(1, 49) =  6.20, p = .016
  Q (kwadratisch): F(1, 49) =  4.10, p = .048
  C (kubisch):    F(1, 49) = 12.50, p < .001

Welke uitspraak past het best?

Het lineaire contrast op de reactietijd is het belangrijkst omdat het altijd voor het kwadratische en kubische uit gaat.
Het kubische contrast is significant en het sterkst — het effect van dosis op de reactietijd vertoont een S-vormig of op-neer-op-neer-patroon dat niet door alleen lineair of kwadratisch gevangen wordt.
Drie significante contrasten zijn onmogelijk; ze moeten orthogonaal zijn dus maar één kan tegelijk significant zijn.
Kubisch is altijd ruis bij vier niveaus en wordt genegeerd.

Antwoord E2 — open na je eigen poging

b) Kubisch significant + sterkst betekent dat het effect van dosis op de reactietijd een S-vorm of zigzag-vorm heeft (op-neer-op-neer of neer-op-neer-op). Lineair vangt een monoton effect van dosis op de reactietijd, kwadratisch een U-vorm; kubisch is wat daarbovenop nog over is. Optie a is onjuist (geen orde-prioriteit; alle drie staan los). Optie c verwart orthogonaliteit met significantie (orthogonaal betekent dat ze onafhankelijke stukken meten, daardoor kunnen ze onafhankelijk significant of niet zijn). Optie d is fout — een significant kubisch patroon is geen ruis, maar een echte vorm-bevinding.

Vraag E3 — Sphericity en GG-correctie (univariate route)

Bij de univariate-route RMA wordt sphericity als aanname getest met Mauchly’s toets:

Mauchly's Test for Sphericity:
W = 0.78, df = 5, p = .021

If sphericity is violated:
  Greenhouse-Geisser ε = 0.83
  Huynh-Feldt ε        = 0.87

Welke uitspraak is correct?

Sphericity is geschonden ($p < .05$); rapporteer de Greenhouse-Geisser-gecorrigeerde univariate $F$, of gebruik de multivariate route die geen sphericity-aanname maakt.
Mauchly’s $p = .021$ betekent dat we de univariate route niet meer mogen gebruiken; alleen Bayesiaanse RMA is veilig.
De Greenhouse-Geisser $\epsilon = 0.83$ betekent dat $83\%$ van de variantie verklaard wordt.
Huynh-Feldt is altijd strenger dan Greenhouse-Geisser.

Antwoord E3 — open na je eigen poging

a) Mauchly significant ($p = .021 < .05$) → sphericity-aanname geschonden bij de univariate route. Twee oplossingen: (1) GG- of HF-correctie op $df$ van de univariate $F$; (2) de multivariate route uit dit thema, die geen sphericity-aanname kent. Optie b is overdreven (Bayes is een optie, maar niet de enige). Optie c verwart $\epsilon$ met variantie-verklaring ($\epsilon$ is een schaalfactor voor $df$, geen $R^2$). Optie d is precies omgekeerd: GG is conservatiever, HF is iets minder conservatief.

Vraag E4 — Mixed-design hoofdeffect plus interactie

Een onderzoeker rapporteert een mixed-design RMA met habitat (bos/woestijn) en moment (drie momenten) op de mate-van-zichzelf:

Multivariate Tests:
  habitat            F(1, 78) =  8.50, p = .005
  moment             F(2, 77) = 31.20, p < .001
  habitat:moment     F(2, 77) =  7.80, p = .001

Welke conclusie is het meest correct?

Het hoofdeffect van habitat op de mate-van-zichzelf, het hoofdeffect van moment op de mate-van-zichzelf en hun interactie zijn alle drie significant; rapporteer naast elkaar zonder verdere stappen.
De interactie habitat $\times$ moment op de mate-van-zichzelf is significant, dus het hoofdeffect van moment mag niet los geïnterpreteerd worden — splits per habitat-niveau en bekijk binnen elk niveau het effect van moment op de mate-van-zichzelf afzonderlijk.
De interactie van $F = 7.80$ is kleiner dan het hoofdeffect van moment $F = 31.20$, dus de interactie kan worden genegeerd.
Bij significante interactie moet de hele MANOVA opnieuw zonder een van de factoren.

Antwoord E4 — open na je eigen poging

b) Qualified by interaction-conventie: bij significante interactie habitat $\times$ moment op de mate-van-zichzelf zegt het hoofdeffect van moment alleen iets over het gemiddelde effect van moment over habitats heen — dat gemiddelde is misleidend als bos- en woestijn-kameleons verschillende moment-patronen in de mate-van-zichzelf hebben. Vervolgstap: split-by-habitat-RMA. Optie a is incompleet (mist de interactie-interpretatie). Optie c is fout (relatieve grootte van $F$’s telt niet voor interpretatie-prioriteit; significantie wel). Optie d is overdreven (factor-eliminatie is geen gangbare reactie).

R-practical opdrachtje

De voorzichtigheid van de muis

Vraag E5 — Mini-RMA bij de muis

Een ecoloog onderzoekt $70$ muizen op drie momenten rond een gebeurtenis: in rust, bij een nieuwe geur (alarm), en na herstel. Per moment is de voorzichtigheid (continu, $0$-$100$) gemeten. De dataset staat in data/muis_voorzichtigheid.RData en bevat het object muis_voorzichtigheid met kolommen id, rust, nieuwe_geur, na_herstel. Sla je R-commando’s op in één scriptbestand: muis.R. Gebruik $\alpha = .05$ voor de multivariate toets en Bonferroni-correctie voor paarsgewijze follow-up.

a) Verken de data: gemiddelden, SD’s, correlaties tussen momenten. Past het patroon bij een psychologisch litteken (rust laag, geur hoog, herstel ergens daartussenin)?

b) Fit een eenwegs-RMA via de multivariate route met lm(cbind(rust, nieuwe_geur, na_herstel) ~ 1) en Anova(., idata, idesign, type = 3). Rapporteer Pillai’s $V$, $F$-approximatie, $df$ en $p$.

c) Voer polynomiale contrasten uit (lineair en kwadratisch). Welk patroon past het sterkst — monotoon stijgend / dalend (lineair) of dip-of-piek (kwadratisch)?

d) Voer alle drie paarsgewijze paired $t$-toetsen uit met Bonferroni-correctie. Welke paren blijven significant? Schrijf één APA-zin met je conclusie.

Antwoord E5 — open na je eigen poging

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

'data.frame':   70 obs. of  4 variables:
 $ id         : Factor w/ 70 levels "m01","m02","m03",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ rust       : num  46.3 54.9 23.7 23.9 16.6 34.7 30.1 22.2 30.4 53.4 ...
 $ nieuwe_geur: num  93.4 91.3 82.9 79.4 71.5 70.8 50.6 53.6 47.1 60.8 ...
 $ na_herstel : num  56.4 72.8 53 53.4 24.1 53.4 39.9 20.9 30.4 44.4 ...

round(colMeans(muis_voorzichtigheid[, c("rust", "nieuwe_geur", "na_herstel")]), 2)

       rust nieuwe_geur  na_herstel 
      29.00       64.41       42.14

round(apply(muis_voorzichtigheid[, c("rust", "nieuwe_geur", "na_herstel")], 2, sd), 2)

       rust nieuwe_geur  na_herstel 
      12.30       14.45       12.40

round(cor(muis_voorzichtigheid[, c("rust", "nieuwe_geur", "na_herstel")]), 2)

            rust nieuwe_geur na_herstel
rust        1.00        0.43       0.60
nieuwe_geur 0.43        1.00       0.46
na_herstel  0.60        0.46       1.00

# Defensief: contrasten ook hier expliciet zetten vlak vóór de Type-III-aanroep.
options(contrasts = c("contr.sum", "contr.poly"))

lm_m <- lm(cbind(rust, nieuwe_geur, na_herstel) ~ 1,
           data = muis_voorzichtigheid)

idata_m <- data.frame(
  moment = factor(c("rust", "nieuwe_geur", "na_herstel"),
                  levels = c("rust", "nieuwe_geur", "na_herstel"))
)

res_m <- Anova(lm_m, idata = idata_m, idesign = ~moment, type = 3)
summary(res_m, multivariate = TRUE)


Type III Repeated Measures MANOVA Tests:

------------------------------------------
 
Term: (Intercept) 

 Response transformation matrix:
            (Intercept)
rust                  1
nieuwe_geur           1
na_herstel            1

Sum of squares and products for the hypothesis:
            (Intercept)
(Intercept)     1286139

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1   0.94847 1270.029      1     69 < 2.22e-16 ***
Wilks             1   0.05153 1270.029      1     69 < 2.22e-16 ***
Hotelling-Lawley  1  18.40622 1270.029      1     69 < 2.22e-16 ***
Roy               1  18.40622 1270.029      1     69 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: moment 

 Response transformation matrix:
            moment1 moment2
rust              1       0
nieuwe_geur       0       1
na_herstel       -1      -1

Sum of squares and products for the hypothesis:
          moment1   moment2
moment1  12096.69 -20484.97
moment2 -20484.97  34689.98

Multivariate Tests: moment
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1  0.860216 209.2332      2     68 < 2.22e-16 ***
Wilks             1  0.139784 209.2332      2     68 < 2.22e-16 ***
Hotelling-Lawley  1  6.153917 209.2332      2     68 < 2.22e-16 ***
Roy               1  6.153917 209.2332      2     68 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

            Sum Sq num Df Error SS den Df F value    Pr(>F)    
(Intercept) 428713      1    23292     69 1270.03 < 2.2e-16 ***
moment       44848      2    12158    138  254.52 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Mauchly Tests for Sphericity

       Test statistic  p-value
moment        0.90122 0.029126


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

       GG eps Pr(>F[GG])    
moment 0.9101  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

          HF eps  Pr(>F[HF])
moment 0.9334227 4.44638e-44

contrasts(idata_m$moment) <- contr.poly(3)
res_m_poly <- Anova(lm_m, idata = idata_m, idesign = ~moment, type = 3)
summary(res_m_poly, multivariate = TRUE)


Type III Repeated Measures MANOVA Tests:

------------------------------------------
 
Term: (Intercept) 

 Response transformation matrix:
            (Intercept)
rust                  1
nieuwe_geur           1
na_herstel            1

Sum of squares and products for the hypothesis:
            (Intercept)
(Intercept)     1286139

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1   0.94847 1270.029      1     69 < 2.22e-16 ***
Wilks             1   0.05153 1270.029      1     69 < 2.22e-16 ***
Hotelling-Lawley  1  18.40622 1270.029      1     69 < 2.22e-16 ***
Roy               1  18.40622 1270.029      1     69 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: moment 

 Response transformation matrix:
                 moment.L   moment.Q
rust        -7.071068e-01  0.4082483
nieuwe_geur -9.073800e-17 -0.8164966
na_herstel   7.071068e-01  0.4082483

Sum of squares and products for the hypothesis:
           moment.L  moment.Q
moment.L   6048.343 -15319.01
moment.Q -15319.014  38799.41

Multivariate Tests: moment
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1  0.860216 209.2332      2     68 < 2.22e-16 ***
Wilks             1  0.139784 209.2332      2     68 < 2.22e-16 ***
Hotelling-Lawley  1  6.153917 209.2332      2     68 < 2.22e-16 ***
Roy               1  6.153917 209.2332      2     68 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

            Sum Sq num Df Error SS den Df F value    Pr(>F)    
(Intercept) 428713      1    23292     69 1270.03 < 2.2e-16 ***
moment       44848      2    12158    138  254.52 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Mauchly Tests for Sphericity

       Test statistic  p-value
moment        0.90122 0.029126


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

       GG eps Pr(>F[GG])    
moment 0.9101  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

          HF eps  Pr(>F[HF])
moment 0.9334227 4.44638e-44

ps_m <- c(
  rust_geur     = t.test(muis_voorzichtigheid$rust,
                         muis_voorzichtigheid$nieuwe_geur, paired = TRUE)$p.value,
  rust_herstel  = t.test(muis_voorzichtigheid$rust,
                         muis_voorzichtigheid$na_herstel,  paired = TRUE)$p.value,
  geur_herstel  = t.test(muis_voorzichtigheid$nieuwe_geur,
                         muis_voorzichtigheid$na_herstel,  paired = TRUE)$p.value
)
round(ps_m, 6)

   rust_geur rust_herstel geur_herstel 
           0            0            0

round(p.adjust(ps_m, method = "bonferroni"), 6)

   rust_geur rust_herstel geur_herstel 
           0            0            0

a) $M_{\text{rust}} = 29.00$ ($SD = 12.30$), $M_{\text{nieuwe\_geur}} = 64.41$ ($SD = 14.45$), $M_{\text{na\_herstel}} = 42.14$ ($SD = 12.40$). De geur tilt voorzichtigheid fors omhoog; na herstel zakt het weer, maar niet helemaal terug naar rust — een psychologisch litteken van ongeveer $13$ punten. Within-subject correlaties matig positief ($r = .43$-$.60$).

b) Pillai’s $V = .86$, $F(2, 68) = 209.23$, $p < .001$ — een zeer sterk multivariaat effect.

c) Kwadratische trend domineert: $t(69) = -18.31$, $p < .001$ — piek in het midden (geur). Lineaire trend ook significant maar veel kleiner: $t(69) = 10.00$, $p < .001$ — een lichte stijging van rust naar herstel, het litteken-effect.

d) Alle drie paren significant na Bonferroni ($\alpha = .05/3 = .017$): rust vs nieuwe geur ($t(69) = -20.53$), rust vs na herstel ($t(69) = -10.00$), nieuwe geur vs na herstel ($t(69) = 13.19$), alle $p < .001$.

Een eenwegs-RMA toonde een sterk effect van moment op voorzichtigheid bij de muis, Pillai’s $V = .86$, $F(2, 68) = 209.23$, $p < .001$. Polynomiale contrasten lieten een sterk kwadratisch patroon zien ($t(69) = -18.31$, $p < .001$): hoge voorzichtigheid bij de nieuwe geur, lager in rust en na herstel. Paarsgewijze paired $t$-toetsen met Bonferroni-correctie ($\alpha = .017$) bevestigden dat alle drie momenten onderling verschilden (alle $p < .001$); voorzichtigheid keerde na herstel niet volledig terug naar het rust-niveau, een effect van $M = 13.14$ punten.

Wat blijft liggen

Onderwerpen die buiten dit thema blijven

RMA via de multivariate route met polynomiale contrasten en Bonferroni-paren is de hoofdingang van herhaalde-metingen-analyse, maar er ligt veel achter. Wat dit werkboek je gaf is genoeg voor de gangbare onderzoeks-rapportage; voor verdere studie:

Sphericity-correcties (Greenhouse-Geisser, Huynh-Feldt) — niet behandeld in dit werkboek omdat we de multivariate route prefereren. Bij kleine $n$ is de univariate route met GG- of HF-correctie krachtiger; check Mauchly’s toets en pas de correctie toe op de univariate $F$.
Mixed-effects models uitgewerkt — random intercepts plus random slopes, kruisgeneste data, niet-balanced designs. Het pakket lme4 met lmer() is de toegangspoort. Sluit aan bij thema’s over multilevel-analyse.
Robuuste RMA — bij niet-normale of outlier-rijke trajecten: getrimde gemiddelden of M-schattingen via Wilcox’ robuuste varianten in WRS2. Niet behandeld.
Bayesiaanse RMA — multivariate posterior-verdelingen via brms::brm(mvbind(Y1, Y2, Y3, Y4) ~ 1 + (1|id)). Toegankelijke ingang als je verder wilt.
Time-series-analyse voor langere reeksen ($k \gg 5$) — autocorrelatie-modellen (AR, MA, ARMA), niet RMA. Een ander hoofdstuk uit de statistiek.
Latent growth curve modeling — herhaalde metingen als latente variabele in een SEM-frame. Pakket lavaan. Verwijst naar de SEM-cursus.
Missing data in RMA — complete-case-deletie versus multiple imputation (mice, Amelia). Mixed-effects-modellen hanteren missing automatisch beter dan de standaard RMA-route.

Voor verdieping: Maxwell, Delaney & Kelley (2018, Designing Experiments and Analyzing Data); Tabachnick & Fidell (2019, Using Multivariate Statistics); Bates, Mächler, Bolker & Walker (2015, Fitting Linear Mixed-Effects Models Using lme4, JSS).

Aan het eind van de dag

Toen de zon onderging boven de tak, klapte de kameleon zijn schrift dicht.

“De vier momenten zijn vier gezichten,” zei hij. “Mijn ochtend ben ik bedacht, mijn middag onbestemd, mijn avond herstellend, mijn nacht volledig. Niet één kameleon, maar een loop.”

“En wat trekt eraan?” vroeg de pad, die langs kwam.

“Een dip in de middag — de zon staat hoog, ik los op. En een piek in de nacht — de wereld is stil, ik ben er weer.” Hij wachtte even. “Pillai zegt: $.88$. Het is bijna alles wat ik ben.”

“En lineair?” vroeg de pad.

“Stijgend,” zei de kameleon. “Maar niet glad. Met een dip die alleen kwadratisch ziet, en een hapering die alleen kubisch ziet. Drie trends, één dier, vier momenten — en zes paren die elk significant verschillen.”

“En morgen?”

“Morgen meet ik weer,” zei de kameleon. “Misschien ben ik dan een ander dier. Of dezelfde kameleon, op vier andere momenten.”

De pad knikte en hopte verder. Boven de tak werd het donker, en de kameleon liet zijn kleur uitvloeien tot wat de schemer ervan maakte.

Verantwoording

Dit werkboek is geschreven voor studenten die repeated measures ANOVA 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 Maxwell, Delaney & Kelley (2018), Tabachnick & Fidell (2019) en Bates et al. (2015) volgen de gebruikelijke citaatpraktijk in dit veld.

Versie: July 2026 — CountCamp Lab

#	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

Bron	Formule	Hier
Within-subjects (binnen-persoon) totaal	\(df_W = (n-1)(k-1)\)	\(89 \cdot 3 = 267\)
Conditie (moment)	\(df_{cond} = k - 1\)	\(3\)
Residu within	\(df_{resW} = (n-1)(k-1)\)	\(267\)
Multivariate \(df_H\)	\(k - 1\)	\(3\)
Multivariate \(df_E\)	\(n - k + 1\)	\(87\)

Trend	Ochtend	Middag	Avond	Nacht
Lineair (L)	\(-3/\sqrt{20}\)	\(-1/\sqrt{20}\)	\(1/\sqrt{20}\)	\(3/\sqrt{20}\)
Kwadratisch (Q)	\(1/2\)	\(-1/2\)	\(-1/2\)	\(1/2\)
Kubisch (C)	\(-1/\sqrt{20}\)	\(3/\sqrt{20}\)	\(-3/\sqrt{20}\)	\(1/\sqrt{20}\)

\(k\)	\(\binom{k}{2}\)	Bonferroni-drempel (\(\alpha = .05\))
2	1	\(.050\)
3	3	\(.0167\)
4	6	\(.0083\)
5	10	\(.0050\)
6	15	\(.0033\)

Het verhaal van dit thema

Werkmaterialen — R-pakketten en functies

Notatie — symbolen voor RMA

6.0 Project- en datavoorbereiding

6.1 Het RMA-model fitten en interpreteren

T1 — Techniek-keuze: RMA, eenwegs-ANOVA, paired t-toets of mixed-effects?

T2 — Waarom RMA krachtiger dan one-way ANOVA?

6.1.a Verkenning — gemiddelden, correlaties en spaghetti

Means-tabel — moment-overzicht met gt

T3 — Multivariate vs univariate route, en de sphericity-aanname

6.1.b RMA fitten — Leiden’s route via factorial_design()

6.1.c De vier teststatistics — Pillai als primary

T4 — Polynomiale contrasten in RMA

6.1.d Multivariate effectgrootte

6.2 Contrasten en follow-up

6.2.0 Contrasten met de hand opbouwen — rustige bijles-volgorde

6.2.a Polynomiale contrasten — eerst hand, dan R

In R — met icontrasts = "contr.poly"

T5 — Geplande contrasten en orthogonaliteit

6.2.b Paarsgewijze vergelijking met Bonferroni

T6 — Bonferroni voor RMA-paren handrekenen

6.2.c Visualisatie — gemiddeld profiel + spaghetti

6.3 Polynomial-trends-surprise — vier dagdeel-profielen

De vier profielen

Visualisatie — vier profielen in één plot

6.3.a Reken de drie contrasten uit

T7 — RMA-aannames

T8 — Missing data in RMA

T9 — RMA versus mixed-effects models

6.A Uitbreiding — mixed-design (between \(\times\) within)

T10 — Between \(\times\) within: wanneer interactie interessant?

R-spiekblad bij RMA

Pakketten en data laden — wide en long format

Verkenning — spaghetti, correlaties, boxplot

RMA fitten — Leiden’s factorial_design() + Anova()

De vier teststatistics — Pillai, Wilks, Hotelling-Lawley, Roy

Polynomiale contrasten — contr.poly()

Geplande contrasten + Bonferroni voor paren

Mixed-effects alternatief (kort, voor wie verder wil)

Visualisatie — spaghetti + gemiddelde lijn

Voorbeeld-tentamenvragen

Theorie en handreken

R-practical opdrachtje

Wat blijft liggen

Aan het eind van de dag

Verantwoording

Means-tabel — moment-overzicht met `gt`

6.1.b RMA fitten — Leiden’s route via `factorial_design()`

In R — met `icontrasts = "contr.poly"`

RMA fitten — Leiden’s `factorial_design()` + `Anova()`

Polynomiale contrasten — `contr.poly()`