Når det kommer til modellering af tællinger (dvs. hele tal større end eller lig med 0), starter vi ofte med Poisson-regression. Dette er en generaliseret lineær model, hvor et svar antages at have en Poisson-fordeling betinget af en vægtet sum af prædiktorer. Vi kan f.eks. modellere antallet af dokumenterede hjernerystelser hos NFL-quarterbacks som en funktion af antallet af spillede snaps og det samlede antal års erfaring hos hans offensive linje. En potentiel ulempe ved Poisson-regression er imidlertid, at den måske ikke nøjagtigt beskriver variabiliteten i tællingerne.
En Poisson-fordeling er parameteriseret ved \(\lambda\), som tilfældigvis er både dens middelværdi og varians. Selv om det er praktisk at huske, er det ikke ofte realistisk. En fordeling af tællinger vil normalt have en varians, der ikke er lig med dens middelværdi. Når vi ser dette ske med data, som vi antager (eller håber) er Poisson-fordelt, siger vi, at vi har under- eller overspredning, afhængigt af om variansen er mindre eller større end middelværdien. Hvis man udfører Poisson-regression på tælledata, der udviser denne adfærd, resulterer det i en model, der ikke passer godt.
En tilgang, der behandler dette problem, er Negative Binomialregression. Den negative binomialfordeling beskriver ligesom Poisson-fordelingen sandsynlighederne for forekomsten af hele tal, der er større end eller lig med 0. I modsætning til Poisson-fordelingen er variansen og middelværdien ikke ækvivalente. Dette tyder på, at den kan tjene som en nyttig tilnærmelse til modellering af tællinger med en variabilitet, der er forskellig fra dens gennemsnit. Variansen af en negativ binomialfordeling er en funktion af dens middelværdi og har en ekstra parameter, k, som kaldes spredningsparameteren. Hvis vores tælling er en tilfældig variabel Y fra en negativ binomialfordeling, er variansen af Y
$$$ var(Y) = \mu + \mu^{2}/k $$
Da sprednings parameteren bliver større og større, konvergerer variansen mod den samme værdi som middelværdien, og den negative binomialfordeling bliver til en Poisson-fordeling.
For at illustrere den negative binomialfordeling, lad os arbejde med nogle data fra bogen Categorical Data Analysis af Alan Agresti (2002). Dataene er præsenteret i tabel 13.6 i afsnit 13.4.3. Dataene stammer fra en undersøgelse af 1308 personer, hvor de blev spurgt om, hvor mange mordofre de kender. Variablerne er resp, antallet af ofre, som respondenten kender, og race, respondentens race (sort eller hvid). Er race med til at forklare, hvor mange mordofre en person kender? Dataene skal først indtastes i R:
> # 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)
Hvor vi går i gang med modelleringen, skal vi undersøge dataene. Først bemærker vi, at de fleste respondenter er hvide:
> table(race)racewhite black 1149 159
Sorte har et højere gennemsnitligt antal end hvide:
> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258
For hver race er stikprøvevariansen ca. dobbelt så stor som gennemsnittet. Det ser ud til, at vi har overspredning.
> with(victim, tapply(resp, race, var)) white black 0.1552448 1.1498288
Sidst ser vi på fordelingen af tællingerne efter race.
> 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 til modeltilpasning. Først prøver vi Poisson-regression ved hjælp af glm()-funktionen og viser en del af det sammenfattende output.
> # 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 er meget signifikant. Det ser ud til, at sorte er meget mere tilbøjelige til at kende nogen, der har været offer for et drab. Men hvad betyder koefficienten 1,73? I denne enkle model med én dikotom prædiktor er det forskellen i logaritmisk forventede antal. Hvis vi eksponerer koefficienten, får vi et forhold mellem stikprøvens gennemsnit:
> exp(coef(pGLM))raceblack 5.658419 > # same thing> mean(victim$resp)/mean(victim$resp) 5.658419
Hvis vi laver en forudsigelse med denne model og eksponerer resultaterne, får vi faktisk stikprøvens gennemsnit:
> 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
Dette siger, at antallet af kendte ofre for hvide er fordelt som en Poisson med middelværdi og varians lig med 0.09, mens antallet af kendte ofre for sorte er fordelt som en Poisson med middelværdi og varians lig med 0,52. Vi ved allerede fra vores eksplorative analyse, at de observerede varianser var meget større, så vi bør ikke være alt for tilfredse med modellens estimerede varianser. Hvis vi undersøger de tilpassede tællinger, vil vi se endnu flere beviser for den manglende tilpasning:
> # 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
Ovenfor gemte vi først de forudsagte middelværdier i et objekt kaldet fmeans. Vi genererede derefter tilpassede tællinger ved at bruge dpois-funktionen sammen med de estimerede middelværdier til at forudsige sandsynligheden for at få 0 til 6. Vi multiplicerede derefter disse sandsynligheder med antallet af respondenter for at opnå tilpassede tællinger. Til sidst kombinerede vi det hele i en dataramme for nemt at kunne sammenligne observerede og tilpassede værdier. To af de mere dramatiske ting, der er værd at bemærke, er, at vi underjusterer 0-tallene og overjusterer 1-tallene.
Vi kan bruge et rodogram til at visualisere tilpasningen af en regressionsmodel for antalsregression. Funktionen rootogram() i pakken countreg gør dette let.
> countreg::rootogram(pGLM)
Den røde buede linje er den teoretiske Poisson-tilpasning. “Hængende” fra hvert punkt på den røde linje er en søjle, hvis højde repræsenterer forskellen mellem det forventede og det observerede antal. En søjle, der hænger under 0, indikerer en undertilpasning. En søjle, der hænger over 0, angiver en overtilpasning. Tællingerne er blevet transformeret med en kvadratrodstransformation for at forhindre, at mindre tællinger bliver overskygget og overvældet af større tællinger. Vi ser en stor undertilpasning for tællinger 2 og højere og massiv overtilpasning for tælling 1.
Lad os nu prøve at tilpasse en negativ binomialmodel. Vi bemærkede, at variabiliteten i tællingerne var større for begge racer. Det ser ud til, at den negative binomialfordeling bedre ville tilnærme sig fordelingen af tællingerne.
For at tilpasse en negativ binomialmodel i R vender vi os til glm.nb()-funktionen i MASS-pakken (en pakke, der er installeret sammen med R). Igen viser vi kun en del af det sammenfattende output:
> # 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
Først bemærker vi, at koefficienterne er de samme som før. Igen kan vi eksponere løbskoefficienten for at få et forhold mellem stikprøvens middelværdier og foretage forudsigelser for at få de oprindelige stikprøvens middelværdier.
> # 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
Men bemærk, at standardfejlen for racekoefficienten er større, hvilket indikerer større usikkerhed i vores skøn (0,24 mod 0,15). Dette giver mening i betragtning af den observerede variabilitet i vores optællinger. Bemærk også estimatet for Theta. Det er vores spredningsparameter. Vi kan få adgang til den direkte fra vores modelobjekt på følgende måde:
> nbGLM$theta 0.2023119
Og vi kan bruge den til at få estimerede varianser for tællingerne:
> fmeans + fmeans^2 * (1/nbGLM$theta) 1 2 0.134322 1.868928
Disse er meget tættere på de observerede varianser end dem, der er givet af Poisson-modellen.
En gang til visualiserer vi tilpasningen ved hjælp af et rodogram:
> countreg::rootogram(nbGLM)
Dette ser meget bedre ud end Poisson-modellens rodogram. Der er en lille undertilpasning/overtilpasning for tællinger 1 til 3, men ellers ser det ret godt ud.
For at få yderligere indsigt i vores negative binomialmodel, skal vi bruge dens parametre til at simulere data og sammenligne de simulerede data med de observerede data. Nedenfor indlæser vi pakken magrittr for at få adgang til operatoren %>%, som giver os mulighed for at “kæde” funktioner.
Først tabulerer vi tællingerne og opretter et søjlediagram for henholdsvis de hvide og sorte deltagere. Derefter bruger vi modelparametrene til at simulere data fra en negativ binomialfordeling. De to parametre er mu og størrelse (dvs. dispersionsparameter). Bemærk, at vi bruger funktionen coef() til at udtrække de relevante koefficienter for hver race. For hvid er det kun interceptet, for sort er det interceptet og hældningen (derfor summerer vi dem). Derefter eksponerer vi for at konvertere fra log-skalaen til den oprindelige skala. Vi simulerer det samme antal observationer, som vi har i vores oprindelige data.
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 simulerede data er meget lig de observerede data, hvilket igen giver os tillid til at vælge negativ binomialregression til at modellere disse data. Disse plotter viser også den betingede karakter af vores model. Den negative binomialfordeling af tællingerne afhænger af, eller er betinget af, race. Hver race har en forskellig middelværdi, men en fælles spredningsparameter.
- 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
For spørgsmål eller afklaringer vedrørende denne artikel, kontakt UVA Library StatLab: [email protected]
Se hele samlingen af UVA Library StatLab artikler.
Clay Ford
Statistical Research Consultant
University of Virginia Library