Solved – Tukey with interaction and extracting letters in multcomp

lme4-nlmepost-hoctukey-hsd-test

I've seen similar questions posted elsewhere with no answer, so I thought I'd put this out into the CrossValidated universe.

I am attempting to do a Tukey's post-hoc test on a model that has an interaction. To do this, I have used the interaction() function to create a new variable that combines my two main effects, Habitat (3 levels) and Detritus (2 levels) into one variable. 

   reg.invert$interaction <- interaction(reg.invert$Habitat, reg.invert$Detritus)

Thus, I have the following model:

   model <- glmer (Species Richness ~ interaction + (1|Site)).

I have been trying to do this using glht() function in the package multcomp.

This is what I ran, with notes about the output:

 inter.test1 <- glht(model,  mcp(interaction = "Tukey"))
 summary(inter.test1)
 cld(inter.test1)

This resulted in 5 warnings that said

In RET$pfunction("adjusted", …) : Completion with error > abseps.

Why am I getting these warnings? Should I not use the results of cld(), given these warnings? How can I avoid the warnings? In looking at the raw data, some of the letters seem a bit off, so I was wondering if I can trust this output.

I have also run the exact same model with a different response variable without any errors, and the extracted letters make sense.

I would be very grateful for any advice for how to successfully use Tukey's method in this model, as well as for how to extract letters to denote significant differences. Either suggestions for adjusting my code or for a new approach entirely would be fine.

Thanks for your help!

Edit: See below for what I hope is a reproducible example using my data…please bear with me- I've never made one of these before:

 library(lme4)
 library(multcomp)


  ####re-create my crazy dataframe:
 new.data<-structure(list(Site = structure(c(13L, 13L, 13L, 13L, 13L, 13L, 
                              13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 12L, 
                              12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 
                              12L, 12L, 12L, 12L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 
                              11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 10L, 10L, 10L, 10L, 
                              10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
                              10L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
                              9L, 9L, 9L, 9L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
                              8L, 8L, 8L, 8L, 8L, 8L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
                              7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
                              6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 5L, 5L, 5L, 5L, 5L, 
                              5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 4L, 4L, 4L, 
                              4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 2L, 2L, 
                              2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                              1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                              3L, 3L), .Label = c("JM", "SL", "AV", "SEB", "CANA", "SC", "NP", 
                                                  "GTM", "NERR", "TAL", "JEK", "CAT", "SKID"), class = "factor"), 
           Habitat = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("MAN", 
                                                                                             "MIX", "SM"), class = "factor"), Detritus = structure(c(2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
                                                                                                                                                     3L, 1L, 2L, 3L, 1L), .Label = c("CON", "AVI", "SPAR"), class = "factor"), 
           Species.Number = c(10L, 7L, 10L, 5L, 7L, 0L, 6L, 6L, 11L, 
                              5L, 4L, 3L, 5L, 5L, 1L, 4L, 7L, 4L, 7L, 8L, 6L, 10L, 11L, 
                              6L, 7L, 5L, 7L, 10L, 10L, 5L, 11L, 8L, 6L, 8L, 10L, 9L, 4L, 
                              6L, 5L, 7L, 5L, 3L, 7L, 8L, 3L, 7L, 5L, 5L, 4L, 3L, 4L, 7L, 
                              6L, 5L, 8L, 3L, 4L, 5L, 7L, 3L, 6L, 6L, 6L, 6L, 7L, 4L, 6L, 
                              6L, 5L, 6L, 8L, 8L, 6L, 4L, 7L, 11L, 6L, 6L, 4L, 3L, 1L, 
                              7L, 5L, 4L, 9L, 6L, 6L, 3L, 4L, 3L, 8L, 11L, 9L, 7L, 10L, 
                              6L, 6L, 11L, 8L, 4L, 6L, 5L, 5L, 6L, 5L, 6L, 5L, 2L, 2L, 
                              3L, 1L, 3L, 1L, 0L, 3L, 5L, 3L, 6L, 4L, 2L, 3L, 6L, 3L, 6L, 
                              10L, 4L, 5L, 6L, 4L, 8L, 5L, 1L, 6L, 4L, 4L, 4L, 4L, 4L, 
                              4L, 3L, 1L, 6L, 3L, 3L, 5L, 4L, 5L, 5L, 5L, 4L, 9L, 9L, 6L, 
                              11L, 4L, 5L, 7L, 8L, 4L, 11L, 17L, 1L, 3L, 2L, 3L, 5L, 2L, 
                              3L, 5L, 2L, 4L, 2L, 3L, 4L, 3L, 4L, 3L, 9L, 4L, 0L, 2L, 2L, 
                              2L, 2L, 6L, 2L, 1L, 3L, 1L, 5L, 3L, 1L, 6L, 4L, 1L, 3L, 4L, 
                              4L, 4L, 8L, 2L, 3L, 5L, 1L, 1L, 3L, 2L, 6L, 8L, 4L, 3L, 3L, 
                              1L, 1L, 1L, 2L)), .Names = c("Site", "Habitat", "Detritus", 
                                                           "Species.Number"), row.names = c(NA, -216L), class = "data.frame")



 ####analysis:


 new.data$interaction<-interaction(new.data$Habitat,new.data$Detritus)
 m6<-glmer(Species.Number~interaction + (1|Site), family = "poisson", data = new.data, control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 20000)))


 inter.test1<-glht(m6, mcp(interaction="Tukey"))
 plot(inter.test1)
 summary(inter.test1)
 cld(inter.test1, level=0.05)

Best Answer

I also cannot seem to figure out what the warning is trying to tell us but one alternative approach would be using the emmeans() function in the emmeans package and then the cld() function:

inter.test1 <- emmeans(m6, "interaction")
#NOTE: I would probably rename "interaction" to something else as this word is often used
#in other function arguments or even as a function itself, i.e. interaction()!
cld(inter.test1, Letter="abcdefg")
#interaction    emmean        SE  df asymp.LCL asymp.UCL .group
#MAN.CON     0.7938065 0.1786641 Inf 0.4436314  1.143982  a    
#MAN.AVI     1.2638100 0.1502522 Inf 0.9693211  1.558299  abc  
#MIX.CON     1.2995688 0.1335088 Inf 1.0378964  1.561241  ab   
#MAN.SPAR    1.3096194 0.1479132 Inf 1.0197148  1.599524  abc  
#SM.CON      1.5982275 0.1094568 Inf 1.3836962  1.812759   bcd 
#MIX.AVI     1.7414015 0.1177577 Inf 1.5106006  1.972202    cd 
#MIX.SPAR    1.8103945 0.1157236 Inf 1.5835804  2.037209    cd 
#SM.SPAR     1.8133390 0.1034105 Inf 1.6106582  2.016020   bcd 
#SM.AVI      1.8908973 0.1014570 Inf 1.6920453  2.089749     d 

#Results are given on the log (not the response) scale. 
#Confidence level used: 0.95 
#P value adjustment: tukey method for comparing a family of 9 estimates 
#significance level used: alpha = 0.05

If you want to see the mean estimates on the response scale, i.e. back-transformed from the log scale, you can add type = response:

inter.test2 <- emmeans(m6, "interaction", type = "response")
cld(inter.test2, Letter="abcdefg")

See also the very helpful emmeans vignettes.

Related Question