Vierde beweging · Verschil onthechten

Oefening 10.1 · Als de moderator een groep is

Interactie met een groep als moderator · categorical moderation

Data bij deze oefening — penguins.csv 342 pinguïns · diamantjes.csv de twaalf steentjes uit het boek · wat elke kolom betekent

prikkel

Eén getal dat uit twee getallen bleek te bestaan

helling per groep

Bij de enkelvoudige regressie trok je een lijn door 342 pinguïns en las je er één getal van af: 49.69 gram per millimeter.¹ Elke millimeter vin erbij, bijna vijftig gram pinguïn erbij. Eén getal, netjes uitgerekend, en er was geen enkele reden om het te wantrouwen.

Splits diezelfde 342 beesten nu eens in tweeën — de Gentoo’s aan de ene kant, de Adélies en Chinstraps aan de andere — en trek de lijn twee keer. In deze steekproef komt er bij de niet-Gentoo’s 29.49 uit, en bij de Gentoo’s 54.62.

Kijk daar even naar. 49.69 ligt tussen die twee in, en het is geen van beide.

En het is ook niet netjes hun gemiddelde — dat zou te makkelijk zijn. Weeg de twee hellingen naar groepsgrootte (219 niet-Gentoo’s tegen 123 Gentoo’s) en je komt uit op 38.53. Je eerste getal ligt daar ruim elf gram bóven.

Dat komt doordat één lijn door alle pinguïns twee dingen tegelijk moet doen: het verband bínnen elke groep volgen, én de afstand tússen de groepen overbruggen. En die tweede taak duwt hem omhoog, want de Gentoo’s zijn niet alleen zwaarder maar hebben óók langere vinnen. Trek je een lijn door niets anders dan de twee groepsgemiddelden, dan is die zelfs 53.73.

Je eerste getal was dus niet fout gerekend. Het was een menging — en niet eens een eerlijke.

En dat roept meteen de vraag op waar deze oefening over gaat. Is die kloof echt, of kan het toeval van 342 beesten hem namaken? En als hij echt is — hoe schrijf je dan een model op dat niet één helling wil, maar twee?

speel het

Twee lijnen: één groep tegenover de rest

dummy · interactieterm · simple slope

Eerst zoals in het boek — de twaalf diamantjes

d <- read.csv("diamantjes.csv")

# het sterretje na de variabelen vraagt om beide effecten PLUS hun product
model.glad.x.D <- lm(glans ~ gladheid * D, data = d)
summary(model.glad.x.D)
# gladheid     b = 6.14
# gladheid:D   b = 2.64   p = .054

lm(), t.test() en aov() zijn in beide talen precies dezelfde functie — daar valt niets te vertalen. Het verschil tussen base R en dplyr zit in het klaarzetten van de data, niet in het model.

D is een kolom met nullen en enen: nul voor golfje, één voor sterretje. Zo’n nul-of-één-kolom heet een dummy, en hij is het eenvoudigste soort moderator die er bestaat — twee groepen, meer niet. De interactie gladheid × merk komt uit op b = 2.64 met p = .054: bij het ene merk telt een punt gladheid harder mee voor de glans dan bij het andere. Eén helling is dan niet meer genoeg.

Nu jij — dezelfde beesten als in oefening 3.2, nieuwe vraag. Is die 49.69 één helling, of zijn het er twee?

De pinguïns hebben geen kant-en-klare dummy in het bestand, dus die maak je zelf: één kolom die zegt of een pinguïn een Gentoo is of niet.

peng <- read.csv("penguins.csv")

# TRUE/FALSE wordt 1/0 zodra je er een getal van maakt
peng$gentoo <- as.numeric(peng$soort == "Gentoo")
table(peng$gentoo)      # → 219 nullen, 123 enen

library(dplyr)
peng <- read.csv("penguins.csv") |>
  mutate(gentoo = as.numeric(soort == "Gentoo"))

count(peng, gentoo)     # → 219 nullen, 123 enen

Nu het model. De ster * tussen vinlengte en gentoo is de toestemming: reken het verband tussen vinlengte en gewicht uit, maar laat de Gentoo’s er een eigen helling op na houden.

In R · trek eraan

model.vl.x.gentoo <- lm(gewicht ~ vinlengte * gentoo, data = peng)
summary(model.vl.x.gentoo)

# vinlengte          b =  29.49   <- de helling bij gentoo = 0
# gentoo             b = -4841.81
# vinlengte:gentoo   b =  25.13   p < .001   <- hoevéél steiler de Gentoo's lopen
# Multiple R-squared: 0.785

Vier getallen, en ze doen alle vier iets anders. Het loont om er even bij stil te staan, want dit is de tabel die je de rest van de oefening blijft lezen:

