30
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • $\begingroup$ @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. $\endgroup$ – lokheart Jan 25 '11 at 15:32
  • $\begingroup$ 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 $\endgroup$ – Mr. Nobody Sep 18 '15 at 13:36
  • 1
    $\begingroup$ epicalc was removed (source). An alternative could be lmtest. $\endgroup$ – Martin Thoma Dec 17 '15 at 13:51
26
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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 $\endgroup$ – lokheart Jan 25 '11 at 8:04
  • 3
    $\begingroup$ @lokheart anova(model1, model0) will work too. $\endgroup$ – chl Jan 25 '11 at 8:05
  • 5
    $\begingroup$ @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. $\endgroup$ – Gavin Simpson Jan 25 '11 at 8:13
  • $\begingroup$ Or you can do it "manually": p-value of the LR test = 1-pchisq(deviance, dof) $\endgroup$ – Umka Jan 24 '19 at 23:40
25
$\begingroup$

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)
Loading required package: lmtest
Loading required package: zoo
> ## 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
Improve this answer
$\endgroup$
4
  • $\begingroup$ +1 It's good to know (and it seems I forgot about that package). $\endgroup$ – chl Jan 25 '11 at 9:28
  • 3
    $\begingroup$ @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? $\endgroup$ – Kerry Jul 17 '14 at 13:27
  • 5
    $\begingroup$ @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). $\endgroup$ – Gavin Simpson Oct 24 '14 at 14:15
  • $\begingroup$ @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? $\endgroup$ – Branden Keck Jan 25 at 19:30

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.