Compare poisson and negative binomial regression with LR test

Question

My question is related to the question Compare negative binomial models. I have some difficulties understanding the UCLA guide (I am using R, too).

In the section Checking model assumption they run

pchisq(2 * (logLik(m1) - logLik(m3)), df = 1, lower.tail = FALSE)

I understand the first part, but why do they set the degrees to freedom equal to 1 (df = 1)?

Their output then is

## 'log Lik.' 2.157e-203 (df=5)

and they write

In this example the associated chi-squared value is 926.03 with one degree of freedom. This strongly suggests the negative binomial model, estimating the dispersion parameter, is more appropriate than the Poisson model.

What does "926.03 with one degree of freedom" refer to?

Answer 1

The Poisson and negative binomial (NB) model are nested: Poisson is a special case with theta = infinity. So a likelihood ratio test comparing the two models is testing the null hypothesis that "theta = infinity" against the alternative that "theta < infinity". Here the two models have the following log-likelihoods

R> logLik(m3)
'log Lik.' -1328.642 (df=4)
R> logLik(m1)
'log Lik.' -865.6289 (df=5)

Thus, the fitted log-likelihood in the NB model is much larger/better using just one additional parameter (4 regression coefficients plus 1 theta). The value of m1$theta is 1.032713, i.e., essentially corresponding to a geometric distribution (theta = 1). So clearly there is overdispersion compared to the Poisson distribution with a theta much lower than infinity. (My personal rule of thumb is that everything beyond 10 is already rather close to infinity.)

The likelihood ratio test statistic is then twice the absolute difference of log-likelihoods (2 * (1328.642 - 865.6289) = 926.0262) and has to be compared with a chi-squared distribution with degrees of freedom (df) equal to the difference in df's between the two models (5 - 4 = 1). This is what the code above does. However, it may be somewhat confusing that the difference of logLik() objects retains the "logLik" class and hence displays it with a label and df. Applying as.numeric() to drop the class might clarify this:

R> 2 * (logLik(m1) - logLik(m3))
'log Lik.' 926.0272 (df=5)
R> as.numeric(2 * (logLik(m1) - logLik(m3)))
[1] 926.0272

And then you could compute the p-value by hand:

R> stat <- as.numeric(2 * (logLik(m1) - logLik(m3)))
R> pchisq(stat, df = 5 - 4, lower.tail = FALSE)
[1] 2.157298e-203

There are also functions that carry out such likelihood ratio tests in a generic way, e.g., lrtest() in lmtest (among others):

R> library("lmtest")
R> lrtest(m3, m1)
Likelihood ratio test

Model 1: daysabs ~ math + prog
Model 2: daysabs ~ math + prog
  #Df   LogLik Df  Chisq Pr(>Chisq)    
1   4 -1328.64                         
2   5  -865.63  1 926.03  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As a final detail: The likelihood ratio test of Poisson vs. NB is non-standard because the parameter theta = infinity is on the boundary of the parameter space. Therefore, the p-value from the chi-squared distribution has to be halved (similar to a one-sided test). Of course, if you already have a significant result without halving, it's still significant afterwards...