vinlengte is de helling bij gentoo = 0, dus bij de niet-Gentoo’s: 29.49 gram per millimeter.
gentoo is het hoogteverschil bij vinlengte = 0 — een vin van nul millimeter, dus dat getal woont op een adres dat niet bestaat. Laat het staan; in oefening 10.2 gaan we er iets aan doen.
vinlengte:gentoo is waar het om gaat: 25.13 gram per millimeter extra voor de Gentoo’s. Niet hun helling — het verschil tussen hun helling en die van de rest.
En die p < .001 zegt: zo’n kloof tussen twee hellingen zie je bijna nooit als er in werkelijkheid geen kloof is.

De helling van de Gentoo’s zelf staat er dus niet; die tel je op. En dat is de hele truc:

29.49 + 25.13      # → 54.62 gram per millimeter bij de Gentoo's

Twee hellingen, allebei uit één tabel: 29.49 en 54.62. Zulke hellingen binnen één groep heten simple slopes, en je hebt er niets voor nodig behalve optellen. Geloof je het niet, reken ze dan los na — het moet op de komma hetzelfde geven:

In R · ter controle, de twee lijnen los

model.niet.gentoo <- lm(gewicht ~ vinlengte, data = subset(peng, gentoo == 0))
model.wel.gentoo  <- lm(gewicht ~ vinlengte, data = subset(peng, gentoo == 1))

coef(model.niet.gentoo)[2]   # → 29.49
coef(model.wel.gentoo)[2]    # → 54.62

En daar ligt je 49.69 uit oefening 3.2 uit elkaar. Eén lijn door twee groepen die elk hun eigen lijn hebben, geeft een getal dat bij geen van beide groepen hoort — te steil voor de Adélies en Chinstraps, te flauw voor de Gentoo’s. Het is een gemiddelde waar niemand woont.

Teken het maar, dan zie je het in één keer:

In R · de twee lijnen tekenen

plot(peng$vinlengte, peng$gewicht, col = peng$gentoo + 1,
     xlab = "vinlengte (mm)", ylab = "gewicht (g)")

# per groep een lijn, in dezelfde kleur als de punten -- de twee modellen
# van hierboven, dus je hoeft ze niet opnieuw te bouwen
abline(model.niet.gentoo, col = 1, lwd = 2)
abline(model.wel.gentoo,  col = 2, lwd = 2)

# en de ene lijn uit oefening 3.2 er als stippellijn doorheen
model.vl <- lm(gewicht ~ vinlengte, data = peng)
abline(model.vl, lty = 2, lwd = 2)

Twee lijnen die uit elkaar lopen, en een stippellijn die er tussendoor probeert te kruipen. Die stippellijn is 49.69.

Nog twee lijnen, maar dan parallel: de sekse

niet-significante interactie · hoogte tegenover steilheid

Voordat je dit doortrekt naar meer groepen, eerst een tegenvoorbeeld — want anders ga je overal waaiers zien. In penguins.csv staat nóg een kolom met twee groepen: sekse. Mannetjes zijn duidelijk zwaarder dan vrouwtjes; in deze steekproef wegen ze gemiddeld 4545.68 gram tegen 3862.27 — een verschil van 683.41 gram. Dan zal een millimeter vin bij mannetjes ook wel zwaarder tellen?

In R · dezelfde zet, andere groep

# de negen pinguïns zonder genoteerde sekse vallen er vanzelf uit
model.vl.x.sekse <- lm(gewicht ~ vinlengte * sekse, data = peng)
summary(model.vl.x.sekse)
# vinlengte:seksemale   b = -0.29   p = .92   <- 47.15 - 0.29 = 46.86

# helling per sekse los
model.vrouw <- lm(gewicht ~ vinlengte, data = subset(peng, sekse == "female"))
model.man   <- lm(gewicht ~ vinlengte, data = subset(peng, sekse == "male"))

coef(model.vrouw)[2]   # → 47.15
coef(model.man)[2]     # → 46.86

47.15 tegen 46.86. Die liggen zo goed als op elkaar, en de interactie is dan ook nergens: p = .92. De mannetjes zijn zwaarder, maar een millimeter vin weegt bij hen niet zwaarder.

Dat is het onderscheid dat de rest van deze oefening draagt. “Mannetjes zijn zwaarder” is een verschil in hoogte — twee wolken, de ene boven de andere, met lijnen die keurig parallel lopen. Een interactie is een verschil in steilheid — lijnen die uit elkaar waaieren. Twee groepen kunnen ver uit elkaar liggen en tóch precies dezelfde helling hebben, en dan is er geen interactie te bekennen.

Drie lijnen: de hele soort erbij

nominale moderator · modelvergelijking

Een dummy kon maar twee groepen aan. Maar soort heeft er drie — Adélie, Chinstrap, Gentoo — en dat verandert meer dan je zou denken.

