GR’s Website

4.A Modelle für überdispergierte Zähldaten

Wir verwenden die Daten von Long (1990) über die Anzahl der Veröffentlichungen von promovierten Biochemikern, um die Anwendung von Poisson-, überdispergierten Poisson-, negativen Binomial- und Null-inflationären Poisson-Modellen zu illustrieren.

Die Variablen im Datensatz sind

kodiert eins, wenn verheiratet

Prestige des Ph.D. Programm

Variable	Beschreibung
art	Artikel in den letzten drei Jahren der Promotion.
fem	kodiert eins für Frauen
mar
kid5	Anzahl der Kinder unter sechs Jahren
phd
ment	Artikel durch Mentor in den letzten drei Jahren

Diese Daten wurden auch von Long und Freese (2001) analysiert und sind auf der Stata-Website verfügbar:

> library(foreign)> ab names(ab) "art" "fem" "mar" "kid5" "phd" "ment"> r c(mean=r, var=r, ratio=r/r) mean var ratio 1.692896 3.709742 2.191358 

Die durchschnittliche Anzahl der Artikel beträgt 1,69 und die Varianz 3,71, also etwas mehr als das Doppelte des Mittelwertes. Die Daten sind überdispers, aber natürlich haben wir noch keine Kovariaten berücksichtigt.

Ein Poisson-Modell

Wir wollen das von Long und Freese (2001) verwendete Modell anpassen, ein einfaches additives Modell, das alle fünf Prädiktoren verwendet.

> mp summary(mp)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = poisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.304617 0.102981 2.958 0.0031 ** femWomen -0.224594 0.054613 -4.112 3.92e-05 ***marMarried 0.155243 0.061374 2.529 0.0114 * kid5 -0.184883 0.040127 -4.607 4.08e-06 ***phd 0.012823 0.026397 0.486 0.6271 ment 0.025543 0.002006 12.733 
We see that the model obviously doesn't fit the data. The five-percent critical value for a chi-squared with 909 d.f. is
>qchisq(0.95, df.residual(mp)) 980.2518> deviance(mp) 1634.371> pr sum(pr^2) 1662.547
und die Abweichung und das Chi-Quadrat von Pearson liegen beide in den 1600ern.
 4.A.1. Extra-Poisson-Variation
Wir nehmen nun an, dass die Varianz proportional und nicht gleich dem Mittelwert ist, und schätzen den Skalenparameter φ, indem wir das Pearson-Chi-Quadrat durch seinen d.f. teilen:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Wir sehen, dass die Varianz etwa 83% größer ist als der Mittelwert. Das bedeutet, dass wir die Standardfehler anpassen sollten, indem wir mit 1,35, der Quadratwurzel aus 1,83, multiplizieren.
R kann diese Berechnung für uns durchführen, wenn wir die quasipoissonFamilie verwenden:
> mq summary(mq)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = quasipoisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.304617 0.139273 2.187 0.028983 * femWomen -0.224594 0.073860 -3.041 0.002427 ** marMarried 0.155243 0.083003 1.870 0.061759 . kid5 -0.184883 0.054268 -3.407 0.000686 ***phd 0.012823 0.035700 0.359 0.719544 ment 0.025543 0.002713 9.415 
The estimates are exactly the same as before, but the standard errorsare about 35% larger. We can verify this fact easily. First we writea useful function to extract standard errors and then use it onour fits:
> se round(data.frame(p=coef(mp), q=coef(mq),+ se.p=se(mp), se.q=se(mq), ratio=se(mq)/se(mp)), 4) p q se.p se.q ratio(Intercept) 0.3046 0.3046 0.1030 0.1393 1.3524femWomen -0.2246 -0.2246 0.0546 0.0739 1.3524marMarried 0.1552 0.1552 0.0614 0.0830 1.3524kid5 -0.1849 -0.1849 0.0401 0.0543 1.3524phd 0.0128 0.0128 0.0264 0.0357 1.3524ment 0.0255 0.0255 0.0020 0.0027 1.3524
Eine alternative Herangehensweise besteht darin, ein Poisson-Modell anzupassen und den robusten oder Sandwich-Schätzer der Standardfehler zu verwenden. Dies führt in der Regel zu Ergebnissen, die dem überdispersen Poisson-Modell sehr ähnlich sind.
 4.A.2 Negative Binomialregression
