[R] Help with glht function and binary data!
Jessica McLachlan
jessica.mclachlan at bristol.ac.uk
Tue Nov 20 22:02:10 CET 2012
Hi,
I carried out an experiment, using a repeated measures design, in which I
broadcast 6 playbacks to 15 nests. I am trying to run a posthoc pairwise
comparison on a binomial GLMM with missing data points to determine which
playback treatments differ.
The response variable (Response) is binary, there is one categorical factor
(Treatment) and Nest is included as a random factor. I have run a binomial
GLMM using the lme4 package and found the effect of Treatment to be highly
significant. I followed this up using the glht function in the multcomp
package to look at the pairwise comparisons.
However, it is giving me strange results for some of the comparisons,
producing p-values of 1 for comparisons between what graphically (plot of
mean and std error/boxplots) appear to be the most different treatments.
For my other binomial GLMM, where I look at proportion rather than binary
data, the glht function yields the expected results without problems.
Can anyone help? Is the glht function not appropriate for binary data? Or
have I got something wrong?
Thank you!
Jessica
My data and code are below.
# My data:
Nest<-c(1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8,9,9,9,9,9,9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12,13,13,13,13,13,14,14,14,14,15,15,15,15,15)
Treatment<-factor(c("CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CRS","CL","CS","SL","SS","CRL","CL","CS","SL","SS","CRL","CS","SL","SS","CRL","CRS","CS","SL","SS"))
Response<-c(0,0,1,1,1,1,0,0,1,1,1,1,0,0,1,1,1,1,0,1,1,1,1,1,0,0,1,0,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,0,1,0,1,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,0,1,1,1,0,0,1,1,1)
binomialGLMM<-data.frame(Nest,Treatment,Response)
# Binomial GLMM
library(lme4)
lmm1<-lmer(Response~Treatment+(1|factor(Nest)),data=binomialGLMM,family=binomial)
lmm2<-lmer(Response~1+(1|factor(Nest)),data=binomialGLMM,family=binomial)
anova(lmm1,lmm2)
Data: binomialGLMM
Models:
lmm2: Response ~ 1 + (1 | factor(Nest))
lmm1: Response ~ Treatment + (1 | factor(Nest))
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
lmm2 2 109.405 114.314 -52.703
lmm1 7 57.013 74.193 -21.506 62.392 5 3.889e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Post-hoc pairwise comparison
library(multcomp)
posthoc<-glht(lmm1, linfct = mcp(Treatment = "Tukey"))
summary(posthoc)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glmer(formula = Response ~ Treatment + (1 | factor(Nest)), data =
binomialGLMM,
family = binomial)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
!!! *CRL - CL == 0 -2.088e+01 4.581e+03 -0.005 1.00000*
!!! *CRS - CL == 0 -2.247e+01 4.581e+03 -0.005 1.00000*
CS - CL == 0 -1.652e+01 4.581e+03 -0.004 1.00000
SL - CL == 0 9.148e-07 6.276e+03 0.000 1.00000
SS - CL == 0 -1.735e+01 4.581e+03 -0.004 1.00000
CRS - CRL == 0 -1.595e+00 1.348e+00 -1.183 0.79608
CS - CRL == 0 4.359e+00 1.359e+00 3.206 0.01107 *
!!! *SL - CRL == 0 2.088e+01 4.290e+03 0.005 1.00000 *
SS - CRL == 0 3.531e+00 1.072e+00 3.294 0.00790 **
CS - CRS == 0 5.953e+00 1.700e+00 3.502 0.00382 **
!!! *SL - CRS == 0 2.247e+01 4.290e+03 0.005 1.00000*
SS - CRS == 0 5.125e+00 1.480e+00 3.464 0.00421 **
SL - CS == 0 1.652e+01 4.290e+03 0.004 1.00000
SS - CS == 0 -8.279e-01 1.465e+00 -0.565 0.98998
SS - SL == 0 -1.735e+01 4.290e+03 -0.004 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
!!! on the left indicates unexpected p-values.
--
View this message in context: http://r.789695.n4.nabble.com/Help-with-glht-function-and-binary-data-tp4650222.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list