4. Logistische regressieanalyse (LRA)

Al dan niet naar jezelf turen in de spiegel

De pad zat aan zijn vijver. Het was een wateroppervlak met geheugen — een vijver die soms iets wist over wie erin keek dat de kijker zelf niet wist.

Sommige dieren bleven staan en tuurden naar zichzelf in het water. Andere liepen door zonder een blik over hun schouder. Het ene dier checkte zijn snorharen alsof er een afspraak op het spel stond. Het andere had het te druk met de geur van de wind, of was te ontspannen om zich druk te maken om zijn eigen kop.

“Twee mogelijkheden per dier,” dacht de pad. “Wel turen, of niet turen. Geen tussenstand.”

Hij vroeg zich af of er iets te voorspellen viel. Misschien dat dieren die net een afwijzing hadden gehad het water opzochten om te checken of er iets aan ze mankeerde. Misschien dat oudere dieren juist minder zelf-bekijkers waren — minder uit te zoeken, of minder bereid het uit te zoeken.

Eén ding wist hij zeker: een gewone rechte regressielijn ging hem niet helpen. Zijn uitkomst was geen continue maat — geen hoeveel, alleen wel of niet. Hij had een ander gereedschap nodig.

De techniek heet logistic regression LRA — en in dit hoofdstuk gaan we hem stap voor stap opbouwen, te beginnen met de allereenvoudigste vraag: als ik niets over een dier weet, hoe goed kan ik dan al gokken of hij tuurt?

Wat je na dit thema moet kunnen

Herkennen dat een binaire $Y$ om LRA vraagt.
Met de hand van $b$ naar kans, odds en odds ratio rekenen, en het kantelpunt $-b_0/b_1$ vinden.
Met de hand het $95\%$ -interval van $b$ en van de odds ratio uitrekenen.
Per voorspeller de Wald-toets lezen; de likelihood-ratio-toets voor het hele model of tussen twee geneste modellen.
$R^2_L$ uitrekenen uit twee deviances.
Een model in blokken opbouwen (hiërarchische LRA): per blok de $\chi^2$ -toets en per model $R^2_L$ .
Uit een classificatietabel het percentage goed, de sensitiviteit, de specificiteit, PPV en NPV halen.
Weten welke aannames je checkt en welke juist niet.
De uitkomst in een APA-zin opschrijven.

Wat dichtgeklapt staat onder Voor wie verder wil hoort hier niet bij. Wat blijft liggen, onderaan, noemt wat niet behandeld wordt.

Twee formules komen in dit thema steeds terug: de logistische vergelijking en die van $R^2_L$ .

4.0 Project- en datavoorbereiding

Aantekeningen netjes, paden kloppen

Download de bron van dit oefenboek (de zip op de voorpagina), pak hem uit en dubbelklik op mvda_werkboek.Rproj — één project voor alle thema’s. Open daarna 04_logistische_regressie/04_logistische_regressie.qmd en draai de code vanuit dat bestand: RStudio voert een chunk uit in de map van het hoofdstuk, dus de paden kloppen. De data van dit thema ligt in 04_logistische_regressie/data/, de helper-functies in common/R/. Niets te installeren als de twee pakketten (car, tidyverse) op je computer staan; anders één keer:

install.packages(c("car", "tidyverse"))

Open de hoofd-dataset.

load("data/pad_spiegel.RData")
str(pad_spiegel)

'data.frame':   152 obs. of  4 variables:
 $ bekeek_zichzelf   : int  0 0 1 0 0 0 1 0 1 1 ...
 $ leeftijd          : num  7.7 2.9 2.7 7 8.9 8 7.2 7.7 7.6 4.9 ...
 $ aantal_afwijzingen: int  2 6 6 0 4 3 3 1 2 3 ...
 $ diersoort         : Factor w/ 3 levels "ree","vos","das": 3 3 1 3 1 1 2 3 2 2 ...

Per dier vier waarden: bekeek_zichzelf (0 = nee, 1 = ja), leeftijd (jaren), aantal_afwijzingen (recente sociale flops) en diersoort (ree/vos/das). De drie voorspellers vormen samen het kandidaat-model.

4.1 De pad telt — eerst met één voorspeller

Van een eenvoudige gok naar een model dat luistert naar één gegeven

“Ik begin klein,” zei de pad. “Eerst tel ik alleen wat ik zie. Dan voeg ik één ding toe. Daarna kijk ik of het beter wordt.”

In deze vijf mini-secties (4.1.a t/m 4.1.e) bouwen we het allereenvoudigste logistische model — één voorspeller, niets anders — en oefenen we kans → odds → odds ratio. Géén Wald-toets en géén likelihood-ratio nog; het woord log-odds komt alleen even langs. Eerst de drie rekenstappen. De rest volgt in 4.1.f en verder.

4.1.a Null-model — hoe goed gokt de pad zonder verdere info?

De pad heeft $152$ dieren in zijn boekje. Voordat hij naar voorspellers kijkt, vraagt hij: als ik niets weet over een specifiek dier, wat is dan mijn beste gok?

load("data/pad_spiegel.RData")

# Hoeveel dieren in totaal? Hoeveel daarvan turen?
N       <- nrow(pad_spiegel)
n_tuurt <- sum(pad_spiegel$bekeek_zichzelf == 1)
p_tuurt <- n_tuurt / N

c(N = N, n_tuurt = n_tuurt, p_tuurt = round(p_tuurt, 3))

      N n_tuurt p_tuurt 
 152.00   73.00    0.48

# Dezelfde telling als één zin: neem de data, en vat hem samen. Het pijltje
# |> geeft wat links staat door aan wat rechts komt. summarise() perst de
# hele tabel samen tot één regel; n() telt de rijen. Een gemiddelde van
# WAAR/ONWAAR is het aandeel WAAR, dus mean(... == 1) is de proportie turers.
pad_spiegel |>
  summarise(N       = n(),
            n_tuurt = sum(bekeek_zichzelf == 1),
            p_tuurt = round(mean(bekeek_zichzelf == 1), 3))

    N n_tuurt p_tuurt
1 152      73    0.48

Zo’n model zonder voorspellers heet het null-model. Het geeft elk dier dezelfde kans, wat je verder ook over dat dier weet. Zijn beste gok is dus ook voor elk dier dezelfde: de uitkomst die het vaakst voorkomt. Elk model dat hierna komt, leg je naast dit null-model: brengt wat je over een dier weet de gok dichter bij wat er echt gebeurt?

Vraag 4.1.a — Pen-en-papier

a) Hoeveel dieren noteerde de pad in totaal? Hoeveel daarvan tuurden er?

b) Bereken handmatig $P(\text{tuurt}) = n_{\text{tuurt}} / N$ . Op twee decimalen.

c) De pad ziet een nieuw dier langs lopen. Hij weet niets over dit dier. Wat is zijn beste gok: tuurt of tuurt niet? Met welke kans gokt hij goed?

d) Bij die strategie — altijd hetzelfde gokken voor elk dier — hoeveel keer per $100$ dieren zou hij naar schatting fout zitten?

Antwoord 4.1.a — open na je eigen poging

a) $N = 152$ dieren in totaal, waarvan $n_{\text{tuurt}} = 73$ turen en $n_{\text{niet}} = 79$ niet.

b) $P(\text{tuurt}) = 73 / 152 \approx .48$ . Bijna de helft tuurt.

c) Omdat niet-tuurt iets vaker voorkomt dan tuurt ( $79$ versus $73$ ), is de beste gok zonder verdere info “tuurt niet”. Hij heeft dan kans $79/152 \approx .52$ om goed te zitten. (Anders zou hij omgekeerd $73/152 \approx .48$ goed zitten — niet veel verschil, hier is het null-model bijna een muntje gooien.)

d) Met de “altijd niet-tuurt”-strategie zit hij ongeveer $48$ keer per $100$ dieren fout. Dat is veel — bijna de helft. De vraag is dus: kan hij beter doen als hij iets over de dieren weet?

4.1.b Een voorspeller toevoegen — `aantal_afwijzingen`

De pad heeft een hypothese: misschien turen dieren die net afgewezen zijn vaker in het water — om te checken of er iets aan ze mankeert.

Hij bouwt een logistisch model met één voorspeller: aantal_afwijzingen (interval, aantal recente sociale flops).

m_afw <- glm(bekeek_zichzelf ~ aantal_afwijzingen,
             data   = pad_spiegel,
             family = binomial(link = "logit"))

summary(m_afw)$coefficients

                     Estimate Std. Error   z value   Pr(>|z|)
(Intercept)        -0.5429478  0.3142876 -1.727551 0.08406876
aantal_afwijzingen  0.1812410  0.1050218  1.725747 0.08439298

Het getal $b_1$ zelf is het effect op de log-odds, en daar heeft bijna niemand gevoel voor. Het teken lees je wel. Een positieve $b_1$ betekent: hoe meer afwijzingen, hoe groter de kans dat een dier tuurt. Een negatieve $b_1$ : hoe meer afwijzingen, hoe kleiner die kans. Hoe gróót het effect van afwijzingen op turen is, lees je straks af aan de odds ratio.

Vraag 4.1.b — Pen-en-papier

a) Lees uit de output de intercept $b_0$ en de coëfficiënt $b_1$ voor aantal_afwijzingen. Noteer ze op drie decimalen.

b) Wat is het teken van $b_1$ ? Wat zegt dat teken over de relatie tussen aantal afwijzingen en de kans op turen?

c) Bedenk in één zin of de hypothese van de pad door dit teken bevestigd of weerlegd wordt.

Antwoord 4.1.b — open na je eigen poging

a) $b_0 = -0.543$ (intercept), $b_1 = 0.181$ (helling voor aantal_afwijzingen).

b) $b_1$ is positief: hoe meer afwijzingen, hoe groter de kans dat een dier in het water staart.

c) De hypothese van de pad wordt bevestigd door het teken: afgewezen dieren zoeken vaker hun spiegelbeeld op. Of het effect ook significant van nul afwijkt is een vraag voor straks (4.1.f).

4.1.c Van $b$ naar kans — de sigmoid-formule

Het teken alleen voorspelt geen concrete kans. Voor “wat is de kans dat een dier met $5$ afwijzingen tuurt?” heb je de logistische regressie-vergelijking nodig:

$P(Y=1) \;=\; \frac{e^{b_0 + b_1 X}}{1 + e^{b_0 + b_1 X}}$

Hier is $e \approx 2.718$ — Eulers getal, zoals $\pi$ een vaste constante. Op een rekenmachine staat de knop als eˣ of exp; in R is dat exp(). Wat er in de exponent ook uitkomt, groot of klein, positief of negatief: de formule maakt er altijd een kans tussen $0$ en $1$ van.

Reken-truc — de gek-ding-vorm

Voor snelheid op je rekenmachine: noem de exponent-uitdrukking even $E$ (hoofdletter; let op — niet hetzelfde als $e$ , kleine letter):

$P = \frac{E}{1 + E} \quad \text{waar} \quad E = e^{b_0 + b_1 X} = \exp(b_0 + b_1 X)$

Eerst reken je $E$ uit (één keer typen). Dan gek-ding gedeeld door $1$ plus gek-ding. Klaar.

In R kun je het ook automatisch via predict():

# Bereken kansen voor een paar X-waarden via R:
nieuw <- data.frame(aantal_afwijzingen = c(0, 2, 3, 5))
nieuw$P_tuurt <- predict(m_afw, newdata = nieuw, type = "response")
round(nieuw, 3)

  aantal_afwijzingen P_tuurt
1                  0   0.368
2                  2   0.455
3                  3   0.500
4                  5   0.590

Vraag 4.1.c — Pen-en-papier

Met $b_0 = -0.543$ en $b_1 = 0.181$ :

a) Reken handmatig $P(\text{tuurt})$ voor een dier met $\text{aantal\_afwijzingen} = 0$ . Toon ook tussenstap $E$ .

b) Idem voor $\text{aantal\_afwijzingen} = 3$ . Wat valt je op aan de uitkomst?

c) Idem voor $\text{aantal\_afwijzingen} = 5$ .

d) Zonder uit te rekenen: bij welke $X$ -waarde ligt $P(\text{tuurt})$ precies op $.5$ ? (Hint: gebruik de drempel-truc $X^* = -b_0/b_1$ uit de mini-callout hieronder.)

De drempel-truc — bij welk

X

P = .5

Met het gek-ding van hierboven: $P = E/(1 + E)$ is precies $.5$ als $E = 1$ , want $1/(1 + 1) = .5$ . En $E = e^{b_0 + b_1 X}$ is precies $1$ als de exponent nul is. Dus: zet $b_0 + b_1 X = 0$ , los op naar $X$ :

$X^* = -\dfrac{b_0}{b_1}$

Voor ons: $X^* = -(-0.543)/0.181 = 3.00$ . Bij precies $3$ afwijzingen kantelt het model van “voorspelt niet-tuurt” naar “voorspelt tuurt”.

Antwoord 4.1.c — open na je eigen poging

a) $E = e^{-0.543 + 0.181 \cdot 0} = e^{-0.543} \approx 0.581$ . Dan $P = 0.581 / (1 + 0.581) \approx 0.367$ . Een dier zonder afwijzingen heeft ~ $37\%$ kans om te turen.

b) $E = e^{-0.543 + 0.181 \cdot 3} = e^{0} = 1$ . Dan $P = 1 / (1 + 1) = 0.50$ . Bij precies $3$ afwijzingen voorspelt het model fifty-fifty. Dat is het kantelpunt.

c) $E = e^{-0.543 + 0.181 \cdot 5} = e^{0.362} \approx 1.436$ . Dan $P = 1.436 / 2.436 \approx 0.589$ . Een dier met $5$ afwijzingen heeft ~ $59\%$ kans om te turen.

d) Bij $X^* = -b_0/b_1 = 0.543/0.181 = 3.0$ afwijzingen. Klopt met antwoord (b).

Inzicht. Het model maakt van aantal_afwijzingen een continue kans tussen $0$ en $1$ . Bij weinig afwijzingen voorspelt het “niet-tuurt”, bij veel afwijzingen “tuurt”, en bij precies $3$ afwijzingen weet het model het niet — die waarde is de drempel.

4.1.d Odds — een eerste laag ratio

Een odds is een verhouding: hoe vaak verwachten we $Y = 1$ ten opzichte van $Y = 0$ ? In formule:

$\text{odds}(Y = 1) \;=\; \frac{P(Y = 1)}{1 - P(Y = 1)} \;=\; \frac{P(Y = 1)}{P(Y = 0)}$

Zie wat hier gebeurt: de odds is al een ratio. Eén dier-situatie, twee kansen tegen elkaar afgezet. Eerste laag: één ratio.

Eén meeuw — odds is al een ratio

Stel een zilvermeeuw vangt vis bij $70\%$ van zijn duik-pogingen. Op een dag van $100$ pogingen:

	Vis gevangen	Niet gevangen	Totaal
Eén dag duiken: de vangst van de zilvermeeuw
Zilvermeeuw	70	30	100

Kans op vis: $P_{\text{vis}} = .70$
Kans op niet-vis: $P_{\text{geen}} = .30$
Odds op vis $= .70 / .30 \approx 2.33$ — “voor elke keer dat hij faalt, slaagt hij $2.33$ keer” of “ $7$ op $3$ ”

Eén meeuw, twee mogelijkheden, één verhouding. Eerste laag.

Odds en kans zeggen hetzelfde, op twee schalen. Een kans van $.5$ is odds $1$ : even vaak wel als niet. Een kans boven $.5$ geeft odds boven $1$ , een kans onder $.5$ odds onder $1$ . Alleen de bereiken verschillen: een kans zit opgesloten tussen $0$ en $1$ , odds lopen van $0$ tot oneindig.

Voor onze pad — terug naar aantal_afwijzingen. Uit de kansen die je in 4.1.c hebt uitgerekend, volgen direct de odds:

# Bouw verder op nieuw uit 4.1.c:
nieuw$odds_tuurt <- nieuw$P_tuurt / (1 - nieuw$P_tuurt)
round(nieuw, 3)

  aantal_afwijzingen P_tuurt odds_tuurt
1                  0   0.368      0.581
2                  2   0.455      0.835
3                  3   0.500      1.001
4                  5   0.590      1.438

# mutate() zet er een kolom bij en geeft de hele tabel terug. Vandaar de
# toewijzing nieuw <- ervoor: zonder die regel reken je het uit en gooi je
# het meteen weer weg.
nieuw <- nieuw |>
  mutate(odds_tuurt = P_tuurt / (1 - P_tuurt))
round(nieuw, 3)

  aantal_afwijzingen P_tuurt odds_tuurt
1                  0   0.368      0.581
2                  2   0.455      0.835
3                  3   0.500      1.001
4                  5   0.590      1.438

Nog één woord voordat je rekent. Wat in de gek-ding-vorm in de exponent stond, $b_0 + b_1 X$ , heet de log-odds of logit: de logaritme van de odds. Een raar woord, maar laat je er niet door blokkeren. Denk gewoon: hoe hoger de log-odds, hoe meer kans. Het kantelpunt ligt bij nul: log-odds $0$ hoort bij odds $1$ en dus bij kans $.5$ . Negatieve log-odds horen bij een kans onder $.5$ , positieve bij een kans erboven.

Het gek-ding $E$ uit 4.1.c ís dus de odds. Bij nul afwijzingen: $E = e^{-0.543} \approx 0.581$ , hetzelfde getal als in de odds-kolom hierboven. Terug gaat met de logaritme: $\ln(0.581) \approx -0.543$ . Gebruik daarvoor de knop ln op je rekenmachine, niet log (die rekent met grondtal $10$ ). In R en in de formules van dit thema heet hij gewoon $\log$ , en dan bedoelen we steeds deze natuurlijke logaritme.

Vraag 4.1.d — Pen-en-papier

a) Reken uit de kansen die je in 4.1.c hebt gevonden de odds voor elk niveau van aantal_afwijzingen ( $0$ , $3$ , $5$ ).

b) Bij aantal_afwijzingen = 3 heb je $P = .50$ . Wat zijn de odds dan? Wat valt je op aan dat getal?

c) Snelle weg: omdat $e^{\ldots}$ de logaritme ongedaan maakt, is $\text{odds} = e^{b_0 + b_1 X}$ . Reken dat uit voor dezelfde drie waarden. Klopt het met je antwoorden bij (a)?

Antwoord 4.1.d — open na je eigen poging

a) Met de drie kansen $.367, .500, .589$ :

$X$	$P$	odds = $P/(1-P)$
$0$	$.367$	$.367/.633 \approx 0.580$
$3$	$.500$	$.500/.500 = 1.000$
$5$	$.589$	$.589/.411 \approx 1.433$

b) Bij $P = .50$ is odds $= 1$ . Even waarschijnlijk wel als niet — fifty-fifty in odds-taal heet één-op-één.

