Derde beweging · Verschil op de proef

Oefening 5.1 · Het betrouwbaarheidsinterval

Het verschil tussen twee gemiddelden, geschat mét zijn onzekerheid · confidence interval

Data bij deze oefening — plantgroei.csv plantgroei onder drie condities · diamantjes.csv de twaalf steentjes uit het boek · wat elke kolom betekent

prikkel

Hoe zeker is die 0.87 gram?

schatting uit een steekproef · onzekerheid

Dertig potten in een kas, dertig dezelfde plantjes.¹ Tien kregen gewoon water, tien kregen plantenvoer één, tien plantenvoer twee. En toen — dit is het onsentimentele deel — werden alle dertig geknipt, gedroogd en gewogen. Droog, want nat gewicht is vooral water, en water zegt niks over groeien.

De controle-potten laten we vandaag even links staan; het gaat ons om de twee soorten voer. Voer één: gemiddeld 4.66 gram droge plant. Voer twee: 5.53 gram. Verschil: 0.87 gram, in het voordeel van voer twee.

En nu níét de vraag die je verwacht. Niet “is er verschil, ja of nee” — dat is de bijvraag. De echte vraag is: hoe zeker is die 0.87? Er stonden maar tien potten per groep. Had de kweker tien ándere potten gepakt, was het dan 0.9 geweest? 0.5? Nul komma niks? Een getal uit een steekproef is een schatting, en een schatting zonder onzekerheids-bandje eromheen is een mening met cijfers. Dat bandje gaan we maken.

speel het

Reken het verschil uit in R — mét het bandje

betrouwbaarheidsinterval aflezen · zit de nul erin

Eerst zoals in het boek — de twaalf diamantjes

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

# t.test geeft het interval er gratis bij
t.test(d$glans)$conf.int    # → 36.19 tot 63.81

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

Het 95%-betrouwbaarheidsinterval voor de gemiddelde glans loopt van 36.19 tot 63.81 — bijna dertig punten breed. Dat is geen slordigheid maar eerlijkheid: twaalf steentjes weten gewoon niet zo veel, en het interval zegt dat hardop.

Nu jij — met nieuwe data: plantjes onder drie condities, en de vraag hoe zeker dat verschil van 0.87 gram eigenlijk is.

Haal de potten binnen, pak de twee voer-groepen eruit, en laat R het verschil uitrekenen — mét het bandje, want dat doet t.test er gratis bij.

In R · trek eraan

pg <- read.csv("plantgroei.csv")
voer <- pg[pg$groep %in% c("behandeling1", "behandeling2"), ]

tapply(voer$gewicht, voer$groep, mean)   # → 4.66 en 5.53
# FALSE is de standaard: Welch, het huismerk van deze boekjes
t.test(gewicht ~ groep, data = voer, var.equal = FALSE)   # → t = -3.01, p = .009
                                         #   95 percent confidence interval:
                                         #      -1.48  -0.25

Twee regels uit die uitvoer verdienen je aandacht. Eerst het min-teken: R rekent voer één mín voer twee, en omdat voer twee zwaarder uitpakt is het verschil −0.87. Niks engs, gewoon de rekenvolgorde.

En dan de regel waar dit blok om draait: 95 percent confidence interval: −1.48 tot −0.25. Dat is het betrouwbaarheidsinterval / confidence interval: niet één schatting maar een héle reeks — alle waarden voor het echte verschil die goed bij deze dertig potten passen. Eerste kijkvraag, altijd: zit de 0 erin? Kijk zelf: van −1.48 tot −0.25, en nul valt erbuiten. “Geen verschil” past dus níét bij deze data. Maar waar komt zo’n bandje eigenlijk vandaan? Dat kun je vóélen, zonder één formule.

trek het

De bootstrap: trek duizend keer opnieuw

trekken met teruglegging · de middelste 95% als interval