Wir passen nun ein negatives Binomialmodell mit denselben Prädiktoren an.
> library(MASS)> mnb summary(mnb)Call:glm.nb(formula = art ~ fem + mar + kid5 + phd + ment, data = ab, init.theta = 2.264387695, link = log)Deviance Residuals: Min 1Q Median 3Q Max -2.1678 -1.3617 -0.2806 0.4476 3.4524 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.256144 0.137348 1.865 0.062191 . femWomen -0.216418 0.072636 -2.979 0.002887 ** marMarried 0.150489 0.082097 1.833 0.066791 . kid5 -0.176415 0.052813 -3.340 0.000837 ***phd 0.015271 0.035873 0.426 0.670326 ment 0.029082 0.003214 9.048 1/mnb$theta 0.4416205> -2*(logLik(mp)-logLik(mnb))'log Lik.' 180.196 (df=6)
Dazu benötigen wir die Funktion glm.nb() im PaketMASS.
> library(MASS)> mnb summary(mnb)Call:glm.nb(formula = art ~ fem + mar + kid5 + phd + ment, data = ab, init.theta = 2.264387695, link = log)Deviance Residuals: Min 1Q Median 3Q Max -2.1678 -1.3617 -0.2806 0.4476 3.4524 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.256144 0.137348 1.865 0.062191 . femWomen -0.216418 0.072636 -2.979 0.002887 ** marMarried 0.150489 0.082097 1.833 0.066791 . kid5 -0.176415 0.052813 -3.340 0.000837 ***phd 0.015271 0.035873 0.426 0.670326 ment 0.029082 0.003214 9.048 1/mnb$theta 0.4416205> -2*(logLik(mp)-logLik(mnb))'log Lik.' 180.196 (df=6)
R's theta ist die Genauigkeit des multiplikativen Zufallseffekts und entspricht 1/σ2 in den Anführungszeichen. Die Schätzung entspricht einer geschätzten Varianz von 0,44 und ist hochsignifikant.
Um die Signifikanz dieses Parameters zu testen, kann man sich vorstellen, die doppelte Differenz der Log-Likelihoods zwischen diesem Modell und dem Poisson-Modell, 180,2, zu berechnen und sie als Chi-Quadrat mit einer d.f. zu behandeln. Die übliche Asymptotik gilt jedoch nicht, da die Nullhypothese an einer Grenze des Parameterraums liegt.
Es gibt einige Arbeiten, die zeigen, dass eine bessere Annäherung darin besteht, die Statistik als 50:50-Mischung aus Null und einem Chi-Quadrat mit einer d.f. zu behandeln. 
Für das Testen von Hypothesen über die Regressionskoeffizienten können wir entweder Wald-Tests oder Likelihood-Ratio-Tests verwenden, die möglich sind, weil wir vollständige Verteilungsannahmen getroffen haben.
Es gibt auch eine negative.binomial-Familie für glm, und diese kann verwendet werden, sofern die Parametertheta gegeben sind. (Dies basiert auf dem Ergebnis, dass das negative Binomial in der glm-Familie für feste Varianz ist.) Hier eine kurze Überprüfung:
> mnbg all(abs(coef(mnb)==coef(mnbg)) 
The standard errors would differ, however, because  allows for the fact that was estimated, whereas  does not.
Unobserved Heterogeneity
R has a function  to compute the density of a gamma distribution with given shape and scale (or its reciprocal the rate). In particular, when the random effect hasvariance  the density is. This can beused to draw the density.
> v = 1/mnb$theta> gd = data.frame(x = seq(0.0, 3, 0.1), g = dgamma(x, shape = 1/v, scale = v))> ggplot(gd, aes(x, g)) + geom_line()> ggsave("gamdenr.png", width=500/72, height=400/72, dpi=72)

Wir können auch Quantile mitqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE) berechnen. Insbesondere, wenn der zufällige Effekt eine Varianz von v hat, sind die Quartile qgamma((1:3)/4, shap e= 1/v, scale = v).
> qgamma((1:3)/4, shape = 1/v, scale = v) 0.5114167 0.8572697 1.3347651
Biochemiker in Q1 der Verteilung der unbeobachteten Heterogenität veröffentlichen 49% weniger Arbeiten als aufgrund ihrer beobachteten Merkmale zu erwarten wäre, während diejenigen im Median 14% weniger und diejenigen in Q3 33% mehr als erwartet veröffentlichen.
Vergleich der Schätzungen und Standardfehler
Vergleichen wir die Parameterschätzungen und Standardfehler des Poisson-Modells, des überdispersen Poisson-Modells und des negativen Binomialmodells:
> round(data.frame(+ p=coef(mp),q=coef(mq),nb=coef(mnb),+ se.p=se(mp),se.q=se(mq),se.nb=se(mnb)),4) p q nb se.p se.q se.nb(Intercept) 0.3046 0.3046 0.2561 0.1030 0.1393 0.1373femWomen -0.2246 -0.2246 -0.2164 0.0546 0.0739 0.0726marMarried 0.1552 0.1552 0.1505 0.0614 0.0830 0.0821kid5 -0.1849 -0.1849 -0.1764 0.0401 0.0543 0.0528phd 0.0128 0.0128 0.0153 0.0264 0.0357 0.0359ment 0.0255 0.0255 0.0291 0.0020 0.0027 0.0032
Die Schätzungen des negativen Binomialmodells unterscheiden sich nicht wesentlich von denen des Poisson-Modells, und beide Schätzungen würden zu denselben Schlussfolgerungen führen. 
Betrachtet man die Standardfehler, so stellt man fest, dass beide Ansätze zur Überstreuung zu sehr ähnlichen geschätzten Standardfehlern führen, und dass die gewöhnliche Poisson-Regression die Standardfehler unterschätzt.
Anpassungsgüte
Wir können die Anpassungsgüte des negativen Binomialmodells anhand der Abweichung
> deviance(mnbg) 1004.281
Das negative Binomialmodell passt besser als das Poissonmodell, hat aber immer noch eine Abweichung über dem kritischen Fünf-Prozent-Wert von980,25.
Die Varianzfunktion
Das überstreute Poissonmodell und das negative Binomialmodell haben unterschiedliche Varianzfunktionen. Eine Möglichkeit zu prüfen, welches Modell besser geeignet ist, besteht darin, Gruppen auf der Grundlage des linearen Prädiktors zu erstellen, den Mittelwert und die Varianz für jede Gruppe zu berechnen und schließlich die Mittelwert-Varianz-Beziehung darzustellen.
Hier sind Gruppen auf der Grundlage des linearen Prädiktors mit negativem Binomialmodell zu sehen, die mit cut mit Unterbrechungen bei den 5(5)95-Perzentilen erstellt wurden, um zwanzig Gruppen von annähernd gleicher Größe zu erzeugen. Beachten Sie, dass predict() einen linearen Prädiktor berechnet, sofern nicht anders angegeben. Für eine Vorhersage auf der Skala der Antwort fügen Sie die Option type="response" hinzu. 
> xb g m v png("c4afig1.png",width=500,height=400)> plot(m, v, xlab="Mean", ylab="Variance", + main="Mean-Variance Relationship")> mtext("Articles Published by Ph.D. Biochemists",padj=-0.5)> x lines(x, x*phi, lty="dashed")> lines(x, x*(1+x/mnb$theta))> legend("topleft", lty=c("dashed","solid"), + legend=c("Q. Poisson","Neg. Binom."), inset=0.05)> dev.off();null device 1 