c) Snelle weg via de logit: als $\log(\text{odds}) = b_0 + b_1 X$ , dan $\text{odds} = e^{b_0 + b_1 X}$ . Voor $X = 0$ : $e^{-0.543} \approx 0.581$ . Voor $X = 3$ : $e^{0} = 1$ . Voor $X = 5$ : $e^{0.362} \approx 1.436$ . Dat zijn dezelfde odds als bij (a), op de derde decimaal na: daar liep de route via de kans, en die was al afgerond ( $0.580$ tegenover $0.581$ ). Rond pas aan het eind af, dan kom je op beide routes op hetzelfde getal uit.

4.1.e Odds ratio — een ratio van ratios

Als odds een ratio is binnen één situatie, dan is de odds ratio (OR) een ratio tussen twee situaties — bv. tussen twee waarden van $X$ .

$\mathit{OR} \;=\; \frac{\text{odds}_{X = k+1}}{\text{odds}_{X = k}} \;=\; e^{b_1} \;=\; \exp(b_1)$

Dat geldt voor elke eenheids-stijging in $X$ — de OR is constant. Dat is de essentie van logistische regressie: per eenheid $X$ omhoog, vermenigvuldigen de odds met dezelfde factor.

Twee meeuwen — odds ratio (een ratio van ratios)

Naast de zilvermeeuw vliegt een stormmeeuw, die in $40\%$ van zijn pogingen vis vangt:

	Vis gevangen	Niet gevangen	Totaal
Eén dag duiken: de vangst van twee meeuwen
Zilvermeeuw	70	30	100
Stormmeeuw	40	60	100

Odds zilvermeeuw $= 70/30 \approx 2.33$ ← eerste laag ratio
Odds stormmeeuw $= 40/60 \approx 0.67$ ← eerste laag ratio
Odds ratio $= (70/30) / (40/60) = 3.50$ ← tweede laag ratio

“De zilvermeeuw heeft $3.5$ keer hogere odds op vis dan de stormmeeuw.” De OR is letterlijk een ratio van ratios — eerst per dier odds, dan die odds op elkaar gedeeld.

OR als factor — de drie ankerwoorden

De OR is een factor: een vermenigvuldigingsgetal waarmee je de odds vermenigvuldigt om bij de nieuwe odds te komen, als je $X$ met één eenheid laat toenemen. Drie woorden, één concept.

Factor — als ik “factor tien” tegen jouw leeftijd gooi, en jij bent $21$ , dan ben je opeens $210$ jaar. Geen plus tien, geen min tien — keer tien. Zo werkt OR. Lineair gedraagt zich additief; exponentieel gedraagt zich multiplicatief. Plus / min versus keer / door.

Het teken van $b$ zie je terug in de OR. Bij $b = 0$ is $\mathit{OR} = e^0 = 1$ : keer één, de odds blijven gelijk, geen effect. Een positieve $b$ geeft een OR boven $1$ , een negatieve $b$ een OR onder $1$ . Wil je de OR in procenten zeggen, dan gaat dat over de odds en niet over de kans: $\mathit{OR} = 1.5$ maakt de odds $50\%$ hoger, $\mathit{OR} = 0.8$ maakt ze $20\%$ lager.

Voor onze pad:

# Odds ratio per extra afwijzing — uit het model:
OR <- exp(coef(m_afw)["aantal_afwijzingen"])
round(OR, 3)

aantal_afwijzingen 
             1.199

Vraag 4.1.e — Pen-en-papier

a) Reken $\mathit{OR} = e^{b_1} = e^{0.181}$ handmatig uit (op je rekenmachine eˣ of exp).

b) Interpreteer in één zin: wat zegt deze OR over het effect van één extra afwijzing op de kans dat een dier tuurt?

c) Voor vijf extra afwijzingen: bereken $\mathit{OR}(5) = \mathit{OR}(1)^5$ . Wat zijn dan de odds van een dier met $5$ afwijzingen, ten opzichte van een dier met $0$ afwijzingen?

d) Vergelijk je antwoord bij (c) met de odds-tabel uit 4.1.d: klopt de berekening?

Antwoord 4.1.e — open na je eigen poging

a) $\mathit{OR} = e^{0.181} \approx 1.198$ . (R rekent met de onafgeronde $b_1$ en geeft $1.199$ .)

b) Per één extra afwijzing worden de odds dat het dier tuurt vermenigvuldigd met factor $1.20$ — een toename van ongeveer $20\%$ in de odds. (Niet in de kans direct! De kans schaalt niet-lineair met de odds.)

c) $\mathit{OR}(5) = 1.198^5 \approx 2.468$ . De odds van een dier met $5$ afwijzingen zijn ongeveer $2.5$ keer zo groot als de odds van een dier met $0$ afwijzingen.

d) Via de odds-tabel: $\text{odds}_{X=5} / \text{odds}_{X=0} = 1.436 / 0.581 \approx 2.472$ . Op de derde decimaal na hetzelfde als bij (c): beide routes rekenen met afgeronde tussenstappen. R, zonder iets af te ronden, komt op $2.475$ .

Inzicht. De OR is constant over $X$ — dezelfde factor $1.20$ geldt of je van $0$ naar $1$ gaat, of van $10$ naar $11$ . Op odds-schaal is logistische regressie lineair-in-de-vermenigvuldiging. Dat is wat $e^{b_1}$ in essentie zegt.

Tot hier: drie rekenstappen — kans → odds → odds ratio — met één voorspeller. Hierna komen toets-vragen: wijkt $b_1$ van nul af in de populatie? (Wald-toets, 4.1.f) en doet het model als geheel iets? (LR-toets, 4.1.g). En daarna pas een uitbreiding naar meerdere voorspellers tegelijk.

4.1.f Wald-toets — wijkt $b_1$ in de populatie af van nul?

We hebben in de steekproef $b_1 = 0.181$ gevonden. Maar dat is een schatting. De vraag is: kunnen we ook in de populatie zeggen dat afwijzingen iets te maken hebben met turen? Daarvoor toetsen we de populatie-coëfficiënt $b^*_1$ tegen nul. De ster ( $^*$ ) maakt zichtbaar dat het over de populatie gaat en niet over de $b_1$ van deze $152$ dieren, net als in thema 1. Een Griekse $\beta$ zie je hier dus niet: die letter is in dit boek het gestandaardiseerde gewicht uit thema 1.

H_0

H_1

— populatie-taal, geen significantie

De hypothesen gaan over de populatie, niet over je steekproef en niet over een $p$ -waarde:

$H_0: b^*_1 = 0$ — in de populatie heeft aantal_afwijzingen geen effect op de log-odds van turen.
$H_1: b^*_1 \neq 0$ — in de populatie heeft aantal_afwijzingen wel een effect op de log-odds van turen.

“De $H_0$ is dat het niet significant is” — die zin klopt nooit. Significantie is een eigenschap van de toetsuitkomst, niet van de coëfficiënt of van de hypothese.

De Wald-toets vergelijkt de steekproef-coëfficiënt $b_1$ met zijn standaardfout $SE_{b_1}$ . Twee gedaanten, mathematisch identiek:

$z = \frac{b_1}{SE_{b_1}} \qquad \text{of equivalent} \qquad \mathit{\chi}^2(1) = z^2$

R’s summary() geeft $z$ (geen $df$ ). In SPSS-output en in veel rapportages staat $\chi^2(1)$ . Beide verwijzen naar dezelfde toets — de ene is het kwadraat van de andere — dezelfde $p$ -waarde.

# Coëfficiënten met SE en z-waarde:
round(summary(m_afw)$coefficients, 3)

                   Estimate Std. Error z value Pr(>|z|)
(Intercept)          -0.543      0.314  -1.728    0.084
aantal_afwijzingen    0.181      0.105   1.726    0.084

Vraag 4.1.f — Pen-en-papier

Met $b_1 = 0.181$ en $SE_{b_1} = 0.105$ :

a) Bereken Wald- $z$ handmatig.

b) Bereken $\mathit{\chi}^2(1) = z^2$ .

c) De $p$ -waarde uit de output is $.084$ . Verwerp je $H_0$ op $\alpha = .05$ ?

d) Welke conclusie trek je over de populatie? Formuleer in één zin volgens de regels uit de $H_0/H_1$ -callout hierboven.

e) De pad’s hypothese was dat afgewezen dieren vaker turen. Heeft de Wald-toets die hypothese ondersteund?

Antwoord 4.1.f — open na je eigen poging

a) $z = 0.181 / 0.105 \approx 1.724$ . (Klein-positief — wijst dezelfde kant op als $b_1$ , maar niet zwaar.) De uitvoer van R zegt $1.726$ , omdat R met de onafgeronde $b$ en $SE$ rekent.

b) $\chi^2(1) = 1.724^2 \approx 2.97$ . (Mathematisch identiek aan de $z$ -toets — zelfde toets, andere vorm.) Rond pas aan het eind af: wie eerst $z$ afrondt op $1.72$ , komt op $2.96$ . En R zelf, met zijn onafgeronde $z$ , komt op $2.98$ : dat getal zie je in SPSS-uitvoer staan.

c) $p = .084 > .05$ . We verwerpen $H_0$ niet. Op $\alpha = .05$ kunnen we niet concluderen dat $b^*_1 \neq 0$ in de populatie.

d) In de huidige steekproef ( $N = 152$ ) is het effect van aantal_afwijzingen op de log-odds van turen niet groot genoeg om met voldoende zekerheid te zeggen dat $b^*_1$ in de populatie van nul afwijkt. Korter, niet-significant betekent niet “geen effect” — het betekent “niet genoeg signaal in deze data om te concluderen dat het effect bestaat”.

e) Strikt genomen: de Wald-toets geeft onvoldoende bewijs om de hypothese te ondersteunen. Maar het teken van $b_1$ is wel positief, en de $p$ -waarde zit dichtbij $.05$ — er zit iets in de data dat de hypothese-richting wijst. Mogelijk wordt het effect van afwijzingen op de log-odds van turen wel significant als we andere voorspellers toevoegen (in het volle model hierna komt aantal_afwijzingen als deel van een groter model wel uit met $p = .032$ ). Een voorspeller alleen vertelt zelden het hele verhaal.

Inzicht. De Wald-toets ( $z$ of $\chi^2(1)$ ) zegt iets over één specifieke coëfficiënt in het model dat je nu hebt gefit. Bij een enkel-voorspeller-model zit veel ruis nog in de residuen — het model verklaart maar een klein stukje. In het volle model hierna breiden we uit met leeftijd en diersoort, en dan komt het beeld scherper.

4.1.g Likelihood-ratio-toets — doet het model als geheel iets?

De Wald-toets kijkt per coëfficiënt. De likelihood-ratio-toets (LR) kijkt naar het model als geheel: brengt onze voorspeller aantal_afwijzingen de voorspellingen dichter bij de echte data dan een model zonder voorspellers (het null-model)?

Verhaalfiguur — drie steden, drie modellen, één lijn

Stel je drie steden voor op een rechte lijn: Leiden, Voorhout, Haarlem. En ergens voorbij Haarlem ligt de data — de werkelijkheid.

Stad	Model	Wat het doet
Leiden	$M_0$ — null-model	gokt $P(Y=1)$ als één algemeen gemiddelde. Geen voorspellers.
Voorhout	$M_1$ — model met één voorspeller	bv. alleen `aantal_afwijzingen`. Iets dichter bij data.
Haarlem	$M_{\text{vol}}$ — volledig model	alle voorspellers (komt straks, in het volle model). Het dichtst bij data.

De deviance ( $-2LL$ ) is de afstand van dat model tot de data. Hoe groter, hoe verder weg. Het null-model in Leiden ligt het verst. Voeg een goede voorspeller toe, en je schuift op naar Voorhout. Een volledig model brengt je naar Haarlem.

De LR-toets vraagt: hoeveel kilometer zijn we opgeschoten door deze voorspeller toe te voegen — en is dat opvallend meer dan we van een willekeurige uitbreiding zouden verwachten?

Concreet voor onze pad:

# Null-model — alleen intercept (= Leiden):
m_nul <- glm(bekeek_zichzelf ~ 1,
             data   = pad_spiegel,
             family = binomial(link = "logit"))

# LR-toets: model met aantal_afwijzingen versus model zonder:
anova(m_nul, m_afw, test = "Chisq")

Analysis of Deviance Table

Model 1: bekeek_zichzelf ~ 1
Model 2: bekeek_zichzelf ~ aantal_afwijzingen
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
1       151     210.48                       
2       150     207.40  1   3.0796  0.07928 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Drie namen voor één ding — deviance,

-2LL

, residual deviance

In formules, R-output en handboeken kom je drie verschillende termen tegen die naar hetzelfde getal verwijzen:

Hoe genoemd	Waar je het ziet
$-2LL$ (of $-2\log L$ )	Formules, SPSS-uitvoer
Deviance	Handboeken, artikelen
Residual deviance	R-output van `summary(glm(...))`

$LL$ staat voor log-likelihood: hoe aannemelijk de data zijn onder het model. Keer $-2$ maakt er een positief getal van dat je als afstand leest. Hoe kleiner, hoe dichter bij de data. Voor onze pad:

    Null deviance: 210.48  on 151  degrees of freedom
Residual deviance: 207.40  on 150  degrees of freedom

$\chi^2_{\text{LR}} = 210.48 - 207.40 = 3.08$ — de afstand-verbetering van Leiden naar Voorhout.

Let op in de uitvoer van anova() hierboven: daar staan de twee afstanden in de kolom Resid. Dev. De kolom die Deviance heet, is hun verschil, de $\chi^2_{\text{LR}}$ zelf.

Vraag 4.1.g — Pen-en-papier

a) Lees uit de R-output (of het kader Drie namen voor één ding hierboven) de twee deviances: null en residual.

b) Bereken $\mathit{\chi}^2_{\text{LR}} = \text{Dev}_{\text{null}} - \text{Dev}_{\text{model}}$ handmatig.

c) Hoeveel vrijheidsgraden heeft de LR-toets hier? (Hint: $df$ = aantal toegevoegde voorspellers.)

d) Uit de R-output is $p \approx .079$ . Verwerp je $H_0$ op $\alpha = .05$ ?

e) Vergelijk je $\chi^2_{\text{LR}}$ met de Wald- $\chi^2(1) \approx 2.97$ uit 4.1.f. Geven beide toetsen dezelfde conclusie? Wat valt je op aan het kleine verschil in getallen?

Antwoord 4.1.g — open na je eigen poging

a) Null deviance $= 210.48$ ; residual deviance van het model met aantal_afwijzingen $= 207.40$ .

b) $\chi^2_{\text{LR}} = 210.48 - 207.40 = 3.08$ — let op: gebruik onafgeronde getallen anders klopt de aftrekking niet (de sudoku-eigenschap; die komt terug bij de modelvergelijking in het volle model).

c) $df = 1$ — we voegen één voorspeller toe (aantal_afwijzingen).

d) $p \approx .079 > .05$ . We verwerpen $H_0$ niet. De LR-toets ondersteunt op $\alpha = .05$ niet de hypothese dat het model significant beter is dan het null-model.

e) Beide toetsen geven $p \approx .08$ — dezelfde conclusie (niet-significant). Wald- $\chi^2(1) = 2.97$ versus LR- $\chi^2(1) = 3.08$ verschillen marginaal. Bij grote $N$ en simpele modellen geven Wald en LR vrijwel dezelfde uitkomst. Het kleine verschil komt doordat ze via andere formules tot een toets-statistiek komen — Wald via $b/SE$ , LR via deviance-aftrekken. In ons geval levert dat zo’n tiende verschil op. Verwaarloosbaar voor de conclusie.

Inzicht. Wald en LR doen op model-niveau vrijwel hetzelfde werk. Het verschil wordt merkbaar bij (1) factoren met meerdere niveaus — daar bundelt LR alle dummies in één toets, terwijl Wald per dummy een aparte toets levert; en (2) kleine steekproeven of heel grote effecten — daar kunnen de twee uiteenlopen. Vuistregel: bij twijfel, en altijd voor factoren als geheel, LR.

We hebben nu de drie rekenstappen (kans → odds → OR) én de twee toetsen (Wald, LR) op één voorspeller doorlopen. Het effect van aantal_afwijzingen alleen blijkt marginaal niet-significant ( $p \approx .08$ ). Dat is niet “geen effect” — dat is “in deze data niet sterk genoeg om met zekerheid te concluderen, gegeven alleen deze ene voorspeller”. Hierna breidt het model zich uit naar leeftijd en diersoort — en daar komt het beeld scherper uit.

4.2 Het volle model — meerdere voorspellers tegelijk

Tellen wat voorbij komt — en wat voorbij gaat

“Eén bit per dier,” zei de pad. “Maar als ik weet bij welke combinatie het spiegelen wél gebeurt, weet ik morgen wie er nieuwsgierig naar zijn eigen smoeltje zal zijn.”

Voor je gaat rekenen — vijf vragen aan jezelf

Het stappenplan dat in elk thema terugkomt, met de techniek-keuze als laatste stap.

Wie of wat wordt er gemeten? De pad heeft $152$ dieren genoteerd.
Wat wordt er gemeten? Vier variabelen per dier: bekeek_zichzelf, leeftijd, aantal_afwijzingen, diersoort.
Onafhankelijk of afhankelijk? Onderzoeksvraag: voorspellen leeftijd, aantal afwijzingen en diersoort of een dier zichzelf bekijkt? Dan is bekeek_zichzelf afhankelijk ( $Y$ ), de andere drie onafhankelijk ( $X$ ).
Meetniveau van elke variabele? bekeek_zichzelf binair, leeftijd interval, aantal_afwijzingen interval, diersoort nominaal (3 niveaus).
Welke techniek? Doorloop de beslisboom:

flowchart TD
    A[Hoeveel afhankelijke<br/>variabelen?] -->|één| B[Y meetniveau?]
    A -->|meerdere| Z[MANOVA<br/><i>thema 5</i>]
    B -->|interval| C[X-en meetniveau?]
    B -->|binair| D[Logistic regression<br/><i>dit thema</i>]
    C -->|alleen interval| E[Multiple regression<br/><i>thema 1</i>]
    C -->|alleen nominaal| F[ANOVA<br/><i>thema 2</i>]
    C -->|gemengd| G[ANCOVA<br/><i>thema 3</i>]
    style D fill:#faf3e2,stroke:#c9a05a,stroke-width:2px

Pijltjes-overzicht — alle MVDA-technieken in één tabel

Dezelfde keuze, anders genoteerd. Lees als predictor(en) (meetniveau) → DV(s) (meetniveau):

#	Predictor(en)	DV(s)	Techniek
	dich \|X\|	\|Y\| int	$t$ -toets
	nom 3+lvl \|X\|	\|Y\| int	éénweg ANOVA
1	int \|X X …\|	\|Y\| int	MRA — thema 1
2	nom \|X X\|	\|Y\| int	factorial ANOVA — thema 2
3	nom + int \|X C\|	\|Y\| int	ANCOVA — thema 3
4	bo(in) \|X\|	\|Y\| bin	LRA — thema 4
5	nom \|X\|	\|Y₁ Y₂ …\| int	MANOVA — thema 5
6	within-subject momenten	\|Y₁ Y₂ Y₃ Y₄\| int	RMA — thema 6
7	indirect via $M$	$X \to M \to Y$ pad	Mediation — thema 7

