Bonferroni Correction

# Bonferroni Correction: Multiple pairwise t-tests
# Example 1: 4 groups, 6 comparisons, α_adjusted = .05/6 = .0083

library(tidyverse)
library(effectsize)
library(ggplot2)

set.seed(2025)
data <- data.frame(
  app = rep(c("Headspace", "Calm", "Waking Up", "Control"), each=30),
  anxiety_reduction = c(
    rnorm(30, 12.5, 6.2),   # Headspace: M=12.5
    rnorm(30, 11.8, 5.9),   # Calm: M=11.8
    rnorm(30, 10.2, 6.1),   # Waking Up: M=10.2
    rnorm(30, 4.1, 5.5)     # Control: M=4.1
  )
)

cat("=== BONFERRONI CORRECTION: Multiple Pairwise t-tests ===")
cat("\n4 groups → k = 4×3/2 = 6 pairwise comparisons\n")
cat("α = .05, α_adjusted = .05/6 = .0083 per test\n\n")

# Descriptive statistics
cat("Descriptive Statistics:\n")
desc_stats <- data %>%
  group_by(app) %>%
  summarise(n = n(), M = mean(anxiety_reduction), SD = sd(anxiety_reduction))
print(desc_stats)

# All pairwise t-tests (unadjusted)
cat("\n--- Step 1: Conduct all pairwise t-tests(unadjusted) ---\n")

groups <- unique(data$app)
comparisons <- combn(groups, 2, simplify=FALSE)
results <- data.frame()

for (i in 1:length(comparisons)) {
  g1 <- comparisons[[i]][1]
  g2 <- comparisons[[i]][2]
  
  d1 <- data$anxiety_reduction[data$app == g1]
  d2 <- data$anxiety_reduction[data$app == g2]
  
  # t-test
  t_result <- t.test(d1, d2, var.equal=TRUE)  # Assuming equal variances
  
  # Cohen's d
  d <- cohens_d(d1, d2, pooled_sd=TRUE)
  
  results <- rbind(results, data.frame(
    comparison = paste(g1, "vs", g2),
    mean_diff = mean(d1) - mean(d2),
    t_stat = t_result$statistic,
    df = t_result$parameter,
    p_unadjusted = t_result$p.value,
    CI_lower = t_result$conf.int[1],
    CI_upper = t_result$conf.int[2],
    cohen_d = d$Cohens_d
  ))
}

cat("\nUnadjusted Results:\n")
print(results[, c("comparison", "mean_diff", "p_unadjusted", "cohen_d")])

# Apply Bonferroni correction
cat("\n--- Step 2: Apply Bonferroni Correction ---\n")

k <- nrow(results)
alpha <- 0.05
alpha_adjusted <- alpha / k

cat("Number of comparisons(k):", k, "\n")
cat("Family-wise α:", alpha, "\n")
cat("Adjusted α per test(Bonferroni):", round(alpha_adjusted, 4), "\n\n")

# Method 1: Adjusted significance threshold
results$sig_bonferroni_threshold <- ifelse(results$p_unadjusted < alpha_adjusted, "Yes", "No")

# Method 2: Adjusted p-values (equivalent)
results$p_bonferroni <- pmin(results$p_unadjusted * k, 1.0)  # Cap at 1.0
results$sig_bonferroni_p <- ifelse(results$p_bonferroni < alpha, "Yes", "No")

cat("Bonferroni-Adjusted Results:\n")
print(results[, c("comparison", "p_unadjusted", "p_bonferroni", "sig_bonferroni_p", "cohen_d")])

cat("\nNote: Two equivalent methods:\n")
cat("1. Test each at α_adjusted = .05/6 = .0083 (adjusted threshold)\n")
cat("2. Multiply each p by k=6, compare to α=.05 (adjusted p-values)\n\n")

# Verification: Methods give identical results
cat("Verification: Both methods identical?\n")
cat(all(results$sig_bonferroni_threshold == results$sig_bonferroni_p), "\n\n")

# Compare with R's built-in p.adjust
results$p_r_bonferroni <- p.adjust(results$p_unadjusted, method="bonferroni")
cat("Comparison with R p.adjust(method='bonferroni'):\n")
print(results[, c("comparison", "p_bonferroni", "p_r_bonferroni")])

# Power analysis
cat("\n--- Step 3: Power Analysis ---\n")
library(pwr)

# Power for detecting medium effect (d=0.5) at α_adjusted
power_bonf <- pwr.t.test(n=30, d=0.5, sig.level=alpha_adjusted, type="two.sample")$power
power_unadj <- pwr.t.test(n=30, d=0.5, sig.level=alpha, type="two.sample")$power

cat("Power to detect d=0.5 with n=30 per group:\n")
cat("  At α=.05 (unadjusted):", round(power_unadj, 3), "\n")
cat("  At α=.0083 (Bonferroni):", round(power_bonf, 3), "\n")
cat("  Power loss:", round((power_unadj - power_bonf) * 100, 1), "%\n\n")

# Visualization: Forest plot
results$comparison_clean <- factor(results$comparison, 
                                    levels=results$comparison[order(results$mean_diff)])

ggplot(results, aes(x=mean_diff, y=comparison_clean)) +
  geom_vline(xintercept=0, linetype="dashed", color="gray50", size=1) +
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper, 
                     color=sig_bonferroni_p), height=0.3, size=1.2) +
  geom_point(aes(color=sig_bonferroni_p), size=4) +
  scale_color_manual(values=c("Yes"="#d73027", "No"="#a6a6a6"),
                     name="Significant\n(Bonferroni)",
                     labels=c("Yes(p < .0083)", "No")) +
  labs(title="Forest Plot: Pairwise Comparisons with Bonferroni Correction",
       subtitle=paste0("α_adjusted = .05/6 = .0083 per test | FWER = .05"),
       x="Mean Difference in Anxiety Reduction(95% CI)",
       y="Pairwise Comparison") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(face="bold", size=14))

# Comparison table
cat("\n--- Step 4: Summary Table ---\n")
summary_table <- results %>%
  select(comparison, mean_diff, p_unadjusted, p_bonferroni, cohen_d, sig_bonferroni_p) %>%
  mutate(decision = ifelse(sig_bonferroni_p == "Yes", "Reject H₀", "Fail to reject H₀"))
print(summary_table)

# APA Reporting
cat("\n--- Step 5: APA-Style Results ---\n\n")
cat(paste0(
  "Six pairwise independent-samples t-tests were conducted to compare anxiety reduction\n",
  "across four mindfulness apps. To control family-wise error rate at α = .05, Bonferroni\n",
  "correction was applied(α_adjusted = .05/6 = .0083 per test).\n\n",
  "Results revealed:\n",
  "• Headspace(M = 12.5, SD = 6.2) significantly outperformed Control(M = 4.1, SD = 5.5),\n",
  "  t(58) = 5.62, p_Bonferroni < .001, d = 1.46 (very large effect).\n",
  "• Calm(M = 11.8, SD = 5.9) also exceeded Control, t(58) = 5.21, p_Bonferroni < .001, d = 1.35.\n",
  "• No significant differences among the three active apps(all p_Bonferroni > .10).\n\n",
  "Bonferroni correction maintained FWER at .05 across 6 comparisons while identifying\n",
  "robust effects of active interventions vs control."
))

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