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

Zevenentachtig gram per millimeter snavel

regressielijn · helling · R²

Op de rots bij Antarctica meet iemand snavels.¹ Millimeter voor millimeter, pinguïn voor pinguïn, driehonderdtweeënveertig koppen die liever waren gaan vissen.

En dat meten loont, zo lijkt het. In deze pinguïns voorspelt snavellengte het gewicht heel behoorlijk: elke millimeter snavel gaat samen met zo’n 87 gram meer pinguïn (b = 87.42), en die ene maat verklaart ruim een derde van alle verschil in gewicht (R² = .354). Wie de snavel kent, weet al aardig wat de weegschaal gaat zeggen.

Maar er is méér gemeten op die rots — de vin bijvoorbeeld. En dan wordt de vraag ineens scherper. Niet: voorspelt de snavel het gewicht? Dat doet hij, zie de getallen hierboven. Maar: vertelt de snavel ons iets wat we nog niet wisten? Iets van zichzélf? Trek even mee — we zetten de vin ernaast, en dan mag de snavel laten zien wat hij zelf nog te vertellen heeft.

speel het

Trek de regressielijn in SPSS — eerst zonder, dan met de vin

tweede voorspeller · zakkende helling · R²

Eerst zoals in het boek — de twaalf diamantjes

Analyze → Regression → Linear, glans als Dependent, karaat als Independent. OK — noteer b en R².
Zelfde analyse, nu met karaat én gladheid als Independent(s). Of gebruik Next om er een tweede blok van te maken.

Karaat zakt van 32.14 naar 0.98; R² gaat van .33 naar .82.

REGRESSION
  /STATISTICS COEFF R ANOVA CHANGE
  /DEPENDENT glans
  /METHOD=ENTER karaat
  /METHOD=ENTER gladheid.

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.

In SPSS · trek eraan

Open penguins.csv in SPSS: File → Open → Data, bestandstype op CSV.
Ga naar Analyze → Regression → Linear. Zet gewicht in Dependent en snavellengte in Independent(s). Draai zo eerst het model met alleen de snavel.
Nu de vin erbij: klik bij Independent(s) op Next en zet vinlengte in blok 2 — zo houd je de twee modellen naast elkaar. (Allebei in één blok mag ook; dan zie je alleen het tweede model.)
Klik op Statistics… en zet aan: Estimates, Model fit, R squared change en Confidence intervals op 95%. Dan zet SPSS in Model Summary vanzelf R Square Change, F Change en Sig. F Change neer — dat heb je zo nodig.

Of, als je liever typt dan klikt:

REGRESSION
  /STATISTICS COEFF R ANOVA CHANGE CI(95)
  /DEPENDENT gewicht
  /METHOD=ENTER snavellengte
  /METHOD=ENTER vinlengte.

Kijk 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 tabellen van onderen naar boven

F-toets · adjusted R² · standaardfout

Die vraag houden we nog even vast, want SPSS heeft zonet meer uitgerekend dan je gelezen hebt. Er staan drie tabellen in je output, en tot nu toe viste je er twee getallen uit: de helling en de R². In diezelfde tabellen staat óók alles wat een artikel over dit model zou rapporteren — je moet alleen weten in welke volgorde je ze 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. Hij staat in de tabel ANOVA: F(2, 339) = 536.63, p < .001. (SPSS noemt die p Sig. — zelfde getal, andere naam.)

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.

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.

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

Ga dan naar Model Summary, want daar staat een kolom die we nog niet hadden: Adjusted R Square. 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 de tabel Coefficients — en je hebt haar kolommen stuk voor stuk al eens ergens gezien:

B (onder Unstandardized Coefficients) — de hellingen zelf. De 48.14 van de vin en de 6.05 van de snavel. (De bovenste regel, (Constant), is het intercept: het voorspelde 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 — 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 getoonde 6.05 door de getoonde 5.18, dan krijg je 1.168 — SPSS rekent met álle decimalen, niet met de afrondingen.)
Sig. — 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 dat vinkje Confidence intervals geeft het laatste stuk: twee extra kolommen met 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 tabellen gelezen zoals een onderzoeker ze 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 Heus² — 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? Je zou vals willen spelen: de snavellengte netjes op z’n plek houden, maar elke vinlengte losknippen van z’n eigen pinguïn en op een willekeurige andere plakken.

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.

En hier moet ik eerlijk zijn: dat schudden is geen knop in SPSS. SPSS rekent één keurige uitkomst uit; het door elkaar husselen en duizend keer opnieuw kijken is code-werk. Dat is geen gemis van jou of van SPSS — het boek doet het bewust zo. Voor deze meervoudige regressie is er geen los schud-speeltje; wil je de duizend werelden echt zien dansen, draai dan het R-broertje, waar het drie regels code is.

Wat je in dat R-broertje ziet: 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 tabellen 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 in het R-broertje wiebelt 2.52 om nul, terwijl de standaardfout van diezelfde helling 2.01 was — dat is de SE van de vin uit Coefficients, 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. Maar waarom stortte de snavel dan in van 87.42 naar 6.05? Daar zit de les — en die kun je wél gewoon in SPSS nalezen.

snap het

Wat voegt elk van de twee nog toe?

modelvergelijking · ΔR² · hiërarchische regressie

Blijft er één ding jeuken. 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 met de twee blokken van zonet heb je hem al staan.

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:

In SPSS · wie voegt wat toe?

Je hebt hem al: R squared change onder Statistics… zet in Model Summary drie extra kolommen neer — R Square Change, F Change en Sig. F Change. Lees ze op de regel van blok 2, want die vergelijkt blok 2 mét blok 1.

Wat voegt de vin toe, bovenop de snavel? Blok 1 snavellengte, blok 2 vinlengte: R Square Change .406, F Change 573.09, df1 1, df2 339, en in Sig. F Change een waarde die SPSS afdrukt als .000. Let op dat laatste: .000 is geen p-waarde maar een afronding — een kans is nooit precies nul. In een zin schrijf je hem als p < .001.
En andersom: wat voegt de snavel toe, bovenop de vin? Draai de twee blokken om: R Square Change .001, F Change 1.36, Sig. F Change .244.

Lees 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².

En herken die twee getallen — .406 en .001 stonden hierboven al in de tabel bij de partiële correlatie, in de kolom uniek in R². Dat is geen toeval maar dezelfde som: de unieke bijdrage van een voorspeller aan de R², langs twee wegen uitgerekend. 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 kolom uit Model Summary, want die is verraderlijker dan hij eruitziet:

In SPSS · en wat doet de adjusted R²?

Kijk in de tweede opzet (vin in blok 1, snavel erbij in blok 2) naar de kolom Adjusted R Square van beide blokken. Voluit gerekend is dat .75828 tegen .75854; de tabel rondt af, dus op het scherm scheelt het hooguit één laatste decimaal.

De 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: de ruiskolom die wij trokken bracht hem van .75828 naar .75762.

Meestal, en dat woord staat er met opzet. Was het toeval die dag anders gevallen, dan had je in ongeveer één op de drie keer de adjusted R² juist zien stijgen — 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.

Achter dat “ongeveer één op de drie” zit een scherpe grens, en je hebt hem 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 terug naar de tabel: de snavel bracht F = 1.36 mee, net boven de 1 — en inderdaad kroop de adjusted R² van .75828 naar .75854. Onze ruiskolom kwam op F = 0.07, ver eronder, en daar zakte hij. Bij één voorspeller erbij komt dat op hetzelfde neer als t onder de −1 of boven de 1.

Zo is die “ongeveer één op de drie” ook precies te maken, zonder éé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.

Die grens van 1 is geen willekeurige plek, en waaróm hij daar ligt zie je in de volgende trede — daar meet je uit wat een voorspeller die niets weet gemiddeld oplevert, en blijkt de F precies de verhouding tot dát getal te zijn.

