R-08: Post-hoc test: Tukey HSD

Post hoc tests are designed for situations in which we have already obtained a significant omnibus F-test with a factor that consists of three or more means and additional exploration of the differences among means is needed to provide specific information on which means are significantly different from each other. In our example, we saw that Sepal.Length was significantly different between the 3 species. Now, we want to find out which species had significantly longer Sepal.Length.

Load packages

library(tidyverse)
library(agricolae)

    ## Warning: package 'agricolae' was built under R version 3.3.3

Read input data

irisdata<-read_csv("irisdata.csv")

    ## Parsed with column specification:
    ## cols(
    ##   Sepal.Length = col_double(),
    ##   Sepal.Width = col_double(),
    ##   Petal.Length = col_double(),
    ##   Petal.Width = col_double(),
    ##   Species = col_character()
    ## )

Convert character into factors

Since Species is not a factor, we will convert it into a factor:

irisdata$Species <- as.factor(irisdata$Species)

Fit a model and view residuals

model1 <- aov(Sepal.Length ~ Species, data=irisdata)
par(mfrow=c(2,2))
plot(model1, pch=19, col="blue")

Get summary of model

summary(model1)

  ##              Df Sum Sq Mean Sq F value Pr(>F)    
  ## Species       2  63.21  31.606   119.3 <2e-16 ***
  ## Residuals   147  38.96   0.265                   
  ## ---
  ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Apply post-hoc test

hsd_result<-HSD.test(model1, "Species", group=TRUE)
hsd_result

## $statistics
##     MSerror  Df     Mean       CV       MSD
##   0.2650082 147 5.843333 8.809859 0.2437727
## 
## $parameters
##    test  name.t ntr StudentizedRange alpha
##   Tukey Species   3         3.348424  0.05
## 
## $means
##            Sepal.Length       std  r Min Max   Q25 Q50 Q75
## setosa            5.006 0.3524897 50 4.3 5.8 4.800 5.0 5.2
## versicolor        5.936 0.5161711 50 4.9 7.0 5.600 5.9 6.3
## virginica         6.588 0.6358796 50 4.9 7.9 6.225 6.5 6.9
## 
## $comparison
## NULL
## 
## $groups
##            Sepal.Length groups
## virginica         6.588      a
## versicolor        5.936      b
## setosa            5.006      c
## 
## attr(,"class")
## [1] "group"

This shows that the Sepal.Length of all 3 species is significantly different than one another.

<– Click here to go to the previous tutorial                         Click here to go to the next tutorial –>