Solved – Post Hoc Pairwise Comparison of Interaction in Mixed Effects (lmer) Model

Question

I have searched other posts and textbooks and found numerous variations of the pairwise comparisions for my mixed effects model using code from the multicomp package as such…How to perform post-hoc test on lmer model?

lsmeansLT(MODEL, pairwise~GROUP, adjust="tukey")

or

summary(glht(MODEL, linfct=mcp(GROUP = "Tukey")))

My data includes four treatments including Burned North, Unburned North, Burned South, and Unburned South facing aspects with snow depth as my response variable. Burned Condition and Aspect are both fixed effects in my model

Depth1 = lmer(Depth~Burn_Con*Aspect+Year+(1|Block), data = snow)

I can compare total burned vs unburned but cannot seem to find how to compare the interaction of burned north and unburned north facing aspects. Trying code like this spits out a number of errors…

summary(glht(Depth1, linfct=mcp(Burn_Con*Aspect = "Tukey")))

Any help would be very much appreciated!

Answer 1

Try the emmeans package. Something like

library(emmeans)
emm = emmeans(Depth1, ~ Burn_Con * Aspect)
pairs(emm)
# or for simple comparisons
pairs(emm, simple = “each”)