Tweede beweging · Verschillen koppelen

Oefening 3.3 · De residuen bekijken

Wat de lijn liet liggen, en waar R² vandaan komt · residuals

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

prikkel

Tweeëntwintig pinguïns van precies 190 millimeter

voorspelling · misser

De pinguïns zijn oude bekenden — in de vorige oefening trok je er zelf een lijn door: het gewicht voorspeld uit de vinlengte, bijna vijftig gram pinguïn per millimeter vin, en R² = .76.¹ Driekwart van alle verschil in gewicht loopt langs de vin mee. Dat is veel, voor levende dieren, en het klinkt alsof de lijn het gewicht daarmee wel zo ongeveer wéét.

En nu iets raars. In het bestand zitten tweeëntwintig pinguïns met een vin van precies 190 millimeter. De lijn kent alleen de vinlengte, dus die tweeëntwintig krijgen alle tweeëntwintig dezelfde voorspelling: 3659.4 gram. Wat de weegschaal werkelijk zei, in grammen:

3050 · 3325 · 3400 · 3450 · 3450 · 3450 · 3575 · 3600 · 3600 · 3650 · 3700 · 3700 · 3725 · 3800 · 3900 · 3900 · 3900 · 3950 · 4250 · 4250 · 4250 · 4600

Van 3050 tot 4600 — anderhalve kilo verschil, bij één en dezelfde vinlengte. De lijn kan daar niets aan doen: hij weet maar één ding van een pinguïn, en dat ene ding was tweeëntwintig keer gelijk.

En die tweeëntwintig zijn geen rariteit die je toevallig trof. Kies een andere vinlengte en je ziet hetzelfde, alleen met minder dieren tegelijk op een kluitje: overal langs de lijn zit de werkelijkheid erboven en eronder. Dat zet een vraag op scherp die we tot nu toe hebben laten liggen: hoe ver zit de lijn er per pinguïn naast — en wat kun je met die missers doen?

Trek even mee. De lijn heeft gegokt; wij gaan nakijken.

speel het

Reken de missers uit op de twaalf steentjes

residu · voorspelde waarde

Zo’n misser heeft een naam, en je hebt hem in het boek al zien staan: DATA = FIT + RESIDU. De fit is wat de lijn gokt, het residu is wat er overblijft — de gokfout, nu niet meer ten opzichte van het gemiddelde maar ten opzichte van de lijn.

Eerst zoals in het boek — de twaalf diamantjes

d <- read.csv("diamantjes.csv")
model.karaat <- lm(glans ~ karaat, data = d)

# wat de lijn gokt, en wat er overblijft
round(fitted(model.karaat), 1)     # → 30.7 33.9 33.9 40.4 50.0 50.0 53.2 56.4 59.6 59.6 66.1 66.1
round(resid(model.karaat),  1)     # → -10.7 -3.9 -3.9 -10.4 0.0 20.0 26.8 3.6 10.4 20.4 -36.1 -16.1

# de kolom die je in het boek al hebt zien optellen tot nul
sum(resid(model.karaat))           # → 0 (op afrondingsgruis na)

Kijk of die som ook bij jou nul is. Doet hij dat niet, dan heb je onderweg iets anders in handen gekregen dan de residuen van dit model.

Dat zijn precies de twee kolommen uit de tabel in het boek: Else valt op nul (de lijn gaat dwars door haar heen), Kees zit er met −36.1 het verst naast. Nieuw is alleen de derde regel hieronder — je zet die twee kolommen tegen elkaar in een plaatje in plaats van ze te lezen:

plot(fitted(model.karaat), resid(model.karaat))    # de gok op de x-as, de misser op de y-as
abline(h = 0)                # de lijn zelf ligt hier plat op nul

Kijk goed wat dat plaatje doet. De regressielijn is plat geworden: alles wat de lijn al wist is eruit gerekend, dus wat overblijft hoort rond de nullijn te zweven. Een residuenplot is een schoongeveegde puntenwolk — je kijkt niet meer naar het verband, je kijkt naar wat het verband liet liggen.

