Vierde beweging · Verschil onthechten

Oefening 10.2 · Als de moderator een getal is

Interactie tussen twee continue voorspellers · continuous moderation

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

prikkel

De moderator heeft geen hokjes meer

continue moderator

In oefening 10.1 kon je het altijd tekenen. Twee groepen, twee lijnen. Drie soorten, drie lijnen. Sekse naast soort, zes hokjes met een punt in elk hokje. Wat de moderator ook was, je kon hem aanwijzen en de hellingen ervan aflezen.

Probeer dat nu eens met de snavellengte.

Dezelfde 342 pinguïns.1 Je wilt weten of het effect van vinlengte op gewicht afhangt van hoe lang de snavel is. Redelijke vraag — een pinguïn met een lange snavel is misschien een ander soort dier dan een met een korte, en bij een ander dier telt een millimeter vin misschien anders.

Maar teken het maar eens. Snavellengte heeft geen twee waarden en geen drie; hij loopt van 32.1 tot 59.6 millimeter, en in deze 342 beesten komen er 164 verschillende waarden voor. Eén lijn per snavellengte worden er dan 164. Er is geen hokje om in te gaan staan.

En tóch is er iets. Knip de pinguïns straks maar eens grof in tweeën bij de middelste snavellengte — dat doen we verderop echt — dan ligt de helling van gewicht op vinlengte bij de korte-snavelhelft rond de 37 en bij de lange-snavelhelft rond de 57. Dus of het verband tussen vinlengte en gewicht meebeweegt met de snavellengte, of het toeval van deze 342 beesten is dat.

Deze oefening gaat over die twee vragen tegelijk: hoe je zoiets tóétst als er geen groepen zijn, en — lastiger — hoe je het opschrijft als er geen groepen zijn om je hellingen aan op te hangen.

speel het

Waar de nul woont, en waarom dat nu wél uitmaakt

centreren · intercept

OpmerkingEerst zoals in het boek — de twaalf diamantjes
d <- read.csv("diamantjes.csv")

# twee getallen als voorspellers, en hun product
model.gl.k <- lm(glans ~ gladheid * karaat, data = d)
summary(model.gl.k)
# gladheid          b = 10.00
# karaat            b =  4.70
# gladheid:karaat   b = -1.21   p = .80

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.

En hier moet ik meteen eerlijk zijn: op twaalf steentjes valt over deze interactie niets te zeggen. p = .80, en dat blijft zo hoe je ook kijkt. Een interactie tussen twee getallen is de hongerigste term die er in een model kan zitten — hij vraagt om veel meer waarnemingen dan een hoofdeffect, want hij moet een verschil in een verschil zichtbaar maken. Twaalf steentjes zijn genoeg om een lijn te trekken en te krap om een waaier te zien. Dat is geen tegenvaller maar een les die je meeneemt.

Kijk wél naar dat eerste getal, want dáár leren de diamantjes je iets wat de pinguïns straks nodig hebben. gladheid staat op 10.00. Maar in oefening 10.1 stond de helling van gladheid op 6.14, en in de meervoudige regressie weer op iets anders. Hoe kan dat?

Omdat een coëfficiënt in een model met een interactie altijd de helling is bij nul van de ander. Die 10.00 is de helling van glans op gladheid bij een steentje van nul karaat. En een steentje van nul karaat bestaat niet — het gemiddelde steentje weegt 1.6 karaat. Je leest dus een helling af op een adres waar niemand woont.

De reparatie is één aftrekking:

# de ruwe variabele min het gemiddelde: dat is centreren, meer niet
d$glad_c <- d$gladheid - mean(d$gladheid)   # gemiddelde gladheid = 5.0
d$kar_c  <- d$karaat   - mean(d$karaat)     # gemiddelde karaat   = 1.6

model.gl.k.c <- lm(glans ~ glad_c * kar_c, data = d)
summary(model.gl.k.c)
# glad_c          b = 8.06    <- helling bij GEMIDDELDE karaat
# glad_c:kar_c    b = -1.21   <- exact hetzelfde als hierboven
library(dplyr)
d <- d |> mutate(glad_c = gladheid - mean(gladheid),
                 kar_c  = karaat   - mean(karaat))

model.gl.k.c <- lm(glans ~ glad_c * kar_c, data = d)
summary(model.gl.k.c)

Twee dingen om vast te houden. De interactie verandert niet — hij was −1.21 en hij blijft −1.21. Centreren verschuift alleen waar de nul ligt; het verzet het model niet. De helling van gladheid verandert wél, van 10.00 naar 8.06, en dat komt doordat je hem nu op een ander adres afleest: niet bij nul karaat maar bij een gemiddeld steentje.

Nul is hier gemiddeld. Kijk je naar een moderatieregressie, dan staat het effect van de voorspeller op de uitkomst in de coëfficiëntentabel altijd bij een nulwaarde van de moderator. De vraag is dus: wat betekent nul?

Onthoud die vraag. Vanaf hier is het de vraag die je bij elke tabel stelt.

De waaier bij de pinguïns

gecentreerde voorspellers · interactieterm

Nu jij — dezelfde beesten als in 10.1, en nu twee getallen als voorspellers. Centreer ze allebei voordat je het model draait.

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

peng$vin_c  <- peng$vinlengte    - mean(peng$vinlengte)     # gemiddelde: 200.92 mm
peng$snav_c <- peng$snavellengte - mean(peng$snavellengte)  # gemiddelde:  43.92 mm

