Calculate Tukey HSD by RStudio
I want to perform Tukey HSD and other tests using RStudio. In this article, I will describe how to use R for Tukey HSD and other tests as a memo.
Install packages
In this article, we will use multcomp and car for testing and tidyverse for graph creation. To use them as libraries, first download the packages to RStudio with the following code:
install.packages("multcomp")
install.packages("car")
install.packages("tidyverse")
In my environment, I encountered an error downloading tidyverse because the R version of RStudio was old. If tidyverse cannot be downloaded properly, consider updating R.
TukeyHSD
Before performing Tukey HSD, we performed Shapiro's test and Levene's test to confirm whether the data to be used for testing follows a normal distribution and has homoscedasticity, respectively.
Input data format
The data used is in the form of a csv file as shown in the following table. Column names can be anything.
| color | values |
|---|---|
| red | 0.237 |
| red | 0.230 |
| red | 0.156 |
| blue | 0.990 |
| blue | 0.967 |
| blue | 1.146 |
| blue | 1.2005 |
| yellow | 0.15093 |
| yellow | 0.0866 |
| yellow | 0.1337 |
| green | 0.686 |
| green | 0.557 |
| green | 0.486 |
| green | 0.5866 |
Load libraries
library(multcomp)
library(car)
library(tidyverse)
Read the csv file
Read the csv file into R. Also, change the column names of the csv data.
file_path <- "C:/Users/hoge/R/your_data.csv"
data <- read.csv(file_path, header = TRUE, col.names = c("treatment", "values"))
Type conversion
Convert the type of the treatment column to a factor type so that it can be used for testing.
data$treatment <- factor(data$treatment)
Attach to make it easier to specify items
attach(data)
Testing
Normality test
tapply(values, treatment, shapiro.test)
Output
$blue
Shapiro-Wilk normality test
data: X[[i]]
W = 0.87592, p-value = 0.3215
$green
Shapiro-Wilk normality test
data: X[[i]]
W = 0.98214, p-value = 0.9144
$red
Shapiro-Wilk normality test
data: X[[i]]
W = 0.81429, p-value = 0.1491
$yellow
Shapiro-Wilk normality test
data: X[[i]]
W = 0.93295, p-value = 0.4998
Homoscedasticity test
leveneTest(values, treatment, center = mean)
Output
Levene's Test for Homogeneity of Variance (center = mean)
Df F value Pr(>F)
group 3 3.7002 0.05027 .
10
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ANOVA
anova(lm(values ~ treatment))
Output
Analysis of Variance Table
Response: values
Df Sum Sq Mean Sq F value Pr(>F)
treatment 3 2.00564 0.66855 100.44 9.133e-08 ***
Residuals 10 0.06656 0.00666
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey HSD
TukeyHSD(aov(values ~ treatment))
Output
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = values ~ treatment)
$treatment
diff lwr upr p adj
green-blue -0.49697500 -0.6734702 -0.3204798 0.0000304
red-blue -0.86820833 -1.0588449 -0.6775717 0.0000004
yellow-blue -0.95213167 -1.1427683 -0.7614951 0.0000002
red-green -0.37123333 -0.5618699 -0.1805967 0.0006723
yellow-green -0.45515667 -0.6457933 -0.2645201 0.0001270
yellow-red -0.08392333 -0.2877224 0.1198758 0.6061560
Create alphabet from significant difference
Adjust the significance level with level (in this case, 5%). decreasing = T labels the alphabet so that the process with the highest value is a.
result <- aov(values ~ treatment, data = data)
tukey <- glht(result, linfct = mcp(treatment = "Tukey"))
cld(tukey, level = 0.05, decreasing = T)
Output
blue green red yellow
"c" "b" "a" "a"
Creating a bar graph
Let's also create a bar graph.
Store the alphabet created based on Tukey's multiple comparison results in a table data
cld_result <- cld(tukey, decreasing = F)
cld_table <- as.data.frame(cld_result[["mcletters"]][["Letters"]])
colnames(cld_table) <- "alphabet"
cld_letters <- as_tibble(rownames_to_column(cld_table, "treatment"))
Calculate the mean and standard deviation and join with the alphabet table
data <- data %>%
group_by(treatment) %>%
summarise(mean = mean(values), sd = sd(values)) %>%
full_join(cld_letters, by = "treatment")
Create a bar graph
ggplot(data, aes(x = treatment, y = mean)) +
geom_bar(stat = "identity", fill = "steelblue") +
geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +
geom_text(aes(label = alphabet), vjust = -1, hjust = -1, size = 6) +
labs(x = "Treatment", y = "Value", title = "Show Significant Difference")
Output
All code
library(multcomp)
library(car)
library(tidyverse)
# Specify file path
file_path <- "C:/Users/hoge/R/your_data.csv"
# Read CSV file
data <- read.csv(file_path, header = TRUE, col.names = c("treatment", "values"))
data$treatment <- factor(data$treatment)
attach(data)
# Test for normality
tapply(values, treatment, shapiro.test)
# Test for homogeneity of variance
leveneTest(values, treatment, center = mean)
# Perform ANOVA analysis
anova(lm(values ~ treatment))
# Perform Tukey HSD test
TukeyHSD(aov(values ~ treatment))
# Create alphabet based on Tukey HSD results
result <- aov(values ~ treatment, data = data)
tukey <- glht(result, linfct = mcp(treatment = "Tukey"))
cld_result <- cld(tukey, level=0.05,decreasing = T)
cld_table <- as.data.frame(cld_result[["mcletters"]][["Letters"]])
colnames(cld_table) <- "alphabet"
cld_letters <- as_tibble(rownames_to_column(cld_table, "treatment"))
# Combine alphabet table with mean and standard deviation
data <- data %>%
group_by(treatment) %>%
summarise(mean = mean(values), sd = sd(values)) %>%
full_join(cld_letters, by = "treatment")
# Create bar graph
ggplot(data, aes(x = treatment, y = mean)) +
geom_bar(stat = "identity", fill = "steelblue") +
geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +
geom_text(aes(label = alphabet), vjust = -1, hjust = -1, size = 6) +
labs(x = "Treatment", y = "Value", title = "Show Significant Difference")