Mikor a számlálások (azaz a 0-nál nagyobb vagy azzal egyenlő egész számok) modellezéséről van szó, gyakran Poisson-regresszióval kezdjük. Ez egy általánosított lineáris modell, ahol a válaszról feltételezzük, hogy Poisson-eloszlású, a prediktorok súlyozott összegének függvényében. Például modellezhetjük az NFL irányítóinak dokumentált agyrázkódásainak számát a játszott labdamenetek és a támadósorának teljes éves tapasztalata függvényében. A Poisson-regresszió egyik lehetséges hátránya azonban az, hogy nem feltétlenül írja le pontosan a számlálások változékonyságát.

A Poisson-eloszlást \(\lambda\) paraméterezi, ami történetesen az átlaga és a szórása is. Bár kényelmes megjegyezni, de gyakran nem reális. Egy számlálási eloszlásnak általában olyan varianciája van, amely nem egyenlő az átlagával. Amikor ezt látjuk olyan adatoknál, amelyekről feltételezzük (vagy reméljük), hogy Poisson-eloszlásúak, akkor alul- vagy felüldiszperzióról beszélünk, attól függően, hogy a variancia kisebb vagy nagyobb, mint az átlag. Az ilyen viselkedést mutató számlálási adatokon Poisson-regresszió végrehajtása olyan modellt eredményez, amely nem illeszkedik jól.

Az egyik megközelítés, amely ezt a problémát kezeli, a negatív binomiális regresszió. A negatív binomiális eloszlás a Poisson-eloszláshoz hasonlóan a 0-nál nagyobb vagy azzal egyenlő egész számok előfordulási valószínűségeit írja le. A Poisson-eloszlástól eltérően a variancia és az átlag nem egyenértékű. Ez arra utal, hogy hasznos közelítésként szolgálhat az átlagtól eltérő változékonyságú számlálások modellezésére. A negatív binomiális eloszlás varianciája az átlag függvénye, és rendelkezik egy további paraméterrel, k-val, amelyet szórási paraméternek nevezünk. Tegyük fel, hogy a számlálásunk egy negatív binomiális eloszlás Y véletlen változója, akkor Y varianciája

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

Amint a szórási paraméter egyre nagyobb lesz, a variancia az átlaggal azonos értékhez konvergál, és a negatív binomiális eloszlás Poisson-eloszlássá alakul.

A negatív binomiális eloszlás szemléltetésére dolgozzunk Alan Agresti (2002) Categorical Data Analysis című könyvéből származó néhány adattal. Az adatokat a 13.4.3. szakaszban található 13.6. táblázatban mutatjuk be. Az adatok egy 1308 ember körében végzett felmérésből származnak, amelyben megkérdezték, hogy hány emberölés áldozatát ismerik. A változók a következők: resp, a válaszadó által ismert áldozatok száma, és race, a válaszadó faji hovatartozása (fekete vagy fehér). Segít-e a faji hovatartozás megmagyarázni, hogy egy személy hány emberölés áldozatát ismeri? Az adatokat először be kell vinni az R-be:

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

Mielőtt rátérnénk a modellezésre, vizsgáljuk meg az adatokat. Először is észrevesszük, hogy a legtöbb válaszadó fehér:

> table(race)racewhite black 1149 159 

A feketék átlagos száma magasabb, mint a fehéreké:

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

Minden faj esetében a minta szórása nagyjából kétszerese az átlagnak. Úgy tűnik, hogy túldiszperzióval állunk szemben.

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

Végül nézzük meg a számok fajonkénti eloszlását.

> 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

Térjünk rá a modellillesztésre. Először a Poisson-regressziót próbáljuk ki a glm() függvény segítségével, és megmutatjuk az összefoglaló kimenet egy részét.

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

A faji hovatartozás nagyon jelentős. Úgy tűnik, hogy a feketék sokkal nagyobb valószínűséggel ismernek valakit, aki emberölés áldozata lett. De mit jelent az 1,73-as együttható? Ebben az egyszerű, egy dichotóm prediktort tartalmazó modellben ez a logaritmikus várható számok különbsége. Ha az együtthatót exponenciáljuk, akkor a mintaátlagok arányát kapjuk:

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

Sőt, ha előrejelzést készítünk ezzel a modellel, és az eredményeket exponenciáljuk, akkor a mintaátlagokat kapjuk:

> 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 

Ez azt mondja, hogy a fehérek ismert áldozatainak száma Poisson-eloszlású, átlaga és szórása 0-val egyenlő.09, míg a feketék esetében az ismert áldozatok száma Poisson-eloszlású, átlaga és szórása 0,52. A feltáró elemzésünkből már tudjuk, hogy a megfigyelt szórások sokkal nagyobbak voltak, így nem lehetünk túlságosan elégedettek a modell becsült szórásával. Ha megvizsgáljuk az illesztett számokat, még több bizonyítékot látunk az illeszkedés hiányára:

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