$Y$ is binair (bekeek_zichzelf), $X$ -en zijn gemengd. Meetniveau van $X$ doet er minder toe — wat de techniek bepaalt is de aard van $Y$ . Antwoord: logistic regression.

T1 — Techniek-keuze: LRA, MRA, ANOVA of ANCOVA?

Vraag T1 — Pen-en-papier

Voor elk vignet: bepaal afhankelijke variabele, meetniveau van elke voorspeller, en kies dan de techniek.

a) Een vleermuis registreert per nachtelijke vlucht of hij prooi heeft gevangen (ja/nee), met de continue voorspellers maan-helderheid en jachtduur.

b) Een psycholoog meet bij $100$ studenten hun studie-uren-per-week als functie van slaapkwaliteit en zelf-gerapporteerde motivatie — beide continu.

c) Een sportcoach vergelijkt drie trainingsregimes (traag opbouwend / interval / krachtgericht) op rondetijd ( $s$ ). Bij $60$ lopers.

d) Een kerkuil meet per nest of een kuiken het overleeft (ja/nee), met als voorspellers nesthabitat (drie typen) en gewicht-bij-uitkomst (continu).

Antwoord T1 — open na je eigen poging

a) $Y$ = prooi-gevangen (BIN), twee continue $X$ -en. $\Rightarrow$ LRA.

b) $Y$ = studie-uren (interval), twee continue $X$ -en. $\Rightarrow$ multiple regressie-analyse (thema 1).

c) $Y$ = rondetijd (interval), één factor (trainingsregime, NOM, 3 niveaus). $\Rightarrow$ eenweg-ANOVA (thema 2).

d) $Y$ = overleeft (BIN), één factor (nesthabitat) + één continue covariaat (gewicht). $\Rightarrow$ LRA. Niet ANCOVA: ANCOVA vereist een continue $Y$ . Bij binaire $Y$ wordt het altijd LRA, ongeacht de mix van $X$ -meetniveaus.

Vuistregel. Het meetniveau van $Y$ bepaalt de techniek-familie. Binaire $Y$ $\Rightarrow$ LRA. Continu $Y$ + alleen continue $X$ $\Rightarrow$ MRA. Continu $Y$ + alleen factor(en) $\Rightarrow$ ANOVA. Continu $Y$ + factor + covariaat $\Rightarrow$ ANCOVA.

T2 — Odds versus probability

Vraag T2 — Pen-en-papier

Odds en kansen zijn twee manieren om hetzelfde feit te beschrijven, maar op verschillende schalen.

a) Bij een kans $P = .50$ : wat is de odds? En de log-odds (logit)?

b) Bij $P = .80$ : wat is de odds? En de log-odds?

c) Stel je leest af dat de odds van slagen voor een tentamen $= 3$ . Wat is de bijbehorende kans $P$ ?

d) Bij negatieve log-odds is de kans… groter of kleiner dan $.50$ ?

Antwoord T2 — open na je eigen poging

a) Odds = $.50/.50 = 1$ . Log-odds = $\log(1) = 0$ .

b) Odds = $.80/.20 = 4$ . Log-odds = $\log(4) \approx 1.39$ .

c) $P = \text{odds}/(1 + \text{odds}) = 3/4 = .75$ .

d) Kleiner dan $.50$ . Bij $\text{logit}(P) < 0$ geldt $P < .50$ ; bij $\text{logit}(P) > 0$ geldt $P > .50$ ; bij $\text{logit}(P) = 0$ geldt $P = .50$ .

4.2.a Verkenning — frequenties en wolken

Voordat de pad fit, kijkt hij eerst.

Algemene vorm.

table(mijn_data$Y)
prop.table(table(mijn_data$Y))
table(mijn_data$factor, mijn_data$Y)

# Continue X tegenover Y — boxplot per uitkomst-categorie.
boxplot(X1 ~ Y, data = mijn_data)

Voor onze dieren.

# Aantallen per uitkomst en proportie zelf-bekijkers.
table(pad_spiegel$bekeek_zichzelf)


 0  1 
79 73

prop.table(table(pad_spiegel$bekeek_zichzelf))


        0         1 
0.5197368 0.4802632

# Kruistabel diersoort x bekeek_zichzelf.
table(pad_spiegel$diersoort, pad_spiegel$bekeek_zichzelf)

# count() telt per waarde; mutate() zet de proportie ernaast.
pad_spiegel |>
  count(bekeek_zichzelf) |>
  mutate(proportie = n / sum(n))

  bekeek_zichzelf  n proportie
1               0 79 0.5197368
2               1 73 0.4802632

# Kruistabel: tel per combinatie, en zet de uitkomst daarna in twee kolommen
# naast elkaar (pivot_wider, uit tidyr), zodat hij leest als table().
pad_spiegel |>
  count(diersoort, bekeek_zichzelf) |>
  pivot_wider(names_from = bekeek_zichzelf, values_from = n)

# A tibble: 3 × 3
  diersoort   `0`   `1`
  <fct>     <int> <int>
1 ree          40    24
2 vos           9    46
3 das          30     3

Toon code (twee boxplots naast elkaar)

par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))

boxplot(leeftijd ~ bekeek_zichzelf, data = pad_spiegel,
        names = c("Niet bekeken", "Wel bekeken"),
        ylab = "Leeftijd (jaren)", xlab = "",
        col = c("#BBBBBB", "#0077BB"))

boxplot(aantal_afwijzingen ~ bekeek_zichzelf, data = pad_spiegel,
        names = c("Niet bekeken", "Wel bekeken"),
        ylab = "Aantal afwijzingen", xlab = "",
        col = c("#BBBBBB", "#0077BB"))

Toon code (twee boxplots naast elkaar)

par(mfrow = c(1, 1))

Vragen 4.2 — verkenning (a–d)

a) Hoeveel dieren heeft de pad genoteerd? Welke proportie bekeek zichzelf?

b) Hoe verschilt de mediane leeftijd tussen “wel bekeken” en “niet bekeken”? In welke richting?

c) Idem voor aantal afwijzingen.

d) Lijkt diersoort uit te maken? Welke richting?

Antwoord 4.2 — open na je eigen poging

In gewone woorden.

$N = 152$ , waarvan $73$ zelf-bekijkers ( $48.0\%$ ) en $79$ niet-bekijkers ( $52.0\%$ ) — vrijwel balanced ( $n_{\max}/n_{\min} = 79/73 \approx 1.08$ ).
Op de boxplot ligt de mediane leeftijd bij bekijkers lager dan bij niet-bekijkers — jonger geeft meer spiegel-gedrag.
De afwijzingen-boxplot wijst de andere kant op: bij bekijkers liggen de aantallen afwijzingen iets hoger — wie afgewezen wordt, zoekt het water op om eens goed te kijken wie er voor schut staat.
De kruistabel laat zien dat vossen relatief vaak in de “wel bekeken”-cel zitten ( $46$ van $55$ vossen), dassen bijna nooit ( $3$ van $33$ ), reeën wisselend ( $24$ van $64$ ). Vossen zijn de narcisten van het bos.

APA-stijl.

Van $N = 152$ dieren bekeken $73$ ( $48.0\%$ ) zichzelf in het water. Beschrijvend lag de mediane leeftijd bij bekijkers lager dan bij niet-bekijkers, lag het aantal afwijzingen bij bekijkers iets hoger, en kwam zelfreflectie veel vaker voor bij vossen dan bij dassen.

4.2.b Het volle model fitten — alle voorspellers tegelijk

Algemene vorm.

mijn_model <- glm(Y ~ X1 + X2 + factor,
                  data = mijn_data,
                  family = binomial(link = "logit"))
summary(mijn_model)

De familie binomial(link = "logit") zegt: “ $Y$ is binair, en model dit op de log-odds-schaal.” glm() schat de coëfficiënten op een andere manier dan lm() bij MRA (met maximum likelihood in plaats van met kleinste kwadraten). Wat je daarvan in de uitvoer merkt: per coëfficiënt een $z$ -toets in plaats van een $t$ -toets.

Wald-

z

in summary() versus Wald-

\chi^2(1)

in rapportage

summary() toont per coëfficiënt de Wald- $z$ : $z_j = \dfrac{b_j}{SE_{b_j}}$

Maar in rapportage (zoals in SPSS-output en in veel APA-stijl-papers) staat meestal de Wald- $\chi^2(1)$ — dezelfde toets, andere vorm:

$\mathit{\chi}^2(1)_j = z_j^2$

Mathematisch identiek. Een $z = 2.5$ wordt $\chi^2(1) = 6.25$ — exact dezelfde $p$ -waarde. De ene is gewoon het kwadraat van de andere. De $df = 1$ omdat we één coëfficiënt toetsen.

Lees-conventie:

Uit summary() lees je z value en Pr(>|z|).
In rapportage schrijf je $\chi^2(1) = \ldots$ , $p = \ldots$ — kwadrateer de $z$ .
Bij twijfel: noteer beide, bijvoorbeeld “ $z = 2.50$ , equivalent $\chi^2(1) = 6.25$ , p = .012”.

Voor onze dieren.

m_pad <- glm(bekeek_zichzelf ~ leeftijd + aantal_afwijzingen + diersoort,
             data = pad_spiegel,
             family = binomial(link = "logit"))
summary(m_pad)


Call:
glm(formula = bekeek_zichzelf ~ leeftijd + aantal_afwijzingen + 
    diersoort, family = binomial(link = "logit"), data = pad_spiegel)

Coefficients:
                   Estimate Std. Error z value    Pr(>|z|)    
(Intercept)        -0.43231    0.58160  -0.743      0.4573    
leeftijd           -0.14707    0.06703  -2.194      0.0282 *  
aantal_afwijzingen  0.29008    0.13499   2.149      0.0316 *  
diersoortvos        2.35041    0.47589   4.939 0.000000785 ***
diersoortdas       -1.74150    0.67515  -2.579      0.0099 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 210.48  on 151  degrees of freedom
Residual deviance: 144.26  on 147  degrees of freedom
AIC: 154.26

Number of Fisher Scoring iterations: 5

Categorische voorspellers met meer dan twee niveaus

diersoort is een factor met drie niveaus (ree, vos, das). Wat doet R daarmee in glm()?

R maakt automatisch dummy-variabelen (treatment-coding): voor $K$ niveaus krijg je $K - 1$ dummies, één per niet-referentieniveau.
De referentiecategorie is het eerste level van de factor (hier: ree). R gebruikt levels(factor)[1] als referentie. Bij factor() zonder expliciet levels=-argument is dat alfabetisch eerste; in onze data wordt de level-volgorde expliciet gezet (ree, vos, das). Wijzigen kan via pad_spiegel$diersoort <- relevel(pad_spiegel$diersoort, ref = "vos").
In summary() zie je per niet-referentieniveau één coëfficiënt + Wald-toets — dat is een vergelijking met de referentie, niet een toets op het niveau zelf.
Eén OR per dummy = factor waarmee de odds vermenigvuldigen vergeleken met de referentie.
Voor de vraag “is de factor als geheel voorspellend?” gebruik je anova(model_zonder, model_met, test = "Chisq") (LR-toets) of car::Anova(model, type = "II", test.statistic = "LR"). Dat bundelt alle dummies in één $\chi^2$ -toets met $df = K - 1$ .

Mini-oefening. Maak een kopie van de data met vos als referentiecategorie (pad_vos <- pad_spiegel; pad_vos$diersoort <- relevel(pad_vos$diersoort, ref = "vos")) en fit het model opnieuw op pad_vos; zo blijft pad_spiegel zelf zoals hij was, en gaan de vragen hieronder nog over hetzelfde model — om te zien hoe Wald-toetsen veranderen ten opzichte van de huidige ree-referentie. Hoe verandert de OR voor ree? En voor das? Wat blijft hetzelfde aan de output van car::Anova(., test = "LR")? Inzicht: dummies + Wald-toetsen veranderen, de LR-toets op de hele factor niet — die is referentie-onafhankelijk.

Vragen 4.2 — volle model (e–g)

e) Welke voorspellers zijn significant volgens de Wald-toets ( $p < .05$ )?

f) Wat is de geschatte log-odds-coëfficiënt voor leeftijd? Voor aantal_afwijzingen? Wat zegt het teken? (Het getal $b$ zelf lees je niet als effect, zijn teken wel.)

g) Wat doet diersoortvos ten opzichte van diersoortree (de referentie)? En diersoortdas?

Antwoord 4.2 — open na je eigen poging

In gewone woorden.

Significant (p < .05): leeftijd (summary(): $z = −2.19$ , equivalent $\chi^2(1) = 4.81$ , p = .028), aantal_afwijzingen ( $z = 2.15$ / $\chi^2(1) = 4.62$ , p = .032), de dummy diersoortvos ( $z = 4.94$ / $\chi^2(1) = 24.39$ , p < .001) en diersoortdas ( $z = −2.58$ / $\chi^2(1) = 6.65$ , p = .010). Alle vier de effecten staan vlaggetjes.
Leeftijd $b = -0.147$ — negatief teken: jonger geeft meer spiegel-gedrag, een puberale narcisten-trek. Aantal_afwijzingen $b = 0.290$ — positief teken: wie afgewezen wordt, zoekt het water op.
diersoortvos-coëfficiënt $2.350$ — positief, dus vossen bekijken zichzelf veel vaker dan reeën (de referentie). diersoortdas $-1.742$ — sterk negatief, dassen tonen weinig zelfreflectie. De pad zag een tweedeling tussen vos en das: vossen, blijkbaar, hebben iets te bewijzen.

APA-stijl (we komen er straks mooier op terug, in de APA-zin aan het eind van deze sectie).

Een logistische regressie van zelf-bekijken op leeftijd, aantal afwijzingen en diersoort toonde een negatief leeftijdseffect ( $\chi^2(1) = 4.81$ , p = .028) en een positief afwijzingseffect ( $\chi^2(1) = 4.62$ , p = .032); vossen onderscheidden zich sterk van reeën ( $\chi^2(1) = 24.39$ , p < .001).

T3 — $\mathit{OR}$ uit $b$ via $\exp()$

Vraag T3 — Pen-en-papier

De odds ratio is $\mathit{OR} = \exp(b)$ . Reken zelf:

a) $b = 0$ . Wat is $\mathit{OR}$ ? Interpreteer.

b) $b = 0.69$ . Wat is $\mathit{OR}$ (afgerond)? Interpreteer.

c) $b = -0.69$ . Wat is $\mathit{OR}$ ? Interpreteer.

d) Voor de pad-data: aantal_afwijzingen $b = 0.290$ . Wat is $\mathit{OR}$ per afwijzing? Hint: $\exp(0.290) \approx 1.336$ .

Antwoord T3 — open na je eigen poging

a) $\mathit{OR} = \exp(0) = 1$ . Geen effect: bij één-eenheid stijging in $X$ blijven de odds gelijk.

b) $\mathit{OR} = \exp(0.69) \approx 2.0$ . Per één-eenheid stijging in $X$ verdubbelen de odds.

c) $\mathit{OR} = \exp(-0.69) \approx 0.50$ . Per één-eenheid stijging halveren de odds. Spiegelbeeld van b.

d) $\mathit{OR}_{\text{afw}} = \exp(0.290) \approx 1.336$ . Per extra afwijzing worden de odds van zelfreflectie met factor $1.336$ vermenigvuldigd — een stijging van ruim $33\%$ per afgewezen vibe.

T4 — $\mathit{OR}(k) = \mathit{OR}(1)^k$ : meerdere eenheden

Vraag T4 — Pen-en-papier

De OR voor één-eenheid stijging is $\mathit{OR}(1) = \exp(b)$ . Voor $k$ -eenheden stijging: $\mathit{OR}(k) = \exp(k \cdot b) = \mathit{OR}(1)^k$

Voor de pad-data: $\mathit{OR}_{\text{afw}}(1) \approx 1.34$ per afwijzing.

a) Wat is de OR voor een verschil van $3$ afwijzingen?

b) Wat is de OR voor een verschil van $5$ afwijzingen?

c) Stel een dier had vorige maand $0$ afwijzingen, deze maand $2$ extra. Hoeveel keer hoger zijn de odds van zelfreflectie, alle andere voorspellers gelijk?

d) Voor leeftijd: $\mathit{OR}(1) \approx 0.86$ per jaar (negatief teken — jonger geeft meer spiegel-gedrag). Wat is de OR voor $5$ jaar verschil?

Antwoord T4 — open na je eigen poging

a) $\mathit{OR}(3) = 1.34^3 \approx 2.41$ . Bij drie afwijzingen meer worden de odds bijna $2.4$ keer zo groot.

b) $\mathit{OR}(5) = 1.34^5 \approx 4.32$ . Vijf afwijzingen meer maakt de odds ruim $4$ keer zo groot — dat lijkt extreem; zo werkt exponentieel rekenen.

c) $\mathit{OR}(2) = 1.34^2 \approx 1.80$ . De odds zijn bijna verdubbeld na twee extra afwijzingen.

d) $\mathit{OR}(5) = 0.86^5 \approx 0.47$ . Bij vijf jaar oudere dieren halveren de odds van spiegel-gedrag ongeveer.

Inzicht. Odds rekenen vermenigvuldigend, niet optellend. Per eenheid telt $b$ op in log-odds, vermenigvuldigt $\mathit{OR}$ in odds. Vandaar de exponentiële groei (of krimp) bij meer eenheden.

Bij welk

x

P = .5

? — de drempel-truc (afleiding in drie stappen)

Dit is de drempel-truc uit 4.1.c nog één keer, nu afgeleid in logit-taal, want hij komt regelmatig terug op tentamens: “voor welke waarde van $X$ is de slagings-kans precies 50%?” Dat punt is de natuurlijke drempel van het model — daar slaat het van “voorspelt zakken” naar “voorspelt slagen” om. Komt straks terug bij de classificatie (de cutoff).

Voor een model met één voorspeller, $P(Y=1) = \dfrac{e^{b_0 + b_1 X}}{1 + e^{b_0 + b_1 X}}$ , leid je dat punt zo af:

Stap 1. $P = .5$ correspondeert met odds $= 1$ , en odds $= 1$ correspondeert met $\text{logit}(P) = \log(1) = 0$ . Dus de logit-vorm wordt: $b_0 + b_1 X = 0$

Stap 2. Los op naar $X$ . Sleep $b_0$ naar de andere kant: $b_1 X = -b_0$

Stap 3. Deel beide kanten door $b_1$ : $\boxed{\,X^* = -\frac{b_0}{b_1}\,}$

Bv. bij een fictief logistic-regressie-model met $b_0 = -4.663$ en $b_1 = 0.707$ : $X^* = -(-4.663)/0.707 = 6.59$ . Bij $X = 6.59$ is $P(Y=1) = .5$ . Voor $X < 6.59$ voorspelt het model “zakken”, voor $X > 6.59$ “slagen”.

Voor een model met meerdere voorspellers werkt het ook, maar dan moet je de andere $X$ ’s vastzetten op een waarde naar keuze (bv. hun gemiddelde, of een specifiek profiel). De formule blijft dezelfde structuur: zet logit op nul, los op naar $X_j$ .

