9
$\begingroup$

I have these data:

set.seed(1)
predictor  <- rnorm(20)
set.seed(1)
counts <- c(sample(1:1000, 20))
df <- data.frame(counts, predictor)

I ran a poisson regression

poisson_counts <- glm(counts ~ predictor, data = df, family = "poisson")

And a negative binomial regression:

require(MASS)
nb_counts <- glm.nb(counts ~ predictor, data = df)

Then I calculated for dispersion statistics for the poisson regression:

sum(residuals(poisson_counts, type="pearson")^2)/df.residual(poisson_counts)

# [1] 145.4905

And the negative binomial regression:

sum(residuals(nb_counts, type="pearson")^2)/df.residual(nb_counts)

# [1] 0.7650289

Is anyone able to explain, WITHOUT USING EQUATIONS, why the dispersion statistic for the negative binomial regression is considerably smaller than the dispersion statistic for the poisson regression?

$\endgroup$

2 Answers 2

10
$\begingroup$

This is rather straightforward, but the "without using equations" is a substantial handicap. I can explain it in words, but those words will necessarily mirror equations. I hope that will be acceptable / still of some value to you. (The relevant equations are not difficult.)

There are several types of residuals. Raw residuals are simply the difference between the observed response values (in your case the counts) and the model's predicted response values. Pearson residuals divide those by the standard deviation (the square root of the variance function for the particular version of the generalized linear model that you are using).

The standard deviation associated with the Poisson distribution is smaller than that of the negative binomial. Thus, when you divide by a larger denominator, the quotient is smaller.

In addition, the negative binomial is more appropriate to your case, because your counts will be distributed as a uniform in the population. That is, their variance will not equal their mean.

$\endgroup$
4
  • 4
    $\begingroup$ Although the O.P. asks for a non-mathematical explanation, it would still be nice to see mathematical (or some equally rigorous and clear) justification for this answer. Upon reading the question my intuition was that "Because the Poisson is a (limiting) special case of the NB and the NB has more parameters, there is more flexibility in fitting, so of course any reasonable measure of residuals ought not to increase when replacing a Poisson GLM by a NB GLM." I am wondering whether such intuition was really correct. $\endgroup$
    – whuber
    May 20, 2014 at 16:01
  • $\begingroup$ If $X\sim\text{Poisson}(\lambda)$, $E[X]=V[X]=\lambda$. If $X\sim\text{NegBin}(r,p)$, $E[X]=pr/(1-p)$ and $V[X]=pr/(1-p)^2$. So a Poisson variance is equal to the mean, a NegBin variance is larger than the mean ($p<1\Rightarrow (1-p)^2<(1-p)$). This is why "the standard deviation associated with the Poisson distribution is smaller than that of the negative binomial." $\endgroup$
    – Sergio
    May 20, 2014 at 20:04
  • 3
    $\begingroup$ @Sergio The crux of the matter, though, is that in the Poisson model we are working with the estimate $\hat\lambda$ rather than $\lambda$ itself and in the NB model we are similarly working with two estimates $\hat{r}$ and $\hat{p}$. Your comparison therefore does not directly apply. Without actually writing down the formulas for the MLEs in both models, it is not at all obvious what the relationships must be between those sets of estimates. Furthermore, the Pearson residual is a ratio and the argument about variances addresses only the denominators, which is just half the story. $\endgroup$
    – whuber
    May 20, 2014 at 20:37
  • $\begingroup$ MLE estimates are consistent. The problem is that when, as gung says, "counts will be distributed as a uniform in the population. That is, their variance will not equal their mean", you'll never be able to get an estimated Poisson variance larger than an estimated Poisson mean, even if your estimates are unbiased and consistent. It's a problem of misspecification. $\endgroup$
    – Sergio
    May 21, 2014 at 5:08
6
$\begingroup$

For the Poisson model, if the expection for the $i$th observation $Y_i$ is $\mu_i$ its variance is $\mu_i$, & the Pearson residual therefore

$$\frac{y_i-\hat\mu_i}{\sqrt{\hat\mu_i}}$$

where $\hat\mu$ is the estimate of the mean. The parametrization of the negative binomial model used in MASS is explained here. If the expection for the $i$th observation $Y_i$ is $\mu_i$ its variance is $\mu_i + \frac{\mu^2}{\theta}$, & the Pearson residual therefore

$$\frac{y_i-\tilde\mu_i}{\sqrt{\tilde\mu_i+\frac{\tilde\mu'^2}{\theta}}}$$

where $\tilde\mu$ is the estimate of the mean. The smaller the value of $\theta$— i.e. the more extra-Poisson variance—, the smaller the residual compared to its Poisson equivalent. [But as @whuber has pointed out, the estimates of the means are not the same, $\hat\mu\neq\tilde\mu$, because the estimation procedure weights observations according to their assumed variance. If you were to make replicate measurements for the $i$th predictor pattern, they'd get closer, & in general adding a parameter should give a better fit across all observations, though I don't know how to demonstrate this rigorously. All the same, the population quantities you're estimating are larger if the Poisson model holds, so it shouldn't be a surprise.]

$\endgroup$
8
  • 1
    $\begingroup$ Thanks for introducing some of the equations. But are the $\mu_i$ in the two models going to have the same values? (I don't think so.) If not, how then is it possible to compare the two Pearson residuals? $\endgroup$
    – whuber
    May 20, 2014 at 20:02
  • $\begingroup$ @whuber In this case, it turns out that the fitted values for both models are nearly identical. After all, the "true" model really just has an intercept and is basically modeling the mean since there is no relationship between x and Y in the simulation. $\endgroup$
    – jsk
    May 20, 2014 at 20:22
  • 1
    $\begingroup$ @jsk Yes, I have looked at the data and ran the code. (BTW, it's possible to change the data and obtain essentially the same dispersion statistic for the two models.) Alas, your point, which is valid, still doesn't settle the specific question nor does it address the (implicit) general question about comparing Poisson residuals to NB residuals, because the estimated variances also could be nearly identical. One potentially confusing aspect about the present answer is the use of the symbol "$\mu_i$" to refer to what (in principle) could be different estimates in two models of the same data. $\endgroup$
    – whuber
    May 20, 2014 at 20:34
  • 1
    $\begingroup$ @whuber Indeed, you have valid points about the use of $\mu_i$. Interestingly, I can't seem to be able to find a way to simulate data that would result in a lower dispersion statistic for Poisson than NB. Perhaps it is not possible? I agree that this makes sense intuitively. Not easy to prove since there do not exist closed form solution for the mle when you have a glm with a link function other than the identity. But yes, it is easy to make the two dispersion statistics very similar. $\endgroup$
    – jsk
    May 20, 2014 at 21:01
  • 1
    $\begingroup$ @jsk - one theoretical argument to suspect that a NB model will always fit better than Poisson, is that you can write NB as a poisson-gamma compound distribution. So you have $( y_i|\lambda, v_i, r)\sim Poisson (\lambda v_i) $ and then $(v_i|\lambda, r)\sim Gamma(r, r) $ gives a negative binomial model $ ( y_i|\lambda, r)\sim NB (r,\frac {\lambda}{r+\lambda} )$. Now the addition of those $ v_i $ parameters allows the model to make the predicted mean closer to the observed value (when $ y_i>\lambda $ you would see $ v_i> 1$ , reducing the residual.) $\endgroup$ May 21, 2014 at 13:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.