round(c(mean(peng$vin_c), mean(peng$snav_c)), 10)   # → 0 en 0, zoals het hoort
library(dplyr)
peng <- read.csv("penguins.csv") |>
  mutate(vin_c  = vinlengte    - mean(vinlengte),
         snav_c = snavellengte - mean(snavellengte))
OpmerkingIn R · trek eraan
model.vl.sn.c <- lm(gewicht ~ vin_c * snav_c, data = peng)
summary(model.vl.sn.c)

# (Intercept)    b = 4141.49   <- gewicht van een dóórsnee pinguïn
# vin_c          b =   45.39   <- helling bij GEMIDDELDE snavellengte
# snav_c         b =   11.82   p = .027
# vin_c:snav_c   b =    1.20   p < .001
# Multiple R-squared: 0.769

Lees die tabel nu met de vraag van hierboven in je hoofd — wat betekent nul?

  • De intercept, 4141.49 gram, is het voorspelde gewicht van een pinguïn met een gemiddelde vin én een gemiddelde snavel. Dat is een pinguïn die kán bestaan. Vergelijk dat met oefening 3.2, waar de intercept op −5780.83 gram stond: een pinguïn met een vin van nul millimeter, die bijna zes kilo negatief weegt. Het adres is nu bewoond.
  • 45.39 is de helling van gewicht op vinlengte bij een pinguïn met een gemiddelde snavel — niet bij álle pinguïns, bij een gemiddelde. En let op: in oefening 3.2 stond hier 49.69. Dat verschil komt niet door het centreren maar doordat de snavellengte nu in het model zit: 45.39 is de helling van gewicht op vinlengte bij pinguïns die even lang van snavel zijn. Precies de verschuiving die je in oefening 8.1 bij de meervoudige regressie zag.
  • En 1.20 is de interactie: per millimeter extra snavellengte komt er 1.20 gram per millimeter bij de helling van de vin bij.

Dat laatste getal is klein en je zou eroverheen kunnen lezen. Maar het is een helling van een helling, en dat is een ander soort getal. Reken maar mee.

Kies je eigen hokjes: de simple slopes

simple slopes · min één SD, gemiddelde, plus één SD

Er zijn geen groepen, dus de hellingen komen er niet vanzelf uit. Je moet ze zelf gaan halen — en dat betekent: je kiest waarden van de moderator waarop je de helling wilt aflezen. De afspraak die je overal tegenkomt is: één standaardafwijking eronder, het gemiddelde, en één standaardafwijking erboven.

De standaardafwijking van de snavellengte is 5.46 millimeter. Dus je leest de helling af bij een snavel van 38.46, van 43.92 en van 49.38 millimeter — kort, gemiddeld, lang.

En dan is het optellen, net als bij de dummy in 10.1:

sd(peng$snav_c)                    # → 5.46

45.39 - 1.20 * 5.46                # → 38.84   bij een korte snavel
45.39                              # → 45.39   bij een gemiddelde snavel
45.39 + 1.20 * 5.46                # → 51.94   bij een lange snavel
OpmerkingIn R · dezelfde drie hellingen, uitgerekend uit het model
k <- coef(model.vl.sn.c)
s <- sd(peng$snav_c)

k["vin_c"] + k["vin_c:snav_c"] * c(-s, 0, s)      # → 38.84  45.39  51.94

38.84 · 45.39 · 51.94. Bij een pinguïn met een korte snavel weegt elke millimeter vin bijna negenendertig gram; bij een lange snavel bijna tweeënvijftig. Dertien gram per millimeter verschil tussen de buitenste twee. Hoe langer de snavel, hoe zwaarder elke millimeter vin telt.

En dat zijn geen echte groepen — er staat nergens een pinguïn met het bordje “lange snavel” om zijn nek. Het zijn drie plekken op een glijdende schaal, door jou gekozen omdat je ergens moest gaan staan. Kies je andere waarden, dan krijg je andere hellingen; het model is hetzelfde gebleven.

Wil je zien dat het niet uit de lucht komt vallen, knip de pinguïns dan even grofweg in tweeën bij de middelste snavellengte en trek twee losse lijnen:

OpmerkingIn R · ter controle, grof in tweeën geknipt
middelste <- median(peng$snavellengte)

model.kort <- lm(gewicht ~ vinlengte, data = subset(peng, snavellengte <  middelste))
model.lang <- lm(gewicht ~ vinlengte, data = subset(peng, snavellengte >= middelste))

coef(model.kort)[2]   # → 37.17
coef(model.lang)[2]   # → 57.27

37.17 tegen 57.27 — dezelfde richting, en forser, want een mediaan-knip zet de uitersten tegenover elkaar in plaats van één standaardafwijking. Doe dat níét als je het echt gaat rapporteren: je gooit de helft van de informatie weg en je verzint groepen die niet bestaan. Als controle bij het leren is het prima.

Waarom centreren hier geen gemak meer is

intercept · onmogelijke nul

In oefening 3.2 was centreren een leuk extraatje. Hier is het dat niet meer, en dat kun je in één regel zien. Draai hetzelfde model nog eens, nu zonder te centreren:

OpmerkingIn R · dezelfde data, ongecentreerd
model.vl.sn <- lm(gewicht ~ vinlengte * snavellengte, data = peng)
summary(model.vl.sn)

# vinlengte                b =   -7.31   p = .63
# snavellengte             b = -229.24   p < .001
# vinlengte:snavellengte   b =    1.20   p < .001   <- exact hetzelfde

Kijk naar de eerste regel. De helling van gewicht op vinlengte is nu −7.31. Negatief. Alsof een langere vin een pinguïn lichter maakt.