Eerst de twaalf diamantjes

Twee vragen voor je verder leest.

a) Vraag de standaardafwijking van de glans op, en die van de residuen. (twee getallen, en één zin over wat het verschil tussen die twee betekent)

b) Kwadrateer allebei die getallen, tel de twaalf gekwadrateerde afwijkingen van de glans op (sum((d$glans - mean(d$glans))^2)) en doe hetzelfde met de residuen. (twee sommen, en zeg welk deel van de eerste door de tweede wordt opgeslokt)

Antwoord — open na je eigen poging

a) De glans zelf spreidt met SD = 21.74 glanspunten, de residuen met SD = 17.75. Zonder lijn is je beste gok voor elk steentje het gemiddelde, en dan zit je er typisch bijna 22 punten naast; met de lijn nog bijna 18. De lijn scheelt je dus wel iets, maar niet veel — bij deze twaalf steentjes is het karaat lang niet alles.

sd(d$glans)                    # → 21.74
sd(resid(model.karaat))        # → 17.75

b) De glans om zijn eigen gemiddelde geeft 5200.0; de glans om de lijn geeft 3464.3. Van alle verschil in glans blijft er dus tweederde over nadat de lijn zijn deel heeft gehad — en dat ene derde dat hij wél pakt, is precies de R² = .33 die je in de vorige oefeningen zag.

sum((d$glans - mean(d$glans))^2)    # → 5200
sum(resid(model.karaat)^2)          # → 3464.3

Waaróm dat zo is, staat hieronder.

snap het

Het residu, en waar R² vandaan komt

residu · residuenplot · R²

Kijk nog eens naar die som van de residuen. Nul — niet ongeveer, maar op afrondingsgruis na exact. Dat is geen toeval en ook geen compliment aan de lijn: het is de bouw. Zo is de lijn gelegd — voor elke misser in de plus staat er evenveel min tegenover, want bleef er iets over, dan kon de hele lijn een stukje opschuiven en was het de best passende lijn niet.

Zoek voor jezelf momentjes om te controleren — het zou nul moeten zijn, klopt dat? Onthoud die nul dus wél als controle-momentje: telt jouw kolom residuen niet op tot nul, dan is er onderweg iets misgegaan.

Maar let op wat het óók betekent. Juist omdát die som altijd nul is — bij een lijn die strak door de wolk loopt én bij een lijn die er van alles naast laat liggen — zegt hij níéts over hoe goed de lijn gokt. Plus en min heffen elkaar toch wel op. Wil je weten hoe goed, dan moet je niet optellen maar kijken: hoe breed liggen de missers om de nul?

En dat is precies wat je in de residuenplot ziet. De tweeëntwintig waar dit blok mee opende staan er als één kolom in: zelfde gok, tweeëntwintig missers, boven en onder de nul.

In R · teken de missers en tel de tuin

peng <- read.csv("penguins.csv")
model.vl <- lm(gewicht ~ vinlengte, data = peng)

# de schoongeveegde wolk: gok op de ene as, misser op de andere
plot(fitted(model.vl), resid(model.vl))
abline(h = 0)

# je gokfout zónder lijn, en wat er mét lijn van overblijft
sd(peng$gewicht)                                 # → 801.95
sd(resid(model.vl))                              # → 393.70

# kwadrateer eerst, deel dan
1 - var(resid(model.vl)) / var(peng$gewicht)     # → 0.759

Zonder lijn was je beste gok voor elk gewicht het gemiddelde — dat spel ken je — en dan was je typische misser 801.95 gram, de standaardafwijking van het gewicht. Mét de lijn blijft er 393.70 gram over: de standaardafwijking van de missers. De lijn scheelt je dus 408.25 gram, ruim de helft van je gokfout.