In andere blokken schudden we: band doorknippen, kijken of het verschil wegvalt. Vandaag maken we een andere beweging — we trekken opnieuw. Het probleem is namelijk: je hebt maar één steekproef, en je zou willen weten wat er gebeurd was met een ándere. Nieuwe potten kweken duurt weken. Maar er is een brutale sluiproute: doe alsof je eigen data de kas is. Trek uit je tien potten opnieuw tien potten, mét teruglegging — sommige potten kom je dubbel tegen, andere sla je net over. Elke trekking is een wereld die nét even anders liep, met een nét even ander verschil. Doe dat duizend keer. Dit heet de bootstrap / bootstrap.

Diezelfde teruglegging zag je in Oefening 4.1, waar je uit een ton met 342 pinguïns graaide. Het verschil zit niet in de handeling maar in waar je je arm in steekt: daar was er een populatie en kenden we de waarheid, hier heb je alleen je eigen tien potten en doe je alsóf die de kas zijn.

In R · trek duizend keer opnieuw

b1 <- voer$gewicht[voer$groep == "behandeling1"]
b2 <- voer$gewicht[voer$groep == "behandeling2"]

# één keer opnieuw trekken: tien potten mét teruglegging, per groep
mean(sample(b1, replace = TRUE)) - mean(sample(b2, replace = TRUE))

# doe het 1000 keer en kijk waar de verschillen landen
set.seed(42)
werelden <- replicate(1000,
  mean(sample(b1, replace = TRUE)) - mean(sample(b2, replace = TRUE)))
hist(werelden)                            # een bergje rond -0.87
abline(v = 0, col = "red", lwd = 2)       # en dáár staat de nul
quantile(werelden, c(.025, .975))         # → -1.39 en -0.29

▶ Zie het dansen: het trek-speeltje

Kijk naar dat plaatje. Duizend werelden die dansen rond je schatting — niet rond nul, zoals bij het schudden, maar rond −0.87, want elke wereld is uit jóúw data getrokken en draagt het verschil in zich mee. De middelste 95% van die dans loopt van −1.39 tot −0.29. Herken je ’m? Bijna: t.test gaf je −1.48 tot −0.25. Het betrouwbaarheidsinterval ís de dans, en de formule is de snelweg ernaartoe.

Bijna, en dat “bijna” is de les. Leg de twee naast elkaar: de formule begint 0.09 gram lager en eindigt 0.04 gram hoger. Aan allebei de kanten plakt hij er iets bij. Pech met jouw duizend trekkingen is het niet — begin met een ander toevalsgetal en de dans blijft de smalste van de twee.

Waar dat bijgeplakte stuk vandaan komt: t.test rekent een toeslag omdat je de spreiding niet ként maar schat, en bij tien potten per groep is die fors. De dans betaalt hem niet; die leest gewoon af waar de middelste 95% van jouw werelden ligt. Smaller is hier dus niet zorgvuldiger — eerder wat brutaal. Bij Z) zie je wat die brutaliteit kost.

En de nul? Die staat als rode streep búíten het feest: 999 van de 1000 werelden landen onder nul, één verdwaalde er net overheen. Zit nul buiten de dans, dan is het verschil er echt. Zo lees je voortaan élk interval: eerst kijken waar het ligt, dan pas of nul nog mee mag doen.

Drie keer duizend werelden, drie verschillende vragen

Dit is de derde keer dat je duizend werelden bouwt. In code lijken ze verdacht veel op elkaar, maar ze beantwoorden verschillende vragen — en het verschil zit niet in de handeling.

	wat je met de data doet	de werelden dansen rond	de vraag die je stelt
schudden (Oefening 3.1)	de twee kolommen van elkaar losknippen	de nul	kan toeval dit namaken?
grabbelen (Oefening 4.1)	grijpen uit de populatie, mét terugleggen	de waarheid (4202 gram)	hoe hard wiebelt een gemiddelde van steekproef tot steekproef?
bootstrappen (hier)	grijpen uit je eigen steekproef, mét terugleggen	jouw schatting (−0.87)	hoe zeker is dit ene getal?

