[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