En toch zei de vorige oefening R² = .76. Driekwart. Wie liegt er — de helft of de driekwart? Geen van beide; ze tellen in een andere maat. Denk aan alle verschil in gewicht als een tuin vol variatie. Een verband verklaren is geen lijntje verklaren maar oppervlakte, en oppervlakte gaat in kwadraten — dus kwadrateer de twee gokfouten vóór je ze op elkaar deelt. Dan staat er .76: het verklaarde stuk van de tuin, gedeeld door de hele tuin. In meters langs de rand wint de lijn ruim de helft; in tuinoppervlakte driekwart. Het is dezelfde winst, twee keer opgemeten.

En R² heeft daar niets nieuws voor nodig: alles wat hij weet, zit al in de missers. Doe dezelfde deling op de twaalf steentjes en er komt .33 uit — precies het kwadraat van de r = .58 die je aan het begin van dit hoofdstuk zelf in elkaar zette. Twee routes, één getal: de correlatie kwadrateren, of de missers boekhouden.

jouw beurt

De boekhouding van de pinguïnlijn

som van de residuen · residuenplot · R² uit de missers

Nu jij, op de lijn uit de prikkel: het gewicht van 342 pinguïns, voorspeld uit hun vinlengte. Je hebt hem in de vorige oefening getrokken en gelezen; nu ga je uitrekenen wat hij liet liggen.

In R · zet klaar

peng <- read.csv("penguins.csv")
model.vl <- lm(gewicht ~ vinlengte, data = peng)

Jouw beurt

Beantwoord de vragen op volgorde; ze bouwen op elkaar. De laatste heet met opzet Z) en niet k): de wolk-vraag staat altijd achteraan.

Het zijn er elf, en dat is te veel voor één zitting. Na e) heb je de missers in handen en heb je de prikkel teruggevonden — een goed moment om te stoppen. Bij f) begint de boekhouding, en dat is het stuk waar R² uit rolt.

a) Vraag van de eerste zes pinguïns de voorspelling en het residu op (head(fitted(model.vl), 6) en head(resid(model.vl), 6)). Tel bij één van de zes de twee op en leg de uitkomst naast zijn echte gewicht. (twee rijtjes van zes, en één som)

b) Tel alle 342 residuen op. Wat komt eruit, en waarom is dat géén compliment aan de lijn? (één getal, en twee zinnen)

c) Teken de residuenplot. Wat hoort daar te staan, en wat zie je? (twee zinnen)

d) Zoek de pinguïn waar de lijn er het verst naast zit, naar boven en naar onder. (twee getallen mét eenheid, en zeg van welke soort ze zijn)

e) Terug naar de tweeëntwintig uit de prikkel. Vraag de residuen op van de pinguïns met een vin van precies 190 millimeter. (de kleinste en de grootste, en één zin over wat er in de residuenplot boven die ene gok staat)

f) Nu de boekhouding. Tel de gekwadrateerde afwijkingen van het gewicht om zijn eigen gemiddelde op, en die van de residuen om nul. Tel daarna de gekwadrateerde afstanden van de lijn tot het gemiddelde op (sum((fitted(model.vl) - mean(peng$gewicht))^2)). (drie getallen — en tel de laatste twee bij elkaar op)

g) Reken R² uit die getallen: één min de residu-som gedeeld door de totaal-som. Leg het naast summary(model.vl)$r.squared. (twee getallen)

h) Vraag de standaardafwijking van het gewicht op en die van de residuen. Hoeveel procent van je gokfout scheelt de lijn — en waarom is dat een ánder getal dan de R² uit g)? (twee getallen, en twee zinnen)

i) Kleur de residuenplot per soort en vraag het gemiddelde residu per soort op. Wat zegt dat over de lijn? (drie getallen, en twee zinnen)

j) Schrijf je bevinding op in één zin. (wat de lijn dekt, en wat hij per pinguïn laat liggen)

Z) Is de wolk recht, en zitten er geen punten tussen die de boel verneuken? Kijk naar je residuenplot uit c) — ligt er een boog in, of hangt er één dier zo ver weg dat het in zijn eentje aan de lijn trekt? (twee korte uitspraken, en zeg welke van je antwoorden hierboven eraan hangen)

