[R] Interaction contrasts or posthoc test for glm (MASS) with ANOVA design
David Winsemius
dwinsemius at comcast.net
Mon Nov 2 16:39:43 CET 2009
On Nov 2, 2009, at 8:40 AM, Maya Pfaff wrote:
> 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.
I'm a bit puzzled by this request, since it appears that you already
have the test result you seek in the form of the line beginning
Site1:Taxon1. (It might be a different story if you had more species.)
In the default glm setup, the contrasts are of type "treatment". Your
Site1:Taxon1 coefficients would then be the mean difference (and se)
from the estimates based on Intercept+Taxon+Site( on the scale being
regressed on). You will notice that the z value does not change when
you use treatment contrasts instead of those you specified.
The only thing further to do would be to construct a reduced model
without the interaction, take the difference in deviances of the
model, and compare to chi-sq(significance level, df=1) and that would
give you the ML test rather than the score test which results from the
by coefficient tests.
When I do that I get
s1: 2 x log-likelihood: -1155.4010
s2: 2 x log-likelihood: -1181.8600
so
> 1-pchisq(26.459,df=1)
[1] 2.691913e-07
(Which is immaterially different than the Pr(>|z|) that you get from
the default output.
Of course that was the way we did it in school 20 years ago and it
would be much faster to do:
> anova(s1,s2)
Likelihood ratio tests of Negative Binomial Models
Response: counts
Model theta Resid. df 2 x log-lik. Test df LR
stat. Pr(Chi)
1 Site + Taxon 0.5263897 151 -1181.860
2 Site * Taxon 0.6126726 150 -1155.401 1 vs 2 1
26.45960 2.691083e-07
--
David
> 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]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
More information about the R-help
mailing list