11

I recently fit 4 multiple regression models for the same predictor/response data. Two of the models I fit with Poisson regression.

model.pois <- glm(Response ~ P1 + P2 +...+ P5, family=poisson(), ...)
model.pois.inter <- glm(Response ~ (P1 + P2 +...+ P5)^2, family=poisson(), ...)

Two of the models I fit with negative binomial regression.

library(MASS)
model.nb <- glm.nb(Response ~ P1 + P2 +...+ P5, ...)
model.nb.inter <- glm.nb(Response ~ (P1 + P2 +...+ P5)^2, ...)

Is there a statistical test I can use to compare these models? I've been using the AIC as a measure of the fit, but AFAIK this doesn't represent an actual test.

MERose
  • 403
  • 1
  • 6
  • 20
Daniel Standage
  • 1,109
  • 3
  • 13
  • 21
  • You want to compare the models' *fit* using a statistical test, right? What kind of hypothesis would you like to test? – Firefeather Dec 16 '10 at 20:56
  • @Firefeather For example, I would like to test whether the fit of `model.nb.inter` is _significantly_ better than that of `model.pois.inter`. Yes, the AIC is lower, but how much lower constitutes _significantly better_? – Daniel Standage Dec 16 '10 at 21:00
  • Note: the answer to this question need not actually include the AIC. – Daniel Standage Dec 16 '10 at 21:07
  • I don't know the answer to this question, but I can give a start. I know you can use an $F$ test to compare `model.pois` against `model.pois.inter` (and similarly compare `model.nb` against `model.nb.inter`), but I can't guarantee that comparisons between a Poisson model and a negative binomial model would work. I wonder if an $F$ test to compare the variances of each pair would be reliable. – Firefeather Dec 16 '10 at 21:16
  • @Daniel, note that if you're doing more than one comparison between pairs of models, then you'll need to adjust for multiple comparisons (e.g., use [Scheffé's method](http://en.wikipedia.org/wiki/Scheffé's_method)). – Firefeather Dec 16 '10 at 21:19
  • 1
    @Firefeather, yes I'm aware of the need to control the familywise confidence level. Would Scheffe be more appropriate here than, say, Bonferroni? – Daniel Standage Dec 16 '10 at 21:40
  • @Daniel, I'd have to do some digging to remind myself whether Scheffé can be applied to Poisson or negative binomial models. However, I do know that other methods are typically preferred over the Bonferroni correction, because it is commonly regarded as too conservative (you give up more power than you should have to). – Firefeather Dec 16 '10 at 21:48

3 Answers3

16

You can compare the negative binomial model to the corresponding Poisson model with a likelihood ratio test. A Poisson model is equivalent to a negative binomial model with an overdispersion parameter of zero. Therefore they are nested models and likelihood ratios are valid. The complication is that the overdispersion parameter is restricted to be non-negative, i.e. it logically can't be less than zero, so the null hypothesis is on the boundary of the parameter space. This means that instead of comparing twice the log-likelihood to a chi-squared distribution with one degree of freedom, you need to compare it to a mixture distribution consisting of equal parts of a chi-squared with 1 d.f. and a point mass at zero (a chi-squared distribution with zero degrees of freedom). What that means in practice is that you can calculate the p-value using the chi-squared with 1 d.f. and then halve it. For more details and background, see Case 5 of Self & Liang JASA 1987; 82:605-610..

Note that some statistical software packages, such as Stata, will do this all for you automatically when you fit a negative binomial model. In fact I've shamelessly cribbed much of the above from the Stata help system -- if you have Stata see help j_chibar.

Scortchi - Reinstate Monica
  • 27,560
  • 8
  • 81
  • 248
onestop
  • 16,816
  • 2
  • 53
  • 83
5

I believe anova() in R can be used for this. Despite its name, it's a likelihood ratio test. Crawley in his The R Book has some examples of usage.

Roman Luštrik
  • 3,338
  • 3
  • 31
  • 39
1

As onestop notes, because the models are nested you can perform a likelihood ratio test.

In general though that isn't true, so if you want to compare non-nested models you could use Vuong's test.

Xodarap
  • 2,448
  • 2
  • 17
  • 24