GLM. Binomial negative gives confusing output

Question

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?

Answer 1

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 interaction 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.