Dat is geen rekenfout en het spreekt oefening 3.2 niet tegen. Het is de helling van gewicht op vinlengte bij een snavellengte van nul millimeter — bij een pinguïn zonder snavel. Zo’n dier bestaat niet, de data bevat er geen enkele in de buurt van, en het model verzint daar dus maar wat. Je leest een echt getal af op een onmogelijk adres.

De interactie blijft ondertussen exact 1.20, gecentreerd of niet. Dát is het patroon: centreren raakt de hoofdeffecten en laat de interactie met rust. En precies daarom kun je het niet overslaan — het getal dat je wél mag aflezen verandert van betekenis, en niets in de uitvoer waarschuwt je.

schud het

Duizend werelden waarin de snavel niets modereert

interactieterm schudden · dansen rond nul

Is die 1.20 echt, of kan het toeval hem namaken? We spelen weer vals. Deze keer husselen we de snavellengtes los: elke pinguïn houdt z’n eigen vinlengte en z’n eigen gewicht, maar krijgt de snavel van een willekeurige andere pinguïn toebedeeld.

Dat husselen koppelt de snavellengte los van het vin-gewicht-verband. Je haalt met eigen handen weg dat een bepaalde snavellengte bij een bepaalde helling hoort. In zo’n losgekoppelde wereld kán de snavellengte niets meer modereren — de ware interactie is er nul. Wat er dan tóch uitkomt, is puur toeval.

OpmerkingIn R · schud de snavels los
echt <- coef(model.vl.sn.c)["vin_c:snav_c"]
echt                            # → 1.1998

set.seed(42)
nep <- replicate(1000, {
  d <- peng
  d$snav_c <- sample(d$snav_c, replace = FALSE)  # snavels lukraak herverdeeld
  coef(lm(gewicht ~ vin_c * snav_c, data = d))["vin_c:snav_c"]
})

mean(nep)                       # → -0.002   de dans zit om nul
sd(nep)                         # →  0.269
range(nep)                      # → -0.76 tot 0.94
sum(abs(nep) >= abs(echt))      # → 0 van de 1000

# het venster moet allebei kunnen tonen: het bergje én de streep
hist(nep, xlim = c(-1.4, 1.4), main = "", xlab = "interactie in een geschudde wereld")
abline(v = echt, col = "red", lwd = 2)

Duizend werelden waarin de snavellengte de vin-gewicht-relatie niet modereert. Het gemiddelde van die duizend is −0.002 — zo goed als nul, precies zoals het hoort, want in elke geschudde wereld ís de waarheid nul. En tóch landt geen enkele wereld precies op nul: ze dansen eromheen, met een spreiding van 0.27, en de wildste haalde 0.94.

En dan die rode streep op 1.1998, buiten het hele bergje. Geen van de duizend kwam er in de buurt. De echte pinguïns weten iets wat de geschudde niet weten: bij een langere snavel telt een millimeter vin écht zwaarder mee.

Kijk voor de aardigheid ook eens wat er in zo’n geschudde wereld met de drie simple slopes gebeurt. Die komen alle drie rond de 49.6 uit te liggen — de gewone helling van gewicht op vinlengte over alle pinguïns, zonder waaier. Als er niets te moduleren valt, valt de waaier dicht tot één lijn.

snap het

Een vlak dat kantelt

continue moderatie · symmetrie van de interactie

Bij een groep als moderator kon je lijnen tekenen. Hier niet, en het is de moeite waard om te weten wat je dan wél voor je moet zien.

Met één voorspeller heb je een lijn. Met twee voorspellers en géén interactie heb je een vlak — een schuin dak boven de wolk, dat overal even schuin is. De interactie is wat dat dak laat kantelen: aan de kant van de lange snavels loopt het steiler op dan aan de kant van de korte. Geen twee lijnen die uiteenwaaieren, maar één oppervlak dat verdraait.

En daar zit iets in verstopt wat makkelijk te missen is. Een interactie is symmetrisch. Wij hebben hem gelezen als “hoe langer de snavel, hoe zwaarder elke millimeter vin telt”, maar je mag hem net zo goed omdraaien: “hoe langer de vin, hoe zwaarder elke millimeter snavel telt”. Het is één getal en het staat op één regel; welke van de twee de moderator heet, is jouw keuze en niet die van het model. Bij een groep als moderator voelde dat vanzelfsprekend — soort modereert vinlengte, niet andersom — maar dat kwam alleen doordat een soort geen helling kán hebben. Tussen twee getallen is die asymmetrie weg.

Kies dus bewust. Wat je de moderator noemt, bepaalt welk verhaal je vertelt, en je moet dat verhaal kunnen verdedigen met iets buiten de statistiek om: theorie, of de volgorde waarin de dingen gebeuren.

Wat dit hoofdstuk bij elkaar houdt. In 10.1 was de moderator een groep en kon je het tekenen; de hellingen rolden uit de tabel omdat elke groep er een had. Hier is de moderator een getal, is er geen enkele groep om een helling van af te lezen, en moet je zélf de plekken kiezen waarop je kijkt. Dat is de hele stap. Alles eromheen — de interactieterm, het product, de vraag of het verschil in verschil groot genoeg is — is hetzelfde gebleven.

En de vraag eronder ook. Denk aan de therapie in het land waar geen eten is: het gaat er niet om óf vinlengte effect heeft op gewicht, maar onder welke conditie, en hoeveel. Alleen is de conditie nu geen hokje meer maar een schaal waarop je zelf een plek moet aanwijzen.

