# R: anova() vs. Anova() for test of categorical predictor from glmer or glm.nb object

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,

####################
# 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,

####################
# 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:

1. Isn't restricted to only mixed or only fixed effects models
2. 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.

anova{stats} is for Type I only, and has no way of doing Type III ANOVA. Anova{car} uses Type II or III tests.
• 2/2: I also found the test statistics are different even when fitting a fixed effects only negative binomial count model with glm.nb() from the MASS package, where the only thing in the model is Species and the offset() option. This is inconsistent with the two functions providing the same test statistic when fitting the Poisson count model using glm(), as was shown in my original post. I also updated my post to show this case. – Meg Jul 14 '16 at 19:34