0

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)))

$emmeans
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
etc.

Any help greatly appreaciated.

guest
  • 1
  • The `multcomp` package provides `glht` which should allow to test for this type of contrasts. See also [this application](https://stats.stackexchange.com/a/125321/930). – chl Nov 22 '20 at 17:38
  • 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...? – guest Nov 22 '20 at 21:50
  • 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. – Russ Lenth Nov 23 '20 at 03:41
  • You need to write up the interaction explicitely, either using `interaction()`, or as suggested by @Russ. – chl Nov 23 '20 at 16:48
  • @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. – Russ Lenth Nov 23 '20 at 18:08
  • @RussLenth Thank you so much! You've explained it perfectly for me and thank you for emmeans! – guest Nov 23 '20 at 19:26

0 Answers0