4.2.c Odds ratios en hun 95%-CI — meerdere voorspellers naast elkaar

Algemene vorm.

exp(coef(mijn_model))           # Odds ratios per coëfficiënt
exp(confint(mijn_model))        # 95%-CI op odds-schaal (profile-likelihood)
exp(confint.default(mijn_model)) # 95%-CI op odds-schaal (Wald, sneller)

confint.default() geeft het Wald-interval: precies $b \pm 1.96 \cdot SE_b$ , dezelfde som die je met de hand maakt, en met exp() eromheen het interval van de odds ratio. confint() rekent anders (via de zogeheten profile likelihood) en geeft daardoor net iets andere grenzen; bij grote $N$ is het verschil klein. In dit hoofdstuk gebruiken we confint.default().

Voor onze dieren.

# Odds ratios.
exp(coef(m_pad))

       (Intercept)           leeftijd aantal_afwijzingen       diersoortvos 
         0.6490102          0.8632319          1.3365355         10.4898740 
      diersoortdas 
         0.1752566

# 95%-CI op odds-schaal (Wald).
exp(confint.default(m_pad))

                        2.5 %     97.5 %
(Intercept)        0.20758278  2.0291384
leeftijd           0.75695243  0.9844336
aantal_afwijzingen 1.02583858  1.7413335
diersoortvos       4.12756519 26.6591687
diersoortdas       0.04666461  0.6582049

OR-tabel — alles op één plek

Toon code (gt-tabel met b, SE, z, p, OR, CI)

term_labels <- c(
  "(Intercept)"        = "Intercept",
  "leeftijd"           = "Leeftijd (jaren)",
  "aantal_afwijzingen" = "Aantal afwijzingen",
  "diersoortvos"       = "Diersoort: vos vs ree",
  "diersoortdas"       = "Diersoort: das vs ree"
)

co  <- summary(m_pad)$coefficients
ci  <- confint.default(m_pad)
or  <- exp(co[, "Estimate"])
or_lo <- exp(ci[, 1])
or_hi <- exp(ci[, 2])

# De getallen worden tekst via de fmt_*-functies van de huisstijl (KAAPA):
# twee decimalen, een echt minteken, en een p zonder voorloopnul.
or_df <- data.frame(
  term = term_labels[rownames(co)],
  b    = fmt_b(co[, "Estimate"]),
  se   = fmt_se(co[, "Std. Error"]),
  z    = fmt_t(co[, "z value"]),
  p    = palign(sapply(co[, "Pr(>|z|)"], fmt_p)),
  OR   = fmt_b(or),
  ci   = fmt_ci(or_lo, or_hi),
  row.names = NULL
)

or_df |>
  gt() |>
  fmt_markdown(columns = p) |>   # fmt_p() levert "&lt; .001"; zo wordt dat weer "<"
  cols_align("right", columns = c(b, se, z, p, OR, ci)) |>
  cols_label(
    term = "",
    b    = md("*b*"),
    se   = md("*SE*"),
    z    = md("Wald *z*"),
    p    = md("*p*"),
    OR   = md("*OR*"),
    ci   = md("95%-CI voor *OR*")
  ) |>
  gt_apa(title = "Wie van de dieren naar zichzelf tuurt in de spiegel: coëfficiënten en odds ratios") |>
  tab_footnote_apa(paste0(
    "Logistische regressie van zelf in de spiegel kijken (1 = ja) op leeftijd, ",
    "aantal afwijzingen en diersoort (referentie: ree), *N* = ", nobs(m_pad), ". ",
    "*b* = coëfficiënt op log-odds-schaal; *SE* = standaardfout; ",
    "*OR* = odds ratio, exp(*b*); CI = betrouwbaarheidsinterval (Wald)."))

	b	SE	Wald z	p	OR	95%-CI voor OR
Wie van de dieren naar zichzelf tuurt in de spiegel: coëfficiënten en odds ratios
Intercept	−0.43	0.58	−0.74	.46	0.65	[0.21, 2.03]
Leeftijd (jaren)	−0.15	0.07	−2.19	.028	0.86	[0.76, 0.98]
Aantal afwijzingen	0.29	0.13	2.15	.032	1.34	[1.03, 1.74]
Diersoort: vos vs ree	2.35	0.48	4.94	< .001	10.49	[4.13, 26.66]
Diersoort: das vs ree	−1.74	0.68	−2.58	.010	0.18	[0.05, 0.66]
Note. Logistische regressie van zelf in de spiegel kijken (1 = ja) op leeftijd, aantal afwijzingen en diersoort (referentie: ree), N = 152. b = coëfficiënt op log-odds-schaal; SE = standaardfout; OR = odds ratio, exp(b); CI = betrouwbaarheidsinterval (Wald).

Met de hand. Voor een $95\%$ -betrouwbaarheidsinterval wandel je vanaf $b$ $1.96$ standaardfouten naar links en $1.96$ naar rechts. Dat geeft het interval van $b$ ; daarna zet je beide grenzen met $e^{\ldots}$ om naar de odds ratio. Voor leeftijd, met $b$ en $SE$ uit de uitvoer van summary() in 4.2.b, op drie decimalen:

stapgrootte: $1.96 \cdot 0.067 = 0.13132$
ondergrens van $b$ : $-0.147 - 0.13132 = -0.27832$
bovengrens van $b$ : $-0.147 + 0.13132 = -0.01568$
interval van de OR: $[e^{-0.27832},\; e^{-0.01568}] = [0.757,\; 0.984]$

Rond pas aan het eind af; dan kom je uit op wat R ook geeft. Het interval van $b$ bevat de $0$ niet, dus het interval van de OR bevat de $1$ niet: het effect van leeftijd op de log-odds van turen wijkt significant van nul af.

Vragen 4.2 — odds ratios (h–j)

h) Wat is de OR voor leeftijd? Voor aantal_afwijzingen? Wat betekent dat in plain Nederlands?

i) Reken het $95\%$ -interval van de OR voor aantal_afwijzingen met de hand uit, zoals hierboven voor leeftijd. Bevatten de intervallen van leeftijd en aantal_afwijzingen de waarde $1$ ? Wat zegt dat over significantie?

j) Wat is de OR voor diersoort vos versus ree? Wat betekent een OR boven $1$ ?

Antwoord 4.2 — open na je eigen poging

In gewone woorden.

$\mathit{OR}_{\text{leef}} = 0.86$ — per jaar dat een dier ouder is, worden de odds van zelf-bekijken met factor $0.86$ vermenigvuldigd, oftewel $14\%$ lager per jaar. $\mathit{OR}_{\text{afw}} = 1.34$ — per extra afwijzing $34\%$ hogere odds van zelfreflectie.
Met de hand voor aantal_afwijzingen:

stapgrootte: $1.96 \cdot 0.135 = 0.2646$
ondergrens van $b$ : $0.290 - 0.2646 = 0.0254$
bovengrens van $b$ : $0.290 + 0.2646 = 0.5546$
interval van de OR: $[e^{0.0254},\; e^{0.5546}] = [1.026,\; 1.741]$

De OR-tabel geeft, afgerond op twee decimalen, $[1.03, 1.74]$ . 95%-CI leeftijd $[0.76, 0.98]$ — bevat $1$ niet, dus significant. CI aantal_afwijzingen — ook niet $1$ , ook significant. Een CI die $1$ bevat zou betekenen “geen effect kan niet uitgesloten worden”.

$\mathit{OR}_{\text{vos vs ree}} = 10.49$ — vossen hebben odds van zelfreflectie die ruim tien keer hoger liggen dan reeën, mooie sterke zelfbevestigings-energie. Een OR boven $1$ betekent: deze categorie heeft hogere odds dan de referentie.

APA-stijl.

Per jaar dat een dier ouder is, lagen de odds van zelf-bekijken $0.86$ keer zo hoog ( $95\%$ -CI $[0.76, 0.98]$ ). Per extra afwijzing waren de odds $1.34$ keer hoger ( $[1.03, 1.74]$ ). Vergeleken met reeën hadden vossen odds van $10.49$ ( $[4.13, 26.66]$ ), terwijl dassen de odds verlaagden tot $0.18$ ( $[0.05, 0.66]$ ).

4.2.d Modelvergelijking — likelihood-ratio-toets bij meerdere voorspellers

Hoe ver is dit model van waar de data zijn?

“Hou eens vast,” zei de pad. “Afstanden. Hoe ver is mijn vijver van Leiden? Een eind. Hoe ver is Leiden van Voorhout? Tien kilometer. Hoe ver is Voorhout van Haarlem? Twintig. Drie steden op één lijn. Dat zijn drie modellen.”

Dit zijn de drie steden uit 4.1.g: Leiden is het null-model, Voorhout het model met één voorspeller, en nu komt Haarlem in zicht, het volle model. De deviance ( $-2LL$ ) is nog steeds de afstand tot de data, en de likelihood-ratio-toets vraagt nog steeds hoeveel kilometer je bent opgeschoten.

De deviances van het volle model

Wat R laat zien (knip uit summary(m_pad)):

    Null deviance: 210.48  on 151  degrees of freedom
Residual deviance: 144.26  on 147  degrees of freedom

Null deviance = $-2LL_0$ (afstand van het null-model tot de data)
Residual deviance = $-2LL_k$ (afstand van het volle model tot de data)
$\chi^2_{\text{LR}}$ = het verschil = $210.48 - 144.26 = 66.22$ — hoeveel kilometer dichterbij data je bent gekomen door voorspellers toe te voegen

In rapportage: schrijf gewoon deviance (of $-2LL$ in formule-vorm). Residual deviance is R-jargon, niet APA.

Wald versus LR — twee toetsen, één doel

Bij MRA gebruik je $t$ -toetsen op coëfficiënten en $F$ -toetsen op (sub)modellen. Bij LRA bestaan beide soorten ook:

Wald- $z$ per coëfficiënt: vergelijkt $b_j$ met $0$ via zijn $SE$ . Snel, direct uit summary().
Likelihood ratio likelihood ratio tussen twee geneste modellen: $\mathit{\chi}^2_{\text{LR}} = -2LL_{\text{simpel}} - (-2LL_{\text{complex}})$ . Theoretisch nauwkeuriger, vooral bij factoren met meerdere niveaus of bij lage $N$ . Bij grote $N$ geven ze meestal vergelijkbare $p$ -waardes.

Voor een factor met $K$ niveaus toetst de LR-toets alle $K-1$ dummies tegelijk — één $\chi^2_{\text{LR}}$ met $df = K-1$ . Dat lees je niet af uit summary(); daarvoor heb je anova(m_zonder, m_met, test = "Chisq") nodig.

Algemene vorm.

# Volledig model versus null-model — toetst of het model als geheel iets doet.
m_null <- glm(Y ~ 1, data = mijn_data, family = binomial(link = "logit"))
anova(m_null, mijn_model, test = "Chisq")

# LR-toets per term (model met- vs zonder die term).
car::Anova(mijn_model, type = "II", test.statistic = "LR")

Voor onze dieren.

# Null-model (alleen intercept).
m_null <- glm(bekeek_zichzelf ~ 1,
              data = pad_spiegel,
              family = binomial(link = "logit"))

# LR-toets: doet het volledige model iets, ten opzichte van het null-model?
anova(m_null, m_pad, test = "Chisq")

Analysis of Deviance Table

Model 1: bekeek_zichzelf ~ 1
Model 2: bekeek_zichzelf ~ leeftijd + aantal_afwijzingen + diersoort
  Resid. Df Resid. Dev Df Deviance           Pr(>Chi)    
1       151     210.48                                   
2       147     144.26  4   66.215 0.0000000000001427 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# LR-toets per term (factor diersoort als geheel).
Anova(m_pad, type = "II", test.statistic = "LR")

Analysis of Deviance Table (Type II tests)

Response: bekeek_zichzelf
                   LR Chisq Df         Pr(>Chisq)    
leeftijd              5.025  1            0.02498 *  
aantal_afwijzingen    5.043  1            0.02472 *  
diersoort            58.097  2 0.0000000000002423 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Geneste modellen — twee voorwaarden voor de LR-toets

Vraag eerst: wanneer is een vogeltje lekker genest?

Als hij er niet doorheen pleurt. Het nestje draagt elk takje waar de vogel op landt. Niet helemaal de juiste vergelijking voor de wiskunde — wel een geheugensteuntje waar je altijd om kunt lachen. Onthou dat: genest = draagvermogen.

De LR-toets vergelijkt twee modellen, $M_{\text{simpel}}$ en $M_{\text{complex}}$ . Hij werkt alleen als de modellen genest zijn:

Op dezelfde data gefit (zelfde aantal observaties — anders zijn de likelihoods niet vergelijkbaar).
$M_{\text{simpel}}$ is een deelverzameling van $M_{\text{complex}}$ qua voorspellers (alle voorspellers in $M_{\text{simpel}}$ zitten ook in $M_{\text{complex}}$ , plus nog één of meer).

Twee modellen met overlappende-maar-niet-geneste voorspellers toets je niet met de LR-toets. Bij anova(m_null, m_pad, test = "Chisq") moet R-output óók kloppen qua $df$ van het verschil: $df_{\text{LR}} = k_{\text{complex}} - k_{\text{simpel}}$ .

Deviances onafgerond rapporteren — anders breekt de sudoku

De LR-toetsstatistiek is een verschil tussen twee deviances:

$\mathit{\chi}^2_{\text{LR}} = (-2LL_{\text{simpel}}) - (-2LL_{\text{complex}}) = \text{Deviance}_{\text{simpel}} - \text{Deviance}_{\text{complex}}$

Als je in een tabel of paper de twee deviances afzonderlijk rapporteert én daarnaast de $\chi^2_{\text{LR}}$ , dan moet de aftrekking kloppen. Bv. bij onze pad-data:

Model	$-2LL$
De twee deviances van de pad en hun verschil
$M_0$ (null)	210.48
$M_{\text{vol}}$ (alle voorspellers)	144.26
$\Delta = \chi^2_{\text{LR}}$	66.22
Note. $M_0$ = nulmodel, alleen een intercept; $M_{\text{vol}}$ = het model met leeftijd, aantal afwijzingen en diersoort. $LL$ = log-likelihood. Het verschil is de likelihood-ratio-toets van $M_{\text{vol}}$ tegen $M_0$ : $\chi^2$ (4) = 66.22, p < .001.

Reken zelf: $210.48 - 144.26 = 66.22$ . ✓ Sudoku-eigenschap.

Dit gaat mis als je rondt voordat je aftrekt. Bij onze pad gaat het toevallig goed: afgerond op gehele getallen staat er $210 - 144 = 66$ , en $66.22$ wordt afgerond ook $66$ . Maar neem twee andere deviances, $198.47$ en $131.56$ : afgerond staat er $198 - 132 = 66$ , terwijl het echte verschil $66.91$ is en afgerond $67$ . Een lezer die de drie getallen ziet en wil controleren krijgt dan een aftrekking die er net naast zit. Voor APA-rapportage in een paper:

Liever: deviances rapporteren op $2$ decimalen, ook al staat in jouw R-output meer.
Beter nog: deviances alleen rapporteren als ze relevant zijn voor wat je beweert; vaak is alleen $\chi^2_{\text{LR}}(df) = \ldots, p = \ldots$ genoeg.

(Uitvoer in een opgave of uitwerking is soms al afgerond; dan kun je de aftrek niet meer precies narekenen. Niet jouw fout, wél een leesmoeilijkheid.)

Vrijheidsgraden bij LRA

Bij Wald per coëfficiënt en LR-toets per term:

Toets	$df$
Wald per coëfficiënt	$1$ (normaal verdeeld onder $H_0$ )
LR-toets voor één continue voorspeller	$1$
LR-toets voor factor met $K$ niveaus	$K - 1$
LR-toets volledig vs null-model	$k$ (totaal aantal voorspellers in het model)

Hoofdrekenbaar voorbeeld. Drie voorspellers (twee continu, één factor met $3$ niveaus) leveren in totaal $1 + 1 + (3-1) = 4$ vrijheidsgraden voor de LR-toets volledig vs null. Voor onze pad: leeftijd ( $1$ ) + aantal_afwijzingen ( $1$ ) + diersoort ( $2$ ) = $4$ — dat is precies het $df$ in de anova()-output.

Vragen 4.2 — modelvergelijking (k–m)

k) Wat is $\chi^2_{\text{LR}}$ voor het volledige model versus null? Met welk $df$ en welke $p$ ?

l) Is diersoort als geheel significant volgens de LR-toets? Vergelijk met de Wald-toetsen op de losse dummies.

m) Klopt $df = 4$ met de hoofdreken-formule?

Antwoord 4.2 — open na je eigen poging

In gewone woorden.

$\chi^2_{\text{LR}}(4) = 66.22$ , p < .001 — het volledige model verklaart de uitkomst significant beter dan het null-model.
Diersoort als geheel is sterk significant: $\chi^2_{\text{LR}}(2) = 58.10$ , p < .001. De losse Wald-toetsen wijzen dezelfde kant op (vos $p < .001$ , das $p = .010$ ); bij factoren is de LR-toets toch principieel beter omdat hij beide dummies in één hypothese bundelt — vergelijking blijft referentie-onafhankelijk.
$df = 4$ : leeftijd $1$ + aantal_afwijzingen $1$ + diersoort $2$ = $4$ — klopt.

APA-stijl.

Het volledige model verklaarde de uitkomst significant beter dan het null-model, $\chi^2_{\text{LR}}(4) = 66.22$ , p < .001. Bij gecontroleerde toetsing per term waren leeftijd ( $\chi^2_{\text{LR}}(1) = 5.02$ , p = .025), aantal afwijzingen ( $\chi^2_{\text{LR}}(1) = 5.04$ , p = .025) en diersoort ( $\chi^2_{\text{LR}}(2) = 58.10$ , p < .001) elk significant.

T5 — Wald versus LR-toets: wanneer welke?

Vraag T5 — Pen-en-papier

De Wald- $z$ -toets (uit summary()) en de LR- $\chi^2$ -toets (uit anova()) toetsen vaak hetzelfde, maar niet altijd op dezelfde manier.

a) Bij één continue voorspeller met één coëfficiënt: zijn de twee toetsen exact hetzelfde, of ongeveer hetzelfde?

b) Bij een factor met $4$ niveaus: hoeveel Wald- $z$ -waardes zie je in summary()? Hoeveel LR-toetsen?

c) Welke toets is geschikter om “doet deze factor er als geheel toe?” te beantwoorden? Waarom?

Antwoord T5 — open na je eigen poging

a) Ongeveer hetzelfde — beide toetsen $H_0: b_j = 0$ en zijn asymptotisch equivalent. Bij grote $N$ en moderate effecten geven ze vrijwel identieke $p$ -waardes; bij kleine $N$ of grote $|b|$ wijken ze uiteen.

b) Bij een factor met $4$ niveaus: $3$ dummies in summary(), dus $3$ Wald- $z$ ’s — elk een afzonderlijke vergelijking met de referentie. Eén LR-toets in Anova(., test = "LR") met $df = 3$ — de hele factor-bijdrage tegelijk.