In R · drie soorten, drie hellingen

model.vl.x.soort <- lm(gewicht ~ vinlengte * soort, data = peng)
summary(model.vl.x.soort)

# vinlengte:soortChinstrap   b =  1.74   p = .82
# vinlengte:soortGentoo      b = 21.79   p = .002

Onderin staan nu twéé interactie-regels in plaats van één. R neemt de Adélie als ijkpunt — alfabetisch de eerste — en legt de andere twee soorten daarnaast; een regel vinlengte:soortAdelie bestaat daarom niet. Wat je leest is dus: de Chinstraps lopen 1.74 gram per millimeter steiler dan de Adélies (en dat is niets, p = .82), en de Gentoo’s 21.79 steiler (en dat is wél iets, p = .002).

En hier houdt “lees de coëfficiënt af” op te werken. Wil je weten of de helling ergens tussen de soorten verschilt, dan heb je één getal nodig over die twee regels samen — en dat staat nergens in de tabel. Je krijgt het door twee modellen naast elkaar te leggen: eentje waarin alle soorten dezelfde helling moeten delen (+), en eentje waarin elke soort een eigen helling mag hebben (*).

In R · de twee modellen naast elkaar

# alle soorten moeten één helling delen -- dat is de + in plaats van de *
model.vl.soort <- lm(gewicht ~ vinlengte + soort, data = peng)

anova(model.vl.soort, model.vl.x.soort)   # model.vl.x.soort is de * van hierboven
# F(2, 336) = 5.53   p = .004

F(2, 336) = 5.53, p = .004. Twee vrijheidsgraden, want er kwamen twee regels bij. Dezelfde trechter als bij de eenweg ANOVA: eerst overkoepelend of er íéts is, en pas daarna per regel waar het zit. Hier zit het bij de Gentoo’s.

De drie hellingen los:

In R · de helling per soort

model.adelie    <- lm(gewicht ~ vinlengte, data = subset(peng, soort == "Adelie"))
model.chinstrap <- lm(gewicht ~ vinlengte, data = subset(peng, soort == "Chinstrap"))
model.gentoo    <- lm(gewicht ~ vinlengte, data = subset(peng, soort == "Gentoo"))

coef(model.adelie)[2]      # → 32.83
coef(model.chinstrap)[2]   # → 34.57
coef(model.gentoo)[2]      # → 54.62

32.83, 34.57 en 54.62 gram per millimeter vin. Twee soorten die vrijwel samenvallen, en één die er ver bovenuit steekt. Trek de flauwste van de steilste af — 54.62 − 32.83 — en je krijgt 21.79. Dat is exact de coëfficiënt vinlengte:soortGentoo uit de tabel hierboven, en dat is geen toeval: de Adélie is het ijkpunt én toevallig ook de flauwste, dus “hoeveel steiler dan de Adélie” en “hoeveel tussen de steilste en de flauwste” zijn hier hetzelfde getal. Onthoud het, want in de volgende trede komt het terug. Merk op dat de Gentoo’s op precies dezelfde 54.62 uitkomen als in de eerste trede — logisch, het is dezelfde groep pinguïns met dezelfde lijn. Wat veranderde is waarmee je ze vergelijkt: eerst met alle andere pinguïns op één hoop, nu met elke soort apart.

Zes hokjes: twee groeperingen tegelijk

tweeweg-ANOVA · celgemiddelde · type III

Tot nu toe modereerde een groep het effect van een getal op het gewicht. Maar twee groeperingen kunnen elkaar ook onderling moduleren, en dan blijft er van de lijnen niets over — je houdt hokjes over.

Zet sekse (twee groepen) naast soort (drie groepen) en je krijgt zes hokjes. Dat heet een tweeweg-ANOVA, of een factoriële ANOVA; dat is hetzelfde. Twee groeperingen, dus tweeweg — net zoals de vorige oefening er één had en eenweg heette. De vraag is dan niet meer hoe steil een lijn loopt, maar of het gemiddelde per hokje zich netjes gedraagt.

In R · eerst de zes hokjes zelf

# de negen pinguïns zonder genoteerde sekse laten we hier expliciet weg
ps <- subset(peng, sekse %in% c("male", "female"))

table(ps$sekse, ps$soort)                                    # hoeveel zitten er in elk hokje
tapply(ps$gewicht, list(ps$sekse, ps$soort), mean)           # en wat weegt zo'n hokje gemiddeld

	Adélie	Chinstrap	Gentoo
vrouwtje	3368.8 g (n = 73)	3527.2 g (n = 34)	4679.7 g (n = 58)
mannetje	4043.5 g (n = 73)	3939.0 g (n = 34)	5484.8 g (n = 61)