jouw beurt

Een tweede moderator, die de waaier de andere kant op zet

continue moderatie · negatieve interactie

Jouw beurt. Je hebt gezien dat de snavellengte het effect van vinlengte op gewicht versterkt: hoe langer de snavel, hoe zwaarder elke millimeter vin telt. Er staat nog een kolom in het bestand die je niet gebruikt hebt: de snaveldiepte — hoe hoog de snavel is, van boven naar onder. Dezelfde vraag, dezelfde uitkomst, dezelfde vin-voorspeller, andere moderator.

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”.

OpmerkingIn R · zet klaar
peng$dp_c <- peng$snaveldiepte - mean(peng$snaveldiepte)
model.vl.dp.c <- lm(gewicht ~ vin_c * dp_c, data = peng)
summary(model.vl.dp.c)
OpmerkingJouw beurt

a) Hoe diep is een gemiddelde pinguïnsnavel, en hoeveel spreiding zit daarop? (rapporteer M en SD, in millimeters)

b) Lees de vier getallen uit de tabel af, en zeg bij elk wát het is — dus niet alleen het getal maar ook op welk adres het gemeten is. (vier b’s, elk met één regel uitleg)

c) Kijk naar het teken van de interactie en zeg in gewone woorden wat er staat, vóórdat je verder rekent. Hij is negatief — wat betekent dat voor de helling van gewicht op vinlengte bij een diepe snavel?

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 b, SE, t, p en het 95%-interval)

f) Op welke drie waarden van de snaveldiepte ga je de helling aflezen, en waarom juist die drie? (drie getallen in millimeters, en één zin)

g) Reken de drie simple slopes uit met k <- coef(model.vl.dp.c); s <- sd(peng$dp_c); k["vin_c"] + k["vin_c:dp_c"] * c(-s, 0, s). (drie hellingen in gram per millimeter)

h) Hoeveel keer zwaarder telt een millimeter vin bij een ondiepe snavel dan bij een diepe? (één verhouding)

i) Leg deze waaier naast die van de snavellengte uit dit blok. Wat is er hetzelfde en wat is er tegengesteld? (twee reeksen van drie, en één zin)

j) Een interactie is symmetrisch. Schrijf jouw bevinding daarom nog een keer op, maar dan met de vinlengte als moderator. (één zin, en hetzelfde getal)

k) Draai het model nu nog eens zonder te centreren, met lm(gewicht ~ vinlengte * snaveldiepte, data = peng). Welk getal blijft precies gelijk, en welk getal verandert onherkenbaar? (twee getallen, en waarom dat zo is)

l) In oefening 10.1 bleek de soort een sterke moderator — en de Gentoo’s hebben nu net een opvallend ondiepe snavel. Zou de snaveldiepte hier zelf de moderator kunnen zijn, of staat hij in voor iets anders? Wat zou je moeten doen om die twee uit elkaar te trekken?

Z) Mag je alles hierboven geloven? Bij een continue moderator komt er een aanname bij die je bij een groep niet had: je hebt de hellingen afgelezen op drie waarden die je zélf gekozen hebt. Wonen daar pinguïns? Kijk naar range(peng$snaveldiepte). En loop daarna de gewone drie na: liggen de residuen in een egale band, liggen ze symmetrisch, en telt elke pinguïn maar één keer mee? (één bereik en drie uitspraken, elk met het getal waar je hem op baseert)

a) M = 17.15 mm, SD = 1.97 mm.

b) Vier getallen, vier adressen:

  • Intercept 4062.79 — het voorspelde gewicht van een pinguïn met een gemiddelde vin én een gemiddelde snaveldiepte. Een dier dat kán bestaan, want je hebt allebei gecentreerd.
  • vin_c 48.63 — de helling van gewicht op vinlengte bij een pinguïn met een gemiddeld diepe snavel. Niet bij álle pinguïns; bij een gemiddelde.
  • dp_c 44.64 — de helling van gewicht op snaveldiepte bij een pinguïn met een gemiddelde vin.
  • Interactie −8.60 — een helling van een helling: per millimeter diepere snavel gaat er 8.60 gram per millimeter af van de helling van de vin.

c) Negatief betekent: hoe diéper de snavel, hoe mínder elke millimeter vin meetelt voor het gewicht. Bij de snavellengte was het andersom. De waaier kantelt hier dus de andere kant op.

d) De nulhypothese: het effect van vinlengte op gewicht hangt niet af van de snaveldiepte. In symbolen: \(H_0:\ \beta_3 = 0\) tegen \(H_1:\ \beta_3 \neq 0\), waarin \(\beta_3\) het gewicht van de interactieterm in de populatie is. Je mag ook zeggen dat de helling van de vin voor elke snaveldiepte dezelfde is; één van de twee is genoeg. Grieks, want het gaat over de populatie en niet over deze 342.

e) b = −8.60, SE = 1.34, t(338) = −6.41, p < .001, 95% CI [−11.23, −5.96]. Het interval blijft ruim aan één kant van de nul.

f) Bij 15.18, 17.15 en 19.13 mm: één standaardafwijking onder het gemiddelde, het gemiddelde zelf, en één erboven. Dat is een afspraak en geen wet — er staat nergens een pinguïn met het bordje “ondiepe snavel” om zijn nek. Je moest ergens gaan staan, en dit is waar iedereen gaat staan.

g) 65.61 · 48.63 · 31.66 gram per millimeter.

h) 65.61 / 31.66 = 2.07 keer. Bij een ondiepe snavel telt een millimeter vin dus ruim tweemaal zo zwaar als bij een diepe.

