4.1 Overdispersion: Erkennung und Handhabung
Es gibt mehrere Möglichkeiten Poisson-overdispersion zu erkennen, wobei die Pearson-Dispersionsstatistik in der Regel der erste Schritt ist.
Im Folgenden werden zusätzlich einige formale Tests aufgeführt.
4.1.1 Lagrange Mulitplier Test
Definition 4.1 (Lagrange Multiplier Test) Nach Hilbe (2014) (p. 85f):
Eine \(\chi^2\)-Teststatistik mit einem Freiheitsgrad.
Eine rudimentäre R-Implementation findet sich in Abschnitt 1.2.2.
4.1.2 (Boundary) Likelihood Ratio Test (BLR)
Zum vergleich von zwei geschachtelten (nested) Modellen kann der Likelihood Ratio Test verwendet werden:
Definition 4.2 (Likelihood Ratio Test) \[\begin{equation} LR = -2 (\mathcal{L}_R - \mathcal{L}_F) \end{equation}\]
Mit \(\mathcal{L}_F\) als log-likelihood des vollen (oder “größeren”) Modells, und \(\mathcal{L}_R\) als log-likelihood des reduzierteren Modells.Eine Variante des Tests kann verwendet werden, um den Dispersionsparameter \(\alpha\) eines NB-Modells zu testen. Da eine NB-Verteilung für \(\alpha = 0\) äquivalent zur Poisson ist (siehe Abschnitt 3.1), kann ein Poisson-Modell als reduzierte Variante eines NB-Modells betrachtet werden. In diesem Fall verwendet man den Boundary Likelhihood Ratio test:
It is important, though, to remember that the BLR test has a lower limiting case for the value of \(\alpha\), which is what is being tested. Given that the standard parameterization of the negative binomial variance function is \(\mu + \alpha \mu^2\) , when \(\alpha = 0\), the variance reduces to \(\mu\).
— Hilbe (2014), p. 115
Der resultierende Wert ist \(\chi^2_{(1)}\)-verteilt. Der resultierende p-Wert muss zusätzlich halbiert werden (siehe Hilbe 2014, p. 115).
Since the distribution being tested can go no lower than 0, that is the boundary. Only one half of the full distribution is used. Therefore the Chi2 (sic!) test is divided by 2
— Hilbe (2014), p. 115
Am Beispiel der rwm1984-Daten:
# Poisson-Modell
mod_p <- glm(docvis ~ outwork + age, family = poisson(), data = rwm1984)
# NB-Modell
mod_nb <- MASS::glm.nb(docvis ~ outwork + age, data = rwm1984)
# BLR: Recht eindeutig.
lmtest::lrtest(mod_p, mod_nb)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 3 | -15636.390 | NA | NA | NA |
| 4 | -8332.762 | 1 | 14607.26 | 0 |
# Bei gegebenem BLR-Wert von 4.2 ließe sich der p-Wert wie folgt berechnen:
pchisq(4.2, df = 1, lower.tail = FALSE) / 2#> [1] 0.02021199#> [1] 0.020211994.1.3 Umgang mit Overdispersion
Zur expliziten Modellierung des Dispersionparameters kann ein NBH Modell (3.5.1) genutzt werden, falls die Quelle der overdispersion von besonderem Interesse ist.
Abseits davon bleiben 2 grobe Strategien :
- Adjustierung der durch die overdispersion verzerrten Standardfehler (Quasipoisson, Quasi-Likelihood, robuste Varianzschätzer)
- Wechsel auf ein Modell, das overdispersion (oder allgemeine extradispersion) erlaubt (e.g. NB)
(Bemerke: Standardfehler für IRRs werden i.A. über die delta method bestimmt, die ich noch recherchieren und irgendwo einbauen muss)
Quasipoisson: Skalierung der Standardfehler
Hier werden lediglich die Standardfehler der Koeffizienten eines Poisson-Modells adjustiert, die bei overdispersion i.d.R. unterschätzt werden – die eigentliche Parameterschätzung wird nicht beeinflusst:
\[\mathrm{SE}(\beta_k) \cdot \sqrt{D}\] Wobei \(D\) der Pearson-Dispersionsindex ist (vgl. Hilbe (2014), p. 92ff), den wir in 2.5 als \(\frac{\chi^2_{\mathrm{Pearson}}}{\mathrm{df}}\) definiert hatten.
Ein möglicher Nachteil jedoch: Es wird zuerst ein reguläres Poisson-Modell gefittet, Standardfehler und Dispersionsindex bestimmt, und dann dasselbe Modell mit skalierten Standardfehlern erneut gefittet – dementsprechend ist diese Methode vermutlich für größere Datensätze eher ineffizient.
Einmal auf unser voriges Beispiel anhand der rwm5yr-Daten angewandt:
# Urprüngliches Model
mod <- glm(docvis ~ outwork + age, family = poisson(), data = rwm1984)
mod_qp <- glm(docvis ~ outwork + age, family = quasipoisson(), data = rwm1984)
dispersion(mod)#> X-squared(3871) = 43909.62
#> Pearson Dispersion = 11.343#> X-squared(3871) = 43909.62
#> Pearson Dispersion = 11.343An der Dispersion hat sich nichts geändert, nur an den Standardfehlern der Koeffizienten (in folgenden Tabellen auf log-Skala):
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -0.03352 | 0.03918 | -0.8554 | 0.3923 |
| outwork | 0.4079 | 0.01884 | 21.65 | 6.481e-104 |
| age | 0.02208 | 0.0008377 | 26.36 | 3.73e-153 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -0.03352 | 0.132 | -0.254 | 0.7995 |
| outwork | 0.4079 | 0.06347 | 6.427 | 1.457e-10 |
| age | 0.02208 | 0.002821 | 7.827 | 6.386e-15 |
In diesem Fall “lohnt” sich dieser Ansatz zwar eher nicht, da wir neben overdispersion noch ein Problem mit zero-inflation haben, aber in manchen Fällen kann Quasipoisson-Methode ausreichen.
“Scaling is an easy, quick-and-dirty method of adjusting standard errors for overdispersion. However, when data are highly correlated or clustered, model slopes or coefficients usually need to be adjusted as well. Scaling does not accommodate that need, but simulations demonstrate that it is quite useful for models with little to moderate overdispersion.”
— Hilbe (2014) (p. 96)
4.1.3.1 Robuste Varianzschätzer (Sandwich Estimators)
(Andere Namen: Huber & White standard errors, empirical standard errors)
Sandwich estimators sind in erster Linie für nicht unabhängige bzw. korrellierte Daten gedacht, zum Beispiel für Daten, die innerhalb mehrere Haushalte, Krankenhäuser, Städte etc. gesammelt wurden.
Robuste Varianzschätzer können (und je nach Quelle: sollten) standardmäßig für Count-response Daten verwendet werden, da die resultierenden Schätzer im Falle tatsächlich unkorrellierter Daten äquivalent zu den Standardfehlern des ursprünglichen Modells sind (i.e. “schadet nicht”).
Bootstrapped SEs wiederum erfordern mehr Aufwand, ähneln aber den robusten SEs sowieso sehr stark, weshalb sich der Aufwand ggf. nicht wirklich lohnt.
Siehe Hilbe (2014) (p. 99f) für eine ausführlichere Beschreibung. Zudem gibt der Autor explizit den Rat:
“Unless your Poisson or negative binomial model is well fitted and meets its respective distributional assumptions, use robust or empirical standard errors as a default”
— Hilbe (2014) (p. 133)
Literatur
Hilbe, Joseph M. 2014. Modeling Count Data. Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9781139236065.