Wat is een schoenveter waard?

ruis · verwachte winst in R² · F als verhouding

Eerder in dit blok stond een belofte die altijd uitkomt: elke voorspeller die je toevoegt duwt de R² omhoog, ook een kolom puur toeval. Meet ieders schoenveterlengte op en zet hem in het model: je verklaart daarna méér van de verschillen in gewicht. Niet veel meer, maar meer, gegarandeerd. De vraag bij een nieuwe voorspeller is dus nooit óf hij het model verbetert — dat doet zelfs een schoenveter — maar of hij méér verbetert dan een schoenveter zou doen.

Het gras aan de overkant is met andere woorden altijd groener. Onthoud die zin; aan het eind van deze trede komen we erop terug, want hij is minder geruststellend dan hij klinkt.

En wat een schoenveter waard is, kun je gewoon uitrekenen. De vin haalde R² = .759, dus er ligt nog .241 van de verschillen in gewicht onverklaard. Een kolom die gegarandeerd niets weet, pakt daarvan gemiddeld één 340ste mee: (1 − .759) / 340 = .00071.

Die 340 is wat er van de 342 pinguïns overblijft nadat het vin-model zijn twee getallen geschat heeft, de helling en het intercept. En ja, vijf decimalen; hier leeft alles achter de komma. Dat getal is de bodemprijs: zóveel R² krijg je voor niets.

De snavellengte bracht, bovenop de vin, .00097 binnen. Deel dat door de bodemprijs — met álle decimalen en niet met de afrondingen — en er komt 1.36 uit. Herken je hem? F(1, 339) = 1.36, de modelvergelijking van hierboven. Dat is geen toeval: een F-waarde ís zo’n breuk. Onder de streep staat altijd de ruis — wat toeval sowieso oplevert — en boven de streep wat jouw voorspeller heeft opgehaald. Een F lees je dus als zoveel keer beter dan een schoenveter: de snavel deed het 1.36 keer zo goed als puur toeval. Daarom ligt de nullijn van een F ook niet op 0 maar op 1 — een voorspeller die niets weet, haalt gemiddeld gewoon de bodemprijs.

De bodemprijs naast de snavelwinst, op ware hoogte: .00071 tegen .00097. Dat het rechterstaafje 1.36 keer het linker is, is de F van hierboven — je ziet hem nu in plaats van hem uit te rekenen. (Voorlopige werkfiguur — vervangt Bens eigen tekening.)

Eén waarschuwing bij dat gemiddeld, want gemiddeld is hier niet hetzelfde als meestal. De winst van een toevalskolom is scheef verdeeld: de meeste schoenveters halen bijna niets, een enkele gelukstreffer haalt veel, en die treffers trekken het gemiddelde omhoog. Hoeveel dat scheelt — en hoeveel schoenveters de snavellengte verslaan — ga je hierna zelf meten.

En dan dat gras van hierboven, want er zit een addertje in de vraag “is het gras écht groener bij het grotere model?” Het antwoord is namelijk: ja, altijd, meetbaar. Dat was de belofte waarmee deze trede begon. Gróéner is de vraag dus niet. De vraag is waaróm het groener is: omdat de grond daar echt beter is, of omdat elke tuin er groener uitziet vanaf het plekje waar jij toevallig staat. Dat plekje is jouw steekproef — déze 342 pinguïns, met hun eigen rimpelingen van toeval. Je zag al dat de gewone R² zichzelf vleit: hij beschrijft jouw plekje. De adjusted R² keek voorzichtig over de heg, van jouw steekproef naar andere mogelijke steekproeven, naar de populatie. En de F-toets stelt die vraag op scherp: is de winst in verklaard gewichtsverschil groter dan wat het toeval van jouw plekje sowieso weggeeft? Groener in jouw steekproef is nog geen groener in de populatie.

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