c) De LR-toets. De vraag “doet deze factor er toe?” gaat over de hele factor, niet over één specifieke dummy. De Wald per dummy zou je $3$ aparte toetsen geven, met multiple-testing-issues. De LR-toets bundelt ze in één getoetste hypothese.

Inzicht. Wald is comfortabel (komt direct uit summary()); LR is robuuster (vooral bij factoren). In rapportage van LRA: gebruik Wald voor losse continue voorspellers, LR voor factoren als geheel en voor model-vergelijking.

Voor wie verder wil — waarom de Wald-toets bij een heel groot effect kan haperen

Bij een heel sterk effect, en zeker in een kleine steekproef, schat glm() naast een grote $b_j$ ook een grote standaardfout $SE_{b_j}$ . De verhouding $z = b_j / SE_{b_j}$ valt dan kleiner uit dan je zou verwachten: een zeer sterk effect van $X_j$ op de log-odds lijkt volgens de Wald-toets ineens minder significant dan volgens de LR-toets. Dat heet het Hauck-Donner-effect. De LR-toets kijkt alleen naar hoeveel de deviance daalt als $X_j$ erbij komt, en heeft er geen last van. Nog een reden om bij twijfel de LR-toets te nemen.

4.2.e Pseudo- $R^2$

De pseudo- $R^2$ van dit thema, $R^2_L$ , is het deel van de afstand dat het model wegwerkt:

$R^2_L \;=\; \frac{-2LL_0 \;-\; (-2LL_1)}{-2LL_0}$

met $-2LL_0$ de deviance van het null-model en $-2LL_1$ die van jouw model. In de code hieronder staat dezelfde breuk twee keer: één keer met de log-likelihoods zelf (de $-2$ valt in teller en noemer tegen elkaar weg) en één keer met de deviances.

Twee dingen volgen meteen uit die breuk. $R^2_L$ ligt altijd tussen $0$ en $1$ : $0$ als jouw model geen stap dichter bij de data komt dan het null-model, $1$ als de afstand helemaal verdwijnt ( $-2LL_1 = 0$ : elk dier perfect voorspeld). En hij komt zelden hoog, lager dan je van de $R^2$ bij MRA gewend bent. Zelfs een model dat de kans van elk dier precies kent, weet niet of juist dít dier vandaag tuurt: wel of niet blijft een gok, en die rest blijft in de deviance staan.

# Pseudo-R^2 van Hosmer en Lemeshow — handberekend.
LL_null  <- as.numeric(logLik(m_null))
LL_full  <- as.numeric(logLik(m_pad))
pseudoR2 <- 1 - (LL_full / LL_null)
round(pseudoR2, 3)

[1] 0.315

# Equivalent via deviances:
1 - (m_pad$deviance / m_pad$null.deviance)

[1] 0.3145922

Wat het getal betekent

$R^2_L = .315$ : het volle model haalt ruwweg $31\%$ van de deviance van het null-model weg. Dat is de effectgrootte van het model als geheel: hoeveel dichter bij de data kom je met de voorspellers, vergeleken met alleen het intercept? Je leest dit vaak als “het model verklaart $31\%$ van de variantie”. Strikt genomen is het verklaarde deviance en geen verklaarde variantie zoals bij MRA; welk van de twee woorden je ook gebruikt, het gaat om hetzelfde getal. Is dat veel? Voor $R^2_L$ bestaat geen vaste grens voor klein of groot. Gebruik het vooral om modellen op dezelfde data met elkaar te vergelijken: het model met de hoogste $R^2_L$ zit het dichtst bij de data. Het model uit 4.1 met alleen afwijzingen (Voorhout) haalt $R^2_L = (210.48 - 207.40)/210.48 = .015$ ; met leeftijd en diersoort erbij (Haarlem) wordt dat $.315$ .

Voor wie verder wil — wanneer is

R^2_L

groot?

Pseudo- $R^2$ is een andere maat dan de $R^2$ van MRA, met een eigen schaal. Een ruwe vuistregel, geen vaste afspraak: vanaf $.10$ acceptabel, vanaf $.20$ goed, vanaf $.30$ sterk. Met $R^2_L = .315$ valt het model van de pad in de trede sterk.

Twee schalen voor twee maten, dus leg ze niet zomaar naast elkaar. Een $R^2_L = .12$ is op zijn eigen schaal acceptabel. Een MRA- $R^2 = .12$ is volgens Cohen nog klein: zijn grenzen zijn klein $.02$ , matig $.13$ , groot $.26$ .

Deviance-staafgrafiek — visualiseer $R^2_L$

Toon code (illustratie deviance-decompositie)

D0 <- m_pad$null.deviance
Dk <- m_pad$deviance
D_uitleg <- D0 - Dk

dev_df <- data.frame(
  bron = factor(c("Verklaard\n(null - residual)", "Residual deviance"),
                levels = c("Verklaard\n(null - residual)", "Residual deviance")),
  waarde = c(D_uitleg, Dk)
)

ggplot(dev_df, aes(x = bron, y = waarde, fill = bron)) +
  geom_col(width = 0.55) +
  geom_text(aes(label = sprintf("%.2f", waarde)),
            vjust = -0.4, size = 3.6) +
  scale_fill_manual(values = c("#0077BB", "#BBBBBB"), guide = "none") +
  labs(x = NULL, y = "−2 log-likelihood",
       title = "Deviance-decompositie",
       # Twee decimalen: met één (66.2/210.5) komt er .314 uit, niet .315.
       subtitle = paste0("Pseudo-R² = ", sprintf("%.2f", D_uitleg), "/",
                         sprintf("%.2f", D0), " = ",
                         sub("^0", "", sprintf("%.3f", D_uitleg / D0)))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  theme_minimal(base_size = 11)

Wat het plaatje zegt

De totale null-deviance ( $-2LL_0 = 210.48$ ) wordt verdeeld in een verklaard deel (links: $66.22$ ) en een residueel deel (rechts: $144.26$ ). De pseudo- $R^2$ is de verhouding tussen “verklaard” en “totaal”: $66.22/210.48 = .315$ . Het is geen variantie-decompositie zoals in MRA, maar het idee is analoog — wat heeft het model toegevoegd ten opzichte van niets weten?

T6 — Pseudo- $R^2$ : interpretatie en grenzen

Vraag T6 — Pen-en-papier

a) Wat zijn de minimum- en maximumwaarden van Hosmer-Lemeshow pseudo- $R^2$ ?

b) Waarom liggen pseudo- $R^2$ -waardes typisch lager dan MRA’s $R^2$ ?

c) Voor welke vergelijking is pseudo- $R^2$ wel zinvol?

Antwoord T6 — open na je eigen poging

a) Minimum $0$ (het model komt geen stap dichter bij de data dan het null-model), maximum $1$ (perfect model: residual deviance $= 0$ ).

b) Bij een binaire $Y$ blijft er altijd gok over: ook wie de kans van een dier precies kent, weet niet of het dier vandaag tuurt. Die rest blijft in de deviance staan, dus komt de teller van $R^2_L$ zelden in de buurt van de hele noemer.

c) Twee logistische modellen op dezelfde data, met verschillende voorspeller-sets — bijvoorbeeld een model dat stap voor stap voorspellers erbij krijgt. Hoger $R^2_L$ = meer verklaard. Voor cross-data-vergelijking is de absolute waarde minder informatief.

Inzicht. Pseudo- $R^2$ is een vergelijkings-maat, geen absolute prestatie-maat. Combineer met inhoudelijke effectgroottes ( $\mathit{OR}$ ’s) en classificatie-prestatie.

4.2.f Hiërarchische LRA — het model in blokken opbouwen

Eerst wat een dier is, dan wat het meemaakte

“Afwijzingen alleen zeiden te weinig,” zei de pad. “Maar toen wist ik verder ook niets van de dieren. Wat als ik eerst opschrijf wat ik al weet? Doen de afwijzingen er dán nog toe?”

Je kent het uit thema 1, Hiërarchische regressie-analyse: je bouwt het model op in blokken, in een volgorde die je vooraf kiest, op inhoud. Eerst wat je al weet, dan waar je eigenlijk naar vraagt. Bij elke stap stel je twee vragen: is die toevoeging significant? en hoeveel komt erbij? Bij LRA beantwoord je ze met twee dingen uit dit thema.

Is de toevoeging significant? De likelihood-ratio-toets tussen het model vóór en het model ná het blok, met als $df$ het aantal coëfficiënten dat het blok toevoegt. Een voorspeller met één getal (interval of 0/1) kost er één; een factor met drie niveaus kost er twee, één per dummy.
Hoeveel komt erbij? De $R^2_L$ van elk model, uit zijn twee deviances. Die zet je naast elkaar.

De toetsgrootheid van een blok is het verschil in deviance, dezelfde aftrekking als bij elke LR-toets:

$\mathit{\chi}^2_{\text{blok}} = -2LL_{\text{vóór}} - (-2LL_{\text{ná}})$

De volgorde. De pad kiest drie blokken, van wat een dier ís naar wat het onlangs meemaakte: eerst diersoort (ree, vos of das), dan leeftijd, en als laatste zijn eigen vermoeden, aantal afwijzingen. Zo vraagt het laatste blok precies wat hij wil weten: voegen afwijzingen iets toe als je soort en leeftijd al kent?

Algemene vorm.

m0 <- glm(Y ~ 1,            data = mijn_data, family = binomial(link = "logit"))
m1 <- glm(Y ~ X1,           data = mijn_data, family = binomial(link = "logit"))
m2 <- glm(Y ~ X1 + X2,      data = mijn_data, family = binomial(link = "logit"))
m3 <- glm(Y ~ X1 + X2 + X3, data = mijn_data, family = binomial(link = "logit"))

# Alle blokken in één tabel; per paar kan ook: anova(m1, m2, test = "Chisq").
anova(m0, m1, m2, m3, test = "Chisq")

# R2_L van één model, uit zijn twee deviances.
(m1$null.deviance - m1$deviance) / m1$null.deviance

Voor onze dieren.

# Vier modellen, elk één blok groter dan het vorige, allemaal op dezelfde
# dieren: ze zijn dus genest.
m_h0 <- glm(bekeek_zichzelf ~ 1,
            data = pad_spiegel, family = binomial(link = "logit"))
m_h1 <- glm(bekeek_zichzelf ~ diersoort,
            data = pad_spiegel, family = binomial(link = "logit"))
m_h2 <- glm(bekeek_zichzelf ~ diersoort + leeftijd,
            data = pad_spiegel, family = binomial(link = "logit"))
m_h3 <- glm(bekeek_zichzelf ~ diersoort + leeftijd + aantal_afwijzingen,
            data = pad_spiegel, family = binomial(link = "logit"))

# Per blok de LR-toets: elk model tegen het model op de regel erboven.
anova(m_h0, m_h1, m_h2, m_h3, test = "Chisq")

Analysis of Deviance Table

Model 1: bekeek_zichzelf ~ 1
Model 2: bekeek_zichzelf ~ diersoort
Model 3: bekeek_zichzelf ~ diersoort + leeftijd
Model 4: bekeek_zichzelf ~ diersoort + leeftijd + aantal_afwijzingen
  Resid. Df Resid. Dev Df Deviance           Pr(>Chi)    
1       151     210.48                                   
2       149     153.81  2   56.672 0.0000000000004941 ***
3       148     149.31  1    4.500            0.03390 *  
4       147     144.26  1    5.043            0.02472 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# R2_L per model.
r2_L <- function(m) (m$null.deviance - m$deviance) / m$null.deviance
round(c(model_1 = r2_L(m_h1), model_2 = r2_L(m_h2), model_3 = r2_L(m_h3)), 3)

model_1 model_2 model_3 
  0.269   0.291   0.315

Lees de tabel regel voor regel. Regel 1 is het null-model, elke volgende regel is één blok erbij. Per kolom:

Resid. Df: die heb je hier niet nodig.
Resid. Dev: de deviance van het model ná dat blok.
Df: hoeveel coëfficiënten er bij dat blok bijkwamen. Dat is de $df$ van de toets.
Deviance: het verschil met de regel erboven. Dat is de $\chi^2$ van het blok.
Pr(>Chi): de $p$ -waarde van het blok.

Let op de nummers in de kop. R telt vanaf 1, dus daar heet het null-model Model 1. In de tekst heet het null-model model 0, en draagt elk model het nummer van zijn laatste blok: model 1 is alleen diersoort, model 2 diersoort en leeftijd, model 3 alle drie (in de code m_h1 tot en met m_h3).

Met de hand, voor het laatste blok. Model 2 (diersoort en leeftijd) zit op deviance $149.31$ , model 3 (met de afwijzingen erbij) op $144.26$ . Het verschil is $149.31 - 144.26 = 5.05$ . R zegt $5.043$ , omdat R aftrekt met de onafgeronde deviances; dat is precies de val uit het kader Deviances onafgerond rapporteren. Neem daarom de $\chi^2$ uit de kolom Deviance: $\chi^2_{\text{LR}}(1) = 5.04$ , p = .025. En $R^2_L$ van model 3: $(210.48 - 144.26)/210.48 = 66.22/210.48 = .315$ .

Kijk nog eens naar die teller, $66.22$ : dat is alle afstand die het model sinds het null-model gewonnen heeft, en dus de kolom Deviance opgeteld ( $56.67 + 4.50 + 5.04$ , op een afronding na). Elk blok legt er zijn eigen stuk bij.

Zo loopt $R^2_L$ op van $.269$ via $.291$ naar $.315$ , en alle drie de toevoegingen zijn significant. Twee dingen om bij te onthouden:

$R^2_L$ kan alleen stijgen. Een blok erbij maakt de afstand tot de data nooit groter, dus $R^2_L$ stijgt of blijft gelijk, ook als het blok niets betekent. Of de stijging meer is dan toeval, zegt de $\chi^2$ -toets van het blok. Vraagt een opgave welk model volgens $R^2_L$ het beste is, dan is dat het model met de hoogste $R^2_L$ , ook als het laatste blok niet significant was.
Dezelfde voorspeller, een andere vraag. In het begin van dit thema stond aantal afwijzingen er alleen voor, en voegde hij niet significant iets toe aan het null-model ( $\chi^2_{\text{LR}}(1) = 3.08$ , p = .079). Als laatste blok wel. Dat is geen tegenspraak, want het zijn twee vragen: wat doen afwijzingen als je verder niets weet, en wat voegen ze toe als je soort en leeftijd al kent?

Vragen 4.2 — hiërarchische LRA (n–s)

Gebruik de anova()-uitvoer hierboven.

n) Voorspelt diersoort (blok 1) significant of een dier zichzelf bekijkt? Rapporteer $\chi^2$ , $df$ en $p$ . Waarom is $df$ hier $2$ ?

o) Reken $R^2_L$ van model 1 (alleen diersoort) met de hand uit, uit twee deviances.

p) Voegt leeftijd (blok 2) significant iets toe bovenop diersoort? Rapporteer $\chi^2$ , $df$ en $p$ .

q) Reken $R^2_L$ van model 2 (diersoort en leeftijd) met de hand uit.

r) Stel dat de pad nog een vierde blok had toegevoegd, de kleur van de vacht (één 0/1-voorspeller), en dat de deviance daarmee nog $0.80$ zakte, met p = .37 voor dat blok. Wat wordt $R^2_L$ ? Is het vierde blok significant? Welk model is dan volgens $R^2_L$ het beste?

s) Schrijf de uitkomst van de drie blokken op in APA-stijl.

Antwoord 4.2 — open na je eigen poging

In gewone woorden.

Ja: $\chi^2_{\text{LR}}(2) = 56.67$ , p < .001. De $df$ is $2$ omdat diersoort drie niveaus heeft en dus als twee dummies in het model komt (vos en das, elk tegen ree): twee coëfficiënten erbij.
$R^2_L = (210.48 - 153.81)/210.48 = 56.67/210.48 = .269$ .
Ja: $\chi^2_{\text{LR}}(1) = 4.50$ , p = .034.
$R^2_L = (210.48 - 149.31)/210.48 = 61.17/210.48 = .291$ .
De deviance wordt $144.26 - 0.80 = 143.46$ , dus $R^2_L = (210.48 - 143.46)/210.48 = .318$ : hoger dan de $.315$ van model 3. Het blok zelf is niet significant: $\chi^2_{\text{LR}}(1) = 0.80$ , p = .37. Toch is model 4 volgens $R^2_L$ het beste, want het heeft de hoogste waarde. Dat is de valkuil: $R^2_L$ stijgt bij elk blok, ook bij een blok dat niets betekent. Of een blok significant iets toevoegt, zegt alleen zijn $\chi^2$ -toets.

APA-stijl.

In een hiërarchische logistische regressie van zelf-bekijken voorspelde diersoort in de eerste stap of een dier zichzelf bekeek, $\chi^2_{\text{LR}}(2) = 56.67$ , p < .001. Leeftijd voegde daar in de tweede stap significant iets aan toe, $\chi^2_{\text{LR}}(1) = 4.50$ , p = .034, en aantal afwijzingen in de derde stap eveneens, $\chi^2_{\text{LR}}(1) = 5.04$ , p = .025. De $R^2_L$ van het model steeg daarbij van $.27$ via $.29$ naar $.31$ .

Let op bij het afronden naar twee decimalen. De $R^2_L$ van model 3 is $.3146$ ; dat wordt $.31$ . Wie vanuit de al afgeronde $.315$ verder rondt, komt op $.32$ . Rond dus af vanuit het getal van R, niet vanuit een getal dat al een keer afgerond is.

4.2.g Classificatie — kansen, cutoff, accuracy

Algemene vorm.

phat <- predict(mijn_model, type = "response")     # Voorspelde kansen
yhat <- ifelse(phat >= 0.5, 1, 0)                  # Cutoff op 0.5
table(observed = mijn_data$Y, predicted = yhat)    # Confusion matrix

Voor onze dieren.

phat <- predict(m_pad, type = "response")
yhat <- ifelse(phat >= 0.5, 1, 0)
ctab <- table(observed = pad_spiegel$bekeek_zichzelf, predicted = yhat)
ctab

        predicted
observed  0  1
       0 66 13
       1 20 53

# Accuracy, sensitivity, specificity.
acc  <- sum(diag(ctab)) / sum(ctab)
sens <- ctab["1", "1"] / sum(ctab["1", ])   # Onder daadwerkelijke bekijkers: hoeveel correct?
spec <- ctab["0", "0"] / sum(ctab["0", ])   # Onder daadwerkelijke niet-bekijkers: hoeveel correct?

round(c(accuracy = acc, sensitivity = sens, specificity = spec), 3)

   accuracy sensitivity specificity 
      0.783       0.726       0.835

# Zet de voorspelde kans en de voorspelde klasse als kolommen naast de echte
# uitkomst. if_else() is de dplyr-broer van ifelse().
pad_klas <- pad_spiegel |>
  mutate(phat = predict(m_pad, type = "response"),
         yhat = if_else(phat >= 0.5, 1, 0))

# De classificatietabel als lijst van combinaties: echt en voorspeld.
pad_klas |>
  count(observed = bekeek_zichzelf, predicted = yhat)

  observed predicted  n
1        0         0 66
2        0         1 13
3        1         0 20
4        1         1 53