Grabbelen en bootstrappen zijn in code niet uit elkaar te houden — allebei sample(..., replace = TRUE). Het verschil is waar je je arm in steekt: in 4.1 in een ton waarvan we de waarheid kennen, een luxe die alleen in een oefening bestaat; hier in de tien potten die de kweker écht heeft, juist omdát de kas onbereikbaar is.

Geeft die dans ook een p? Nee. Jouw werelden dansen rond −0.87, jouw eigen antwoord — dat zag je net in het histogram. Een p vraagt iets anders — hoe ver zou het toeval komen als er niks was? — en dat is een dans rond nul. Die heb je hier niet gebouwd. De bootstrap levert je dus een interval, geen toets. Geen gemis: de nul-vraag heb je al beantwoord, want het interval loopt van −1.39 tot −0.29 en de nul valt erbuiten.

Let wel goed op wát je met die werelden telt, want twee dingen lijken op elkaar en zijn het niet. Kijken of de nul buiten de middelste 95% valt: dat mag — dat deed je hierboven, en dat is het interval aflezen. De staartfractie een p noemen: dat mag niet. Want dat is de verleiding: tel hoeveel werelden aan de andere kant van nul landen, en noem dat getal je p. Bij ons was dat er één van de duizend, dus dat leest als p = .001 — terwijl t.test op dezelfde data p = .009 gaf. Bijna tien keer zo veel, en dat is geen toeval: die telling meet hoe ver jóúw dans over de nul heen reikt, niet hoe vaak het toeval jouw uitkomst zou halen. Was er géén enkele wereld overheen gegaan, dan had de telling “kleiner dan .001” geroepen — een getal dat alleen verklapt hoeveel werelden je hebt gebouwd.

In Oefening 6.1 tel je wél werelden om een p te krijgen. Dat zijn geschudde werelden, en die dansen rond nul. Dat is het hele verschil.

snap het

Wat het interval wél en niet belooft

de 95% gaat over de methode · breedte en precisie

Wat zegt dat interval nou eigenlijk? Dit: het echte verschil tussen de twee voeren — dat van álle mogelijke potten, niet alleen deze dertig — kan van alles zijn tussen een kwart gram en bijna anderhalve gram. Dat is eerlijk: het verschil is echt (nul doet niet mee), maar over hoe gróót het is, zijn deze data nog behoorlijk ruimhartig. En zie je meteen waarom het interval rijker is dan de toets: t = −3.01, p = .009 zegt alleen “niet nul”. Het interval zegt “niet nul, én ergens tussen klein en fors”. Daarom kijken we schattend eerst naar het interval; de toets is de bijvraag die er al in zit.

Nu de zin die íédereen opschrijft en die tóch fout is: “er is 95% kans dat het echte verschil tussen −1.48 en −0.25 ligt.” Klinkt redelijk. Klopt niet. Het echte verschil is geen kansspel — het ligt al ergens vast, wij weten alleen niet waar. En dít interval staat ook al vast: het vangt het echte verschil, of het mist het. Daar valt niks meer aan te verloten. De 95% gaat over de methode: doe het hele ritueel — steekproef, interval — honderd keer, en ongeveer 95 van die honderd intervallen vangen de echte waarde. Denk aan een visser met een net dat 95 van de 100 worpen raak is. Over déze ene worp kan hij maar één ding zeggen: het net komt uit een goede familie. Jouw bootstrap liet precies dat zien — élke nieuwe trekking een iets ander interval, en het vertrouwen zit in de hele dans, niet in één wereld.

