Dear R experts

I am running a negative-binomial GLM (glm.nb) to test the null hypotheses
that species 1 and 2 are equally abundant between site 1 and site2, and
between each other. So, I have a 2x2 factorial design with factors Site
(1,2) and Taxon (1,2).
Since the Site:Taxon interaction is significant, I need to do the equivalent
to a "post-hoc test" for ANOVA, however, the same tests (e.g. Tukey HSD) do
not seem to be applicable for the GLM. I tried specifying orthogonal
contrasts, but could not figure out what the interaction contrast (see
Site1:Taxon1 in below example) means.

Could you please advise me how to specify a meaningful interaction contrast
(i.e. contrast species within sites)? Alternatively, could you recommend a
way to do posthoc comparisons?

Thanks for your time and kind regards
Maya

> library(MASS)
> counts <- c(1, 4, 9, 2, 1, 4, 2, 4, 1, 3, 2, 2, 1, 3, 1, 1, 2, 1, 113, 83,
49, 46, 13, 52, 4, 10, 14, 10, 3, 19, 8, 21, 151, 186, 99, 11, 29, 24, 24,
62, 15, 98, 30, 21, 63, 29, 48, 11, 16, 35, 21, 17, 6, 2, 2, 3, 3, 4, 4, 2,
1, 2, 2, 3, 4, 8, 10, 3, 14, 3, 11, 23, 3, 51, 8, 8, 7, 1, 13, 8, 1, 0, 1,
0, 1, 1, 1, 0, 1, 1, 1, 5, 8, 1, 1, 20, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,
0, 3, 0, 1, 0, 5, 1, 1, 9, 0, 34, 4, 1, 17, 0, 7, 33, 86, 73, 67, 79, 109,
27, 37, 23, 12, 17, 41, 8, 38, 4, 23, 14, 49, 64, 39, 31, 156, 110, 97, 33,
170, 137, 72, 28, 54)
> Site <- factor(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2))
> Taxon <- factor(c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1))

> contrasts(Site) <- cbind(c(1,-1))
> contrasts(Taxon) <- cbind(c(1,-1))
> s1 <- glm.nb(counts ~ Site*Taxon)
> summary(s1)

Call:
glm.nb(formula = counts ~ Site * Taxon, init.theta = 0.612672617555492,
    link = log)

Deviance Residuals:
      Min         1Q     Median         3Q        Max
-2.269021  -1.215727  -0.454719   0.003515   2.517288

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    2.6139     0.1097  23.831  < 2e-16 ***
Site1         -0.8664     0.1097  -7.899 2.81e-15 ***
Taxon1        -0.3814     0.1097  -3.477 0.000506 ***
Site1:Taxon1  -0.5977     0.1097  -5.449 5.06e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.6127) family taken to be 1)

    Null deviance: 242.94  on 153  degrees of freedom
Residual deviance: 176.20  on 150  degrees of freedom
AIC: 1165.4

Number of Fisher Scoring iterations: 1


              Theta:  0.6127
          Std. Err.:  0.0690

 2 x log-likelihood:  -1155.4010
>
> TukeyHSD(s1)
Error in UseMethod("TukeyHSD") : no applicable method for "TukeyHSD"

	[[alternative HTML version deleted]]


