Vierde beweging · Verschil onthechten
Oefening 8.1 · Meervoudige regressie
Twee continue voorspellers in één model, en wat elk uniek toevoegt · multiple regression
Data bij deze oefening — penguins.csv 342 pinguïns · diamantjes.csv de twaalf steentjes uit het boek · wat elke kolom betekent
prikkel
Trek de regressielijn in R — eerst zonder, dan met de vin
tweede voorspeller · zakkende helling · R²
d <- read.csv("diamantjes.csv")
# eerst karaat alleen
lm(glans ~ karaat, data = d) # → b = 32.14
# en nu met gladheid erbij — let op wat er met karaat gebeurt
lm(glans ~ karaat + gladheid, data = d) # → b = 0.98lm(), 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 dat klaarzetten is waar het echte werk zit. Alle handelingen die eraan voorafgaan — filteren, kolommen maken, groepen samenvatten, koppelen — staan bij elkaar in de gereedschapskist achterin. Voor dit blok in het bijzonder: centreren, want zodra je twee voorspellers naast elkaar zet, gaat het intercept iets betekenen dat je wilt kunnen voorlezen.
Los voorspelt karaat de glans met b = 32.14, maar zet gladheid ernaast en er blijft b = 0.98 van over — terwijl R² van .33 naar .82 springt. Niet het gewicht maakte de steentjes mooi; de slijper deed dat, en karaat liftte mee.
Nu jij — bij de pinguïns, die je goed kent inmiddels. De snavellengte voorspelt het gewicht keurig — tot je de vinlengte ernaast zet.
Haal de pinguïns binnen en trek eerst die ene simpele lijn: gewicht voorspeld uit snavellengte, meer niet. Kijk dan wat er gebeurt als je de vinlengte erbij zet.
peng <- read.csv("penguins.csv")
# eerst: gewicht uit alleen snavellengte
model.sn <- lm(gewicht ~ snavellengte, data = peng)
summary(model.sn)
coef(model.sn)[["snavellengte"]] # → 87.42 ruim 87 gram per millimeter snavel
summary(model.sn)$r.squared # → .354
# nu: gooi de vinlengte erbij
model.sn.vl <- lm(gewicht ~ snavellengte + vinlengte, data = peng)
summary(model.sn.vl)
coef(model.sn.vl)[["snavellengte"]] # → 6.05 de helling zakt door de bodem
summary(model.sn.vl)$r.squared # → .760Kijk naar dat ene getal bij snavellengte. In het eerste model staat er 87.42: ruim 87 gram per millimeter, en een derde van alle verschil in gewicht verklaard. Zet je de vinlengte ernaast, dan zakt diezelfde helling naar 6.05 — nog geen veertiende van wat hij was. En de R² doet precies het omgekeerde van de helling: die springt van .35 naar .76.
Dat is dezelfde beweging als bij de steentjes, alleen harder. Daar viel karaat van 32.14 naar 0.98 zodra gladheid binnenkwam; hier valt de snavel van 87.42 naar 6.05 zodra de vin binnenkomt. De vraag is dus dezelfde: wat hoorde er nou eigenlijk bij de snavel?
Lees de uitvoer van onderen naar boven
F-toets · adjusted R² · standaardfout
Die vraag houden we nog even vast, want je hebt zonet meer uitgerekend dan je gelezen hebt. summary() drukte een hele lap af, en tot nu toe viste je er twee getallen uit: de helling en de R². In diezelfde lap staat óók alles wat een artikel over dit model zou rapporteren — je moet alleen weten in welke volgorde je hem leest. En die volgorde is niet van boven naar beneden.
Eigenlijk wil je eerst het model als geheel toetsen: is er iets, of is er niets aan de hand? De nulhypothese is hier lekker groots: dit model — snavellengte en vinlengte samen — verklaart níéts van de verschillen in gewicht. De alternatieve zegt dat hij er wél iets van verklaart. En dát toetsen we met de F-toets: één getal dat afweegt hoeveel verschil in gewicht het model verklaart tegenover hoeveel er als ruis overblijft.
In symbolen ziet diezelfde nulhypothese er zo uit:
\[H_0:\ \rho^2 = 0 \qquad H_1:\ \rho^2 > 0\]
of, wat hetzelfde zegt maar dan via de twee hellingen:
\[H_0:\ \beta_1 = \beta_2 = 0 \qquad H_1:\ \text{minstens één } \beta_j \neq 0\]
Dat die twee hetzelfde zeggen, kun je zien: zijn béide hellingen nul, dan voorspelt het model overal hetzelfde getal en verklaart het dus niets — en dan is de verklaarde variantie in de populatie (ρ²) ook nul. Andersom net zo. Die \(j\) is niet meer dan een aftelnummer: hier staat hij voor 1 of 2, en in een model met vijf voorspellers voor 1 tot en met 5.
Eén van de twee is genoeg, en dat is de vraag die je je nu stelt. Het zijn niet twee hypothesen die je allebei moet kennen maar één hypothese in twee notaties. Welke je opschrijft hangt af van waar je zin daarna heen gaat: praat je over hoevéél het model van het gewicht verklaart, neem dan de ρ²-vorm; ga je daarna per voorspeller door, neem dan de β-vorm, want die sluit aan op de regels van de tabel die je zo gaat lezen.
Let op de letters. Er staan Griekse, want de hypothese gaat over de populatie: die ρ², en de hellingen dáár (β). Je uitvoer geeft R² en b — Latijnse letters, de schattingen uit déze 342 pinguïns. Een nulhypothese over R² zou onzin zijn: die heb je al, je kunt hem gewoon aflezen.
En \(H_1\) kijkt hier maar één kant op — groter dan nul, niet ongelijk aan nul — want een model kan niet minder dan niets verklaren.
# R zet dit onderaan zijn uitvoer, en je begint er juist
summary(model.sn.vl)$fstatistic # → 536.63 2 339
summary(model.sn.vl)$r.squared # → .7600
summary(model.sn.vl)$adj.r.squared # → .7585
# en de p die bij die F hoort
f <- summary(model.sn.vl)$fstatistic
pf(f[1], f[2], f[3], lower.tail = FALSE) # ruim kleiner dan .001R zet die F op de állerlaatste regel van de uitvoer, maar je begint er dus: F(2, 339) = 536.63, p < .001. De 2 telt de hellingen in het model — één voor de snavel, één voor de vin — en de 339 is de rest: 342 pinguïns min de drie getallen die het model schat, die twee hellingen plus het intercept. Die p ken je uit de oefening over de p-waarde: in werelden waarin het model werkelijk niets verklaart, maakt toeval zó’n F vrijwel nooit na. De nulhypothese dat dit model niets van de verschillen in gewicht verklaart, mag dus de deur uit. En pas als dát vaststaat, heeft de rest van de tabel zin — is er niks, dan valt er ook niks te verdelen.
Let op wat die F wél en niet zegt, want dat wordt zo belangrijk. Hij toetst het model als geheel. Hij zegt dus: hier zit iets in. Hij zegt niet dat allebei de voorspellers dat iets leveren.
Vlak boven de F staat een regel die we nog niet hadden: Adjusted R-squared. De gewone R² is uitgerekend op déze 342 pinguïns, en hij vleit zichzelf: elke voorspeller die je toevoegt kan hem alleen maar omhoog duwen, nooit omlaag — ook een kolom met puur toeval erin duwt hem nog een fractie op. De adjusted R² haalt die vleierij eraf: hij generaliseert de R² van jouw steekproef naar andere mogelijke steekproeven — naar de populatie dus. Hier schelen de twee elkaar bijna niets, .760 tegenover .759, en dát is het nieuws: met 342 pinguïns en maar twee hellingen valt er nauwelijks iets af te halen. Lopen die twee ver uiteen, dan heb je te veel voorspellers op te weinig data gezet, en schept je R² vooral op over jouw eigen steekproef.
Pas als het geheel zegt dat er iets is, ga je uitzoeken waar het aan ligt: aan de snavel, aan de vin, of aan allebei. Daarvoor is het blok Coefficients middenin de uitvoer — vier kolommen, en je hebt ze stuk voor stuk al eens ergens gezien:
tabel <- coef(summary(model.sn.vl))
tabel # de vier kolommen, voor elke voorspeller
tabel["vinlengte", 1:4] # → 48.14 2.01 23.94 0.0000
tabel["snavellengte", 1:4] # → 6.05 5.18 1.17 0.2438
confint(model.sn.vl)["vinlengte", ] # → 44.19 52.10
confint(model.sn.vl)["snavellengte", ] # → -4.14 16.24- Estimate — de hellingen zelf. De 48.14 van de vin en de 6.05 van de snavel. (De bovenste regel, het intercept, voorspelt het gewicht bij een vin én een snavel van nul millimeter — laat die even voor wat hij is; die pinguïn bestaat niet.)
- Std. Error — de standaardfout van elke helling: hoe ver zo’n helling van steekproef tot steekproef zou dansen. Dezelfde wiebel die je in de oefening over de steekproevenverdeling voor het gemiddelde uitrekende, nu voor een helling. Voor de vin 2.01, voor de snavel 5.18.
- t value — de helling gedeeld door zijn standaardfout: hoeveel standaardfouten ligt hij van nul af? Voor de vin 23.94, voor de snavel 1.17. (Deel je de gedrukte 6.05 door de gedrukte 5.18, dan krijg je 1.168 — R rekent met álle decimalen, niet met de afrondingen.)
- Pr(>|t|) — de p die daarbij hoort.
Elke regel toetst zo zijn eigen kleine nulhypothese: déze helling is in de populatie nul, met de andere voorspeller vastgehouden. En hier lopen de twee regels ver uiteen. De vin ligt bijna vierentwintig standaardfouten van nul, p < .001: die levert, náást de snavel, een significante unieke bijdrage aan het voorspellen van het gewicht. De snavel ligt er nog geen anderhalve standaardfout vandaan, p = .24. Voor de snavel houd je de nulhypothese dus gewoon overeind: naast de vin is er geen bewijs dat hij nog íéts van zichzelf bijdraagt aan het voorspellen van het gewicht.
Kijk daar even naar, want het is precies wat de prikkel vroeg. Dezelfde kolom die in zijn eentje ruim een derde van het gewichtsverschil verklaarde, haalt naast de vin de drempel niet eens meer.
En confint() geeft het laatste stuk: het 95%-betrouwbaarheidsinterval van elke helling. Voor de vin loopt het van 44.19 tot 52.10. Lees het zoals je élk interval leest: eerst waar het loopt — alle waarden voor de echte helling die nog goed bij deze data passen, van krap 44 tot ruim 52 gram per millimeter — en dan of nul meedoet. Nul valt er ruim buiten.
Voor de snavel loopt het van −4.14 tot 16.24, en daar zit nul dus gewoon tússen. Dat is dezelfde uitspraak als die p van .244, in andere kleren: het interval sluit een helling van nul niet uit. Sterker nog, het sluit ook een negatieve helling niet uit. Het interval is het schat-antwoord, de p de toets-bijvraag die er al in zit; het is hetzelfde feit, twee keer bekeken. Valt nul buiten het 95%-interval, dan is de p kleiner dan .05 — en valt nul erbinnen, zoals hier, dan is hij groter.
Dat interval zegt trouwens iets wat de p níét zegt, en het is het waard er even bij stil te staan. p = .24 betekent niet “de snavel doet niets”. Het interval laat zien wat er nog past bij deze data: van vier gram eraf tot zestien gram erbij. Wat je hebt vastgesteld is niet dat de bijdrage nul ís, maar dat je hem met 342 pinguïns niet van nul kunt onderscheiden.
Daarmee heb je de uitvoer gelezen zoals een onderzoeker hem leest: eerst het geheel, dan pas de delen. Studeer je psychologie, dan loop je deze trap in het tweede jaar gewoon opnieuw, in het MVDA-werkboek van Peter de Heus2 — alleen in een andere volgorde: daar keur je éérst de residuen en ga je pas daarna het model in, hier kwam het model eerst en deed je die keuring al eerder, in de oefening over de residuen. Geen van beide volgordes is de juiste; je moet ze allebei gedaan hebben voordat je je conclusie opschrijft. Maar twee dingen staan nog open: of je die p van de vin zomaar op zijn woord moet geloven — en hoevéél elk van de twee nou eigenlijk zelf inbrengt. Dat tweede eerst, want daar is een getal voor.
Schud de vinlengte los van het gewicht
schud-toets · husselen · helling uit puur toeval
Voordat we dat uitzoeken, eerst de eerlijkheidscheck: is die 48.14 van de vin écht, of kan zo’n helling ook uit het niets ontstaan? We spelen vals. We houden de snavellengte netjes op z’n plek, maar we knippen elke vinlengte los van z’n eigen pinguïn en plakken ’m op een willekeurige andere.
Dat loskoppelen is de hele truc, en het is het waard er even bij stil te staan: door de vinlengtes te husselen haal je met eigen handen weg dat een bepaalde vin bij een bepaald gewicht hoort. In zo’n losgekoppelde wereld kán de vin niks meer bijdragen — elk beetje helling dat je er dan nóg in meet, is per definitie puur toeval. Doe het duizend keer en kijk waar de vin-helling belandt.
# één keer: vinlengte door elkaar husselen, snavellengte houden we vast
geschud <- transform(peng, vinlengte = sample(vinlengte, replace = FALSE))
coef(lm(gewicht ~ snavellengte + vinlengte, data = geschud))["vinlengte"] # ~ rond 0
# doe het 1000 keer en kijk waar de vin-helling dan landt
set.seed(42)
nep <- replicate(1000, {
geschud <- transform(peng, vinlengte = sample(vinlengte, replace = FALSE))
coef(lm(gewicht ~ snavellengte + vinlengte, data = geschud))[["vinlengte"]]
})
hist(nep, xlim = c(-10, 50)) # een bergje rond nul
abline(v = 48.14, col = "red", lwd = 2) # en dáár staat de echte
sd(nep) # → 2.52 zo breed wiebelt het bergje
range(nep) # → -7.77 7.74
sum(nep >= 48.14) # → 0 van de duizendKijk naar dat plaatje. Duizend werelden die dansen rond nul — want in elke geschudde wereld is de band tussen vin en gewicht doorgeknipt, maar toeval laat de helling nooit precies op nul landen. Ze wiebelen allemaal dicht om nul, van zo’n −8 tot +8 gram per millimeter, niet verder. En dan die rode streep op 48.14, mijlenver buiten het hele bergje. In duizend keer schudden haalde geen enkele wereld het. Dát is wat “geen toeval” betekent — geen p-waarde die je moet geloven, maar een afstand die je zíet. De 48.14 is echt van de vin.
Nog één ding, omdat je de tabel al gelezen hebt. Het schudden en de p stelden dezelfde vraag — kan toeval dit namaken? — en gaven hetzelfde antwoord: nee. Maar ze meten met een andere lat, en dat zie je meteen: het geschudde bergje wiebelt 2.52 om nul (sd(nep)), terwijl de standaardfout van diezelfde helling 2.01 was — dat is de SE van de vin uit de coëfficiëntentabel, een paar secties terug. Dat is geen slordigheid, en het is de moeite waard om te weten waaróm ze verschillen: in een wereld waarin de vin is losgeknipt, is er niemand meer die het gewicht goed voorspelt — alleen de snavel blijft over, en die kan het niet alleen. Er blijft dus veel meer onverklaarde ruis liggen, en in meer ruis wiebelt elke helling breder. Onthoud dus het antwoord van het schudden, en de lat van de tabel.
Dus de vin draagt echt bij, óók naast de snavel. Blijft de andere kant nog over: hoevéél brengt ieder van de twee nou eigenlijk zelf mee?
snap het
Wat voegt elk van de twee nog toe?
modelvergelijking · ΔR² · hiërarchische regressie
Je hebt nu twee losse antwoorden. Het schudden zei: die 48.14 van de vin is echt. De tabel zei: de snavel is naast de vin niet van nul te onderscheiden. Maar allebei die antwoorden zijn ja-of-nee, en de vraag die overblijft is hoeveel. Als de vin zó dominant is, verklaart de snavel er dan écht niets meer bovenop? Dat kun je scherp krijgen door de twee modellen zélf tegen elkaar te laten toetsen. Zo’n toets heet een modelvergelijking, en anova() doet hem voor je. Schrik niet van die naam: je kent anova() van het vergelijken van groepsgemiddelden, maar onder de motorkap is het dezelfde F-afweging — verklaarde ruimte tegen overgebleven ruis. Geef je hem twee modellen in plaats van één, dan weegt hij die tegen elkaar af.
Wat je hier aan het doen bent, heeft trouwens een eigen naam. Telkens het model een stap groter maken en per stap kijken hoeveel er extra verklaard wordt van de verschillen in gewicht: dat heet een hiërarchische regressie (hierarchical regression). Meer is het niet; het woord klinkt zwaarder dan het werk.
En hier kun je hem twee kanten op draaien, want er zijn twee volgordes:
# wat voegt de VIN toe, bovenop de snavel?
anova(model.sn, model.sn.vl) # → F = 573.09 op 1 en 339 p < .001
summary(model.sn.vl)$r.squared - summary(model.sn)$r.squared # → .406
# en nu andersom: wat voegt de SNAVEL toe, bovenop de vin?
model.vl <- lm(gewicht ~ vinlengte, data = peng)
anova(model.vl, model.sn.vl) # → F = 1.36 op 1 en 339 p = .244
summary(model.sn.vl)$r.squared - summary(model.vl)$r.squared # → .001Lees die twee naast elkaar, want dat is het hele blok in vier getallen.
De vin voegt, bovenop de snavel, F(1, 339) = 573.09, p < .001 toe, en tilt de R² met .406 omhoog. De snavel voegt, bovenop de vin, F(1, 339) = 1.36, p = .24 toe: één duizendste R².
Die twee getallen hebben een naam: de unieke bijdrage van elke voorspeller aan de R². Onthoud ze, want je komt ze in de volgende oefening langs een heel andere weg opnieuw tegen — en dan blijkt er nog een derde getal te bestaan dat van geen van beide is. Merk ook op dat die p van .244 exact de p uit de coëfficiëntentabel is. Voor één voorspeller erbij zijn de modelvergelijking en de t-toets van die helling hetzelfde ding.
Nog één regel uit de uitvoer, want die is verraderlijker dan hij eruitziet:
summary(model.vl)$adj.r.squared # → .75828 alleen de vin
summary(model.sn.vl)$adj.r.squared # → .75854 snavel erbij
# en nu een kolom pure ruis in plaats van de snavel
set.seed(7) # zodat jij dezelfde ruis trekt als wij
peng$ruis <- rnorm(nrow(peng))
summary(lm(gewicht ~ vinlengte + ruis, data = peng))$adj.r.squared # → .75762De gewone R² moest omhoog — die kan niet anders. De adjusted R² hoefde niet mee, en ging tóch een piepklein beetje omhoog: van .75828 naar .75854, twee tienduizendsten. Zet er een kolom pure ruis in plaats van de snavel in, dan zakt hij meestal wél: onze ruiskolom bracht hem van .75828 naar .75762.
Meestal, en dat woord staat er met opzet. Haal die set.seed(7) weg, draai het opnieuw, en in ongeveer één op de drie keer stijgt de adjusted R² juist — van een kolom die per ongeluk niets weet. Hij is dan ook geen slot op de deur maar een voorzichtiger gok naar wat je model in de populatie waard zou zijn, en zo’n gok gaat soms mis. Een voorspeller die de adjusted R² niet laat zakken, heeft daarmee dus nog niets bewezen — zijn lat ligt véél lager dan die van de significantie. De snavel haalde die lat, en is alsnog niet van nul te onderscheiden.
Voor wie verder wil — wanneer stijgt hij dan wél?
De oefening is compleet zonder dit kader; sla het gerust over.
Dat “ongeveer één op de drie” is geen vuistregel maar volgt uit een scherpe grens, en die grens heb je hierboven al zien werken. De adjusted R² gaat omhoog precies wanneer de F van de nieuwe voorspeller boven de 1 uitkomt, en omlaag zodra hij eronder blijft. Kijk maar terug: de snavel bracht F = 1.36 mee, net boven de 1 — en inderdaad kroop de adjusted R² van .75828 naar .75854. De ruiskolom van set.seed(7) bracht F = 0.07 mee, ver eronder, en daar zakte hij. Bij één voorspeller erbij komt dat op hetzelfde neer als t onder de −1 of boven de 1.
Daarmee is die “ongeveer één op de drie” ook precies te maken, zonder ook maar één keer te schudden. Een kolom pure ruis levert een t op die rond de nul zwabbert, en zo’n t komt in .318 van de trekkingen voorbij de −1 of de 1. Dat getal komt uit de t-verdeling; je hoeft het niet zelf te kunnen halen, wel te weten dat het er ligt.
De les van deze oefening, in drie regels:
- De schud-toets en de toetsen in de tabel zeiden, elk met hun eigen meetlat: die 48.14 is echt, de vin draagt bij aan het gewicht — óók naast de snavel.
- Het tweede model zei: van de “87 gram per millimeter snavel” bleef, met de vin ernaast, 6.05 over — en dat is niet van nul te onderscheiden.
- De modelvergelijking zei hoe groot dat “niets” is: ΔR² = .001, oftewel ongeveer wat een willekeurig getal ook had opgeleverd.
Echt en hoe groot zijn dus twee verschillende vragen, en je hebt ze allebei beantwoord. Maar er is een derde die nog open ligt, en het is de interessantste: waar is die 87 gram gebléven? Een helling die instort verdwijnt niet in het niets — hij gaat ergens heen. Dat is de volgende oefening.
jouw beurt
Doe het hele rondje op de twaalf steentjes
modeluitvoer lezen · modelvergelijking · ΔR²
Aan het begin van dit blok kwamen de diamantjes in twee regels langs: karaat voorspelde de glans met b = 32.14, en met gladheid ernaast bleef er b = 0.98 van over. Daar keek je naar. Nu loop je de hele route zelf, op diezelfde twaalf steentjes — dezelfde volgorde als hierboven: eerst het model als geheel, dan per voorspeller, dan wat elk er nog bovenop legt.
Twaalf steentjes is weinig, en dat is precies waarom ze hier staan. Een paar schermen terug beloofde dit blok je iets over de adjusted R²: lopen die twee ver uiteen, dan heb je te veel voorspellers op te weinig data gezet. Bij de pinguïns liepen ze niet uiteen. Hier wel.
Draai eerst de drie modellen. Tussen haakjes staat telkens wat er in een compleet antwoord hoort, zodat je jezelf kunt nakijken zonder ons antwoord te openen. De laatste heet met opzet Z) en niet l): de aannames staan altijd achteraan.
d <- read.csv("diamantjes.csv")
model.ka <- lm(glans ~ karaat, data = d)
model.gl <- lm(glans ~ gladheid, data = d)
model.ka.gl <- lm(glans ~ karaat + gladheid, data = d)a) Lees model.ka en model.gl los af. (per model b mét eenheid, R² en p)
b) Toets model.ka.gl als geheel. Schrijf de nulhypothese eerst voluit op, in woorden én in symbolen. (R², F(df1, df2) en p, plus H0 en H1 — één notatie is genoeg, en welke mag je zelf kiezen)
c) Die F draagt twee vrijheidsgraden: een 2 en een 9. Waar komen ze vandaan? (reken ze allebei na uit het aantal steentjes)
d) Zet de gewone R² en de adjusted R² van model.ka.gl naast elkaar, en leg ze naast de .760 en .759 van de pinguïns hierboven. (vier getallen op drie decimalen, en één zin waarom het verschil hier zoveel groter is)
e) Lees de coëfficiëntentabel van model.ka.gl. (per voorspeller b, SE, t en p)
f) Vraag de betrouwbaarheidsintervallen op met confint(model.ka.gl). Bij welke van de twee voorspellers zit nul ertussen, en wat stel je daarmee vast? (twee intervallen, en één zin)
g) Wat is er met de helling van karaat gebeurd nu gladheid ernaast staat? Leg het naast wat de snavellengte deed toen de vinlengte ernaast kwam. (twee keer twee hellingen, en één zin over de gelijkenis)
h) Hoeveel voegt de gladheid uniek toe aan de verklaarde glans, bovenop karaat? Gebruik anova(model.ka, model.ka.gl). (ΔR² en de F-toets)
i) En andersom: hoeveel voegt karaat toe bovenop gladheid? Vergelijk die p met die uit je coëfficiëntentabel bij e). (ΔR², de F-toets, en één zin)
j) Reken de overlap tussen de twee voorspellers uit met cor(d$karaat, d$gladheid), en leg hem naast de .656 van vinlengte en snavellengte. (één getal, en één zin over wat overlap doet met wat een voorspeller nog uniek kan bijdragen)
k) Geef een inhoudelijke conclusie in niet meer dan drie zinnen. (over steentjes en slijpers, niet over getallen)
Z) Mag je alles hierboven geloven? Loop de aannames na. Je hebt er drie regels voor nodig en je kent ze alle drie al: complete.cases() uit de oefening over de residuen, resid() en fitted() van diezelfde plek, en cor() van nog eerder.
- heeft elk steentje een karaat, een glans én een gladheid, of vallen er rijen weg?
- liggen de residuen ongeveer symmetrisch, en is hun spreiding gelijkmatig over het bereik van de voorspelde glans? (knip bij de mediaan van de voorspelde waarden, en deel de brede helft door de smalle)
- zitten de twee voorspellers elkaar in de weg?
En dan de vraag die er altijd achteraan hoort: gaat er iets mis, welke antwoorden hierboven doe je dan over en welke mogen blijven staan?
a) Karaat: b = 32.14 glanspunten per karaat, R² = .33, p = .049 — nét onder de .05. Gladheid: b = 8.03 glanspunten per gladheidspunt, R² = .82, p < .001.
Onthoud die p van .049 even. Los is de karaat van deze steentjes dus significant; bij e) is daar niets van over.
b) R² = .82, F(2, 9) = 20.32, p < .001. De nulhypothese die je verwerpt is dat dit model — karaat en gladheid samen — niets verklaart van de verschillen in glans tussen de steentjes. In symbolen schrijf je hem zoals bij de pinguïns: H0: ρ² = 0 tegen H1: ρ² > 0, met Griekse letters omdat de hypothese over de populatie gaat. Het aantal hellingen is hetzelfde als daar, dus de β-vorm ziet er ook hetzelfde uit; alleen de steekproef eronder is een andere.
c) De 2 telt de hellingen in het model: één voor karaat, één voor gladheid. De 9 is wat er van de twaalf steentjes overblijft nadat het model zijn drie getallen geschat heeft — die twee hellingen plus het intercept. Dus 12 − 3 = 9.
d) De steentjes: R² = .819 tegen adjusted R² = .778. De pinguïns: R² = .760 tegen .759. Bij de pinguïns kost een tweede helling bijna niets, want er zijn 342 dieren om hem op te schatten; bij twaalf steentjes eet diezelfde helling een merkbaar deel van je vrijheidsgraden op, en dát is wat de adjusted R² eraf haalt.
e) Karaat: b = 0.98, SE = 10.14, t(9) = 0.10, p = .93. Gladheid: b = 7.93, SE = 1.62, t(9) = 4.91, p < .001.
coef(summary(model.ka.gl)) # → karaat 0.98 10.14 0.10 0.9251
# → gladheid 7.93 1.62 4.91 0.0008Kijk wat er met karaat gebeurd is. Los haalde hij p = .049; naast de gladheid staat hij op .93. Dat is niet “iets minder significant” — dat is de drempel niet eens meer in zicht.
f) Karaat: 95% CI [−21.95, 23.91] — nul zit er ruim tussen, en een flink negatieve helling is niet eens uitgesloten. Gladheid: 95% CI [4.27, 11.59] — nul valt er ruim buiten. Voor karaat heb je dus níét vastgesteld dat hij niets bijdraagt aan de glans; je hebt vastgesteld dat je het met twaalf steentjes niet van nul kunt onderscheiden.
g) Karaat zakt van 32.14 naar 0.98; de snavellengte zakte van 87.42 naar 6.05. Twee keer dezelfde beweging: een voorspeller die er in zijn eentje toe lijkt te doen, houdt naast een sterkere buurman vrijwel niets van zichzelf over.
h) ΔR² = .485, F(1, 9) = 24.06, p < .001. De gladheid brengt in haar eentje bijna de halve verklaarde glans mee.
i) ΔR² = .0002, F(1, 9) = 0.01, p = .93. Diezelfde p van .93 stond al in je coëfficiëntentabel bij e), en dat hoort zo: voor één voorspeller erbij zijn de modelvergelijking en de t-toets van die helling hetzelfde ding.
j) r(karaat, gladheid) = .627, tegen .656 voor vinlengte en snavellengte. Vrijwel evenveel overlap, en dezelfde afloop: wat twee voorspellers samen weten, kan geen van beide nog als eigen inbreng opeisen.
k) Bijvoorbeeld: bij deze twaalf steentjes gaat een beter geslepen steen samen met meer glans, en die samenhang houdt stand als je het karaat ernaast legt. Het karaat zelf voegt aan de glans niets toe wat van nul te onderscheiden is. Wie een steen op zijn karaat uitkiest, betaalt dus voor de slijper zonder het te weten — en ziet dat pas als de slijpkwaliteit mee het model in gaat.
Let op dat dit een uitspraak over déze twaalf steentjes is, en geen uitspraak dat karaat er niet toe doet: bij f) zag je hoe breed het interval nog openstaat.
Z) Hier past geen kale “ja”. Twee van de drie checks geven een hard antwoord; de derde is precies degene waar twaalf steentjes te weinig voor zijn, en dát is het antwoord.
Draai de drie checks. resid() en fitted() ken je uit de oefening over de residuen; de rest is optellen en delen.
sum(complete.cases(d[, c("karaat", "glans", "gladheid")])) # → 12 van de 12
e <- resid(model.ka.gl)
f <- fitted(model.ka.gl)
(mean(e) - median(e)) / sd(e) # → -0.30 scheefheid
laag <- f <= median(f) # de laagst VOORSPELDE helft
sd(e[!laag]) / sd(e[laag]) # → 1.44 hoog gedeeld door laag
cor(d$karaat, d$gladheid) # → .63
sqrt(1 / (1 - cor(d$karaat, d$gladheid)^2)) # → 1.28De knip gaat over de voorspelling, niet over de spreiding. Dat de hoog voorspelde helft hier ook de bréédste is, is de uitkomst en niet de definitie — bij een trechter die andersom loopt komt er een getal onder de 1 uit.
- Geen rij valt weg: alle 12 steentjes hebben een karaat, een glans én een gladheid, dus de drie modellen draaien op dezelfde steentjes. Dat is hier extra belangrijk, want twee modellen op verschillende steekproeven mag je niet van elkaar aftrekken. Dit is een telling, dus hier valt niets te twijfelen.
- De residuen liggen iets scheef — gemiddelde min mediaan is −0.30 standaardafwijking, dus het gemiddelde ligt ónder de mediaan — en de hoogst voorspelde helft spreidt 1.44 keer zo ver als de laagste. Allebei die getallen zeggen bij twaalf steentjes bijna niets: de gewone wiebel is daar zó breed dat er nauwelijks iets buiten valt. En de knip zelf wiebelt mee, want met twaalf punten hangt het van één steentje af aan welke kant van de mediaan het valt. Dat maakt deze steentjes niet verdacht; het maakt je oordeel erover ongegrond — en dat is iets anders.
- De twee voorspellers zitten elkaar niet dodelijk in de weg: r = .63. Let op welke maat je daarbij pakt, want er zijn er twee en ze schelen een wortel. De variantie van de hellingen wordt 1/(1 − r²) = 1.65 keer zo groot, maar de standaardfout is een standaardafwijking en gaat dus met de wortel daarvan: 1.28 keer zo breed. Merkbaar, niet dodelijk. (Dat getal heet de variance inflation factor, VIF; je hoeft het niet uit je hoofd te kennen, wel te weten dat het over de variantie gaat en niet over de fout zelf. En vul je de gedrukte .63 in, dan krijg je 1.66 en 1.29 — reken met de volle
cor()zoals in het codeblok, niet met het afgeronde getal op deze bladzij.) Ook dit is een som op de twaalf zelf, dus ook hier valt niets te twijfelen — anders dan bij de residuen hierboven, en dat verschil is de hele les van deze vraag.
Wat dat betekent voor je antwoorden. Niet “overdoen” — er is niets aangetoond dat je zou moeten repareren. Wel: het rijtje valt in tweeën, en de tweede helft draag je met een slag om de arm mee.
- Staat vast, wat er ook gebeurt: de hellingen bij a), e) en g), de gewone R²’s, de twee ΔR²’s bij h) en i), en de overlap bij j). Dat zijn beschrijvingen van déze twaalf steentjes; daar verandert geen aanname iets aan.
- Leunt op de aanname: alles wat de stap naar buiten maakt — de F en de p bij b), e), h) en i), de standaardfouten en de t’s bij e), de twee intervallen bij f), en daarmee ook de conclusie bij k). Die komen uit een verdeling die aanneemt dat de residuen zich netjes gedragen, en dat is precies wat je met twaalf steentjes niet kunt nakijken. Ook de adjusted R² hoort in dit rijtje en niet in het vorige: die doet een gok over de populatie, en zo’n gok is niet beter dan de aannames eronder.
Zie waar de scheidslijn ligt, want die is dezelfde als bij de letters een paar schermen terug: wat over déze twaalf gaat, staat vast; wat over de populatie gaat, staat op een aanname. Bij 342 pinguïns kun je die aanname nog nakijken. Bij twaalf steentjes moet je hem lenen.
rapportage
Schrijf het op als een onderzoeker
rapportagezin · F, b en R² · twee decimalen
Bij een meervoudige regressie is het boodschappenlijstje voor het eerst echt volgroeid: je rapporteert het model als geheel, dan wat de tweede voorspeller toevoegde, en dan de voorspellers apart. De eerste en de derde laag las je in de uitvoer; de middelste kwam uit de modelvergelijking. Het klinkt als veel, maar het zijn gewone zinnen die elk één laag dragen. Zo:
Om te onderzoeken of de snavellengte iets eigens bijdraagt aan het gewicht, voorspelden we bij 342 pinguïns het gewicht uit hun snavellengte en hun vinlengte. Het model als geheel verklaarde het grootste deel van de verschillen in gewicht, R² = .76, F(2, 339) = 536.63, p < .001. Bij gelijke snavellengte ging elke millimeter langere vin samen met ruim achtenveertig gram meer gewicht, b = 48.14, SE = 2.01, t(339) = 23.94, p < .001, 95% CI [44.19, 52.10]. De snavellengte droeg daarnaast niet significant bij aan het voorspellen van het gewicht, b = 6.05, SE = 5.18, t(339) = 1.17, p = .24, 95% CI [−4.14, 16.24]; bovenop de vinlengte voegde zij .001 aan de verklaarde variantie toe, F(1, 339) = 1.36, p = .24.
Kijk wat er cursief staat en wat niet. De symbolen die zelf een waarde dragen leunen schuin — R², F, p, b, SE, t — en de gewone woorden staan rechtop. De letters CI staan ook rechtop, met een spatie vóór de haak: ze benoemen wat er tussen de haken staat en dragen zelf geen waarde. De getallen staan in de zin zelf op twee decimalen. De p is de uitzondering: die krijgt er drie zolang hij ónder de .10 ligt, en twee zodra hij erboven komt — vandaar p = .24 hierboven. In een uitleg of een antwoord mag een getal best een derde decimaal dragen als het verschil er anders af valt (een ΔR² van .001 bestaat niet op twee decimalen); in de rapportagezin houd je het strak. En een p die onder de .001 duikt schrijf je niet uit: niet p = .0004, maar p < .001, zonder nul voor de punt. Alles staat in de verleden tijd, want je vertelt wat er gevónden werd, en de zinnen staan op eigen benen: wie ze uit dit blok tilt, weet nog steeds wie er gewogen werd, wat er voorspeld werd, en met welke medespeler vastgehouden.
En let op hoe die tweede voorspeller wordt opgeschreven, want daar gaat het vaak mis. Er staat p = .24 — het getal zelf. Niet “n.s.”, niet “geen effect”, en al helemaal niet weggelaten. Een niet-significante uitkomst is een uitkomst: hij hoort er met zijn eigen getallen bij te staan, mét het interval, zodat je lezer zelf kan zien hoeveel er nog paste. Wie de snavel uit de zin laat omdat hij “niks deed”, vertelt niet wat er gevonden is maar wat er goed uitkwam.
Eén getal uit de uitvoer laat je hier bewust thuis: de adjusted R². Niet omdat hij niets waard is — hij liet je net zien hoe laag zijn lat eigenlijk ligt — maar kijken en rapporteren zijn twee verschillende klussen. Het MVDA-werkboek van De Heus geeft een advies dat we overnemen: rapporteer alleen de R², en laat de adjusted weg — behalve wanneer je het model echt als voorspelformule op nieuwe gevallen gaat loslaten (zíjn voorbeeld: een selectieformule voor nieuwe piloten). Zolang jouw model geen volgende lichting pinguïns hoeft te voorspellen, is de adjusted R² een instrument om mee te kíjken, geen getal voor in de zin.
Jouw beurt: schrijf dezelfde drie lagen op voor het steentjes-model uit de taak hierboven. Het model als geheel (R² = .82, F(2, 9) = 20.32, p < .001), dan wat de gladheid bovenop karaat toevoegde (ΔR² = .48, F(1, 9) = 24.06, p < .001), en dan de twee voorspellers apart (gladheid b = 7.93, SE = 1.62, t(9) = 4.91, p < .001, 95% CI [4.27, 11.59]; karaat b = 0.98, SE = 10.14, t(9) = 0.10, p = .93, 95% CI [−21.95, 23.91]). Noem de twaalf steentjes, wat je voorspelde en waaruit, zodat je zin uit deze opdracht te tillen is.
Let bij het opschrijven op twee dingen. Ten eerste: dat interval van karaat loopt van ruim tweeëntwintig glanspunten eraf tot bijna vierentwintig erbij, en die breedte is het eerlijke antwoord — je hebt niet vastgesteld dat karaat niets doet met de glans, je hebt vastgesteld dat je het met twaalf steentjes niet kunt uitsluiten. En ten tweede: die R² van .82 ziet er prachtig uit en is vrijwel helemaal van de gladheid. Een model als geheel kan glansrijk slagen terwijl een van zijn twee voorspellers niets van zichzelf inbrengt. Hardop voorlezen mag.
◠ Zelfde vorm in het boek
Hier: twee voorspellers naast elkaar in één model, wat elk van de twee er nog bovenop legt, en de lat waaraan je dat afmeet — wat een willekeurig getal ook al had opgeleverd.
In het boek: W8 — Meer dan één verschil
Daar draait het niet om méér knoppen, maar om dezelfde vraag als bij die correlatie-wolk: wat blijft er van een verband over als je vasthoudt wie er nog meer meepraat?
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.↩︎Peter de Heus, MVDA Exercise book 2025-2026, Universiteit Leiden — het oefenboek waarmee tweedejaars psychologie daar multivariate data-analyse (MVDA) leren. Daar komt ook het rapporteer-advies over de adjusted R² in dit blok vandaan.↩︎