Jeśli chodzi o modelowanie liczebności (tzn. liczb całkowitych większych lub równych 0), często zaczynamy od regresji Poissona. Jest to uogólniony model liniowy, w którym zakłada się, że odpowiedź ma rozkład Poissona w zależności od ważonej sumy predyktorów. Na przykład, możemy modelować liczbę udokumentowanych wstrząsów mózgu u rozgrywających NFL jako funkcję rozegranych snapów i całkowitego wieloletniego doświadczenia linii ofensywnej. Jednak jedną z potencjalnych wad regresji Poissona jest to, że może ona niedokładnie opisywać zmienność zliczeń.

Rozkład Poissona jest parametryzowany przez \(\), która jest zarówno jego średnią, jak i wariancją. Choć jest to wygodne do zapamiętania, nie jest często realistyczne. Rozkład liczebności zwykle ma wariancję, która nie jest równa jego średniej. Kiedy widzimy, że dzieje się tak z danymi, które zakładamy (lub mamy nadzieję), że mają rozkład Poissona, mówimy, że mamy do czynienia z niedostatecznym lub nadmiernym rozproszeniem, w zależności od tego, czy wariancja jest mniejsza czy większa od średniej. Przeprowadzenie regresji Poissona na danych zliczeniowych, które wykazują takie zachowanie, skutkuje modelem, który nie pasuje dobrze.

Jednym z podejść, które rozwiązuje ten problem, jest ujemna regresja dwumianowa. Ujemny rozkład dwumianowy, podobnie jak rozkład Poissona, opisuje prawdopodobieństwo wystąpienia liczb całkowitych większych lub równych 0. W przeciwieństwie do rozkładu Poissona, wariancja i średnia nie są równoważne. To sugeruje, że może on służyć jako użyteczne przybliżenie do modelowania liczb o zmienności różnej od średniej. Wariancja rozkładu dwumianowego ujemnego jest funkcją jego średniej i ma dodatkowy parametr, k, zwany parametrem dyspersji. Powiedzmy, że naszym licznikiem jest zmienna losowa Y z rozkładu dwumianowego ujemnego, to wariancja Y wynosi

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

Gdy parametr dyspersji staje się coraz większy, wariancja zbiega do tej samej wartości co średnia, a dwumian ujemny zamienia się w rozkład Poissona.

Aby zilustrować ujemny rozkład dwumianowy, popracujmy z danymi z książki Categorical Data Analysis, autorstwa Alana Agresti (2002). Dane te są przedstawione w tabeli 13.6 w rozdziale 13.4.3. Dane te pochodzą z ankiety przeprowadzonej wśród 1308 osób, w której zapytano ich, ile ofiar zabójstw znają. Zmiennymi są resp, liczba ofiar, które zna respondent, oraz race, rasa respondenta (czarny lub biały). Czy rasa pomaga wyjaśnić, ile ofiar zabójstw zna dana osoba? Dane trzeba najpierw wprowadzić do programu 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)

Zanim przejdziemy do modelowania, zbadajmy dane. Najpierw zauważamy, że większość respondentów to biali:

> table(race)racewhite black 1149 159 

Czarni mają wyższą średnią liczebność niż biali:

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

Dla każdej rasy wariancja próby jest mniej więcej dwukrotnie większa od średniej. Wygląda na to, że mamy nadmiar dyspersji.

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

W końcu przyjrzyjmy się rozkładowi liczebności według 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

Przejdziemy do dopasowania modelu. Najpierw wypróbowujemy regresję Poissona przy użyciu funkcji glm() i pokazujemy część podsumowania danych wyjściowych.

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

Rasa jest bardzo znacząca. Wygląda na to, że czarni są znacznie bardziej skłonni znać kogoś, kto był ofiarą zabójstwa. Ale co oznacza współczynnik 1,73? W tym prostym modelu z jednym dychotomicznym predyktorem, jest to różnica w log spodziewanych liczbach. Jeśli wykładamy współczynnik, otrzymujemy stosunek średnich próbkowych:

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

W rzeczywistości, jeśli dokonamy predykcji za pomocą tego modelu i wykładamy wyniki, otrzymujemy średnie próbkowe:

> 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 

To mówi, że liczba znanych ofiar dla białych jest rozłożona jako rozkład Poissona ze średnią i wariancją równą 0.09, podczas gdy liczba znanych ofiar dla czarnych jest rozłożona jako rozkład Poissona ze średnią i wariancją równą 0.52. Wiemy już z naszej analizy eksploracyjnej, że obserwowane wariancje były znacznie większe, więc nie powinniśmy być zbyt zadowoleni z oszacowanych wariancji modelu. Jeśli zbadamy dopasowane liczebności, zobaczymy jeszcze więcej dowodów na brak dopasowania:

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