i) Hetzelfde: allebei de moderatoren zijn een getal, allebei op precies dezelfde 342 pinguïns, met dezelfde uitkomstmaat en dezelfde voorspeller. Tegengesteld: bij de snavellengte loopt de reeks omhoog (38.84 → 45.39 → 51.94), hier omlaag (65.61 → 48.63 → 31.66). Twee moderatoren die de waaier de tegengestelde kant op kantelen.

j) “Per millimeter langere vin gaat er 8.60 gram per millimeter af van de helling van de snaveldiepte.” Precies hetzelfde getal, precies dezelfde regel in de tabel — alleen jouw verhaal is omgedraaid. Welke van de twee je de moderator noemt, is jouw keuze en niet die van het model.

k) Gelijk blijft de interactie: −8.60, tot op de decimaal. Onherkenbaar verandert vinlengte: van 48.63 naar 196.07. Dat komt doordat je hem nu afleest bij een snaveldiepte van nul millimeter — bij een pinguïn zónder snavel. Zo’n dier bestaat niet, de data bevat er geen enkele in de buurt van, en het model verzint daar maar wat. Dat is het patroon van dit blok: centreren verzet alleen waar de nul ligt, dus het raakt de hoofdeffecten en laat de interactie met rust.

l) Je weet het op dit punt nog niet, en dat is het eerlijke antwoord. De Gentoo’s hebben gemiddeld een snaveldiepte van 14.98 mm tegen 18.37 bij de rest, dus “ondiepe snavel” en “Gentoo” zijn bijna hetzelfde bericht. Om ze uit elkaar te trekken moet je één van de twee vasthouden: draai de interactie binnen één soort, want daar valt aan de soort niets meer te variëren. Dat doe je in de Eindopdracht hieronder, bij k). Dit is de vraag uit oefening 8.2 in een nieuw jasje: een verband dat verdwijnt zodra je een derde variabele in het model zet die met allebei samenhangt.

Z) Eerst de nieuwe aanname, want die is de belangrijkste. De snaveldiepte loopt in deze data van 13.1 tot 21.5 mm, en je leest af op 15.18, 17.15 en 19.13 — alle drie ruim binnen dat bereik. Er zitten zelfs 23 pinguïns binnen een kwart millimeter van je onderste afleespunt en 29 binnen een kwart millimeter van je bovenste. Je leest je hellingen dus af op adressen waar écht dieren wonen. Was een van de drie buiten het bereik gevallen, dan had je een helling gerapporteerd voor een pinguïn die niet bestaat — precies de fout die de ongecentreerde 196.07 bij k) maakt.

De gewone drie houden ook stand. De residuen liggen in een egale band: bij de lage helft van de voorspelde waarden is hun standaardafwijking 359.85 gram, bij de hoge helft 381.21 — een verhouding van 1.06, geen trechter. Ze liggen bijna symmetrisch (scheefheid 0.37). En elke pinguïn is één keer gewogen, dus de waarnemingen zijn onafhankelijk.

Eén ding om op te merken en niet af te doen als niets: vinlengte en snaveldiepte hangen onderling samen, r = −.58. Dat is genoeg om de standaardfouten wat op te blazen, maar niet genoeg om het model onbruikbaar te maken.

De weg terug, voor als de residuen wél getrechterd hadden:

  • Blijft staan: de vier b’s bij b) en de drie simple slopes bij g). Dat zijn beschrijvingen van deze 342 pinguïns.
  • Doe je over: de SE, de t, de p en het interval bij e). Die rusten alle vier op één gedeelde restspreiding, en die is er dan niet.

eindopdracht

Van getal naar groep — en weer terug

dummy · simple slopes · modelvergelijking · R²

De afsluiter van het hele hoofdstuk. In de vorige oefening was de moderator een groep en rolden de hellingen uit de tabel; in deze oefening was hij een getal en moest je zelf kiezen waar je ging staan. Die twee hebben elkaar nog niet ontmoet. Dat gaan we nu regelen, op dezelfde 342 pinguïns en met dezelfde vraag — hangt het effect van vinlengte op gewicht af van de snavellengte? Want een getal kun je in groepen knippen, en een groep kun je soms als getal meten. Deze oefening zei dat dat eerste zonde is. Dat gaan we niet geloven maar narekenen. En let op de verrassing die eraan zit te komen.

Een moderator met maar twee waarden is de meest simpele situatie, want dan kun je weer in twee groepen denken. Wat altijd fijn is. Gun jezelf dat plezier — en kijk daarna wat het kost.

En let op de vorm: een artikel valt in twee gelabelde stukken, en die volgorde is niet dezelfde als de volgorde waarin je het leerde. Toen mocht de aannamecontrole achteraf, want je moest eerst zien wát je nou eigenlijk controleerde. Nu je het opschrijft voor een lezer gaat hij vóórop, want die wil weten of hij je getallen mag geloven vóórdat hij ze leest.

OpmerkingIn R · zet de vier gedaanten klaar
# 1. het getal (staat er al uit de romp van dit blok)
model.vl.sn.c <- lm(gewicht ~ vin_c * snav_c, data = peng)

# 2. geknipt in tweeën, bij de middelste snavel
peng$lang <- as.numeric(peng$snavellengte >= median(peng$snavellengte))
model.vl.lang <- lm(gewicht ~ vin_c * lang, data = peng)

# 3. geknipt in drieën, in vrijwel even grote klassen
grens <- quantile(peng$snavellengte, c(1/3, 2/3))
peng$klasse <- cut(peng$snavellengte, c(-Inf, grens, Inf),
                   labels = c("kort", "midden", "lang"))