# mean(yhat == bekeek_zichzelf) is het aandeel goed. Binnen de echte
# bekijkers ([bekeek_zichzelf == 1]) is hetzelfde aandeel de sensitiviteit.
pad_klas |>
  summarise(accuracy    = mean(yhat == bekeek_zichzelf),
            sensitivity = mean(yhat[bekeek_zichzelf == 1] == 1),
            specificity = mean(yhat[bekeek_zichzelf == 0] == 0)) |>
  round(3)

  accuracy sensitivity specificity
1    0.783       0.726       0.835

Twee kanten van dezelfde tabel. De rijen zijn wat er echt gebeurde, de kolommen wat het model voorspelt.

Langs de rijen: sensitiviteit en specificiteit (in de code hierboven sensitivity en specificity). Gegeven dat een dier echt tuurt, hoe vaak zegt het model dat ook?
Langs de kolommen: PPV en NPV. Zet alle dieren waarvan het model zegt “die tuurt” bij elkaar aan de vijver en pak er één uit: hoe groot is de kans dat die ook echt tuurt? Dat is de positief voorspellende waarde (PPV). Bij de pad: $53$ van de $66$ dieren met de voorspelling tuurt, dus PPV $= .80$ .
De negatief voorspellende waarde (NPV) doet hetzelfde met de kolom tuurt niet: $66$ van de $86$ , dus NPV $= .77$ .

Is dat goed? Leg het aandeel goed naast de gok zonder model. Het null-model van het begin van dit thema zat bij $.52$ van de dieren goed, door altijd tuurt niet te zeggen; het volle model komt op $.78$ . Je ziet die vergelijking ook wel “beter dan toeval (50%)” genoemd worden. Dat is hetzelfde getal zolang de twee uitkomsten even vaak voorkomen. Bij een scheve uitkomst scheelt het veel: zijn $90$ van de $100$ dieren niet-turers, dan zit je met altijd tuurt niet zeggen bij $.90$ goed, zonder iets over de dieren te weten.

De grens van $.50$ is een afspraak, geen wet. Leg je hem lager, dan zegt het model vaker tuurt. Je mist dan minder echte turers (de sensitiviteit stijgt), maar je zegt ook vaker tuurt bij een dier dat niet tuurt (de specificiteit daalt). Hoger leggen werkt andersom. Dit oefenboek houdt de grens op $.50$ : dan kies je per dier gewoon de uitkomst die het model het waarschijnlijkst vindt.

PPV en NPV hangen ook af van hoe vaak de uitkomst voorkomt base rate. Sensitiviteit en specificiteit zijn eigenschappen van het model. PPV en NPV niet alleen: is turen zeldzaam, dan zitten er onder de dieren met de voorspelling tuurt naar verhouding veel die het niet doen, en zakt de PPV, ook bij een goed model.

Vragen 4.2 — classificatie (t–w)

t) Hoeveel dieren zijn correct geclassificeerd? Welk percentage?

u) Wat is de sensitivity? Wat zegt dat?

v) Wat is de specificity? Wat zegt dat?

w) De cutoff staat op $0.5$ . Waarom is dat een keuze, geen wet?

Antwoord 4.2 — open na je eigen poging

In gewone woorden.

Diagonaal: $66 + 53 = 119$ correct uit $152$ , accuracy $78.3\%$ .
Sensitivity $.73$ : van de $73$ daadwerkelijke bekijkers heeft het model er $53$ correct herkend ( $73\%$ ).
Specificity $.84$ : van de $79$ niet-bekijkers er $66$ correct als niet-bekijker voorspeld ( $84\%$ ) — het model is iets beter in “die kijkt zichzelf níet” zeggen dan andersom.
De cutoff van $0.5$ is een conventie, niet een wet. Als je sensitivity wilt maximaliseren (bv. een dier dat dringend zelfreflectie nodig heeft niet missen), zet je de cutoff lager. Andersom hoger. ROC-curves vergelijken alle cutoffs tegelijk — zie Wat blijft liggen.

APA-stijl.

Op basis van een cutoff van $\hat{p} \geq .50$ classificeerde het model $119$ van $152$ dieren correct (accuracy $= .78$ ). Sensitivity $= .73$ , specificity $= .84$ .

T7 — Classificatie-accuratesse + PPV/NPV

Vraag T7 — Pen-en-papier

Naast accuracy, sensitivity en specificity zijn ook PPV en NPV relevant (zie Twee kanten van dezelfde tabel hierboven). Stel je leest:

        predicted
observed  0  1
       0 60 20
       1 15 55

a) Bereken accuracy, sensitivity en specificity.

b) PPV $= P(Y=1 \mid \hat{Y}=1)$ : van degenen voor wie het model $1$ voorspelt, welk deel heeft echt $Y = 1$ ? Bereken.

c) NPV $= P(Y=0 \mid \hat{Y}=0)$ : van degenen voor wie het model $0$ voorspelt, welk deel heeft echt $Y = 0$ ? Bereken.

d) Waarom hangt PPV/NPV af van de base rate (het aandeel met $Y = 1$ in de populatie)?

Antwoord T7 — open na je eigen poging

a) Accuracy = $(60+55)/150 = .767$ . Sensitivity = $55/(15+55) = .786$ . Specificity = $60/(60+20) = .750$ .

b) PPV = $55/(20+55) = .733$ . Van degenen voor wie het model $1$ voorspelt, klopt het bij $73\%$ .

c) NPV = $60/(60+15) = .800$ . Van degenen voor wie het model $0$ voorspelt, klopt het bij $80\%$ .

d) Sensitivity en specificity zijn eigenschappen van het model — onafhankelijk van hoe vaak $Y=1$ in de populatie voorkomt. PPV en NPV daarentegen hangen af van de base rate. In een populatie waar $Y=1$ heel zeldzaam is, zal PPV laag zijn — zelfs een goed model zal veel false positives geven, simpelweg omdat er weinig echte positives zijn om voorspeld te worden. Vandaar dat een test-met-hoge-sensitivity in een lage-prevalentie-context (bv. zeldzame ziekte) toch tot veel valse alarmen kan leiden.

4.2.h APA-zin — alles in één

Toon code (de getallen voor de APA-zin ophalen)

# Een paar getallen rechtstreeks ophalen voor de zin hieronder.
chi <- anova(m_null, m_pad, test = "Chisq")[2, "Deviance"]
df_full <- m_pad$df.null - m_pad$df.residual
p_full  <- anova(m_null, m_pad, test = "Chisq")[2, "Pr(>Chi)"]

cat("LR overall: chi^2(", df_full, ") = ",
    round(chi, 2), ", p = ", signif(p_full, 3),
    ", pseudo-R^2 = ", round(1 - m_pad$deviance/m_pad$null.deviance, 3),
    sep = "")

LR overall: chi^2(4) = 66.22, p = 0.000000000000143, pseudo-R^2 = 0.315

Een nette APA-rapportage

Een logistische regressie van zelf-bekijken op leeftijd, aantal afwijzingen en diersoort verklaarde de uitkomst significant beter dan een null-model, $\chi^2_{\text{LR}}(4) = 66.22$ , p < .001, pseudo- $R^2_L = .31$ . Per jaar dat een dier ouder was, lagen de odds van zelf-bekijken $0.86$ keer zo hoog ( $95\%$ -CI $[0.76, 0.98]$ , p = .028); per extra afwijzing waren de odds $1.34$ keer hoger ( $[1.03, 1.74]$ , p = .032). Diersoort droeg sterk bij ( $\chi^2_{\text{LR}}(2) = 58.10$ , p < .001): vossen hadden odds van $10.49$ ten opzichte van reeën, dassen slechts $0.18$ . Op basis van een cutoff $\hat{p} \geq .50$ classificeerde het model $78\%$ van de dieren correct (sensitivity $= .73$ , specificity $= .84$ ).

Voor wie verder wil — interactie tussen leeftijd en afwijzingen

“Misschien,” zei de pad. “Misschien dat een afgewezen oud dier net zo hard in het water staart als een afgewezen jong dier. Dan is leeftijd niet meer interessant. Dan is afwijzing alles.”

Een logistisch model kun je net als een gewone regressie uitbreiden met een interactieterm. Een interactie tussen leeftijd en aantal afwijzingen zegt: het effect van leeftijd op de log-odds van turen verandert in grootte of richting met het aantal afwijzingen.

m_int <- glm(bekeek_zichzelf ~ leeftijd * aantal_afwijzingen + diersoort,
             data = pad_spiegel,
             family = binomial(link = "logit"))

# Vergelijk met het model zonder interactie: voegt de interactie iets toe?
anova(m_pad, m_int, test = "Chisq")

Analysis of Deviance Table

Model 1: bekeek_zichzelf ~ leeftijd + aantal_afwijzingen + diersoort
Model 2: bekeek_zichzelf ~ leeftijd * aantal_afwijzingen + diersoort
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       147     144.26                     
2       146     143.95  1  0.31881   0.5723

Wat zegt de output? De LR-toets vergelijkt het model met en zonder de interactieterm: $\chi^2_{\text{LR}}(1) = 0.32$ , $p = .57$ . Geen bewijs dus dat het effect van leeftijd op de log-odds van turen afhangt van het aantal afwijzingen, of omgekeerd. De pad blijft bij het model zonder interactie.

Was de interactie wél significant geweest, dan rapporteer je het model mét interactie en lees je de losse effecten voorwaardelijk, net als in MRA en ANOVA:

Coëfficiënten: $b_{\text{leeftijd}}$ is dan het effect van leeftijd op de log-odds van turen bij nul afwijzingen, en niet meer “het” effect van leeftijd.
Odds ratios van die losse effecten gelden dan ook alleen bij die ene waarde van de andere voorspeller.
Vrijheidsgraden: de interactieterm van twee continue voorspellers kost er één.
Toetsen doe je de interactie met de LR-toets: model met tegen model zonder.
Centreren van de continue voorspellers helpt bij het lezen: de losse effecten gelden dan bij het gemiddelde van de andere voorspeller in plaats van bij nul. Thema 3 legt centreren uit.

4.3 Aannamechecks

“Lopen mijn voorspellers wel mooi op de log-odds-schaal?”

LRA leunt op andere aannames dan MRA: vier aannames, plus genoeg dieren.

Log-odds-lineariteit — continue voorspellers moeten lineair samenhangen met de logit van $Y$ (zie Log-odds-lineariteit hieronder).
Onafhankelijke observaties — geen geclusterde data, geen herhaalde metingen op hetzelfde dier.
Voorspellers zonder meetfout — meetfouten in $X$ vertekenen $b$ .
Geen sterke multicollineariteit — de VIF-check (zie Multicollineariteit hieronder).

Plus genoeg dieren per voorspeller: $N/k \geq 30$ (zie Sample-size hieronder).

Onafhankelijkheid en meetfout check je niet in R: die volgen uit hoe de data verzameld zijn (elk dier één keer geteld, een betrouwbaar meetinstrument). Twee aannames die je van MRA kent, vallen hier juist weg: normaliteit en homoscedasticiteit.

Normaliteit en homoscedasticiteit hoeven NIET (klassieke trick-question)

Op tentamens komt het terug: “Welke aannames moet je checken bij LRA?” en dan staat in de antwoord-opties iets als “normaal verdeelde residuen” en “gelijke varianties tussen groepen”. Dat zijn de afleiders. Wie reflexmatig de MRA- of ANOVA-aannames opdreunt, valt erin.

Waarom werken die aannames hier niet?

$Y$ is binair. Een binaire variabele is per definitie binomiaal verdeeld, niet normaal. En de variantie van $Y$ hangt vast aan de kans $P$ dat $Y = 1$ : $\text{Var}(Y) = P(1 - P)$ . Een $P = .5$ geeft variantie $.25$ ; een $P = .9$ geeft variantie $.09$ . Heteroscedasticiteit is dus geen schending — het is de structuur. Daarom ook geen $F$ -toets bij LRA, en daarom ook geen Levene of Shapiro-Wilk.

Tentamen-veiligheid: als je bij een LRA-vraag een aanname-checklist moet opnoemen, bouw hem op uit de vier hierboven, plus $N/k \geq 30$ , niet uit de MRA-aannames. Wie “normaliteit” of “homoscedasticiteit” in zijn LRA-antwoord opneemt, verliest punten.

4.3.a Log-odds-lineariteit

Let op wát er recht moet lopen: niet de kans, want die volgt altijd een S-curve, maar de log-odds. En de aanname geldt alleen voor continue voorspellers. Een factor als diersoort gaat via dummies, met een eigen coëfficiënt per niveau; daar is geen lijn die recht hoeft te zijn.

Algemene vorm. Je bekijkt het: deel een continue voorspeller op in een paar groepjes (bins), reken per groepje de log-odds van $Y = 1$ uit, en zet die uit tegen $X$ . Liggen de punten ongeveer op een rechte lijn, dan houdt de aanname. Is de samenhang duidelijk krom, dan kun je een kwadratische term toevoegen (I(X^2) in de formule) of de voorspeller in een paar groepen indelen. Dat laatste is eenvoudig, maar gooit informatie weg.

Voor onze dieren.

Toon code (visuele log-odds-check)

# Bin de continue voorspellers en bereken de empirische log-odds per bin.
bin_logit_plot <- function(x, y, k = 8, xlab = deparse(substitute(x))) {
  qs <- quantile(x, probs = seq(0, 1, length.out = k + 1), na.rm = TRUE)
  # Maak unieke breaks om bins te garanderen.
  qs <- unique(qs)
  bins <- cut(x, breaks = qs, include.lowest = TRUE)
  agg <- aggregate(y, by = list(bin = bins), FUN = mean)
  # Bereken bin-midpoint
  mids <- aggregate(x, by = list(bin = bins), FUN = mean)$x
  ps <- pmin(pmax(agg$x, 0.01), 0.99)
  logits <- log(ps / (1 - ps))
  plot(mids, logits,
       pch = 19, col = "#0077BB",
       xlab = xlab, ylab = "Empirische log-odds",
       main = paste("Log-odds-lineariteit:", xlab))
  abline(lm(logits ~ mids), lty = 2, col = "#EE7733")
}

par(mfrow = c(1, 2), mar = c(4, 4, 2.5, 1))
bin_logit_plot(pad_spiegel$leeftijd, pad_spiegel$bekeek_zichzelf,
               k = 6, xlab = "Leeftijd")
bin_logit_plot(pad_spiegel$aantal_afwijzingen, pad_spiegel$bekeek_zichzelf,
               k = 5, xlab = "Aantal afwijzingen")

Toon code (visuele log-odds-check)

par(mfrow = c(1, 1))

Vragen 4.3

a) Lopen de empirische log-odds redelijk lineair met leeftijd?

b) Idem met aantal_afwijzingen — let op het bereik aan de rechterkant.

c) Wat zou je doen als één van de twee duidelijk gekromd is?

Antwoord 4.3 — open na je eigen poging

In gewone woorden.

Voor leeftijd lopen de bin-punten redelijk dalend: jongere dieren hogere log-odds, oudere lagere — een rechte daling past goed.
Voor aantal_afwijzingen iets meer wisselvallig — bij hoge afwijzings-bins is de schatting onstabiel doordat er minder dieren zijn, maar de gemiddelde trend is nog stijgend.
Bij duidelijke kromming: opties zijn (1) een kwadratische term toevoegen (I(X^2)); (2) categoriseren in groepen; (3) een niet-lineaire link gebruiken (cubic spline, GAM). Categoriseren is laagdrempelig maar verliest informatie; spline/GAM houdt continuïteit (zie Wat blijft liggen).

APA-stijl.

Visuele inspectie van empirische log-odds tegen de continue voorspellers ondersteunde de aanname van log-odds-lineariteit; afwijkingen waren beperkt tot bins met weinig dieren.

Voor wie verder wil — de aanname toetsen in plaats van bekijken

Naast het plaatje bestaat er een toets, de Box-Tidwell-check: zet per continue voorspeller een extra term $X \cdot \log(X)$ in het model. Is die term significant, dan buigt de samenhang tussen $X$ en de log-odds af van een rechte lijn.

m_check <- glm(Y ~ X1 + X1:log(X1) + ..., data = mijn_data,
               family = binomial(link = "logit"))
summary(m_check)

Hij werkt alleen voor een $X$ met uitsluitend positieve waarden, want $\log(0)$ bestaat niet. Bij de pad kan het dus niet zomaar voor aantal_afwijzingen: daar komen dieren met nul afwijzingen in voor.

4.3.b Multicollineariteit — `vif()`

Dezelfde check als bij MRA: een VIF boven de $10$ betekent dat een voorspeller te sterk samenhangt met de andere voorspellers. Alleen ziet de uitvoer er hier anders uit. Omdat diersoort een factor is, geeft vif() geen rijtje getallen maar een tabel met drie kolommen. Lees de eerste, GVIF (de gegeneraliseerde VIF). Voor een voorspeller met één vrijheidsgraad (Df is $1$ ), zoals leeftijd en aantal_afwijzingen, is dat gewoon de VIF. Ligt alles in die kolom onder de $10$ , dan ben je klaar, ook voor de factor.

vif(m_pad)

                       GVIF Df GVIF^(1/(2*Df))
leeftijd           1.032157  1        1.015951
aantal_afwijzingen 1.050872  1        1.025120
diersoort          1.075163  2        1.018283

Vragen 4.3

d) Wat zijn de VIF’s? Liggen ze onder de grens van $10$ ?

Antwoord 4.3 — open na je eigen poging

In gewone woorden.

Leeftijd en aantal_afwijzingen hebben elk één vrijheidsgraad, dus hun getal in de eerste kolom ís de VIF: $1.03$ en $1.05$ . Diersoort komt op $1.08$ . Alles ligt rond de $1$ , ruim onder de grens van $10$ . Geen multicollineariteit-zorg.

APA-stijl.

Multicollineariteit was geen probleem: de VIF’s van leeftijd en aantal afwijzingen waren $1.03$ en $1.05$ , de gegeneraliseerde VIF van diersoort $1.08$ , alle ruim onder $10$ .

Voor wie verder wil — de GVIF bij een factor

Voor een factor rekent R een gegeneraliseerde VIF (GVIF) uit over alle dummies samen; de kolom Df zegt over hoeveel dummies. Ligt die GVIF onder de $10$ , dan is er niets aan de hand. Komt hij erboven, dan is dat bij een factor nog geen alarm: een GVIF groeit mee met het aantal dummies. Kijk dan naar de derde kolom, $\text{GVIF}^{1/(2 \cdot Df)}$ . Die is te vergelijken met de wortel uit een gewone VIF, dus de grens van $10$ wordt daar $\sqrt{10} \approx 3.16$ . Bij de pad komt die kolom niet boven de $1.03$ .

4.3.c Sample-size — Pedhazur’s $N/k \geq 30$

N_obs   <- nrow(pad_spiegel)
k_terms <- length(coef(m_pad)) - 1   # exclusief intercept
ratio   <- N_obs / k_terms
round(c(N = N_obs, k = k_terms, ratio = ratio), 2)

    N     k ratio 
  152     4    38