Das Diagramm stellt den Mittelwert der Varianz gegenüber und überlagert die Kurven, die dem überstreuten Poisson-Modell, bei dem die Varianz φμ ist, und dem negativen Binomialmodell, bei dem die Varianz μ(1+μσ2) ist, entsprechen.
Die Poisson-Varianzfunktion ist für den Großteil der Daten recht gut geeignet, erfasst jedoch nicht die hohen Varianzen der produktivsten Wissenschaftler. Die negative binomiale Varianzfunktion unterscheidet sich nicht allzu sehr von der Poisson-Funktion, kann aber als quadratische Funktion schneller ansteigen und leistet am oberen Ende bessere Arbeit. Wir kommen zu dem Schluss, dass das negative Binomialmodell die Daten besser beschreibt als das überstreute Poisson-Modell.
 4.A.3 Null-Inflations-Poisson-Modelle
 Ein häufiges Phänomen bei Zähldaten ist ein Übermaß an Nullen im Vergleich zu dem, was bei einem Poisson-Modell erwartet wird. Dies ist bei unseren Daten tatsächlich ein Problem: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Wir sehen, dass 30,0 % der Wissenschaftler in der Stichprobe in den letzten drei Jahren nach ihrer Promotion keine Artikel veröffentlicht haben, aber das Poisson-Modell sagt voraus, dass nur 20,9 % keine Veröffentlichungen haben würden. Das Modell unterschätzt eindeutig die Wahrscheinlichkeit von Nullzählungen.
