Shapiro-Wilk Normality Test

# Shapiro-Wilk Normality Test: Exam Scores Example
# ================================================================
# Research Question: Are exam scores normally distributed?

library(tidyverse)
library(moments)  # skewness(), kurtosis()

# Create exam score data (n=35)
set.seed(101)
exam_scores <- c(
  rnorm(30, mean = 75, sd = 12),  # Most students
  rnorm(5, mean = 90, sd = 5)     # Few high achievers
)
exam_scores <- round(pmin(100, pmax(0, exam_scores)))  # Bound 0-100

n <- length(exam_scores)

cat("Sample size:", n, "\n")
cat("Mean:", round(mean(exam_scores), 2), "\n")
cat("SD:", round(sd(exam_scores), 2), "\n\n")

# -------- DESCRIPTIVE STATISTICS --------

cat("=== DESCRIPTIVE STATISTICS ===\n")
cat("Min:", min(exam_scores), "\n")
cat("Q1:", quantile(exam_scores, 0.25), "\n")
cat("Median:", median(exam_scores), "\n")
cat("Q3:", quantile(exam_scores, 0.75), "\n")
cat("Max:", max(exam_scores), "\n\n")

# Skewness and kurtosis
skew <- skewness(exam_scores)
kurt_excess <- kurtosis(exam_scores) - 3  # Excess kurtosis

cat("Skewness:", round(skew, 3), 
    ifelse(abs(skew) < 0.5, "(approximately symmetric)",
           ifelse(skew > 0, "(right-skewed)", "(left-skewed)")), "\n")
cat("Excess kurtosis:", round(kurt_excess, 3),
    ifelse(abs(kurt_excess) < 1, "(approximately mesokurtic)",
           ifelse(kurt_excess > 0, "(heavy tails)", "(light tails)")), "\n\n")

# -------- SHAPIRO-WILK TEST --------

shapiro_result <- shapiro.test(exam_scores)

cat("=== SHAPIRO-WILK TEST ===\n")
cat("W statistic:", round(shapiro_result$statistic, 4), "\n")
cat("p-value:", format.pval(shapiro_result$p.value, digits = 3), "\n\n")

# Interpret W statistic
if (shapiro_result$statistic > 0.95) {
  cat("W > 0.95: Good fit to normal distribution\n")
} else if (shapiro_result$statistic > 0.90) {
  cat("0.90 < W < 0.95: Moderate fit\n")
} else {
  cat("W < 0.90: Poor fit to normal distribution\n")
}

# Interpret p-value
if (shapiro_result$p.value >= 0.05) {
  cat("p >= 0.05: Fail to reject H₀; data consistent with normality\n\n")
} else {
  cat("p < 0.05: Reject H₀; significant departure from normality\n\n")
}

# -------- GRAPHICAL DIAGNOSTICS --------

par(mfrow = c(2, 2))

# 1. Histogram with normal overlay
hist(exam_scores, breaks = 10, freq = FALSE, 
     main = "Histogram with Normal Curve",
     xlab = "Exam Score", col = "lightblue", border = "white")
curve(dnorm(x, mean = mean(exam_scores), sd = sd(exam_scores)),
      add = TRUE, col = "red", lwd = 2)

# 2. Q-Q plot
qqnorm(exam_scores, main = "Q-Q Plot", pch = 19, col = "steelblue")
qqline(exam_scores, col = "red", lwd = 2)

# 3. Boxplot
boxplot(exam_scores, main = "Boxplot", ylab = "Exam Score", 
        col = "lightgreen", horizontal = FALSE)

# 4. Density plot
plot(density(exam_scores), main = "Kernel Density Estimate",
     xlab = "Exam Score", lwd = 2, col = "darkblue")
rug(exam_scores, col = "gray")

par(mfrow = c(1, 1))

# -------- OUTLIER CHECK --------

z_scores <- scale(exam_scores)
outliers <- which(abs(z_scores) > 3)

cat("=== OUTLIER CHECK ===\n")
if (length(outliers) > 0) {
  cat("Outliers detected(|z| > 3):\n")
  print(data.frame(
    Index = outliers,
    Score = exam_scores[outliers],
    Z = round(z_scores[outliers], 2)
  ))
  cat("\n")
} else {
  cat("No extreme outliers(|z| > 3) detected\n\n")
}

# -------- SENSITIVITY ANALYSIS --------

if (length(outliers) > 0) {
  scores_no_outliers <- exam_scores[-outliers]
  shapiro_no_outliers <- shapiro.test(scores_no_outliers)
  
  cat("=== SENSITIVITY ANALYSIS(without outliers) ===\n")
  cat("W statistic:", round(shapiro_no_outliers$statistic, 4), "\n")
  cat("p-value:", format.pval(shapiro_no_outliers$p.value, digits = 3), "\n")
  cat("Change in W:", 
      round(shapiro_no_outliers$statistic - shapiro_result$statistic, 4), "\n\n")
}

# -------- INTERPRETATION & RECOMMENDATIONS --------

cat("=== INTERPRETATION ===\n")
if (shapiro_result$p.value >= 0.05 && abs(skew) < 0.5 && abs(kurt_excess) < 1) {
  cat("Conclusion: Data are approximately normally distributed.\n")
  cat("Recommendation: Parametric tests(t-test, ANOVA) are appropriate.\n")
} else if (shapiro_result$p.value < 0.05) {
  cat("Conclusion: Data show significant departure from normality.\n")
  cat("Primary issue:",
      ifelse(abs(skew) > 1, "Skewness",
             ifelse(abs(kurt_excess) > 2, "Kurtosis", "Other")), "\n")
  cat("Recommendation: Consider(1) transformation(log, sqrt),\n")
  cat("(2) nonparametric tests, or(3) robust methods.\n")
  if (n >= 30) {
    cat("Note: With n ≥ 30, CLT may justify parametric tests for means\n")
    cat("despite mild non-normality(check Q-Q plot).\n")
  }
} else {
  cat("Conclusion: Borderline normality(non-significant but W < 0.95).\n")
  cat("Recommendation: Examine Q-Q plot closely. Consider robust methods.\n")
}