there are a lot of questions about post-hoc tests for GLMMs on this site and thanks to the replies I almost have my question solved. I want to see if there is a difference in treatment groups over time but for all pairwise comparisons.

My model has a count response (count) with a categorical predictor treatment (as factor with 3 levels: A,B,C) and year (repeated measures of count over time). I also include an interaction to test whether the treatment levels differ over time as follows (simplified here, removed random effects for brevity):

fit <- glmmTMB(count ~ treatment + year + year:treatment)

Using the posts here and here (along with the emmeans vignettes) I have contrasts between groups for each year separately:

emmeans(fit, pairwise ~ treatment | year, at = list(year = c(1, 2, 3,4,5,6)))

year = 1:
 treatment emmean    SE  df lower.CL upper.CL
 A   13.55 0.276 423    13.00     14.1
 B   13.51 0.276 423    12.96     14.0
 C   13.24 0.276 423    12.69     13.8
year = 6:
 treatment emmean    SE  df lower.CL upper.CL
 A   12.90 0.241 423    12.43     13.4
 B   12.89 0.241 423    12.42     13.4
 C   12.58 0.241 423    12.11     13.1

But I want all pairwise contrasts for treatment:year to compare treatments within and across all years. This is easy if it were a linear model followed by Tukey's test. The result looks something like:

    A:year1 - B:year1
    A:year1 - B:year2
    A:year1 - B:year3
    A:year2 - B:year2
    A:year2 - B:year3

Any help greatly appreaciated.

  • $\begingroup$ The multcomp package provides glht which should allow to test for this type of contrasts. See also this application. $\endgroup$
    – chl
    Commented Nov 22, 2020 at 17:38
  • $\begingroup$ Thanks for the reply. Hmmm...I'm still not able to figure out how glht provides contrasts by treatment:year. I again only get contrasts by treatment, or whatever the factor levels are...? $\endgroup$
    – guest
    Commented Nov 22, 2020 at 21:50
  • $\begingroup$ How does this differ from pairwise ~ treatment * year? (In other words, don't make year a by variable). That will do all pairs from those 18 combinations (warning -- that is 171 comparisons, and the Tukey correction is very conservative. $\endgroup$
    – Russ Lenth
    Commented Nov 23, 2020 at 3:41
  • $\begingroup$ You need to write up the interaction explicitely, either using interaction(), or as suggested by @Russ. $\endgroup$
    – chl
    Commented Nov 23, 2020 at 16:48
  • $\begingroup$ @chl @guest the approach using interaction()' requires starting from scratch: defining that variable, fitting a new model with that variable as the one predictor, and running glht() or emmeans(). glht() is really not very easy to use except for one-factor models, and that's one of the main reasons I wrote emmeans. BTW you can also use glht but specify an emmeans::emm() call for its linfct argument. $\endgroup$
    – Russ Lenth
    Commented Nov 23, 2020 at 18:08


Your Answer

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