5. Multivariate Variantieanalyse (MANOVA)

Eindelijk mogen we in zinvolle optellingen denken

Het verhaal van dit thema

De leeuw lag op een rotsplateau en keek over zijn troep.

“Een leeuw voelt nooit één ding tegelijk,” zei hij. “Trots, jaloezie, eenzaamheid, woede — die lopen door elkaar. Wie ze los toetst, mist de leeuw.”

Honderdtien leeuwen, vier gevoelens per dier. Per leeuw noteerde hij trots, jaloezie, eenzaamheid en woede (elk op een $0$-$100$ schaal), plus zijn sociale_status in de troep (alpha, middenrang, outsider). Vraag: verschilt het gevoels-profiel tussen statusgroepen — niet één gevoel apart, maar de vier samen?

De techniek heet multivariate analysis of variance MANOVA.

Hehe, eindelijk multivariaat.

Vier weken lang werkten we met één $Y$ tegelijk. MRA: één voorspelde uitkomst. ANOVA: één gemiddelde over groepen. ANCOVA: één DV met een covariaat. LRA: één binaire $Y$. Voor elke vraag een keurig opgeschoonde projectie naar één afhankelijke variabele — alsof de wereld zich zo netjes in losse $Y$’s laat ontleden.

Maar één leeuw voelt niet drie keer apart-iets. Een mens met depressie heeft slaap-, eet-, energie- én moed-problemen tegelijk, en die hangen samen. Een opleiding met goede docenten heeft vaak ook goede roosters en goede gebouwen. Honderden variabelen zijn er in toegepast werk — en veel ervan zeggen ongeveer hetzelfde. Multivariate analyse begint precies waar je dat erkent: je gooit niet je honderd metingen los over de tafel, je reduceert ze tot een paar zinvolle gewogen optellingen (de $D$’s die nog komen) die het hele profiel vangen. Vanaf hier mag het.

Bij MANOVA werk je in vier stappen — in het Engels heet de route ook wel protected F: drie stappen, met een voorbereidende stap ervoor:

0. Klopt het allemaal wel? (preliminary) — multivariate normaliteit, gelijke covariantie-matrices (Box’s M), groepsgrootte-balans. Eerst checken, dan rekenen.
1. Voorspelt het samen iets? — de multivariate $F$-toets (Pillai’s trace en varianten) op het gezamenlijke effect van de factor over alle DV’s.
2. Welke DV’s dragen bij? — univariate $F$ per DV, met Bonferroni- of Holm-Bonferroni-correctie op $\alpha$.
3. Tussen welke groepen op die DV? — pairwise comparisons (Tukey HSD), bij $I$ groepen zijn dat $I(I-1)/2$ paren. Voor onze drie statusgroepen: drie paren (alpha–mid, alpha–outsider, mid–outsider).

TipWat je na dit thema moet kunnen

Na dit thema kun je:

1. herkennen wanneer MANOVA de juiste techniek is: groepen die je op twee of meer afhankelijke variabelen tegelijk vergelijkt;
2. de drie aannames kennen — gelijke covariantiematrices, multivariate normaliteit en onafhankelijke metingen — en de eerste toetsen met Box’s $M$, bij $\alpha = .001$;
3. met de vuistregels ($n \geq 20$ per groep, een gebalanceerd design) beoordelen of de MANOVA robuust is als een aanname niet klopt;
4. de uitvoer van een MANOVA lezen: de vier multivariate toetsen (Pillai, Wilks, Hotelling-Lawley en Roy), elk met zijn $F$, vrijheidsgraden en $p$;
5. beide manieren kennen om een significante MANOVA uit te werken:
   • per afhankelijke variabele een univariate $F$, met de Bonferroni- of de Holm-Bonferroni-correctie, en daarna Tukey binnen die variabele;
   • een discriminantanalyse;
6. uitleggen wat er schort aan de eerste manier, de protected F-route. Je volgt hem wel, maar de naam belooft te veel: de MANOVA vooraf beschermt de losse toetsen maar half tegen een vals alarm, en de losse toetsen kijken niet naar de samenhang tussen de afhankelijke variabelen (het antwoord bij T2 laat zien hoe dat zit);
7. de uitvoer van een discriminantanalyse lezen:
   • hoeveel functies er hooguit zijn: $\min(I - 1, p)$, met $I$ groepen en $p$ afhankelijke variabelen;
   • per functie de eigenvalue en het percentage verklaarde groepsvariantie;
   • wat een functie betekent, uit de structure matrix: een afhankelijke variabele telt mee als zijn correlatie met de functie boven $.40$ of onder $-.40$ ligt;
   • waar elke groep op een functie ligt, afgelezen uit de uitvoer;
8. er een conclusie bij schrijven, in gewone woorden en in APA-stijl.

In R doe je dat met options(contrasts = …), lm(cbind(...) ~ factor), Manova(., type = 3), summary(., multivariate = TRUE, univariate = FALSE), box_m(), summary.aov(), p.adjust() en TukeyHSD(aov(...)).

Wat dichtgeklapt staat onder Voor wie verder wil, hoort niet bij dit lijstje. De dichtgeklapte antwoorden horen er wel bij.

Twee dingen horen ook bij discriminantanalyse, al oefent dit hoofdstuk ze niet:

• waar een groep op een discriminantfunctie ligt, uitrekenen: je vult de groepsgemiddelden in de gewogen optelling in, met de gewichten die de opgave geeft, en soms nog een vast getal erbij, de constante. De tussenstap met de gewogen optellingen, verderop, laat het idee zien met zelfgekozen gewichten;
• uit een classificatietabel aflezen hoe goed een discriminantanalyse de gevallen in hun eigen groep indeelt: welk deel goed ingedeeld is, en bij welke groep het best.

Factoriële MANOVA, met twee factoren tegelijk, staat in dit hoofdstuk alleen in algemene vorm, zonder eigen data.

Wat blijft liggen, onderaan, noemt wat buiten dit thema blijft.

OpmerkingA. Wat is MANOVA?

MANOVA test of meerdere uitkomsten tegelijk verschillen tussen groepen — multivariate ANOVA. Tot nu toe was er één afhankelijke variabele $Y$; bij MANOVA zijn er meerdere ($Y_1, Y_2, \ldots, Y_p$) die je gezamenlijk vergelijkt.

Waarom multivariate? (1) Sommige constructen zijn meerdimensionaal (stemming = trots + jaloezie + eenzaamheid + woede; intelligentie = verbal + spatial + numeric). (2) Correlaties tussen DV’s bevatten zelf informatie. (3) Eén multivariate toets stelt de gezamenlijke vraag in één keer, met één kans van $.05$ op een vals alarm; vier losse ANOVA’s stellen vier vragen, en dan loopt de kans op minstens één vals alarm op — de family-wise error rate FWER. De losse toetsen die je ná een significante MANOVA doet, beschermt de MANOVA maar half; dat werk doet de Bonferroni-correctie (het antwoord bij T2 laat zien hoe dat zit).

De H₀ in twee niveaus

Het slimme aan MANOVA is dat de $H_0$ eigenlijk twee stellingen tegelijk toetst.

Stelling A — per DV zijn de groepsgemiddelden gelijk. Per variabele $j$ schrijf je dat uit als één rijtje:

$$\mathit{\mu}_{j,1} = \mathit{\mu}_{j,2} = \ldots = \mathit{\mu}_{j,I}$$

waar $\mathit{\mu}_{j,g}$ het gemiddelde is van variabele $j$ in groep $g$ (afspraak in dit oefenboek: eerste index = positie-binnen-set, tweede index = groep). Bij $p$ DV’s krijg je $p$ van die rijtjes, één per variabele. Dit is wat $p$ losse univariate ANOVA’s afzonderlijk toetsen.

Kleine noot: in informele uitleg zie je dit ook zonder subscripts (“alle gemiddelden gelijk”). Wij houden ze er altijd bij — kost één teken, levert helderheid op.

Stelling B — voor élke manier waarop je de DV’s met gewichten optelt, is de gemiddelde uitkomst over groepen gelijk. Concreter: kies een willekeurige gewichten-vector, bijvoorbeeld $\mathbf{w} = (3, 1, -2, 0)$. Bereken voor élke leeuw $D = 3 \cdot \text{trots} + 1 \cdot \text{jaloezie} - 2 \cdot \text{eenzaamheid} + 0 \cdot \text{woede}$ — één $D$-score per leeuw. Dan moet onder stelling B gelden: $\bar{D}_{\text{alpha}} = \bar{D}_{\text{middenrang}} = \bar{D}_{\text{outsider}}$. En dat moet gelden voor elke keuze van gewichten, niet alleen deze.

De gezamenlijke $H_0$ is: stelling A én stelling B zijn beide waar.

Als H₀ niet waar is

In de populatie staan of vallen stelling A en stelling B samen. Verschilt er één DV tussen de groepen, dan verschilt ook de gewogen optelling die alleen naar die DV kijkt — gewicht $1$ op die DV, $0$ op de rest. Dan is stelling B dus ook niet waar. En andersom: verschilt er een gewogen optelling tussen de groepen, dan verschilt ook minstens één DV zélf tussen de groepen. Alleen A niet waar of alleen B niet waar bestaat dus niet. Waarom dan toch twee stellingen? Omdat ze laten zien waar elke toets kijkt: losse ANOVA’s kijken alleen naar stelling A, MANOVA zoekt ook alle gewogen optellingen van stelling B af.

Het verschil tussen MANOVA en $p$ losse ANOVA’s zit niet in wát ze toetsen, maar in wat ze in een steekproef kunnen zien:

1. Eén of meer DV’s verschillen duidelijk tussen de groepen. Dan vinden de losse ANOVA’s het, en de multivariate toets ook.
2. Elke DV verschilt een beetje tussen de groepen — te weinig om in zijn eentje op te vallen. Geen enkele losse ANOVA is significant, en tóch kan een gewogen optelling de groepen duidelijk scheiden. Zo’n optelling telt de kleine verschillen bij elkaar op. En hangen de DV’s binnen de groepen samen, dan kan een goed gekozen combinatie bovendien een deel van de ruis wegstrepen: die correlatie-structuur bepaalt hoeveel ruis er in een combinatie overblijft. Dit geval ziet alleen MANOVA.

Belangrijk: dat tweede geval is precies waarom MANOVA niet hetzelfde is als $p$ losse ANOVA’s. Bij de tulpen, aan het eind van de follow-up, zie je twee oogsten met precies dezelfde groepsgemiddelden en precies dezelfde losse ANOVA’s — en tóch is de MANOVA in de ene oogst significant en in de andere niet. Het verschil zit in hoe de maten binnen de groepen samenhangen.

OpmerkingVoor wie verder wil — dezelfde $H_0$ in vectoren

Hierboven staat de $H_0$ uitgeschreven per variabele. In tekstboeken (Stevens, Tabachnick & Fidell, Field) zie je een compactere notatie met vectoren:

$$H_0: \boldsymbol{\mu}_1 = \boldsymbol{\mu}_2 = \ldots = \boldsymbol{\mu}_I$$

waar $\boldsymbol{\mu}_g$ de vector van $p$ DV-gemiddelden in groep $g$ is — alle rijtjes van stelling A samengevat in één regel. Stelling B krijgt dan:

$$H_0: \mathbf{a}^\top \boldsymbol{\mu}_1 = \mathbf{a}^\top \boldsymbol{\mu}_2 = \ldots = \mathbf{a}^\top \boldsymbol{\mu}_I \quad \text{voor élke gewichten-vector } \mathbf{a}$$

Dezelfde stelling, andere schrijfwijze. Handig om te herkennen in tijdschriftartikelen en methodeboeken. In de hoofdtekst houden we het gezellig in losse rijtjes.

OpmerkingB. De vier teststatistics

Bij MANOVA bestaan vier multivariate teststatistieken. Allemaal worden ze omgezet naar een $F$-approximatie waarop je een $p$-waarde leest. Je rapporteert ze alle vier: Pillai’s trace eerst, omdat die het robuustst is als een aanname niet helemaal klopt, en de andere drie ernaast, zodat je lezer ziet of ze het eens zijn.

Statistic	Symbool	Range	Richting	Voorkeur bij
Wilks’ Lambda	$\Lambda$	$0$ tot $1$	lager = sterker effect (tegenovergesteld aan $F$)	klassieke conventie; gelijke covariantie-matrices
Pillai’s trace	$V$	$0$ tot $\min(I-1, p)$	hoger = sterker	moderne standaard, meest robuust tegen schendingen
Hotelling-Lawley trace	$T^2$	$0$ tot $\infty$	hoger = sterker	hoogste power bij grote effecten + kleine sample
Roy’s largest root	$\theta$	$0$ tot $\infty$	hoger = sterker	kan misleidend zijn (alleen op grootste eigenvalue)

In R: summary(res_Manova, multivariate = TRUE) geeft alle vier. Bij grote effecten zeggen ze meestal hetzelfde; bij borderline kan het verschillen — dan is Pillai het veiligst.

Let op de Wilks-richting: $\Lambda$ is een ratio van error-variantie tot totale variantie. Klein $\Lambda$ = weinig error relatief tot totaal = sterk discriminerend. Tegenovergesteld aan hoe $F$, $V$, $T^2$ en $\theta$ werken (waar hoger = sterker). SPSS-tijd-veteranen wisten dit uit hun hoofd; in R is het soms even schakelen. De $p$-waarde zelf is wél altijd “lager = significanter” — zoals overal.

OpmerkingC. De leeuw voelt nooit één ding tegelijk

Een leeuw kan tegelijk trots zijn op zijn pluk-en-tand, jaloers op een rivaal, eenzaam in zijn rust, en woedend over een grensovertreding. Vier gevoelens, niet onafhankelijk: jaloezie kan woede aanwakkeren, eenzaamheid kan trots ondermijnen, en alpha’s voelen anders dan middenrangers of outsiders. Een aparte ANOVA per gevoel mist die samenhang.

MANOVA pakt het hele gevoels-spectrum tegelijk: zien sociale alpha’s, middenrangers en outsiders andere patronen over deze vier dimensies samen? Niet “verschilt trots”, maar “verschilt het gevoels-profiel als geheel”. Pas als die multivariate vraag bevestigend beantwoord wordt, ga je per DV verder kijken.

OpmerkingVoor wie verder wil — drie statusgroepen, of drie hechtingsstijlen met manen?

Een leeuw is geen mens, maar het is verleidelijk. Sociale alpha lijkt op wat hechtings-theoretici secure noemen: bovenaan de hiërarchie, gerust in eigen kunnen, mist ’m soms zijn moeder maar bekent het niet hardop. De middenrang zit dichter bij anxious-resistant: kijkt voortdurend omhoog naar de alpha en omlaag naar de outsider, hangt en klampt aan rangorde. De outsider doet avoidant: houdt afstand, lijkt onverschillig, voelt het wel.

Voor MANOVA-doeleinden gaat het hier om drie groepen met vier gevoelens als DV. Voor wie nieuwsgierig is naar de echte theorie: Bowlby (1969) is de pionier, Ainsworth (1970s, Strange Situation) deed de empirische uitwerking. Latere lijn loopt via Hazan & Shaver naar het moderne tweedimensionale model — avoidance × anxiety, Brennan-Clark-Shaver ECR. Hier is het een knipoog.

OpmerkingTussenstap — wat doe je met vier scores per leeuw?

Stel je voor: je moet een collega vertellen hoe een leeuw zich voelt. Vier scores per leeuw oplepelen voor élke leeuw is onoverzichtelijk. Bij honderdtien leeuwen wordt het 440 getallen die niemand meer overziet.

Dan maar even terugbrengen tot één getal per leeuw. We pakken even drie van de vier gevoelens — trots, jaloezie en eenzaamheid — om het idee overzichtelijk te houden. Een gewogen optelling: $$D = w_1 \cdot \text{trots} + w_2 \cdot \text{jaloezie} + w_3 \cdot \text{eenzaamheid}$$

De gewichten $(w_1, w_2, w_3)$ kies je zelf, omdat het mag. We gaan twee verschillende keuzes proberen en kijken wat er gebeurt. Eerste keuze: $\mathbf{w}_1 = (4, 1, 2)$. De uitkomst per groep noemen we $D_{1, g}$ — eerste index = gewichten-vector-nummer (dezelfde afspraak als bij $\mathit{\mu}_{j,g}$: positie-binnen-set eerst, groep tweede):

	trots	jaloezie	eenzaamheid	$D_{1, g} = 4\cdot\text{trots} + 1\cdot\text{jaloezie} + 2\cdot\text{eenzaamheid}$
alpha	90	30	15	$D_{1, \text{alpha}} = 4 \cdot 90 + 1 \cdot 30 + 2 \cdot 15 = 420$
midden	50	80	50	$D_{1, \text{midden}} = 4 \cdot 50 + 1 \cdot 80 + 2 \cdot 50 = 380$
outsider	15	40	80	$D_{1, \text{outsider}} = 4 \cdot 15 + 1 \cdot 40 + 2 \cdot 80 = 260$

Onder gewichten-vector $\mathbf{w}_1$ scoort alpha het hoogst, midden tweede, outsider laagst. Logisch: $\mathbf{w}_1$ legt de meeste nadruk op trots (gewicht $4$), en alpha is daar de duidelijke koploper.

Maar we hadden de gewichten ook anders kunnen kiezen. Tweede keuze: $\mathbf{w}_2 = (1, 4, 2)$. Uitkomst per groep is $D_{2, g}$:

	$D_{2, g} = 1\cdot\text{trots} + 4\cdot\text{jaloezie} + 2\cdot\text{eenzaamheid}$
alpha	$D_{2, \text{alpha}} = 1 \cdot 90 + 4 \cdot 30 + 2 \cdot 15 = 240$
midden	$D_{2, \text{midden}} = 1 \cdot 50 + 4 \cdot 80 + 2 \cdot 50 = 470$
outsider	$D_{2, \text{outsider}} = 1 \cdot 15 + 4 \cdot 40 + 2 \cdot 80 = 335$

Onder $\mathbf{w}_2$ scoort opeens midden het hoogst, outsider tweede, alpha laatst. De ranking is omgedraaid. Dezelfde drie groepen, dezelfde gemiddelden, andere gewichten — andere conclusie. Met $\mathbf{w}_2$ ligt de nadruk op jaloezie (gewicht $4$), en daar is midden het sterkst.

Dat is precies waar MANOVA over gaat. Niet één gewichten-vector kiezen, maar alle mogelijke gewichten-vectoren tegelijk doortesten en kijken of er ergens een combinatie is waarin de groepen verschillen. Dat is stelling B uit callout A. Een gewogen optelling waarbij de groepen écht verschillen, heet later — als we hem door de data laten kiezen in plaats van zelf — een discriminant-functie. Daar komen we in 5.2 op terug.

TipSjonge sjonge — al die woorden voor één ding

In dit hoofdstuk en in latere thema’s kom je een rits termen tegen die allemaal hetzelfde gereedschap aanwijzen. Een paar:

Term	Waar je ’m tegenkomt
gewogen optelling	hier, didactische ingang
lineaire combinatie	de wiskundige naam; ook in vector-context
som van producten	zo lees je het bij hand-rekenen
composite score	psychometrie, als items tot een schaal-score worden samengevoegd
schaalscore / itemtotaalscore	toegepast (vragenlijsten)
contrast	speciaal geval: gewichten tellen op tot $0$ — dan meet je een verschil i.p.v. een totaal
discriminant-functie / discriminant-score	in DDA (sectie 5.2), als de data zelf de gewichten kiest
(canonical) variate / functie-variaat	zelfde ding in MANOVA- en kanonische-correlatie-context
projectie	wiskundig: pijl-op-as
weighted sum	Engelse tekstboeken

Onder de motorkap is het allemaal dezelfde operatie: pak het rijtje scores van één leeuw (één persoon, één geval), vermenigvuldig elke score met een eigen gewicht, tel ze op, krijg één getal terug. Dat getal heet gewogen optelling als jij de gewichten verzint, discriminant-functie als MANOVA ze voor je kiest, contrast als je gewichten optellen tot nul, composite score als jij vragenlijst-items optelt tot een schaal. Verschillende contexten, verschillende namen, één gereedschap.

Als jij straks “discriminant-functie” leest of “canonical variate”, denk: “dat is gewoon een gewogen optelling, andere achternaam.” Zodra je dat verbindt, valt een aantal kwartjes tegelijk. Sjonge.

OpmerkingVoor wie verder wil — wat is een vector eigenlijk?

Toen we zeiden “gewichten $(w_1, w_2, w_3)$”, was dat al een vector — een rijtje van getallen onder elkaar geschreven. Compact:

$$\mathbf{w} = \begin{pmatrix} 4 \\ 1 \\ 2 \end{pmatrix}$$

Vet gedrukt symbool $\mathbf{w}$ zegt: dit is geen los getal, dit is een rijtje. Bij twee getallen onder elkaar kun je het tekenen — als pijl uit de oorsprong naar een punt op een assenkruis.

library(ggplot2)

ggplot() +
  geom_segment(aes(x = 0, y = 0, xend = 3, yend = 2),
               arrow = arrow(length = unit(0.35, "cm"), type = "closed"),
               linewidth = 1, colour = "#4477AA") +
  geom_point(aes(x = 3, y = 2), size = 2.5, colour = "#4477AA") +
  annotate("text", x = 3.2, y = 2.25, label = "(3, 2)",
           hjust = 0, size = 4) +
  annotate("text", x = 1.4, y = 1.3, label = "v",
           fontface = "bold.italic", size = 5, colour = "#4477AA") +
  coord_cartesian(xlim = c(0, 4.5), ylim = c(0, 4.5)) +
  scale_x_continuous(breaks = 0:4) +
  scale_y_continuous(breaks = 0:4) +
  labs(x = expression(Y[1]), y = expression(Y[2])) +
  theme_minimal(base_size = 11) +
  theme(panel.grid.minor = element_blank())

