dispersion of a negative binomial model

Question

In R's, glm.nb summary, it says dispersion parameter $\phi$ is set to 1. When the model is

$Y \sim \text{Negbin}(\mu,\theta)$

where $E(Y)=\mu$ and $V(Y)=\mu+\mu^2/\theta$. Is it correct that dispersion $\phi$ here means $V(Y)=\phi(\mu+\mu^2/\theta)$ in the exponential family formulation? Given that $\theta$ can be freely adjusted, does it make sense to test, e.g., $H_0: \phi=1$? I know negative binomial models may not fit well to data for various reasons (e.g., relationship between many zeros and other values), but if we say it is overdispersion, we also need to consider overdispersion for a normal distribution, but it is often stated there is no overdispersion in gaussian models.

set.seed(1)
x <- round(rep(seq(3,30,by=3),each=10))
y <- rnbinom(length(x),mu=exp(1+0.1*x),size=5)

plot(x,y)
library(MASS)
model <- glm.nb(y~x)
r <- resid(model,type="pearson")
(phi <- sum(r^2)/(length(x)-3)) # dispersion

Similar post: link

Answer 1

Your primary question seems to be this one:

Given that $\theta$ can be freely adjusted, does it make sense to test, e.g., $H_0: \phi = 1$?

Two methods one can use for the Poisson are the score test and Lagrange multiplier tests. I might be wrong, but I imagine these can be generalized to a degree to an negative binomial model. A score test is a post-hoc test which tests whether the hypothesis of no overdispersion can be rejected:

$$ z = \frac{(y - \mu)^2 - y}{\mu\sqrt{2}} $$

where $\mu$ is the predicted response and $y$ is the raw response. You basically fit a regular NB model, take the predictions from that model, and then solve for $z$:

#### Load Data ####
library(COUNT)
data(rwm5yr)
rwm1984 <- subset(rwm5yr, year==1984) 

#### Fit NB Model ####
nb <- glm.nb(
  docvis ~ outwork + age, 
  data=rwm1984
)

#### Get Mu ####
mu <- predict(nb, type="response")

#### Estimate Z ####
z <- ((rwm1984$docvis - mu)^2 - rwm1984$docvis)/ (mu * sqrt(2))
summary(zscore <- lm(z ~ 1))

The t-test from the intercept ($t = 7.082$, $p < .001$) shows the "no dispersion" hypothesis can be rejected. However, this test typically runs with the assumption that the dataset is large. One can alternatively use the Lagrange multiplier test instead, which is defined as:

$$ \chi^2 = \frac{\left(\sum_{i=1}^{n} \mu^2_i-n \bar{y}_i\right)^2}{2 \sum_{i=1}^{n} \mu^2_i} $$

Estimated as so:

#### Setup Formula for Chi ####
obs <- nrow(rwm1984)
mmu <- mean(mu)
nybar <- obs*mmu
musq <- mu*mu
mu2 <- mean(musq)*obs 

#### Chi Square Estimation ####
chival <- (mu2 - nybar)^2/(2*mu2)
chival

#### P Value ####
pchisq(chival,1,lower.tail = FALSE)

Once again, the hypothesis is rejected, with a fairly large chi-square score. Keep in mind that if you just want to get the dispersion from the model, this function can do that quite easily, and was created specifically for glm.nb fits in the msme package but written explicitly here:

#### Create Dispersion Function ####
p_disp <- function(x) {
  pr <- sum(residuals(x, type = "pearson")^2)
  dispersion <- pr/x$df.residual
  return(c(pearson.chi2 = pr, dispersion = dispersion))
}

#### Check ####
p_disp(nb)

Where you can see the dispersion is clearly above 1. More about these tests can be found in the below reference.

Reference

Hilbe, J. M. (2014). Modeling count data. Cambridge Univ. Press.

Answer 2

we also need to consider overdispersion for a normal distribution, but it is often stated there is no overdispersion in gaussian models.

This can be related to the following difference:

• For the normal distribution, there is intrinsically a dispersion parameter, the squared standard deviation $\sigma^2$.
• For the negative binomial distribution there is not such a parameter. The plain negative binomial distribution doesn't have a dispersion parameter.

The addition of a dispersion parameter is artificial. In the framework of generalized linear models if the score function (which is used to compute the maximum likelihood) does not yet have a dispersion parameter, then we can add an additional dispersion parameter.

With this adjustment, the resulting score function does not need to relate to an actual existing distribution and likelihood function. It is called a quasi-likelihood function.

For those cases, when this dispersion parameter is not normally part of the actual distribution (e.g. the binomial distribution, Poisson distribution and negative binomial distribution), there we call it over-dispersion when the variance is larger than expected. A good example case of the opposite, under-dispersion is: A chart of daily cases of COVID-19 in a Russian region looks suspiciously level to me - is this so from the statistics viewpoint?

---

Where does this dispersion parameter fit in?

For many GLM models, the score is a simple function of the difference between the observation and the mean, divided by a function for the variance that is a function solely of the mean.

$$s(y;\mu) = \frac{\partial \log[\mathcal{L}(y,\mu]}{\partial \mu} = \frac{y-\mu}{V(\mu)}$$

Several GLM models can have an additional dispersion parameter, then the score looks like

$$s(y;\mu) = \frac{\partial \log[\mathcal{L}(y,\mu]}{\partial \mu} = \frac{y-\mu}{\phi V(\mu)}$$

For distributions like Poisson or Binomial distribution this parameter is not in the model and we can say $\phi = 1$. For example,

• with a Bernoulli distribution $V(\mu) = \mu(1-\mu)$ and $$s(y;\mu) = \frac{y-\mu}{\mu(1-\mu)}$$
• with a Poisson distribution $V(\mu) = \mu$ and $$s(y;\mu) = \frac{y-\mu}{\mu}$$

And we could come up with quasi-likelihood models like $s(y;\mu) = \frac{y-\mu}{\phi\mu(1-\mu)}$ and $s(y;\mu) = \frac{y-\mu}{\phi\mu}$, then the parameter $\phi$ relates to under or over-dispersion and we are modelling as if the distribution is something that resembles Poisson or Binomial/Bernoulli but with a different (larger or smaller) variance.

---

does it make sense to test, e.g., $H_0: \phi=1$?

You can do that. The GLM model can be fitted independently from estimating the dispersion parameter. After estimating the GLM model, the dispersion parameter is fitted based on the sum of squared residuals divided by the function $V(\mu)$

$$\hat\phi = \frac{1}{n-m} \sum_{i=1}^n \frac{(x_i-\hat\mu_i)^2}{V(\hat\mu_i)}$$

See also: Is there a Relationship Between Variance and Chi-Square?

For large $n$ this is approximately normal distributed and the variance might be estimated with the kurtosis of the distribution of the residuals. (F-test for equality of variance for truncated normal distributions)

Answer 3

Ignore what summaries with $\theta = 1$ say and with a reasonable amount of data you can estimate $\theta$ from the data. The dispersion parameter $\theta$ in this parameterization is describing how the model deviates from a Poisson distribution (how much it is overparameterized). For $\theta \to \infty$ the distribution increasingly becomes like a Poisson distribution.

Doing a test vs. "$\theta = \infty$" is a bit tricky, because it's not a finite parameter value and it's also on the edge of the parameter space, which makes many things like confidence intervals etc. tricky (and that doesn't change if you use the paramterization used by SAS, which uses $\kappa := 1/\theta$ instead - although that makes it easier to deal with $\hat{\kappa}=0$). However, even with R's parameterization you could e.g. test for it being so high a value that it's almost a Poisson distribution anyway. I.e. null hypothesis it's a negative binomial with meaningful overdispersion vs. Poisson, say, $\theta < \theta_0$ for some relatively large $\theta_0$ and alternative that it's pretty close to a Poisson $\theta \geq \theta_0$.

There's of course other existing tests for overdispersion.