[R] glht function in multcomp gives unexpected p=1 for all comparisons
Katie Boothroyd-Roberts
katieboothroyd at yahoo.ca
Mon Nov 19 03:34:21 CET 2012
Hi, I have data with binomial response variable (survival) and 2 categorical independent variables (site and treatment) (see below). I have run a binomial GLM and found that both IVs and the interaction are significant. Now I want to do a post-hoc test for all pairwise comparisons to see which treatment groups differ. I tried the glht function in the multcomp package, but I get surprising results, with p=1 for all comparisons.
# My data:
surv.data <- data.frame(
Site=c(rep("Site1", 9), rep("Site2", 9)),
Treatment=rep(c(rep("Treat1", 3), rep("Treat2", 3), rep("Treat3", 3)), 2),
survival=c(0.9, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.8, 0.4, 0.6, 0.5, 1.0, 0.7, 0.2, 0.2, 0.4)
)
# Binomial GLM:
glm.out <- glm(survival~Site*Treatment, data=surv.data, family="binomial", weights=rep(10, nrow(surv.data)))
anova(glm.out, test="Chisq") # Treatment effect: p=0.001291
# Post-hoc test
library(multcomp)
summary(glht(glm.out, mcp(Treatment="Tukey", interaction_average=TRUE)))
This gives me the following results:
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glm(formula = survival ~ Site * Treatment, family = "binomial",
data = surv.small, weights = rep(10, nrow(surv.small)))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Deltoide - Baumier == 0 -9.179 2133.341 -0.004 1
Forest - Baumier == 0 -1.012 3016.999 0.000 1
Forest - Deltoide == 0 8.167 2133.341 0.004 1
(Adjusted p values reported -- single-step method)
Can anyone explain to me why p=1 for all comparisons even though the GLM showed that the treatment main effect was highly significant?
I searched in the archives and found two similar questions, but without helpful answers:
http://r.789695.n4.nabble.com/glht-problem-tt890671.html
http://r.789695.n4.nabble.com/glht-multicomparisons-with-a-binomial-response-variable-tt4360898.html
Thank-you in advance!
Katie
More information about the R-help
mailing list