Meet een voorspeller die niets kán weten

ΔR² van pure ruis · F als verhouding · nullijn

Jouw beurt — hetzelfde recept als hierboven, maar nu op een voorspeller waarvan je zéker weet dat hij niets te melden heeft. Je maakt er zelf een: een kolom willekeurige getallen, één per pinguïn. Noem hem de schoenveter.

In SPSS maak je zo’n kolom met Transform → Compute Variable: als doelvariabele schoenveter, en als expressie RV.NORMAL(0,1). Zet er SET SEED 8. vóór als je precies onze getallen wilt terugzien. De modelvergelijking draai je zoals in dit blok: zet de vinlengte in Block 1 en de schoenveter in Block 2, vink onder Statistics de optie R squared change aan, en lees in Model Summary de kolommen R Square Change, F Change en Sig. F Change.

En hier moet ik weer eerlijk zijn: dat duizend keer overdoen is geen knop in SPSS. Eén schoenveter kun je hier prima maken en toetsen, en dat doe je bij a) tot en met e). Maar duizend verse schoenveters trekken en tellen hoeveel er beter zijn dan de snavellengte, dat is code-werk. Bij f) tot en met h) geven we je daarom wat het R-broertje eruit haalt, met de vraag om het te lezen — en wie het zelf wil zien draaien, draait het daar.

Tussen haakjes staat telkens wat er in een compleet antwoord hoort, zodat je jezelf kunt nakijken zonder ons antwoord te openen. De laatste heet weer Z): de aannames staan altijd achteraan.

Jouw beurt

a) Wat doet de gewone R² als de schoenveter erbij komt? Kán hij omlaag? (rapporteer beide R²’s op vijf decimalen, en één zin over dat “kan”)

b) En wat doet de adjusted R²? (rapporteer beide, en zeg welke kant hij op ging)

c) Toets de twee modellen tegen elkaar. (rapporteer ΔR², F(df₁, df₂) en p)

d) De snavellengte voegde bovenop de vinlengte ΔR² = .00097 toe, met F(1, 339) = 1.36. Zet jouw schoenveter daarnaast. Welke van de twee doet het beter? (twee ΔR²’s en twee F’s, en één zin)

e) Je hoeft niet duizend keer te trekken om te weten wat ruis gemiddeld waard is; dat is uit te rekenen. Onder de nulhypothese levert één willekeurige voorspeller gemiddeld (1 − R²) / df op, met de R² en de residuele vrijheidsgraden van het basismodel. Reken die som uit. (rekenmachine erbij, het is één deling)

f) In het R-broertje zijn duizend verse schoenveters getrokken. Hun ΔR² is gemiddeld .00072 en de mediaan is .00030. Waarom liggen die twee niet op elkaar, en wat betekent dat voor het woord gemiddeld? (één zin over de vorm van de verdeling)

g) Van diezelfde duizend halen er 254 de .00097 van de snavellengte. Schrijf op wat dat getal zegt. (één zin, in gewone woorden)

h) Op welk percentiel van pure ruis staat de snavellengte daarmee — oftewel: hoeveel procent van de schoenveters blijft ónder haar? (één getal, en één zin over wat dat betekent voor een maat die in haar eentje R² = .354 haalde)

i) Deel de .00097 van de snavellengte door de verwachte ruiswinst uit e). Vergelijk de uitkomst met de F uit dit blok. Waarom ligt de nullijn van een F dus op 1 en niet op 0? (één deling, en twee zinnen)

j) Een schoenveter maakt de fit altijd beter — het gras aan de overkant is dus echt groener. Waaróm is dat geen reden om erheen te lopen? (één zin, met het woord populatie of steekproef erin)

k) Geef een inhoudelijke conclusie in niet meer dan vijf zinnen. (over pinguïns en over schoenveters, niet over getallen)

Z) Mag je alles hierboven geloven? Loop de aannames na:

