5
$\begingroup$

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:

  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.

$\endgroup$
2
$\begingroup$

anova{stats} is for Type I only, and has no way of doing Type III ANOVA. Anova{car} uses Type II or III tests.

You might find other helpful bits of information in these two threads:

Choice between Type-I, Type-II, or Type-III ANOVA

Difference between anova and Anova function

$\endgroup$
  • 1
    $\begingroup$ True, though here there's only one predictor variable. $\endgroup$ – Aaron - Reinstate Monica Jul 14 '16 at 18:53
  • $\begingroup$ 1/2: Thank you for your response. However, I don't know if that's at play here since the nested models I'm comparing only differ by one term (Species), so the type I SS for Species should be the same as Type III here: What is the effect of Species in a model that already includes all other terms? I also played with the models a bit (not shown in original post), and no matter whether I put Species or the random effect (Group) first in the model statement, I get the same results as shown in my post. I edited my post to show this. $\endgroup$ – Meg Jul 14 '16 at 19:34
  • $\begingroup$ 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. $\endgroup$ – Meg Jul 14 '16 at 19:34

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.