De pijl heeft een richting (waar wijst-ie?) en een lengte (hoe ver kom je?) — allebei samengevat in de twee getallen $(3, 2)$. Drie naar rechts, twee omhoog.

Bij twee DV’s kun je de gemiddelden van een groep tekenen als zo’n pijl: van de oorsprong naar het punt $(\bar{Y}_{1g}, \bar{Y}_{2g})$. Voor élke groep een eigen pijl. Wijzen ze allemaal in dezelfde richting en zijn ze even lang? Dan multivariaat geen verschil. Wijken ze af in richting of lengte? Dan zit er ergens een verschil in het gemiddelden-profiel.

Handig om te herkennen wanneer je het in een tekstboek tegenkomt — zoals in de schrijfwijze van stelling B met vectoren, in het kader daarover bij Wat is MANOVA?.

OpmerkingDe boomstructuur — vier-stappen-flow op een compact voorbeeld

Drie groepen ($A$, $B$, $C$), twee DV’s ($Y_1$, $Y_2$). Hoe loopt de MANOVA-procedure door?

flowchart TD
    S0["**Stap 0** — Assumpties checken<br/>Box's *M*, MV-normaliteit, n-balans"]
    S1{"**Stap 1** — Multivariate F<br/>Pillai's trace"}
    Stop1["**STOP** — geen verschil"]
    S2{"**Stap 2** — Univariate F per DV<br/>Bonferroni α = .05/2 = .025"}
    Y1sig{"Y₁ significant?"}
    Y2sig{"Y₂ significant?"}
    S3a["**Stap 3 voor Y₁** — Tukey HSD<br/>3 paren: A-B, A-C, B-C"]
    S3b["**Stap 3 voor Y₂** — Tukey HSD<br/>3 paren: A-B, A-C, B-C"]
    Klaar1["klaar voor Y₁"]
    Klaar2["klaar voor Y₂"]

    S0 --> S1
    S1 -- niet sig --> Stop1
    S1 -- significant --> S2
    S2 --> Y1sig
    S2 --> Y2sig
    Y1sig -- ja --> S3a
    Y1sig -- nee --> Klaar1
    Y2sig -- ja --> S3b
    Y2sig -- nee --> Klaar2

Lees van boven naar beneden. Stap 0 staat los — die check je vóór het rekenen, niet in de inferentie-keten. De inhoudelijke kern is de driehoek $1 \to 2 \to 3$. Bij stap 1 nee zit je klaar zonder enige univariate uitspraak; bij stap 1 ja loopt het via stap 2 (per DV met Bonferroni-correctie) naar stap 3 (alleen op DV’s die in stap 2 nog significant zijn).

OpmerkingVraag T-vooraf — Hoe ziet deze boom eruit bij 4 groepen en 6 DV’s?

Zelfde vier stappen, maar de takken breder. Twee dingen veranderen:

• Bij stap 2: de Bonferroni-grens wordt $.05/6 \approx .0083$ (zes DV’s in plaats van twee).
• Bij stap 3: per significante DV $4 \cdot 3/2 = 6$ pairwise comparisons (in plaats van 3).

Probeer zelf de boom te schetsen — zonder zes aparte takken vol te tekenen mag je “per DV in stap 2: zie subboom” schrijven. Het patroon blijft hetzelfde.

OpmerkingMANOVA en DDA in sets met pijlen — eerst los, dan samen

Een tweede manier om de techniek vast te leggen: sets met pijlen ertussen. Per set een rechthoek (meetniveau bovenin, variabel-namen in vet, aantal onder); pijl met techniek-label tussen de sets.

Eerst MANOVA alleen — wat doet de techniek?

Eén nominale IV op de vier continue DV’s tegelijk — de multivariate $F$ vraagt of er ergens in dat gevoels-profiel een groeps-verschil zit (stelling A of B uit callout A).

Dan DDA alleen — wat doet die in 5.2?

DDA vertrekt vanuit de DV’s en maakt er een kleiner aantal gewogen optellingen van — de discriminant-functies $D_1, D_2, \ldots$. De pijl wijst dus de andere kant op: niet “groepen voorspellen DV’s” maar “DV’s worden samengevouwen tot een paar dimensies die groepen scheiden”. Bij DDA staat de pijl omgekeerd ten opzichte van MANOVA.

Allebei samen — zo passen ze:

MANOVA test óf er ergens een verschil zit in het gevoels-profiel. DDA vouwt het profiel daarna samen tot één of twee dimensies waarop het verschil het meest zichtbaar is. “Waar MANOVA stopt, daar gaat DDA verder” — zoals je het in 5.2 zelf zal ervaren.

Dezelfde set-pijl-stijl gebruiken we straks ook in andere thema’s (MRA, ANCOVA, LRA, RMA, Mediation) — herken het als overzicht-conventie.

OpmerkingHoe dit hoofdstuk leest

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

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

De antwoorden staan dichtgeklapt onder hun vraag: klap ze pas open na je eigen poging. T-vragen vlechten zich tussen de hoofdvragen door om je rekenintuïtie te oefenen.

Werkmaterialen — R-pakketten en functies

De leeuw pakt zijn gereedschap

Gereedschap	Waar het voor is	Pakket
lm(cbind(Y1, Y2, Y3) ~ factor)	Multivariaat lineair model — de basis	base
car::Manova(lm_obj, type = 3)	Type-III MANOVA met alle vier teststatistics	car
summary(res, multivariate = TRUE, univariate = FALSE)	De vier multivariate toetsen	car
rstatix::box_m()	Box’s M-toets op gelijke covariantie-matrices	rstatix
summary.aov(lm_obj)	Univariate ANOVA per DV — hier haal je de univariate $F$’s vandaan, niet uit Manova() (zie de waarschuwing in 5.1.c)	base
p.adjust(p, method = "bonferroni")	Multiple-testing-correctie op $p$-waardes	base
TukeyHSD(aov(Y ~ factor))	Tukey HSD post-hoc binnen één DV	base
candisc::candisc(lm_obj)	Descriptive Discriminant Analysis (DDA)	candisc

OpmerkingWaarom lm() + Manova() in plaats van manova()?

R kent een base-functie manova() die direct werkt, maar standaard Type-I sums of squares gebruikt. car::Manova() op een lm-object met type = 3 geeft Type-III SS — past beter bij factorial designs met onbalans en is consistent met onze ANOVA-aanpak in thema 2 en 3 (zie de contr.sum callout daar). Eerst lm(cbind(...) ~ factor) voor het multivariate lineair model, dán Manova() op dat object.

Notatie — symbolen voor MANOVA

OpmerkingSleutelsymbolen in dit hoofdstuk

In dit oefenboek zie je telkens dezelfde notatie:

• $\Lambda$ — Wilks’ Lambda. Klein = sterk effect.
• $V$ — Pillai’s trace.
• $T^2$ — Hotelling-Lawley trace.
• $\theta$ — Roy’s largest root: de grootste eigenvalue (zie de discriminantanalyse verderop).
• $F$ — APA-italic, hier als $F$-approximatie bij elk van de vier teststatistics.
• $df_H$ en $df_E$ — de twee vrijheidsgraden van de $F$-approximatie, in de uitvoer num Df en den Df. Die rapporteer je: bij onze leeuw voor Pillai $F(8, 210)$. De Df-kolom vooraan in de uitvoer is iets anders: daar staat het aantal groepen min één.
• $p$ — $p$-waarde van de $F$-approximatie.
• $N$ — totale steekproef-grootte.

Statistische symbolen met Latijnse letters staan cursief, ook als ze uit twee letters bestaan ($F$, $t$, $p$, $r$, $N$, $df$, $SD$, $SS$). Griekse letters staan rechtop als ze een toetsgrootheid of een drempel noemen ($\Lambda$, $\theta$, $\eta^2$, $\alpha = .05$), zoals in een APA-zin. Rekent een kleine Griekse letter mee in een formule, dan staat hij cursief ($\mathit{\alpha} / k$); dat geldt ook als hij een index draagt die naar zo’n formule wijst, zoals de populatiegemiddelden $\mathit{\mu}_{j,g}$, de drempels $\mathit{\alpha}_h$ van Holm en de eigenvalues $\mathit{\lambda}_1, \mathit{\lambda}_2$. Staat er naast zo’n letter alleen een getal, dan staat hij weer rechtop: in $\lambda_1 = 0.80$ en $\alpha_1 = .05/5 = .010$ krijgt hij een waarde en rekent hij niet. Hoofdletters blijven altijd rechtop, en namen van technieken en toetsen ook (MANOVA, Pillai).

OpmerkingVoor wie verder wil — de vier toetsgrootheden in matrixtaal

In tekstboeken staan de vier toetsgrootheden in matrixtaal:

• $\boldsymbol{\mu}_g$ — vector van DV-gemiddelden in groep $g$. Bij onze vier DV’s: $\boldsymbol{\mu}_g = (\mathit{\mu}_{1,g}, \mathit{\mu}_{2,g}, \mathit{\mu}_{3,g}, \mathit{\mu}_{4,g})^\top$ — eerste index de DV, tweede de groep, net als bij stelling A.
• $\boldsymbol{H}$ — hypothese-matrix (between-groups sums-of-squares-and-cross-products, SSCP).
• $\boldsymbol{E}$ — error-matrix (within-groups SSCP).
• $\Lambda$ $= |\boldsymbol{E}| / |\boldsymbol{H} + \boldsymbol{E}|$.
• $V$ = som van de eigenvalues van $\boldsymbol{H}(\boldsymbol{H}+\boldsymbol{E})^{-1}$.
• $T^2$ = som van de eigenvalues van $\boldsymbol{HE}^{-1}$.
• $\theta$ = grootste eigenvalue van $\boldsymbol{HE}^{-1}$.

De eigenvalues van $\boldsymbol{HE}^{-1}$ zijn dezelfde getallen als de eigenvalues van de discriminant-functies verderop. De SSCP-matrices zie je ook terug in de uitvoer van summary(Manova(...)), als “Sum of squares and products for the hypothesis”.

Voor inferentie spreek je over de populatie-vector $\boldsymbol{\mu}_g^*$ (asterisk-notatie); voor rapportage over de steekproef-vector $\bar{\boldsymbol{Y}}_g$.

TipVuistregels zijn afspraken, geen wetten

De drempels in dit oefenboek:

• Cel-grootte: $n_i \geq 20$ per groep — dan is de $F$-approximatie robuust tegen lichte normaliteit-schendingen.
• Box’s M: standaard $\alpha$ voor Box’s M is $.001$, niet $.05$. De toets is bekend overgevoelig — bij $\alpha = .05$ verwerpt hij vrijwel altijd. Vandaar deze strenge drempel.
• Welke teststatistic: bij gelijke groepen en moderate samples zijn alle vier ongeveer gelijk; bij twijfel of ongelijke covariantie-matrices is Pillai het veiligst.
• Bonferroni voor follow-up: $\mathit{\alpha}_{\text{per DV}} = .05/p$. Bij onze $p = 4$ DV’s: $.05/4 = .0125$. Strenger dan unprotected ANOVA’s.
• MANOVA $F$-robuustheid: bij gebalanceerde celgroottes en $n_i \geq 20$ is de $F$-approximatie robuust tegen normaliteit- en covariantie-matrix-schendingen. Onbalans + $n_{\max}/n_{\min} > 1.5$ + Box’s M $p < .001$ samen is wel een rood vlaggetje.

Met deze drempels voor groepsgrootte, balans en Box’s $M$ rekent dit oefenboek. Handboeken kiezen soms iets andere getallen; dat is een keuze, geen natuurwet.

OpmerkingCode in dit hoofdstuk — wat moet je kunnen typen?

Code die je zelf moet kunnen typen staat open: options(contrasts = c("contr.sum", "contr.poly")), lm(cbind(Y1, Y2, Y3) ~ factor), Manova(., type = 3), summary(., multivariate = TRUE, univariate = FALSE), box_m(), summary.aov(), p.adjust(), TukeyHSD(aov(...)). Code die alleen een plaatje of een tabel mooi maakt — kleuren, opmaak, een tabel met gt — of die je eerder in dit hoofdstuk al open zag staan, zit achter een knopje “Toon code”. Klap hem open als je nieuwsgierig bent. Achter zo’n knopje zit ook de code van de discriminantanalyse, candisc(): bij die analyse oefen je hier het lezen van de uitvoer, en die uitvoer staat gewoon open.

5.0 Project- en datavoorbereiding

Aantekeningen netjes, paden kloppen

Pak het meegestuurde werkboek uit en dubbelklik in de hoofdmap op mvda_werkboek.Rproj. Open daarna dit hoofdstuk, 05_multivariate_variantieanalyse/05_multivariate_variantieanalyse.qmd. RStudio draait de codeblokken vanuit de map van het hoofdstuk, dus load("data/...") vindt de data in 05_multivariate_variantieanalyse/data/. Typ je liever in een eigen script, zet dan eerst de werkmap goed: setwd("05_multivariate_variantieanalyse"). Niets te installeren als de pakketten (car, tidyverse, rstatix, candisc) op je computer staan; anders één keer:

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

1. Open de hoofd-dataset.

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

'data.frame':   110 obs. of  5 variables:
 $ sociale_status: Factor w/ 3 levels "alpha","midden",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ trots         : num  94 66.9 87.5 72.3 81.7 ...
 $ jaloezie      : num  45.4 49.6 49.3 43.8 42.7 ...
 $ eenzaamheid   : num  23.4 29.8 21.4 23.2 43 ...
 $ woede         : num  40.1 43.5 41.2 36.9 42.5 ...

Per leeuw vijf waarden: sociale_status (factor: alpha / midden / outsider), en vier gevoelens — trots, jaloezie, eenzaamheid en woede (elk continu, schaal $0$-$100$). De vier DV’s vormen samen het gevoels-profiel.

OpmerkingVier DV’s of vier aparte ANOVA’s? Het multivariate alternatief

Je zou vier aparte ANOVA’s kunnen runnen — één voor trots, één voor jaloezie, één voor eenzaamheid, één voor woede. Wat is daar mis mee?

• FWER inflation: bij $\alpha = .05$ per ANOVA en vier onafhankelijke toetsen wordt de kans op minstens één type-I-fout $1 - (1-.05)^4 \approx .185$ — bijna één op de vijf. Bonferroni-correctie (een drempel van $.05/4 = .0125$ per toets) repareert dit, maar verliest power.
• Geen gebruik van DV-correlaties: de vier gevoelens hangen samen; die samenhang draagt informatie. MANOVA gebruikt de gezamenlijke covariantie-structuur.
• Eén multivariate hypothese: “verschilt het gevoels-profiel?” is conceptueel één vraag. Vier aparte ANOVA’s beantwoorden vier vragen.

Vandaar de volgorde: eerst MANOVA (multivariate $F$, één hypothese), pas bij significantie univariate follow-up met Bonferroni — de zogenoemde protected F-route. Als de multivariate $F$ niet significant is, stop je: je hebt geen bewijs voor profiel-verschil, dus losse DV-toetsen zijn dan niet meer geloofwaardig.

5.1 Het MANOVA-model fitten en interpreteren

Vier gevoelens tegelijk, één hypothese

“Trots, jaloezie, eenzaamheid, woede,” mompelde de leeuw. “Als ik ze los toets, mis ik het patroon. Als ik ze samen toets, vraag ik of de troep echt drie gevoels-werelden bewoont.”

OpmerkingVoor je gaat rekenen — zes vragen aan jezelf

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

1. Wie of wat wordt er gemeten? De leeuw heeft $110$ leeuwen genoteerd.

2. Wat wordt er gemeten? Vijf variabelen per leeuw: sociale_status, trots, jaloezie, eenzaamheid, woede.

3. Onafhankelijk of afhankelijk? Onderzoeksvraag: voorspelt sociale_status het gevoels-profiel? Dan zijn trots, jaloezie, eenzaamheid en woede afhankelijk ($Y_1, Y_2, Y_3, Y_4$), en sociale_status onafhankelijk (factor).

4. Welke factor(en) en hoeveel niveaus? Eén factor (sociale_status) met drie niveaus (alpha, middenrang, outsider) — dus eenwegs-MANOVA. Bij twee factoren had het tweewegs-MANOVA geheten (zie de sectie over factoriële MANOVA verderop); bij drie factoren drieweg, enzovoort.

5. Meetniveau van elke variabele? Sociale_status nominaal ($3$ niveaus), vier DV’s interval ($0$-$100$).

6. Welke techniek? Doorloop de beslisboom:

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

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

Dezelfde keuze, anders genoteerd. Lees als predictor(en) (meetniveau) → DV(s) (meetniveau):

#	Predictor(en)	DV(s)	Techniek
	dich |X|	|Y| int	$t$-toets
	nom 3+lvl |X|	|Y| int	éénweg ANOVA
1	int |X X …|	|Y| int	MRA — thema 1
2	nom |X X|	|Y| int	factorial ANOVA — thema 2
3	nom + int |X C|	|Y| int	ANCOVA — thema 3
4	bo(in) |X|	|Y| bin	LRA — thema 4
5	nom |X|	|Y₁ Y₂ …| int	MANOVA — thema 5
6	within-subject momenten	|Y₁ Y₂ Y₃ Y₄| int	RMA — thema 6
7	indirect via $M$	$X \to M \to Y$ pad	Mediation — thema 7

Vier continue $Y$’s, één factor $X$. Het aantal DV’s is wat de techniek bepaalt — de mix van DV-meetniveaus moet wel uniform interval zijn (geen mengsels van interval en nominaal als DV). Antwoord: MANOVA.

T1 — Techniek-keuze: MANOVA, ANOVA, ANCOVA of LRA?

OpmerkingVraag T1 — Pen-en-papier

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

a) Een orka onderzoekt of drie pod-types (familiepod / mannenpod / mengpod) verschillen op vier emotionele dimensies samen: rust, opwinding, zorgzaamheid en agressie — alle continu.

b) Een gezondheidspsycholoog meet of stress-niveau (continu) verschilt tussen vier therapie-types (factor) na correctie voor leeftijd (continu). Bij $200$ cliënten.

c) Een onderwijskundige registreert per student of een tentamen wordt gehaald (ja/nee), met als voorspeller aantal studie-uren-per-week (continu). Bij $300$ studenten.

d) Een uil vergelijkt drie jachtperioden (factor, 3 niveaus) op het aantal succesvolle vluchten per nacht (continu).

OpgeletAntwoord T1 — open na je eigen poging

a) Vier continue DV’s (rust, opwinding, zorgzaamheid, agressie), één factor (pod-type, $3$ niveaus). $\Rightarrow$ MANOVA.

b) Eén continue $Y$ (stress-niveau), één factor + één continue covariaat. $\Rightarrow$ ANCOVA (thema 3).

c) Eén binair $Y$ (tentamen-gehaald ja/nee), één continue $X$. $\Rightarrow$ logistische regressie (thema 4).

d) Eén continue $Y$ (succes-aantal), één factor (jachtperiode, NOM, $3$ niveaus). $\Rightarrow$ eenweg-ANOVA (thema 2).

Vuistregel. Aantal DV’s bepaalt of het MANOVA wordt. Eén continue DV + factor(en) + (eventueel) covariaat = ANOVA-familie. Meerdere continue DV’s tegelijk + factor = MANOVA. Binair $Y$ blijft LRA, ongeacht aantal voorspellers.

T2 — Waarom multivariaat in plaats van $p$ aparte ANOVA’s?

OpmerkingVraag T2 — Pen-en-papier

Een onderzoeker heeft drie continue DV’s en één factor met drie niveaus. Hij overweegt drie aparte eenweg-ANOVA’s in plaats van één MANOVA.

a) Bij $\alpha = .05$ per ANOVA en drie aparte ANOVA’s: hoe groot is de family-wise type-I-fout?

b) Welke informatie bevat een MANOVA dat drie ANOVA’s missen?

c) Wanneer mag je tóch drie aparte ANOVA’s runnen? (Hint: na een significante MANOVA, met Bonferroni-correctie.)

OpgeletAntwoord T2 — open na je eigen poging

a) $1 - (1 - .05)^3 = 1 - .857 = .143$ — ruim $14\%$ kans op minstens één vals-positief, ongeveer drie keer zo hoog als de nominale $.05$.

b) De correlatie-structuur tussen de DV’s. Een MANOVA leest dat de drie DV’s samen een patroon vormen — bv. dat alpha’s tegelijk hoog op trots en laag op eenzaamheid scoren. Drie ANOVA’s missen dat samen-patroon.

c) Na een significante multivariate $F$-toets: dan rapporteer je univariate $F$’s als follow-up, met Bonferroni-correctie (per DV de drempel $.05/p$). Dit heet de protected F-route, maar de naam belooft te veel: de MANOVA vooraf beschermt maar half. Verschillen de groepen op geen enkele DV, dan houdt de MANOVA je meestal tegen. Maar verschilt er één DV écht, dan komt de MANOVA al door de drempel dankzij die ene — en daarna toets je óók de DV’s die niet verschillen, elk met zijn eigen kans van $.05$ op een vals alarm. Daarom doet de Bonferroni-correctie hier het beschermende werk, niet de MANOVA. De route blijft wel staan: eerst de MANOVA, dan per DV met Bonferroni. De MANOVA is daarin niet overbodig — hij beantwoordt de gezamenlijke vraag, en vindt ook het geval dat losse ANOVA’s missen. Alleen het beschermen van de losse toetsen laat hij aan Bonferroni over.

