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
Best Answer
First, it's worth considering whether your interaction truly is significant. Looking at the
summary()of your model: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:
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:If you did have an interaction, it would be due to an effect of
grouponly whenperiod = within.