model.klasse.zonder <- lm(gewicht ~ vin_c + klasse, data = peng)
model.klasse.met    <- lm(gewicht ~ vin_c * klasse, data = peng)

# 4. de echte soort, uit de vorige oefening
model.vl.soort <- lm(gewicht ~ vin_c * soort, data = peng)
OpmerkingEindopdracht

Methode en data-analyse

a) Wat onderzoek je, waarop, en in welke vier gedaanten laat je de moderator meedoen? Schrijf het op zonder één uitkomst te noemen. (twee of drie zinnen; wie of wat er gemeten is, hoeveel, en de vier modellen)

b) Tel per model hoeveel parameters het mag invullen — getallen die het model zelf uit de data haalt. Welke twee paren mag je daarna tegen elkaar uitkomen, en welke niet? (vier getallen en twee paren)

c) Waarom mag je R² niet tussen modellen met ongelijke parameteraantallen vergelijken? Zeg het in gewone woorden. (één zin; zie het als draaiknoppen)

d) Kijk of elke gedaante genoeg dieren per groep overhoudt: table(peng$lang) en table(peng$klasse). En hoeveel pinguïns doen er in alle vier de modellen mee? (twee plus drie groepsaantallen, en één totaal)

Resultaten

e) Waar ligt de middelste snavellengte, en wat geeft het dummy-model? Draai summary(model.vl.lang). (rapporteer de mediaan, de helling bij de korte helft, de extra helling met p, en R²)

f) Tel die twee hellingen bij elkaar op. Welke twee getallen uit de romp van deze oefening krijg je terug? (twee getallen)

g) Nu de eerste aanwijzing. Het getal-model zei dat de helling met 1.20 per millimeter snavel meebeweegt, en de twee helften liggen gemiddeld diff(tapply(peng$snavellengte, peng$lang, mean)) millimeter uit elkaar. Hoeveel kloof vóórspelt het getal-model dus tussen de helften, en hoeveel meet de knip er? (één som en twee getallen)

h) Knip nu in drieën en toets zoals bij de drie soorten: eerst één toets over allebei de interactieregels samen, met anova(model.klasse.zonder, model.klasse.met), dan de drie hellingen los. (rapporteer de twee grenzen, F(df1, df2), p en de drie hellingen)

i) Zet de vier R²’s op een rij, en vergelijk alleen wat je bij b) mocht vergelijken. Waar zit de verrassing? (vier R²’s en één zin)

j) Die verrassing is niet gratis. Tel na wie er in die groepen zit met table(peng$lang, peng$soort) en table(peng$klasse, peng$soort). Wat blijkt de knip eigenlijk te zijn? (twee percentages en één zin)

k) Draai het nu om en geef de groep een getal. Vergelijk model.vl.gentoo met het diepte-model uit je Jouw beurt, en reken uit wat de gemiddelde snaveldiepte van de Gentoo’s bij een interactie van −8.60 zou voorspellen. (twee R²’s, twee gemiddelde diepten, en één voorspelde helling naast één gemeten)

l) Houd tot slot één soort tegelijk vast en draai dezelfde interactie binnen elke soort. Wat blijft er van de snaveldiepte over? (per soort b, SE en p, en daarna twee of drie zinnen met je eerlijke antwoord)

Z) Mag je alles hierboven geloven? Je hebt bij b) en c) al gecontroleerd of de vergelijking eerlijk was — dus nu de vraag die erachteraan hoort. Waar kwam die mediaan van 44.45 mm vandaan, en waarom is de R² van .7762 daardoor eerder te mooi dan te lelijk? (één zin over waar het mes vandaan kwam, en in welke richting dat je uitkomst duwt)

Methode en data-analyse

a) Bijvoorbeeld: “Bij 342 pinguïns onderzochten we of het effect van vinlengte op gewicht afhangt van de snavellengte. We lieten die moderator in vier gedaanten meedoen: als continu getal, geknipt bij de mediaan, geknipt in drie even grote klassen, en vervangen door de soort waartoe het dier behoort. Van elk model bepaalden we de verklaarde variantie, en we vergeleken alleen modellen met gelijke aantallen parameters.” Geen enkel getal uit de uitkomsten.

b) Het getal-model heeft er vier (intercept, twee hellingen, één interactie), de knip in tweeën ook vier, de knip in drieën zes, en het soort-model zes. Je mag dus twee paren tegen elkaar uitkomen: het getal tegen de knip in tweeën, en de knip in drieën tegen de soort. De andere combinaties niet.

c) Een parameter erbij kan R² nooit laten zakken. Zie het als draaiknoppen: een model met zes knoppen kan er altijd twee op nul zetten en ís dan het model met vier knoppen — verliezen kan het niet, hoogstens gelijkspelen. Een model met meer knoppen naast een model met minder leggen is dus doorgestoken kaart.

d) De knip in tweeën geeft 171 en 171, de knip in drieën 115, 113 en 114 — allemaal ruim genoeg. En alle 342 pinguïns doen in alle vier de modellen mee; er valt niemand weg, dus de vier R²’s gaan over precies dezelfde dieren. Dat is geen bijzaak: draaiden twee modellen op verschillende pinguïns, dan mag je hun R² sowieso niet naast elkaar leggen.

Resultaten

e) De mediane snavellengte is 44.45 mm. Bij de korte helft is de helling b = 37.17 g/mm, en de lange helft loopt daar b = 20.10 g/mm bovenop, p < .001. R² = .7762.

