You have been running a repeated measures ANOVA with lme() and anova(), and the F-test has revealed the existence of a significant difference between some of the tested groups. But which groups? To answer that question, you will need to run the appropriate post-hoc tests to assess the significance of differences between pairs of group means. The functions emmeans() and glht() will help you do this.

We will reuse the example introduced here (repeated measures ANOVA). To create the corresponding dataframe, use the following code:

# response variable
rat.weight <- c(164,164,158,159,155,220,230,226,227,222,261,275,264,280,272,306,326,320,330,312)

# predictor variable
time.point <- as.factor(c(rep("week08",5), rep("week12",5), rep("week16",5), rep("week20",5)))

# individual ID
rat.ID <- as.factor(rep(c("rat1","rat2","rat3", "rat4", "rat5"),4))

# dataframe
my.dataframe <- data.frame(rat.ID,time.point,rat.weight)



For recall, here is the corresponding plot:

and the model that was fitted with lme():

# linear mixed effect model stored in the object 'model'
model <- lme(rat.weight ~ time.point, random=~1|rat.ID, data=my.dataframe)



Alternative 1: 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 ~ predictor) where predictor is the predictor variable which levels have to be compared pairwise:

# running emmeans()
emmeans(model, pairwise ~ time.point)
## $emmeans
##  time.point emmean   SE df lower.CL upper.CL
##  week08        160 3.08  4      151      169
##  week12        225 3.08  4      216      234
##  week16        270 3.08  4      262      279
##  week20        319 3.08  4      310      327
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast        estimate   SE df t.ratio p.value
##  week08 - week12    -65.0 3.39 12 -19.159 <.0001 
##  week08 - week16   -110.4 3.39 12 -32.541 <.0001 
##  week08 - week20   -158.8 3.39 12 -46.807 <.0001 
##  week12 - week16    -45.4 3.39 12 -13.382 <.0001 
##  week12 - week20    -93.8 3.39 12 -27.648 <.0001 
##  week16 - week20    -48.4 3.39 12 -14.266 <.0001 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 4 estimates

The output shows first some statistics for each group in the upper part, while the lower part displays the six pairwise comparisons and p-values. Here, all comparisons are associated with a p-value under 0.05.

Alternative 2: running the test with glht()

glht() is part of the package multcomp which we first need to activate:

library(multcomp)


Running the test with glht() involves two steps:

Let’s apply this to our example:

# running glht()
post.hoc <- glht(model, linfct = mcp(time.point = 'Tukey'))

# displaying the result table with summary()
summary(post.hoc)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = rat.weight ~ time.point, data = my.dataframe, 
##     random = ~1 | rat.ID)
## 
## Linear Hypotheses:
##                      Estimate Std. Error z value Pr(>|z|)    
## week12 - week08 == 0   65.000      3.393   19.16   <2e-16 ***
## week16 - week08 == 0  110.400      3.393   32.54   <2e-16 ***
## week20 - week08 == 0  158.800      3.393   46.81   <2e-16 ***
## week16 - week12 == 0   45.400      3.393   13.38   <2e-16 ***
## week20 - week12 == 0   93.800      3.393   27.65   <2e-16 ***
## week20 - week16 == 0   48.400      3.393   14.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

The output shows all six pairwise comparisons, p-values and corresponding stars, and all pairwise comparisons are accompanied by a p-value less than 0.05.