Nog één blik, op de breedte. Breed interval = onzekere schatting; smal interval = precieze. Dit interval is aan de brede kant — logisch, tien potten per groep is een handjevol. Wil je het smaller? Meer potten. De vraag “hoe groot is het effect, en hoe precies weten we dat?” is bijna altijd interessanter dan “significant, ja of nee” — en het interval beantwoordt ze allebei in één regel.

jouw beurt

De controle-potten doen weer mee

t-toets · bootstrap · een interval dat de nul schampt

De controle-potten mogen eindelijk van de zijlijn af. Voer twee won het van voer één, maar de echte kwekersvraag is natuurlijk: doet voer twee het ook beter dan gewoon water? Zet eerst de twee groepen klaar, en beantwoord daarna de vragen. De laatste heet Z): de aannames staan altijd achteraan, en die letter is met opzet geen volgletter.

In R · de twee groepen klaarzetten

pg <- read.csv("plantgroei.csv")
ez <- pg[pg$groep %in% c("controle", "behandeling2"), ]

ct <- ez$gewicht[ez$groep == "controle"]
v2 <- ez$gewicht[ez$groep == "behandeling2"]

nrow(ez) hoort 20 te geven — tien potten water, tien potten voer twee. Staat er 30, dan is je filter niet aangekomen.

Jouw beurt

a) Hoeveel droge plant leverde elke groep gemiddeld op, en hoe groot is het verschil? (rapporteer beide M’s en het verschil, in grammen)

b) Laat t.test(gewicht ~ groep, data = ez) het bandje erbij geven en lees het af. (rapporteer het 95% CI met allebei de grenzen)

c) Zit de nul in dat interval? Zeg wat dat betekent voor de kweker. (over planten, niet over getallen)

d) In de romp van dit blok stond een min-teken voor het verschil en hier niet. Wat is er tussen die twee vergelijkingen veranderd? (één zin over de rekenvolgorde)

e) Maak de dans nu zelf: trek per groep duizend keer opnieuw, mét teruglegging, en reken elke wereld als voer twee mín controle. Zet eerst set.seed(42). Waar ligt de middelste 95%? (rapporteer de twee grenzen)

f) Hoeveel van je duizend werelden landen boven nul? Dat getal ziet eruit als een p, en dat is het niet. (het aantal, en één zin over wat een p wél telt)

g) Leg de twee intervallen naast elkaar — dat van t.test en dat van je bootstrap. Hoe dicht zitten ze op elkaar, en welke van de twee is het smalst? (vier grenzen en één zin)

h) Vergelijk de breedte van het interval dat t.test je gaf met die van het voer-één-tegen-voer-twee-interval hierboven. Allebei tien potten per groep, en tóch is de ene breder. Waar komt dat verschil vandaan? (twee breedtes in grammen, en wat je in de spreiding per groep ziet)

i) Wat betekent die 95% precies? (over de methode, niet over dit ene interval)

j) Schrijf het advies aan de kweker, in maximaal vijf zinnen: “voer twee werkt beter dan water”, of “kweek eerst nog een rij potten”? (een advies waarin de onzekerheid zelf ook voorkomt)

Z) Mag je alles hierboven geloven? Waar leunt dat bandje op — kijk naar de vorm van de twintig gewichten (hist() per groep, of boxplot(gewicht ~ groep, data = ez)), naar de twee spreidingen (tapply(ez$gewicht, ez$groep, sd)), en naar de vraag of elke pot maar één keer meetelt.

Antwoord — open na je eigen poging

a) Gewoon water M = 5.03 gram, voer twee M = 5.53. Verschil: 0.49 gram in het voordeel van voer twee.

b) 95% CI [0.01, 0.98]. Het echte verschil ligt volgens deze data ergens tussen een haartje en bijna een gram.

c) Nét niet. De ondergrens is 0.01 — de nul valt er met een millimeter buiten. Voor de kweker: op deze twintig potten past “voer twee doet niks meer dan water” nog maar net níét bij wat hij zag. Dat is iets anders dan bewijs dat het werkt.

