GLM. Binomial negative gives confusing output

Given the following data:

Treatment <- as.factor(c(rep("Ploughed",8),rep("Mowed",8),
rep("Ploughed",9),rep("Mowed",10)))
Strata <- as.factor(c(rep("Canopy",16), rep("Ground",19)))
Herbivores <- as.integer(c(4,3,39,15,22,10,61,11,3,3,1,6,7,3,3,4,10,37,44,6,
77,222,49,101,74,112,23,47,86,42,104,203,111,310,41))
data <- data.frame(Treatment,Strata,Herbivores)


I am trying to fit a Negative Binomial generalized linear model, to see if Treatment affects the number of Herbivores, while controlling for strata. The data is overdispersed because we sampled two very different strata (canopy and ground cover).

First, I use this, to get the "theta" parameter:

library(MASS)
m1 <- glm.nb(formula=Herbivores~Treatment+Strata,link=log, data=data)
summary(m1)


Then, I use this to fit the final model:

m2 <- glm(formula=Herbivores~Treatment+Strata,
family=negative.binomial (theta=m1$theta, link=log),data=data) summary(m2)  Questions: 1. Why am I getting different outputs when using these models (m1 gives a$z$-value, m2 gives a$t$-value)? As far as I know they should do the same. 2. Why don't the results show the differences in abundance when there is quite a big difference? Am I using a wrong model? I see very different abundances using this code: aggregate(data$Herbivores, by=list(Category=data$Treatment), FUN=sum)  3. Why is the$t$-value higher (positive) for Ploughed treatment than for Mowed treatment (negative$t$-value) when ploughed treatment has less abundance of Herbivores? • I don't see why this should be closed. Questions 2 and 3 are statistical, not R-related. Close-voters, please reconsider. – Stephan Kolassa Oct 17 '17 at 7:52 • ... plus I'm working on an answer... – Stephan Kolassa Oct 17 '17 at 7:54 • Thank you Stephan, very clarifying. Regarding question 3 and your m3 model. Why is the z-value 3.993 times greater for "TreatmentPloughed" than "TreatmentMowed"? The TreatmentMowed level should have a "steeper" slope and also a higher abundance of Hervivores. Thank you again, Charly – Charly Oct 17 '17 at 10:10 • I edited my answer. You cannot interpret the coefficients on a main effect without taking the interaction into account. The effect of going from Ploughed to Mowed is$+1.70$in the Canopy, but$+1.70-2.15=-0.45$on the Ground. Plug a few values into the model to see the fits. (Incidentally, when you need clarification on a specific CV answer, it's better to comment at that specific answer, rather than at the question.) – Stephan Kolassa Oct 17 '17 at 13:11 1 Answer 1. The two commands fit similar models, but report different statistics. After all, they come from different R packages by different authors. No surprise that different authors choose to report different statistics. 2. (& 3.) Let's look at a plot. Always look at plots. plot(c(.5,2.5),range(data$Herbivores),type="n",xlab="",xaxt="n",ylab="")
axis(1,1:2,levels(data$Treatment)) with(data,points(jitter(as.numeric(Treatment)),Herbivores,pch=19,col=Strata)) legend("topright",pch=19,col=1:2,legend=levels(data$Strata))


We see that the abundance per stratum differs at the Mowed treatment - but not at the Ploughed treatment. This strongly suggests an between stratum and treatment. Let's fit this:

m3 <- glm.nb(formula=Herbivores~Treatment*Strata,link=log,data=data)
summary(m3)


Result:

                               Estimate Std. Error z value Pr(>|z|)
(Intercept)                      1.3218     0.3237   4.084 4.43e-05 ***
TreatmentPloughed                1.7047     0.4269   3.993 6.52e-05 ***
StrataGround                     3.3594     0.4035   8.325  < 2e-16 ***
TreatmentPloughed:StrataGround  -2.1535     0.5527  -3.897 9.76e-05 ***


That is, the model is

$$\text{Herbivores} = 1.32 + 1.70\times\text{Ploughed} + 3.36\times\text{Ground} -2.15\times\text{Ploughed}\times\text{Ground}.$$

And we indeed see a highly significant interaction. Note that Akaike's Information Criterion also very much prefers the model including an interaction:

> AIC(m1)
[1] 327.7352
> AIC(m3)
[1] 317.0077


This should also answer your 3rd question. Your model was misspecified.

Incidentally, I'm not sure why you'd get one NB parameter from one model, then fit a new model constraining that parameter. Better to do all in one approach with MASS:glm.nb() as we did for m3 here.