Wenn es um die Modellierung von Zählungen (d. h. ganze Zahlen größer oder gleich 0) geht, beginnen wir oft mit der Poisson-Regression. Dabei handelt es sich um ein verallgemeinertes lineares Modell, bei dem davon ausgegangen wird, dass eine Antwort eine Poisson-Verteilung in Abhängigkeit von einer gewichteten Summe von Prädiktoren aufweist. So könnte man beispielsweise die Anzahl der dokumentierten Gehirnerschütterungen bei NFL-Quarterbacks als Funktion der gespielten Snaps und der Gesamtzahl der Jahre an Erfahrung der Offensive Line modellieren. Ein potenzieller Nachteil der Poisson-Regression ist jedoch, dass sie die Variabilität der Zählungen möglicherweise nicht genau beschreibt.
Eine Poisson-Verteilung wird durch \(\lambda\) parametrisiert, das zufällig sowohl ihr Mittelwert als auch ihre Varianz ist. Das ist zwar bequem zu merken, aber oft nicht realistisch. Eine Verteilung von Zählungen hat in der Regel eine Varianz, die nicht gleich ihrem Mittelwert ist. Wenn dies bei Daten der Fall ist, von denen wir annehmen (oder hoffen), dass sie Poisson-verteilt sind, sprechen wir von einer Unter- oder Überdispersion, je nachdem, ob die Varianz kleiner oder größer als der Mittelwert ist. Die Durchführung einer Poisson-Regression auf Zähldaten, die dieses Verhalten aufweisen, führt zu einem Modell, das nicht gut passt.
Ein Ansatz, der dieses Problem angeht, ist die Negative Binomialregression. Die negative Binomialverteilung beschreibt wie die Poisson-Verteilung die Wahrscheinlichkeiten für das Auftreten ganzer Zahlen größer oder gleich 0. Im Gegensatz zur Poisson-Verteilung sind Varianz und Mittelwert nicht äquivalent. Dies deutet darauf hin, dass sie als nützlicher Näherungswert für die Modellierung von Zählungen mit einer vom Mittelwert abweichenden Variabilität dienen kann. Die Varianz einer negativen Binomialverteilung ist eine Funktion ihres Mittelwerts und hat einen zusätzlichen Parameter, k, der als Dispersionsparameter bezeichnet wird. Angenommen, unsere Zählung ist eine Zufallsvariable Y aus einer negativen Binomialverteilung, dann ist die Varianz von Y
$$ var(Y) = \mu + \mu^{2}/k $$
Wenn der Dispersionsparameter immer größer wird, konvergiert die Varianz zum gleichen Wert wie der Mittelwert, und die negative Binomialverteilung wird zu einer Poissonverteilung.
Um die negative Binomialverteilung zu veranschaulichen, arbeiten wir mit einigen Daten aus dem Buch Categorical Data Analysis von Alan Agresti (2002). Die Daten sind in Tabelle 13.6 in Abschnitt 13.4.3 dargestellt. Die Daten stammen aus einer Umfrage unter 1308 Personen, in der sie gefragt wurden, wie viele Mordopfer sie kennen. Die Variablen sind resp, die Anzahl der Opfer, die der Befragte kennt, und race, die Rasse des Befragten (schwarz oder weiß). Lässt sich anhand der Rasse erklären, wie viele Mordopfer eine Person kennt? Die Daten müssen zunächst in R eingegeben werden:
> # 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)
Bevor wir mit der Modellierung beginnen, wollen wir die Daten untersuchen. Zunächst stellen wir fest, dass die meisten Befragten weiß sind:
> table(race)racewhite black 1149 159
Schwarze haben eine höhere durchschnittliche Anzahl als Weiße:
> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258
Für jede Rasse ist die Stichprobenvarianz ungefähr doppelt so hoch wie der Mittelwert. Es scheint, dass wir eine Überdispersion haben.
> with(victim, tapply(resp, race, var)) white black 0.1552448 1.1498288
Schließlich betrachten wir die Verteilung der Zählungen nach Rasse.
> 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
Weiter zur Modellanpassung. Zunächst versuchen wir die Poisson-Regression mit der Funktion glm() und zeigen einen Teil der zusammenfassenden Ausgabe.
> # 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 ***
Die Rasse ist sehr signifikant. Es scheint, dass Schwarze viel wahrscheinlicher jemanden kennen, der Opfer eines Tötungsdelikts war. Aber was bedeutet der Koeffizient 1,73? In diesem einfachen Modell mit einem dichotomen Prädiktor ist er die Differenz der logarithmischen erwarteten Zählungen. Wenn wir den Koeffizienten potenzieren, erhalten wir ein Verhältnis der Stichprobenmittelwerte:
> exp(coef(pGLM))raceblack 5.658419 > # same thing> mean(victim$resp)/mean(victim$resp) 5.658419
Wenn wir mit diesem Modell eine Vorhersage machen und die Ergebnisse potenzieren, erhalten wir die Stichprobenmittelwerte:
> 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
Dies besagt, dass die Anzahl der bekannten Opfer für Weiße als Poisson mit Mittelwert und Varianz gleich 0 verteilt ist.09, während die Anzahl der bekannten Opfer für Schwarze als Poisson mit Mittelwert und Varianz gleich 0,52 verteilt ist. Wir wissen bereits aus unserer explorativen Analyse, dass die beobachteten Varianzen viel größer waren, so dass wir mit den geschätzten Varianzen des Modells nicht allzu zufrieden sein sollten. Wenn wir die angepassten Zählungen untersuchen, sehen wir noch mehr Beweise für die mangelnde Anpassung:
> # 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
Oben haben wir zuerst die vorhergesagten Mittelwerte in einem Objekt namens fmeans gespeichert. Dann haben wir angepasste Zählungen erzeugt, indem wir die Funktion dpois zusammen mit den geschätzten Mittelwerten verwendet haben, um die Wahrscheinlichkeit des Erreichens von 0 bis 6 vorherzusagen. Anschließend haben wir diese Wahrscheinlichkeiten mit der Anzahl der Befragten multipliziert, um die angepasste Anzahl zu erhalten. Schließlich haben wir alles in einem Datenrahmen zusammengefasst, um die beobachteten und angepassten Werte leicht vergleichen zu können. Besonders auffällig ist, dass die 0-Werte unterangepasst und die 1-Werte überangepasst sind.
Wir können ein Rootogramm verwenden, um die Anpassung eines Zählungsregressionsmodells zu visualisieren. Die Funktion rootogram() im countreg-Paket macht dies einfach.
> countreg::rootogram(pGLM)
Die rote gebogene Linie ist die theoretische Poisson-Anpassung. An jedem Punkt der roten Linie “hängt” ein Balken, dessen Höhe die Differenz zwischen erwarteten und beobachteten Zählungen darstellt. Ein Balken, der unter 0 hängt, bedeutet eine Unteranpassung. Ein Balken, der über 0 hängt, bedeutet eine Überanpassung. Die Zählungen wurden mit einer Quadratwurzeltransformation umgewandelt, um zu verhindern, dass kleinere Zählungen durch größere Zählungen verdeckt und überlagert werden. Wir sehen eine starke Unteranpassung für Zählungen ab 2 und eine massive Überanpassung für die Zählung 1.
Wir wollen nun versuchen, ein negatives Binomialmodell anzupassen. Wir stellen fest, dass die Variabilität der Zählungen bei beiden Rassen größer ist. Es scheint, dass die negative Binomialverteilung die Verteilung der Zählungen besser wiedergibt.
Um ein negatives Binomialmodell in R anzupassen, verwenden wir die Funktion glm.nb() im Paket MASS (ein Paket, das mit R installiert wird). Auch hier zeigen wir nur einen Teil der zusammenfassenden Ausgabe:
> # 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
Zunächst fällt auf, dass die Koeffizienten dieselben sind wie zuvor. Auch hier können wir den Rennkoeffizienten potenzieren, um ein Verhältnis der Stichprobenmittelwerte zu erhalten, und Vorhersagen machen, um die ursprünglichen Stichprobenmittelwerte zu erhalten.
> # 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
Allerdings ist der Standardfehler für den Rassenkoeffizienten größer, was auf eine größere Unsicherheit in unserer Schätzung hinweist (0,24 gegenüber 0,15). Dies ist angesichts der beobachteten Variabilität in unseren Zählungen sinnvoll. Beachten Sie auch die Schätzung von Theta. Dies ist unser Dispersionsparameter. Wir können direkt von unserem Modellobjekt wie folgt darauf zugreifen:
> nbGLM$theta 0.2023119
Und wir können ihn verwenden, um geschätzte Varianzen für die Zählungen zu erhalten:
> fmeans + fmeans^2 * (1/nbGLM$theta) 1 2 0.134322 1.868928
Diese liegen viel näher an den beobachteten Varianzen als die, die das Poisson-Modell liefert.
Wir veranschaulichen die Anpassung noch einmal mit einem Rootogramm:
> countreg::rootogram(nbGLM)
Dieses sieht viel besser aus als das Rootogramm des Poisson-Modells. Es gibt eine leichte Unteranpassung/Überanpassung für die Zählungen 1 bis 3, aber ansonsten sieht es ziemlich gut aus.
Um weitere Erkenntnisse über unser negatives Binomialmodell zu gewinnen, verwenden wir seine Parameter, um Daten zu simulieren und die simulierten Daten mit den beobachteten Daten zu vergleichen. Im Folgenden laden wir das Paket magrittr für den Zugriff auf den %>%-Operator, mit dem wir Funktionen “verketten” können.
Zunächst werden die Zählungen tabelliert und ein Balkendiagramm für die weißen bzw. schwarzen Teilnehmer erstellt. Dann verwenden wir die Modellparameter, um Daten aus einer negativen Binomialverteilung zu simulieren. Die beiden Parameter sind mu und size (d. h. der Dispersionsparameter). Beachten Sie, dass wir die Funktion coef() verwenden, um die entsprechenden Koeffizienten für jede Rasse zu extrahieren. Für Weiße ist es nur der Achsenabschnitt, für Schwarze sind es der Achsenabschnitt und die Steigung (wir addieren sie also). Dann wird exponentiert, um von der logarithmischen Skala in die ursprüngliche Skala umzuwandeln. Wir simulieren die gleiche Anzahl von Beobachtungen wie in unseren Originaldaten.
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)
Die simulierten Daten sind den beobachteten Daten sehr ähnlich, was uns wiederum zuversichtlich stimmt, dass die negative Binomialregression zur Modellierung dieser Daten geeignet ist. Diese Diagramme zeigen auch die bedingte Natur unseres Modells. Die negative Binomialverteilung der Zählungen hängt von der Rasse ab bzw. ist von dieser abhängig. Jede Rasse hat einen anderen Mittelwert, aber einen gemeinsamen Streuungsparameter.
- 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
Bei Fragen oder Unklarheiten zu diesem Artikel wenden Sie sich bitte an das StatLab der UVA-Bibliothek: [email protected]
Die gesamte Sammlung der StatLab-Artikel der UVA-Bibliothek ansehen.
Clay Ford
Statistical Research Consultant
University of Virginia Library