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?
@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?
by = NULLI 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.nbto avoid thinking too hard about the simulation :P $\endgroup$