Antwoord — open na je eigen poging

a) De zes voorspellingen zijn 3212.3 · 3460.7 · 3907.9 · 3808.5 · 3659.4 · 3212.3 gram, en de zes residuen 537.7 · 339.3 · −657.9 · −358.5 · −9.4 · 412.7. Tel ze per dier op en je hebt het gewicht terug: de eerste pinguïn weegt 3212.3 + 537.7 = 3750 gram, en dat staat ook in het bestand. DATA = FIT + RESIDU, letterlijk.

round(head(fitted(model.vl), 6), 1)   # → 3212.3 3460.7 3907.9 3808.5 3659.4 3212.3
round(head(resid(model.vl), 6), 1)    # → 537.7 339.3 -657.9 -358.5 -9.4 412.7
head(peng$gewicht, 6)                 # → 3750 3800 3250 3450 3650 3625

b) Nul, op afrondingsgruis na. En dat is geen compliment: het is de manier waarop de lijn gelegd is. Elke misser in de plus heeft een misser in de min tegenover zich staan, en dat geldt voor élke best passende lijn — ook voor een die er compleet naast zit. Een getal dat altijd hetzelfde is, kan je niets vertellen over hoe goed het gaat.

sum(resid(model.vl))   # → 0 (op afrondingsgruis na)

c) Er hoort een vormeloze wolk te staan rond de nullijn: alles wat de lijn al wist is eruit gerekend, dus wat overblijft mag geen patroon meer dragen. En dat is ook wat je ziet — een band van punten die van links naar rechts ongeveer even breed blijft en nergens duidelijk omhoog of omlaag loopt. Je ziet de verticale strepen wél: pinguïns zijn op hele millimeters opgemeten, dus de gokken vallen op een beperkt aantal plekken.

d) Naar boven zit de lijn er 1288.69 gram naast (een Adélie, rij 39) en naar onder 1058.80 gram (een Chinstrap, rij 313). Ruim een kilo, allebei — bij dieren die gemiddeld ruim vier kilo wegen.

max(resid(model.vl))              # → 1288.69
min(resid(model.vl))              # → -1058.80
peng$soort[which.max(resid(model.vl))]   # → "Adelie"

e) Van −609 tot +941 gram. In de residuenplot staat boven die ene gok van 3659.4 dus één verticale streep van tweeëntwintig punten, ruim anderhalve kilo lang. Dat is de prikkel, nu als plaatje: dezelfde voorspelling, tweeëntwintig verschillende missers.

round(range(resid(model.vl)[peng$vinlengte == 190]))   # → -609  941

f) De totaal-som is 219 307 697, de residu-som 52 854 796 en de fit-som 166 452 902. Tel die laatste twee op en je krijgt de eerste terug, tot op de eenheid. Dat is DATA = FIT + RESIDU, maar dan in gekwadrateerde afwijkingen — en dát is de vorm waarin de boekhouding klopt.

sum((peng$gewicht - mean(peng$gewicht))^2)        # → 219307697
sum(resid(model.vl)^2)                            # → 52854796
sum((fitted(model.vl) - mean(peng$gewicht))^2)    # → 166452902

g) .76, allebei.

1 - sum(resid(model.vl)^2) / sum((peng$gewicht - mean(peng$gewicht))^2)   # → 0.759
summary(model.vl)$r.squared                                               # → 0.759

En dat is de pointe van deze oefening. R² is geen aparte maat die R ergens vandaan tovert: het is de tuin die de lijn pakt, gedeeld door de hele tuin. Je hebt hem net zelf uit de missers gerekend.

h) De standaardafwijking van het gewicht is 801.95 gram, die van de residuen 393.70. De lijn scheelt je dus 408.25 gram, en dat is 51 procent van je gokfout — ruim de helft, niet driekwart.

Twee zinnen waarom die twee getallen allebei kloppen: de 51 procent telt in grammen, langs de rand van de tuin, en de .76 telt in oppervlakte. Kwadrateren maakt van ruim de helft driekwart — reken maar na, 393.70 gedeeld door 801.95 is .49, en .49 in het kwadraat is .24, en één min .24 is .76.