Lees die tabel eens per kolom. Bij de Adélies zijn de mannetjes 674.7 gram zwaarder dan de vrouwtjes, bij de Chinstraps 411.8, bij de Gentoo’s 805.1. Het sekseverschil is er bij alle drie — maar het is niet overal even groot. Dát is de interactie: het verschil in verschil. Niet of sekse uitmaakt, maar of het bij elke soort evenveel uitmaakt.

En let op wat je nog steeds kunt: je kunt dit tekenen. Zes punten, twee lijntjes die elkaar volgen of niet. Dat is wat de hele 10.1 bij elkaar houdt — of de moderator nu twee, drie of zes hokjes maakt, je kunt het plaatje voor je zien.

In R · het plaatje bij de tabel

# soort op de horizontale as, één lijn per sekse -- lijnen boven staafjes,
# want alleen bij lijnen zie je of ze parallel lopen
interaction.plot(ps$soort, ps$sekse, ps$gewicht,
                 xlab = "soort", ylab = "gemiddeld gewicht (g)", trace.label = "sekse")

Nu de toets. En hier komt een keuze om de hoek kijken die je bij één factor nooit tegenkwam. Hij kost je twee regels code, maar je wilt weten waaróm ze er staan.

Wie krijgt het overlapstuk? Type I en type III

kwadratensommen · type III · contrasten

Teken twee cirkels die elkaar deels overlappen. De ene is alles wat sekse van het gewicht kan verklaren, de andere alles wat soort kan verklaren. In het overlapstuk zit variantie waar ze allebei aanspraak op maken — en iemand moet beslissen wie hem krijgt.

Type I geeft hem aan wie het eerst genoemd wordt. Schrijf je sekse * soort, dan pakt sekse het overlapstuk en houdt soort alleen wat overblijft; schrijf je soort * sekse, dan is het andersom. Je uitkomst hangt dus af van je typvolgorde.

Type III geeft hem aan niemand. Elke factor krijgt alleen wat híj als enige kan verklaren. Dat is onafhankelijk van de volgorde, en dat is wat we hier aanhouden — het is ook wat je bij de meeste opleidingen en in de meeste tijdschriften tegenkomt.

Dan de tweede regel, die met contr.sum. Die gaat niet over je data maar over hoe R je groepen in getallen omzet. Standaard neemt R één groep als ijkpunt en legt de andere daarnaast — precies de nul-of-één-codering die je bij gentoo zelf met de hand maakte. Type III heeft een andere codering nodig: eentje waarin de groepen symmetrisch om nul liggen in plaats van allemaal ten opzichte van één ijkgroep. contr.sum is die codering. (contr.poly hoort er in dezelfde regel bij voor geordende factoren; die heb je hier niet, maar R wil dat je allebei noemt.)

Je hoeft niet te snappen hóé die codering werkt. Wat je moet onthouden is de volgorde van drie handelingen: aanzetten, rekenen, terugzetten.

In R · de tweeweg-ANOVA

# 1. aanzetten -- de codering die type III nodig heeft
options(contrasts = c("contr.sum", "contr.poly"))

# 2. rekenen
model.sekse.x.soort <- lm(gewicht ~ sekse * soort, data = ps)
car::Anova(model.sekse.x.soort, type = 3)
# sekse         F(1, 327) = 311.84   p < .001
# soort         F(2, 327) = 746.92   p < .001
# sekse:soort   F(2, 327) =   8.76   p < .001

lsr::etaSquared(model.sekse.x.soort, type = 3)
# sekse         eta.sq = .139
# soort         eta.sq = .664
# sekse:soort   eta.sq = .008

# 3. terugzetten voor de volgende analyse
options(contrasts = c("contr.treatment", "contr.poly"))

Vergeet je die eerste regel, dan rekent R gewoon door en geeft je andere getallen, zonder een kik. Kijk maar wat er met de F voor sekse gebeurt onder vier varianten van dezelfde data:

hoe je het uitrekent	F voor sekse	F voor de interactie
type I, `sekse` eerst opgeschreven	406.15	8.76
type I, `soort` eerst opgeschreven	387.46	8.76
type III met `contr.sum`	311.84	8.76
type III zonder `contr.sum`	173.55	8.76

Vier verschillende antwoorden op “heeft sekse effect op gewicht”, en één ervan is minder dan de helft van een andere. Maar kijk naar de rechterkolom: de interactie is in alle vier precies 8.76. Dat is geen toeval — de interactie is de term die als laatste komt, wat je ook doet, dus er valt voor hem niets te verdelen.

Reken het maar na, met drie verschillende functies naast elkaar — dezelfde interactie, drie wegen erheen:

In R · drie wegen, één interactie

# alle drie geven voor sekse:soort exact F(2, 327) = 8.757
anova(lm(gewicht ~ sekse * soort, data = ps))["sekse:soort", ]
summary(aov(gewicht ~ sekse * soort, data = ps))
car::Anova(lm(gewicht ~ sekse * soort, data = ps), type = 3)["sekse:soort", ]