f) 37.17 + 20.10 = 57.27 — precies de twee losse lijnen (37.17 en 57.27) uit de mediaan-knip eerder in deze oefening. Dezelfde twee getallen, nu uit één tabel in plaats van uit twee losse modellen.

g) De twee helften liggen 9.40 mm uit elkaar, gemiddeld 39.22 tegen 48.62 mm. Het getal-model voorspelt dus een kloof van 1.20 × 9.40 = 11.28 g/mm. En de knip meet er 20.10. Bijna het dubbele — onthoud dat, het is de eerste aanwijzing.

h) De grenzen liggen op 40.80 en 46.73 mm. De toets over de twee interactieregels samen: F(2, 336) = 20.89, p < .001 — dezelfde vrijheidsgraden als bij de soort-vergelijking in de vorige oefening, want ook hier kwamen er twee regels bij. De drie hellingen: 23.02, 42.22 en 60.26 g/mm. De waaier is wéér breder geworden: het getal gaf 38.84 tot 51.94, de knip in tweeën 37.17 tot 57.27, en deze drie klassen lopen van 23.02 tot 60.26.

i) .7694 voor het getal, .7762 voor twee groepen, .7863 voor drie groepen en .7896 voor de soort. Vergelijk alleen de twee toegestane paren, en daar staat de verrassing in het bovenste: de knip wint van het getal, 77.6% tegen 76.9%, met evenveel parameters. Terwijl deze oefening je net op het hart drukte dat knippen informatie weggooit. In het onderste paar wint de soort van de geknipte klassen, 79.0% tegen 78.6% — de vermomming verliest van het origineel.

j) Van de korte helft is 87% Adélie (148 van de 171), en omgekeerd zit 98% van de Adélies (148 van de 151) in de korte helft. De kortste van de drie klassen is zelfs zuiver: 115 Adélies en niemand anders. De knip op snavellengte is dus bijna dezelfde knip als “Adélie of niet”. Hij gooit wel degelijk informatie weg; hij vindt er alleen per ongeluk iets belangrijkers voor terug — de soort, de moderator uit de vorige oefening. Dat verklaart je aanwijzing uit g): de waaier kantelt niet glad met de snavellengte mee, hij springt waar de soorten wisselen, en een knip kan zo’n sprong volgen terwijl een rechte kanteling eroverheen middelt.

k) De interactie met de groep komt op 25.13 met p < .001 en R² = .7846; die met het getal gaf −8.60 met R² = .7870. De gemiddelde snaveldiepte is 18.37 mm bij de niet-Gentoo’s en 14.98 bij de Gentoo’s, dus 3.39 mm ertussen. Het diepte-model maakt daarvan −8.60 × (14.98 − 18.37) = 29.15 g/mm. Het getal voorspelt dus dat de Gentoo’s ruim 29 g/mm steiler lopen dan de rest, en de dummy meet zelf 25.13 — het getal doet de groep bijna na. (Die 29.15 rekent met afgeronde getallen; met alle decimalen die het model zelf heeft komt er 29.12 uit. Vandaar ruim 29 en geen twee decimalen — zoveel precisie draagt deze som niet.) Schrik niet van de mintekens: de interactie van de diepte is negatief én de Gentoo-snavel ligt 3.39 mm ónder die van de rest, en min keer min is plus. De diepte spreekt de dummy niet tegen; hij vertelt bijna hetzelfde verhaal, in spiegelschrift. En de twee R²’s ontlopen elkaar amper, 78.5% tegen 78.7% — allebei vier parameters, dus deze vergelijking mág.

l) Binnen één soort valt er aan de soort niets meer te variëren, dus wat daar nog modereert kán de soort niet zijn. Adélies: b = −0.65, SE = 4.36, p = .88 — het getal is kleiner dan zijn eigen ruis. Chinstraps: b = 8.09, SE = 4.75, p = .094, ook niets om over naar huis te schrijven. Gentoo’s: b = −12.55, SE = 4.84, p = .011 — daar blijft hij staan, en hij is er zelfs stérker dan over het hele bestand (−8.60).

Het eerlijke antwoord heeft dus twee delen, en dat mag: bij twee van de drie soorten stond de snaveldiepte inderdaad alleen maar in voor de soort, en bij de Gentoo’s doet hij zelf iets. Pas op voor de valkuil aan béíde kanten: “de diepte stond alleen maar in voor de soort” verklaart de Gentoo’s niet, en “de diepte is gewoon zelf de moderator” verklaart de Adélies en de Chinstraps niet.

Z) Die mediaan van 44.45 mm kwam uit dezelfde 342 pinguïns waarop je het model daarna hebt beoordeeld. Je hebt het mes dus geslepen op het hout dat je ermee ging zagen, en dat duwt maar één kant op: de R² van .7762 is eerder te mooi dan te lelijk. Op een nieuwe groep pinguïns zou de knip het slechter doen dan hier, want de mediaan van díé groep ligt ergens anders.

Dat maakt je conclusie bij i) niet ongeldig — het verschil tussen .7762 en .7694 is echt, en de verklaring bij j) is dat ook: 148 van de 151 Adélies zitten werkelijk in de korte helft, en dat verzin je niet met een mediaan. Maar het maakt de winst kleiner dan hij eruitziet. Hoe je zoiets fatsoenlijk meet — een model beoordelen op dieren die niet meededen aan het bouwen ervan — staat niet in dit oefenboek, en dat zeg je er dan ook eerlijk bij.

