# Likelihood ratio test in R

Suppose I am going to do a univariate logistic regression on several independent variables, like this:

mod.a <- glm(x ~ a, data=z, family=binominal("logistic"))
mod.b <- glm(x ~ b, data=z, family=binominal("logistic"))


I did a model comparison (likelihood ratio test) to see if the model is better than the null model by this command

1-pchisq(mod.a$null.deviance-mod.a$deviance, mod.a$df.null-mod.a$df.residual)


Then I built another model with all variables in it

mod.c <- glm(x ~ a+b, data=z, family=binomial("logistic"))


In order to see if the variable is statistically significant in the multivariate model, I used the lrtest command from epicalc

lrtest(mod.c,mod.a) ### see if variable b is statistically significant after adjustment of a
lrtest(mod.c,mod.b) ### see if variable a is statistically significant after adjustment of b


I wonder if the pchisq method and the lrtest method are equivalent for doing loglikelihood test? As I dunno how to use lrtest for univate logistic model.

• @Gavin thanks for reminding me, as comparing with stackoverflow, I need to spend more time to "digest" the answer before deciding whether the answer is appropriate or not, anyway, thanks again. Commented Jan 25, 2011 at 15:32
• I would not recommend using waldtest from lmtest. Use the aod package for model testing. Its far more straightforward. cran.r-project.org/web/packages/aod/aod.pdf Commented Sep 18, 2015 at 13:36
• epicalc was removed (source). An alternative could be lmtest. Commented Dec 17, 2015 at 13:51

Basically, yes, provided you use the correct difference in log-likelihood:

> library(epicalc)
> model0 <- glm(case ~ induced + spontaneous, family=binomial, data=infert)
> model1 <- glm(case ~ induced, family=binomial, data=infert)
> lrtest (model0, model1)
Likelihood ratio test for MLE method
Chi-squared 1 d.f. =  36.48675 , P value =  0
> model1$deviance-model0$deviance
[1] 36.48675


and not the deviance for the null model which is the same in both cases. The number of df is the number of parameters that differ between the two nested models, here df=1. BTW, you can look at the source code for lrtest() by just typing

> lrtest


at the R prompt.

• thanks, and I just found that i can use glm(output ~ NULL, data=z, family=binomial("logistic")) for creating a NULL model, and so i can use the lrtest afterwards. FYI, thanks again Commented Jan 25, 2011 at 8:04
• @lokheart anova(model1, model0) will work too.
– chl
Commented Jan 25, 2011 at 8:05
• @lokheart glm(output ~ 1, data=z, family=binomial("logistic")) would be a more natural null model, which says that output is explained by a constant term (the intercept)/ The intercept is implied in all your models, so you are testing for the effect of a after accounting for the intercept. Commented Jan 25, 2011 at 8:13
• Or you can do it "manually": p-value of the LR test = 1-pchisq(deviance, dof) Commented Jan 24, 2019 at 23:40

An alternative is the lmtest package, which has an lrtest() function which accepts a single model. Here is the example from ?lrtest in the lmtest package, which is for an LM but there are methods that work with GLMs:

> require(lmtest)
> ## with data from Greene (1993):
> ## load data and compute lags
> data("USDistLag")
> usdl <- na.contiguous(cbind(USDistLag, lag(USDistLag, k = -1)))
> colnames(usdl) <- c("con", "gnp", "con1", "gnp1")
> fm1 <- lm(con ~ gnp + gnp1, data = usdl)
> fm2 <- lm(con ~ gnp + con1 + gnp1, data = usdl)
> ## various equivalent specifications of the LR test
>
> ## Compare two nested models
> lrtest(fm2, fm1)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ gnp + gnp1
#Df  LogLik Df  Chisq Pr(>Chisq)
1   5 -56.069
2   4 -65.871 -1 19.605  9.524e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> ## with just one model provided, compare this model to a null one
> lrtest(fm2)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ 1
#Df   LogLik Df  Chisq Pr(>Chisq)
1   5  -56.069
2   2 -119.091 -3 126.04  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• +1 It's good to know (and it seems I forgot about that package).
– chl
Commented Jan 25, 2011 at 9:28
• @GavinSimpson This may seem silly, but how would you interpret the 'lrtest(fm2,fm1)' results? Model 2 is significantly different than model 1 and therefore the addition of the con1 variable was useful? Or the lrtest(fm2) is saying that model 2 is significantly different from model 1? But which model is better? Commented Jul 17, 2014 at 13:27
• @Kerry fm1 has a lower log likelihood and hence a poorer fit than fm2. The LRT is telling us that the degree to which we made fm1 a poorer model than fm2 is unexpectedly large if the terms that are different between the models were useful (explained the response). lrtest(fm2) isn't compared with fm1 at all, the model fm2 is compared with in that case if, as stated in the output, this: con ~ 1. That model, the null model, says that the best predictor of con is the sample mean of con (the intercept/constant term). Commented Oct 24, 2014 at 14:15
• @GavinSimpson , I know this is an old post, but I want to confirm my understanding: if Pr(>Chisq) were to be a large value (close to 1), is it then the case that we have rejected a null hypothesis stating that the models are significantly different? Then, we can claim that the "con1" term does not have a statistically significant impact on the model? Commented Jan 25, 2021 at 19:30
• @BrandenKeck If Pr(>Chisq) is large, we fail to reject the null hypothesis that the models are equally good. It's not evidence that the models are the same, but it's lack of evidence that they are different. For example, if you don't have a lot of data, you will fail to reject the NULL but you also should not be confident that the models are not different. Commented Jul 31, 2022 at 12:51