# Interpreting output of analysis of deviance table from anova() model comparison

0 down vote favorite

I have a large multivariate abundance data and I am interested in comparing multiple models that fit different combinations of three categorical predictor variables to my species matrix response variable. I have been using anova() to compare my different models, but I am having difficulty interpreting the output. Below, I have given my code as well as the corresponding R output.

invert.mvabund <- mvabund(mva.dat)
null<-manyglm(mva.dat~1, family='negative.binomial')
m1 <- manyglm(mva.dat~Habitat+Detritus, family='negative.binomial')
m2 <- manyglm(mva.dat~Habitat*Detritus, family='negative.binomial')
m3 <- manyglm(mva.dat~Habitat*Detritus+Block, family='negative.binomial')
anova(null,m1,m2,m3)

Analysis of Deviance Table

null: mva.dat ~ 1
m1: mva.dat ~ Habitat + Detritus
m2: mva.dat ~ Habitat * Detritus
m3: mva.dat ~ Habitat * Detritus + Block

Multivariate test:
Res.Df Df.diff   Dev Pr(>Dev)
null     99
m1       94       5 257.2    0.001 ***
m2       90       4  87.7    0.003 **
m3       81       9 173.5    0.003 **
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


How do I interpret these results? Is m2 the best-fitting model because it has the lowest deviance, even though it has a higher p-value than m1? Is this because the p-value is suggesting that there is a significant level of deviance, so the optimal model will have a higher p-value? Any suggestions on how to interpret these results would be much appreciated- I haven't been able to find a clear answer in my Google searches. Thanks!

You misunderstood some values. Multivariate test table's Dev is decrement from upper model (When a model has a interaction term, become a litte more complex). So I think each p-value indicates the difference between the model and upper one is statistically significant or not.

Here is my example code.

library(mvabund)
iris2 <- cbind(iris[,1:4]*10, Species=iris[,5]) # make integer df

null <- manyglm(Petal.Length ~ 1, data=iris2, family="poisson")
m1 <- manyglm(Petal.Length ~ Sepal.Length + Species, data=iris2, family="poisson")
m2 <- manyglm(Petal.Length ~ Sepal.Length * Species, data=iris2, family="poisson")
m3 <- manyglm(Petal.Length ~ Sepal.Length * Species + Petal.Width, data=iris2, family="poisson")
anova(null, m1, m2, m3)

# Analysis of Deviance Table

# null: Petal.Length ~ 1
# m1: Petal.Length ~ Sepal.Length + Species
# m2: Petal.Length ~ Sepal.Length * Species
# m3: Petal.Length ~ Sepal.Length * Species + Petal.Width

# Multivariate test: (I used anova(null, m1, m2, m3)$table to increase digit number) # Res.Df Df.diff Dev Pr(>Dev) # null 149 NA NA NA # m1 146 3 1355.3139486 0.001 # m2 144 2 3.0393296 0.001 # m3 143 1 0.1406748 0.361 anova(null, m1)$table[2,3]                           # 1355.314
anova(m1, m3)$table[2,3] - anova(m2, m3)$table[2,3]  # 3.03933
anova(m2, m3)$table[2,3] # 0.1406748  • Okay, so each model is being compared to the model above it, give or take the interactions, and not to the null model. So if I'm interested in understanding which of these models best represents my species matrix, this analysis does not provide this information. Any suggestions of alternate methods to determine which model best fits my data? – user2096647 Jun 10 '16 at 14:08 • Or do I just need to compare each model individually to one another one at a time? – user2096647 Jun 10 '16 at 14:22 • In the manual p.44 of package ‘mvabund’ or ?manyglm, author(s) used drop1() function, i.e., they did AICs based model comparison. It indicates the lowest AIC model is the best. – cuttlefish44 Jun 10 '16 at 23:05 • Or if you want to test, I think you can do multiple comparison test using all pair's p-values and adjusting them. For example, p <- c(anova(null, m1)$table[2,4], anova(null, m2)$table[2,4], anova(null, m3)$table[2,4], anova(m1, m2)$table[2,4], anova(m1, m3)$table[2,4], anova(m2, m3)\$table[2,4]) mat <- matrix(p.adjust(p, "holm")[c(1, 2, 3, NA, 4, 5, NA, NA, 6)], byrow=T, ncol=3) rownames(mat) <- c("null", "m1", "m2") colnames(mat) <- c("m1", "m2", "m3") mat – cuttlefish44 Jun 10 '16 at 23:08