Eine Möglichkeit, diese Art von Situation zu modellieren, besteht darin, anzunehmen, dass die Daten aus einer Mischung von zwei Populationen stammen, von denen die eine immer Null zählt und die andere eine Poisson-Verteilung mit Mittelwert μ aufweist. In diesem Modell können Nullzählungen aus beiden Populationen stammen, während positive Zählungen nur aus der zweiten Population stammen.
Im Zusammenhang mit den Veröffentlichungen von promovierten Biochemikern können wir uns vorstellen, dass einige eine Stelle im Auge hatten, bei der Veröffentlichungen nicht wichtig sind, während andere eine akademische Stelle anstrebten, bei der eine Vielzahl von Veröffentlichungen erwartet wurde. Die Mitglieder der ersten Gruppe würden null Artikel veröffentlichen, während die Mitglieder der zweiten Gruppe 0,1,2,... veröffentlichen würden, eine Anzahl, die als Poisson-Verteilung angenommen werden kann.
Die Verteilung des Ergebnisses kann dann in Form von zwei Parametern modelliert werden, π, der Wahrscheinlichkeit von "immer Null", und μ, der mittleren Anzahl von Veröffentlichungen für diejenigen, die nicht in der "immer Null"-Gruppe sind. Ein natürlicher Weg, um Kovariaten einzuführen, ist die Modellierung des Logits der Wahrscheinlichkeit π von "immer Null" und des Logs des Mittelwerts μ für diejenigen, die nicht in der "immer Null"-Klasse sind.
Diese Art von Modell kann in R mit der Funktion zeroinfl() im Paket pscl angepasst werden. Die Modellformel kann wie üblich angegeben werden, wenn die gleichen Variablen in beiden Gleichungen enthalten sein sollen. Andernfalls kann man zwei Sätze von Prädiktoren angeben, die durch einen vertikalen Balken getrennt sind. Geben Sie ?zeroinfl ein, um mehr darüber zu erfahren.
Hier ist ein null-inflatiertes Modell mit allen Kovariaten in beiden Gleichungen:
> library(pscl)> mzip summary(mzip)Call:zeroinfl(formula = art ~ fem + mar + kid5 + phd + ment, data = ab)Pearson residuals: Min 1Q Median 3Q Max -2.3253 -0.8652 -0.2826 0.5404 7.2976 Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.640838 0.121307 5.283 1.27e-07 ***femWomen -0.209145 0.063405 -3.299 0.000972 ***marMarried 0.103751 0.071111 1.459 0.144565 kid5 -0.143320 0.047429 -3.022 0.002513 ** phd -0.006166 0.031008 -0.199 0.842378 ment 0.018098 0.002294 7.888 3.07e-15 ***Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.577059 0.509386 -1.133 0.25728 femWomen 0.109746 0.280082 0.392 0.69518 marMarried -0.354014 0.317611 -1.115 0.26502 kid5 0.217097 0.196482 1.105 0.26919 phd 0.001274 0.145263 0.009 0.99300 ment -0.134114 0.045243 -2.964 0.00303 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 21 Log-likelihood: -1605 on 12 Df
Betrachten wir die inflate-Gleichung, so sehen wir, dass der einzige signifikante Prädiktor für die Zugehörigkeit zur Klasse "immer Null" die Anzahl der vom Mentor veröffentlichten Artikel ist, wobei jeder Artikel des Mentors mit einer um 12,6 % geringeren Wahrscheinlichkeit verbunden ist, nie zu veröffentlichen.
Betrachten wir die Gleichung für die durchschnittliche Anzahl der Artikel unter denjenigen, die nicht in der "Immer-Null"-Klasse sind, so finden wir signifikante Nachteile für Frauen und Wissenschaftler mit Kindern unter fünf Jahren und einen großen positiven Effekt der Anzahl der Veröffentlichungen durch den Mentor, wobei jeder Artikel mit einer 1,8 %igen Erhöhung der erwarteten Anzahl der Veröffentlichungen verbunden ist.
Um zu überprüfen, ob das Modell das Problem der überzähligen Nullen löst, sagen wir π und μ voraus und berechnen die kombinierte Wahrscheinlichkeit, dass keine Veröffentlichungen stattfinden. Es gibt Optionen in der predict()Funktion namens "zero" und "count", um diese zu erhalten. Es gibt auch eine Option "prob" zur Berechnung der vorhergesagten Dichte, aber das ist übertrieben, da wir nur die Wahrscheinlichkeit von Null haben wollen.
> pr mu zip mean(zip) 0.2985679
Das Modell löst also das Problem der überzähligen Nullen und sagt voraus, dass 29,9 % der Biochemiker keine Artikel veröffentlichen werden, was dem beobachteten Wert von 30,0 % sehr viel näher kommt.
Es gibt auch ein negatives Binomialmodell mit Nullen, das ein negatives Binomial für die Anzahl in der Klasse "nicht immer Null" verwendet. Dieses Modell kann mit zeroinfl() und dem Parameter dist="negbin" angepasst werden. Zu den alternativen Verknüpfungen für die Inflationsgleichung gehört das Probit-Modell, das mit link="probit" angegeben werden kann. 
Modellvergleich mit AIC
Auch für diese Daten löst das negative Binomial das Problem! Hier ist die Wahrscheinlichkeit von Nullartikeln im negativen Binomialmodell. Wir gehen von den ersten Grundsätzen aus, aber man könnte auch die eingebaute negative Binomialdichte verwenden
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Das Modell sagt voraus, dass 30,4 % der Biochemiker in den letzten drei Jahren ihrer Promotion keine Artikel veröffentlichen würden, was dem beobachteten Wert von 30,0 % sehr nahe kommt.
Um zwischen dem negativen Binomialmodell und dem auf Null aufgeblasenen Modell zu wählen, müssen wir auf andere Kriterien zurückgreifen. Eine sehr einfache Möglichkeit, Modelle mit einer unterschiedlichen Anzahl von Parametern zu vergleichen, ist die Berechnung des Akaike-Informationskriteriums (AIC), das wir wie folgt definieren: 
AIC = -2logL + 2p
wobei p die Anzahl der Parameter im Modell ist. Der erste Term ist im Wesentlichen die Abweichung und der zweite eine Strafe für die Anzahl der Parameter, und er kann für die meisten Modelle mit der Funktion AIC berechnet werden. Ich ermittle ihn auch "von Hand", damit wir die Berechnung überprüfen können. Für unsere Daten
> c(nb=AIC(mnb), zip=AIC(mzip)) nb zip 3135.917 3233.546 > mzip$rank aic sapply(list(mnb,mzip), aic) 3133.917 3233.546
 Für diesen Datensatz ist das negative Binomialmodell ein klarer Gewinner in Bezug auf Parsimonie und Anpassungsgüte. 