Fentebb először a megjósolt átlagokat mentettük el egy fmeans nevű objektumba. Ezután a dpois függvénnyel és a becsült átlagokkal együtt illesztett számokat generáltunk, hogy megjósoljuk a 0-tól 6-ig terjedő valószínűségeket. Ezután megszoroztuk ezeket a valószínűségeket a válaszadók számával, hogy megkapjuk az illesztett számokat. Végül mindent egy adatkeretben egyesítettünk, hogy könnyen összehasonlíthassuk a megfigyelt és az illesztett értékeket. A két drámaibb dolog, amit meg kell jegyeznünk, hogy alulillesztettük a 0-s számokat és túlillesztettük az 1-es számokat.

A számlálási regressziós modell illeszkedésének szemléltetésére használhatjuk a rootogramot. A countreg csomag rootogram() függvénye megkönnyíti ezt.

> countreg::rootogram(pGLM)

A piros görbe vonal az elméleti Poisson-illesztés. A piros vonal minden egyes pontjáról egy sáv “lóg” le, amelynek magassága a várható és a megfigyelt számok közötti különbséget jelöli. A 0 alá lógó sáv alulillesztést jelez. A 0 fölött lógó sáv túlilleszkedést jelez. A számlálásokat négyzetgyök-transzformációval alakítottuk át, hogy a kisebb számlálásokat ne fedjék el és ne nyomják el a nagyobb számlálások. A 2-es és magasabb számlálásoknál nagymértékű alulilleszkedést, az 1-es számlálásnál pedig masszív túlilleszkedést látunk.

Most próbáljunk meg egy negatív binomiális modellt illeszteni. Észrevettük, hogy a számlálások variabilitása mindkét faj esetében nagyobb. Úgy tűnik, hogy a negatív binomiális eloszlás jobban közelítené a számlálások eloszlását.

Negatív binomiális modell illesztéséhez R-ben a MASS csomag (az R-rel együtt telepített csomag) glm.nb() függvényéhez fordulunk. Ismét csak az összefoglaló kimenet egy részét mutatjuk meg:

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

Először is vegyük észre, hogy az együtthatók ugyanazok, mint korábban. Ismét exponenciálhatjuk a fajlagos együtthatót, hogy megkapjuk a mintaátlagok arányát, és előrejelzéseket végezhetünk, hogy megkapjuk az eredeti mintaátlagokat.

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

De vegyük észre, hogy a faji együttható standard hibája nagyobb, ami nagyobb bizonytalanságot jelez a becslésünkben (0,24 a 0,15-tel szemben). Ennek van értelme, tekintve a számlálásunkban megfigyelt változékonyságot. Figyeljük meg a Theta becslését is. Ez a szórási paraméterünk. Közvetlenül a modellobjektumunkból érhetjük el a következőképpen:

> nbGLM$theta 0.2023119

És használhatjuk arra, hogy becsült szórásokat kapjunk a számlálásokhoz:

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

Ezek sokkal közelebb állnak a megfigyelt szórásokhoz, mint a Poisson-modell által adottak.

Még egyszer vizualizáljuk az illeszkedést egy gyökogram segítségével:

> countreg::rootogram(nbGLM)

Ez sokkal jobban néz ki, mint a Poisson-modell gyökogramja. Az 1-3. számok esetében enyhe alul-/túlilleszkedés tapasztalható, de egyébként elég jól néz ki.

Hogy további betekintést nyerjünk a negatív binomiális modellünkbe, használjuk a paramétereit az adatok szimulálására, és hasonlítsuk össze a szimulált adatokat a megfigyelt adatokkal. Az alábbiakban betöltjük a magrittr csomagot, hogy hozzáférjünk a %>% operátorhoz, amely lehetővé teszi számunkra a függvények “láncolását”.

Először táblázatba foglaljuk a számlálásokat, és barplotot készítünk a fehér, illetve a fekete résztvevők számára. Ezután a modellparaméterek segítségével szimuláljuk a negatív binomiális eloszlásból származó adatokat. A két paraméter a mu és a méret (azaz a szórás paramétere). Vegyük észre, hogy a coef() függvényt használjuk az egyes fajokra vonatkozó megfelelő együtthatók kinyerésére. A fehérek esetében ez csak a metszéspont, a feketék esetében a metszéspont és a meredekség (tehát összegezzük őket). Ezután exponenciáljuk, hogy a log-skáláról az eredeti skálára konvertáljuk. Ugyanannyi megfigyelést szimulálunk, mint amennyi eredeti adatunk van.

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)

A szimulált adatok nagyon hasonlóak a megfigyelt adatokhoz, ami ismét bizalmat ad nekünk abban, hogy a negatív binomiális regressziót válasszuk ezen adatok modellezésére. Ezek az ábrák a modellünk feltételes jellegét is szemléltetik. A számlálások negatív binomiális eloszlása a faji hovatartozástól függ, vagy attól függ. Minden fajnak más az átlaga, de közös a szórásparamétere.

• 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

A cikkel kapcsolatos kérdések vagy pontosítások esetén forduljon az UVA Library StatLab-hoz: [email protected]

Nézze meg az UVA Library StatLab teljes cikkgyűjteményét.

Clay Ford
Statisztikai kutatási tanácsadó
University of Virginia Library