När det gäller modellering av antal (dvs. hela tal som är större än eller lika med 0) börjar vi ofta med Poisson-regression. Detta är en generaliserad linjär modell där ett svar antas ha en Poissonfördelning betingad av en viktad summa av prediktorer. Vi kan till exempel modellera antalet dokumenterade hjärnskakningar hos NFL-quarterbacks som en funktion av antalet spelade ögonblick och det totala antalet år av erfarenhet hos hans offensiva linje. En potentiell nackdel med Poisson-regressionen är dock att den kanske inte beskriver variabiliteten i räkningarna på ett korrekt sätt.

En Poissonfördelning parametreras av \(\lambda\), som råkar vara både dess medelvärde och varians. Även om det är bekvämt att komma ihåg är det inte ofta realistiskt. En fördelning av antalstillfällen kommer vanligtvis att ha en varians som inte är lika stor som dess medelvärde. När vi ser detta hända med data som vi antar (eller hoppas) är Poissonfördelad, säger vi att vi har under- eller överspridning, beroende på om variansen är mindre eller större än medelvärdet. Att utföra Poisson-regression på räkneuppgifter som uppvisar detta beteende resulterar i en modell som inte passar bra.

En metod som behandlar detta problem är negativ binomialregression. Den negativa binomialfördelningen, liksom Poissonfördelningen, beskriver sannolikheten för förekomsten av hela tal som är större än eller lika med 0. Till skillnad från Poissonfördelningen är variansen och medelvärdet inte likvärdiga. Detta tyder på att den kan fungera som en användbar approximation för att modellera tal med en variabilitet som skiljer sig från dess medelvärde. Variansen för en negativ binomialfördelning är en funktion av dess medelvärde och har ytterligare en parameter, k, som kallas spridningsparameter. Säg att vår räkning är en slumpmässig variabel Y från en negativ binomialfördelning, så är variansen för Y

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

När spridningsparametern blir större och större konvergerar variansen mot samma värde som medelvärdet, och den negativa binomialfördelningen förvandlas till en Poissonfördelning.

För att illustrera den negativa binomialfördelningen kan vi arbeta med några data från boken Categorical Data Analysis av Alan Agresti (2002). Uppgifterna presenteras i tabell 13.6 i avsnitt 13.4.3. Uppgifterna kommer från en undersökning med 1308 personer där de tillfrågades om hur många mordoffer de känner till. Variablerna är resp, antalet offer som respondenten känner till, och race, respondentens ras (svart eller vit). Hjälper ras till att förklara hur många mordoffer en person känner till? Uppgifterna måste först föras in 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)

För att komma till modellering ska vi utforska uppgifterna. Först märker vi att de flesta respondenterna är vita:

> table(race)racewhite black 1149 159 

Svarta har ett högre medelantal än vita:

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

För varje ras är urvalets varians ungefär dubbelt så stor som medelvärdet. Det verkar som om vi har överspridning.

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

Slutligt tittar vi på fördelningen av antalen per 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

Om till modellanpassning. Först prövar vi Poisson-regression med hjälp av funktionen glm() och visar en del av sammanfattningen.

> # 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 är mycket signifikant. Det verkar som om svarta har mycket större sannolikhet att känna någon som varit offer för ett mord. Men vad betyder koefficienten 1,73? I denna enkla modell med en dikotom prediktor är det skillnaden i log förväntat antal. Om vi exponerar koefficienten får vi ett förhållande mellan provmedlen:

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

Om vi gör en förutsägelse med den här modellen och exponerar resultaten får vi provmedlen:

> 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 

Detta säger att antalet kända offer för vita är fördelat som en Poisson med medelvärde och varians lika med 0.09, medan antalet kända offer för svarta är fördelat som en Poisson med medelvärde och varians lika med 0,52. Vi vet redan från vår explorativa analys att de observerade varianserna var mycket större, så vi bör inte vara alltför nöjda med modellens uppskattade varianser. Om vi undersöker de inpassade räknevärdena ser vi ännu fler bevis för bristen på anpassning:

> # 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

