Likelihood ratio test vs p value for Poisson regression

Question

I have a Poisson regression model, from its summary table, I could see the p-value for a certain variable, e.g. gender. Since the p-value is testing the hypothesis whether the coefficient of gender equals to 0 (and it is comparing a model with all variables vs a model without this variable gender), seeing a p-value of 0.000 gives me the information that this variable is significant and should be included. Now consider a likelihood ratio test where the null hypothesis is the reduced model without gender, and the alternative hypothesis is the full model with gender. From likelihood ratio test, if the p-value is small than a significance value, then we reject the null (reduced model), and therefore gender as a variable should be included.

I am a little bit confused about whether these two tests are equivalent in this scenario. Are they testing the same thing and will always produce same decision in this specific case? And if not, where are their differences?

Answer 1

Assuming that you're using glm(), which is the standard/most-used tool for fitting Poisson regressions in R: the $p$-values in the table printed by summary() are so-called Wald p-values, which make a specific assumption about the sampling distribution of the parameter estimates (specifically that they are [multivariate] Normal, or equivalently that the log-likelihood surface is quadratic). The $p$-values that you get from a likelihood ratio test (e.g., with anova(full_model, restricted_model) where restricted_model omits the focal parameter, or drop1(full_model)) also makes an assumption (that 2 x the log-likelihood difference between the models follows a chi-squared distribution when the null hypothesis is true), but the LRT assumption is much weaker than the one assumed by the Wald p-values, so the LRT $p$-values are in general more accurate.

As a side note, your statement that "... this variable is significant and should be included ..." makes me a bit nervous, as it suggests that you are using hypothesis testing to decide on the structure of your final model, which may be problematic (see this answer) ...

As an example from UCLA OARC:

p <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
p <- within(p, {
  prog <- factor(prog, levels=1:3, labels=c("General", "Academic", 
                                                     "Vocational"))
  id <- factor(id)
})
m1 <- glm(num_awards ~ prog + math, family="poisson", data=p)

Wald (summary())

              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -5.24712    0.65845  -7.969 1.60e-15 ***
progAcademic    1.08386    0.35825   3.025  0.00248 ** 
progVocational  0.36981    0.44107   0.838  0.40179    
math            0.07015    0.01060   6.619 3.63e-11 ***

We get the same p-value for math (for example) using aod::wald.test:

wald.test(b=coef(m1), Sigma=vcov(m1), Term=4)
Wald test:
----------

Chi-squared test:
X2 = 43.8, df = 1, P(> X2) = 3.6e-11

(The X2 value is the square of the z-value listed in the glm summary)

LRT (drop1)

drop1(m1, test= "Chisq")
Single term deletions

Model:
num_awards ~ prog + math
       Df Deviance    AIC    LRT  Pr(>Chi)    
<none>      189.45 373.50                     
prog    2   204.02 384.08 14.572 0.0006852 ***
math    1   234.46 416.51 45.010  1.96e-11 ***

You'll note a few differences here ...

• for a categorical (factor) variable, summary() gives you a test for every level of the factor (except the baseline/first), representing differences between the baseline and other levels; drop1() gives you a single p-value for the combined effects (you can use car::Anova() to get a single Wald-based p-value per factor variable)
• for a numeric predictor, or a factor with only two levels (such as gender in most but not all cases), you'll get a single p-value either way, but they won't be identical. Both methods strongly reject the null hypothesis that there is no effect of math, but the p-values aren't the same (3.6e-11 for Wald and 1.96e-11 for LRT)