Inzicht. MANOVA is geen “betere” ANOVA — het is een ander type vraag (gezamenlijk profiel) met andere aannames (multivariate normaliteit, gelijke covariantie-matrices) en eigen interpretatie.

5.1.a Verkenning — DV-means per groep en correlatiematrix

Voordat de leeuw fit, kijkt hij eerst.

Algemene vorm.

# Gemiddelden per groep voor elke DV.
aggregate(cbind(Y1, Y2, Y3, Y4) ~ factor, data = mijn_data, FUN = mean)

# Correlatiematrix tussen de DV's.
cor(mijn_data[, c("Y1", "Y2", "Y3", "Y4")])

# Boxplot per DV.
boxplot(Y1 ~ factor, data = mijn_data)

Voor onze dieren. Waar je data samenvat — per groep middelen, kolommen kiezen — staan in dit hoofdstuk twee schrijfwijzen naast elkaar, in tabbladen. Base R staat vooraan en is wat je nodig hebt; dplyr is er om te kijken. Beide geven dezelfde getallen.

OpmerkingHet pijpje |> in elk dplyr-blok

|> geeft door wat links staat aan de functie rechts. Lees het als “en dan”.

# Deze twee doen precies hetzelfde.
round(mean(c(1, 2, 3)), 1)
c(1, 2, 3) |> mean() |> round(1)

Base R

# Gemiddelden per groep voor elke DV.
aggregate(cbind(trots, jaloezie, eenzaamheid, woede) ~ sociale_status,
          data = leeuw_gevoelens, FUN = mean)

  sociale_status trots jaloezie eenzaamheid woede
1          alpha    80       40          25    45
2         midden    55       75          50    70
3       outsider    25       45          75    30

# Correlatiematrix tussen de vier DV's.
round(cor(leeuw_gevoelens[, c("trots", "jaloezie", "eenzaamheid", "woede")]), 2)

            trots jaloezie eenzaamheid woede
trots        1.00    -0.12       -0.91  0.32
jaloezie    -0.12     1.00        0.16  0.79
eenzaamheid -0.91     0.16        1.00 -0.28
woede        0.32     0.79       -0.28  1.00

dplyr

# group_by() knipt de tabel in groepen; summarise() vat elke groep samen tot
# één regel. across() zegt: doe dat voor deze vier kolommen tegelijk.
leeuw_gevoelens |>
  group_by(sociale_status) |>
  summarise(across(c(trots, jaloezie, eenzaamheid, woede), mean))

# A tibble: 3 × 5
  sociale_status trots jaloezie eenzaamheid woede
  <fct>          <dbl>    <dbl>       <dbl> <dbl>
1 alpha             80       40          25    45
2 midden            55       75          50    70
3 outsider          25       45          75    30

# select() kiest de vier DV's; cor() en round() krijgen ze doorgegeven.
leeuw_gevoelens |>
  select(trots, jaloezie, eenzaamheid, woede) |>
  cor() |>
  round(2)

            trots jaloezie eenzaamheid woede
trots        1.00    -0.12       -0.91  0.32
jaloezie    -0.12     1.00        0.16  0.79
eenzaamheid -0.91     0.16        1.00 -0.28
woede        0.32     0.79       -0.28  1.00

par(mfrow = c(1, 4), mar = c(4, 4, 2.5, 1))
boxplot(trots ~ sociale_status, data = leeuw_gevoelens,
        ylab = "Trots", xlab = "", ylim = c(0, 100),
        col = c("#0077BB", "#33BBEE", "#EE7733"),
        main = "Trots")
boxplot(jaloezie ~ sociale_status, data = leeuw_gevoelens,
        ylab = "Jaloezie", xlab = "", ylim = c(0, 100),
        col = c("#0077BB", "#33BBEE", "#EE7733"),
        main = "Jaloezie")
boxplot(eenzaamheid ~ sociale_status, data = leeuw_gevoelens,
        ylab = "Eenzaamheid", xlab = "", ylim = c(0, 100),
        col = c("#0077BB", "#33BBEE", "#EE7733"),
        main = "Eenzaamheid")
boxplot(woede ~ sociale_status, data = leeuw_gevoelens,
        ylab = "Woede", xlab = "", ylim = c(0, 100),
        col = c("#0077BB", "#33BBEE", "#EE7733"),
        main = "Woede")

par(mfrow = c(1, 1))

DV-means-tabel — kruispunt met gt

means_df <- aggregate(cbind(trots, jaloezie, eenzaamheid, woede) ~ sociale_status,
                      data = leeuw_gevoelens, FUN = mean)
sds_df   <- aggregate(cbind(trots, jaloezie, eenzaamheid, woede) ~ sociale_status,
                      data = leeuw_gevoelens, FUN = sd)
ns_df    <- aggregate(trots ~ sociale_status,
                      data = leeuw_gevoelens, FUN = length)

tab_df <- data.frame(
  # als tekst, niet als factor: gt centreert een factorkolom, en rijnamen horen links.
  status      = as.character(means_df$sociale_status),
  n           = ns_df$trots,
  trots_m     = means_df$trots,
  trots_sd    = sds_df$trots,
  jaloez_m    = means_df$jaloezie,
  jaloez_sd   = sds_df$jaloezie,
  eenz_m      = means_df$eenzaamheid,
  eenz_sd     = sds_df$eenzaamheid,
  woede_m     = means_df$woede,
  woede_sd    = sds_df$woede
)

# Voeg overall-rij toe
overall <- data.frame(
  status   = "Totaal",
  n        = nrow(leeuw_gevoelens),
  trots_m  = mean(leeuw_gevoelens$trots),
  trots_sd = sd(leeuw_gevoelens$trots),
  jaloez_m = mean(leeuw_gevoelens$jaloezie),
  jaloez_sd= sd(leeuw_gevoelens$jaloezie),
  eenz_m   = mean(leeuw_gevoelens$eenzaamheid),
  eenz_sd  = sd(leeuw_gevoelens$eenzaamheid),
  woede_m  = mean(leeuw_gevoelens$woede),
  woede_sd = sd(leeuw_gevoelens$woede)
)
tab_df <- rbind(tab_df, overall)

# In huisstijl (KAAPA, via ../common/R/kaapa.R). Variabelenamen in de
# kolomkoppen rechtop: cursief is voor statistische symbolen (M, SD, n).
tab_df |>
  gt() |>
  tab_spanner(label = "Trots",       columns = c(trots_m, trots_sd)) |>
  tab_spanner(label = "Jaloezie",    columns = c(jaloez_m, jaloez_sd)) |>
  tab_spanner(label = "Eenzaamheid", columns = c(eenz_m, eenz_sd)) |>
  tab_spanner(label = "Woede",       columns = c(woede_m, woede_sd)) |>
  fmt_number(columns = c(trots_m, trots_sd, jaloez_m, jaloez_sd,
                         eenz_m, eenz_sd, woede_m, woede_sd),
             decimals = 2) |>
  cols_label(
    status    = "Sociale status",
    n         = md("*n*"),
    trots_m   = md("*M*"),  trots_sd  = md("*SD*"),
    jaloez_m  = md("*M*"),  jaloez_sd = md("*SD*"),
    eenz_m    = md("*M*"),  eenz_sd   = md("*SD*"),
    woede_m   = md("*M*"),  woede_sd  = md("*SD*")
  ) |>
  gt_apa(title = "Vier gevoelens per statusgroep") |>
  tab_footnote_apa(paste0(
    "Elk gevoel is gemeten op een schaal van 0 tot 100. *n* = aantal leeuwen; ",
    "*M* = gemiddelde; *SD* = standaardafwijking."))

Vier gevoelens per statusgroep
Sociale status	n	Trots		Jaloezie		Eenzaamheid		Woede
		M	SD	M	SD	M	SD	M	SD
alpha	36	80.00	9.36	40.00	8.54	25.00	8.79	45.00	9.15
midden	38	55.00	8.54	75.00	7.65	50.00	7.73	70.00	7.84
outsider	36	25.00	8.80	45.00	10.08	75.00	8.72	30.00	10.51
Totaal	110	53.36	24.06	53.73	17.92	50.00	21.96	48.73	19.01
Note. Elk gevoel is gemeten op een schaal van 0 tot 100. n = aantal leeuwen; M = gemiddelde; SD = standaardafwijking.

OpmerkingVragen 5.1

a) Hoeveel leeuwen heeft de leeuw genoteerd? Hoe groot is elke statusgroep?

b) In welke groep is het gemiddelde trots het hoogst? Het laagst?

c) In welke groep is jaloezie het hoogst? Past dat bij de verwachting?

d) In welke groep is eenzaamheid het hoogst? In welke groep is woede het hoogst? Welke combinatie van DV’s onderscheidt outsiders van alpha’s?

e) Wat zegt de correlatiematrix? Welke twee DV’s hangen het sterkst samen? Welke twee zijn elkaars tegenpolen?

OpgeletAntwoord 5.1 — open na je eigen poging

In gewone woorden.

1. $N = 110$ leeuwen, met $36$ alpha’s, $38$ middenrangers en $36$ outsiders — vrijwel balanced ($n_{\max}/n_{\min} = 38/36 \approx 1.06$, ruim onder de $1.5$-drempel).

2. Trots is het hoogst bij alpha’s ($M = 80.00$) en het laagst bij outsiders ($M = 25.00$).

3. Jaloezie is het hoogst bij middenrangers ($M = 75.00$) — past op het Sapolsky-baboon-patroon: middenrangers kijken omhoog naar alpha’s en zijn bedreigd van onder, dat genereert de meeste sociale-vergelijking-stress.

4. Eenzaamheid is het hoogst bij outsiders ($M = 75.00$) en het laagst bij alpha’s ($M = 25.00$). Woede is het hoogst bij middenrangers ($M = 70.00$) en het laagst bij outsiders ($M = 30.00$); alpha’s zitten ertussen ($M = 45.00$) — defensief, niet plat nul. Outsiders verschillen van alpha’s vooral op trots (laag versus hoog) en eenzaamheid (hoog versus laag); middenrangers verschillen van de andere twee vooral op jaloezie en woede.

5. De correlatiematrix laat een sterk positieve $r$ tussen jaloezie en woede zien ($r = .79$, wedijver-cluster) en een sterk negatieve $r$ tussen trots en eenzaamheid ($r = -.91$, tegenpolen).

APA-stijl.

Van $N = 110$ leeuwen verschilden de vier gevoelens beschrijvend tussen statusgroepen: alpha’s scoorden het hoogst op trots ($M = 80.00$, $SD = 9.36$), middenrangers het hoogst op zowel jaloezie ($M = 75.00$, $SD = 7.65$) als woede ($M = 70.00$, $SD = 7.84$), en outsiders het hoogst op eenzaamheid ($M = 75.00$, $SD = 8.72$). De samenhang tussen de vier gevoelens liep sterk uiteen; het sterkst hingen jaloezie en woede samen ($r = .79$) en trots en eenzaamheid ($r = -.91$).

OpmerkingVoor wie verder wil — als twee DV’s bijna hetzelfde meten

Hangen twee DV’s heel sterk samen, zeg $|r| > .90$, dan meten ze bijna hetzelfde. MANOVA wint daar niets mee: twee keer bijna dezelfde meting is geen tweede blik op het profiel. En het rekenwerk wordt wankel. De vier toetsgrootheden delen als het ware door de samenhang tussen de DV’s binnen de groepen; vallen twee DV’s bijna samen, dan deel je door bijna niets, en dan schieten de eigenvalues en de gewichten van de discriminant-functies alle kanten op.

Een veelgebruikte vuistregel uit de handboeken: MANOVA heeft de meeste meerwaarde bij een samenhang tussen de DV’s van ongeveer $.20$ tot $.70$. Daaronder zeggen de DV’s elkaar weinig en leveren losse ANOVA’s bijna hetzelfde op; daarboven meten ze bijna hetzelfde. Wat je dan doet: maak van de twee één score (een gemiddelde, een som, of een hoofdcomponent), of laat er op inhoudelijke gronden één vallen.

Let op welke samenhang je bekijkt. Het gaat om de samenhang binnen de groepen. De correlatiematrix hierboven is berekend over alle leeuwen samen, en daar telt het verschil tussen de statusgroepen mee: trots en eenzaamheid komen daar uit op $r = -.91$. Binnen de groepen is dat $r = -.35$ — een gewone, matige samenhang, geen dubbele meting.

5.1.b Aannames van MANOVA

OpmerkingDe drie MANOVA-aannames in één oogopslag

Dit hoofdstuk werkt met drie aannames van MANOVA (uitgebreidere handboeken noemen er meer, zoals uitbijters en lineariteit):

1. Multivariate normaliteit — de errors per DV zijn normaal verdeeld, óók per deelgroep van personen met identieke scores op de andere DV’s. Te toetsen? In de praktijk niet. Robuust? Ja, mits $n \geq 20$ per groep.
2. Homogeniteit van variantie-covariantie-matrices — alle groepen hebben dezelfde $\boldsymbol{\Sigma}$. Te toetsen? Met Box’s M, maar zie de strenge-$\alpha$-vuistregel hieronder. Robuust? Ja, mits $n_{\max} / n_{\min} < 1.5$.
3. Onafhankelijke errors — error van persoon $i$ is onafhankelijk van error van persoon $j$. Te toetsen? Nee, dat zit in het onderzoeksdesign (bv. geen herhaalde metingen, geen geneste structuur).

In dit hoofdstuk besteden we de meeste aandacht aan Box’s M (de enige formele toets); MV-normaliteit en onafhankelijkheid behandelen we via vuistregels en design-overwegingen.

Geldigheid van deze check: de controle van deze drie aannames geldt voor de hele MANOVA-procedure inclusief follow-up. Je hoeft Box’s M of de normaliteits-vuistregel niet opnieuw te doen per univariate ANOVA in 5.2. Dat komt doordat de multivariate aannames de univariate al bevatten: gelijke covariantie-matrices betekent ook gelijke varianties per DV (die staan op de diagonaal), multivariaat normaal betekent ook normaal per DV, en onafhankelijke errors zitten in het design — dat is voor elke DV hetzelfde.

Box’s M — gelijke covariantie-matrices

OpmerkingVan varianties naar covariantie-matrices — Box’s M opbouwen

Bij ANOVA in thema 2 had je per groep één variantie te checken (Levene). Bij MANOVA hebben we meerdere DV’s per persoon, dus per groep:

• De variantie per DV (hoe gespreid is de score binnen die groep?).
• De covariantie tussen elke twee DV’s (gaat een hoge trots in deze groep samen met een hoge jaloezie, of juist niet?).

Voor onze leeuw met vier DV’s: per groep krijg je $4$ varianties (één per DV) en $6$ covarianties (één per DV-paar: trots–jaloezie, trots–eenzaamheid, trots–woede, jaloezie–eenzaamheid, jaloezie–woede, eenzaamheid–woede). Die tien getallen leg je vast in één $4 \times 4$-tabel: varianties op de diagonaal, covarianties off-diagonal (symmetrisch — de covariantie van DV $j$ met DV $k$ is dezelfde als die van $k$ met $j$).

Voor groep $g$ ziet die tabel er schematisch zo uit:

$$\begin{pmatrix} \text{var}(\text{trots})_g & \text{cov}(\text{trots, jal.})_g & \text{cov}(\text{trots, eenz.})_g & \text{cov}(\text{trots, woede})_g \\ \text{cov}(\text{trots, jal.})_g & \text{var}(\text{jal.})_g & \text{cov}(\text{jal., eenz.})_g & \text{cov}(\text{jal., woede})_g \\ \text{cov}(\text{trots, eenz.})_g & \text{cov}(\text{jal., eenz.})_g & \text{var}(\text{eenz.})_g & \text{cov}(\text{eenz., woede})_g \\ \text{cov}(\text{trots, woede})_g & \text{cov}(\text{jal., woede})_g & \text{cov}(\text{eenz., woede})_g & \text{var}(\text{woede})_g \end{pmatrix}$$

Bij $p$ DV’s wordt dat een $p \times p$-tabel. Dezelfde vorm als een correlatie-matrix uit thema 1 (MRA), maar met ruwe verbanden in de cellen (covarianties), niet gestandaardiseerde (correlaties). Een covariantie deel je door de SD’s van beide variabelen en je krijgt een correlatie — dezelfde informatie, andere schaal.

Deze tabel heet de variantie-covariantie-matrix (kortweg covariantie-matrix). Het wiskundige symbool is $\boldsymbol{\Sigma}$ — hoofdletter sigma, vetgedrukt om aan te geven dat het een matrix is, niet één getal. Per groep dus één eigen $\boldsymbol{\Sigma}_g$ (met $g$ = alpha, midden, outsider).

Wat toetst Box’s M dan? Of álle groepen dezelfde covariantie-matrix hebben:

$$H_0: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2 = \ldots = \boldsymbol{\Sigma}_I$$

In gewone woorden: ook de variantie van trots is in elke statusgroep gelijk; ook de covariantie tussen jaloezie en eenzaamheid; en zo voor élke cel van de matrix. De multivariate analoog van Levene’s $F$ — die toetste alleen de varianties per DV; Box’s M toetst tegelijk varianties én covarianties.

Box’s M is bekend overgevoelig — bij grote samples verwerpt hij vrijwel altijd, ook bij irrelevante schendingen. Vandaar de vuistregel: hanteer $\alpha = .001$ voor Box’s M (niet de gewone $.05$). Pas bij $p < .001$ ben je serieus bezorgd. Bij gebalanceerde celgroottes en $n_i \geq 20$ is de MANOVA-$F$ bovendien robuust tegen lichte schendingen — dan ga je gewoon door met de analyse.

Wat als Box’s M wél significant is op $\alpha = .001$? Kijk dan naar je celgroottes. Als $n_{\max} / n_{\min} < 1.5$ en alle $n_i \geq 20$ is de MANOVA-$F$ nog steeds robuust — vermeld de schending in je rapportage, maar stop niet met analyseren. Alleen als zowel $p < .001$ én zware onbalans (of kleine groepen) speelt, wordt het problematisch, en dan is er geen simpele oplossing. Bij die combinatie overweeg je robuuste varianten of bootstrap, of je herziet je design.

Kleine noot: in lopende tekst heet het concept meestal gewoon “covariance matrices”; het symbool $\boldsymbol{\Sigma}$ zie je vooral in formules. Wij introduceren het hier expliciet — kost één teken, levert helderheid op.

Algemene vorm.

# Box's M — drie of meer DV's, één factor.
rstatix::box_m(data = mijn_data[, c("Y1", "Y2", "Y3")],
               group = mijn_data$factor)

Voor onze dieren.

box_m(data  = leeuw_gevoelens[, c("trots", "jaloezie", "eenzaamheid", "woede")],
      group = leeuw_gevoelens$sociale_status)

# A tibble: 1 × 4
  statistic p.value parameter method                                            
      <dbl>   <dbl>     <dbl> <chr>                                             
1      14.4   0.808        20 Box's M-test for Homogeneity of Covariance Matric…

OpmerkingVragen 5.1

f) Wat is de Box’s M-statistiek, met welke $df$ en welke $p$-waarde?

g) Verwerp je $H_0$ van gelijke covariantie-matrices op $\alpha = .001$? Op $\alpha = .05$?

h) Wat is je conclusie over robuustheid: kun je doorgaan met MANOVA?

OpgeletAntwoord 5.1 — open na je eigen poging

In gewone woorden.

6. Box’s $M = 14.43$, $df = 20$, $p = .81$.

7. Niet significant op $\alpha = .001$, en ook niet op de gewone $\alpha = .05$: geen aanwijzing dat de covariantie-matrices tussen de groepen verschillen, geen reden tot zorg.

8. De celgroottes zijn vrijwel balanced ($n_i = 36$-$38$, ratio $1.06$) en alle $\geq 20$, dus zelfs als Box’s M wél significant was geweest, was de MANOVA-$F$ robuust geweest. We gaan met vertrouwen door.

APA-stijl.

Box’s $M$-test op homogeniteit van de covariantie-matrices was niet significant op $\alpha = .001$, $M = 14.43$, $df = 20$, $p = .81$. Gegeven de gebalanceerde celgroottes ($n_{\max}/n_{\min} = 1.06$) en $n_i \geq 36$ is de MANOVA-$F$-toets bovendien robuust tegen lichte schendingen.

T3 — Box’s M en zijn strenge $\alpha$

OpmerkingVraag T3 — Pen-en-papier

a) Wat toetst Box’s M precies, in eigen woorden?

b) Waarom hanteren we $\alpha = .001$ voor Box’s M en niet de gewone $.05$?

c) Een onderzoeker rapporteert Box’s $M = 122.9$, $df = 75$, $p = .0004$. Op $\alpha = .001$ wel of niet significant? Wat is zijn vervolgstap?

d) Welke twee andere kenmerken van het design helpen MANOVA robuust te maken bij borderline Box’s M?

OpgeletAntwoord T3 — open na je eigen poging

a) Box’s M test of de covariantie-matrices van de DV’s gelijk zijn over de groepen — dus zowel varianties per DV als covarianties tussen DV-paren binnen elke groep. De multivariate uitbreiding van Levene’s $F$ uit thema 2.