En hoe weet je nu of je die contrast-regel vergeten bent? Niet aan de uitvoer. Daar is niets aan te zien — geen waarschuwing, geen rare waarde, gewoon andere getallen. Dat is precies het gemene eraan, en het is de reden dat die regel een gewoonte moet worden in plaats van iets wat je erbij bedenkt. Wat wél werkt is één regel die je vóór de analyse draait:

options("contrasts")     # de eerste moet "contr.sum" zijn, niet "contr.treatment"

En vergelijk je je uitkomst ooit met die van een medestudent, kijk dan éérst naar de interactie. Komt die overeen terwijl de hoofdeffecten verschillen, dan heeft niemand zich verrekend — dan heeft een van jullie tweeën een andere kwadratensom gebruikt. Wélke van de twee, zegt die vergelijking je niet; daarvoor moeten jullie allebei opschrijven wat je gedaan hebt. Vandaar dat het in de rapportage hoort.

schud het

Duizend werelden waarin de soort niets modereert

soort-labels husselen · kloof tussen hellingen

Terug naar de drie lijnen. Is die kloof tussen 32.83 en 54.62 echt, of kan het toeval hem namaken? We spelen weer vals. Deze keer husselen we de soort-labels los van de pinguïns: elke pinguïn houdt z’n eigen vinlengte en gewicht, maar krijgt een willekeurig soort-kaartje opgeplakt. Adélie, Gentoo, Chinstrap — lukraak rondgestrooid.

Dat husselen koppelt de soort los van het vin-gewicht-verband. Je haalt met eigen handen weg dat een bepaalde soort bij een bepaalde helling hoort. In zo’n losgekoppelde wereld kán de ene soort geen steilere lijn meer hebben dan de andere — als de hellingen dan tóch verschillen, is dat puur toeval van de verdeling.

In R · schud de soort-labels los

# de helling per soort, in één functietje
hellingen <- function(d) sapply(c("Adelie", "Chinstrap", "Gentoo"),
  function(s) coef(lm(gewicht ~ vinlengte, data = subset(d, soort == s)))[2])

echt <- hellingen(peng)
max(echt) - min(echt)          # → 21.79: de kloof tussen steilste en flauwste

# doe het 1000 keer met geschudde soort-labels
set.seed(42)
nep <- replicate(1000, {
  d <- peng
  d$soort <- sample(d$soort, replace = FALSE)   # soort-kaartjes lukraak rondgestrooid
  h <- hellingen(d)
  max(h) - min(h)
})

# het venster moet allebei kunnen tonen: het bergje én de streep
hist(nep, xlim = c(0, max(c(nep, 21.79)) * 1.1),
     main = "", xlab = "kloof tussen steilste en flauwste helling")
abline(v = 21.79, col = "red", lwd = 2)
sum(nep >= 21.79)              # → 0 van de 1000

Kijk naar dat plaatje. Duizend werelden waarin de soort de vin-gewicht-relatie niet modereert — want in elke geschudde wereld is de soort losgeknipt van de helling. De kloof tussen de hellingen danst rond een paar gram: soms toevallig wat groter, meestal klein. En dan die rode streep op 21.79, ver buiten het hele bergje. In duizend keer schudden haalde geen enkele geschudde wereld de echte kloof. De echte pinguïns weten iets wat de geschudde niet weten: bij de Gentoo’s telt een langere vin écht zwaarder mee.

snap het

Onder welke conditie doet iets ertoe

interactie · moderatie · verschil in verschil

Wat je hier drie keer zag, heeft één naam: een interactie, in de wandelgangen ook moderatie. En de nette formulering is de moeite waard om precies over te nemen, want de slordige versie kost je later punten: het effect van vinlengte op gewicht verandert in grootte en/of richting met het niveau van de soort. Niet “vinlengte interacteert met soort” — dat zegt niet wát er met wát gebeurt.

De kortste weg naar waar het echt over gaat, loopt even bij de pinguïns vandaan. Zou wel of geen therapie uitmaken voor het welzijn van mensen? Ja, hè. Maar ga nu eens naar een moeilijk land waar geen eten is, waar mensen verhongeren. Geef sommige van die hongerende mensen therapie en sommige geen therapie — zou dat verschil maken? Nee: als je geen eten hebt, heeft therapie ook geen zin.

Therapie geeft dus soms een verschil in welzijn, maar alleen in bepaalde situaties. Bij een interactie-effect is dát de vraag: wanneer, onder welke conditie, heeft een variabele een effect op de uitkomst?

