# Posthoc analysis LMM for interaction containing a factor with 3 levels

I am fitting a Linear Mixed Model using the lmerTest package. I have two fixed effects: group (categorical, 2 levels) and period (categorical, 3 levels). I have an interaction term between the two fixed effects period*group The response variable is continuous. I have one random effect (id). Please see the data below:

> moddf <- dput(mydata)
structure(list(id = c("id3", "id5", "id3", "id5", "id3", "id5",
"id1", "id2", "id4", "id6", "id7", "id8", "new", "id1", "id2",
"id4", "id6", "id7", "id8", "new", "id1", "id2", "id4", "id6",
"id7", "id8", "new"), area = c(28.28, 11.04, 31, 10.28, 18.4,
18.12, 17.56, 18.28, 13.04, 14.24, 21.44, 26.48, 22.56, 7.92,
16.28, 8.52, 16.8, 7.36, 11.56, 17.64, 4.56, 5.04, 8.16, 6.2,
3.24, 5.28, 6.6), period = structure(c(1L, 1L, 2L, 2L, 3L, 3L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L,
3L, 3L, 3L, 3L, 3L), levels = c("night", "outside", "within"), class = "factor"),
group = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L), levels = c("groupN", "groupP"), class = "factor")), row.names = c(NA,
-27L), class = "data.frame")


I specified the following model:

lmR_g <- lmerTest::lmer(area ~ group + period + period*group + (1|id), data = moddf, REML = F)


The summary shows that the interaction is significant. And visual representation asks for a post hoc analysis as the effects of period clearly differs in the two groups (see plot). While there are different options described online, I am unsure which one is most appropriate. I have used lsmeans so far to test for each period at the levels of the group (see below). However, any advice regarding post hoc options for this situation as well as advice if the $lsmeans or$contrasts should be used if lsmeans() is applied would be highly appreciated.

lsmeans(lmR_g, pairwise ~ period | group)

$lsmeans group = groupN: period lsmean SE df lower.CL upper.CL night 19.66 3.89 29.9 11.71 27.61 outside 20.64 3.89 29.9 12.69 28.59 within 18.26 3.89 29.9 10.31 26.21 group = groupP: period lsmean SE df lower.CL upper.CL night 19.09 2.08 29.9 14.83 23.34 outside 12.30 2.08 29.9 8.05 16.55 within 5.58 2.08 29.9 1.33 9.83 Degrees-of-freedom method: kenward-roger Confidence level used: 0.95$contrasts
group = groupN:
contrast         estimate   SE   df t.ratio p.value
night - outside     -0.98 4.65 23.1  -0.211  0.9759
night - within       1.40 4.65 23.1   0.301  0.9515
outside - within     2.38 4.65 23.1   0.511  0.8666

group = groupP:
contrast         estimate   SE   df t.ratio p.value
night - outside      6.79 2.49 23.1   2.729  0.0309
night - within      13.50 2.49 23.1   5.427  <.0001
outside - within     6.71 2.49 23.1   2.699  0.0329

Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 3 estimates


First, it's worth considering whether your interaction truly is significant. Looking at the summary() of your model:

> summary(lmR_g)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: area ~ group + period + period * group + (1 | id)
Data: moddf

AIC      BIC   logLik deviance df.resid
176.0    186.3    -80.0    160.0       19

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.6776 -0.6796 -0.1875  0.8329  1.6776

Random effects:
Groups   Name        Variance Std.Dev.
id       (Intercept)  6.729   2.594
Residual             16.849   4.105
Number of obs: 27, groups:  id, 9

Fixed effects:
Estimate Std. Error       df t value Pr(>|t|)
(Intercept)                19.6600     3.4335  23.2179   5.726 7.58e-06 ***
groupgroupP                -0.5743     3.8933  23.2179  -0.148   0.8840
periodoutside               0.9800     4.1048  18.0000   0.239   0.8140
periodwithin               -1.4000     4.1048  18.0000  -0.341   0.7370
groupgroupP:periodoutside  -7.7686     4.6544  18.0000  -1.669   0.1124
groupgroupP:periodwithin  -12.1029     4.6544  18.0000  -2.600   0.0181 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) grpgrP prdtsd prdwth grpgrpP:prdt
groupgroupP  -0.882
periodoutsd  -0.598  0.527
periodwithn  -0.598  0.527  0.500
grpgrpP:prdt  0.527 -0.598 -0.882 -0.441
grpgrpP:prdw  0.527 -0.598 -0.441 -0.882  0.500


You can see that your model adds two interaction terms, one of which is p<0.05. To test whether the interaction as a whole is significant, you can consider adopting more of an ANOVA framework and testing whether the addition of the interaction term improves model fit:

> lmR_g_null = lmerTest::lmer(area ~ group + period  + (1|id), data = moddf, REML = F)
> anova(lmR_g_null,lmR_g)
Data: moddf
Models:
lmR_g_null: area ~ group + period + (1 | id)
lmR_g: area ~ group + period + period * group + (1 | id)
npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
lmR_g_null    6 177.84 185.61 -82.920   165.84
lmR_g         8 175.97 186.33 -79.984   159.97 5.8721  2    0.05308 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


The chisquared test of the log-likelihood ratio isn't quite significant, and while the AIC is smaller in the full model, the BIC is not. So you may not have an interaction here.

Even so, if you did have an interaction, the way your lsmeans() call is structured only indirectly tells you about it's source. In your example the results reflect the effect of effect of period within each group, instead you want the effect of group at each level of period:

> lsmeans(lmR_g, pairwise ~ group | period)
$lsmeans period = night: group lsmean SE df lower.CL upper.CL groupN 19.66 3.89 29.9 11.71 27.61 groupP 19.09 2.08 29.9 14.83 23.34 period = outside: group lsmean SE df lower.CL upper.CL groupN 20.64 3.89 29.9 12.69 28.59 groupP 12.30 2.08 29.9 8.05 16.55 period = within: group lsmean SE df lower.CL upper.CL groupN 18.26 3.89 29.9 10.31 26.21 groupP 5.58 2.08 29.9 1.33 9.83 Degrees-of-freedom method: kenward-roger Confidence level used: 0.95$contrasts
period = night:
contrast        estimate   SE   df t.ratio p.value
groupN - groupP    0.574 4.41 29.9   0.130  0.8974

period = outside:
contrast        estimate   SE   df t.ratio p.value
groupN - groupP    8.343 4.41 29.9   1.890  0.0685

period = within:
contrast        estimate   SE   df t.ratio p.value
groupN - groupP   12.677 4.41 29.9   2.872  0.0074

Degrees-of-freedom method: kenward-roger


If you did have an interaction, it would be due to an effect of group only when period = within.

• The anova function in lmerTest constructs the anova-like table, so you can just use anova(lmR_g) for that step. Commented Jan 17, 2023 at 19:52
• This is very helpful - thank you very much for the detailed answer David! Commented Jan 19, 2023 at 17:28