sd(peng$gewicht)          # → 801.95
sd(resid(model.vl))       # → 393.70

i) Adélie +43.5, Chinstrap −215.7, Gentoo +65.8 gram.

plot(fitted(model.vl), resid(model.vl), col = factor(peng$soort), pch = 16)
abline(h = 0)
tapply(resid(model.vl), peng$soort, mean)    # → 43.5  -215.7  65.8

Twee zinnen: over alle 342 heen is het gemiddelde residu nul, maar per soort is het dat niet — de lijn zit voor de Chinstraps stelselmatig ruim twee ons te hoog en voor de andere twee een beetje te laag. Dat wil zeggen dat er ook hier nog een derde speler in de wolk zit, precies zoals in de eerste oefening van dit hoofdstuk; alleen is hij deze keer klein genoeg om de lijn niet om te gooien.

j) Bijvoorbeeld:

De lijn van gewicht op vinlengte dekte ruim driekwart van de verschillen in gewicht bij 342 pinguïns (R² = .76), maar liet per dier nog altijd bijna vierhonderd gram spreiding om zich heen liggen (SD van de residuen = 393.70).

Z) Ja — allebei komen ze goed uit.

Er ligt geen boog in. De residuen lopen niet stelselmatig omhoog of omlaag met de voorspelling; sterker nog, de samenhang tussen de gok en de misser is exact nul, en dat is geen bevinding maar de bouw — dat is precies wat “de lijn is plat geworden” betekent. Wat je wél kunt nakijken is of er een kromme in zit, en dat doe je door de lijn een knik toe te staan. De R² was .76, op vier decimalen .7590; met de knik erbij wordt het .7753. Ruim anderhalve procent erbij — maar dat is geen bocht, dat is een parameter erbij, en die kan de R² alleen maar omhoog duwen.

En geen pinguïn trekt in zijn eentje aan de lijn. De grootste misser is 1288.69 gram, en dat is ruim drie keer de standaardafwijking van de residuen — een dier dat opvalt, maar niet een dier dat de lijn kantelt. Laat je hem weg, dan beweegt de helling nauwelijks.

cor(fitted(model.vl), resid(model.vl))   # → 0 (op machineruis na)

Wat eraan hángt, als het wél mis was gegaan: de hele boekhouding blijft staan — de sommen bij f), de R² bij g), de standaardafwijkingen bij h) en de soortgemiddelden bij i) zijn optellingen op deze 342 dieren, en die kloppen wat er ook krom is. Wat je overdoet is elke zin waarin je de lijn aanbeveelt: bij j) staat dat hij driekwart dekt, en als het verband krom loopt is dat geen leugen maar wel een lijn die je op de verkeerde plek in de steek laat.

eindopdracht

En nu bij mensen — dezelfde boekhouding op 2634 vragenlijsten

schaalscore · R² uit de missers · steekproefzin

Eindopdracht

De afsluiter van het hele hoofdstuk, en hij eindigt waar dit boek altijd eindigt: bij mensen.

In bigfive.csv staan 2634 mensen die vijfentwintig uitspraken over zichzelf beoordeelden op een schaal van 1 tot 6. Je gebruikt er hier tien — vijf over neuroticisme, hoe snel je van slag raakt, en vijf over openheid, de smaak voor nieuwe dingen. Je hebt ze in de vorige twee oefeningen al ontmoet; nu doe je er de boekhouding van dit blok op.

b <- read.csv("bigfive.csv")

b$neuroticisme <- rowMeans(b[, c("N1", "N2", "N3", "N4", "N5")])
b$openheid     <- rowMeans(b[, c("O1", "O2", "O3", "O4", "O5")])

model.mens <- lm(neuroticisme ~ openheid, data = b)

a) Vraag de residuen op en tel ze op. (één getal, en één zin over waarom je dat had kunnen voorspellen)