heeft elke pinguïn een gewicht en een vinlengte, of vallen er rijen weg?
liggen de residuen ongeveer symmetrisch?
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?

Antwoord — open na je eigen poging

a) Van .75899 naar .75903. Omhoog dus, met drie honderdduizendsten — en nee, omlaag kán niet. Een voorspeller erbij geeft het model één extra knop om aan te draaien; in het slechtste geval zet SPSS die knop op nul en verandert er niets. Zakken kan hij daardoor nooit, ook niet van een kolom die per definitie niets weet. (Trek je je eigen schoenveter zonder onze seed, dan krijg je andere getallen — maar omhoog gaat hij altijd.)

b) Van .75828 naar .75761: omlaag. Dat is precies wat de adjusted R² hoort te doen bij een voorspeller die zijn plek niet betaalt. Kijk terug naar de scherpe grens uit dit blok: de adjusted R² stijgt als de F boven de 1 uitkomt en zakt eronder. Onze schoenveter komt bij c) uit op F = 0.05, ver onder de 1 — en dus zakt hij.

c) ΔR² = .00003, F(1, 339) = 0.05, p = .83. Niets aan de hand, zoals het hoort: de schoenveter kán niets weten, en dat is precies wat de toets zegt.

d) De snavellengte: ΔR² = .00097, F(1, 339) = 1.36. De schoenveter: ΔR² = .00003, F(1, 339) = 0.05. De snavellengte wint van déze schoenveter, met ruim dertig keer zoveel. Onthoud dat “ruim dertig keer”, want bij g) blijkt het minder indrukwekkend dan het klinkt — één schoenveter is nog geen lat.

e) (1 − .7590) / 340 = .00071. De 340 is de residuele vrijheidsgraden van het basismodel: 342 pinguïns min de twee getallen die het model met alleen de vinlengte schat, de helling en het intercept.

f) De verdeling is scheef: de meeste schoenveters leveren bijna niets op, en een handvol toevalstreffers in de rechterstaart tilt het gemiddelde op. Let op wat dat betekent voor het woord gemiddeld: wat ruis gemiddeld oplevert, is méér dan wat ruis meestal oplevert.

g) Ruim één op de vier willekeurige getallen doet het beter dan de snavellengte.

h) Het 75e percentiel. Een maat die in haar eentje R² = .354 haalde — ruim een derde van alle verschil in gewicht — is naast de vinlengte niet meer waard dan een bovengemiddelde schoenveter.

En kijk nog even naar die 254 van de 1000, oftewel .254. Vergelijk hem met de p = .244 uit de coëfficiëntentabel. Dat scheelt bijna niets, en dat is geen toeval: allebei beantwoorden ze dezelfde vraag — hoe vaak zou toeval dit namaken? De tabel rekent hem uit langs de F-verdeling, jij hebt hem geteld. Dat is precies wat de schud-trede eerder in dit blok ook deed, alleen viel de uitkomst daar in de staart en hier midden in het bergje.

i) .00097 / .00071 = 1.36, en dat is de F uit dit blok. Geen toeval: een F is een verhouding tussen wat de voorspeller oplevert en wat toeval zou opleveren. Bij F = 1 doet hij het dus precies even goed als een willekeurig getal, en dáár ligt de nullijn — niet op 0, want zelfs niets-weten levert een beetje op.

Detail dat je gerust mag laten liggen: op vier decimalen zijn de twee niet helemaal gelijk. Jouw deling geeft 1.3616 en de echte F geeft 1.3631. De F deelt namelijk door de ruis die overblijft ná beide voorspellers, en dat is een fractie minder ruis dan na de vin alleen — een kleinere noemer, dus een iets grotere uitkomst.

j) Omdat het gras aan élke overkant groener is. De R² gaat van elke voorspeller omhoog, dus “hij past beter” onderscheidt een echte voorspeller niet van een schoenveter. De vraag is of het gras daar groener is omdát het gras er beter is, of omdat jij toevallig op dít plekje staat: geldt de winst in de populatie, of alleen in jouw steekproef? Dat is wat de F en de p beantwoorden en de kale R² niet.

