clean

rm(list = ls()); graphics.off()

Set R packages path and load packages

.libPaths("C:/Rpack")

library(openxlsx)
library(ggpubr)
## Loading required package: ggplot2
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(multcompView)

Import data

dtaAOV <- read.xlsx("dtaHedonic.xlsx", sheet = 1)

ANOVA model

Set factors

dtaAOV$assessor <- as.factor(dtaAOV$assessor)
dtaAOV$product <- as.factor(dtaAOV$product)

levels(dtaAOV$product)[1:7] <- paste0("P",1:7)

Run ANOVA

res.aov2 <- aov(liking~assessor+product, data = dtaAOV)
summary(res.aov2)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## assessor     49  302.8    6.18   1.952 0.000396 ***
## product       6  702.2  117.04  36.963  < 2e-16 ***
## Residuals   294  930.9    3.17                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post-hoc test

TukeyHSD(res.aov2, which = "product")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = liking ~ assessor + product, data = dtaAOV)
## 
## $product
##        diff        lwr         upr     p adj
## P2-P1 -1.00 -2.0565836  0.05658362 0.0769798
## P3-P1 -4.70 -5.7565836 -3.64341638 0.0000000
## P4-P1 -1.48 -2.5365836 -0.42341638 0.0008235
## P5-P1 -0.70 -1.7565836  0.35658362 0.4380568
## P6-P1 -1.18 -2.2365836 -0.12341638 0.0176121
## P7-P1 -2.32 -3.3765836 -1.26341638 0.0000000
## P3-P2 -3.70 -4.7565836 -2.64341638 0.0000000
## P4-P2 -0.48 -1.5365836  0.57658362 0.8281611
## P5-P2  0.30 -0.7565836  1.35658362 0.9801587
## P6-P2 -0.18 -1.2365836  0.87658362 0.9987655
## P7-P2 -1.32 -2.3765836 -0.26341638 0.0046044
## P4-P3  3.22  2.1634164  4.27658362 0.0000000
## P5-P3  4.00  2.9434164  5.05658362 0.0000000
## P6-P3  3.52  2.4634164  4.57658362 0.0000000
## P7-P3  2.38  1.3234164  3.43658362 0.0000000
## P5-P4  0.78 -0.2765836  1.83658362 0.3031554
## P6-P4  0.30 -0.7565836  1.35658362 0.9801587
## P7-P4 -0.84 -1.8965836  0.21658362 0.2194678
## P6-P5 -0.48 -1.5365836  0.57658362 0.8281611
## P7-P5 -1.62 -2.6765836 -0.56341638 0.0001570
## P7-P6 -1.14 -2.1965836 -0.08341638 0.0250691
tukey.plot <- TukeyHSD(res.aov2, which = "product")
windows(10,10)
plot(tukey.plot, las = 1, cex.axis = 0.7)

## Mean plot

windows(10,10)
ggbarplot(dtaAOV, x = "product", y = "liking",
          add = "mean_se", label = TRUE, lab.vjust = -2,
          lab.nb.digits = 1)

Create letter display

data frame with factors, means, se

std.error <- function(x) sd(x)/sqrt(length(x))

data_summary <- group_by(dtaAOV, product) %>%
  summarise(mean = mean(liking), sd = sd(liking), se = std.error(liking)) %>%
  arrange(desc(mean))

create letters

res.aov1 <- aov(liking~product, data = dtaAOV)
tukey <- TukeyHSD(res.aov1)
tukey.cld <- multcompLetters4(res.aov1,tukey)
print(tukey.cld)
## $product
##   P1   P5   P2   P6   P4   P7   P3 
##  "a" "ab" "ab"  "b" "bc"  "c"  "d"

create data frame with means, letters

cld <- as.data.frame.list(tukey.cld$'product')
data_summary$Tukey <- cld$Letters

plot

windows(10,10)
ggplot(data_summary, aes(x=product, y = mean)) + 
  geom_bar(stat = "identity", position = "dodge", alpha = 0.5,
           fill = "lightblue", color = "darkblue") +
  geom_errorbar(aes(ymin = mean-se, ymax = mean+se),
                position = position_dodge(0.9), width = 0.1,
                color="darkblue")+
  labs(x="Product", y="Liking") +
  theme_classic() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  # theme(legend.position = c(0.1, 0.75)) +
  geom_text(aes(label=Tukey), position = position_dodge(0.9), size=3,
            vjust=-2.5, hjust=-0.5, color="darkblue")