b) Vraag de standaardafwijking van het neuroticisme op en die van de residuen. Hoeveel scheelt de lijn je? (twee getallen, en één zin)

c) Reken R² uit de missers, net als bij de pinguïns. (één getal, en leg het naast summary(model.mens)$r.squared)

d) Nu de mensen-versie van de tweeëntwintig pinguïns. Hoeveel mensen scoren precies 4.0 op openheid, welke voorspelling krijgen die allemaal, en hoe ver loopt hun échte neuroticisme uiteen? (drie antwoorden, en één zin)

e) Schrijf op wat je gevonden hebt, in twee zinnen en in deze vaste vorm:

Binnen deze n mensen geldt: … (R² = …). Of dat ook buiten deze n opgaat, weet je nog niet.

(vul de n, de bevinding in gewone taal en de R² in — en laat die tweede zin staan zoals hij er staat)

Antwoord — open na je eigen poging

a) Nul, op afrondingsgruis na — en dat had je kunnen voorspellen, want dat geldt voor élke best passende lijn, hoe goed of slecht hij ook past. Bij de diamantjes was het zo, bij de pinguïns was het zo, en hier dus ook.

sum(resid(model.mens))   # → 0 (op afrondingsgruis na)

b) Het neuroticisme spreidt met SD = 1.1946 schaalpunten, de residuen met SD = 1.1789. De lijn scheelt je dus 0.0156 schaalpunt op een schaal van 1 tot 6 — oftewel: bijna niets. (Vier decimalen hier met opzet: op twee zie je het verschil niet eens meer.) Wie de openheid van iemand kent, weet daarmee nauwelijks meer over zijn neuroticisme dan wie alleen het gemiddelde kent.

sd(b$neuroticisme)          # → 1.1946
sd(resid(model.mens))       # → 1.1789

c) .03, allebei. Van de hele tuin pakt deze lijn nog geen drie procent.

1 - var(resid(model.mens)) / var(b$neuroticisme)   # → 0.026
summary(model.mens)$r.squared                      # → 0.026

d) 470 mensen scoren precies 4.0 op openheid. Die krijgen alle 470 dezelfde voorspelling: 3.21 schaalpunten neuroticisme. En hun échte neuroticisme loopt van 1.0 tot 6.0 — van de laagste score die op deze schaal bestaat tot de hoogste, allemaal binnen die ene gok.

sum(b$openheid == 4)                                  # → 470
round(predict(model.mens, data.frame(openheid = 4)), 2)   # → 3.21
range(b$neuroticisme[b$openheid == 4])                # → 1  6

Dat is dezelfde figuur als de tweeëntwintig pinguïns waar dit blok mee opende, maar dan harder: daar liep het echte gewicht over anderhalve kilo uiteen bij één gok, hier loopt het echte neuroticisme over de hele schaal uiteen bij één gok. Bij mensen is dat de regel en niet de uitzondering.

e) Bijvoorbeeld:

Binnen deze 2634 mensen geldt: wie opener is, scoort gemiddeld iets hoger op neuroticisme, maar openheid en neuroticisme delen nog geen drie procent van hun verschillen (R² = .03). Of dat ook buiten deze 2634 opgaat, weet je nog niet.

Kijk wat die tweede zin doet. Alles wat je in dit hoofdstuk hebt uitgerekend — de covariantie, de correlatie, de lijn, de residuen, de R² — beschrijft de mensen die je vóór je hebt, en verder niemand. Je hebt geen enkele reden gekregen om aan te nemen dat de volgende 2634 mensen hetzelfde zouden laten zien, en ook geen manier om te zeggen hoe waarschijnlijk dat is.

Dat is precies waar de volgende beweging mee begint.

◠ Zelfde vorm in het boek

Hier: de missers zelf uitgerekend, en R² uit de boekhouding in plaats van uit een knop.

In het boek: W4 — Hangen ze samen?

Daar valt DATA = FIT + RESIDU voor het eerst; hier zie je wat er in die derde term zit, en dat het verklaarde deel er gewoon uit te tellen is.

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