Dáár gaat het om. Een interactie is geen ingewikkeldere versie van een effect — het is de erkenning dat de vraag “hoeveel effect heeft X op Y?” soms geen antwoord heeft zonder de vervolgvraag “bij wie, en wanneer?” Bij de pinguïns: “hoeveel extra gewicht per millimeter vin?” heeft geen antwoord zonder “bij welke soort?”

En daarom loog die 49.69 uit de prikkel niet, maar antwoordde hij op een vraag die niet goed gesteld was.

Drie dingen om vast te houden, want ze komen alle drie terug.

Eén: hoogte is geen steilheid. De sekse maakte de pinguïns 683.41 gram zwaarder en veranderde niets aan de helling. Dat is een hoofdeffect, geen interactie. Twee wolken op verschillende hoogte met parallelle lijnen — verschil, maar geen verschil in verschil.

Twee: bij meer dan twee groepen stopt de coëfficiënt met antwoorden. Twee groepen geven één interactieregel die je gewoon kunt aflezen. Drie groepen geven er twee, en dan heb je een modelvergelijking nodig om te weten of er überhaupt iets is. Dat is geen extra formaliteit; het is dezelfde reden waarom je bij drie gemiddelden een eenweg ANOVA doet en niet drie t-toetsen.

Drie: een hoofdeffect lees je aan de randen af. Bij de zes hokjes vroeg je niet naar één cel maar naar het gemiddelde van een hele rij of kolom — de marginalen. Doe het even met lengtes: in Nederland zijn de mannen 1.80 en de vrouwen 1.70, in Japan 1.70 en 1.65. Vraag je “er komt een man binnen, hoe lang is die?”, dan gebruik je het marginaal van de mannen — (1.80 + 1.70) / 2 = 1.75 — en bij een vrouw (1.70 + 1.65) / 2 = 1.675. Die twee getallen gaan enkel en alleen over het hoofdeffect van geslacht op lengte, niet over land, en het verschil ertussen is 7.5 centimeter. Zo lees je ook de pinguïn-tabel: de rijen voor sekse, de kolommen voor soort, en de manier waarop de hokjes van hun eigen randen afwijken is de interactie.

jouw beurt

Dezelfde zes hokjes, een andere meetlat

tweeweg-ANOVA · type III · parallel

Jouw beurt. Je hebt net gezien dat sekse en soort samen iets doen met het gewicht van een pinguïn: het sekseverschil is groter bij de Gentoo’s dan bij de Chinstraps. Nu dezelfde twee groeperingen, dezelfde zes hokjes — maar je meet iets anders: de snavellengte. De vraag: is het sekseverschil in snavellengte ook niet overal even groot?

Zet eerst klaar, en beantwoord daarna de vragen. De laatste heet met opzet Z) en niet l): de aannames staan altijd achteraan, en die vreemde letter zegt “en ook, elke keer, deze”.

In R · zet klaar

ps <- subset(peng, sekse %in% c("male", "female"))

options(contrasts = c("contr.sum", "contr.poly"))              # 1. aanzetten
model.snavel <- lm(snavellengte ~ sekse * soort, data = ps)    # 2. rekenen
car::Anova(model.snavel, type = 3)
lsr::etaSquared(model.snavel, type = 3)
options(contrasts = c("contr.treatment", "contr.poly"))        # 3. terugzetten

Jouw beurt

a) Hoeveel pinguïns zitten er in de zes hokjes, en hoeveel per hokje? Gebruik table(ps$sekse, ps$soort). Waarom is het totaal geen 342? (zes getallen, één totaal, en één zin)

b) Reken het gemiddelde per hokje uit met tapply(ps$snavellengte, list(ps$sekse, ps$soort), mean) en zet ze in een tabel van twee rijen bij drie kolommen. (zes getallen)

c) Trek per soort de vrouwtjes van de mannetjes af. Zijn die drie kloven ongeveer even groot? (drie getallen, en hoeveel de grootste en de kleinste schelen)

d) Schrijf de nulhypothese voor de interactie voluit op, in woorden én in symbolen. (één zin, en die ook in symbolen)

e) Toets de interactie. (rapporteer F(df₁, df₂) en p)

f) Rapporteer nu ook de twee hoofdeffecten. (per hoofdeffect F(df₁, df₂) en p)

g) Reken de drie effectgroottes uit met lsr::etaSquared(model.snavel, type = 3). Welke van de drie doet het meeste werk? (drie η²’s)

h) Laat de contrast-regel nu eens weg en draai car::Anova(..., type = 3) opnieuw. Welke F verandert er, en welke niet? (twee getallen voor sekse, één voor de interactie, en één zin over waaróm de interactie stilstaat)

i) Teken het met interaction.plot(ps$soort, ps$sekse, ps$snavellengte, xlab = "soort", ylab = "gemiddelde snavellengte (mm)", trace.label = "sekse") en leg het naast het eerdere plaatje van het gewicht. (twee plaatjes, en wat je aan de lijnen ziet)

