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')

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.

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
  • $\begingroup$ 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? $\endgroup$ – user2096647 Jun 10 '16 at 14:08
  • $\begingroup$ Or do I just need to compare each model individually to one another one at a time? $\endgroup$ – user2096647 Jun 10 '16 at 14:22
  • 1
    $\begingroup$ 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. $\endgroup$ – cuttlefish44 Jun 10 '16 at 23:05
  • $\begingroup$ 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 $\endgroup$ – cuttlefish44 Jun 10 '16 at 23:08

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.