4
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ There seems to be some typos in your code, e.g.: In the model you use Interaction while here reg.invert$interaction so one is capitalized the other is not. I tried on an example dataset and for me everything works fine. Maybe provide a reproducible example. $\endgroup$
    – Stefan
    Commented Feb 13, 2018 at 22:13
  • $\begingroup$ Awesome, thanks so much for your help! I just added what I hope is a reproducible example using my dataset to the end of the post. Any suggestions would be great! $\endgroup$ Commented Feb 13, 2018 at 22:59

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Hi Stefan- thanks for this suggestion! Any ideas on why the df = Inf in the emmeans output? Also, from reading one of the EMM vignettes, they state that they "really don’t recommend this method, though, as it imposes a stark difference between P values slightly less and slightly more than alpha." Any additional suggestions on this? I really want to know what that warning means from the multcomp package... $\endgroup$ Commented Feb 15, 2018 at 14:56
  • $\begingroup$ Here is an explanation what the df=Inf means in the output: cran.r-project.org/web/packages/emmeans/vignettes/models.html#G Regarding the other part of your question GLMMs by default provide Wald Z-statistics, which are asymptotic approximations. Also have a look here and search for degrees of freedom bbolker.github.io/mixedmodels-misc/glmmFAQ.html I am sure multcomp uses the same asymptotic approximations for its p-values. See also answer+comments here: stats.stackexchange.com/questions/204741/… $\endgroup$
    – Stefan
    Commented Feb 15, 2018 at 21:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.