j) Waarom moet je in je rapportagezin schrijven dát het type III-kwadratensommen zijn? (één zin, met het getal uit h) erin)

k) Bij het gewicht wás er een interactie, bij de snavellengte niet — op precies dezelfde 333 pinguïns, met precies dezelfde twee groeperingen. Wat zegt dat over waar een interactie in zit?

Z) Mag je alles hierboven geloven? Loop de aannames na: is de spreiding in snavellengte ongeveer gelijk over de zes hokjes (tapply(ps$snavellengte, list(ps$sekse, ps$soort), sd)), en telt elke pinguïn maar één keer mee? En dan de vraag die bij dít ontwerp hoort: kijk nog eens naar je zes celaantallen uit a) — wat zou er met de hele type I/type III-kwestie gebeuren als die zes gelijk waren geweest?

Antwoord — open na je eigen poging

a) Vrouwtjes 73, 34 en 58; mannetjes 73, 34 en 61 — samen 333. Geen 342, want bij negen pinguïns is de sekse nooit vastgesteld, en subset() laat die hier expliciet weg. Onthoud die zes getallen; bij Z) blijken ze ertoe te doen.

	Adélie	Chinstrap	Gentoo
vrouwtje	37.26	46.57	45.56
mannetje	40.39	51.09	49.47

c) 3.13, 4.52 en 3.91 mm. Het sekseverschil is er bij alle drie de soorten, en het scheelt tussen de grootste en de kleinste maar 1.39 mm. Vergelijk dat met het gewicht, waar diezelfde drie kloven 674.7, 411.8 en 805.1 gram waren — daar was de grootste bijna twee keer de kleinste.

d) De nulhypothese: het sekseverschil in snavellengte is bij alle drie de soorten even groot. In symbolen: $H_0:\ \mu_{\text{man},j} - \mu_{\text{vrouw},j}$ is hetzelfde voor elke $j$, waarin $j$ een aftelnummer voor de soort is. Je mag het ook schrijven als alle interactietermen in de populatie zijn nul; één van de twee is genoeg.

e) F(2, 327) = 2.28, p = .103. Niet significant: deze keer lopen de lijnen zo goed als parallel.

f) Sekse F(1, 327) = 208.06, p < .001. Soort F(2, 327) = 650.21, p < .001. Allebei ruimschoots.

g) η² = .11 voor sekse, .70 voor soort en .00 voor de interactie. De soort doet vrijwel al het werk; dat mannetjes langere snavels hebben komt daar als een keurig parallelle laag bovenop.

h) Voor sekse: 208.06 mét de contrast-regel en 66.81 zonder — minder dan een derde. Voor de interactie in allebei de gevallen 2.28. De interactie staat stil omdat hij als laatste komt, wat je ook doet: er valt voor hem niets meer te verdelen. Precies wat je bij het gewicht in de vier-varianten-tabel zag.

i) Bij de snavellengte lopen de twee lijnen zo goed als parallel: overal ongeveer vier millimeter tussen de vrouwtjes en de mannetjes. Bij het gewicht liepen ze zichtbaar uiteen — bij de Gentoo’s stond het sekseverschil veel wijder open dan bij de Chinstraps. Hoogte tegenover steilheid, en hier is het alleen hoogte.

j) Omdat dezelfde data zonder die vermelding niet na te rekenen is: je kreeg bij h) 208.06 en 66.81 op precies dezelfde 333 pinguïns. Wie niet weet welke van de twee je gebruikte, kan je F niet reproduceren en dus ook niet controleren.

k) Dat een interactie niet in de twee groeperingen zit maar in wat je meet. Sekse en soort zijn in allebei de analyses dezelfde twee kolommen, op dezelfde 333 dieren. Het gewicht van een Gentoo-mannetje loopt anders uit de pas dan zijn snavel. “Sekse en soort interacteren” is dus een onaffe zin — de vraag is altijd: interacteren ze op welke uitkomst?

Z) De spreiding per hokje loopt van 1.56 mm (Chinstrap-mannetjes) tot 3.11 mm (Chinstrap-vrouwtjes) — de grootste is 1.99 keer de kleinste. Dat is krap onder de twee, dus het mag, maar het is geen ruime marge; bij een grensgeval als p = .103 is dat het vermelden waard. Elke pinguïn is één keer gemeten, dus de waarnemingen zijn onafhankelijk.

En dan de tweede helft van de vraag, want dáár zit de winst. Waren de zes hokjes even groot geweest, dan had de hele type I/type III-kwestie niet bestaan. De twee cirkels overlappen alleen doordat de aantallen ongelijk zijn; bij gelijke aantallen valt er niets te verdelen en geven alle varianten hetzelfde getal. Op deze data is dat na te meten: trek 34 dieren uit elk hokje en type I met sekse eerst, type I met soort eerst en type III komen alle drie op dezelfde F uit. Op de échte, scheve verdeling geven ze 219.23, 211.81 en 208.06.

