Anderson-Darling Test for Normality

# Anderson-Darling Test for Normality
# Reaction time data with Q-Q diagnostics

library(nortest)      # ad.test
library(ggplot2)
library(gridExtra)
library(car)          # qqPlot
library(moments)      # skewness, kurtosis
library(boot)

set.seed(42)

# Generate reaction time data (slightly right-skewed)
rt <- c(rnorm(65, 350, 50), rgamma(10, shape = 15, scale = 25))  # Mixture for slight skew

cat("ANDERSON-DARLING TEST FOR NORMALITY\n")
cat("====================================\n\n")

# ====================
# 1. CHECK ASSUMPTIONS
# ====================

cat("ASSUMPTION CHECKS\n")
cat("==================\n\n")

# A1: Continuous data
cat("1. Continuous data:\n")
cat("   Sample size: n =", length(rt), "\n")
cat("   Unique values:", length(unique(rt)), "\n")
cat("   Data type:", class(rt), "\n")
cat("   ✓ Continuous measurement\n\n")

# A2: Independence
cat("2. Independence:\n")
cat("   ✓ Verified by design(separate trials, different stimuli)\n\n")

# A3: Sample size
cat("3. Sample size:\n")
cat("   n =", length(rt), "(minimum n ≥ 7)\n")
cat("   ✓ Adequate sample size\n\n")

# A4: Single population check
cat("4. Single population check:\n")
cat("   Examining histogram for multimodality...\n")
cat("   (Visual inspection in plots)\n\n")

# ====================
# 2. DESCRIPTIVE STATISTICS
# ====================

cat("DESCRIPTIVE STATISTICS\n")
cat("=======================\n\n")

cat("  Mean:", round(mean(rt), 2), "ms\n")
cat("  SD:", round(sd(rt), 2), "ms\n")
cat("  Median:", round(median(rt), 2), "ms\n")
cat("  IQR:", round(IQR(rt), 2), "ms\n")
cat("  Range: [", round(min(rt), 2), ",", round(max(rt), 2), "]\n\n")

cat("  Skewness:", round(skewness(rt), 3), "\n")
cat("  Kurtosis:", round(kurtosis(rt), 3), "(normal = 3)\n")
cat("  Excess kurtosis:", round(kurtosis(rt) - 3, 3), "\n\n")

if (abs(skewness(rt)) > 0.5) {
  cat("  ⚠ Moderate skewness detected\n")
}
if (abs(kurtosis(rt) - 3) > 0.5) {
  cat("  ⚠ Excess kurtosis detected\n")
}
cat("\n")

# ====================
# 3. ANDERSON-DARLING TEST
# ====================

cat("ANDERSON-DARLING TEST\n")
cat("======================\n\n")

ad_result <- ad.test(rt)

cat("Test Results:\n")
cat("  A² statistic:", round(ad_result$statistic, 4), "\n")
cat("  p-value:", round(ad_result$p.value, 4), "\n\n")

alpha <- 0.05
if (ad_result$p.value < alpha) {
  cat("Decision: REJECT H₀ (α = 0.05)\n")
  cat("Interpretation: Data significantly departs from normality\n\n")
} else {
  cat("Decision: FAIL TO REJECT H₀ (α = 0.05)\n")
  cat("Interpretation: Data consistent with normal distribution\n\n")
}

# Interpretation guidelines
cat("Interpretation Guidelines for A²:\n")
cat("  A² < 0.5:  Close to normal\n")
cat("  0.5 ≤ A² < 1.0:  Moderate departure\n")
cat("  A² ≥ 1.0:  Strong departure from normality\n\n")

if (ad_result$statistic < 0.5) {
  cat("  Current A² = ", round(ad_result$statistic, 3), " → Close to normal\n\n")
} else if (ad_result$statistic < 1.0) {
  cat("  Current A² = ", round(ad_result$statistic, 3), " → Moderate departure\n\n")
} else {
  cat("  Current A² = ", round(ad_result$statistic, 3), " → Strong departure\n\n")
}

# ====================
# 4. COMPARISON WITH OTHER TESTS
# ====================

cat("COMPARISON WITH OTHER NORMALITY TESTS\n")
cat("======================================\n\n")

# Shapiro-Wilk
sw_result <- shapiro.test(rt)
cat("Shapiro-Wilk test:\n")
cat("  W statistic:", round(sw_result$statistic, 4), "\n")
cat("  p-value:", round(sw_result$p.value, 4), "\n\n")

# Kolmogorov-Smirnov (Lilliefors)
library(nortest)
lillie_result <- lillie.test(rt)
cat("Lilliefors(KS) test:\n")
cat("  D statistic:", round(lillie_result$statistic, 4), "\n")
cat("  p-value:", round(lillie_result$p.value, 4), "\n\n")

cat("Test Comparison Notes:\n")
cat("  • Shapiro-Wilk: Generally most powerful overall\n")
cat("  • Anderson-Darling: Most sensitive to TAIL deviations\n")
cat("  • Lilliefors/KS: Least powerful; use for general distribution testing\n\n")

# ====================
# 5. BOOTSTRAP CI FOR A²
# ====================

cat("Bootstrap 95% CI for A² Statistic\n")
cat("===================================\n\n")

boot_ad <- function(data, indices) {
  d <- data[indices]
  test <- ad.test(d)
  return(test$statistic)
}

boot_results <- boot(rt, boot_ad, R = 2000)
boot_ci <- boot.ci(boot_results, type = "perc")

cat("  95% CI for A²: [", round(boot_ci$percent[4], 4), ",", 
    round(boot_ci$percent[5], 4), "]\n\n")

