Tukey HSD

# Tukey HSD: All pairwise diet comparisons following ANOVA
# Example 1: 4 diet groups, 6 pairwise comparisons

# Load packages
library(car)           # Levene's test
library(effectsize)    # Effect sizes
library(tidyverse)     # Data manipulation
library(multcomp)      # Advanced contrasts
library(emmeans)       # Estimated marginal means

# Simulate realistic data (or load: data <- read.csv("diet_study.csv"))
set.seed(2025)
data <- data.frame(
  diet = rep(c("Low-carb", "Mediterranean", "Low-fat", "Control"), each=28),
  weight_loss = c(
    rnorm(28, mean=7.2, sd=3.1),   # Low-carb: M=7.2 kg
    rnorm(28, mean=6.5, sd=2.8),   # Mediterranean: M=6.5 kg
    rnorm(28, mean=4.3, sd=2.9),   # Low-fat: M=4.3 kg
    rnorm(28, mean=2.1, sd=2.5)    # Control: M=2.1 kg
  )
)

cat("=== TUKEY HSD POST-HOC ANALYSIS ===")
cat("\nAfter 12-week diet intervention: Which diets differ?\n\n")

# === STEP 1: Check Assumptions ===

cat("--- STEP 1: Assumption Checks ---\n")

# 1. Normality per group
cat("\n1. Normality(Shapiro-Wilk per group):\n")
by(data$weight_loss, data$diet, shapiro.test)
# All p > .05 → normality OK

# Q-Q plots
par(mfrow=c(2,2))
for (diet_type in c("Low-carb", "Mediterranean", "Low-fat", "Control")) {
  qqnorm(data$weight_loss[data$diet == diet_type], main=diet_type)
  qqline(data$weight_loss[data$diet == diet_type], col="red")
}
par(mfrow=c(1,1))

# 2. Homogeneity of variance (Levene's test - CRITICAL for Tukey HSD)
cat("\n2. Homogeneity of Variance(Levene's test):\n")
levene_result <- leveneTest(weight_loss ~ diet, data=data)
print(levene_result)
if (levene_result$`Pr(>F)`[1] < 0.05) {
  cat("\n*** WARNING: Unequal variances detected! Use Games-Howell instead. ***\n")
} else {
  cat("\n✓ Equal variances assumption met(p > .05). Tukey HSD is appropriate.\n")
}

# 3. Outliers
cat("\n3. Outliers Check:\n")
ggplot(data, aes(x=diet, y=weight_loss, fill=diet)) +
  geom_boxplot(alpha=0.7) +
  labs(title="Weight Loss by Diet Type\n(Check for outliers)",
       x="Diet Group", y="Weight Loss(kg)") +
  theme_classic() +
  theme(legend.position="none")

# === STEP 2: Descriptive Statistics ===

cat("\n--- STEP 2: Descriptive Statistics ---\n")
desc_stats <- data %>%
  group_by(diet) %>%
  summarise(n = n(),
            M = mean(weight_loss),
            SD = sd(weight_loss),
            SE = SD/sqrt(n),
            CI_lower = M - 1.96*SE,
            CI_upper = M + 1.96*SE)
print(desc_stats)

# === STEP 3: Run Omnibus ANOVA (prerequisite) ===

cat("\n--- STEP 3: Omnibus One-Way ANOVA ---\n")
anova_model <- aov(weight_loss ~ diet, data=data)
anova_summary <- summary(anova_model)
print(anova_summary)

# Extract F-statistic and p-value
f_stat <- anova_summary[[1]]$`F value`[1]
p_value <- anova_summary[[1]]$`Pr(>F)`[1]

if (p_value >= 0.05) {
  cat("\n*** STOP: ANOVA non-significant(p ≥ .05). Do NOT proceed to Tukey HSD. ***\n")
  cat("Conducting post-hoc tests after non-significant ANOVA inflates Type I error.\n")
  stop("ANOVA must be significant before post-hoc testing.")
} else {
  cat("\n✓ ANOVA is significant(p < .05). Proceed to Tukey HSD post-hoc tests.\n")
}

# Effect size
omega_sq <- omega_squared(anova_model)
cat("\nEffect size: ω² =", round(omega_sq$Omega2, 3))
if (omega_sq$Omega2 < 0.06) {
  cat(" (small)")
} else if (omega_sq$Omega2 < 0.14) {
  cat(" (medium)")
} else {
  cat(" (large)")
}
cat("\n")

