Wanneer het aankomt op het modelleren van tellingen (d.w.z. gehele getallen groter dan of gelijk aan 0), beginnen we vaak met Poisson regressie. Dit is een veralgemeend lineair model waarbij een respons verondersteld wordt een Poisson-verdeling te hebben, afhankelijk van een gewogen som van voorspellers. Wij zouden bijvoorbeeld het aantal gedocumenteerde hersenschuddingen bij NFL quarterbacks kunnen modelleren als een functie van de gespeelde snaps en het totale aantal jaren ervaring van zijn offensive line. Een potentieel nadeel van Poisson regressie is echter dat het de variabiliteit van de tellingen niet accuraat beschrijft.

Een Poisson-verdeling wordt geparametriseerd door het gemiddelde en de variantie. Hoewel handig om te onthouden, is het niet vaak realistisch. Een verdeling van tellingen zal meestal een variantie hebben die niet gelijk is aan het gemiddelde. Wanneer we dit zien gebeuren met gegevens waarvan we aannemen (of hopen) dat ze Poisson verdeeld zijn, zeggen we dat we onder- of overdispersie hebben, afhankelijk van of de variantie kleiner of groter is dan het gemiddelde. Het uitvoeren van Poisson regressie op telgegevens die dit gedrag vertonen resulteert in een model dat niet goed past.

Een benadering die dit probleem aanpakt is Negatieve Binomiale Regressie. De negatieve binomiale verdeling beschrijft, net als de Poisson-verdeling, de kans op het voorkomen van gehele getallen groter dan of gelijk aan 0. In tegenstelling tot de Poisson-verdeling zijn de variantie en het gemiddelde niet gelijkwaardig. Dit suggereert dat deze verdeling kan dienen als een nuttige benadering voor het modelleren van tellingen met een variabiliteit die afwijkt van het gemiddelde. De variantie van een negatieve binomiale verdeling is een functie van het gemiddelde en heeft een extra parameter, k, die de spreidingsparameter wordt genoemd. Stel dat onze telling een willekeurige variabele Y is uit een negatieve binomiale verdeling, dan is de variantie van Y

$$ var(Y) = \mu + \mu^{2}/k $$

Naarmate de spreidingsparameter groter wordt, convergeert de variantie naar dezelfde waarde als het gemiddelde, en verandert de negatieve binomiaal in een Poisson-verdeling.

Om de negatieve binomiale verdeling te illustreren, laten we werken met enkele gegevens uit het boek, Categorical Data Analysis, van Alan Agresti (2002). De gegevens staan in tabel 13.6 in paragraaf 13.4.3. De gegevens zijn afkomstig uit een enquête onder 1308 mensen waarin hun werd gevraagd hoeveel slachtoffers van moord zij kennen. De variabelen zijn resp, het aantal slachtoffers dat de respondent kent, en race, het ras van de respondent (zwart of blank). Kan ras helpen verklaren hoeveel slachtoffers van moord iemand kent? De gegevens moeten eerst in R worden ingevoerd:

> # Table 13.6> # Agresti, p. 561> black <- c(119,16,12,7,3,2,0)> white <- c(1070,60,14,4,0,0,1)> resp <- c(rep(0:6,times=black), rep(0:6,times=white))> race <- factor(c(rep("black", sum(black)), rep("white", sum(white))),+ levels = c("white","black"))> victim <- data.frame(resp, race)

Voordat we aan de modellering beginnen, laten we de gegevens eerst verkennen. Eerst zien we dat de meeste respondenten blank zijn:

> table(race)racewhite black 1149 159 

Zwarte mensen hebben een hoger gemiddeld aantal dan blanke mensen:

> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

Voor elk ras is de steekproefvariantie ruwweg het dubbele van het gemiddelde. Het lijkt erop dat er sprake is van overdispersie.

> with(victim, tapply(resp, race, var)) white black 0.1552448 1.1498288 

Ten slotte kijken we naar de verdeling van de tellingen naar ras.

> table(resp, race) raceresp white black 0 1070 119 1 60 16 2 14 12 3 4 7 4 0 3 5 0 2 6 1 0

Op naar de modelaanpassing. Eerst proberen we Poisson-regressie met behulp van de glm()-functie en tonen we een deel van de samenvattende uitvoer.