Weitere diagnostische Kriterien, die wir uns ansehen könnten, sind die Randverteilung der vorhergesagten und beobachteten Zählungen und die Varianzfunktionen. 
Zero-Truncated- und Hurdle-Modelle
Andere Modelle, die wir noch nicht behandelt haben, sind das null-trunkierte Poisson- und das negative Binomialmodell, die für Daten entwickelt wurden, die keine Nullen enthalten. Ein gängiges Beispiel ist die Aufenthaltsdauer in einem Krankenhaus, die mindestens einen Tag beträgt. Ein sinnvoller Ansatz ist die Anpassung eines Poisson- oder negativen Binomialmodells, das Nullen ausschließt und die anderen Wahrscheinlichkeiten so skaliert, dass sie in der Summe eins ergeben. Bei der Interpretation dieser Modelle ist Vorsicht geboten, denn μ ist nicht das erwartete Ergebnis, sondern der Mittelwert einer zugrundeliegenden Verteilung, in der die Nullen enthalten sind. Diese Modelle sind in der Funktion vglm() des Pakets VGAM unter Verwendung der Familien pospoisson und posnegbinomial implementiert.
Eine alternative Herangehensweise an einen Überschuss (oder einen Mangel) an Nullen ist die Verwendung eines zweistufigen Prozesses, mit einem Logit-Modell zur Unterscheidung zwischen Nullen und positiven Zählungen und dann einem Poisson- oder Negativ-Binomial-Modell mit Nullabbruch für die positiven Zählungen. In unserem Beispiel könnten wir ein Logit-Modell verwenden, um zwischen denjenigen zu unterscheiden, die veröffentlichen, und denjenigen, die nicht veröffentlichen, und dann ein Poisson- oder Negativ-Binomial-Modell mit Abbruch für die Anzahl der Artikel derjenigen, die mindestens einen Artikel veröffentlichen. Diese Modelle werden oft als Hürdenmodelle bezeichnet. Sie können in R unter Verwendung der separaten Logit- und Null-abgeschnittenen Poisson- oder Negativ-Binomial-Modelle und einfacher Addition der Log-Likelihoods oder unter Verwendung der Funktion hurdle() im Paket pscl angepasst werden.
Vergleicht man Hurdle- und Zero-inflated-Modelle, so finde ich, dass die Unterscheidung zwischen Null und Eins oder mehr bei Hurdle-Modellen klarer ist, aber die Interpretation des Mittelwerts ist bei Zero-inflated-Modellen klarer.
Das folgende JSS-Papier enthält eine nützliche Diskussion über alle diese Modelle: Regression Models for Count Data in R.