Pairwise comparisons via emmeans

Question

I have a question about emmeans and mixed effect model. I fit a complex model using lmer() with the following variables:

A: a binary categorical predictor, between-subject

B: a binary categorical predictor, within-subject

C: a categorical predictor with 3 levels, within-subject

The model is as follows:

library(lme4)
model= lmer(d ~ A * B * C + (1|subject), data = ddata, REML = FALSE) 

The results show a significant interaction between the predictors. I need to know the pairwise comparisons so I use emmeans as follow:

library(emmeans)
emmeans(model, pairwise ~ A * B * C)

emmeans use the Tukey method for the pairwise comparisons. A reviewer of my paper does not like Tukey method as she said:

"Tukey’s HSD assumes all measures are independent from each other and is not appropriate for repeated measures. Please consider performing a different test even if that changes the results."

My questions are:

Is the claim that Tukey in this case is not appropriate? if no can you cite any source?

if yes what is the alternative method?

Answer 1

Your reviewer has a good point, because you have both between- and within-subject comparisons in your collection, and those have different standard errors.

But the other thing I notice is that there are a total of $2\times2\times3=12$ cell means, and the emmeans() call shown produces all ${12\choose2}=66$ possible comparisons. I truly wonder if all of those are really of interest. If you do, I suggest trying the "mvt" adjustment, which is the same idea as Tukey but is based on the multivariate $t$ distribution with the actual correlation structure in your model (which in this case is not the same as the correlation structure assumed by the Tukey method).

However, often when there is an interaction, people opt for "simple" comparisons -- in this case, comparing the levels of one factor while holding the other two fixed. Those are easily done via

emm <- emmeans(model, ~ A * B * C)
simp <- pairs(emm, simple = "each")
simp

This will yield 6 comparisons of the levels of A, 6 comparisons of the two levels of B, and 4 sets of 3 comparisons among the levels of C, for a total of 24 comparisons instead of 66. Moreover, the issues of Tukey being inappropriate go away, because each set of simple comparisons is homogeneous.

Some additional comments:

• While REML = FALSE is often a good idea for testing your model against other models, I recommend refitting your model with REML = TRUE before proceding to post hoc comparisons, because the REML method reduces bias in the estimates.
• Consider using meaningful names in place of A, B, and C. We're not doing a generic math problem here; we are trying to tell a useful story about an actual experiment. Meaningful names help everybody understand what you're talking about.
• I suggest doing things in steps, as shown above, over trying to get every result you want in one R call. That makes it more natural to focus on particular results or go in different directions; e.g., summary(emm) shows the 12 cell means and pairs(emm, by = “C”) could be used to compare the four A:B combinations at each level of C.
• You might want to do a stronger multiplicity correction rather than a separate one for each of the comparisons. For example, test(simp[[1]], by = NULL, adjust = "mvt") puts all 6 of the A comparisons in one family and applies the multivariate $t$ adjustment. (The Tukey adjustment is completely inappropriate for that because it is not a set of pairwise comparisons; in fact, the software doesn't even allow that adjustment.)