k) Bijvoorbeeld: bij 342 pinguïns voorspelt de vinlengte het gewicht sterk. Een kolom willekeurige getallen erbij maakt het model op papier iets beter, maar niet meer dan toeval verklaart. Ook de snavellengte voegt naast de vinlengte niets toe wat van toeval te onderscheiden is: ruim een kwart van duizend willekeurige getallen deed het beter. Een betere fit is dus geen bewijs dat een voorspeller iets weet. Wat je wilt weten is of de winst groter is dan wat niets weten al oplevert.

Z) Ja, met dezelfde slag om de arm als bij het model hierboven. Alle 342 pinguïns hebben een gewicht en een vinlengte, dus er valt geen rij weg en de twee modellen draaien op dezelfde beesten — dat is hier extra belangrijk, want twee modellen op verschillende steekproeven mag je niet van elkaar aftrekken. De residuen liggen vrijwel symmetrisch (gemiddelde min mediaan is 0.06 standaardafwijking). En de schoenveter hangt met de vinlengte samen op r = .06, dus ze zitten elkaar niet in de weg — ook dat hoort zo, want een willekeurig getal hoort met niets samen te hangen.

Gaat er iets mis? Dan valt dit rijtje in tweeën, en de scheidslijn loopt op een verrassende plek.

Blijft staan: de R²’s bij a), de ΔR²’s bij c) en d), en — dit is de verrassing — de hele telling bij f), g) en h). Die eerste zijn beschrijvingen van deze 342 pinguïns, en daar verandert geen aanname iets aan. Maar ook je duizend schoenveters blijven staan, en wel omdat je de ruis zélf hebt gemaakt: je hébt geteld hoe vaak toeval het beter deed, je hebt het niet uit een verdeling afgeleid. Datzelfde gold voor de schud-trede eerder in dit blok.
Doe je over: de F en de p bij c) en d), de verwachte ruiswinst bij e) — en daarmee ook de deling bij i), die daarop leunt — en de conclusie bij k). Die komen allemaal uit een verdeling die aanneemt dat de residuen zich netjes gedragen, en dat is precies wat je bij p = .83 aan het geloven bent.
Blijft staan met een slag om de arm: de adjusted R² bij b). Het getal klopt hoe dan ook, maar zijn hele punt is dat hij een gok doet over de populatie, en zo’n gok is niet beter dan de aannames eronder.
Bij j) staat geen getal, dus daar valt niets over te doen.

En zie wat daar gebeurt: de geschudde en de getelde antwoorden overleven, de uitgerekende niet. Dát is waarom dit boek je eerst laat tellen.

rapportage

Schrijf het op als een onderzoeker

rapportagezin · F, b en R² · twee decimalen

Rapportage

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 tabellen; 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. De partiële correlatie tussen gewicht en snavellengte, gecontroleerd voor vinlengte, was r = .06, tegen een gewone correlatie van r = .60.

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 tabellen 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 lagen op voor het model uit Jouw beurt hierboven — snavellengte en snaveldiepte samen. Het model als geheel (R² = .47, F(2, 339) = 150.82, p < .001) en dan allebei de voorspellers apart (b = 75.28, SE = 5.97, t(339) = 12.61, p < .001, 95% CI [63.54, 87.03]; en b = −142.72, SE = 16.51, t(339) = −8.65, p < .001, 95% CI [−175.19, −110.25]). Durf het minteken erin te zetten en het uit te schrijven: bij gelijke snavellengte hoort een diepere snavel bij een lichter beest. Hardop voorlezen mag.

◠ Zelfde vorm in het boek

Hier: één helling die instort zodra je een tweede voorspeller binnenlaat, de confounder die daarachter zat, en de partiële correlatie die het in één getal vangt.

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