I was looking for ways to do a likelihood ratio test in R to compare model fits. I first coded it myself, then found both the default
anova() function and also
lrtest() in the
lmtest package. When I checked, though,
anova() always produces a slightly different p-value from the other two even though the 'test' parameter is set to "LRT". Is
anova() actually performing some subtly different test, or am I not understanding something?
Platform: R 3.2.0 running on Linux Mint 17,
lmtest version 0.9-33
set.seed(1) # Reproducibility n=1000 y = runif(n, min=-1, max=1) a = factor(sample(1:5, size=n, replace=T)) b = runif(n) # Make y dependent on the other two variables y = y + b * 0.1 + ifelse(a==1, 0.25, 0) mydata = data.frame(y,a,b) # Models base = lm(y ~ a, data=mydata) full = lm(y ~ a + b, data=mydata) # Anova anova(base, full, test="LRT") # lrtest library(lmtest) lrtest(base, full) # Homebrew log-likelihood test like.diff = logLik(full) - logLik(base) df.diff = base$df.residual - full$df.residual pchisq(as.numeric(like.diff) * 2, df=df.diff, lower.tail=F)
When I run it,
anova() gives a p-value of 0.6071, whereas the other two give 0.60599. A small difference, but consistent, and too big to be imprecision in how floating point numbers are stored. Can someone explain why
anova() gives a different answer?