b) Box’s M is overgevoelig voor kleine, irrelevante schendingen, vooral bij grote samples. Bij $\alpha = .05$ verwerpt hij vrijwel altijd, zonder dat de schending praktisch problematisch is voor de MANOVA-$F$. De drempel $.001$ vangt alleen serieuze schendingen.

c) $p = .0004 < .001$, dus wel significant op de strenge drempel. Vervolgstap: kijken of celgroottes balanced en $\geq 20$ zijn. Zo ja, doorgaan met MANOVA en in de rapportage vermelden dat de schending bekend is. Zo nee, overweeg robuuste varianten of pseudo-bootstrap.

d) Gebalanceerde celgroottes ($n_{\max}/n_{\min} \leq 1.5$) en voldoende grote groepen ($n_i \geq 20$). Onder die twee voorwaarden wordt de MANOVA-$F$ robuust tegen schendingen van zowel normaliteit als covariantie-matrix-gelijkheid.

Inzicht. Box’s M is een waarschuwingslampje, niet een go/no-go-toets. Bij rood lampje plus klein design plus zware onbalans: pas op. Bij rood lampje plus balanced + grote $n$: ga door.

Multivariate normaliteit en onafhankelijke errors — niet toetsen, wel meedragen

OpmerkingDe andere twee aannames — hoe ga je ermee om?

Multivariate normaliteit. De errors per DV zijn normaal verdeeld, en — dit is het multivariate stuk — óók voor elke deelgroep van personen met identieke scores op de andere DV’s. Toetsen doe je dit niet, en je hoeft er niet wakker van te liggen: er is geen praktisch toetsbare procedure die robuust schaalt, en de MANOVA-$F$ is in de praktijk robuust tegen lichte schendingen mits $n \geq 20$ per groep.

Vuistregel: als alle $n_i \geq 20$, niet toetsen, gewoon doorgaan. Bij kleinere groepen: histogrammen of QQ-plots per DV inspecteren (univariate normaliteit als minimale check).

Onafhankelijke errors. De error van persoon $i$ is onafhankelijk van die van persoon $j$. Ook dit toets je niet; het zit in je onderzoeksontwerp, en daar is het wél belangrijk. Praktisch: geen herhaalde metingen, geen geneste structuur (broers/zussen, dezelfde school, dezelfde therapeut), geen tijdsvolgorde. Als je die wél hebt, is MANOVA niet de juiste keuze — dan wil je repeated measures (thema 6) of een mixed model.

Voor ons leeuw-design ($N = 110$ verschillende leeuwen uit losse troepen, één meting per dier): aanname plausibel, geen verdere check.

OpmerkingVoor wie verder wil — multivariate normaliteit en uitbijters zelf bekijken

Multivariaat normaal is meer dan “elke DV is los normaal verdeeld”. Ook elke gewogen optelling van de DV’s is dan normaal verdeeld, en de DV’s vormen samen een ellipsvormige puntenwolk. Bij twee DV’s kun je dat zien: een scatterplot per groep, en dan hoop je op een ellips — geen banaan, geen losse puntjes ver buiten de wolk.

Bij vier DV’s kun je die wolk niet meer tekenen. Het midden ervan kun je nog wel uitrekenen: het gemiddelde op elke DV. De Mahalanobis-afstand $D^2$ is de afstand van één leeuw tot dat midden, en hij telt de samenhang tussen de DV’s mee. Trots en eenzaamheid lopen binnen de groepen tegen elkaar in; een leeuw die op allebei hoog scoort, ligt daardoor verder buiten de wolk dan je aan de twee scores los zou zien.

Een gangbare grens: een $D^2$ boven de $.999$-grens van een $\chi^2$-verdeling met $p$ vrijheidsgraden is een uitbijter die je één voor één bekijkt. Bij onze vier DV’s ligt die grens op $18.47$.

En als het niet klopt? Bij gebalanceerde groepen van minstens twintig is de MANOVA-$F$ robuust tegen lichte schendingen. Bij een sterkere schending kun je een scheve DV transformeren (log, wortel), of uitwijken naar een robuuste MANOVA (zie Wat blijft liggen). Gebalanceerde, grote groepen blijven de beste verzekering — beter dan elk diagnostisch plaatje.

5.1.c MANOVA fitten

Algemene vorm.

# Stap 1: multivariaat lineair model.
lm_man <- lm(cbind(Y1, Y2, Y3) ~ factor, data = mijn_data)

# Stap 2: MANOVA met type-III SS (consistent met factorial design en onbalans).
res_Manova <- car::Manova(lm_man, type = 3)

# Stap 3: de vier multivariate toetsen.
summary(res_Manova, multivariate = TRUE, univariate = FALSE)

# Stap 4: de univariate F's per DV, uit summary.aov() — niet uit Manova().
summary.aov(lm_man)

cbind() plakt de DV’s tot één matrix-uitkomst — dat is hoe R weet dat je een multivariaat lineair model wilt. Manova() (hoofdletter $M$, uit car) doet de Type-III sums-of-squares-and-cross-products en levert alle vier teststatistics in één output.

Voor onze dieren.

# Defensief: zet sum-to-zero contrasten ook hier expliciet, vlak vóór de
# Type-III-aanroep — globale state kan tussentijds gewijzigd zijn door
# package-load of een eerdere chunk.
options(contrasts = c("contr.sum", "contr.poly"))

# Stap 1: lm-object met cbind() over de vier DV's.
lm_leeuw <- lm(cbind(trots, jaloezie, eenzaamheid, woede) ~ sociale_status,
               data = leeuw_gevoelens)

# Stap 2: MANOVA met Type-III SS (consistent met factorial designs en onbalans).
res_leeuw <- Manova(lm_leeuw, type = 3)

# Stap 3: samenvatting met multivariate output.
summary(res_leeuw, multivariate = TRUE, univariate = FALSE)

Type III MANOVA Tests:

Sum of squares and products for error:
                 trots  jaloezie eenzaamheid     woede
trots        8479.0327 -2666.984  -2774.7276 -885.9486
jaloezie    -2666.9836  8273.166   2401.9944 4452.3347
eenzaamheid -2774.7276  2401.994   7580.1544  908.4074
woede        -885.9486  4452.335    908.4074 9070.9927

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

Sum of squares and products for the hypothesis:
               trots jaloezie eenzaamheid    woede
trots       312685.7 312685.7    293142.9 283371.4
jaloezie    312685.7 312685.7    293142.9 283371.4
eenzaamheid 293142.9 293142.9    274821.4 265660.7
woede       283371.4 283371.4    265660.7 256805.4

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1   0.99389 4230.963      4    104 < 2.22e-16 ***
Wilks             1   0.00611 4230.963      4    104 < 2.22e-16 ***
Hotelling-Lawley  1 162.72933 4230.963      4    104 < 2.22e-16 ***
Roy               1 162.72933 4230.963      4    104 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: sociale_status 

Sum of squares and products for the hypothesis:
                 trots  jaloezie eenzaamheid     woede
trots        54605.455 -2929.091      -49500  16870.91
jaloezie     -2929.091 26721.818        4500  24921.82
eenzaamheid -49500.000  4500.000       45000 -13500.00
woede        16870.909 24921.818      -13500  30321.82

Multivariate Tests: sociale_status
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            2  1.699512 148.4656      8    210 < 2.22e-16 ***
Wilks             2  0.017492 170.5844      8    208 < 2.22e-16 ***
Hotelling-Lawley  2 15.178253 195.4200      8    206 < 2.22e-16 ***
Roy               2 11.664087 306.1823      4    105 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Stap 4: univariate F's apart via summary.aov() — zie waarschuwing hieronder.
summary.aov(lm_leeuw)

 Response trots :
                Df Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2  54605 27302.7  344.54 < 2.2e-16 ***
Residuals      107   8479    79.2                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response jaloezie :
                Df  Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2 26721.8 13360.9   172.8 < 2.2e-16 ***
Residuals      107  8273.2    77.3                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response eenzaamheid :
                Df Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2  45000 22500.0  317.61 < 2.2e-16 ***
Residuals      107   7580    70.8                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response woede :
                Df Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2  30322 15160.9  178.84 < 2.2e-16 ***
Residuals      107   9071    84.8                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hoe je deze uitvoer leest: de multivariate uitvoer begint met een blok Term: (Intercept). Dat toetst of de vier gemiddelden samen nul zijn — in de psychologie vrijwel nooit een vraag, dus dat blok sla je over. Jouw blok is Term: sociale_status, en daarin de regels onder Multivariate Tests.

WaarschuwingUnivariate $F$’s uit summary.aov(), niet uit de Manova-tabel

Je kunt summary() op een Manova()-object ook om de univariate $F$’s vragen, met univariate = TRUE. Doe dat niet: car rekent een deel van die $F$’s verkeerd uit zodra de termen in het model ongelijke $df$ hebben. Bij type = 3 is dat al zo bij één factor met drie of meer groepen: het intercept staat dan als term met $df = 1$ in de tabel, naast de factor met $df = 2$. In E5 aan het eind van dit hoofdstuk geeft car zo voor voorzichtigheid $F = 62.40$, waar $F = 31.20$ klopt — precies twee keer te groot. Welke $F$’s het treft, hangt af van het aantal DV’s: bij onze leeuw zijn het trots en eenzaamheid, in E5 alleen voorzichtigheid. Vertrouw die tabel dus in zijn geheel niet. Vandaar dat we hierboven univariate = FALSE zetten en de univariate $F$’s via summary.aov(lm_leeuw) halen. Die kloppen, ook met TukeyHSD() verderop.

Ook in ander lesmateriaal. Uitvoer die met univariate = TRUE is gemaakt, kom je ook in lesmateriaal en uitwerkingen tegen, met dezelfde fout erin — en soms verandert daardoor welke DV significant is. Neem zo’n $F$, en de $p$ ernaast, dus niet over zonder na te rekenen: de kwadratensom van de factor gedeeld door zijn $df$, gedeeld door de kwadratensom van de residuen gedeeld door hun $df$. Die kwadratensommen staan in dezelfde uitvoer, en die kloppen wel. Heb je de data zelf, dan geeft summary.aov() de goede $F$ meteen. Komt er een andere $F$ uit dan in de uitwerking, dan heb jij goed gerekend.

TipType-III SS — contr.sum voor multivariate analyses

Manova(., type = 3) werkt het cleanst met sum-to-zero contrasten (contr.sum) in plaats van de R-default contr.treatment. Daarom de options(contrasts = ...) regel boven het model. Vergelijk thema 2 (ANOVA) waar dezelfde keuze gemaakt werd voor consistente Type-III SS bij interactietermen — dezelfde logica geldt hier voor het multivariate geval. Bij eenwegs-MANOVA met balanced design maakt het cijfermatig nauwelijks uit; bij factoriële MANOVA (zie Factoriële MANOVA verderop) wordt het belangrijk.

Na de analyse kun je terugzetten met options(contrasts = c("contr.treatment", "contr.poly")).

5.1.d De vier teststatistics — Pillai eerst, de andere drie ernaast

OpmerkingVragen 5.1

i) Voor de term sociale_status: welke vier teststatistics rapporteert de output? Welke waardes hebben ze?

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

k) Welke teststatistic rapporteer je eerst? Waarom? En wat doe je met de andere drie?

OpgeletAntwoord 5.1 — open na je eigen poging

In gewone woorden.

1. Voor sociale_status: Pillai’s $V = 1.70$, Wilks’ $\Lambda = .017$, Hotelling-Lawley $T^2 = 15.18$, Roy’s $\theta = 11.66$.

10. $df_H = 8$ voor Pillai/Wilks/Hotelling (twee niveau-vergelijkingen × vier DV’s), $df_E$-approximaties iets verschillend per statistic; alle vier $p < .001$ — allemaal in dezelfde richting, sterk multivariaat effect.

11. Pillai’s trace rapporteer je eerst: hij is het robuustst tegen lichte schendingen van multivariate normaliteit en covariantie-matrix-gelijkheid, en hij staat in de uitvoer ook bovenaan. De andere drie zet je ernaast, zoals in de APA-zin hieronder: bij twee groepen komen alle vier op dezelfde $F$ en $p$ uit, maar bij drie of meer groepen kunnen ze uiteenlopen, en wie ze alle vier ziet, ziet of de conclusie stabiel is.

APA-stijl.

Een eenwegs-MANOVA toonde een zeer sterk multivariaat effect van sociale status op het gevoels-profiel: Pillai’s $V = 1.70$, $F(8, 210) = 148.47$, $p < .001$. De andere drie teststatistics wezen in dezelfde richting (Wilks’ $\Lambda = .017$, $F(8, 208) = 170.58$, $p < .001$; Hotelling-Lawley $T^2 = 15.18$, $F(8, 206) = 195.42$, $p < .001$; Roy’s $\theta = 11.66$, $F(4, 105) = 306.18$, $p < .001$).

OpmerkingVoor wie verder wil — hoe groot is het multivariate effect?

Pillai’s $V$ ligt tussen $0$ en $s = \min(I-1, p)$. Bij twee groepen is dat maximum $1$, maar bij drie of meer groepen kan $V$ boven de $1$ uitkomen — bij onze leeuw tot $\min(2, 4) = 2$. Daarom deel je $V$ eerst door dat maximum. Wat je dan krijgt, $V/s$, ligt weer tussen $0$ en $1$ en lees je als de $\eta^2$ uit de univariate ANOVA: het is de multivariate partiële $\eta^2$ (zo rapporteert SPSS hem ook). Cohens vuistregels gelden voor dát getal:

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

Voor onze leeuw: $V = 1.70$ en $s = 2$, dus $V/s = .85$ — ver boven de grote-effect-grens, een zeer sterk effect.

Per DV kun je de $\eta^2$ uit de univariate ANOVA’s uitrekenen: de kwadratensom van de groepen gedeeld door de totale kwadratensom van die DV. In een eenwegs-ontwerp is dat hetzelfde getal als de partiële $\eta^2$, $SS_{H, j}/(SS_{H, j} + SS_{E, j})$, want daar vormen de groepen en de error samen precies het totaal.

# eta^2 per DV: SS van de groepen gedeeld door de totale SS van die DV.
sa <- summary.aov(lm_leeuw)
eta_per_dv <- function(aov_tab) {
  ss <- aov_tab[["Sum Sq"]]
  ss[1] / (ss[1] + ss[2])
}
eta_p <- setNames(sapply(sa, eta_per_dv), dv_leeuw)
round(eta_p, 2)   # twee decimalen, zoals in de tekst

      trots    jaloezie eenzaamheid       woede 
       0.87        0.76        0.86        0.77 

Sociale status verklaart bij trots $.87$ van de variantie, bij jaloezie $.76$, bij eenzaamheid $.86$ en bij woede $.77$. Trots en eenzaamheid het meest, jaloezie en woede iets minder — dezelfde volgorde als de univariate $F$’s verderop.

OpmerkingVoor wie verder wil — waar de vrijheidsgraden van de multivariate $F$ vandaan komen

Voor de eenwegs-MANOVA met factor sociale_status en vier DV’s:

Bron	Formule	Hier
Hypothesis $df$ (Pillai/Wilks)	$df_H = (I-1) \cdot p$	$(3-1) \cdot 4 = 8$
Error $df$ ($F$-approximatie Pillai)	$df_E = s \cdot (N - I - p + s)$, met $s = \min(I-1, p)$	$2 \cdot (110 - 3 - 4 + 2) = 210$ (varieert per statistic — Wilks $208$, Hotelling $206$, Roy $105$)
Univariaat $df$ per DV (de $F$ van één ANOVA)	teller $I-1$, noemer $N-I$	$2$ en $107$

Bij $I = 3$, $N = 60$, $p = 2$: $df_H = 4$, univariaat $F(2, 57)$ per DV — kleinere data, hetzelfde patroon.

T4 — De vier teststatistics: wanneer welke?

OpmerkingVraag T4 — Pen-en-papier

a) Wat is de relatie tussen Pillai’s $V$ en Wilks’ $\Lambda$? Waarom geeft een lage $\Lambda$ een hoog $V$?

b) In welke situatie geeft Pillai $p > .05$ terwijl Roy $p < .05$? (Hint: Roy ziet alleen de grootste eigenvalue.)

c) Welke teststatistic rapporteer je eerst, en wat doe je met de andere drie?

OpgeletAntwoord T4 — open na je eigen poging

a) Ze kijken naar hetzelfde, van twee kanten. Wilks’ $\Lambda$ is het deel van de spreiding dat de groepen níét verklaren: hoe kleiner, hoe sterker het effect. Pillai’s $V$ telt op wat de groepen wél verklaren: hoe groter, hoe sterker. Wordt het effect sterker, dan gaat $\Lambda$ omlaag en $V$ omhoog. (Wie de eigenvalues $\mathit{\lambda}_i$ van de discriminant-functies verderop al kent: $V = \sum \frac{\mathit{\lambda}_i}{1 + \mathit{\lambda}_i}$ en $\Lambda = \prod \frac{1}{1 + \mathit{\lambda}_i}$.)

b) Als de groepen voornamelijk verschillen langs één discriminant-dimensie (één grote eigenvalue, andere klein), is Roy’s $\theta$ hoog — gevoelig voor die ene as. Pillai middelt over alle dimensies en is dan minder uitgesproken. Roy kan zo “te enthousiast” significant rapporteren bij smalle effecten; Pillai blijft conservatief.

c) Eerst: Pillai’s $V$ — robuustst tegen aanname-schendingen, balanced en kleine samples. Ernaast de andere drie: als ze ruwweg dezelfde $p$-waarde geven, is je conclusie stabiel. Als Pillai en Roy fors verschillen, zit je effect in één dimensie en moet je dat in de discussie melden.

Inzicht. De vier teststatistics zijn niet “vier verschillende toetsen”, maar vier samenvattingen van dezelfde eigenvalues — die van de discriminant-functies. Bij sterke effecten convergeren ze; bij borderline of dimensionaal-smalle effecten lopen ze uiteen.

OpmerkingVoor wie verder wil — waarom $T^2$ en $\theta$ bij het intercept gelijk zijn

In de uitvoer van 5.1.c zijn Hotelling-Lawley en Roy gelijk voor het intercept ($T^2 = 162.73$, $\theta = 162.73$), maar verschillend voor sociale_status ($T^2 = 15.18$, $\theta = 11.66$). Bij het intercept is $df_H = 1$, dus er is maar één eigenvalue die niet nul is. $T^2$ telt de eigenvalues op en $\theta$ neemt de grootste — bij één eigenvalue is dat allebei diezelfde ene. Pillai en Wilks rekenen die eigenvalue om, $V = \mathit{\lambda}/(1+\mathit{\lambda})$ en $\Lambda = 1/(1+\mathit{\lambda})$, en komen dus níét op dezelfde waarde uit. Bij sociale_status ($df_H = 2$) zijn er twee eigenvalues; $T^2$ telt ze op, $\theta$ neemt alleen de grootste.

BelangrijkDon’t: na MANOVA klakkeloos elke DV univariaat doortoetsen

De multivariate $F$ heeft je toegang gegeven tot meerdere afzonderlijke ANOVA’s. Mooi. Niet familywise corrigeren is je eigen Type-I-foutkans opwaarderen — bij drie DV’s en $\alpha = .05$ per stuk wordt familywise $\alpha \approx .14$. Vier DV’s: $\approx .19$. Vijf: $\approx .23$.

Bonferroni $\mathit{\alpha} / k$, Holm-stepdown of equivalent — minimum. Wie drie univariate $F$’s rapporteert na een significante Pillai zonder enige correctie, hoopt erop dat de reviewer het over het hoofd ziet. Dat is geen wetenschap, dat is gokken op de gunfactor van een ander.

5.2 Follow-up na significante MANOVA

“Welke gevoelens zijn de echte motoren?”

“Het profiel verschilt,” zei de leeuw. “Maar verschilt elk gevoel afzonderlijk? En als ze allemaal verschillen, zit dat verschil tussen alpha en outsider, of binnen de middenrang?”

De multivariate $F$ zegt alleen dát het profiel verschilt. Voor welke DV’s en tussen welke groepen heb je follow-up nodig. De gangbare route bij MANOVA:

1. Univariate $F$ per DV met Bonferroni-correctie op $\alpha$.
2. Bij significante DV: Tukey HSD of vergelijkbare post-hoc binnen die DV om de groepsverschillen te lokaliseren.
3. Alternatief: descriptive discriminant analysis DDA — welke DV’s dragen het meest bij aan de groepsscheiding (zie 5.2.c).

TipDe drie tabellen in 5.2 — herken welke je leest

In de R-output zie je in 5.2 drie aparte tabellen. Eén van de meest voorkomende fouten bij follow-up-vragen is: in de verkeerde tabel kijken. Voor je een uitspraak doet over een groepsverschil, check eerst welke tabel je voor je hebt — lees de kolomtitels hardop voor jezelf.

Tabel	Wat je ziet	Wat het toetst	Hoe lees je ’m?
(1) Multivariate test (uit 5.1.c)	Pillai’s trace, Wilks $\Lambda$, $F$-approximatie, $p$	Verschilt het profiel over alle DV’s ergens tussen groepen?	Eén regel per teststatistic. Pillai eerst, de andere drie ernaast.
(2) Univariate ANOVA per DV (5.2.a)	$F$, $df_1$, $df_2$, $p$ per rij = DV	Verschilt deze ene DV univariaat tussen groepen?	Eén regel per DV. Vergelijk $p$ met Bonferroni-drempel $\mathit{\alpha} / k$.
(3) Pairwise comparisons (5.2.b, Tukey HSD)	Verschil, CI, $p$ per rij = groepenpaar binnen één DV	Verschilt deze specifieke twee groepen op deze DV?	Eén regel per paar (alpha-mid, alpha-out, mid-out). Per significante DV apart.

