20
$\begingroup$

I'm trying to fit generalized linear models to some sets of count data that might or might not be overdispersed. The two canonical distributions that apply here are the Poisson and Negative Binomial (Negbin), with E.V. $\mu$ and variance

$Var_P = \mu $

$Var_{NB} = \mu + \frac{\mu^2}{\theta} $

which can be fitted in R using glm(..,family=poisson) and glm.nb(...), respectively. There is also the quasipoisson family, which in my understanding is an adjusted Poisson with the same E.V. and variance

$Var_{QP} = \phi\mu$,

i.e. falling somewhere in between Poisson and Negbin. The main problem with the quasipoisson family is that there is no corresponding likelihood for it, and hence a lot of extremely useful statistical tests and fit measures (AIC, LR etcetera) are unavailable.

If you compare the QP and Negbin variances, you might notice that you could equate them by putting $\phi = 1 + \frac{\mu}{\theta}$. Continuing on this logic, you could try to express the quasipoisson distribution as a special case of the Negbin:

$QP\,(\mu,\phi) = NB\,(\mu,\theta = \frac{\mu}{\phi-1})$,

i.e. a Negbin with $\theta$ linearly dependent on $\mu$. I tried to verify this idea by generating a random sequence of numbers according to the above formula and fitting it with glm:

#fix parameters

phi = 3
a = 1/50
b = 3
x = 1:100

#generating points according to an exp-linear curve
#this way the default log-link recovers the same parameters for comparison

mu = exp(a*x+b) 
y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator

#fit a generalized linear model y = f(x)  
glmQP = glm(y~x, family=quasipoisson) #quasipoisson
glmNB = glm.nb(y~x) #negative binomial

> glmQP

Call:  glm(formula = y ~ x, family = quasipoisson)

Coefficients:
(Intercept)            x  
    3.11257      0.01854  
(Dispersion parameter for quasipoisson family taken to be 3.613573)

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      2097 
Residual Deviance: 356.8    AIC: NA

> glmNB

Call:  glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)

Coefficients:
(Intercept)            x  
    3.10182      0.01873  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      578.1 
Residual Deviance: 107.8    AIC: 824.7

Both fits reproduce the parameters, and the quasipoisson gives a 'reasonable' estimate for $\phi$. We can now also define an AIC value for the quasipoisson:

df = 3 # three model parameters: a,b, and phi
phi.fit = 3.613573 #fitted phi value copied from summary(glmQP)
mu.fit = glmQP$fitted.values 

#dnbinom = negbin density, log=T returns log probabilities
AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T))
> AIC
[1] 819.329