d) In de romp stond voer één tegenover voer twee, en t.test zet de groepsnamen op alfabet: behandeling1 mín behandeling2 is negatief. Hier is het behandeling2 mín controle, en dat is positief. Er is niets aan de planten veranderd, alleen aan de volgorde waarin R aftrekt.

e) De middelste 95% van de duizend werelden loopt van 0.06 tot 0.90.

f) 985 van de 1000. Dat is een controlegetal — zit je daar ver vandaan, dan is er iets mis met je code. Het is géén p: jouw werelden dansen rond 0.49, jouw eigen antwoord, en een p telt hoe ver het toeval komt in werelden die rond nul dansen. Die heb je hier niet gebouwd.

g) t.test gaf [0.01, 0.98], de bootstrap [0.06, 0.90] — dezelfde streek van de getallenlijn, allebei met de nul er nét buiten. De formule is de snelweg, de bootstrap de landweg. Kijk nu naar de uiteinden: de formule begint lager én eindigt hoger, dus zijn bandje is aan allebei de kanten ruimer. Net als in de romp is de dans de smalste, en om dezelfde reden — de formule rekent een toeslag voor de geschatte spreiding, de dans niet. Speling is dat niet: begin met een ander toevalsgetal en de dans blijft smaller. Minder aannemen over de vórm van je data maakt de dans dus niet vanzelf nauwkeuriger; bij Z) zie je wanneer die soberheid wél in je voordeel werkt.

h) Voer twee tegen water is 0.98 gram breed, voer één tegen voer twee 1.23 gram. (Trek je de twee grenzen af zoals ze bij b) staan, dan kom je bij de eerste op 0.97 uit. Dáár is het wél afronding, en die is van dit boekje: bij b) staan de grenzen op twee decimalen, terwijl de breedte uit de hele getallen komt — 0.98287 − 0.00513 = 0.97774, en dat is 0.98.) Toch dezelfde tien potten per groep. Het verschil zit in de spreiding bínnen de groepen: voer twee is de rustigste bak van de drie (SD = 0.44), water zit op SD = 0.58, en voer één is de wildste met SD = 0.79. Een breed interval komt dus niet alleen van weinig potten, maar ook van potten die onderling veel verschillen.

i) Niet: “er is 95% kans dat het echte verschil hierin ligt.” Dit interval staat al vast, en het echte verschil ook — het vangt hem, of het mist hem. De 95% gaat over de methode: doe het hele ritueel honderd keer, en ongeveer 95 van die honderd intervallen vangen de echte waarde. Jouw bootstrap liet dat zien: elke trekking een iets ander bandje.

j) Bijvoorbeeld:

Voer twee leverde in dit experiment gemiddeld 0.49 gram meer droge plant op dan water, 95% CI [0.01, 0.98]. De nul valt er net buiten, dus “water doet net zoveel” past nog maar nauwelijks bij deze twintig potten. Maar de ondergrens ligt zó dicht bij nul dat het voordeel ook vrijwel niets kan zijn. Voor een kweker die hierop zijn inkoop baseert is dat te wankel. Een tweede rij potten kost minder dan een jaar lang het verkeerde voer.

Wat erin moet staan: dát er een verschil gemeten is, én dat de grootte ervan wankel staat. Een advies dat maar één van die twee noemt, is de helft.

Z) Ja, hier mag het. Drie dingen om te bekijken, en alle drie vallen ze goed uit. De twintig gewichten liggen redelijk symmetrisch — de scheefheid is 0.23 bij water en 0.48 bij voer twee, en met tien potten per groep is dat gewoon ruis, geen scheve verdeling. De twee spreidingen liggen niet ver uiteen: 0.58 tegen 0.44, een verhouding van 1.32, en het huis houdt onder de 2 aan. En elke pot is één keer gewogen, dus de waarnemingen zijn onafhankelijk.

En als het wél was misgegaan, dan is dit de weg terug:

Blijft staan: de twee gemiddelden en het verschil van 0.49 gram bij a). Dat is optellen en delen; daar verandert geen aanname iets aan.
Doe je over: het interval bij b) en daarmee de conclusie bij c) en j). Die leunen op de aanname dat de gemiddelden netjes rond de waarheid wiebelen.

En het mooie is dat je die weg terug al gelopen hebt: je bootstrap-interval bij e) leunt daar veel minder op. Waren de gewichten flink scheef geweest, dan was [0.06, 0.90] het eerlijkere antwoord en [0.01, 0.98] het verdachte. Dat de twee hier op dezelfde plek liggen, is zelf een geruststelling. Dat de dans smaller uitvalt is dat niet: dat komt van de toeslag die de formule rekent, niet van scheefheid in jouw potten.

En die toeslag niet betalen kost iets — dat is wat de romp met “brutaal” bedoelde. Wil je het zien: kweek duizenden nagemaakte kassen van tien potten per groep waarvan jij het echte verschil kent, en tel bij elke kas of het bandje dat verschil vangt. De formule zit dan op zijn beloofde 95 van de 100; de dans haalt er drie tot vier minder. Vandaar dat je bij een echte analyse geen kale percentielen neemt maar een bijgestelde bootstrap — het BCa-interval, dat corrigeert voor precies deze twee dingen: de scheefheid en de te krappe band.

rapportage

Schrijf het op als een onderzoeker

CI tussen blokhaken · leesrichting van het verschil

Rapportage

Je kent het boodschappenlijstje-in-zinsvorm al: wie, wat, waar lag het. Vanaf dit hoofdstuk komt er één ding bij — het bandje, tussen blokhaken. Kijk maar:

Om te onderzoeken of twee soorten plantenvoer verschillend uitpakten, kweekten we twintig plantjes — tien per voer — en wogen we ze gedroogd. De potten met voer twee leverden gemiddeld 0.87 gram meer droge plant op dan die met voer één (M = 5.53 om M = 4.66 gram), 95% CI van het verschil [0.25, 1.48].

Nog steeds één adem, en de zin staat op eigen benen. Twee dingen vallen op. Eén: het interval schrijf je als CI met de twee grenzen tussen blokhaken, klein naar groot — en CI blijft rechtop. Schuin gaat wat zélf een waarde draagt (M = 5.53, en straks t en p); CI draagt er geen, die benoemt alleen wat er tussen de haken staat. Twee: R gaf ons −1.48 tot −0.25, en hier staat 0.25 tot 1.48 — we vertellen het verschil in de leesrichting (“voer twee leverde méér op”), en dan draait het bandje gewoon mee. Het min-teken was toch al alleen rekenvolgorde. En zie je wat de zin stiekem óók al vertelt? Nul zit niet tussen die blokhaken — de lezer kan de bijvraag zelf beantwoorden. De t en de p netjes leren opschrijven bewaren we voor de opdracht over de p-waarde.

Jouw beurt: schrijf zo’n zelfstandige zin voor de taak hierboven — voer twee tegenover gewoon water (verschil 0.49 gram, M = 5.53 om M = 5.03, 95% CI [0.01, 0.98]). En durf je ’m aan, met die ondergrens van 0.01? Hardop voorlezen mag.

◠ Zelfde vorm in het boek

Hier: één verschil van 0.87 gram, en een bandje eromheen dat vertelt hoe stevig dat getal staat.

In het boek: W6 — Bestaat het verschil echt?

Daar komt het interval vóór de toets — eerst hoe groot en hoe precies, en dan pas de bijvraag of nul nog mee mag doen.

Voetnoten

PlantGrowth — Dobson (1983), An Introduction to Statistical Modelling; het klassieke groei-experiment dat met R wordt meegeleverd (datasets::PlantGrowth). Alle getallen in dit blok zijn op die échte data nagerekend.↩︎