# ====================
# 6. VISUALIZATIONS (6 plots)
# ====================

# Plot 1: Q-Q plot with confidence envelope
p1 <- ggplot(data.frame(sample = rt), aes(sample = sample)) +
  stat_qq() +
  stat_qq_line(color = "red", linetype = "dashed", size = 1) +
  labs(title = "Normal Q-Q Plot",
       subtitle = paste0("A² = ", round(ad_result$statistic, 3), 
                         ", p = ", round(ad_result$p.value, 3)),
       x = "Theoretical Quantiles", y = "Sample Quantiles(ms)") +
  theme_minimal()

# Plot 2: Histogram with fitted normal
p2 <- ggplot(data.frame(rt = rt), aes(rt)) +
  geom_histogram(aes(y = after_stat(density)), bins = 20, 
                 fill = "lightblue", color = "black", alpha = 0.7) +
  stat_function(fun = dnorm, args = list(mean = mean(rt), sd = sd(rt)), 
                color = "red", size = 1.2) +
  labs(title = "Histogram with Fitted Normal Density",
       subtitle = paste0("Mean = ", round(mean(rt), 1), ", SD = ", round(sd(rt), 1)),
       x = "Reaction Time(ms)", y = "Density") +
  theme_minimal()

# Plot 3: ECDF vs Normal CDF
p3 <- ggplot(data.frame(rt = rt), aes(rt)) +
  stat_ecdf(geom = "step", color = "blue", size = 1) +
  stat_function(fun = pnorm, args = list(mean = mean(rt), sd = sd(rt)), 
                color = "red", linetype = "dashed", size = 1) +
  labs(title = "ECDF vs Fitted Normal CDF",
       subtitle = "Blue = Empirical, Red = Theoretical",
       x = "Reaction Time(ms)", y = "Cumulative Probability") +
  theme_minimal()

# Plot 4: Detrended Q-Q plot (residuals from Q-Q line)
qq_data <- qqnorm(rt, plot.it = FALSE)
residuals_qq <- qq_data$y - qq_data$x * sd(rt) - mean(rt)

p4 <- ggplot(data.frame(theoretical = qq_data$x, residual = residuals_qq), 
             aes(theoretical, residual)) +
  geom_point(alpha = 0.6, size = 2) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed", size = 1) +
  geom_smooth(se = FALSE, color = "blue", method = "loess") +
  labs(title = "Detrended Q-Q Plot",
       subtitle = "Highlights systematic deviations from normality",
       x = "Theoretical Quantiles", y = "Residuals(ms)") +
  theme_minimal()

# Plot 5: Bootstrap distribution of A²
p5 <- ggplot(data.frame(A2 = boot_results$t), aes(A2)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "black", alpha = 0.7) +
  geom_vline(xintercept = ad_result$statistic, color = "red", 
             linetype = "dashed", size = 1.2) +
  geom_vline(xintercept = boot_ci$percent[4:5], color = "orange", 
             linetype = "dotted", size = 1) +
  labs(title = "Bootstrap Distribution of A²",
       subtitle = paste0("Observed A² = ", round(ad_result$statistic, 3)),
       x = "A² Statistic", y = "Frequency") +
  theme_minimal()

# Plot 6: Test comparison
test_results <- data.frame(
  Test = c("Anderson-Darling", "Shapiro-Wilk", "Lilliefors"),
  p_value = c(ad_result$p.value, sw_result$p.value, lillie_result$p.value),
  statistic = c(ad_result$statistic, sw_result$statistic, lillie_result$statistic)
)

p6 <- ggplot(test_results, aes(x = Test, y = p_value, fill = Test)) +
  geom_col(alpha = 0.7, show.legend = FALSE) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "red", size = 1) +
  geom_text(aes(label = round(p_value, 3)), vjust = -0.5) +
  labs(title = "Normality Test Comparison",
       subtitle = "Red line = α = 0.05",
       y = "p-value") +
  ylim(0, max(test_results$p_value) * 1.2) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)

# ====================
# 7. INTERPRETATION
# ====================

cat("\nDETAILED INTERPRETATION\n")
cat("=======================\n\n")

cat("The Anderson-Darling test assessed normality of reaction times.\n\n")

cat("Test Results:\n")
cat("  • A² = ", round(ad_result$statistic, 3), "\n")
cat("  • p-value = ", round(ad_result$p.value, 3), "\n")
cat("  • Skewness = ", round(skewness(rt), 3), "\n")
cat("  • Excess kurtosis = ", round(kurtosis(rt) - 3, 3), "\n\n")

if (ad_result$p.value < 0.05) {
  cat("Conclusion: Data significantly deviates from normality.\n\n")
  cat("Recommendations:\n")
  cat("  1. Examine Q-Q plot for nature of departure(tails, skewness)\n")
  cat("  2. Consider transformations: log(right skew), sqrt, inverse\n")
  cat("  3. Use nonparametric alternatives(Mann-Whitney, Kruskal-Wallis)\n")
  cat("  4. Use robust methods(trimmed means, bootstrap)\n")
else {
  cat("Conclusion: Data consistent with normal distribution.\n\n")
  cat("Recommendations:\n")
  cat("  1. Parametric methods appropriate(t-test, ANOVA, regression)\n")
  cat("  2. Still check Q-Q plot for minor deviations\n")
  cat("  3. Consider sample size when interpreting results\n")
}

cat("\nWhy Anderson-Darling?\n")
cat("  • More powerful than KS test(especially for tail deviations)\n")
cat("  • Weights tails more heavily(important for outliers)\n")
cat("  • Adjusts for parameter estimation(unlike standard KS)\n")
cat("  • Comparable power to Shapiro-Wilk\n")