De weg terug, voor als de spreiding wél ver uiteen was gelopen:

Blijft staan: de zes celgemiddelden bij b) en de drie kloven bij c). Dat zijn beschrijvingen van deze 333 dieren.
Doe je over: alle drie de F’s bij e) en f), en daarmee je conclusie bij k). Die rusten op één gedeelde restspreiding over de zes hokjes, en die is er dan niet.

rapportage

Schrijf het op als een onderzoeker

rapportagezin · F en p · nul-vondst

Rapportage

Een interactie opschrijven klinkt als iets voor de gevorderden, maar het blijft een boodschappenlijstje in zinsvorm: wat vergeleek je, verschilde het, en waar zat het. Twee modelzinnen, één voor elke vorm die je vandaag zag.

Een helling die per groep verschilt:

Om te onderzoeken of het effect van vinlengte op gewicht afhangt van de soort, vergeleken we bij 342 pinguïns twee modellen: één waarin alle soorten dezelfde helling moesten delen, en één waarin elke soort een eigen helling mocht hebben. De hellingen verschilden significant tussen de soorten, F(2, 336) = 5.53, p = .004: bij de Gentoo’s ging een millimeter extra vin samen met meer extra gewicht (54.62 g/mm) dan bij de Adélies (32.83) en de Chinstraps (34.57). Het effect van vinlengte op gewicht werd dus gemodereerd door de soort.

Twee groeperingen tegelijk:

Om te onderzoeken of het effect van sekse op gewicht afhangt van de soort, voerden we bij 333 pinguïns een tweeweg-variantieanalyse uit met sekse en soort als factoren (type III-kwadratensommen). De interactie was significant, F(2, 327) = 8.76, p < .001, η² = .01: het sekseverschil was groter bij de Gentoo’s (805.09 g) dan bij de Chinstraps (411.76 g). Beide hoofdeffecten waren eveneens significant, voor sekse F(1, 327) = 311.84, p < .001, η² = .14, en voor soort F(2, 327) = 746.92, p < .001, η² = .66.

Kijk wat er cursief staat en wat niet. Alléén de symbolen F en p hellen; η² is een Griekse letter en blijft rechtop, en de gewone woorden en de gewichtsgetallen ook. De statistieken krijgen twee decimalen, APA-precies, terwijl de lopende tekst hierboven gerust rondde (“ruim zeshonderd gram”) voor de intuïtie. Een p die kleiner is dan .001 schrijf je niet uit maar als p < .001, en zonder nul ervoor. Alles staat in de verleden tijd, want je vertelt wat er gevónden werd, en beide zinnen staan op eigen benen: wie er één uit dit blok tilt, weet nog steeds wie er gemeten werd, wat er gemeten werd, en waar de interactie zat.

In geen van beide zinnen staat welk programma je gebruikt hebt. Dat is met opzet: een methodezin zegt wát je gedaan hebt, niet welke knop je hebt aangeklikt of welke regel je hebt getypt. Zag jij zelf maar één tabel in plaats van twee modellen? Dan legde je programma die twee voor je naast elkaar, en het getal is hetzelfde. Twee onderzoekers die dezelfde pinguïns dezelfde vraag stellen, hebben niet twee verschillende methoden gebruikt — en precies díé verwarring hoort een methodesectie weg te nemen.

Merk ook op dat de tweede zin zégt dat het type III-kwadratensommen zijn. Dat is geen pronken met vaktaal — je zag hierboven dat dezelfde data vier verschillende F-waarden kan opleveren, dus zonder die vermelding kan niemand jouw getal narekenen.

Jouw beurt: schrijf zelf de rapportage-zin voor de taak hierboven — de snavellengte. Die interactie was juist níét significant (F(2, 327) = 2.28, p = .10), en ook een nul-vondst rapporteer je netjes op: dat het sekseverschil in snavellengte bij alle drie de soorten ongeveer even groot was, is een echt resultaat. De modelzinnen hierboven mag je lenen; hardop voorlezen mag ook.

◠ In het boek

Hier: als de moderator een groep is, kun je het tekenen — twee lijnen, drie lijnen, of zes hokjes.

In het boek: W10 — Verschil in verschil

Daar draait het om de lijnen die niet parallel lopen — het moment waarop je merkt dat “het effect” een tweede vraag nodig heeft: bij wie?

Voetnoten

Palmer Penguins — Gorman, Williams & Fraser (2014), Palmer Station Antarctica LTER; via het R-pakket palmerpenguins (CC0). Alle getallen in dit blok zijn op die échte data nagerekend.↩︎