3
$\begingroup$

I simulate data from a Gamma-Poisson model in R as follows. The mean and variance of the negative binomial distributed counts are $a b=10$ and $a b (1+b)=60$.

a=2
b=5
n=100000
set.seed(123)
cnt=rpois(n,rgamma(n,a,scale=b))
c(mean(cnt),var(cnt))
c(a*b,a*b*(1+b))

Next, I use two different implementations of negative binomial regression to fit parameters to the data. One is glm.nb from the MASS package and the other uses glm and the "nbinom2" family from the glmmTMB package. The documentation says the dispersion parameter is called phi and it is defined so that the variance is $V=\mu*(1+\mu/\phi) = \mu+\mu^2/\phi$. Given that the mean and variance are already known, I can set $ab(1+b) = \mu+\mu^2/\phi=ab+(ab)^2/\phi$ to find $\phi=a$.

I fit the models like this:

library(mixlm)
library(glmmTMB)
library(MASS)
glm1=glm.nb(cnt~1)
summary(glm1)
c(glm1$theta,exp(glm1$coeff))
glm2=glm(cnt~1,family=nbinom2)
summary(glm2)
c(summary(glm2)$dispersion,exp(glm2$coefficients[1]))

Then, I get these results:

Call:
glm.nb(formula = cnt ~ 1, init.theta = 2.003857425, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6795  -0.9961  -0.2781   0.4621   3.8360  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 2.304218   0.002447   941.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(2.0039) family taken to be 1)

    Null deviance: 109942  on 99999  degrees of freedom
Residual deviance: 109942  on 99999  degrees of freedom
AIC: 654400

Number of Fisher Scoring iterations: 1


              Theta:  2.0039 
          Std. Err.:  0.0108 

 2 x log-likelihood:  -654395.7440 



Call:
glm(formula = cnt ~ 1, family = nbinom2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.1906  -0.7457  -0.2058   0.3397   2.7879  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.304218   0.002443   943.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for nbinom2 family taken to be 0.5427019)

    Null deviance: 63740  on 99999  degrees of freedom
Residual deviance: 63740  on 99999  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

glm.nb gives the dispersion parameter estimate of $\hat{\theta}=2.0039$ which is close to $a$. The coefficient estimated by glm.nb is actually an estimate of the log of the mean $log(ab)\approx2.302$. On the other hand, glm estimates the same coefficient but a different dispersion parameter. It looks close to 0.5 (the inverse of 2). What is the parameter glm with nbinom2 family estimates?
In order to try to figure it out, I tried many other values of $a$ from 0.25 to 2. Then, I plotted $a$ versus $1/dispersion$

as=disps=c(25:200)/100
b=5
n=100000
set.seed(123)
for (i in 1:length(as)) {
  cnt=rpois(n,rgamma(n,as[i],scale=b))
  glm2=glm(cnt~1,family=nbinom2)
  disps[i]=summary(glm2)$dispersion
}
plot(as,1/disps)
lines(c(0,3),c(0,3))

The plot does look like the point lie on a straight line, but it's not the line $y=x$.

enter image description here

Next, I used lm to find the slope and intercept of the line and I got dispersion=0.165+0.835*a as the best fitting line. If anybody understands what the dispersion parameter estimated by this procedure is, please advise.

$\endgroup$

1 Answer 1

1
$\begingroup$

You confuse the shape $\theta$ of the Negative Binomial distribution, the dispersion parameter $\phi$ and the dispersion statistic $\hat{\phi}$. So you end up comparing two different models.

Note: The glmmTMB package can fit a Negative Binomial GLM. There is no need to mix and match (incorrectly) stats::glm with glmmTMB::nbinom2.

# It's harder to notice what's doing on with a larger sample
n <- 1000; data <- data.frame(cnt = cnt[1:n])

m.negbin <- glm.nb(cnt ~ 1, data = data)
m.glm <- glm(cnt ~ 1, family = nbinom2, data = data)
m.glmmTMB <- glmmTMB(cnt ~ 1, family = nbinom2, data = data)

Here is (part of) the summary. m.glm is different than the other two versions: the logLik function doesn't know how to compute the log likelihood, so the AIC is reported as NA, and the standard error of the intercept estimate is different.

summary(m.negbin)
#> (Dispersion parameter for Negative Binomial(2.0401) family taken to be 1)
#>
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  2.24823    0.02441   92.11   <2e-16 ***
#>
#> AIC: 6439

summary(m.glm)
#> (Dispersion parameter for nbinom2 family taken to be 0.5301128)
#>
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  2.24823    0.02421   92.87   <2e-16 ***
#>
#> AIC: NA

summary(m.glmmTMB)
#> Dispersion parameter for nbinom2 family (): 2.04 
#> 
#> Conditional model:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  2.24824    0.02441   92.11   <2e-16 ***
#>
#>      AIC      BIC   logLik deviance df.resid 
#>   6439.0   6448.8  -3217.5   6435.0      998 

The information about the dispersion is confusing. summary.negbin says the dispersion is 1 and the NB shape is 2.04; summary.glm says the dispersion is 0.53; and summary.glmmTMB says the dispersion is 2.04. What is it?

Actually there are two parameters: the shape $\theta$ of the Negative Binomial and the dispersion $\phi$. gml.nb and glmmTMB assume $\phi = 1$ and estimate $\hat{\theta} = 2.04$ while glm assumes $\theta = 1$. To verify that glm assumes $\theta = 1$, substitute family = negbin2 with family = negative.binomial(1).

The dispersion statistic on the other hand is the sum of the Pearson residuals squared divided by the residual degrees of freedom. It's an estimate of $\phi$.

See also dispersion in summary.glm().

# The shape theta of the Negative Binomial distribution
c(
  m.negbin$theta,
  1,
  sigma(m.glmmTMB)
)
#> [1] 2.040085 1.000000 2.040087

# The residual degrees of freedom
c(
  df.residual(m.negbin),
  df.residual(m.glm),
  df.residual(m.glmmTMB)
)
#> [1] 999 999 998

# Compute the dispersion statistic
dispersion <- function(m) {
  sum(residuals(m, type = "pearson")^2) / df.residual(m)
}

c(
  dispersion(m.negbin),
  dispersion(m.glm), # same as summary(m.glm)$dispersion
  dispersion(m.glmmTMB) * (n - 2) / (n - 1)
)
#> [1] 0.9837584 0.5301128 0.9837538
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.