3.1 Negativ Binomial (NB)
Die Negative Binomialverteilung (NB) kann als zweiparametrische Erweiterung der Poisson betrachtet werden, wobei die Varianz (separat vom Erwartungswert) als Gamma-verteilte Zufallsvariable betrachtet wird (Poisson-Gamma Mischmodell).
Für eine ausführliche Beschreibung, siehe Hilbe (2014) (p. 126ff).
Die Wahrscheinlichkeitsfunktion mit der \(\alpha\)-Parametriserung:
Definition 3.1 (Negative Binimialverteilung 2, alpha) Nach (Hilbe 2014 p. 130):
Für \(Y \sim NB(\mu, \alpha)\) gilt
Mit \(\mathbb{E}(Y) = \mu\) und \(\mathrm{Var}(Y) = \mu + \alpha \mu^2\)
Eine Alternative unter Verwendung von \(\theta = \frac{1}{\alpha}\):
Definition 3.2 (Negative Binomialverteilung 2, theta) Nach Perumean-Chaney u. a. (2013) (p. 1675):
Für \(Y \sim NB(\mu, \theta)\) gilt
Mit \(\mathbb{E}(Y) = \mu\) und \(\mathrm{Var}(Y) = \mu + \frac{\mu^2}{\theta}\)Der gängigste Anwendungsfall der NB findet sich bei Daten mit nicht korrigierbarer overdispersion (siehe Abschnitt 2.2.1), da der Parameter \(\theta\) bzw. auch \(\alpha = \frac 1 \theta\) als Dispersionsparameter dient:
Definition 3.3 (NB-Dispersionsparameter) Meist wird der Parameter als \(\alpha\) bezeichnet, mit der Interpretation (relativ zum Poisson-Modell):
\[\mathrm{Var}(Y) = \mu + \alpha \mu^2\]
Je nach Quelle (und unter Anderem in R (z.B. MASS::glm.nb)) wird die inverse Variante \(\theta = \frac 1 \alpha\) verwendet, womit gilt:
\[\mathrm{Var}(Y) = \mu + \frac{\mu^2}{\theta}\]
Siehe dazu auch Hilbe (2014) (p. 131).
Das GLM sieht keinen weiteren zu schätzenden Parameter vor, weshalb \(\theta\) bzw. \(\alpha\) in der Praxis separat geschätzt wird.3
Der wohl wichtigste Aspekt des Dispersionsparameters, unabhängig davon ob \(\alpha\) oder \(\theta\) verwendet wird, ist sein Vorzeichen: Er ist immer positiv, das heißt die Varianz der NB ist entweder größer oder gleich der Poisson, oder mit anderen Worten:
“The negative binomial model adjusts for Poisson overdispersion; it cannot be used to model underdispersed Poisson data”
— (Hilbe 2014 (p. 11), eigene Hervorhebung)
Je nach Software/Algorithmus ist es auch möglich, dass ein NB-fit nicht möglich ist, wenn \(\alpha \approx 0\) (die Verteilung zu nah an Poisson) bzw. die Daten Poisson-underdispersed sind.
Weiterhin gibt es zwei verschiedene Formulierungen der Negativen Binomialverteilung, NB1 und NB2, deren Nummerierung auf den Exponenten im zweiten Term ihrer Varianzen zurückzuführen ist:
Definition 3.4 (NB1 und NB2) Nach Hilbe (2014) (p. 126f):
Die NB1, oder auch lineare Negative Binomialverteilung hat die Varianz
…und die gängigere NB2, oder auch quadratische Negative Binomialverteilung hat (wie oben beschrieben) die Varianz
Beide Varianten können mittels MLE geschätzt werden, wobei NB2 zusätzlich im Kontext des GLM via IRLS geschätzt werden kann.
Als link wird \(\log \mu\) verwendet, wobei der canonical link eigentlich \(\log\left(\frac{\alpha\mu}{1 + \alpha\mu}\right)\) ist.4
Für ein gut angepasstes NB-Modell gilt analog eines Poisson-Modells, dass die Dispersionsstatistik approximativ 1 ist (siehe 2.5).
Zur Anwendung in R gibt es zwei Möglichkeiten: MASS::glm.nb oder msme::nbinomial, wobei letztere mehrere Optionen zur Parametrisierung hat und zusätzlich die Dispersionsstatistik ausgibt:
#>
#> Call:
#> MASS::glm.nb(formula = docvis ~ outwork + age, data = rwm1984,
#> init.theta = 0.4353451729, link = log)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -1.5316 -1.2816 -0.4895 0.1043 5.0468
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.039550 0.106669 -0.371 0.711
#> outwork 0.414662 0.055451 7.478 7.55e-14 ***
#> age 0.022146 0.002396 9.242 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for Negative Binomial(0.4353) family taken to be 1)
#>
#> Null deviance: 4100.8 on 3873 degrees of freedom
#> Residual deviance: 3909.6 on 3871 degrees of freedom
#> AIC: 16674
#>
#> Number of Fisher Scoring iterations: 1
#>
#>
#> Theta: 0.4353
#> Std. Err.: 0.0135
#>
#> 2 x log-likelihood: -16665.5250#>
#> Call:
#> ml_glm2(formula1 = formula1, formula2 = formula2, data = data,
#> family = family, mean.link = mean.link, scale.link = scale.link,
#> offset = offset, start = start, verbose = verbose)
#>
#> Deviance Residuals:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -1.5320 -1.2816 -0.4896 -0.4553 0.1035 5.0497
#>
#> Pearson Residuals:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -0.637147 -0.607737 -0.376257 0.000055 0.108975 21.659334
#>
#> Coefficients (all in linear predictor):
#> Estimate SE Z p LCL UCL
#> (Intercept) -0.0443 0.10570 -0.419 0.675 -0.2514 0.1629
#> outwork 0.4141 0.05532 7.485 7.15e-14 0.3056 0.5225
#> age 0.0223 0.00237 9.401 5.42e-21 0.0176 0.0269
#> (Intercept)_s 2.2971 0.07098 32.364 8.85e-230 2.1580 2.4362
#>
#> Null deviance: 4100.786 on 3872 d.f.
#> Residual deviance: 3909.522 on 3870 d.f.
#> Null Pearson: 5835.302 on 3872 d.f.
#> Residual Pearson: 5458.2 on 3870 d.f.
#> Dispersion: 1.410388
#> AIC: 16673.53
#>
#> Number of optimizer iterations: 82Literatur
Hilbe, Joseph M. 2014. Modeling Count Data. Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9781139236065.
Perumean-Chaney, Suzanne E., Charity Morgan, David McDowall, und Inmaculada Aban. 2013. „Zero-Inflated and Overdispersed: What’s One to Do?“ Journal of Statistical Computation and Simulation 83 (9): 1671–83. https://doi.org/10.1080/00949655.2012.668550.
Das ist auch der Grund, warum in R nicht
glm()verwendet werden kann, sondern eine separate Implementation mit expliziter Schätzung von \(\theta\) benötigt wird (MASS::glm.nb)↩Daher können Standardfehler nicht wie für canonical models üblich via OIM (observed information matrix) geschätzt werden, sondern via EIM (expected information matrix). Die EIM-Methode liefert weniger genaue Resultate für kleine Stichproben (n < 30). So jedenfalls Hilbe (2014).↩