In R, I'm wondering how the functions anova() (stats package) and Anova() (car package) differ when being used to compare nested models fit using the glmer() (generalized linear mixed effects model; lme4 package) and glm.nb (negative binomial; MASS package) functions.
I've found the two ANOVA functions do not produce the same results for tests of fixed effects in a Poisson mixed model, or a negative binomial fixed effects model (no random effects). Results from both are shown below.
My goal: Correctly test the overall significance of a multi-level categorical predictor (fixed; Species). I'm looking for a type III SS-type p-value.
First: If one fits a fixed effects generalized linear model (Poisson here) using glm(), then these two functions do produce the same results given the arguments as in the following dummy example:
mod01 <- glm(Count ~ Species + offset(log(Area)), data=data01, family=poisson)
####################
# Anova() function #
####################
library(car)
Anova(mod01, type=3)
# Analysis of Deviance Table (Type III Wald chisquare tests)
# Response: Count
# LR Chisq Df Pr(>Chisq)
# Species 255.44 8 < 2.2e-16 ***
####################
# anova() function #
####################
mod01x <- update(mod01, . ~ . - Species)
anova(mod01x, mod01, test="Chisq")
# Model 1: Count ~ offset(log(Area))
# Model 2: Count ~ Species + offset(log(Area))
# Resid. Df Resid. Dev Df Deviance Pr(>Chi)
# 1 1063 1456.4
# 2 1055 1201.0 8 255.44 < 2.2e-16 ***
# Test statistics are the SAME (255.44) for the fixed effects model
However: For a generalized linear mixed effects model (using glmer() with random effect for Group), analogous code gives a different test statistic across the two functions:
library(lme4)
mod02 <- glmer(Count ~ 1 + Species + (1 | Group) + offset(log(Area)), data=data01,
family=poisson(link="log"), nAGQ=100)
####################
# Anova() function #
####################
Anova(mod02, type=3)
# Analysis of Deviance Table (Type III Wald chisquare tests)
# Response: Count
# Chisq Df Pr(>Chisq)
# (Intercept) 4.0029 1 0.04542 *
# Species 197.9012 8 < 2e-16 ***
####################
# anova() function #
####################
mod02x <- update(mod02, . ~ . - Species)
anova(mod02x, mod02, test="Chisq")
# mod02x: Count ~ (1 | Group) + offset(log(Area))
# mod02: Count ~ 1 + Species + (1 | Group) + offset(log(Area))
# Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
# mod02x 2 1423.9 1433.8 -709.95 1419.9
# mod02 10 1191.7 1241.4 -585.85 1171.7 248.21 8 < 2.2e-16 ***
# Now the test statistics are DIFFERENT (197.9012 vs. 248.21)
#####################################################################
# Not a matter of type I vs. III SS since whether the fixed or random
# effect is fit first in the model does not affect results:
# List random effect (Group) before fixed (Species):
mod03 <- glmer(Count ~ 1 + (1 | Group) + Species + offset(log(Area)), data=data01,
family=poisson(link="log"), nAGQ=100)
####################
# Anova() function #
####################
Anova(mod03, type=3)
# Analysis of Deviance Table (Type III Wald chisquare tests)
# Response: Count
# Chisq Df Pr(>Chisq)
# (Intercept) 4.0029 1 0.04542 *
# Species 197.9012 8 < 2e-16 ***
####################
# anova() function #
####################
mod03x <- update(mod03, . ~ . - Species)
anova(mod03x, mod03, test="Chisq")
# mod03x: Count ~ (1 | Group) + offset(log(Area))
# mod03: Count ~ 1 + (1 | Group) + Species + offset(log(Area))
# Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
# mod03x 2 1423.9 1433.8 -709.95 1419.9
# mod03 10 1191.7 1241.4 -585.85 1171.7 248.21 8 < 2.2e-16 ***
# Respective test statistics are the same as above case where order of fixed
# and random effects was reversed
Another example of inconsistent test statistics: Fixed effects negative binomial model:
library(MASS)
mod04 <- glm.nb(Count ~ Species + offset(log(Area)), data=data01)
####################
# Anova() function #
####################
Anova(mod04, type=3)
# Analysis of Deviance Table (Type III tests)
# Response: Spiders_Tree
# LR Chisq Df Pr(>Chisq)
# Species 101.08 8 < 2.2e-16 ***
####################
# anova() function #
####################
mod04x <- update(mod04, . ~ . - Species)
anova(mod04x, mod04)
# Likelihood ratio tests of Negative Binomial Models
# Response: Count
# Model theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
# 1 offset(log(Area_M2)) 0.2164382 1063 -1500.688
# 2 Species + offset(log(Area_M2)) 0.3488095 1055 -1413.651 1 vs 2 8 87.03677 1.887379e-15
# Test statistics are also DIFFERENT here (101.08 vs. 87.03677)
In summary: The problem:
- Isn't restricted to only mixed or only fixed effects models
- Isn't a matter of type I or III SS, since an example with only one predictor (negative binomial fixed effects model) showed the same problem, and even in the case of more than one predictor (mixed model example), the test is only for the removal of one predictor (Species), so I believe the two types of SS should be equivalent in this case.
Could it have to do with the offset? Maybe the functions were written to "behave well" with the glm() function, but process others (such as glmer() and glm.nb()) inconsistently? Something else I'm not thinking of?
I'm not providing data for my example code above, as I'm assuming someone can comment on the differing theories of each function without a minimal working example. However, if you would like to verify the results really do differ (as shown above), I will add a dummy dataset.
