# Load packages
library(dplyr)
library(ggplot2)
library(agricolae)

# Use ChickWeight dataset
data("ChickWeight")

# 1. Fit ANOVA
model <- aov(weight ~ Diet, data = ChickWeight)
summary(model)

# 2. Tukey HSD with post-hoc letter generation
HSD <- HSD.test(model, "Diet", group = TRUE)
HSD$groups     # displays letters

# 3. Create summary table (mean + SE per Diet)
summary_data <- ChickWeight %>%
  group_by(Diet) %>%
  summarise(
    mean_weight = mean(weight),
    se = sd(weight) / sqrt(n())
  ) %>%
  ungroup()

# Add letters to summary_data
summary_data$letters <- HSD$groups[match(summary_data$Diet,
                                         rownames(HSD$groups)), "groups"]

summary_data

# 4. Plot bar chart with error bars + post-hoc letters
ggplot(summary_data, aes(x = Diet, y = mean_weight, fill = Diet)) +
  geom_col(width = 0.7) +
  geom_errorbar(aes(ymin = mean_weight - se,
                    ymax = mean_weight + se),
                width = 0.15) +
  geom_text(aes(label = letters,
                y = mean_weight + se + 5),  # push letters above bars
            vjust = 0,
            size = 5) +
  labs(x = "Diet",
       y = "Mean Weight (g)") +
  theme_classic() +
  theme(legend.position = "none") +
  expand_limits(y = max(summary_data$mean_weight + summary_data$se) * 1.20)