En de weg terug, voor de gewone aannames: gingen de residuen trechteren, dan blijven staan alle hellingen bij e), g), h) en k), en alle vier de R²’s — dat zijn beschrijvingen van deze 342 dieren. Doe je over: de p bij e), de F bij h) en de drie p’s bij l), want die rusten op één gedeelde restspreiding.

rapportage

Schrijf het op als een onderzoeker

rapportagezin · simple slopes · b en p

OpmerkingRapportage

Bij een continue moderator moet je rapportagezin iets extra’s doen wat bij een groep niet nodig was: hij moet zeggen op welke waarden je de hellingen hebt afgelezen. Anders kan niemand jouw getallen narekenen — ze hangen immers af van een keuze die jij gemaakt hebt. Zo:

Om te onderzoeken of het effect van vinlengte op gewicht afhangt van de snavellengte, voerden we bij 342 pinguïns een regressieanalyse uit met beide voorspellers gecentreerd en hun interactieterm in het model. De interactie was significant, b = 1.20, p < .001; het model verklaarde 76.9% van de variantie in gewicht. Bij een gemiddelde snavellengte (43.92 mm) ging een millimeter extra vin samen met 45.39 gram extra gewicht. Simple slopes op één standaardafwijking onder en boven het gemiddelde lieten zien dat dit effect toenam met de snavellengte: 38.84 g/mm bij een korte snavel (38.46 mm) en 51.94 g/mm bij een lange snavel (49.38 mm).

Kijk wat er cursief staat en wat niet. Alléén de symbolen b en p hellen; de gewone woorden, de millimeters en de grammen staan rechtop. De statistieken krijgen twee decimalen, en een p onder de .001 schrijf je als p < .001, zonder nul ervoor. Eén ding wijkt daarvan af, en dat is met opzet: R² staat hier als percentage met één decimaal (76.9%) in plaats van als .77. In deze oefening zit het verschil tussen de modellen in de derde decimaal, en op twee decimalen vallen precies de twee samen die je uit elkaar moet kunnen houden — dat zie je in de eindopdracht gebeuren. Alles staat in de verleden tijd. En de zin staat op eigen benen: wie ’m uit dit blok tilt, weet wie er gemeten werd, wat er gemeten werd, dat er gecentreerd is, en op welke drie waarden de hellingen zijn afgelezen.

Dat laatste is de nieuwe eis. Zou je alleen “de simple slopes waren 38.84 en 51.94” opschrijven, dan staat er een getal zonder adres — en dan is het net zo onnavolgbaar als de −7.31 uit het ongecentreerde model.

Jouw beurt: schrijf zelf de rapportage-zin voor je jouw-beurt-taak — de snaveldiepte. Let op twee dingen: de interactie is daar negatief, dus de reeks loopt omláág, en de standaardafwijking is een andere (1.97 in plaats van 5.46), dus je leest af op andere waarden. De modelzin hierboven mag je lenen; hardop voorlezen mag ook.

En dan de eindopdracht van het hoofdstuk, want die brengt zijn eigen rapportage-plicht mee. Je hebt daar dezelfde moderator in twee gedaanten gebruikt — als getal en als knip — en allebei leveren ze een nette, eerlijke zin op, met verschillende getallen erin. Je zin moet dus zeggen in welke gedaante de moderator meedeed, anders kan niemand zien welke analyse jij gedaan hebt, laat staan hem narekenen. Zo klinkt de knip-versie:

Om te onderzoeken of het effect van vinlengte op gewicht afhangt van de snavellengte, deelden we een steekproef pinguïns (N = 342) bij de mediane snavellengte (44.45 mm) in een korte en een lange helft (n = 171 per helft) en voerden we een regressieanalyse uit met vinlengte, deze tweedeling en hun interactieterm als voorspellers. De interactie was significant, b = 20.10, SE = 3.95, 95% CI [12.34, 27.86], p < .001; het model verklaarde 77.6% van de variantie in gewicht. Bij de pinguïns met een korte snavel ging een millimeter extra vin samen met 37.17 gram extra gewicht, bij die met een lange snavel met 57.27 gram.

Leg hem naast de modelzin bovenaan deze rapportage. Die moest zeggen op welke waarden je de hellingen hebt afgelezen; deze moet zeggen waar het mes viel. Twee keer dezelfde plicht — een helling zonder adres is niet na te rekenen — alleen het soort adres verschilt. En ja: eerder in deze oefening stond dat je níét moet knippen als je echt gaat rapporteren. Dat advies staat nog steeds. Maar wie tóch knipt, zegt het er in elk geval bij, want het is een keuze die de getallen verandert — en dat is precies wat je in de eindopdracht zag gebeuren.

Jouw beurt, de laatste van dit hoofdstuk: schrijf de rapportage-zin voor de knip in drieën. Je hebt eraan: de modelvergelijking F(2, 336) = 20.89, p < .001, de drie hellingen (23.02, 42.22 en 60.26 g/mm) en de twee grenzen (40.80 en 46.73 mm). De eerste modelzin uit de rapportage van de vorige oefening — die over de drie soorten — mag je als voorbeeld ernaast leggen, want dit is precies dat model, alleen met geknipte klassen in plaats van echte soorten. En keur je eigen zin daarna op de nieuwe eis: kan een lezer aan jouw woorden zien dát de snavellengte geknipt is, en waar de grenzen lagen?

In het boek

Hier: als de moderator een getal is, kun je hem niet meer tekenen — je moet zelf kiezen waar je gaat staan om de helling af te lezen.

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

  1. 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.↩︎