I have a rookie question about emmeans in R.

I fit a complex model using lmer() with the following variables:

  • A: a binary categorical predictor, within-subject
  • B: a binary categorical predictor, within-subject
  • C: a categorical predictor with 4 levels, between-subject
  • X & Y: control variables of no interest, one categorical, one continuous.

The model is as follows:

fit1 <- lmer(rt ~ 1 + A*B*C + X + Y + (1+A*B|Subject))

Now I'm mostly interested in how the A*B interaction differs across different levels of C (i.e., is the interaction different across the four groups I have). I was trying to use emmeans to get to the bottom of this, and I have found some very useful threads here on CrossValidated, but I cannot seem to find one that I can generalize easily to my case.

Here's what I did: I created a new model with an interaction term (AB = A*B).

fit1b <- lmer(rt ~ 1 + A*C + B*C + AB*C + X + Y + (1+A*B|Subject))

Then used emmeans like this:

emms <-  emmeans(fit1b, ~ AB*C)
contrast(emms, interaction = "pairwise")

This results in an output that seems to make sense, however, I'm really uncertain about whether this setup makes any sense to begin with. Essentially my goal is to be able to determine whether the A*B interaction in greater in group x compared to group y, while controlling for all the other stuff in the model.

Is this a good way of doing this? Is there an easier/nicer way to do this?

EDIT: I created a simulated data set - here it is: https://osf.io/4cr8x A is actually congruency in a conflict task, and B is previous trial congruency that's why the first trials have no value in that column. C is the group variable, just like in the example above. X and Y are the controls, with X being trial number and Y being sex.

EDIT 2: And here's the exact code I ran on the simulated data:


Data <- read.csv("simdat.csv",header=TRUE, sep=",", na.strings="-999", dec=".", strip.white=TRUE)

Data$A <- as.factor(Data$A)
Data$B <- as.factor(Data$B)
Data$C <- as.factor(Data$C)
Data$Y <- as.factor(Data$Y)
Data$Subject <- as.factor(Data$subject)

fit1 <- lmer(rt ~ 1 + A*B*C + X + Y + (1|Subject), data = Data, verbose = 0, REML = F) #I simplified the random structure as the original wouldn't converge with the simulated data

interaction_term <- (as.numeric(levels(Data$A))[Data$A])*(as.numeric(levels(Data$B))[Data$B])
Data$AB <- as.factor(interaction_term) 

fit2 <- lmer(rt ~ 1 + A*C + B*C + AB*C + X + Y + (1|Subject), data = Data, verbose = 0, REML = F)

emms <-  emmeans(fit2, ~ AB*C)
contrast(emms, interaction = "pairwise")  
  • 1
    $\begingroup$ can you share the data (or some of it)? It would be easier to demonstrate with data $\endgroup$
    – Mark White
    Jul 11, 2018 at 14:29
  • 2
    $\begingroup$ FWIW this shouldn't be mixed-model-specific; I would think the question would be equally relevant (and the solution would be the same) for lm(rt ~ 1 + A*C + B*C + AB*C + X + Y) $\endgroup$
    – Ben Bolker
    Jul 11, 2018 at 14:40
  • $\begingroup$ Thanks for the suggestion. I quickly simulated some data and added an osf link to the xls file. I updated the question. Can you see it? $\endgroup$
    – MGy
    Jul 11, 2018 at 17:26
  • $\begingroup$ I added the exact code I ran too. $\endgroup$
    – MGy
    Jul 12, 2018 at 9:03

1 Answer 1


It shouldn't be necessary to fit a separate model just to do the post-hoc comparisons you want. You had tried:

emms <-  emmeans(fit1b, ~ AB*C)
contrast(emms, interaction = "pairwise")

but you can get the same results from the original model using by variables judiciously:

emms1 <- emmeans(fit1, ~ A*B | C)
con1 <- contrast(emms1, interaction = "pairwise")
pairs(con1, by = NULL)

The con1 results are the desired 1-d.f. interaction effects for each level of C (the by factor is remembered). Then we compare them pairwise, no longer using the by grouping. By default, a Tukey adjustment is made to the family of comparisons, but you may use a different method via adjust.

  • $\begingroup$ Thank you, this is a fantastic reply, this looks like exactly what I need. Before I accept it, could you clarify how to read the output? E.g., the first line is: A0 - A1,B0 - B1,C1 - A0 - A1,B0 - B1,C2 - is this then, the difference in the A*B interaction between groups C1 and C2? (To make the output clearer, I labelled each factor level.) Similarly, line two, A0 - A1,B0 - B1,C1 - A0 - A1,B0 - B1,C3 is the difference between the first and the third group, is that correct? $\endgroup$
    – MGy
    Jul 14, 2018 at 14:21
  • 1
    $\begingroup$ Yes, that’s correct. If you look at con1, it’s pairwise comparisons of those. I agree the labeling is confusing; a result of multiple steps. $\endgroup$
    – Russ Lenth
    Jul 14, 2018 at 14:58

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.