Volgorde: eerst (1) — niet sig → klaar. (1) sig → (2). Per DV in (2) sig → (3) op die DV. Niet andersom; en niet door elkaar.

Korte controle vóór je een conclusie schrijft: “Welke tabel heb ik voor me? Wat is hier de eenheid van rij? Wat heb ik nodig voor mijn uitspraak?”

5.2.a Univariate ANOVA per DV met Bonferroni

TipDit ken je al van thema 2 — één klein verschil

De univariate ANOVA’s hier zijn gewoon de techniek uit thema 2: één DV per analyse, $F$-toets, de robuustheid-vuistregels uit thema 2. Bij univariate ANOVA is er geen covariantie-matrix tussen DV’s — dus geen Box’s M nodig. De ANOVA-$F$ is robuust onder de bekende voorwaarden: balanced design en voldoende $n$ per groep.

Het enige verschil met thema 2: je doet hier meerdere ANOVA’s tegelijk (één per DV), en dat verhoogt de familywise foutkans. Vandaar Bonferroni of Holm-Bonferroni op $\alpha$ — strenger dan één losse ANOVA, maar precies dezelfde mechaniek per ANOVA. Alle andere aannames-checks zijn al gedaan in stap 0 (zie 5.1.b).

Algemene vorm.

# Univariate F per DV uit het multivariate lm-object.
summary.aov(lm_man)

# Of per DV apart:
summary(aov(Y1 ~ factor, data = mijn_data))
summary(aov(Y2 ~ factor, data = mijn_data))
summary(aov(Y3 ~ factor, data = mijn_data))

# Bonferroni-correctie: alpha-per-DV = .05 / aantal DV's.
# Vergelijk de p-waardes met .017 (bij 3 DV's), niet .05.

Voor onze dieren.

summary.aov(lm_leeuw)

 Response trots :
                Df Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2  54605 27302.7  344.54 < 2.2e-16 ***
Residuals      107   8479    79.2                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response jaloezie :
                Df  Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2 26721.8 13360.9   172.8 < 2.2e-16 ***
Residuals      107  8273.2    77.3                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response eenzaamheid :
                Df Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2  45000 22500.0  317.61 < 2.2e-16 ***
Residuals      107   7580    70.8                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response woede :
                Df Sum Sq Mean Sq F value    Pr(>F)    
sociale_status   2  30322 15160.9  178.84 < 2.2e-16 ***
Residuals      107   9071    84.8                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Verzamel de p-waardes en pas Bonferroni toe.
ps <- sapply(summary.aov(lm_leeuw), function(x) x[["Pr(>F)"]][1])
round(ps, 5)

       Response trots     Response jaloezie  Response eenzaamheid 
                    0                     0                     0 
       Response woede 
                    0 

# Bonferroni-corrected (vermenigvuldigd met p, gecapped op 1).
round(p.adjust(ps, method = "bonferroni"), 5)

       Response trots     Response jaloezie  Response eenzaamheid 
                    0                     0                     0 
       Response woede 
                    0 

# Drempel handmatig: .05 gedeeld door het aantal DV's (hier vier).
.05 / length(ps)

[1] 0.0125

OpmerkingVragen 5.2

a) Welke univariate $F$-waardes zie je per DV? Met welke $df$?

b) Welke DV’s zijn significant op de gewone $\alpha = .05$? Welke ook na Bonferroni-correctie ($\alpha = .0125$)?

c) Wat zegt het feit dat alle vier nog significant blijven na Bonferroni?

OpgeletAntwoord 5.2 — open na je eigen poging

In gewone woorden.

1. $F_{\text{trots}}(2, 107) = 344.54$, $F_{\text{jaloezie}}(2, 107) = 172.80$, $F_{\text{eenzaamheid}}(2, 107) = 317.61$, $F_{\text{woede}}(2, 107) = 178.84$ — alle vier extreem sterk significant met $p < .001$.

2. Op $\alpha = .05$: alle vier significant. Bonferroni-drempel $.05/4 = .0125$: ruim overleefd voor alle vier.

3. Dat de statusgroepen op elk van de vier gevoelens apart verschillen, ook onder de strengere drempel. Het multivariate verschil hangt dus niet aan één gevoel. Wat het niet zegt: dat de vier los van elkaar bijdragen. Een univariate $F$ kijkt naar één DV tegelijk en ziet de andere drie niet — en die hangen samen: jaloezie en woede met $r = .79$, trots en eenzaamheid met $r = -.91$ (de correlatiematrix uit 5.1.a). Wat elk gevoel toevoegt bovenop de andere, laat de DDA in 5.2.c beter zien.

APA-stijl.

De $\eta^2$ achter elke $F$ is dezelfde als in thema 2: de kwadratensom van de groepen gedeeld door de totale kwadratensom van die DV (Sum Sq van sociale_status gedeeld door de som van de twee Sum Sq’s in de summary.aov()-uitvoer).

Univariate follow-up met Bonferroni-correctie ($\alpha_{\text{per DV}} = .0125$) liet zien dat de statusgroepen op elk van de vier gevoelens significant verschilden: trots, $F(2, 107) = 344.54$, $p < .001$, $\eta^2 = .87$; jaloezie, $F(2, 107) = 172.80$, $p < .001$, $\eta^2 = .76$; eenzaamheid, $F(2, 107) = 317.61$, $p < .001$, $\eta^2 = .86$; woede, $F(2, 107) = 178.84$, $p < .001$, $\eta^2 = .77$.

T5 — Bonferroni hand-rekenen

OpmerkingVraag T5 — Pen-en-papier

Een onderzoeker heeft $5$ DV’s en doet voor elke DV een aparte ANOVA na een significante MANOVA.

a) Bij gewone $\alpha = .05$ en $5$ aparte ANOVA’s: hoe groot is de family-wise type-I-fout (zonder correctie)?

b) Bonferroni-correctie: wat is de drempel per DV?

c) Een DV heeft $p = .024$. Significant op de Bonferroni-corrected drempel?

d) Een DV heeft $p = .008$. Significant op de Bonferroni-corrected drempel?

OpgeletAntwoord T5 — open na je eigen poging

a) $1 - (1-.05)^5 = 1 - .774 = .226$ — bijna $23\%$ kans op minstens één vals-positief.

b) De drempel per DV is $.05 / 5 = .010$.

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

d) $p = .008 < .010$ — wel significant na Bonferroni.

Inzicht. Bonferroni is conservatief: je verlaagt de drempel evenredig met het aantal toetsen. Voordeel: simpel en gegarandeerd. Nadeel: bij veel DV’s wordt de drempel zo streng dat je power inboet — vandaar Holm-Bonferroni als minder strenge variant (zie T6).

T6 — Holm-Bonferroni als krachtigere variant

OpmerkingHet Holm-Bonferroni-algoritme — sequentieel, één-voor-één

Holm-Bonferroni is strikt sequentieel: je beoordeelt per stap één $p$-waarde, niet meerdere tegelijk. De koppeling tussen stap-nummer en drempel is vast: bij stap $h$ vergelijk je $p_{(h)}$ met $\mathit{\alpha}_h$ — nooit kruislings. $p_{(1)}$ gaat altijd met $\mathit{\alpha}_1$, $p_{(2)}$ altijd met $\mathit{\alpha}_2$, enzovoort.

Het algoritme, stap voor stap:

• Stap 0 (preparatie): Sorteer de $k$ $p$-waardes van klein naar groot: $p_{(1)} \leq p_{(2)} \leq \ldots \leq p_{(k)}$. (In veel beschrijvingen is dit geen eigen genummerde stap — het is een preparatie. Wij zetten het er als nul-stap bij om de eerste actie scherp te markeren.)
• Stap 1: vergelijk $p_{(1)}$ met $\mathit{\alpha}_1 = \mathit{\alpha} / k$. Als $p_{(1)} < \mathit{\alpha}_1$, $p_{(1)}$ significant — door naar stap 2. Anders: stop, niets significant.
• Stap 2: vergelijk $p_{(2)}$ met $\mathit{\alpha}_2 = \mathit{\alpha} / (k - 1)$. Als $p_{(2)} < \mathit{\alpha}_2$, significant — door. Anders: stop.
• Stap $h$: vergelijk $p_{(h)}$ met $\mathit{\alpha}_h = \mathit{\alpha} / (k - (h-1))$. Doorgaan of stoppen.

Hier is $k$ het aantal toetsen, en in de noemer staat $k$ minus het aantal stappen dat al is afgerond ($= h-1$). Je ziet ook wel $\mathit{\alpha} / (k - h + 1)$ — algebraïsch identiek, twee schrijfwijzen, één betekenis.

Kleine noot: in veel beschrijvingen krijgt de $p$-waarde geen nummer (er staat dan “de kleinste $p$”, “de op één na kleinste”). Wij gebruiken $p_{(h)}$ — de $h$-de kleinste $p$-waarde — omdat dat met één teken extra de koppeling tussen stap-nummer en $p$-waarde zichtbaar maakt.

De drempel wordt per stap ruimer (de noemer wordt kleiner). Reden: na een geconcludeerd-significante test heb je voor de overige tests een minder strenge correctie nodig — er resteren minder hypotheses om gelijktijdig te controleren.

Veelvoorkomende verwarring 1 — “$p_{(2)}$ ook eerst met $\mathit{\alpha}_1$ vergelijken?” Nee. Elke $p$-waarde heeft één eigen drempel die hoort bij de stap waarin hij beoordeeld wordt. $p_{(2)}$ wordt nooit met $\mathit{\alpha}_1$ vergeleken. Pas in stap 2, met $\mathit{\alpha}_2$.

Veelvoorkomende verwarring 2 — “Als bij stap 1 al twee p-waardes onder de drempel $\mathit{\alpha}_1$ liggen, mag ik er dan twee tegelijk uithalen en bij stap 2 direct $\mathit{\alpha} / (k - 2)$ pakken?” Nee. Het algoritme kijkt per stap alleen naar de kleinste resterende $p$-waarde, en gebruikt het stap-nummer (niet het aantal-al-significant) om de drempel te bepalen. Pas in stap 2 beoordeel je de tweede-kleinste, met $\mathit{\alpha}_2 = \mathit{\alpha} / (k - 1)$.

OpmerkingVraag T6.1 — Pen-en-papier

Vijf univariate toetsen, $p$-waardes gesorteerd: $.009, .011, .038, .087, .110$.

a) Bonferroni op $\alpha = .05$: welke zijn significant?

b) Holm-Bonferroni: loop de stappen één voor één door. Tot waar kom je?

c) Welk verschil ontstaat tussen Bonferroni en Holm op deze data?

OpgeletAntwoord T6.1 — open na je eigen poging

a) Bonferroni: vergelijk alle vijf met $.05 / 5 = .010$. Alleen $p_1 = .009 < .010$ ✓. De andere vier $> .010$ ✗. Eén test significant.

b) Holm-Bonferroni, sequentieel:

• Stap 1: $p_{(1)} = .009$ vs $\alpha_1 = .05 / 5 = .010$. Significant ✓. Door.
• Stap 2: $p_{(2)} = .011$ vs $\alpha_2 = .05 / 4 = .0125$. Significant ✓. Door.
• Stap 3: $p_{(3)} = .038$ vs $\alpha_3 = .05 / 3 \approx .0167$. Niet significant ✗. Stop. Ook $p_{(4)}, p_{(5)}$ blijven niet-significant.

Twee tests significant na Holm.

c) Holm vangt de tweede test ($p = .011$) die Bonferroni niet pakt — Bonferroni is uniform strenger op de eerste vergelijking, Holm is ruimer op latere stappen.

OpmerkingVraag T6.2 — Pen-en-papier (waar Holm écht verschil maakt)

Vijf univariate toetsen, $p$-waardes gesorteerd: $.009, .012, .016, .021, .049$.

a) Bonferroni op $\alpha = .05$: welke zijn significant?

b) Holm-Bonferroni: loop alle stappen door. Tot waar kom je?

OpgeletAntwoord T6.2 — open na je eigen poging

a) Bonferroni-drempel $.05 / 5 = .010$. Alleen $p_1 = .009$ haalt het. Eén significant.

b) Holm:

• $p_{(1)} = .009$ vs $.05/5 = .010$: sig ✓.
• $p_{(2)} = .012$ vs $.05/4 = .0125$: sig ✓.
• $p_{(3)} = .016$ vs $.05/3 \approx .0167$: sig ✓.
• $p_{(4)} = .021$ vs $.05/2 = .025$: sig ✓.
• $p_{(5)} = .049$ vs $.05/1 = .05$: sig ✓.

Alle vijf significant. Bonferroni vond er één. Holm vindt er vijf. Klassiek voorbeeld waar de sequentiële regel veel power terugwint.

Inzicht. Bonferroni is uniform conservatief (één drempel voor allemaal, simpel maar streng). Holm-Bonferroni is sequentieel slim (drempel verruimt onderweg, gegarandeerd minstens zoveel ontdekkingen als Bonferroni). Allebei controleren ze familywise error rate. Holm verliest minder power. Bij MANOVA-follow-up is één van beide de standaard — kies en blijf consistent.

In R: p.adjust(ps, method = "holm") past Holm toe en geeft gecorrigeerde $p$-waardes terug — vergelijk met de gewone $\alpha = .05$. Output gemakkelijker dan handmatig de drempel verschuiven.

5.2.b Tukey HSD per significante DV

TipHoe lees je het verschil-getal? Richting van $A - B$

In Tukey-output staat per groepenpaar een diff-kolom (het sample-verschil), een lwr/upr-CI, en een p adj. Belangrijke leesregel:

• Naam van de rij zegt welke twee groepen worden vergeleken, met een min-teken: bv. Catering-Sales betekent $\bar{Y}_{\text{Catering}} - \bar{Y}_{\text{Sales}}$.
• Teken van diff geeft de richting: als diff > 0, dan scoort de groep vóór het min-teken hoger. Als diff < 0, dan scoort de groep ná het min-teken hoger.
• 95%-CI bevat nul? → niet significant op $\alpha = .05$. CI helemaal positief óf helemaal negatief → significant, met de richting die het teken aangeeft.

Voorbeeld: Catering-Sales: diff = 1.17, lwr = 0.20, upr = 2.14, p adj = .013. Lees: Catering scoort gemiddeld $1.17$ punt hoger dan Sales; CI ligt helemaal boven nul; verschil is significant.

Dit is precies waar in een bijles soms een knoop ontstaat: vóór je een uitspraak doet, vertaal het min-teken bewust naar groepsnamen. “Het teken zegt wie hoger zit, het CI zegt of het verschil betrouwbaar is.” Diezelfde regel komt terug in MRA (dummy-codering), ANCOVA (adjusted means) en in alle thema’s met emmeans of pairs.

Algemene vorm.

# Voor elke significante DV: Tukey HSD om groepsverschillen te lokaliseren.
TukeyHSD(aov(Y1 ~ factor, data = mijn_data))
TukeyHSD(aov(Y2 ~ factor, data = mijn_data))
TukeyHSD(aov(Y3 ~ factor, data = mijn_data))

Voor onze dieren.

TukeyHSD(aov(trots       ~ sociale_status, data = leeuw_gevoelens))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = trots ~ sociale_status, data = leeuw_gevoelens)

$sociale_status
                diff       lwr       upr p adj
midden-alpha     -25 -29.92086 -20.07914     0
outsider-alpha   -55 -59.98691 -50.01309     0
outsider-midden  -30 -34.92086 -25.07914     0

TukeyHSD(aov(jaloezie    ~ sociale_status, data = leeuw_gevoelens))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = jaloezie ~ sociale_status, data = leeuw_gevoelens)

$sociale_status
                diff          lwr        upr     p adj
midden-alpha      35  30.13924668  39.860753 0.0000000
outsider-alpha     5   0.07399875   9.926001 0.0458218
outsider-midden  -30 -34.86075332 -25.139247 0.0000000

TukeyHSD(aov(eenzaamheid ~ sociale_status, data = leeuw_gevoelens))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = eenzaamheid ~ sociale_status, data = leeuw_gevoelens)

$sociale_status
                diff      lwr      upr p adj
midden-alpha      25 20.34728 29.65272     0
outsider-alpha    50 45.28483 54.71517     0
outsider-midden   25 20.34728 29.65272     0

TukeyHSD(aov(woede       ~ sociale_status, data = leeuw_gevoelens))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = woede ~ sociale_status, data = leeuw_gevoelens)

$sociale_status
                diff       lwr        upr p adj
midden-alpha      25  19.91027  30.089734     0
outsider-alpha   -15 -20.15806  -9.841944     0
outsider-midden  -40 -45.08973 -34.910266     0

OpmerkingVragen 5.2

d) Voor trots: tussen welke statusgroepen is het verschil significant?

e) Voor jaloezie: hetzelfde patroon of anders?

f) Voor eenzaamheid en woede: hetzelfde patroon of anders?

g) Welke groep is op alle vier de DV’s afwijkend van de andere twee?

OpgeletAntwoord 5.2 — open na je eigen poging

In gewone woorden.

4. Voor trots: alle drie de paarsgewijze verschillen significant — alpha > middenrang > outsider, met ruime gaps.

5. Voor jaloezie: ook alle drie paarsgewijze verschillen significant, in de volgorde alpha < outsider < middenrang — middenrangers het meest jaloers, alpha’s het minst, outsiders ertussenin. Het verschil tussen outsiders en alpha’s haalt de gewone $.05$-drempel van Tukey maar net ($p = .046$).

6. Voor eenzaamheid: outsider scoort significant hoger dan zowel alpha als middenrang; alpha en middenrang verschillen ook onderling, met outsider als duidelijkste afwijker. Voor woede: alle drie paarsgewijze verschillen significant, in de volgorde outsider < alpha < middenrang — middenrangers het meest woedend, outsiders het minst.

7. Geen enkele groep is op alle vier identiek aan de andere twee — alle drie hebben een eigen profiel: alpha (hoog trots), middenrang (hoog jaloezie en woede), outsider (hoog eenzaamheid). Dat is precies wat de DDA in de volgende stap zal bevestigen.

APA-stijl.

Tukey HSD-vergelijkingen lieten verschillende patronen per DV zien. Voor trots verschilden alle drie groepen onderling significant (alle $p < .001$), met alpha’s hoog en outsiders laag. Voor jaloezie verschilden alle drie groepen eveneens significant (middenrang vs. alpha en middenrang vs. outsider beide $p < .001$; outsider vs. alpha $p = .046$), met middenrangers het hoogst en alpha’s het laagst. Voor eenzaamheid scoorden outsiders significant hoger dan beide andere groepen (beide $p < .001$). Voor woede verschilden alle drie groepen onderling significant (alle $p < .001$), met middenrangers het hoogst en outsiders het laagst.

5.2.c Descriptive Discriminant Analysis (DDA)

OpmerkingWat is DDA en wanneer gebruik je het?

Descriptive Discriminant Analysis is een alternatieve follow-up: in plaats van DV-voor-DV te kijken, vraag je welke combinaties van DV’s de groepen het beste van elkaar onderscheiden. DDA produceert discriminant-functies — gewogen sommen van de DV’s, gekozen om groepsverschillen te maximaliseren.

Per discriminant-functie krijg je:

• Eigenvalue $\mathit{\lambda}_k$ — hoeveel variantie de functie tussen groepen verklaart.
• Canonical $R^2$ — analoog aan ANOVA-$R^2$ binnen de discriminant-dimensie.
• Standardized coefficients — de gewichten voor elke DV in de functie.
• Structure matrix — correlaties tussen elke DV en elke discriminant-functie. Standaard de hoofdmaat voor “welke DV draagt bij aan welke functie” (analoog aan factor-loadings in factor-analyse). Dit oefenboek legt de grens bij $.40$: een DV telt mee bij een functie als zijn correlatie ermee boven $.40$ of onder $-.40$ ligt.
• Group means on variates — waar elke statusgroep ligt op elke discriminant-as.

Aantal mogelijke discriminant-functies. Niet onbeperkt: $\min(I-1, p)$. Voor onze leeuw $\min(2, 4) = 2$.

Waarom dat minimum? Stel: je gaat lunchen met een vriend, en je spreekt af dat ieder evenveel uitgeeft. Jij hebt €10 op zak, je vriend €5. Hoeveel geeft ieder dan uit? Niet meer dan €5 — wat er samen kan, hangt af van het kleinste budget. Bij DDA werkt het hetzelfde: je hebt $I-1$ “groep-budgetten” (zoveel onafhankelijke richtingen waarin $I$ groepen uit elkaar kunnen vallen) en $p$ “DV-budgetten” (zoveel dimensies waarin je überhaupt kunt meten). Het aantal discriminant-functies kan nooit boven het kleinste van die twee uitkomen — dus $\min(I-1, p)$.

DDA versus PDA. Descriptive DA verklaart waar de groepsverschillen zitten. Predictive DA gebruikt de functies om gevallen toe te wijzen aan groepen — een classificatie-doel. Hier kijken we hóé de groepen op de functies verschillen.

OpmerkingDDA-structuur in één plaatje — gewichten $b_{jk}$ op de pijlen

Je kunt de DDA-structuur tekenen als path-diagram: $X$ (de groepen-IV, $k$ niveaus) links, $D_1$ en $D_2$ (twee discriminant-functies) in het midden, en $Y_1, Y_2, \ldots, Y_p$ (de DV’s) rechts. Op elke pijl van $Y_j$ naar $D_k$ staat het gewicht $b_{j,k}$ — eerste index = DV-positie, tweede index = D-functie (afspraak in dit oefenboek: positie-binnen-set eerst, set-niveau tweede).