Pedhazur’s drempel — strenger dan MRA

Bij MRA hanteren we vaak $N/k \geq 20$ als sample-size-vuistregel. Voor LRA is de aanbeveling strenger: $N/k \geq 30$ (Pedhazur, 1997). Reden: een binaire $Y$ bevat minder informatie per observatie dan een continue $Y$ .

Wat telt als $k$ ? Het aantal coëfficiënten naast het intercept. Een factor telt dus voor zijn aantal dummies: diersoort heeft $3$ niveaus en dus $2$ dummies. Voor de pad is dat $k = 1 + 1 + 2 = 4$ , voor drie voorspellers. In veel opgaven heeft elke voorspeller één coëfficiënt (hij is interval of 0/1); dan zijn de twee tellingen hetzelfde, en tel je gewoon de voorspellers. Zit er toch een factor met meer dan twee niveaus in, tel dan zijn dummies.

Voor onze pad: $N/k = 152/4 = 38$ — boven Pedhazur’s drempel.

Haal je de drempel niet, dan zijn er twee nette uitwegen: minder voorspellers (de inhoudelijk belangrijkste) of meer data.

Voor wie verder wil — tellen in uitkomsten in plaats van in dieren

Een tweede vuistregel kijkt niet naar het totale aantal waarnemingen, maar naar de kleinste van de twee uitkomstgroepen: het aantal events per voorspeller EPV, met meestal $10$ als ondergrens. Bij een scheve uitkomst is dat strenger dan $N/k$ . Neem $160$ dieren waarvan er maar $12$ turen, en $4$ voorspellers: $N/k = 40$ , ruim boven de $30$ , maar EPV $= 12/4 = 3$ , ver onder de $10$ . De schattingen worden dan wankel en de standaardfouten groot. Bij de pad is de kleinste groep $73$ dieren, dus EPV $= 73/4 = 18.25$ .

Voor modellen met veel voorspellers en weinig uitkomsten bestaat naast minder voorspellers of meer data nog een derde weg: penalized regression, zoals de lasso, die de gewichten bewust kleiner schat.

Voor wie verder wil — invloedrijke dieren opsporen met Cook’s distance

Net als bij MRA kun je per dier uitrekenen hoeveel het model verschuift als je dat ene dier weglaat: Cook’s distance.

plot(m_pad, which = 4, sub.caption = "",
     caption = "Cook's distance per observatie")

De hoogste waarde bij de pad is $0.10$ . Een gangbare grens voor “kijk hier eens naar” is $1$ ; daar komt geen enkel dier in de buurt. Zie je wel een twijfelgeval, laat het dan niet stilletjes weg: draai het model met en zonder dat dier, en meld of de conclusies veranderen.

T8 — Log-odds-lineariteit: wat is het, hoe checken?

Vraag T8 — Pen-en-papier

a) Wat zegt de aanname van log-odds-lineariteit precies? Voor welk type voorspeller geldt hij?

b) Hoe check je hem?

c) Stel je vindt dat een continue voorspeller $X$ niet-lineair samenhangt met de log-odds — een U-vormige relatie. Welke aanpassing aan het model zou dat oplossen?

Antwoord T8 — open na je eigen poging

a) De relatie tussen elke continue voorspeller $X_j$ en de logit van $Y$ moet lineair zijn: $\text{logit}(P(Y=1)) = b_0 + b_j X_j + \ldots$ . Voor factoren geldt de aanname niet — dummies coderen automatisch elk niveau apart, dus geen lineariteit nodig. Alleen continue voorspellers.

b) Visueel: deel $X$ op in groepjes, bereken per groepje de empirische logit, zet die uit tegen $X$ en kijk of het ongeveer een rechte lijn is.

c) Voeg een kwadratische term toe: $\text{logit}(P) = b_0 + b_1 X + b_2 X^2 + \ldots$ Of categoriseer $X$ in groepen. Of gebruik een spline/GAM (zie Wat blijft liggen).

T9 — $N/k \geq 30$ : vergelijken met MRA

Vraag T9 — Pen-en-papier

a) Bij MRA hanteert het oefenboek $N/k \geq 20$ ; bij LRA $N/k \geq 30$ . Waarom is de LRA-drempel strenger?

b) Een onderzoeker heeft $N = 80$ en wil $4$ voorspellers in een logistisch model stoppen. $N/k = 20$ . Voldoet dit aan Pedhazur’s vuistregel?

c) Wat is een gangbare oplossing als je sample-size te klein is voor de gewenste voorspeller-set?

Antwoord T9 — open na je eigen poging

a) Een binaire $Y$ bevat per observatie maar één bit informatie. Een continue $Y$ veel meer (zo veel bits als de meetschaal toelaat). De maximum-likelihood-schatter heeft dus meer observaties nodig om dezelfde precisie te bereiken voor logistische coëfficiënten als OLS voor MRA-coëfficiënten.

b) Nee — Pedhazur’s drempel is $30$ per voorspeller, hier zou je $N \geq 120$ willen hebben voor $4$ voorspellers.

c) Verminder het aantal voorspellers: kies de inhoudelijk belangrijkste. Of verzamel meer data.

Inzicht. Sample-size-vuistregels zijn vuistregels, geen wetten. Ze zeggen wanneer je schattingen wankel worden, niet wanneer ze ineens fout zijn.

T10 — Niet-significante voorspeller weglaten?

Een voorspeller die niet significant is, haal je niet zomaar uit het model. Drie dingen om te weten voordat je schrapt:

De dummies van een factor horen bij elkaar. Haal je één dummy weg, dan voeg je dat niveau stilletjes samen met de referentie. Een factor laat je dus heel staan of je haalt hem heel weg. Of hij als geheel iets doet, zegt de LR-toets op de hele factor.
Niet significant is niet hetzelfde als onbelangrijk. Hangt een voorspeller sterk samen met een andere voorspeller in het model, dan kan zijn eigen toets niet significant uitvallen terwijl hij samen met die ander wel iets bijdraagt. Kijk eerst naar de VIF’s.
Kiezen op de $p$ -waarde, achteraf, maakt vals alarm waarschijnlijker. Wie voorspellers weglaat of toevoegt tot het mooi uitkomt, met de hand of met een automatische stapsgewijze procedure, vindt vaker een effect dat er in de populatie niet is. Kies je voorspellers vooraf, op inhoud, en meld wat je allemaal geprobeerd hebt. Wil je weten wat één voorspeller toevoegt bovenop de rest, zet hem dan als laatste blok in een hiërarchische LRA.

Vraag T10 — Pen-en-papier

Stel: in een uitgebreidere versie van de pad-data is een vierde voorspeller seizoen (factor: lente/zomer/herfst/winter) toegevoegd, en blijkt seizoenherfst niet significant ( $p = .41$ ) terwijl andere niveaus van seizoen dat wel zijn.

a) Mag je die ene dummy weglaten uit het model? Wat zou er met het model gebeuren?

b) Wat heeft multicollineariteit met deze keuze te maken?

c) Wat zegt de LR-toets voor de hele factor seizoen? Is die uitslag relevant voor het wel/niet weglaten?

d) Algemenere vraag: wanneer is “weglaten op basis van $p$ -waarde” een verdedigbare strategie, en wanneer niet?

Antwoord T10 — open na je eigen poging

a) Je kunt een dummy niet zomaar weglaten: dummies bij een factor horen samen — als je seizoenherfst schrapt, codeer je impliciet herfst samen met de referentie (bv. lente), wat inhoudelijk niet klopt. Wel kun je de hele factor weglaten, of niveaus samenvoegen op theoretische gronden.

b) Bij MRA is “een niet-significante voorspeller weglaten” verdedigbaar als hij geen unieke variantie verklaart. Maar bij sterke multicollineariteit kan hij niet-significant lijken terwijl hij wél bijdraagt — alleen samen met een collineaire collega. Eerst VIF checken, dan beslissen.

c) Als de LR-toets op de hele factor significant is (bv. $\chi^2_{\text{LR}}(3) = 9.20$ , $p = .027$ ), dan zit er als geheel signaal in seizoen, ook al ziet één dummy er los niet-significant uit. Sterk argument om de hele factor in het model te houden.

d) Wel verdedigbaar: als je vooraf een datagestuurde model-selectie-strategie (bv. backward elimination) hebt aangekondigd, of als de inhoud zegt “deze voorspeller hoorde er sowieso niet thuis”. Niet verdedigbaar: als je achteraf gaat shoppen tot het mooi uitkomt — dat is HARKing. Hou je vooraf vast aan een theoretische voorspeller-set en rapporteer alles wat je probeerde.

Inzicht. Een $p$ -waarde is geen weeg-maat voor “model-fitness van een voorspeller”. Het is een toetsingsuitkomst onder veronderstelde nul-hypothese. Modelreductie verdient een doordachte strategie, niet ad-hoc-knippen.

R-spiekblad bij LRA

Alle commando’s op één plek

Pakketten en data laden

# Pakketten activeren — één keer per sessie.
library(car)         # voor vif() en Anova(., test = "LR")
library(dplyr)       # alleen voor de dplyr-tabbladen
library(tidyr)       # voor pivot_wider() in de dplyr-kruistabel

# Hoofd-dataset laden — object 'pad_spiegel' verschijnt vanzelf.
load("data/pad_spiegel.RData")
str(pad_spiegel)

# Mini-dataset voor E5.
load("data/kerkuil_nacht.RData")

Verkenning — frequenties en boxplots

# Aantallen per uitkomst-categorie en proportie.
table(pad_spiegel$bekeek_zichzelf)
prop.table(table(pad_spiegel$bekeek_zichzelf))

# Kruistabel factor x Y.
table(pad_spiegel$diersoort, pad_spiegel$bekeek_zichzelf)

# Aantallen per uitkomst-categorie en proportie.
pad_spiegel |>
  count(bekeek_zichzelf) |>
  mutate(proportie = n / sum(n))

# Kruistabel factor x Y.
pad_spiegel |>
  count(diersoort, bekeek_zichzelf) |>
  pivot_wider(names_from = bekeek_zichzelf, values_from = n)

# Boxplot continue X per uitkomst.
boxplot(leeftijd           ~ bekeek_zichzelf, data = pad_spiegel)
boxplot(aantal_afwijzingen ~ bekeek_zichzelf, data = pad_spiegel)

Model fitten — `glm(..., family = binomial(link = "logit"))`

# Het LRA-model — let op family-argument.
m_pad <- glm(bekeek_zichzelf ~ leeftijd + aantal_afwijzingen + diersoort,
             data   = pad_spiegel,
             family = binomial(link = "logit"))

# Coëfficiënten en Wald-z per coëfficiënt.
summary(m_pad)

Coëfficiënten en odds ratios

# Log-odds-coëfficiënten.
coef(m_pad)

# Odds ratios.
exp(coef(m_pad))

# 95%-CI op log-odds-schaal.
confint.default(m_pad)

# 95%-CI op odds-schaal.
exp(confint.default(m_pad))

Modelvergelijking via likelihood ratio

# Null-model (alleen intercept).
m_null <- glm(bekeek_zichzelf ~ 1,
              data = pad_spiegel,
              family = binomial(link = "logit"))

# Volledig model versus null.
anova(m_null, m_pad, test = "Chisq")

# LR-toets per term (factor als geheel).
car::Anova(m_pad, type = "II", test.statistic = "LR")

Pseudo- $R^2$ van Hosmer en Lemeshow

# Via log-likelihoods.
LL_null <- as.numeric(logLik(m_null))
LL_full <- as.numeric(logLik(m_pad))
1 - LL_full / LL_null

# Equivalent via deviances.
1 - m_pad$deviance / m_pad$null.deviance

Hiërarchische LRA — blok voor blok

# Elk model één blok groter dan het vorige, op dezelfde data.
m_h0 <- glm(bekeek_zichzelf ~ 1, data = pad_spiegel, family = binomial(link = "logit"))
m_h1 <- glm(bekeek_zichzelf ~ diersoort, data = pad_spiegel, family = binomial(link = "logit"))
m_h2 <- glm(bekeek_zichzelf ~ diersoort + leeftijd, data = pad_spiegel, family = binomial(link = "logit"))
m_h3 <- glm(bekeek_zichzelf ~ diersoort + leeftijd + aantal_afwijzingen,
            data = pad_spiegel, family = binomial(link = "logit"))

# Per blok: chi2 (kolom Deviance), df en p.
anova(m_h0, m_h1, m_h2, m_h3, test = "Chisq")

# R2_L per model.
r2_L <- function(m) (m$null.deviance - m$deviance) / m$null.deviance
sapply(list(m_h1, m_h2, m_h3), r2_L)

Classificatietabel — predicted vs observed

# Voorspelde kansen.
phat <- predict(m_pad, type = "response")

# Cutoff op 0.5.
yhat <- ifelse(phat >= 0.5, 1, 0)

# Confusion matrix.
ctab <- table(observed = pad_spiegel$bekeek_zichzelf, predicted = yhat)
ctab

# Accuracy, sensitivity, specificity.
sum(diag(ctab)) / sum(ctab)
ctab["1", "1"] / sum(ctab["1", ])
ctab["0", "0"] / sum(ctab["0", ])

# Voorspelde kans en klasse (cutoff 0.5) als kolommen erbij.
pad_klas <- pad_spiegel |>
  mutate(phat = predict(m_pad, type = "response"),
         yhat = if_else(phat >= 0.5, 1, 0))

# Confusion matrix, als lijst van combinaties.
pad_klas |>
  count(observed = bekeek_zichzelf, predicted = yhat)

# Accuracy, sensitivity, specificity.
pad_klas |>
  summarise(accuracy    = mean(yhat == bekeek_zichzelf),
            sensitivity = mean(yhat[bekeek_zichzelf == 1] == 1),
            specificity = mean(yhat[bekeek_zichzelf == 0] == 0))

Aannamechecks — vif en log-odds-lineariteit

# Multicollineariteit: alle VIF's onder de 10?
vif(m_pad)

# Log-odds-lineariteit: bin de continue X, plot empirische logit.
qs   <- quantile(pad_spiegel$leeftijd, probs = seq(0, 1, length.out = 7))
bins <- cut(pad_spiegel$leeftijd, breaks = qs, include.lowest = TRUE)
agg  <- aggregate(bekeek_zichzelf ~ bins, data = pad_spiegel, FUN = mean)
agg$logit <- log(pmin(pmax(agg$bekeek_zichzelf, 0.01), 0.99) /
                 (1 - pmin(pmax(agg$bekeek_zichzelf, 0.01), 0.99)))
plot(seq_along(agg$logit), agg$logit, type = "b")

# Log-odds-lineariteit: per groepje leeftijd het aandeel turers, en dat
# omgezet naar log-odds. group_by() + summarise() = per groep samenvatten.
per_bin <- pad_spiegel |>
  mutate(bin = cut(leeftijd,
                   breaks = quantile(leeftijd, probs = seq(0, 1, length.out = 7)),
                   include.lowest = TRUE)) |>
  group_by(bin) |>
  summarise(p = mean(bekeek_zichzelf)) |>
  mutate(p_begrensd = pmin(pmax(p, 0.01), 0.99),
         logit      = log(p_begrensd / (1 - p_begrensd)))
plot(seq_along(per_bin$logit), per_bin$logit, type = "b")

Voorbeeld-tentamenvragen

Even oefenen op tentamen-toon

Onderaan dit hoofdstuk staan vier korte theorievragen in tentamen-stijl en één mini R-opdracht. Geen verhaal-frame meer: hier oefen je zoals op het tentamen. Hou de vuistregels paraat: $\mathit{OR} = \exp(b)$ als vermenigvuldigingsfactor, het interval $\exp(b \pm 1.96 \cdot SE)$ , $R^2_L$ uit twee deviances, per blok de $\chi^2$ van het verschil in deviance met $df$ = het aantal coëfficiënten dat erbij komt, sensitiviteit en specificiteit langs de rijen en PPV en NPV langs de kolommen, VIF onder $10$ , $N/k \geq 30$ , Wald versus LR. Het hele lijstje staat bovenaan het thema, onder Wat je na dit thema moet kunnen.

Theorie en handreken

Vraag E1 — OR aflezen en interpreteren

Een vleermuis registreert per nachtelijke vlucht of er prooi gevangen is (succes ja/nee), met als voorspellers maan-helderheid en windsnelheid. Hij rapporteert:

                Estimate Std.Error z value Pr(>|z|)
(Intercept)       0.812     0.620    1.31    .190
maan_helderheid  -2.140     0.820   -2.61    .009
windsnelheid     -0.085     0.041   -2.07    .038

Wat is de odds ratio voor maan-helderheid, en hoe interpreteer je hem?

$\mathit{OR} = -2.14$ ; per eenheid maan-helderheid daalt de kans op succes met $2.14$ .
$\mathit{OR} = \exp(-2.14) \approx 0.12$ ; per eenheid maan-helderheid (volle maan = 1) worden de odds van succes met factor $0.12$ vermenigvuldigd — donker maakt jagen veel succesvoller.
$\mathit{OR} = \exp(2.14) \approx 8.50$ ; volle maan verhoogt de kans op succes drastisch.
$\mathit{OR} = 1 + (-2.14) = -1.14$ ; geen geldige interpretatie.

Antwoord E1 — open na je eigen poging

b) $\mathit{OR} = \exp(b) = \exp(-2.14) \approx 0.118$ . Negatief $b$ betekent OR onder $1$ — donkerder geeft meer succes. Optie a verwart $b$ met OR, c heeft het teken verkeerd, d is geen geldige formule.

Vraag E2 —

b

versus

\exp(b)

uit summary

Een marterspecht meet of een nestvogel een ei laat liggen na verstoring (laat-vallen ja/nee), met als voorspeller “verstoring-intensiteit” (continu, schaal 0-10):

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)
(Intercept)              -3.20      0.85    -3.76    < .001
verstoring_intensiteit    0.62      0.18     3.44    < .001

Welke uitspraak is correct?

Per eenheid verstoring stijgt de kans op laat-vallen met $0.62$ .
Per eenheid verstoring stijgt de log-odds met $0.62$ ; de odds vermenigvuldigen met $\exp(0.62) \approx 1.86$ .
De odds ratio is $0.62$ — een effect onder $1$ , dus verstoring verlaagt de kans.
De $p$ -waarde is significant; daarom is de OR per definitie $1$ .

Antwoord E2 — open na je eigen poging

b) $b = 0.62$ is op de log-odds-schaal; de OR is $\exp(0.62) \approx 1.86$ . Optie a verwart kans met odds, c verwart $b$ met OR (negatief $b$ zou onder $1$ geven, niet $0.62$ ), d zegt iets onzinnigs.

Vraag E3 — Classificatie-tabel beoordelen

Een uil bouwt een logistisch model voor jacht-succes. De confusion matrix:

        predicted
observed   0   1
       0  72  18
       1  20  60

Wat is de accuracy en de sensitivity?

Accuracy $= .78$ ; sensitivity $= .75$ .
Accuracy $= .77$ ; sensitivity $= .80$ .
Accuracy $= .78$ ; sensitivity $= .60$ .
Accuracy $= .60$ ; sensitivity $= .80$ .