# === STEP 4: Tukey HSD Post-hoc Test ===

cat("\n--- STEP 4: Tukey HSD Pairwise Comparisons ---\n")
cat("Number of pairwise comparisons: k(k-1)/2 = 4×3/2 = 6\n")
cat("Family-wise error rate(FWER) controlled at α = .05\n\n")

tukey_result <- TukeyHSD(anova_model, conf.level=0.95)
print(tukey_result)

# Extract pairwise comparisons
tukey_df <- as.data.frame(tukey_result$diet)
tukey_df$comparison <- rownames(tukey_df)
tukey_df <- tukey_df %>%
  mutate(sig = ifelse(`p adj` < 0.001, "***",
                      ifelse(`p adj` < 0.01, "**",
                             ifelse(`p adj` < 0.05, "*", "ns"))),
         cohen_d = diff / sqrt(sum(anova_model$residuals^2) / anova_model$df.residual))

cat("\nPairwise Comparisons Summary:\n")
print(tukey_df[, c("comparison", "diff", "lwr", "upr", "p adj", "sig")])

# === STEP 5: Effect Sizes for Pairwise Comparisons ===

cat("\n--- STEP 5: Cohen's d for Each Pairwise Comparison ---\n")

# Calculate pooled SD from ANOVA
MSE <- sum(anova_model$residuals^2) / anova_model$df.residual
pooled_sd <- sqrt(MSE)

# Cohen's d for each comparison
tukey_df <- tukey_df %>%
  mutate(cohen_d = diff / pooled_sd,
         d_interpretation = case_when(
           abs(cohen_d) < 0.2 ~ "negligible",
           abs(cohen_d) < 0.5 ~ "small",
           abs(cohen_d) < 0.8 ~ "medium",
           TRUE ~ "large"
         ))

cat("\nEffect Sizes:\n")
print(tukey_df[, c("comparison", "diff", "cohen_d", "d_interpretation", "p adj")])

# === STEP 6: Compact Letter Display ===

cat("\n--- STEP 6: Compact Letter Display(Homogeneous Groups) ---\n")
cat("Groups sharing a letter are NOT significantly different.\n\n")

library(multcomp)
cld_result <- cld(glht(anova_model, linfct=mcp(diet="Tukey")))
print(cld_result)

# Add letters to descriptive stats
cld_letters <- data.frame(
  diet = names(cld_result$mcletters$Letters),
  letter = cld_result$mcletters$Letters
)
desc_stats_cld <- left_join(desc_stats, cld_letters, by="diet")
print(desc_stats_cld)

# === STEP 7: Bootstrap Confidence Intervals ===

cat("\n--- STEP 7: Bootstrap Confidence Intervals(1500 iterations) ---\n")

set.seed(2025)
n_boot <- 1500
boot_diffs <- matrix(NA, nrow=n_boot, ncol=6)
comparisons <- combn(unique(data$diet), 2, simplify=FALSE)

for (i in 1:n_boot) {
  # Resample within each group
  boot_sample <- data %>%
    group_by(diet) %>%
    slice_sample(n=n(), replace=TRUE) %>%
    ungroup()
  
  # Calculate pairwise differences
  means_boot <- boot_sample %>%
    group_by(diet) %>%
    summarise(M = mean(weight_loss), .groups="drop")
  
  for (j in 1:6) {
    g1 <- comparisons[[j]][1]
    g2 <- comparisons[[j]][2]
    m1 <- means_boot$M[means_boot$diet == g1]
    m2 <- means_boot$M[means_boot$diet == g2]
    boot_diffs[i, j] <- m1 - m2
  }
}

# Calculate bootstrap CIs
boot_cis <- apply(boot_diffs, 2, quantile, probs=c(0.025, 0.975))
colnames(boot_cis) <- sapply(comparisons, function(x) paste(x, collapse=" - "))

cat("\nBootstrap 95% CIs(compare to Tukey CIs):\n")
print(round(boot_cis, 2))

# === STEP 8: Visualizations ===

cat("\n--- STEP 8: Visualizations ---\n")