> # Poisson model> pGLM <- glm(resp ~ race, data=victim, family = poisson)> summary(pGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.38321 0.09713 -24.54 <2e-16 ***raceblack 1.73314 0.14657 11.82 <2e-16 ***

Race is zeer significant. Zwarten blijken veel vaker iemand te kennen die het slachtoffer van een moord is geweest. Maar wat betekent de coëfficiënt 1,73? In dit eenvoudige model met één dichotome voorspeller, is het het verschil in de logische verwachte tellingen. Als we de coëfficiënt exponenteren krijgen we een verhouding van steekproefgemiddelden:

> exp(coef(pGLM))raceblack 5.658419 > # same thing> mean(victim$resp)/mean(victim$resp) 5.658419

In feite, als we met dit model een voorspelling doen en de resultaten exponenteren, krijgen we de steekproefgemiddelden:

> exp(predict(pGLM, newdata = data.frame(race=c("white","black")))) 1 2 0.09225413 0.52201258 > # same thing> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

Dit betekent dat het aantal bekende slachtoffers voor blanken is verdeeld als een Poisson met gemiddelde en variantie gelijk aan 0,09, terwijl het aantal bekende slachtoffers voor blanken is verdeeld als een Poisson met gemiddelde en variantie gelijk aan 0,00.09, terwijl het aantal bekende slachtoffers voor zwarten is verdeeld als een Poisson met gemiddelde en variantie gelijk aan 0,52. Uit onze verkennende analyse weten we al dat de waargenomen varianties veel groter waren, dus moeten we niet al te blij zijn met de geschatte varianties van het model. Als we de gepaste tellingen bekijken, zien we nog meer bewijs voor het gebrek aan fit:

> # fitted counts for Poisson GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans 1 2 0.09225413 0.52201258 > fittedW <- dpois(0:6,lambda = fmeans) * sum(victim$race=="white") > fittedB <- dpois(0:6,lambda = fmeans) * sum(victim$race=="black") > data.frame(Response=0:6,BlackObs=black, BlackFit=round(fittedB,1), + WhiteObs=white, WhiteFit=round(fittedW,1)) Response BlackObs BlackFit WhiteObs WhiteFit1 0 119 94.3 1070 1047.72 1 16 49.2 60 96.73 2 12 12.9 14 4.54 3 7 2.2 4 0.15 4 3 0.3 0 0.06 5 2 0.0 0 0.07 6 0 0.0 1 0.0

Hierboven hebben we eerst de voorspelde gemiddelden opgeslagen in een object met de naam fmeans. Vervolgens hebben we passende tellingen gegenereerd door de dpois-functie te gebruiken samen met de geschatte gemiddelden om de kans op 0 tot 6 te voorspellen. Vervolgens vermenigvuldigen we die kansen met het aantal respondenten om aangepaste tellingen te verkrijgen. Tenslotte hebben wij alles gecombineerd in een dataframe om de waargenomen waarden gemakkelijk te kunnen vergelijken met de aangepaste waarden. Twee van de meer dramatische dingen om op te merken is dat we de 0-tellingen onderfitten en de 1-tellingen overfitten.

We kunnen een wortellogram gebruiken om de fit van een count regression model te visualiseren. De functie rootogram() in het countreg-pakket maakt dit gemakkelijk.

> countreg::rootogram(pGLM)

De rode gebogen lijn is de theoretische Poisson fit. “Aan elk punt van de rode lijn hangt een staaf waarvan de hoogte het verschil tussen de verwachte en de waargenomen telling weergeeft. Een staaf die onder 0 hangt, duidt op een underfitting. Een staaf boven 0 wijst op overfitting. De tellingen zijn getransformeerd met een vierkantsworteltransformatie om te voorkomen dat kleinere tellingen worden overschaduwd door grotere. We zien een grote mate van underfitting voor tellingen van 2 en hoger en een enorme overfitting voor de telling van 1.

Nu gaan we proberen een negatief binomiaal model te passen. We zagen dat de variabiliteit van de tellingen groter was voor beide rassen. Het lijkt erop dat de negatieve binomiale verdeling de verdeling van de tellingen beter benadert.

Om een negatief binomiaal model in R te passen, wenden we ons tot de glm.nb() functie in het MASS pakket (een pakket dat met R wordt meegeleverd). Ook hier tonen we slechts een deel van de samenvattende uitvoer:

> # Fit a negative binomial model> library(MASS)> nbGLM <- glm.nb(resp ~ race, data=victim)> summary(nbGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.3832 0.1172 -20.335 < 2e-16 ***raceblack 1.7331 0.2385 7.268 3.66e-13 *** Theta: 0.2023 

Merk eerst op dat de coëfficiënten dezelfde zijn als voorheen. Ook hier kunnen we de coëfficiënt exponentiëren om een verhouding van steekproefgemiddelden te krijgen en voorspellingen te doen om de oorspronkelijke steekproefgemiddelden te krijgen.

> # fitted counts for Negative Binomial GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans # same as pGLM 1 2 0.09225413 0.52201258 

Maar merk op dat de standaardfout voor de rascoëfficiënt groter is, wat wijst op meer onzekerheid in onze schatting (0,24 tegenover 0,15). Dit is logisch gezien de waargenomen variabiliteit in onze tellingen. Let ook op de schatting van Theta. Dat is onze spreidingsparameter. We kunnen hem als volgt rechtstreeks vanuit ons modelobject benaderen:

> nbGLM$theta 0.2023119

En we kunnen hem gebruiken om geschatte varianties voor de tellingen te krijgen:

> fmeans + fmeans^2 * (1/nbGLM$theta) 1 2 0.134322 1.868928 

Deze liggen veel dichter bij de waargenomen varianties dan die van het Poisson-model.

Wederom visualiseren we de fit met een rootogram:

> countreg::rootogram(nbGLM)

Dit ziet er veel beter uit dan het rootogram van het Poisson-model. Er is een lichte underfitting/overfitting voor de tellingen 1 tot en met 3, maar verder ziet het er vrij goed uit.

Om meer inzicht te krijgen in ons negatief binomiaal model, gebruiken we de parameters om gegevens te simuleren en vergelijken we de gesimuleerde gegevens met de waargenomen gegevens. Hieronder laden we het magrittr pakket voor toegang tot de %>% operator die ons in staat stelt om functies te “chainen”.

Eerst tabelleren we de tellingen en maken we een barplot voor respectievelijk de blanke en zwarte deelnemers. Vervolgens gebruiken we de modelparameters om gegevens van een negatieve binomiale verdeling te simuleren. De twee parameters zijn mu en size (d.w.z. dispersieparameter). Merk op dat we de coef() functie gebruiken om de juiste coëfficiënten voor elk ras te extraheren. Voor blank is dat alleen het intercept, voor zwart is dat het intercept en de helling (dus tellen we ze op). Dan exponentiëren we om van de logschaal naar de oorspronkelijke schaal te gaan. We simuleren hetzelfde aantal waarnemingen als we in onze oorspronkelijke gegevens hebben.

library(magrittr) # for %>% operatorop <- par(mfrow=c(1,2))set.seed(1)victim$resp %>% `)) %>% table() %>% barplot(main = "Simulated White")victim$resp %>% `[`(victim$race=="black") %>% table() %>% barplot(main = "Observed Black")rnbinom(n = 159, size = nbGLM$theta, mu = exp(sum(coef(nbGLM)))) %>% table() %>% barplot(main = "Simulated Black")par(op)

De gesimuleerde gegevens lijken erg op de waargenomen gegevens, wat ons opnieuw vertrouwen geeft in de keuze voor negatieve binomiale regressie om deze gegevens te modelleren. Deze plots tonen ook de voorwaardelijke aard van ons model aan. De negatieve binomiale verdeling van de tellingen hangt af van, of is geconditioneerd door, het ras. Elk ras heeft een ander gemiddelde maar een gemeenschappelijke spreidingsparameter.

• Agresti, Alan (2002), Categorical Data Analysis, Wiley.
• Kleiber, Christian & Zeileis, Achim (2016): Visualizing Count Data Regressions Using Rootograms, The American Statistician, DOI: 10.1080/00031305.2016.1173590

Voor vragen of verduidelijkingen over dit artikel kunt u contact opnemen met het UVA Library StatLab: [email protected]

Bekijk de gehele collectie UVA Library StatLab artikelen.

Clay Ford
Statistical Research Consultant
University of Virginia Library