Antwoord E3 — open na je eigen poging

a) Accuracy $= (72 + 60) / 170 = 132/170 = .776 \approx .78$ . Sensitivity $= 60/(20 + 60) = 60/80 = .750 \approx .75$ . Optie b verwart sensitivity met specificity, c heeft sensitivity verkeerd, d wisselt accuracy en sensitivity om.

Vraag E4 — Pseudo-

R^2

aflezen en interpreteren

Een onderzoek-coördinator runt een logistisch model en rapporteert:

Null deviance: 220.4 on 159 degrees of freedom
Residual deviance: 168.6 on 154 degrees of freedom

Welke uitspraak past het best?

De pseudo- $R^2_L = 168.6/220.4 \approx .76$ ; een uitzonderlijk slecht model.
De pseudo- $R^2_L = (220.4 - 168.6)/220.4 \approx .235$ ; het model haalt zo’n $24\%$ van de deviance van het null-model weg.
De pseudo- $R^2_L = 168.6/(220.4 - 168.6) \approx 3.25$ ; geen geldige waarde.
Pseudo- $R^2_L$ kan niet berekend worden uit deviances; je hebt de log-likelihoods nodig.

Antwoord E4 — open na je eigen poging

b) $R^2_L = 1 - D_{\text{res}}/D_{\text{null}} = 1 - 168.6/220.4 = 1 - 0.765 = 0.235$ . Equivalent: $(220.4 - 168.6)/220.4 = 51.8/220.4 = .235$ . Het model haalt dus ruim een vijfde van de deviance van het null-model weg. Optie a verwart de breuk-kant, c is rekenkundig onzinnig, d is fout (deviance $= -2LL$ , dus pseudo- $R^2$ valt direct uit deviances af te leiden).

R-practical opdrachtje

De nacht van de kerkuil

Vraag E5 — Mini-LRA bij de kerkuil

Een ecoloog onderzoekt het nachtelijke gedrag van $118$ verschillende dieren in het revier van een kerkuil — elk dier op één eigen nacht geobserveerd. Per dier wordt genoteerd of het onrustig sliep (sliep_onrustig, $0$ / $1$ ) en drie voorspellers: aantal_zorgen (continu count, hoeveel zorgen het dier de afgelopen week had), eet_voor_slapen (continu, hapjes-frequentie ’s avonds, schaal $0$ - $5$ ) en sociale_isolatie (factor: nee / ja / twijfelend). De dataset staat in data/kerkuil_nacht.RData en bevat het object kerkuil_nacht. Sla je R-commando’s op in één scriptbestand: kerkuil.R. Gebruik $\alpha = .05$ .

a) Fit het volledige LRA-model met glm(sliep_onrustig ~ aantal_zorgen + eet_voor_slapen + sociale_isolatie, family = binomial(link = "logit")). Welke voorspellers zijn significant volgens de Wald-toets? Rapporteer de toetsstatistiek en de p-waarde voor aantal_zorgen en eet_voor_slapen.

b) Bereken de odds ratios met $\exp(\text{coef}(\cdot))$ en hun $95\%$ -CI. Wat is de OR voor aantal_zorgen, en wat zegt de richting?

c) Toets de hele factor sociale_isolatie met car::Anova(model, type = "II", test.statistic = "LR"). Is sociale_isolatie als geheel significant? Rapporteer $\chi^2_{\text{LR}}$ , $df$ en $p$ .

d) Bereken de pseudo- $R^2$ via $1 - D_{\text{res}}/D_{\text{null}}$ en de classificatie-accuratesse via een tabel met cutoff $\hat{p} \geq 0.5$ . Schrijf één APA-zin.

e) Bouw het model nu in drie blokken op: eerst sociale_isolatie, dan aantal_zorgen, dan eet_voor_slapen. Voegt elk blok significant iets toe? Rapporteer per blok $\chi^2$ , $df$ en $p$ , en per model $R^2_L$ .

Antwoord E5 — open na je eigen poging

load("data/kerkuil_nacht.RData")
str(kerkuil_nacht)

'data.frame':   118 obs. of  4 variables:
 $ sliep_onrustig  : int  0 1 1 1 0 1 1 0 0 0 ...
 $ aantal_zorgen   : int  1 8 5 5 2 2 2 2 4 2 ...
 $ eet_voor_slapen : num  4.8 2.6 1.4 2 0.3 0.7 3.2 1.3 3 1.4 ...
 $ sociale_isolatie: Factor w/ 3 levels "nee","ja","twijfelend": 1 3 3 1 3 3 2 1 1 3 ...

m_uil <- glm(sliep_onrustig ~ aantal_zorgen + eet_voor_slapen + sociale_isolatie,
             data   = kerkuil_nacht,
             family = binomial(link = "logit"))
summary(m_uil)


Call:
glm(formula = sliep_onrustig ~ aantal_zorgen + eet_voor_slapen + 
    sociale_isolatie, family = binomial(link = "logit"), data = kerkuil_nacht)

Coefficients:
                           Estimate Std. Error z value    Pr(>|z|)    
(Intercept)                 -5.1464     1.0271  -5.011 0.000000542 ***
aantal_zorgen                1.0349     0.2100   4.929 0.000000826 ***
eet_voor_slapen              0.6367     0.1936   3.290     0.00100 ** 
sociale_isolatieja           2.5115     0.6730   3.732     0.00019 ***
sociale_isolatietwijfelend   0.8619     0.5667   1.521     0.12827    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 161.41  on 117  degrees of freedom
Residual deviance: 108.51  on 113  degrees of freedom
AIC: 118.51

Number of Fisher Scoring iterations: 5

round(exp(coef(m_uil)), 3)

               (Intercept)              aantal_zorgen 
                     0.006                      2.815 
           eet_voor_slapen         sociale_isolatieja 
                     1.890                     12.324 
sociale_isolatietwijfelend 
                     2.368

round(exp(confint.default(m_uil)), 3)

                           2.5 % 97.5 %
(Intercept)                0.001  0.044
aantal_zorgen              1.865  4.248
eet_voor_slapen            1.294  2.762
sociale_isolatieja         3.295 46.086
sociale_isolatietwijfelend 0.780  7.188

car::Anova(m_uil, type = "II", test.statistic = "LR")

Analysis of Deviance Table (Type II tests)

Response: sliep_onrustig
                 LR Chisq Df      Pr(>Chisq)    
aantal_zorgen      40.875  1 0.0000000001623 ***
eet_voor_slapen    12.774  1       0.0003514 ***
sociale_isolatie   17.030  2       0.0002005 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

phat <- predict(m_uil, type = "response")
yhat <- ifelse(phat >= 0.5, 1, 0)
ctab <- table(observed = kerkuil_nacht$sliep_onrustig, predicted = yhat)
ctab

        predicted
observed  0  1
       0 55 12
       1 15 36

# Pseudo-R^2
1 - m_uil$deviance / m_uil$null.deviance

[1] 0.3276982

# Accuracy
sum(diag(ctab)) / sum(ctab)

[1] 0.7711864

a) Significant volgens Wald: aantal_zorgen ( $z = 4.93$ , p < .001) en eet_voor_slapen ( $z = 3.29$ , p = .001). Een $z$ -toets heeft geen vrijheidsgraden; wie liever de Wald- $\chi^2$ rapporteert, kwadrateert $z$ en schrijft voor aantal_zorgen $\chi^2(1) = 24.30$ . De factor sociale_isolatie heeft één significante dummy (ja vs nee: $z = 3.73$ , p < .001); twijfelend is niet-significant (p = .13).

b) $\mathit{OR}_{\text{zorgen}} = \exp(1.035) \approx 2.82$ met de hand; R rekent met de onafgeronde $b$ en geeft $2.81$ , en dat getal gaat in de APA-zin. 95%-CI $[1.87, 4.25]$ . Per extra zorg worden de odds van onrustig slapen bijna verdrievoudigd — meer zorgen, slechter slapen, weinig verrassend. Voor eet_voor_slapen $\mathit{OR} = 1.89$ per hapjes-eenheid.

c) Sociale_isolatie als geheel is significant: $\chi^2_{\text{LR}}(2) = 17.03$ , p < .001. Geïsoleerde dieren onderscheiden zich duidelijk van niet-geïsoleerde (positief), twijfelende dieren verschillen niet significant.

d) Pseudo- $R^2_L \approx .33$ ; accuracy $\approx .77$ .

Een logistische regressie van onrustig-slapen op aantal zorgen, eten voor het slapen en sociale isolatie verklaarde de uitkomst significant beter dan een null-model, $\chi^2_{\text{LR}}(4) = 52.89$ , p < .001, pseudo- $R^2_L = .33$ . Meer zorgen waren geassocieerd met onrustiger slaap ( $\mathit{OR} = 2.81$ , $95\%$ -CI $[1.87, 4.25]$ , p < .001); meer hapjes voor het slapen eveneens ( $\mathit{OR} = 1.89$ , $[1.29, 2.76]$ , p = .001). Sociale isolatie droeg significant bij ( $\chi^2_{\text{LR}}(2) = 17.03$ , p < .001): geïsoleerde dieren hadden hogere odds van onrustig slapen dan niet-geïsoleerde ( $\mathit{OR} = 12.32$ , $[3.30, 46.09]$ , p < .001). Het model classificeerde $77\%$ van de dieren correct.

u_h1 <- glm(sliep_onrustig ~ sociale_isolatie,
            data = kerkuil_nacht, family = binomial(link = "logit"))
u_h2 <- glm(sliep_onrustig ~ sociale_isolatie + aantal_zorgen,
            data = kerkuil_nacht, family = binomial(link = "logit"))
u_h3 <- glm(sliep_onrustig ~ sociale_isolatie + aantal_zorgen + eet_voor_slapen,
            data = kerkuil_nacht, family = binomial(link = "logit"))

# m0_uil is het null-model; per regel de toets van één blok.
anova(m0_uil, u_h1, u_h2, u_h3, test = "Chisq")

Analysis of Deviance Table

Model 1: sliep_onrustig ~ 1
Model 2: sliep_onrustig ~ sociale_isolatie
Model 3: sliep_onrustig ~ sociale_isolatie + aantal_zorgen
Model 4: sliep_onrustig ~ sociale_isolatie + aantal_zorgen + eet_voor_slapen
  Resid. Df Resid. Dev Df Deviance     Pr(>Chi)    
1       117     161.41                             
2       115     154.09  2    7.315    0.0257932 *  
3       114     121.29  1   32.803 0.0000000102 ***
4       113     108.51  1   12.774    0.0003514 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# R2_L per model, met r2_L() uit de hiërarchische LRA van de pad.
round(sapply(list(u_h1, u_h2, u_h3), r2_L), 3)

[1] 0.045 0.249 0.328

e) Ja, alle drie de blokken voegen significant iets toe. Sociale isolatie (blok 1, twee dummies): $\chi^2_{\text{LR}}(2) = 7.32$ , p = .026, $R^2_L = .045$ . Aantal zorgen erbij (blok 2): $\chi^2_{\text{LR}}(1) = 32.80$ , p < .001, $R^2_L = .249$ . Eten voor het slapen erbij (blok 3): $\chi^2_{\text{LR}}(1) = 12.77$ , p < .001, $R^2_L = .328$ . Het laatste model is hetzelfde als het volle model uit a: zijn $R^2_L$ van $.328$ is de $.33$ uit d, alleen op drie decimalen.

Wat blijft liggen

Onderwerpen die buiten dit thema blijven

LRA met logit-link en Hosmer-Lemeshow pseudo- $R^2$ is de hoofdingang van logistische analyse, maar er ligt veel achter. Wat dit oefenboek je gaf is genoeg voor de gangbare onderzoeks-rapportage; voor verdere studie:

Probit, multinomial en ordinal logistic regression — alternatieve link-functies (probit) of meerdere uitkomst-categorieën (multinomial, ordinal). Niet behandeld; dit oefenboek blijft bij binaire LRA met logit.
Penalized regression (ridge, lasso, elastic net) — voor predictieve modellen met veel voorspellers en/of weinig events. Lasso doet automatische voorspeller-selectie via $L_1$ -shrinkage.
ROC-curve en AUC — alternatieve evaluatie voorbij accuracy/PPV/NPV. Toont sensitivity tegen $1 - \text{specificity}$ over alle cutoff-keuzes; AUC vat het in één getal samen, ongevoelig voor prevalentie.
Mixed-effects logistische regressie — voor geclusterde binaire data (bv. herhaalde metingen per persoon, leerlingen genest in scholen). Pakket lme4::glmer().
Cox-Snell en Nagelkerke pseudo- $R^2$ — varianten van pseudo- $R^2$ die we hier niet behandelen. Cox-Snell heeft een maximumwaarde onder $1$ ; Nagelkerke is daar voor gecorrigeerd. In SPSS-uitvoer staan ze er wel, naast de $-2LL$ ; de $R^2_L$ in dit hoofdstuk is die van Hosmer en Lemeshow, en die reken je zelf uit die $-2LL$ ’s.
Bayesiaanse logistische regressie — alternatief inferentie-frame met posterior-distributies in plaats van $p$ -waardes; toegankelijke ingangen via brms of rstanarm.
Class-imbalance handling — bij sterk scheve uitkomstverdelingen (bv. $5\%$ events): oversampling, undersampling, threshold-tuning, of cost-sensitive learning.
Splines en GAM voor niet-lineaire voorspellers — mgcv::gam(family = binomial) als log-odds-lineariteit niet houdt.

Voor verdieping: Hosmer, Lemeshow & Sturdivant (2013, Applied Logistic Regression); Agresti (2018, An Introduction to Categorical Data Analysis); Long (1997, Regression Models for Categorical and Limited Dependent Variables).

Aan het eind van de dag

Toen het laat werd, klapte de pad zijn schrift dicht.

“Het zijn vooral de vossen,” zei hij. “Tien keer hogere odds dan de reeën. Dassen vrijwel nooit. Een tweedeling waar ik wat aan heb.”

“En de jongere dieren,” vervolgde hij. “Ieder jaar dat ze ouder worden, dalen de odds met zo’n 14 procent. Volwassenen hebben blijkbaar wel betere dingen te doen dan in het water staren.”

“En de afgewezen dieren?” vroeg de bever, die langs zwom om zijn aantekeningen te halen.

“Die ook,” zei de pad. “Per afwijzing een derde meer odds. Wie net is gepasseerd, wil even controleren of hij nog wel een gezicht heeft.”

“Driekwart van mijn schattingen klopt,” voegde hij toe. “Niet alles. Maar genoeg om te weten wie er morgen langs het water zal komen om zichzelf eens goed in de ogen te kijken.”

“Driekwart is meer dan ik vroeger had,” zei de bever.

Ze zaten bij het water. Iemand had ergens een tak laten vallen. Hij dreef voorbij. Een ree liep in de verte langs, schichtig, en bleef even staan bij de waterspiegel. De pad keek niet meer op.

Verantwoording

Dit oefenboek is geschreven voor studenten die logistische regressieanalyse leren via R. De didactische lijnen volgen de gangbare opbouw van Nederlandse universitaire MVDA-cursussen; alle voorbeelden, datasets, vragen en formuleringen in dit hoofdstuk zijn origineel. De gebruikte drempels en formules zijn standaard statistische conventies; verwijzingen naar Pedhazur (1997) en Hosmer-Lemeshow-Sturdivant (2013) volgen de gebruikelijke citaatpraktijk in dit veld.

Versie: October 2026 — CountCamp Lab

Terug naar boven

Al dan niet naar jezelf turen in de spiegel

4.0 Project- en datavoorbereiding

4.1 De pad telt — eerst met één voorspeller

4.1.a Null-model — hoe goed gokt de pad zonder verdere info?

4.1.b Een voorspeller toevoegen — aantal_afwijzingen

4.1.c Van bb naar kans — de sigmoid-formule

4.1.d Odds — een eerste laag ratio

4.1.e Odds ratio — een ratio van ratios

4.1.f Wald-toets — wijkt b1b_1 in de populatie af van nul?

4.1.g Likelihood-ratio-toets — doet het model als geheel iets?

4.2 Het volle model — meerdere voorspellers tegelijk

T1 — Techniek-keuze: LRA, MRA, ANOVA of ANCOVA?

T2 — Odds versus probability

4.2.a Verkenning — frequenties en wolken

4.2.b Het volle model fitten — alle voorspellers tegelijk

T3 — 𝑂𝑅\mathit{OR} uit bb via exp⁡()\exp()

T4 — 𝑂𝑅(k)=𝑂𝑅(1)k\mathit{OR}(k) = \mathit{OR}(1)^k: meerdere eenheden

4.2.c Odds ratios en hun 95%-CI — meerdere voorspellers naast elkaar

OR-tabel — alles op één plek

4.2.d Modelvergelijking — likelihood-ratio-toets bij meerdere voorspellers

T5 — Wald versus LR-toets: wanneer welke?

4.2.e Pseudo-R2R^2

Deviance-staafgrafiek — visualiseer RL2R^2_L

T6 — Pseudo-R2R^2: interpretatie en grenzen

4.2.f Hiërarchische LRA — het model in blokken opbouwen

4.2.g Classificatie — kansen, cutoff, accuracy

T7 — Classificatie-accuratesse + PPV/NPV

4.2.h APA-zin — alles in één

4.3 Aannamechecks

4.3.a Log-odds-lineariteit

4.3.b Multicollineariteit — vif()

4.3.c Sample-size — Pedhazur’s N/k≥30N/k \geq 30

T8 — Log-odds-lineariteit: wat is het, hoe checken?

T9 — N/k≥30N/k \geq 30: vergelijken met MRA

T10 — Niet-significante voorspeller weglaten?

R-spiekblad bij LRA

Pakketten en data laden

Verkenning — frequenties en boxplots

Model fitten — glm(..., family = binomial(link = "logit"))

Coëfficiënten en odds ratios

Modelvergelijking via likelihood ratio

Pseudo-R2R^2 van Hosmer en Lemeshow

Hiërarchische LRA — blok voor blok

Classificatietabel — predicted vs observed

Aannamechecks — vif en log-odds-lineariteit

Voorbeeld-tentamenvragen

Theorie en handreken

R-practical opdrachtje

Wat blijft liggen

Aan het eind van de dag

Verantwoording

4.1.b Een voorspeller toevoegen — `aantal_afwijzingen`

4.1.c Van $b$ naar kans — de sigmoid-formule

4.1.f Wald-toets — wijkt $b_1$ in de populatie af van nul?

T3 — $\mathit{OR}$ uit $b$ via $\exp()$

T4 — $\mathit{OR}(k) = \mathit{OR}(1)^k$ : meerdere eenheden

4.2.e Pseudo- $R^2$

Deviance-staafgrafiek — visualiseer $R^2_L$

T6 — Pseudo- $R^2$ : interpretatie en grenzen

4.3.b Multicollineariteit — `vif()`

4.3.c Sample-size — Pedhazur’s $N/k \geq 30$

T9 — $N/k \geq 30$ : vergelijken met MRA

Model fitten — `glm(..., family = binomial(link = "logit"))`

Pseudo- $R^2$ van Hosmer en Lemeshow