1
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ emm_first still "remembers" f2 as the by variable. try adding by = NULL to the contrast call. $\endgroup$
    – Russ Lenth
    Jun 1, 2022 at 17:12
  • 1
    $\begingroup$ BTW, I wonder why you named the factors "f1" and "f2" rather than "treatment" and "group". I don't think it improves an analysis to make the variable names abstract. In fact, it puts up a barrier in interpreting the results. $\endgroup$
    – Russ Lenth
    Jun 1, 2022 at 17:14
  • $\begingroup$ @RussLenth Thanks! When adding by = NULL I get an output that looks right! Is this a sensible p-value for a planned contrast? FWIW I agree on variable names. I basically copy-pasted the example given in ?lme4::glmer.nb to avoid thinking too hard about the simulation :P $\endgroup$ Jun 1, 2022 at 19:13

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Thanks for this answer, which fully addresses how to use emmeans to obtain the relevant contrast. Could you also comment on the correct p-value computation (which method is preferred for this and why)? $\endgroup$ Jun 2, 2022 at 15:26
  • $\begingroup$ Really appreciate all the support you give on this site to help people learn your package, it's super. Thank you! $\endgroup$ Jun 2, 2022 at 15:30
  • $\begingroup$ In the vignette, it seems contrasts are considered post-hoc, I'd like to leverage any power I can from having a well-formulated question :-). What is the least-conservative, appropriate p-value for a pre-hoc planned comparison of an interaction term? Thanks!! $\endgroup$ Jun 2, 2022 at 15:36
  • 3
    $\begingroup$ The one you have. It is just one test, so there is no multiplicity adjustment. There is no law that says you have to do an ANOVA first -- those really are most useful for model selection, and their importance for testing (as opposed to just summarizing the contributions of model terms) is over-rated, IMO. $\endgroup$
    – Russ Lenth
    Jun 2, 2022 at 15:49

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.