Ovan sparade vi först de förutspådda medelvärdena i ett objekt som heter fmeans. Vi genererade sedan anpassade räkningar genom att använda funktionen dpois tillsammans med de uppskattade medelvärdena för att förutsäga sannolikheten för att få 0 till 6. Vi multiplicerade sedan dessa sannolikheter med antalet respondenter för att få fram anpassade antal. Slutligen kombinerade vi allt i en dataram för att enkelt kunna jämföra observerade och anpassade värden. Två av de mer dramatiska sakerna att notera är att vi underanpassar 0-räkningarna och överanpassar 1-räkningarna.

Vi kan använda ett rotogram för att visualisera anpassningen av en regressionsmodell för antal. Funktionen rootogram() i paketet countreg gör detta enkelt.

> countreg::rootogram(pGLM)

Den röda krökta linjen är den teoretiska Poissonanpassningen. Från varje punkt på den röda linjen “hänger” en stapel vars höjd representerar skillnaden mellan förväntat och observerat antal. En stapel som hänger under 0 tyder på underanpassning. En stapel över 0 indikerar en överanpassning. Räknarna har omvandlats med en kvadratrotsomvandling för att förhindra att mindre räkningar skyms och överskuggas av större räkningar. Vi ser en hel del underanpassning för antal 2 och högre och massiv överanpassning för antal 1.

Nu ska vi försöka anpassa en negativ binomialmodell. Vi märkte att variabiliteten i räkningarna var större för båda raserna. Det verkar som om den negativa binomialfördelningen bättre skulle approximera fördelningen av räkningarna.

För att anpassa en negativ binomialmodell i R vänder vi oss till glm.nb()-funktionen i MASS-paketet (ett paket som installeras med R). Återigen visar vi bara en del av den sammanfattande utgången:

> # 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 noterar vi att koefficienterna är desamma som tidigare. Återigen kan vi exponera kappkoefficienten för att få ett förhållande mellan provmedlen och göra förutsägelser för att få fram de ursprungliga provmedlen.

> # 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 lägg märke till att standardfelet för raskoefficienten är större, vilket tyder på större osäkerhet i vår uppskattning (0,24 jämfört med 0,15). Detta är logiskt med tanke på den observerade variabiliteten i våra räkningar. Lägg också märke till uppskattningen av Theta. Det är vår spridningsparameter. Vi kan få tillgång till den direkt från vårt modellobjekt på följande sätt:

> nbGLM$theta 0.2023119

Och vi kan använda den för att få skattade varianser för räkningarna:

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

Dessa ligger mycket närmare de observerade varianserna än de som ges av Poissonmodellen.

En gång till visualiserar vi anpassningen med hjälp av ett rotogram:

> countreg::rootogram(nbGLM)

Detta ser mycket bättre ut än Poissonmodellens rotogram. Det finns en liten underanpassning/överanpassning för antal 1 till 3, men i övrigt ser det ganska bra ut.

För att få ytterligare insikter om vår negativa binomialmodell kan vi använda dess parametrar för att simulera data och jämföra de simulerade data med de observerade data. Nedan laddar vi paketet magrittr för att få tillgång till operatören %>% som gör det möjligt för oss att “kedja” funktioner.

Först tabellerar vi antalen och skapar ett stapeldiagram för de vita respektive svarta deltagarna. Därefter använder vi modellparametrarna för att simulera data från en negativ binomialfördelning. De två parametrarna är mu och storlek (dvs. spridningsparameter). Lägg märke till att vi använder funktionen coef() för att extrahera lämpliga koefficienter för varje ras. För vita är det bara interceptet, för svarta är det interceptet och lutningen (därför summerar vi dem). Sedan exponerar vi för att konvertera från logaritmisk skala till den ursprungliga skalan. Vi simulerar samma antal observationer som vi har i våra originaldata.

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)

Den simulerade datan är mycket lik den observerade datan, vilket återigen ger oss förtroende för att välja negativ binomialregression för att modellera denna data. Dessa diagram visar också att vår modell är villkorlig. Den negativa binomialfördelningen av antalen beror på, eller är betingad av, ras. Varje ras har ett annat medelvärde men en gemensam spridningsparameter.

• 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

För frågor eller förtydliganden angående denna artikel, kontakta UVA Library StatLab: [email protected]

Se hela samlingen av artiklar från UVA Library StatLab.

Clay Ford
Statistisk forskningskonsult
University of Virginia Library