I'm confused about using an appropriate post hoc pairwise comparison test in the presence of a significant 4-way interaction effect in the model.

Here are the details of the experiment: I want to test the interaction effect of two learning conditions on the two groups of subjects and I have an additional within-subject factor, age.

My data includes 1 dependent variable- reaction time and 4 independent variables of which three(Block, SeqTyp, Age) are within-subject factors and 1(Group)is a between-subject factor.

I used R for data analysis and ran lmer model:

lmer(RTNormalized~ Block* SeqType* Group* Age + (1+Block*SeqType|ParticipantID)

Results showed a significant 4-way interaction effect. Further, I used em-means for the pair-wise comparison

> emm1 <- emmeans(fnl.confir.mod1, pairwise ~ Block* SeqType* Group* Age)

But, R gives a warning message that it's a complicated computation with a large number of observations and to adjust the limits by adding argument 'pbkrtest.limit = 23054' and 'lmerTest.limit = 23054' or larger) and warns of large computation time.

I wonder how to conduct a pairwise comparison. I'm interested in testing which among Blocks (2 levels) and Seq.type (2 levels) interactions per group are significant.

I request your suggestions on resolving em mean complex computation error in this situation and to help me answer if it's appropriate to use marginal means estimates of emmeans in the presence of significant 4way interaction effects


1 Answer 1


I would say first that doing pairwise comparisons of all combinations of those four factors is certainly ill-advised. Even if each of those factors has only 2 levels, you are already proposing to test 16x15/2 = 120 comparisons. Most likely, it is hugely more than that. You should first do

emm <- emmeans(fnl.confir.mod1, ~ Block* SeqType* Group* Age,
               lmer.df = "asymp")

(omitting pairwise -- and I advise you to always omit that and instead do contrasts/comparisons in separate step(s)). (The additional lmer.df argument is discussed below.) I would suggest that you do well-chosen simple comparisons based on what is needed for interpretation of effects. e.g.,

contrast(emm, "pairwise", simple = "Group")
# and probably other contrast statements

There is also such a thing as interaction contrasts, e.g.

contrast(emm, interaction = "consec")
# or
contrast(emm, interaction = "consec", by = c("Group", "Age"))

The first will provide 4-way interaction contrasts equal in number to the number of degrees of freedom for the four-way interaction. The second one obtains 2-way interaction contrasts of Block and SeqType for each combination of the other factors. Again, the above are examples, and you should choose what contrasts to test based on considerations like which interactions are not so important.

You can also get a graphical overview of the model predictions, which may p[rove helpful for understanding what is going on:

emmip(emm, Group ~ Age | Block * SeqType)

(though you may want to organize these differently by switching factors around)

The message you got is another issue altogether; it has to do with the computations needed for degrees of freedom based on the Kenward-Roger method. The computational burden for this method becomes prohibitive for large datasets -- that's why you get the warning message. I suggest just adding lmer.df = "asymp" (as is done above) to bypass those computations. See a discussion of d.f. issues for these models and also the Models reference.

  • $\begingroup$ Also I wonder about the model, and if you need a model with so many random slopes. I might guess that Block is a random effect, as blocks are typically a restriction on randomization. -- if so, it should not be in the fixed part of the model, and typically, blocks are not modeled as interacting with fixed effects. An error term of (1|ParticipantID) + (1 | Block/SeqType) would model participants, blocks, and block*sequences as random intercepts, which may well be enough. $\endgroup$
    – Russ Lenth
    Commented Jan 11, 2023 at 0:44
  • $\begingroup$ Thank you for the detailed explanation and for suggesting corrections to the model, Russ Lenth $\endgroup$
    – Arpitha
    Commented Jan 11, 2023 at 9:00
  • $\begingroup$ I tried--RT ~ Block * SeqType* Group*Age + (1 | ParticipantID) + (1 | Block/SeqType)-- model. However, this specification turns out to have a singularity issue. I'm not sure how best to resolve it. The reason we need to include random slopes of blocks in this experiment is that the "learning" is predicted by the difference between 2 block conditions and I hypothesise this would vary by SeqType per participant. But you rightly pointed out the wrong model specification with the block to be a fixed random slope. $\endgroup$
    – Arpitha
    Commented Jan 11, 2023 at 9:23
  • $\begingroup$ I have a further follow-up question on your suggestion to use the interaction contrasts, the code- contrast(emm, interaction = "consec") resulted in comparisons of two levels of blocks with two levels of SeqTyp (e.g., rand2 - pattern5 SeqTyp1 - SeqTyp2 estimate 0.2355; SE= 0.0361; df= Inf; z.ratio=6.527; p.value <.0001), how to interpret this? I wonder is there a way to extract the interaction of blocks effect for each SeqType condition seperatley? $\endgroup$
    – Arpitha
    Commented Jan 11, 2023 at 9:37
  • $\begingroup$ The rand5-pattern5, SeqTyp1-SeqTyp2 interaction contrast is a contrast of contrasts: (rand5-pattern5 | SeqTyp1) - (rand5-pattern5 | SeqTyp2). On 2nd Q, use `contrast(emmeans(..., ~Block*whatever|SeqType), interaction = "consec")`` $\endgroup$
    – Russ Lenth
    Commented Jan 11, 2023 at 20:54

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.