Hoe lees je dit:

• De pijlen van $X$ naar $D_1$ en $D_2$ zijn de MANOVA-laag: het effect van de groepen op de discriminant-functies (analoog aan het effect van groepen op de DV’s bij MANOVA).
• De pijlen van $Y_j$ naar $D_k$ zijn de DDA-laag: elke DV draagt met een eigen gewicht $b_{j,k}$ bij aan elke functie. Voor functie $D_1$ zijn dat de blauwe pijlen $b_{1,1}, b_{2,1}, \ldots, b_{p,1}$; voor $D_2$ de rode $b_{1,2}, b_{2,2}, \ldots, b_{p,2}$.
• Het aantal discriminant-functies is $\min(I-1, p)$ — voor onze leeuw met $I = 3$ statusgroepen en $p = 4$ DV’s: twee functies, dus dit hele plaatje met $D_1$ en $D_2$ is letterlijk wat je krijgt.

Concreet voor de leeuw: $Y_1$ = trots, $Y_2$ = jaloezie, $Y_3$ = eenzaamheid, $Y_4$ = woede. Elke discriminant-functie is dan:

$$D_1 = b_{1,1} \cdot \text{trots} + b_{2,1} \cdot \text{jaloezie} + b_{3,1} \cdot \text{eenzaamheid} + b_{4,1} \cdot \text{woede}$$

$$D_2 = b_{1,2} \cdot \text{trots} + b_{2,2} \cdot \text{jaloezie} + b_{3,2} \cdot \text{eenzaamheid} + b_{4,2} \cdot \text{woede}$$

De gewichten $b_{j,k}$ kiezen we niet zelf — candisc() rekent ze uit zodat de groepsverschillen op $D_1$ maximaal zijn, en daarna op $D_2$ maximaal onder de voorwaarde dat $D_2$ onafhankelijk is van $D_1$. Dat is wat “discriminant-functie” betekent: een gewogen optelling, met optimaal gekozen gewichten.

TipNaN in de toets van een discriminant-functie — wat betekent dat?

Soms staat er in de toetstabel van candisc() bij een functie $\text{NaN}$ (“Not a Number”) in plaats van een $F$ en een $p$. Dat zegt niets over die functie. candisc() rekent de toets per rij om naar een $F$ met een benaderingsformule, en die formule deelt nul door nul zodra er in die rij precies twee vrijheidsgraden in de teller over zijn. Bij drie groepen en drie DV’s is dat de tweede rij.

Je ziet het zelf in E5 aan het eind van dit hoofdstuk: daar heeft de tweede functie een eigenvalue van $1.03$ en een canonical $R^2$ van $.51$, en tóch staat er NaN. Schrijf dan niet op dat functie 2 niet significant is. Wat je wél opschrijft: dat de toets voor die functie in deze uitvoer niet te berekenen is, en dat hij met die eigenvalue en die canonical $R^2$ een flink deel van de groepsverschillen draagt. Een functie die bijna niets verklaart, heeft een eigenvalue dicht bij $0$.

In onze leeuw-data, met vier DV’s, gebeurt het niet: daar krijgen beide functies een gewone toets.

OpmerkingVoor wie verder wil — de NaN tóch uitrekenen

Voor precies dit geval spreekt de benaderingsformule af dat je in plaats van nul-gedeeld-door-nul een $1$ neemt. Dan komt functie 2 in E5 uit op $F(2, 80) = 41.03$, $p < .001$. Met die afspraak is de toets in dit geval zelfs exact.

Algemene vorm.

# Eerst het lm-object hergebruiken; dan candisc() erop.
library(candisc)
DA <- candisc(lm_man)
DA                  # eigenvalues, canonical R^2, percentages
summary(DA)         # plus structure matrix en group means

Voor onze dieren.

library(candisc)

DA_leeuw <- candisc(lm_leeuw)
DA_leeuw

Canonical Discriminant Analysis for sociale_status:

   CanRsq Eigenvalue Difference Percent Cumulative
1 0.92104    11.6641     8.1499  76.847     76.847
2 0.77848     3.5142     8.1499  23.153    100.000

Test of H0: The canonical correlations in the 
current row and all that follow are zero

  LR test stat approx F numDF denDF   Pr(> F)    
1     0.017492   170.58     8   208 < 2.2e-16 ***
2     0.221525   123.00     3   105 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(DA_leeuw)

Canonical Discriminant Analysis for sociale_status:

   CanRsq Eigenvalue Difference Percent Cumulative
1 0.92104    11.6641     8.1499  76.847     76.847
2 0.77848     3.5142     8.1499  23.153    100.000

Class means:

            Can1     Can2
alpha    -2.9918  2.08083
midden   -1.6823 -2.37152
outsider  4.7676  0.42244

 std coefficients:
                Can1     Can2
trots       -0.64696  0.11081
jaloezie    -0.34631 -0.65883
eenzaamheid  0.58033 -0.19557
woede       -0.31406 -0.34541

# De structure matrix expliciet — DV's tegen functies.
round(DA_leeuw$structure, 3)

              Can1   Can2
trots       -0.921  0.330
jaloezie    -0.217 -0.962
eenzaamheid  0.898 -0.380
woede       -0.620 -0.731

plot(DA_leeuw,
     col   = c("#0077BB", "#33BBEE", "#EE7733"),
     pch   = 16,
     ellipse = TRUE,
     scale = 6,
     main = "DDA — leeuwen op de twee discriminant-functies")

OpmerkingVragen 5.2

h) Hoeveel discriminant-functies zijn er bij $I = 3$ groepen en $p = 4$ DV’s?

i) Wat zijn de eigenvalues van Can1 en Can2? Welk percentage van de groepsvariantie verklaren ze samen?

j) Volgens de structure matrix: welke DV correleert het sterkst met Can1? Met Can2? Wat zegt dat?

k) Hoe liggen de drie statusgroepen op Can1 en Can2? Welke groep ligt het verst weg van de anderen op Can1?

OpgeletAntwoord 5.2 — open na je eigen poging

In gewone woorden.

8. Aantal discriminant-functies = $\min(I-1, p) = \min(2, 4) = 2$ — twee functies.

9. Eigenvalue Can1 $= 11.66$ ($76.8\%$ van variantie); Can2 $= 3.51$ ($23.2\%$). Samen $100\%$ — definitief, omdat er maar twee zijn.

10. Structure matrix: Can1 wordt het sterkst gedragen door trots ($r = -.92$) en eenzaamheid ($r = +.90$) — Can1 is de isolatie-as: hoge waarden = hoge eenzaamheid + lage trots. Can2 wordt het sterkst gedragen door jaloezie ($r = -.96$) en woede ($r = -.73$) — Can2 is de rust-van-wedijver-as: lage waarden = hoge jaloezie + woede. Met de grens uit het kader hierboven (een correlatie boven $.40$ of onder $-.40$) tellen op Can1 trots, eenzaamheid en woede mee, en op Can2 jaloezie en woede. Eén DV telt op allebei mee: woede.

11. Group means: alpha’s liggen laag op Can1 (Can1 $= -2.99$ — hoge trots, lage eenzaamheid) en hoog op Can2 (Can2 $= +2.08$ — weinig wedijver). Outsiders liggen hoog op Can1 (Can1 $= +4.77$ — lage trots, hoge eenzaamheid). Middenrangers liggen laag op Can2 (Can2 $= -2.37$ — veel jaloezie + woede). Op Can1 is alpha versus outsider de grootste afstand (monotone rang-as); op Can2 is midden versus alpha-outsider de grootste afstand (wedijver-piek bij midden).

APA-stijl.

Descriptive discriminant analysis bevestigde dat de vier DV’s twee onderscheidende dimensies vormen. Functie 1 (eigenvalue $= 11.66$, $76.8\%$ verklaarde groepsvariantie) werd vooral gedragen door trots ($r = -.92$) en eenzaamheid ($r = +.90$) en separeerde alpha’s (laag) van outsiders (hoog) — de isolatie-as. Functie 2 (eigenvalue $= 3.51$, $23.2\%$) werd vooral gedragen door jaloezie ($r = -.96$) en woede ($r = -.73$) en separeerde middenrangers (laag) van alpha’s en outsiders (hoog) — de wedijver-piek-as.

T7 — Eigenvalues en structure matrix in DDA

OpmerkingVraag T7 — Pen-en-papier

a) Hoeveel discriminant-functies zijn maximaal mogelijk bij $I = 4$ groepen en $p = 6$ DV’s?

b) Stel: eigenvalues $\lambda_1 = 0.80$, $\lambda_2 = 0.10$, $\lambda_3 = 0.02$. Welk percentage variantie verklaart functie 1?

c) Wat zegt de structure matrix dat de standardized coefficients niet zeggen?

d) Een DV heeft een hoge structure-correlatie ($r = .85$) maar een lage standardized coefficient ($b = 0.10$). Wat suggereert dat?

OpgeletAntwoord T7 — open na je eigen poging

a) $\min(I-1, p) = \min(3, 6) = 3$ functies.

b) $\mathit{\lambda}_1 / \sum \mathit{\lambda} = 0.80 / 0.92 = .870$, oftewel $87\%$ — functie 1 is overheersend.

c) Structure matrix geeft de bivariate samenhang tussen elke DV en elke discriminant-functie — analoog aan factor-loadings in EFA. Geeft het totale belang van een DV voor groepsscheiding. Standardized coefficients geven het unieke gewicht in de discriminant-functie, na controle voor andere DV’s. Bij multicollineaire DV’s kunnen ze sterk verschillen — coëfficiënten lopen leeg, structure-correlaties blijven hoog.

d) Multicollineariteit met andere DV’s. De DV draagt op zich veel informatie (hoge structure $r$), maar verliest dat unieke gewicht in de coefficient zodra de andere DV’s er ook zijn. Een vergelijkbare situatie als in MRA bij multicollineaire voorspellers (zie thema 1, T3 — Multicollineariteit zonder R). De structure-correlatie is dan de meer betrouwbare interpretatie-bron.

Inzicht. Voor het interpreteren van een discriminant-functie (“welke DV’s dragen bij?”) kies je standaard de structure matrix, niet de gestandaardiseerde coëfficiënten. Bij niet-multicollineaire DV’s geven beide hetzelfde antwoord; bij collineariteit geeft de structure matrix het stabielere beeld.

T8 — MANOVA versus repeated-measures ANOVA

OpmerkingVraag T8 — Pen-en-papier

In thema 6 komt repeated-measures ANOVA RMA voor meerdere metingen op dezelfde proefpersoon (bv. drie tijdstippen). MANOVA werkt ook met meerdere uitkomsten — wanneer is het wel/geen RMA?

a) Bij drie tijdstippen op dezelfde dieren: MANOVA of RMA?

b) Bij drie verschillende constructen op dezelfde dieren (trots, jaloezie, eenzaamheid): MANOVA of RMA?

OpgeletAntwoord T8 — open na je eigen poging

a) Drie tijdstippen op dezelfde dieren: in principe beide kunnen, maar RMA is het traditionele werkpaard — het exploiteert dat de drie metingen hetzelfde construct meten en behandelt tijd als within-subjects factor.

b) Drie verschillende constructen: MANOVA. De drie DV’s meten verschillende dingen die niet op één schaal vergelijkbaar zijn. Een tijds-as bestaat hier niet.

Inzicht. MANOVA is breed (verschillende DV’s), RMA is gespecialiseerd (zelfde DV op verschillende tijden).

5.2.d Voorzet — twee oogstjaren tulpen

“De tulp groeit niet op één maat alleen — steel en kelk samen bepalen welke soort het is.”

“De leeuw werd nieuwsgierig naar de tulpenkweker. Drie soorten in de bloembollen, twee maten per tulp. De ene oogst pakte MANOVA wél, de andere niet. Hoe kon dat? Hetzelfde aantal tulpen, dezelfde gemiddelden — en toch een ander verhaal.”

In een tulpenkwekerij worden drie soorten gecultiveerd: Tulipa univariegata, Tulipa multivariegata en Tulipa variabilis. Twee oogstjaren beschikbaar: 2023 en 2024. Per tulp twee maten: steellengte (cm) en kelkdiameter (mm). Per soort zijn er $30$ tulpen per oogstjaar, totaal $N = 90$ per jaar.

De kwekerij vraagt: “Verschillen onze drie soorten in hun bouw? En zo ja, op welke combinatie van maten zit het onderscheid?”

Verkenning

Base R

# Gemiddelde en SD per soort, eerst oogst 2023, dan oogst 2024.
aggregate(cbind(steellengte, kelkdiameter) ~ soort, data = tulpen_2023, FUN = mean)

           soort steellengte kelkdiameter
1   univariegata        39.6         66.4
2 multivariegata        41.4         65.0
3     variabilis        42.0         62.6

aggregate(cbind(steellengte, kelkdiameter) ~ soort, data = tulpen_2023, FUN = sd)

           soort steellengte kelkdiameter
1   univariegata           6            8
2 multivariegata           6            8
3     variabilis           6            8

aggregate(cbind(steellengte, kelkdiameter) ~ soort, data = tulpen_2024, FUN = mean)

           soort steellengte kelkdiameter
1   univariegata        39.6         66.4
2 multivariegata        41.4         65.0
3     variabilis        42.0         62.6

aggregate(cbind(steellengte, kelkdiameter) ~ soort, data = tulpen_2024, FUN = sd)

           soort steellengte kelkdiameter
1   univariegata           6            8
2 multivariegata           6            8
3     variabilis           6            8

dplyr

# Per soort samenvatten: gemiddelde en SD van beide maten, oogst 2023.
tulpen_2023 |>
  group_by(soort) |>
  summarise(across(c(steellengte, kelkdiameter), mean))

# A tibble: 3 × 3
  soort          steellengte kelkdiameter
  <fct>                <dbl>        <dbl>
1 univariegata          39.6         66.4
2 multivariegata        41.4         65  
3 variabilis            42           62.6

tulpen_2023 |>
  group_by(soort) |>
  summarise(across(c(steellengte, kelkdiameter), sd))

# A tibble: 3 × 3
  soort          steellengte kelkdiameter
  <fct>                <dbl>        <dbl>
1 univariegata             6            8
2 multivariegata           6            8
3 variabilis               6            8

# Hetzelfde voor oogst 2024.
tulpen_2024 |>
  group_by(soort) |>
  summarise(across(c(steellengte, kelkdiameter), mean))

# A tibble: 3 × 3
  soort          steellengte kelkdiameter
  <fct>                <dbl>        <dbl>
1 univariegata          39.6         66.4
2 multivariegata        41.4         65  
3 variabilis            42           62.6

tulpen_2024 |>
  group_by(soort) |>
  summarise(across(c(steellengte, kelkdiameter), sd))

# A tibble: 3 × 3
  soort          steellengte kelkdiameter
  <fct>                <dbl>        <dbl>
1 univariegata             6            8
2 multivariegata           6            8
3 variabilis               6            8

OpmerkingVragen 5.2.d — verkenning

a) Hoeveel tulpen per soort? Hoeveel totaal per oogstjaar?

b) Welke soort heeft de langste gemiddelde steel in oogst 2023? In oogst 2024? Verschillen de groepsgemiddelden voor steellengte en kelkdiameter tussen de twee jaren?

c) Vergelijk de SD’s per soort over de twee jaren. Wat valt op?

OpgeletAntwoord 5.2.d — open na je eigen poging

a) $30$ tulpen per soort, $N = 90$ per oogstjaar.

b) In beide jaren heeft variabilis de langste steel ($M = 42.00$ cm) en univariegata de kortste ($M = 39.60$ cm). Bij de kelk is het andersom: univariegata het wijdst ($M = 66.40$ mm), variabilis het smalst ($M = 62.60$ mm). Tussen de jaren verschillen de gemiddelden niet: ze zijn per soort precies gelijk. Tussen de soorten verschillen ze maar een beetje — hooguit 2.4 cm in steellengte en 3.8 mm in kelkdiameter.

c) Ook de SD’s zijn gelijk: in elke soort en in elk jaar $SD = 6.00$ cm voor steellengte en $SD = 8.00$ mm voor kelkdiameter. Vergeleken met die spreiding zijn de verschillen tussen de soorten uit (b) klein. Deze twee oogsten zijn met opzet zo gemaakt: alles wat je per maat los kunt bekijken, is in beide jaren hetzelfde. Wat er tussen de jaren wél verschilt, zie je pas als je de twee maten samen bekijkt.

Steellengte en kelkdiameter samen

beide <- rbind(cbind(oogst = "Oogst 2023", tulpen_2023),
               cbind(oogst = "Oogst 2024", tulpen_2024))
ggplot(beide, aes(x = steellengte, y = kelkdiameter,
                  colour = soort, fill = soort, shape = soort)) +
  geom_point(size = 1.8, alpha = 0.7) +
  stat_ellipse(geom = "polygon", alpha = 0.15, level = 0.68) +
  facet_wrap(~ oogst) +
  scale_colour_manual(values = c("#0077BB", "#33BBEE", "#EE7733")) +
  scale_fill_manual(values   = c("#0077BB", "#33BBEE", "#EE7733")) +
  labs(x = "Steellengte (cm)", y = "Kelkdiameter (mm)",
       colour = "Soort", fill = "Soort", shape = "Soort") +
  theme_minimal() +
  theme(legend.position = "bottom")

Steellengte tegen kelkdiameter, per soort en per oogst. Elke stip is een tulp; elke ellips omsluit ongeveer 68% van de tulpen van die soort.

Base R

# Samenhang tussen steellengte en kelkdiameter, binnen elke soort apart.
# split() knipt de data per soort; sapply() rekent op elk stuk de correlatie.
sapply(split(tulpen_2023, tulpen_2023$soort),
       function(d) cor(d$steellengte, d$kelkdiameter))

  univariegata multivariegata     variabilis 
     0.8942223      0.7896318      0.8821254 

sapply(split(tulpen_2024, tulpen_2024$soort),
       function(d) cor(d$steellengte, d$kelkdiameter))

  univariegata multivariegata     variabilis 
    -0.7318754     -0.7924560     -0.7524337 

dplyr

# Dezelfde correlaties: per soort één regel, met r erin.
tulpen_2023 |>
  group_by(soort) |>
  summarise(r = cor(steellengte, kelkdiameter))

# A tibble: 3 × 2
  soort              r
  <fct>          <dbl>
1 univariegata   0.894
2 multivariegata 0.790
3 variabilis     0.882

tulpen_2024 |>
  group_by(soort) |>
  summarise(r = cor(steellengte, kelkdiameter))

# A tibble: 3 × 2
  soort               r
  <fct>           <dbl>
1 univariegata   -0.732
2 multivariegata -0.792
3 variabilis     -0.752

Een tibble drukt korter af dan base R: 0.894 waar base R 0.8942223 laat zien. Het getal eronder is hetzelfde.

OpmerkingVragen 5.2.d — samen bekeken

d) Vergelijk de richting van de ellipsen in oogst 2023 en oogst 2024. Wat zie je binnen de soorten? Wat zeggen de correlaties binnen elke soort (de uitvoer onder het plaatje)?

e) In welk jaar liggen de drie soorten het duidelijkst uit elkaar? Waarom, als hun gemiddelden in beide jaren dezelfde zijn?

OpgeletAntwoord 5.2.d — open na je eigen poging

d) In 2023 lopen de ellipsen van linksonder naar rechtsboven: binnen elke soort heeft een tulp met een lange steel meestal ook een wijde kelk ($r$ tussen $.79$ en $.89$). In 2024 lopen ze van linksboven naar rechtsonder: een lange steel gaat samen met een smalle kelk ($r$ tussen $-.79$ en $-.73$). Binnen één jaar hebben de drie soorten ongeveer dezelfde samenhang; tussen de jaren draait hij om.

e) In 2023. De soorten liggen in beide jaren op dezelfde plek: univariegata linksboven (korte steel, wijde kelk), variabilis rechtsonder, multivariegata ertussen. Van univariegata naar variabilis ga je dus schuin naar rechtsonder. In 2023 is dat dwars op de ellipsen: elke soort is een smalle, langgerekte wolk, en de drie wolken liggen naast elkaar. In 2024 is dat in de lengte van de ellipsen: de drie wolken schuiven grotendeels over elkaar heen.

MANOVA + univariate ANOVA’s per oogstjaar

Net als bij de leeuw: sla in de uitvoer het blok Term: (Intercept) over. Jouw blok is Term: soort, met daaronder de vier toetsen.

cat("=== MANOVA — oogst 2023 ===\n")

=== MANOVA — oogst 2023 ===

lm_tulp_23 <- lm(cbind(steellengte, kelkdiameter) ~ soort, data = tulpen_2023)
man_tulp_23 <- car::Manova(lm_tulp_23, type = 3)
print(summary(man_tulp_23, multivariate = TRUE, univariate = FALSE))

Type III MANOVA Tests:

Sum of squares and products for error:
             steellengte kelkdiameter
steellengte     3132.000     3571.843
kelkdiameter    3571.843     5568.000

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

Sum of squares and products for the hypothesis:
             steellengte kelkdiameter
steellengte       151290       238620
kelkdiameter      238620       376360

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1   0.98543 2907.591      2     86 < 2.22e-16 ***
Wilks             1   0.01457 2907.591      2     86 < 2.22e-16 ***
Hotelling-Lawley  1  67.61840 2907.591      2     86 < 2.22e-16 ***
Roy               1  67.61840 2907.591      2     86 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: soort 