Powyżej najpierw zapisaliśmy przewidywane średnie do obiektu o nazwie fmeans. Następnie wygenerowaliśmy dopasowane zliczenia, używając funkcji dpois wraz z oszacowanymi środkami, aby przewidzieć prawdopodobieństwo uzyskania wartości od 0 do 6. Następnie pomnożyliśmy te prawdopodobieństwa przez liczbę respondentów, aby uzyskać dopasowane liczebności. Na koniec połączyliśmy wszystko w ramkę danych, aby łatwo porównać wartości obserwowane i dopasowane. Dwie z bardziej dramatycznych rzeczy, które należy zauważyć, to fakt, że jesteśmy niedopasowani do liczebności 0 i przepasowani do liczebności 1.

Możemy użyć rootogramu do wizualizacji dopasowania modelu regresji zliczeń. Funkcja rootogram() w pakiecie countreg ułatwia to zadanie.

> countreg::rootogram(pGLM)

Czerwona zakrzywiona linia to teoretyczne dopasowanie Poissona. Z każdego punktu na czerwonej linii “zwisa” słupek, którego wysokość reprezentuje różnicę między liczbą oczekiwaną a obserwowaną. Słupek wiszący poniżej 0 wskazuje na niedopasowanie. Słupek zwisający powyżej 0 wskazuje na przepasowanie. Liczby zostały przekształcone za pomocą pierwiastka kwadratowego, aby zapobiec przesłonięciu i przytłoczeniu mniejszych liczb przez większe. Widzimy duże niedopasowanie dla liczebników 2 i wyższych oraz masywne przepasowanie dla liczebnika 1.

Spróbujmy teraz dopasować model dwumianowy ujemny. Zauważyliśmy, że zmienność zliczeń była większa dla obu ras. Wydaje się, że ujemny rozkład dwumianowy lepiej przybliżałby rozkład zliczeń.

Aby dopasować ujemny model dwumianowy w R zwróćmy się do funkcji glm.nb() w pakiecie MASS (pakiet, który jest zainstalowany z R). Ponownie pokażemy tylko część podsumowania:

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

Po pierwsze zauważ, że współczynniki są takie same jak poprzednio. Po raz kolejny możemy wykładać współczynnik wyścigu, aby uzyskać stosunek środków próbki i dokonać predykcji, aby uzyskać oryginalne środki próbki.

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

Zauważ jednak, że błąd standardowy dla współczynnika rasy jest większy, wskazując na większą niepewność w naszym oszacowaniu (0,24 w porównaniu z 0,15). Ma to sens, biorąc pod uwagę obserwowaną zmienność w naszych liczeniach. Zauważmy również oszacowanie współczynnika Theta. Jest to nasz parametr dyspersji. Możemy uzyskać do niego dostęp bezpośrednio z naszego obiektu modelu w następujący sposób:

> nbGLM$theta 0.2023119

I możemy go użyć do uzyskania szacunkowych wariancji dla zliczeń:

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

Są one znacznie bliższe obserwowanym wariancjom niż te, które daje model Poissona.

Po raz kolejny wizualizujemy dopasowanie za pomocą rootogramu:

> countreg::rootogram(nbGLM)

Wygląda to znacznie lepiej niż rootogram modelu Poissona. Istnieje niewielkie niedopasowanie / przepasowanie dla zliczeń od 1 do 3, ale poza tym wygląda całkiem dobrze.

Aby uzyskać dalszy wgląd w nasz ujemny model dwumianowy, użyjmy jego parametrów do symulacji danych i porównajmy symulowane dane z danymi obserwowanymi. Poniżej ładujemy pakiet magrittr, aby uzyskać dostęp do operatora %>%, który pozwala nam na “łańcuchowanie” funkcji.

Najpierw tabulujemy zliczenia i tworzymy wykres słupkowy odpowiednio dla białych i czarnych uczestników. Następnie używamy parametrów modelu do symulowania danych z rozkładu dwumianowego ujemnego. Te dwa parametry to mu i wielkość (tzn. parametr dyspersji). Zauważ, że używamy funkcji coef() do wyodrębnienia odpowiednich współczynników dla każdej rasy. Dla białych jest to tylko punkt przecięcia, dla czarnych jest to punkt przecięcia i nachylenie (więc je sumujemy). Następnie wykładamy, aby przekonwertować ze skali logarytmicznej do skali oryginalnej. Symulujemy taką samą liczbę obserwacji, jaką mamy w naszych oryginalnych danych.

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)

Symulowane dane są bardzo podobne do danych obserwowanych, co ponownie daje nam pewność co do wyboru ujemnej regresji dwumianowej do modelowania tych danych. Te wykresy pokazują również warunkową naturę naszego modelu. Ujemny rozkład dwumianowy zliczeń zależy od rasy lub jest od niej uwarunkowany. Każda rasa ma inną średnią, ale wspólny parametr dyspersji.

• 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

W przypadku pytań lub wyjaśnień dotyczących tego artykułu, skontaktuj się z UVA Library StatLab: [email protected]

Zobacz całą kolekcję artykułów UVA Library StatLab.

Clay Ford
Statistical Research Consultant
University of Virginia Library
.