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.

    I would not recommend using waldtest from lmtest. Use the aod package for model testing. Its far more straightforward.

    `epicalc` was removed (source). An alternative could be `lmtest`.

  • chl

    chl Correct answer

    10 years ago

    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

    @lokheart `anova(model1, model0)` will work too.

    @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.

    Or you can do it "manually": p-value of the LR test = 1-pchisq(deviance, dof)

License under CC-BY-SA with attribution

Content dated before 6/26/2020 9:53 AM