planned contrast with **emmeans** on interaction term in mixed model

Question

I think this is simple: I want to test if a treatment (that I know a priori does something) is stronger in one group vs. the other (in the presence of random effects). This seems like it should be straightforward with emmeans, but I'm struggling to find the example that matches my use case and would love a pointer. Below is a simulated dataset and where I got to so far...

How do I compare the strength of treatment between my two groups f2 == "A" and f2 == "B"?

I am partly asking for the right emmeans syntax, but also not sure which p-value correction is appropriate in a situation like this. Thanks!

# load libraries
library(tidyverse)
library(lme4)
library(emmeans)

# simulate data
# this is a slight modification from glmer.nb test data from Ben Bolker
dd <- expand.grid(f1 = factor(1:2), # changed to 2 from 3; this is treatment 
                  f2 = LETTERS[1:2] # this is the "group"
                  , g = 1:9 # random effect
                  , rep = 1:600
                  , KEEP.OUT.ATTRS = FALSE)
# in documentation, BB didn't simulate `mu` as a function of `g` but I will
# this avoids singular fits
summary(mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2) + g/5)))
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
str(dd)

# define and fit model
my_formula <- "y ~ f1:f2 + f1 + f2 + (1|g)"

my_mod <- glmer.nb(formula = my_formula
                   , data = dd
                   # , control = glmerControl(maxit = 1e6)
)

# what I want to say is whether the treatment("f1") is stronger for one group 
# vs. the other (i.e `f2 == "A"` vs. `f2 == "B"`)

# I thought the way to do this (in a planned contrast framework) might be to 
# estimate marginal means (with `emmeans::emmeans`) and then calculate adjusted
# p-values for the interaction term of interest here

# my first stab was this:

emm_first<-emmeans(my_mod
                   ,  trt.vs.ctrl~ f1|f2)

# then, to get the contrast
contrast(emm_first
         , method = "trt.vs.ctrl")

# which returns 0s. Is this right (and I simulated wrong), or am I using emmeans 
# wrong?

UPDATE

@RussLenth suggested modifying the last line to

contrast(emm_first
        , method = "trt.vs.ctrl"
         , by = NULL)

Which returns great-looking ouptput:

$emmeans
 contrast      estimate     SE  df z.ratio p.value
 f12 A - f11 A    0.377 0.0277 Inf  13.603  <.0001
 f11 B - f11 A    1.124 0.0276 Inf  40.692  <.0001
 f12 B - f11 A    1.252 0.0276 Inf  45.385  <.0001

Results are given on the log (not the response) scale. 
P value adjustment: dunnettx method for 3 tests 

$contrasts
 contrast                      estimate     SE  df z.ratio p.value
 (f12 - f11 B) - (f12 - f11 A)   -0.248 0.0389 Inf  -6.374  <.0001

Results are given on the log (not the response) scale. 

Final question: is this "dunnettx" p-value the best choice for asking my question?

Answer 1

I can't reproduce the results here exactly because the random-number seed is not given. But there is some confusion here because, due to the two-sided formula, emm_first is a list of two objects. If I understand the question, we want the contrasts only of the second one -- i.e., a contrast of contrasts:

> contrast(emm_first[[2]]
+          , method = "trt.vs.ctrl"
+          , by = NULL)
 contrast                      estimate     SE  df z.ratio p.value
 (f12 - f11 B) - (f12 - f11 A)   -0.221 0.0392 Inf  -5.631  <.0001

Results are given on the log (not the response) scale. 

This is an example of an interaction contrast. Below are the means and the required interaction contrast:

> (EMM <- emmeans(my_mod, ~ f1 * f2))
 f1 f2 emmean     SE  df asymp.LCL asymp.UCL
 1  A    2.28 0.0580 Inf      2.17      2.40
 2  A    2.69 0.0580 Inf      2.57      2.80
 1  B    3.37 0.0579 Inf      3.25      3.48
 2  B    3.55 0.0579 Inf      3.43      3.66

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

> contrast(EMM, interaction = "trt.vs.ctrl")
 f1_trt.vs.ctrl f2_trt.vs.ctrl estimate     SE  df z.ratio p.value
 2 - 1          B - A            -0.221 0.0392 Inf  -5.631  <.0001

Results are given on the log (not the response) scale. 

Note that the same result is obtained.

This is discussed more fully in the vignette on interactions in the emmeans package.

Side note: Due to fairly widespread confusion resulting from follow-up analyses when a two-sided formula was used, versions after 1.7.4-1 will just do the contrast on the first one in the list (in this case, not the one desired), with a message.