Sum of squares and products for the hypothesis:
             steellengte kelkdiameter
steellengte         93.6       -130.8
kelkdiameter      -130.8        221.6

Multivariate Tests: soort
                 Df test stat  approx F num Df den Df         Pr(>F)    
Pillai            2 0.3156025  8.150517      4    174 0.000004777860 ***
Wilks             2 0.6849291  8.957190      4    172 0.000001355603 ***
Hotelling-Lawley  2 0.4592290  9.758617      4    170 0.000000395203 ***
Roy               2 0.4575326 19.902670      2     87 0.000000076314 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Univariate ANOVA's — oogst 2023 ---\n")

--- Univariate ANOVA's — oogst 2023 ---

print(summary.aov(lm_tulp_23))

 Response steellengte :
            Df Sum Sq Mean Sq F value Pr(>F)
soort        2   93.6    46.8     1.3 0.2778
Residuals   87 3132.0    36.0               

 Response kelkdiameter :
            Df Sum Sq Mean Sq F value Pr(>F)
soort        2  221.6   110.8  1.7312 0.1831
Residuals   87 5568.0    64.0               

cat("=== MANOVA — oogst 2024 ===\n")

=== MANOVA — oogst 2024 ===

lm_tulp_24 <- lm(cbind(steellengte, kelkdiameter) ~ soort, data = tulpen_2024)
man_tulp_24 <- car::Manova(lm_tulp_24, type = 3)
print(summary(man_tulp_24, multivariate = TRUE, univariate = FALSE))

Type III MANOVA Tests:

Sum of squares and products for error:
             steellengte kelkdiameter
steellengte     3132.000    -3169.257
kelkdiameter   -3169.257     5568.000

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

Sum of squares and products for the hypothesis:
             steellengte kelkdiameter
steellengte       151290       238620
kelkdiameter      238620       376360

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1    0.9979 20547.78      2     86 < 2.22e-16 ***
Wilks             1    0.0021 20547.78      2     86 < 2.22e-16 ***
Hotelling-Lawley  1  477.8554 20547.78      2     86 < 2.22e-16 ***
Roy               1  477.8554 20547.78      2     86 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: soort 

Sum of squares and products for the hypothesis:
             steellengte kelkdiameter
steellengte         93.6       -130.8
kelkdiameter      -130.8        221.6

Multivariate Tests: soort
                 Df test stat approx F num Df den Df  Pr(>F)
Pillai            2 0.0505365 1.127664      4    174 0.34509
Wilks             2 0.9499302 1.118691      4    172 0.34935
Hotelling-Lawley  2 0.0522177 1.109625      4    170 0.35369
Roy               2 0.0399060 1.735912      2     87 0.18229

cat("\n--- Univariate ANOVA's — oogst 2024 ---\n")

--- Univariate ANOVA's — oogst 2024 ---

print(summary.aov(lm_tulp_24))

 Response steellengte :
            Df Sum Sq Mean Sq F value Pr(>F)
soort        2   93.6    46.8     1.3 0.2778
Residuals   87 3132.0    36.0               

 Response kelkdiameter :
            Df Sum Sq Mean Sq F value Pr(>F)
soort        2  221.6   110.8  1.7312 0.1831
Residuals   87 5568.0    64.0               

OpmerkingVragen 5.2.d — MANOVA + univariate

f) Wat zijn Pillai’s $V$, de $F$-benadering en de $p$-waarde voor oogst 2023 en voor oogst 2024?

g) Welke univariate $F$’s per DV zie je? Zijn die significant na Bonferroni-correctie (drempel $.05/2 = .025$)?

h) In welk oogstjaar is het MANOVA-effect duidelijker? Komt dat overeen met wat je in de scatterplots zag (vraag e)?

i) Reflectie: in welk jaar kom je met alleen de losse ANOVA’s tot een andere conclusie dan met de MANOVA? Waarom?

OpgeletAntwoord 5.2.d — open na je eigen poging

f) Oogst 2023: Pillai’s $V = 0.32$, $F(4, 174) = 8.15$, $p < .001$. Oogst 2024: Pillai’s $V = 0.05$, $F(4, 174) = 1.13$, $p = .35$. In 2023 verschillen de soorten multivariaat, in 2024 niet.

g) In beide jaren precies dezelfde: steellengte $F(2, 87) = 1.30$, $p = .28$; kelkdiameter $F(2, 87) = 1.73$, $p = .18$. Geen van beide is significant, ook niet op de gewone $.05$, en dus zeker niet op de Bonferroni-drempel $.025$. (Volgens de route van dit hoofdstuk kijk je in 2024 niet eens naar de losse ANOVA’s, want daar was de MANOVA niet significant. Hier doen we het toch, om de twee jaren naast elkaar te leggen.)

Het verschil tussen de jaren staat zelfs in de MANOVA-uitvoer, in de twee blokken getallen boven de toetsen. Leg ze voor 2023 en 2024 naast elkaar. Sum of squares and products for the hypothesis, het verschil tussen de soorten, is in beide jaren precies hetzelfde. Sum of squares and products for error, de spreiding binnen de soorten, heeft ook dezelfde getallen op de diagonaal. Alleen het getal erbuiten verschilt: $3571.843$ in 2023 en $-3169.257$ in 2024. Dat getal hoort bij de samenhang tussen steel en kelk binnen de soorten — en dat ene minteken is het hele verschil tussen de twee oogsten.

h) In 2023, en dat klopt met de scatterplots: daar liggen de soorten als smalle wolken naast elkaar, in 2024 schuiven ze over elkaar heen.

i) In 2023. Wie alleen de losse ANOVA’s draait, concludeert in beide jaren dat de soorten niet verschillen. In 2023 is dat fout. De soorten verschillen wél, alleen niet op één maat apart: ze verschillen in de verhouding tussen de twee.

In 2023 horen binnen een soort een lange steel en een wijde kelk bij elkaar. Een tulp met een korte steel en tóch een wijde kelk valt dan op — en zo’n tulp is meestal een univariegata. Net als bij mensen: dat een langere groep ook zwaarder is, verrast niemand, maar een groep die langer is en tóch lichter, valt meteen op.

MANOVA maakt van die verhouding één getal per tulp, een gewogen optelling. Tel de kelk mee en trek de steel ervan af. Een grote tulp heeft een lange steel én een wijde kelk, en bij het aftrekken valt dat tegen elkaar weg. Wat overblijft, is alleen de verhouding — en daarin verschillen de soorten.

In 2024 kan dat niet. Daar gaat een korte steel binnen elke soort al samen met een wijde kelk, dus een univariegata met die bouw valt niet op tussen de andere twee.

De figuur laat zien wat elke toets ziet. In de eerste twee kolommen staan de soorten op één maat apart; dat is wat een losse ANOVA bekijkt, en daar liggen de drie gemiddelden (de zwarte ruiten) in beide jaren dicht bij elkaar. In de derde kolom staan ze op de beste gewogen optelling van dat jaar (de eerste discriminantfunctie, Can1, die je hieronder voor 2023 uitrekent): in 2023 schuiven de soorten daar duidelijk uit elkaar, in 2024 blijven ze op een kluitje. De beste optelling van 2024 is bovendien vrijwel alleen de kelkdiameter (gestandaardiseerde gewichten $-0.10$ voor de steel en $0.93$ voor de kelk): in dat jaar valt er niets te winnen door de twee maten te combineren, en de derde kolom lijkt daar dan ook op de tweede.

# Per jaar: de scores op Can1 uit een eigen discriminantanalyse. Het teken van
# zo'n functie is willekeurig; we kiezen het zo dat een wijde kelk rechts ligt,
# dan staan de soorten in beide jaren in dezelfde volgorde.
blik <- function(d, jaar) {
  can1 <- candisc::candisc(lm(cbind(steellengte, kelkdiameter) ~ soort, data = d))$scores$Can1
  can1 <- can1 * sign(cor(can1, d$kelkdiameter))
  data.frame(
    oogst = paste("Oogst", jaar),
    soort = rep(d$soort, 3),
    blik  = factor(rep(c("Alleen steellengte", "Alleen kelkdiameter",
                         "Gewogen optelling (Can1)"), each = nrow(d)),
                   levels = c("Alleen steellengte", "Alleen kelkdiameter",
                              "Gewogen optelling (Can1)")),
    score = c(scale(d$steellengte), scale(d$kelkdiameter), scale(can1))
  )
}
schaduw <- rbind(blik(tulpen_2023, 2023), blik(tulpen_2024, 2024))

# Eén rij stippen per soort; de zwarte ruit is het gemiddelde van die soort.
set.seed(5)
ggplot(schaduw, aes(x = score, y = soort, colour = soort)) +
  geom_jitter(height = 0.15, width = 0, size = 1.2, alpha = 0.6) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, colour = "black") +
  facet_grid(oogst ~ blik) +
  scale_colour_manual(values = c("#0077BB", "#33BBEE", "#EE7733"), guide = "none") +
  scale_y_discrete(limits = rev(levels(schaduw$soort))) +
  labs(x = "Gestandaardiseerde score", y = NULL) +
  theme_minimal() +
  theme(panel.grid.minor = element_blank())

Dezelfde tulpen, drie keer bekeken: per maat apart en op de beste gewogen optelling van dat jaar. Elke stip is een tulp, de zwarte ruit het gemiddelde van die soort. De as staat in standaardafwijkingen: 0 is het gemiddelde van alle tulpen van dat jaar.

Inzicht. Dit is het tweede geval uit kader A, Wat is MANOVA?: elke maat verschilt een beetje, te weinig om in zijn eentje op te vallen, en een goed gekozen combinatie strept de ruis weg die de maten delen. Of dat lukt, hangt niet af van de gemiddelden — die zijn in beide jaren gelijk — maar van hoe de maten binnen de soorten samenhangen.

DDA — welke gewogen optelling onderscheidt de soorten?

library(candisc)
cat("=== DDA — oogst 2023 ===\n")

=== DDA — oogst 2023 ===

DA_tulp_23 <- candisc(lm_tulp_23)
print(DA_tulp_23)

Canonical Discriminant Analysis for soort:

     CanRsq Eigenvalue Difference Percent Cumulative
1 0.3139090  0.4575326    0.45584 99.6306     99.631
2 0.0016935  0.0016964    0.45584  0.3694    100.000

Test of H0: The canonical correlations in the 
current row and all that follow are zero

  LR test stat approx F numDF denDF     Pr(> F)    
1      0.68493   8.9572     4   172 0.000001356 ***
2      0.99831   0.1476     1    87      0.7018    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Structure matrix ---\n")

--- Structure matrix ---

print(round(DA_tulp_23$structure, 3))

               Can1   Can2
steellengte  -0.296 -0.955
kelkdiameter  0.342 -0.940

cat("\n--- Gestandaardiseerde gewichten ---\n")

--- Gestandaardiseerde gewichten ---

print(round(DA_tulp_23$coeffs.std, 3))

               Can1   Can2
steellengte  -1.848 -0.558
kelkdiameter  1.870 -0.480

cat("\n--- Group means on variates ---\n")

--- Group means on variates ---

print(round(DA_tulp_23$means, 3))

                 Can1   Can2
univariegata    0.836  0.026
multivariegata -0.045 -0.057
variabilis     -0.791  0.031

OpmerkingVragen 5.2.d — DDA

j) Hoeveel discriminant-functies bestaan er voor deze data? ($\min(I-1, p)$ met $I = 3$ en $p = 2$)

k) Welke DV laadt het sterkst op Can1, en haalt hij de grens van $.40$? Welke DV laadt het sterkst op Can2?

l) Hoe liggen de drie tulpensoorten op Can1 en Can2? Beschrijf het patroon in eigen woorden.

OpgeletAntwoord 5.2.d — open na je eigen poging

j) $\min(I - 1, p) = \min(2, 2) = 2$ functies. Can1 draagt bijna alles: eigenvalue $0.46$, $99.6\%$ van het verschil tussen de soorten. Can2 heeft een eigenvalue van $0.002$ en doet vrijwel niets.

k) Op Can1 is dat kelkdiameter, met $r = .34$; steellengte heeft $r = -.30$. Geen van beide komt boven de $.40$ (of onder de $-.40$), dus met de vuistregel van $.40$ telt hier géén maat mee. Dat is geen fout van jou en geen fout in de uitvoer.

Wat is Can1 dan? Kijk naar de gestandaardiseerde gewichten: kelkdiameter $1.87$, steellengte $-1.85$. Tegengestelde tekens: Can1 telt de kelk en trekt de steel ervan af. Het is de functie wijde kelk voor zijn steellengte — dezelfde verhouding als in antwoord (i).

Waarom zijn de structure-correlaties dan zo laag? Een structure-correlatie zegt hoe goed je de plek van een tulp op de functie kent als je alleen die ene maat weet. Hier: slecht. Een tulp met een lange steel heeft meestal ook een wijde kelk, en dan valt het verschil tussen die twee — zijn plek op Can1 — gewoon in het midden. Eén maat vertelt je dus weinig; pas de twee samen vertellen het.

T7 zegt dat je een functie standaard uit de structure matrix leest. Dat werkt als de functie maten bij elkaar optelt. Is de functie een verschil tussen twee maten die binnen de groepen sterk samenhangen, zoals hier, dan blijven de structure-correlaties laag en vertellen de gewichten wat de functie is.

Op Can2 hebben beide gewichten hetzelfde teken ($-0.56$ en $-0.48$), en beide maten laden sterk ($r = -.96$ en $r = -.94$): Can2 telt steel en kelk bij elkaar op. Dat is de grootte van de tulp, en daarop verschillen de soorten nauwelijks.

l) Op Can1 ligt univariegata aan de ene kant ($0.84$), variabilis aan de andere ($-0.79$) en multivariegata ertussen ($-0.05$). In gewone woorden: univariegata heeft voor zijn steellengte de wijdste kelk, variabilis de smalste. Op Can2 liggen alle drie bij nul: groter of kleiner zijn ze niet.

TipTake-home van 5.2.d

1. Multivariate analyse is niet “univariate × meerdere” — het pakt iets dat losse ANOVA’s missen: de informatie in combinaties van DV’s.
2. De samenhang tussen de DV’s binnen de groepen bepaalt of MANOVA iets ziet wat de losse ANOVA’s missen. Loopt het verschil tussen de groepen dwars op die samenhang, dan ziet MANOVA het scherp; loopt het er in de lengte langs, dan verdrinkt het in de gewone variatie.
3. Zonder verschil in gemiddelden geen multivariaat verschil. Hebben de soorten binnen één oogst exact dezelfde gemiddelden op beide maten, dan vindt geen enkele multivariate toets iets — hoe de samenhang binnen de soorten er ook uitziet. Die samenhang bepaalt hoe scherp MANOVA een verschil in gemiddelden ziet, niet óf er een verschil is. Daarom zijn de gemiddelden hier gelijk tussen de jaren, niet tussen de soorten.
4. DDA = lineaire combinatie waarvan de gewichten door de data zijn gekozen om groepen maximaal te scheiden. Bij twee DV’s en drie groepen krijg je twee discriminant-functies (theoretisch maximum) — maar zoals hier kan één ervan bijna al het werk doen.

5.3 Factoriële MANOVA — twee factoren tegelijk

“Wat als ik ook geslacht meeneem?”

“Misschien,” zei de leeuw, “verschilt het gevoels-profiel niet alleen tussen statusgroepen, maar ook tussen mannetjes en vrouwtjes — en is het effect van status anders bij mannen dan bij vrouwen. Een interactie-effect, multivariate uitgevoerd.”

MANOVA is uit te breiden met meerdere factoren + interactie-termen — net als ANOVA in thema 2. Bij factoriële MANOVA kun je hoofdeffecten en interactie-effecten allemaal multivariaat toetsen.

Algemene vorm.

# Twee factoren + interactie, p DV's.
lm_facman <- lm(cbind(Y1, Y2, Y3, Y4) ~ factor1 * factor2, data = mijn_data)
res <- car::Manova(lm_facman, type = 3)
summary(res, multivariate = TRUE, univariate = FALSE)

# Univariate F's per DV: niet univariate = TRUE (zie de waarschuwing in 5.1.c),
# en hier ook niet summary.aov() -- die geeft type-I SS, en bij een
# ongebalanceerd factorieel design is dat iets anders dan type III.
car::Anova(lm(Y1 ~ factor1 * factor2, data = mijn_data), type = 3)
car::Anova(lm(Y2 ~ factor1 * factor2, data = mijn_data), type = 3)
# ... en zo verder voor Y3 en Y4.

OpmerkingWat verandert er bij factoriële MANOVA?

• Drie multivariate $F$-toetsen: voor factor1, factor2 en factor1:factor2 apart.
• Bij significante interactie: het hoofd-effect heeft minder betekenis en mag niet los geïnterpreteerd worden — qualified by interaction. Splits dan de analyse per niveau van de andere factor en bekijk binnen elk subset het multivariate effect afzonderlijk.
• Univariate follow-up per term: per significante term en per DV, met Bonferroni-correctie over het aantal DV’s binnen die term.
• Type-III SS wordt belangrijk bij onbalans + interactie. Vandaar Manova(., type = 3) met contr.sum-contrasten.
• Aannames: net als eenwegs-MANOVA, maar Box’s M wordt vaak per cel-combinatie problematisch.

R-spiekblad bij MANOVA

Alle commando’s op één plek

Pakketten en data laden

# Pakketten activeren — één keer per sessie.
library(car)         # Manova(), type-III SS
library(rstatix)     # box_m()
library(candisc)     # candisc() voor DDA

# Hoofd-dataset laden.
load("data/leeuw_gevoelens.RData")
str(leeuw_gevoelens)

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

Verkenning — DV-means en correlatiematrix

Base R

# Gemiddelden per groep voor elke DV.
aggregate(cbind(trots, jaloezie, eenzaamheid, woede) ~ sociale_status,
          data = leeuw_gevoelens, FUN = mean)

# Correlatiematrix tussen de DV's.
cor(leeuw_gevoelens[, c("trots", "jaloezie", "eenzaamheid", "woede")])

# Boxplot per DV.
boxplot(trots       ~ sociale_status, data = leeuw_gevoelens)
boxplot(jaloezie    ~ sociale_status, data = leeuw_gevoelens)
boxplot(eenzaamheid ~ sociale_status, data = leeuw_gevoelens)
boxplot(woede       ~ sociale_status, data = leeuw_gevoelens)

dplyr

# Gemiddelden per groep voor elke DV.
leeuw_gevoelens |>
  group_by(sociale_status) |>
  summarise(across(c(trots, jaloezie, eenzaamheid, woede), mean))

# Correlatiematrix tussen de DV's.
leeuw_gevoelens |>
  select(trots, jaloezie, eenzaamheid, woede) |>
  cor()

# De boxplots blijven base R: dat is tekenen, geen data bewerken.

Box’s M — gelijke covariantie-matrices

box_m(data  = leeuw_gevoelens[, c("trots", "jaloezie", "eenzaamheid", "woede")],
      group = leeuw_gevoelens$sociale_status)

# Beoordeel op alpha = .001 (Box's M is overgevoelig op gewone alpha).

MANOVA fitten — lm(cbind(...) ~ factor) + Manova()

# Sum-to-zero contrasten zetten — voorwaarde voor correcte Type-III SS.
# Defensief patroon: ook al staat deze regel in je setup-chunk, herhaal
# hem vlak boven elke Manova(., type = 3)- of Anova(., type = 3)-aanroep.
# Globale state kan tussentijds gewijzigd zijn door package-load of een
# eerdere chunk; deze regel kost niets en sluit één bron van rare output uit.
options(contrasts = c("contr.sum", "contr.poly"))

# Stap 1: multivariaat lineair model.
lm_leeuw <- lm(cbind(trots, jaloezie, eenzaamheid, woede) ~ sociale_status,
               data = leeuw_gevoelens)

# Stap 2: MANOVA met Type-III SS.
res_leeuw <- Manova(lm_leeuw, type = 3)

# Stap 3: de vier multivariate toetsen. univariate = FALSE: de univariate
# F's uit Manova() kunnen fout zijn (ongelijke df) -- haal ze uit summary.aov().
summary(res_leeuw, multivariate = TRUE, univariate = FALSE)

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

# Alle vier komen tegelijk uit summary(), met test stat, F-approximatie, df, p.
# Rapporteer ze alle vier: Pillai's V eerst (het robuustst), de andere drie
# ernaast.

Univariate follow-up + Bonferroni

# Univariate F per DV.
summary.aov(lm_leeuw)

# Verzamel de p-waardes en Bonferroni-corrigeer.
ps <- sapply(summary.aov(lm_leeuw), function(x) x[["Pr(>F)"]][1])
p.adjust(ps, method = "bonferroni")

# Of Holm-Bonferroni (krachtiger, controleert ook FWER).
p.adjust(ps, method = "holm")

# Tukey HSD per significante DV.
TukeyHSD(aov(trots       ~ sociale_status, data = leeuw_gevoelens))
TukeyHSD(aov(jaloezie    ~ sociale_status, data = leeuw_gevoelens))
TukeyHSD(aov(eenzaamheid ~ sociale_status, data = leeuw_gevoelens))
TukeyHSD(aov(woede       ~ sociale_status, data = leeuw_gevoelens))

Descriptive Discriminant Analysis (DDA)

# Hergebruik het lm-object.
DA <- candisc(lm_leeuw)
DA                # eigenvalues + canonical R^2 + percentages
summary(DA)       # plus structure matrix + group means