(I had to manually copy the fitted $\phi$ value from summary(glmQP), as I couldn't find it in the glmQP object)

Since $AIC_{QP} < AIC_{NB}$, this would indicate that the quasipoisson is, unsurprisingly, the better fit; so at least $AIC_{QP}$ does what it should do, and hence it might be a reasonable definition for the AIC (and by extension, likelihood) of a quasipoisson. The big questions I am left with are then

  1. Does this idea make sense? Is my verification based on circular reasoning?
  2. The main question for anyone that 'invents' something that seems to be missing from a well established topic: if this idea makes sense, why isn't it already implemented in glm ?

Edit: figure added

glm fits and +-1 sigma bands

$\endgroup$
  • 1
    $\begingroup$ (+1) Welcome to Cross Validated! And thank you for an excellent question (though a few comments in the code might be nice for people who don't use R). I think you may have re-invented the NB1 model (though I haven't followed it in detail yet). Note also that there isn't a quasi-Poisson distribution - which is why there's no likelihood or AIC - it just refers to a way of fitting means & variances. $\endgroup$ – Scortchi - Reinstate Monica Jun 18 '15 at 21:03
  • 2
    $\begingroup$ Thanks! I've added some comments in the meantime, I hope that clears things up. I understand that the quasi-Poisson distribution does not exist per se - what I was really trying to figure out is why QP is even a thing at all, considering that the NB1 distribution exists and doesn't have any of the QP's quasi-problems (see Achims answer for an apparent resolution). $\endgroup$ – user28400 Jun 19 '15 at 9:51
  • 1
    $\begingroup$ @Scortchi --- actually, there is such a distribution... If $X\sim\text{Pois}(\lambda)$, and $Y=kX$, then $Y$ is exponential family with mean $\mu=k\lambda$ and variance $k\mu$. If $k\neq 1$. It's not necessarily suitable for count data (except as an approximation), since it's defined on $0,k,2k,...$. $\endgroup$ – Glen_b Jun 19 '15 at 10:00
  • 1
    $\begingroup$ @Glen_b: Do people really call that the quasi-Poisson? In any case it's a nice illustration - when you use a "quasiPoisson" model you're not really assuming that distribution, or the NB1, or any other, just a relationship between mean & variance that makes your estimates of coefficients & their standard errors better as the sample gets larger. $\endgroup$ – Scortchi - Reinstate Monica Jun 19 '15 at 14:58
  • 1
    $\begingroup$ @Scortchi It's the only exponential family distribution that satisfies the assumptions of the quasi-Poisson, so sort of -- on occasion I've seen people point out that it's the distribution that the assumption implies. Of course when people use it, they almost* never intend that their data are from that specific distribution -- it's just intended as a rough description of how their mean and variance relate. (It may make sense under very simple assumptions in some insurance applications - total claims cost, where the number of claims is Poisson and the cost per claim is effectively constant.) $\endgroup$ – Glen_b Jun 20 '15 at 2:53
23
$\begingroup$

The quasi-Poisson is not a full maximum likelihood (ML) model but a quasi-ML model. You just use the estimating function (or score function) from the Poisson model to estimate the coefficients, and then employ a certain variance function to obtain suitable standard errors (or rather a full covariance matrix) to perform inference. Hence, glm() does not supply and logLik() or AIC() here etc.

As you correctly point out, a model with the same expectation and variance function can be embedded into the negative binomial (NB) framework if the size parameter $\theta_i$ varies along with the expectation $\mu_i$. In the literature, this is typically called the NB1 parametrization. See for example the Cameron & Trivedi book (Regression Analysis of Count Data) or "Analysis of Microdata" by Winkelmann & Boes.

If there are no regressors (just an intercept) the NB1 parametrization and the NB2 parametrization employed by MASS's glm.nb() coincide. With regressors they differ. In the statistical literature the NB2 parametrization is more frequently used but some software packages also offer the NB1 version. For example in R, you can use the gamlss package to do gamlss(y ~ x, family = NBII). Note that somewhat confusingly gamlss uses NBI for the NB2 parametrization and NBII for NB1. (But jargon and terminology is not unified across all communities.)

Then you could ask, of course, why use quasi-Poisson if there is NB1 available? There is still a subtle difference: The former uses quasi-ML and obtains the estimate from the dispersion from the squared deviance (or Pearson) residuals. The latter uses full ML. In practice, the difference often isn't large but the motivations for using either model are slightly different.

$\endgroup$
  • 1
    $\begingroup$ Thanks! Very helpful answer, I'm experimenting with gamlss now and it looks like it's exactly what I needed. Could you elaborate on the motivations for using quasi-likelihood versus full ML? $\endgroup$ – user28400 Jun 19 '15 at 9:54
  • 2
    $\begingroup$ You assume less: You just assume (1) a log-linear relationship between the expectation and the regressors (2) a linear relationship between variance and expectation. The rest of the likelihood is left completely unspecified. As an alternative to (2), practitioners sometimes employ so-called "robust" sandwich standard errors which would allow for more general heteroskedasticity patterns. Of course, one could also employ the NB1 with sandwich standard errors... A few more comments are in our vignette("countreg", package = "pscl"). $\endgroup$ – Achim Zeileis Jun 19 '15 at 12:00

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.