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

\[\begin{equation*} P(Y = y_i) = \begin{pmatrix} y_i + \frac{1}{\alpha} - 1 \\ \frac{1}{\alpha} - 1\end{pmatrix} \left(\frac{1}{a + \alpha \mu_i}\right)^{\frac{1}{\alpha}} \left(\frac{\alpha \mu_i}{1+\alpha\mu_i}\right) \end{equation*}\]

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

\[\begin{equation*} P(Y = y) = \frac{\Gamma(\theta + y)}{\Gamma(\theta) \Gamma(y+1)} \frac{\theta^\theta \mu^y}{(\theta + \mu)^{(\theta + y)}}, \quad y = 0, 1, 2, \ldots \end{equation*}\]

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\]

\[\begin{align*} \alpha &= 0 &&\Longrightarrow \text{ Äquivalent zur Poissonverteilung} \\ \alpha &> 0 &&\Longrightarrow \text{ Overdispersion} \end{align*}\]

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}\]

\[\begin{align*} \theta &> 0 &&\Longrightarrow \text{ Overdispersion} \\ \theta &\to \infty &&\Longrightarrow \text{ Äquivalent zur Poissonverteilung} \end{align*}\]

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

\[\begin{equation*} \mathrm{Var}(Y) = \mu + \alpha \mu \end{equation*}\]

…und die gängigere NB2, oder auch quadratische Negative Binomialverteilung hat (wie oben beschrieben) die Varianz

\[\begin{equation*} \mathrm{Var}(Y) = \mu + \alpha \mu^2 \end{equation*}\]

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:

mod_nb <- MASS::glm.nb(docvis ~ outwork + age,
                       data = rwm1984)
summary(mod_nb)
#> 
#> 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
mod_nb2 <- msme::nbinomial(docvis ~ outwork + age,
                           data = rwm1984)
summary(mod_nb2)
#> 
#> 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:  82

Literatur

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.