# Structure matrix expliciet — DV's tegen functies.
DA$structure

# Standardized coefficients.
DA$coeffs.std

# Group means op de discriminant-assen.
DA$means

Visualisatie — discriminant-plot

# Group positions op de eerste twee discriminant-functies.
plot(DA, ellipse = TRUE, col = c("#0077BB", "#33BBEE", "#EE7733"))

Voorbeeld-tentamenvragen

Even oefenen op tentamen-toon

Onderaan dit hoofdstuk staan vier korte theorievragen in tentamen-stijl en één mini R-opdracht. Hou de vuistregels paraat: vier teststatistics, Box’s M op $\alpha = .001$, Bonferroni $.05/p$, Pillai eerst en de andere drie ernaast, $\min(I-1, p)$ discriminant-functies.

Theorie en handreken

OpmerkingVraag E1 — Vier teststatistics interpreteren

Een orka-bioloog onderzoekt drie pod-types op vier emotionele DV’s. Hij rapporteert:

                  Df test stat approx F num Df den Df    Pr(>F) 

Pillai             2     0.419     5.96      8    180    < .001 
Wilks              2     0.609     6.27      8    178    < .001 
Hotelling-Lawley   2     0.598     6.58      8    176    < .001 
Roy                2     0.510    11.47      4     90    < .001 

Welke uitspraak is correct?

1. De vier teststatistics geven sterk verschillende $p$-waardes; we moeten Roy eerst rapporteren.
2. De vier teststatistics zijn allemaal significant; we rapporteren eerst Pillai’s $V = 0.42$, $F(8, 180) = 5.96$, $p < .001$, en de andere drie ernaast.
3. Wilks’ $\Lambda = .609$ betekent dat het effect klein is; minder dan $40\%$ verklaard.
4. Roy’s grote $F = 11.47$ verraadt een schending van multivariate normaliteit.

OpgeletAntwoord E1 — open na je eigen poging

b) Allemaal significant; Pillai rapporteer je eerst, de andere drie ernaast. Optie a is fout (alle vier wijzen dezelfde kant op, geen reden Roy te prefereren). Optie c verwart Wilks (lager = sterker effect) — $\Lambda = .609$ duidt juist op een matig effect ($1 - \Lambda = .39$ als ruwe $\eta^2$-analoog). Optie d is onzin: Roy’s hoge $F$ komt door zijn focus op de grootste eigenvalue, geen aanname-schending.

OpmerkingVraag E2 — Box’s M significant op $\alpha = .001$

Een orca-onderzoeker vindt:

Box's M-test for Homogeneity of Covariance Matrices
M = 122.92, df = 75, p = 4.06e-04

De celgroottes zijn $n_1 = 70$, $n_2 = 70$, $n_3 = 70$. Wat is de juiste vervolgstap?

1. Stop met MANOVA; de aanname is geschonden, gebruik aparte ANOVA’s.
2. Box’s M is significant op $\alpha = .001$, maar door balanced design + $n_i = 70$ is MANOVA-$F$ robuust; ga door en vermeld de bevinding in de discussie.
3. Gebruik bootstrap-MANOVA als enige geldige analyse.
4. Pas Holm-Bonferroni toe op de Box’s M $p$-waarde.

OpgeletAntwoord E2 — open na je eigen poging

b) $p = 4 \times 10^{-4} < .001$ dus Box’s M is significant op de strenge $\alpha = .001$. Maar gebalanceerde celgroottes met $n_i \geq 20$ maken MANOVA-$F$ robuust — doorgaan en het in de discussie melden. Optie a is overdreven (en de FWER-discussie van T2 vervalt dan); c is technisch mogelijk maar niet de standaard route; d is onzinnig (Holm-Bonferroni werkt niet op enkele $p$-waardes).

OpmerkingVraag E3 — Bonferroni-correctie hand-rekenen

Na een significante MANOVA met $4$ DV’s krijg je vier univariate $p$-waardes:

DV1: F(2, 87) = 8.42,  p < .001 
DV2: F(2, 87) = 4.51,  p = .014 
DV3: F(2, 87) = 2.10,  p = .13 
DV4: F(2, 87) = 5.83,  p = .004 

Welke DV’s blijven significant na Bonferroni-correctie?

1. DV1 alleen.
2. DV1 en DV4.
3. DV1, DV2 en DV4.
4. Alle vier.

OpgeletAntwoord E3 — open na je eigen poging

b) De Bonferroni-drempel is $.05/4 = .0125$. DV1 ($p < .001$) ✓; DV2 ($p = .014$, groter dan $.0125$) ✗; DV3 ($p = .13$) ✗; DV4 ($p = .004$) ✓. Alleen DV1 en DV4 overleven Bonferroni. (Onder Holm-Bonferroni zou DV2 het wel halen: DV2 heeft de derde-kleinste $p$, dus zijn drempel is $.05/2 = .025$, en $.014$ ligt daaronder. Maar de vraag specificeert Bonferroni.)

OpmerkingVraag E4 — DDA-eigenvalues interpreteren

Een onderzoeker rapporteert na een MANOVA met $I = 4$ groepen en $p = 5$ DV’s:

Canonical Discriminant Analysis:

   CanRsq   Eigenvalue   Difference   Percent   Cumulative 

1  0.682    2.143        1.901        86.7      86.7 
2  0.195    0.242        0.155        9.8       96.5 
3  0.080    0.087                     3.5       100.0 

4  ---      ---          ---          ---       --- 

Welke uitspraak past het best?

1. Er zijn vier discriminant-functies; alle vier moeten geïnterpreteerd worden.
2. Het maximum aantal discriminant-functies is $\min(I-1, p) = \min(3, 5) = 3$; de eerste verklaart $87\%$ van de groepsvariantie en is dominant.
3. Een eigenvalue van $2.14$ betekent dat $214\%$ van de variantie verklaard wordt — dat kan niet, dus de tabel is fout.
4. De ontbrekende rij 4 betekent dat $p$-waardes voor functie 4 missen.

OpgeletAntwoord E4 — open na je eigen poging

b) $\min(I-1, p) = 3$ — er zijn er nooit meer dan $\min(\text{aantal groepen} - 1, \text{aantal DV's})$. De eerste functie domineert ($87\%$); functies 2 en 3 spelen kleine rollen. Optie a is fout (er zijn er drie, niet vier); c verwart eigenvalue met variantie-percentage (eigenvalue is een schaal-vrije matrix-grootheid die boven $1$ kan); d ook fout (er zijn maar drie functies).

R-practical opdrachtje

De zeehonden van drie eilanden

OpmerkingVraag E5 — Mini-MANOVA bij de zeehond

Een ecoloog onderzoekt $84$ zeehonden op drie typen eilanden (klif / strand / binnenmeer). Per zeehond meet hij drie psychologische dimensies: nieuwsgierigheid, voorzichtigheid en sociale_speelsheid (alle continu, schaal $0$-$100$). De dataset staat in data/zeehond_eilanden.RData en bevat het object zeehond_eilanden. Sla je R-commando’s op in één scriptbestand: zeehond.R. Gebruik $\alpha = .05$ voor de multivariate toets en Bonferroni-correctie voor univariate follow-up.

a) Toets de aanname van gelijke covariantie-matrices met box_m(). Rapporteer de toetsstatistiek, de vrijheidsgraden en de $p$-waarde. Ga je door met MANOVA?

b) Fit het multivariaat lineair model met lm(cbind(nieuwsgierigheid, voorzichtigheid, sociale_speelsheid) ~ eilandtype, data = zeehond_eilanden) en run Manova(., type = 3) met summary(., multivariate = TRUE). Rapporteer eerst Pillai’s $V$ met zijn $F$-approximatie, $df$ en $p$, en zet de andere drie toetsen ernaast.

c) Voer univariate follow-up uit per DV (summary.aov()) en pas Bonferroni-correctie toe op de drie $p$-waardes. Welke DV’s blijven significant na correctie?

d) Run een DDA met candisc::candisc(). Hoeveel discriminant-functies zijn er? Welke DV correleert het sterkst met de eerste functie volgens de structure matrix? Schrijf één APA-zin.

OpgeletAntwoord E5 — open na je eigen poging

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

'data.frame':   84 obs. of  4 variables:
 $ eilandtype        : Factor w/ 3 levels "klif","strand",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ nieuwsgierigheid  : num  26.6 44.6 53.4 26.6 21 36.5 49.2 30.2 11.6 29.1 ...
 $ voorzichtigheid   : num  71.5 53.4 72.8 60.5 82.7 72.6 86.8 84.6 93 57 ...
 $ sociale_speelsheid: num  30.7 67.3 43.8 11.8 33.7 57.4 37.1 34.5 22 21.8 ...

box_m(data  = zeehond_eilanden[, c("nieuwsgierigheid", "voorzichtigheid", "sociale_speelsheid")],
      group = zeehond_eilanden$eilandtype)

# A tibble: 1 × 4
  statistic p.value parameter method                                            
      <dbl>   <dbl>     <dbl> <chr>                                             
1      3.55   0.990        12 Box's M-test for Homogeneity of Covariance Matric…

# Defensief: contrasten ook hier expliciet zetten vlak vóór de Type-III-aanroep.
options(contrasts = c("contr.sum", "contr.poly"))
lm_z <- lm(cbind(nieuwsgierigheid, voorzichtigheid, sociale_speelsheid) ~ eilandtype,
           data = zeehond_eilanden)
res_z <- Manova(lm_z, type = 3)
summary(res_z, multivariate = TRUE, univariate = FALSE)

Type III MANOVA Tests:

Sum of squares and products for error:
                   nieuwsgierigheid voorzichtigheid sociale_speelsheid
nieuwsgierigheid          11846.951      -5178.9300          4003.8521
voorzichtigheid           -5178.930      16235.3871          -917.5786
sociale_speelsheid         4003.852       -917.5786         12354.8329

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

Sum of squares and products for the hypothesis:
                   nieuwsgierigheid voorzichtigheid sociale_speelsheid
nieuwsgierigheid           239242.8        255834.8           240657.0
voorzichtigheid            255834.8        273577.6           257347.2
sociale_speelsheid         240657.0        257347.2           242079.6

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1   0.98515 1747.472      3     79 < 2.22e-16 ***
Wilks             1   0.01485 1747.472      3     79 < 2.22e-16 ***
Hotelling-Lawley  1  66.35969 1747.472      3     79 < 2.22e-16 ***
Roy               1  66.35969 1747.472      3     79 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: eilandtype 

Sum of squares and products for the hypothesis:
                   nieuwsgierigheid voorzichtigheid sociale_speelsheid
nieuwsgierigheid           21012.75       -12520.20           14419.44
voorzichtigheid           -12520.20        12508.09          -15873.07
sociale_speelsheid         14419.44       -15873.07           20397.78

Multivariate Tests: eilandtype
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            2  1.215023 41.27590      6    160 < 2.22e-16 ***
Wilks             2  0.143815 43.10579      6    158 < 2.22e-16 ***
Hotelling-Lawley  2  3.458252 44.95727      6    156 < 2.22e-16 ***
Roy               2  2.432500 64.86665      3     80 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Univariate F's per DV -- uit summary.aov(), niet uit Manova().
summary.aov(lm_z)

 Response nieuwsgierigheid :
            Df Sum Sq Mean Sq F value    Pr(>F)    
eilandtype   2  21013 10506.4  71.834 < 2.2e-16 ***
Residuals   81  11847   146.3                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response voorzichtigheid :
            Df Sum Sq Mean Sq F value           Pr(>F)    
eilandtype   2  12508  6254.0  31.202 0.00000000008972 ***
Residuals   81  16235   200.4                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response sociale_speelsheid :
            Df Sum Sq Mean Sq F value    Pr(>F)    
eilandtype   2  20398 10198.9  66.865 < 2.2e-16 ***
Residuals   81  12355   152.5                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ps <- sapply(summary.aov(lm_z), function(x) x[["Pr(>F)"]][1])
round(ps, 5)

   Response nieuwsgierigheid     Response voorzichtigheid 
                           0                            0 
 Response sociale_speelsheid 
                           0 

round(p.adjust(ps, method = "bonferroni"), 5)

   Response nieuwsgierigheid     Response voorzichtigheid 
                           0                            0 
 Response sociale_speelsheid 
                           0 

DA_z <- candisc(lm_z)
DA_z

Canonical Discriminant Analysis for eilandtype:

   CanRsq Eigenvalue Difference Percent Cumulative
1 0.70867     2.4325     1.4067  70.339     70.339
2 0.50636     1.0258     1.4067  29.661    100.000

Test of H0: The canonical correlations in the 
current row and all that follow are zero

  LR test stat approx F numDF denDF   Pr(> F)    
1      0.14381   43.106     6   158 < 2.2e-16 ***
2      0.49364      NaN     2   NaN       NaN    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(DA_z$structure, 3)

                     Can1   Can2
nieuwsgierigheid    0.813  0.581
voorzichtigheid    -0.775  0.134
sociale_speelsheid  0.907 -0.282

a) Box’s $M = 3.55$, $df = 12$, $p = .99$ — verre van significant op zelfs de gewone $\alpha = .05$. Aanname comfortabel gehaald, ga door.

b) Pillai’s $V = 1.22$, $F(6, 160) = 41.28$, $p < .001$ — sterk multivariaat effect. Wilks, Hotelling-Lawley en Roy wijzen dezelfde kant op (alle vier $p < .001$); die zet je ernaast, zoals in de APA-zin hieronder. De vrijheidsgraden zijn num Df en den Df uit de Pillai-rij, niet de Df-kolom vooraan. Dat $V$ boven de $1$ uitkomt mag: bij drie groepen ligt het maximum op $2$ (zie de tabel met de vier teststatistics bovenaan). Je commando uit de vraag, summary(res_z, multivariate = TRUE), laat de univariate $F$’s standaard al weg; hier staat univariate = FALSE er voor de duidelijkheid bij.

c) Univariate $F$’s uit summary.aov(): nieuwsgierigheid $F(2, 81) = 71.83$, $p < .001$; voorzichtigheid $F(2, 81) = 31.20$, $p < .001$; sociale_speelsheid $F(2, 81) = 66.87$, $p < .001$. Vergelijk elke $p$ met de Bonferroni-drempel $.05/3 \approx .017$. Alle drie liggen ver onder $.001$, dus alle drie blijven significant. round(ps, 5) toont $0$: dat is afronding op vijf decimalen, geen echte nul — je rapporteert $p < .001$. (De regel met p.adjust() komt op hetzelfde uit, langs de andere kant: die vermenigvuldigt elke $p$ met $3$, en dan vergelijk je met $.05$.) Haal deze $F$’s niet uit summary(res_z, univariate = TRUE): daar staat voor voorzichtigheid $F = 62.40$, twee keer te groot, en welke $F$’s het treft hangt af van het aantal DV’s — vertrouw die tabel dus niet (zie de waarschuwing in 5.1.c).

d) $\min(I-1, p) = \min(2, 3) = 2$ discriminant-functies. Eigenvalues $\lambda_1 = 2.43$, $\lambda_2 = 1.03$. Structure matrix: met Can1 correleert sociale_speelsheid het sterkst ($r = .91$), nieuwsgierigheid ook hoog ($r = .81$), en voorzichtigheid met het andere teken ($r = -.78$) — Can1 is de “bedachtzaamheid versus uitbundigheid”-as. Staat er bij de toets van functie 2 NaN? Dat zegt niets over die functie; zie het NaN-kader in 5.2.c.

APA-stijl. Eerst de MANOVA met de follow-up, dan de DDA uit vraag d.

Een eenwegs-MANOVA toonde een sterk multivariaat effect van eilandtype op het gedragsprofiel van zeehonden, Pillai’s $V = 1.22$, $F(6, 160) = 41.28$, $p < .001$. De andere drie teststatistics wezen in dezelfde richting (Wilks’ $\Lambda = .144$, $F(6, 158) = 43.11$, $p < .001$; Hotelling-Lawley $T^2 = 3.46$, $F(6, 156) = 44.96$, $p < .001$; Roy’s $\theta = 2.43$, $F(3, 80) = 64.87$, $p < .001$). Univariate follow-up met Bonferroni-correctie ($\alpha = .017$) liet zien dat alle drie DV’s significant verschilden tussen eilandtypes (alle $p < .001$). Klif-zeehonden waren voorzichtiger en minder speels; strand-zeehonden waren juist speelser en minder voorzichtig; binnenmeer-zeehonden combineerden hoge nieuwsgierigheid met gemiddelde voorzichtigheid.

Een discriminantanalyse leverde twee functies op, waarvan de eerste (eigenvalue $= 2.43$, $70.3\%$ van de verklaarde groepsvariantie) vooral werd gedragen door sociale speelsheid ($r = .91$) en nieuwsgierigheid ($r = .81$), tegenover voorzichtigheid ($r = -.78$).

Klaar? Zet de codering terug

Je hebt dit hele thema gerekend met options(contrasts = c("contr.sum", "contr.poly")). Die instelling blijft staan zolang je R-sessie openstaat — ook als je straks aan een ander thema begint. Zet hem daarom na afloop terug op R’s standaard, net als aan het eind van thema 2 en 3:

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

Vergeet je dit, dan rekent elke volgende analyse met een categorische voorspeller in dezelfde sessie met de verkeerde codering: een glm() met een factor, zoals de diersoort in thema 4, geeft dan andere coëfficiënten dan in de uitwerkingen.

Wat blijft liggen

Onderwerpen die buiten dit thema blijven

MANOVA met Manova(., type = 3) en Bonferroni-univariate-follow-up is de hoofdingang van multivariate variantieanalyse, maar er ligt veel achter. Wat dit oefenboek je gaf is genoeg voor de gangbare onderzoeks-rapportage; voor verdere studie:

• Predictive Discriminant Analysis PDA — gebruikt discriminant-functies om nieuwe gevallen te classificeren in groepen. Pakketten MASS::lda() en MASS::qda(). Hier gebruikten we descriptive DA als follow-up, niet als classifier. Eén stukje ervan ligt dichterbij: een classificatietabel lezen — het kader Wat je na dit thema moet kunnen bovenaan noemt het al.
• MANCOVA — MANOVA met continue covariaten erbij. Combineert thema 3 (ANCOVA) en thema 5: lm(cbind(Y1, Y2, Y3) ~ factor + covariate). Aanname-pakket wordt zwaarder (gelijke regressie-hellingen per groep, multivariaat).
• Robuuste MANOVA — bij niet-normale of outlier-rijke data: Wilcox’ robuuste MANOVA met getrimde gemiddelden of M-schattingen. Niet behandeld; pakket WRS2.
• Bayesiaanse MANOVA — multivariate posterior-verdelingen in plaats van $p$-waardes. Toegankelijke ingang via brms::brm(mvbind(Y1, Y2, Y3) ~ factor).
• Mixed multivariate models — voor geclusterde of herhaalde-metingen-data met meerdere DV’s tegelijk. Combineert MANOVA met multilevel-structuur. Pakket MCMCglmm, of brms met multivariate response.
• Structural Equation Modeling SEM — multivariate analyse als generalisatie: pad-modellen, latent-variable-modellen, MIMIC-modellen voor groepsverschillen op een latente $Y$. Pakket lavaan. Veel rijker dan MANOVA, ook complexer.
• Repeated-measures als multivariate approach — RMA via MANOVA-formulering: hetzelfde model, andere interpretatie. Komt impliciet terug in thema 6.

Voor verdieping: Tabachnick & Fidell (2019, Using Multivariate Statistics); Maxwell, Delaney & Kelley (2018, Designing Experiments and Analyzing Data); Stevens (2009, Applied Multivariate Statistics for the Social Sciences).

Aan het eind van de dag

Toen de zon laag stond op de savanne, klapte de leeuw zijn schrift dicht.

“De drie groepen wonen elk in een eigen gevoels-wereld,” zei hij. “Alpha’s met hun trots, middenrangers met hun jaloezie, outsiders met hun eenzaamheid. Niet één gevoel, maar een heel patroon.”

“En het is niet zo dat de ene groep méér voelt en de andere minder,” vervolgde hij. “Het is een ander mengsel. Pillai weet het, Wilks weet het, en zelfs Roy met zijn smalle blik op één eigenvalue zegt hetzelfde.”

“Drie functies?” vroeg de bever, die op de rotsen kwam zitten.

“Twee,” verbeterde de leeuw. “Eén voor trots versus eenzaamheid, één voor jaloezie en woede. Die tweede staat los — middenrangers herkennen zichzelf erin, alpha’s en outsiders niet.”

“Dus de troep heeft een geometrie,” zei de bever.

“Ja,” zei de leeuw. “Drie statussen, twee dimensies, één multivariate $F$ die alles samenvat. Morgen zal ik weer kijken — dezelfde dieren, andere stemming, kijken of het profiel mee verschuift.”

Een leeuwin liep langs, schichtig, en keek even naar het rotsplateau. De leeuw keek niet meer op. Hij wist al onder welke discriminant-as ze paste.

Verantwoording

Dit oefenboek is geschreven voor studenten die multivariate variantieanalyse leren via R. De didactische lijnen volgen de gangbare opbouw van Nederlandse universitaire MVDA-cursussen; alle voorbeelden, datasets, vragen en formuleringen in dit hoofdstuk zijn origineel. De gebruikte drempels en formules zijn standaard statistische conventies; verwijzingen naar Tabachnick & Fidell (2019) en Stevens (2009) volgen de gebruikelijke citaatpraktijk in dit veld.

Versie: October 2026 — CountCamp Lab