# 8a. Bar plot with compact letters
ggplot(desc_stats_cld, aes(x=reorder(diet, -M), y=M, fill=diet)) +
  geom_bar(stat="identity", width=0.6, alpha=0.8) +
  geom_errorbar(aes(ymin=CI_lower, ymax=CI_upper), width=0.2) +
  geom_text(aes(label=letter), vjust=-0.5, size=6, fontface="bold") +
  labs(title="Weight Loss by Diet Type(12-week intervention)",
       subtitle="Bars with different letters differ significantly(Tukey HSD, α=.05)",
       x="Diet Group", y="Mean Weight Loss(kg) ± 95% CI") +
  scale_fill_brewer(palette="Set2") +
  theme_classic() +
  theme(legend.position="none",
        plot.title = element_text(face="bold", size=14))

# 8b. Forest plot of pairwise differences
tukey_df$comparison_clean <- gsub("-", " vs ", tukey_df$comparison)

ggplot(tukey_df, aes(x=diff, y=reorder(comparison_clean, diff))) +
  geom_vline(xintercept=0, linetype="dashed", color="gray50") +
  geom_errorbarh(aes(xmin=lwr, xmax=upr), height=0.2, size=1) +
  geom_point(aes(color=sig), size=4) +
  scale_color_manual(values=c("***"="red", "**"="orange", "*"="yellow", "ns"="gray"),
                     name="Significance") +
  labs(title="Forest Plot: All Pairwise Comparisons(Tukey HSD)",
       subtitle="Error bars = 95% confidence intervals",
       x="Mean Difference in Weight Loss(kg)",
       y="Pairwise Comparison") +
  theme_minimal() +
  theme(plot.title = element_text(face="bold"))

# 8c. Interaction plot
ggplot(desc_stats, aes(x=diet, y=M, group=1)) +
  geom_line(size=1.2, color="steelblue") +
  geom_point(size=4, color="darkblue") +
  geom_errorbar(aes(ymin=CI_lower, ymax=CI_upper), width=0.2) +
  labs(title="Group Means with 95% Confidence Intervals",
       x="Diet Type", y="Mean Weight Loss(kg)") +
  theme_classic()

# === STEP 9: APA-Style Reporting ===

cat("\n--- STEP 9: APA-Style Results ---\n\n")
cat(paste0(
  "A one-way ANOVA was conducted to compare weight loss across four diet interventions.\n",
  "Data met assumptions of normality(Shapiro-Wilk p > .05 for all groups) and homogeneity\n",
  "of variance(Levene's test, F(", anova_summary[[1]]$Df[1], ", ", anova_summary[[1]]$Df[2], ") = ",
  round(levene_result$`F value`[1], 2), ", p = ", round(levene_result$`Pr(>F)`[1], 3), ").\n\n",
  "There was a significant effect of diet type on weight loss, F(", anova_summary[[1]]$Df[1], ", ",
  anova_summary[[1]]$Df[2], ") = ", round(f_stat, 2), ", p < .001, ω² = ",
  round(omega_sq$Omega2, 2), " (large effect).\n\n",
  "Post-hoc comparisons using Tukey HSD indicated:\n",
  "• Low-carb(M = ", round(desc_stats$M[desc_stats$diet=="Low-carb"], 1), ", SD = ",
  round(desc_stats$SD[desc_stats$diet=="Low-carb"], 1), ") produced significantly greater weight loss than\n",
  "  Low-fat(M = ", round(desc_stats$M[desc_stats$diet=="Low-fat"], 1), ", SD = ",
  round(desc_stats$SD[desc_stats$diet=="Low-fat"], 1), "), p < .001, d = ",
  round(abs(tukey_df$cohen_d[grepl("Low-carb.*Low-fat", tukey_df$comparison)]), 2), ",\n",
  "  and Control(M = ", round(desc_stats$M[desc_stats$diet=="Control"], 1), ", SD = ",
  round(desc_stats$SD[desc_stats$diet=="Control"], 1), "), p < .001, d = ",
  round(abs(tukey_df$cohen_d[grepl("Low-carb.*Control", tukey_df$comparison)]), 2), ".\n",
  "• Mediterranean diet showed similar patterns.\n",
  "• No significant difference between Low-carb and Mediterranean diets(p = ",
  round(tukey_df$`p adj`[grepl("Mediterranean.*Low-carb", tukey_df$comparison)], 3), ").\n\n",
  "These findings suggest low-carb and Mediterranean diets are superior to low-fat and\n",
  "control conditions for weight loss, with large effect sizes(d > 0.8)."
))

cat("\n\n=== ANALYSIS COMPLETE ===")