You have been running a factorial (two-way) repeated measures ANOVA with lme() and anova(), and the F-test has revealed a main effect of at least one of the factors, and possibly an interaction. You would like to investigate further to find out which group means are significantly different from each other. To do so, you will need to run the appropriate post-hoc tests to assess the significance of differences between pairs of group means. The function emmeans() will help you do this.
We will reuse the example introduced here (factorial (two-way) repeated measures ANOVA). To create the corresponding dataframe, use the following code:
# response variable
rat.weight <- c(166,168,155,159,151,166,170,160,162,153,220,230,223,233,229,262,274,267,283,274,261,275,264,280,282,343,354,351,359,349,297,311,305,315,303,375,399,388,405,395)
# predictor variables
rat.strain <- as.factor(rep(c(rep("strainA",5),rep("strainB",5)),4))
time.point <- as.factor(c(rep("week08",10), rep("week12",10), rep("week16",10), rep("week20",10)))
# individual ID
rat.ID <- as.factor(rep(c("rat01","rat02","rat03","rat04","rat05","rat06","rat07","rat08","rat09","rat10"),4))
#dataframe
my.dataframe <- data.frame(rat.ID,rat.strain,time.point,rat.weight)For recall, here is the corresponding plot: 
and the model that was fitted with lme() was:
# linear mixed effect model stored in the object 'model'
model <- lme(rat.weight ~ time.point * rat.strain, random=~1|rat.ID, data=my.dataframe)The analysis of variance had shown that there is a main effect of both rat.strain and time.point on rat.weight, as well as an interaction rat.strain:time.point
Running the test with emmeans()
emmeans() is part of the package emmeans, which we first need to activate:
library(emmeans)The next step consists in “feeding” the linear mixed effect model to emmeans(). The syntax is emmeans(model, pairwise ~ predictor1 * predictor2) where predictor1 and predictor2 are the predictor variables which levels have to be compared pairwise:
# running emmeans()
emmeans(model, pairwise ~ rat.strain * time.point)## $emmeans
## rat.strain time.point emmean SE df lower.CL upper.CL
## strainA week08 160 3.5 9 152 168
## strainB week08 162 3.5 8 154 170
## strainA week12 227 3.5 9 219 235
## strainB week12 272 3.5 8 264 280
## strainA week16 272 3.5 9 264 280
## strainB week16 351 3.5 8 343 359
## strainA week20 306 3.5 9 298 314
## strainB week20 392 3.5 8 384 400
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## strainA,week08 - strainB,week08 -2.4 4.96 8 -0.484 0.9994
## strainA,week08 - strainA,week12 -67.2 3.93 24 -17.082 <.0001
## strainA,week08 - strainB,week12 -112.2 4.96 8 -22.636 <.0001
## strainA,week08 - strainA,week16 -112.6 3.93 24 -28.622 <.0001
## strainA,week08 - strainB,week16 -191.4 4.96 8 -38.614 <.0001
## strainA,week08 - strainA,week20 -146.4 3.93 24 -37.214 <.0001
## strainA,week08 - strainB,week20 -232.6 4.96 8 -46.925 <.0001
## strainB,week08 - strainA,week12 -64.8 4.96 8 -13.073 <.0001
## strainB,week08 - strainB,week12 -109.8 3.93 24 -27.910 <.0001
## strainB,week08 - strainA,week16 -110.2 4.96 8 -22.232 <.0001
## strainB,week08 - strainB,week16 -189.0 3.93 24 -48.042 <.0001
## strainB,week08 - strainA,week20 -144.0 4.96 8 -29.051 <.0001
## strainB,week08 - strainB,week20 -230.2 3.93 24 -58.515 <.0001
## strainA,week12 - strainB,week12 -45.0 4.96 8 -9.078 0.0003
## strainA,week12 - strainA,week16 -45.4 3.93 24 -11.540 <.0001
## strainA,week12 - strainB,week16 -124.2 4.96 8 -25.056 <.0001
## strainA,week12 - strainA,week20 -79.2 3.93 24 -20.132 <.0001
## strainA,week12 - strainB,week20 -165.4 4.96 8 -33.368 <.0001
## strainB,week12 - strainA,week16 -0.4 4.96 8 -0.081 1.0000
## strainB,week12 - strainB,week16 -79.2 3.93 24 -20.132 <.0001
## strainB,week12 - strainA,week20 -34.2 4.96 8 -6.900 0.0018
## strainB,week12 - strainB,week20 -120.4 3.93 24 -30.605 <.0001
## strainA,week16 - strainB,week16 -78.8 4.96 8 -15.897 <.0001
## strainA,week16 - strainA,week20 -33.8 3.93 24 -8.592 <.0001
## strainA,week16 - strainB,week20 -120.0 4.96 8 -24.209 <.0001
## strainB,week16 - strainA,week20 45.0 4.96 8 9.078 0.0003
## strainB,week16 - strainB,week20 -41.2 3.93 24 -10.473 <.0001
## strainA,week20 - strainB,week20 -86.2 4.96 8 -17.390 <.0001
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 8 estimatesThe lower part of the table shows all the pairwise comparisons and the corresponding estimates and p-values. This table is rather confusing due to the large amount of information; in addition, not all pairs are meaningful for interpretating the data.
An option to obtain a list of more meaningful pairs is to run the result of emmeans(model, ~ predictor1 * predictor2) through the function pairs() with the argument simple = "each". This creates a table with pairs of means grouped by each single factor. The syntax is pairs(result, simple = "each"):
# running pairs()
emm <- emmeans(model, ~ rat.strain * time.point)
pairs(emm, simple = "each")## $`simple contrasts for rat.strain`
## time.point = week08:
## contrast estimate SE df t.ratio p.value
## strainA - strainB -2.4 4.96 8 -0.484 0.6412
##
## time.point = week12:
## contrast estimate SE df t.ratio p.value
## strainA - strainB -45.0 4.96 8 -9.078 <.0001
##
## time.point = week16:
## contrast estimate SE df t.ratio p.value
## strainA - strainB -78.8 4.96 8 -15.897 <.0001
##
## time.point = week20:
## contrast estimate SE df t.ratio p.value
## strainA - strainB -86.2 4.96 8 -17.390 <.0001
##
## Degrees-of-freedom method: containment
##
## $`simple contrasts for time.point`
## rat.strain = strainA:
## contrast estimate SE df t.ratio p.value
## week08 - week12 -67.2 3.93 24 -17.082 <.0001
## week08 - week16 -112.6 3.93 24 -28.622 <.0001
## week08 - week20 -146.4 3.93 24 -37.214 <.0001
## week12 - week16 -45.4 3.93 24 -11.540 <.0001
## week12 - week20 -79.2 3.93 24 -20.132 <.0001
## week16 - week20 -33.8 3.93 24 -8.592 <.0001
##
## rat.strain = strainB:
## contrast estimate SE df t.ratio p.value
## week08 - week12 -109.8 3.93 24 -27.910 <.0001
## week08 - week16 -189.0 3.93 24 -48.042 <.0001
## week08 - week20 -230.2 3.93 24 -58.515 <.0001
## week12 - week16 -79.2 3.93 24 -20.132 <.0001
## week12 - week20 -120.4 3.93 24 -30.605 <.0001
## week16 - week20 -41.2 3.93 24 -10.473 <.0001
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 4 estimatesIf you are not interested in the full table, but only the part showing the breakdown based on the levels of a specific factor, run pairs(emm, simple = "predictor") instead:
# running pairs()
emm <- emmeans(model, ~ rat.strain * time.point)
pairs(emm, simple = "time.point")## rat.strain = strainA:
## contrast estimate SE df t.ratio p.value
## week08 - week12 -67.2 3.93 24 -17.082 <.0001
## week08 - week16 -112.6 3.93 24 -28.622 <.0001
## week08 - week20 -146.4 3.93 24 -37.214 <.0001
## week12 - week16 -45.4 3.93 24 -11.540 <.0001
## week12 - week20 -79.2 3.93 24 -20.132 <.0001
## week16 - week20 -33.8 3.93 24 -8.592 <.0001
##
## rat.strain = strainB:
## contrast estimate SE df t.ratio p.value
## week08 - week12 -109.8 3.93 24 -27.910 <.0001
## week08 - week16 -189.0 3.93 24 -48.042 <.0001
## week08 - week20 -230.2 3.93 24 -58.515 <.0001
## week12 - week16 -79.2 3.93 24 -20.132 <.0001
## week12 - week20 -120.4 3.93 24 -30.605 <.0001
## week16 - week20 -41.2 3.93 24 -10.473 <.0001
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 4 estimates