# ============================================================================ # 04-09 Hypothesis Testing Toolkit # Chapters 10-01, 10-02, 11-01, 11-02, 11-03 # # Run from the project folder: Rscript r/analysis.R # ============================================================================ fill <- read.csv("data/fill_volume.csv")$ml methods <- read.csv("data/two_methods.csv", stringsAsFactors = TRUE) pairs <- read.csv("data/before_after.csv") report <- function(name, H0, H1, stat_name, stat, df, pval, crit, alpha, effect, ci, conclusion) { cat("\n", strrep("=", 74), "\n", sep = "") cat(name, "\n") cat(strrep("=", 74), "\n", sep = "") cat("H0:", H0, "\n") cat("H1:", H1, "\n") cat(sprintf("alpha = %.2f\n\n", alpha)) cat(sprintf("%-10s = %.5f df = %s\n", stat_name, stat, ifelse(is.na(df), "-", format(round(df, 3))))) cat(sprintf("critical = %s\n", paste(round(crit, 5), collapse = " / "))) cat(sprintf("p-value = %.6f\n", pval)) cat(sprintf("decision = %s\n", ifelse(pval <= alpha, "REJECT H0", "FAIL TO REJECT H0"))) cat(sprintf("effect = %s\n", effect)) cat(sprintf("95%% CI = (%.4f, %.4f)\n", ci[1], ci[2])) cat("\n", conclusion, "\n", sep = "") } interp_d <- function(d) { a <- abs(d) if (a < 0.2) "negligible" else if (a < 0.5) "small" else if (a < 0.8) "medium" else "large" } # ════════════════════════════════════════════════════════════════════════════ # TEST 1 — ONE-SAMPLE t : is the machine filling to 250 ml? # ════════════════════════════════════════════════════════════════════════════ cat("\n>>> ASSUMPTION CHECK — fill volume\n") cat(sprintf("n = %d, Shapiro-Wilk p = %.4f, outliers = %d\n", length(fill), shapiro.test(fill)$p.value, length(boxplot.stats(fill)$out))) n1 <- length(fill); se1 <- sd(fill) / sqrt(n1) t1 <- (mean(fill) - 250) / se1 p1 <- 2 * pt(-abs(t1), n1 - 1) d1 <- (mean(fill) - 250) / sd(fill) ci1 <- mean(fill) + c(-1, 1) * qt(0.975, n1 - 1) * se1 report("TEST 1 — ONE-SAMPLE t : fill volume vs the 250 ml target", "mu = 250", "mu != 250 (two-tailed: we care about over- AND under-filling)", "t", t1, n1 - 1, p1, c(-1, 1) * qt(0.975, n1 - 1), 0.05, sprintf("Cohen's d = %.4f (%s)", d1, interp_d(d1)), ci1, sprintf(paste0("The mean fill volume is %.2f ml. There %s sufficient evidence at\n", "the 5%% level that the machine differs from its 250 ml target."), mean(fill), ifelse(p1 <= 0.05, "IS", "is NOT"))) cat("\nnonparametric check (Wilcoxon signed-rank): p =", round(wilcox.test(fill, mu = 250)$p.value, 5), "\n") # ════════════════════════════════════════════════════════════════════════════ # TEST 2 — TWO-SAMPLE t : method A vs method B # ════════════════════════════════════════════════════════════════════════════ A <- methods$score[methods$method == "A"] B <- methods$score[methods$method == "B"] cat("\n>>> ASSUMPTION CHECK — two methods\n") cat(sprintf("n_A = %d (Shapiro p = %.4f), n_B = %d (Shapiro p = %.4f)\n", length(A), shapiro.test(A)$p.value, length(B), shapiro.test(B)$p.value)) vt <- var.test(A, B) cat(sprintf("variance ratio = %.3f, F-test p = %.4f -> %s\n", var(A) / var(B), vt$p.value, ifelse(vt$p.value > 0.05, "variances are compatible; either test is defensible", "variances DIFFER; use Welch"))) w <- t.test(A, B) # Welch — R's default and the safe choice sp2 <- ((length(A)-1)*var(A) + (length(B)-1)*var(B)) / (length(A)+length(B)-2) d2 <- (mean(A) - mean(B)) / sqrt(sp2) report("TEST 2 — TWO-SAMPLE t (Welch) : method A vs method B", "mu_A = mu_B", "mu_A != mu_B", "t", w$statistic, w$parameter, w$p.value, c(-1, 1) * qt(0.975, w$parameter), 0.05, sprintf("Cohen's d = %.4f (%s)", d2, interp_d(d2)), w$conf.int, sprintf(paste0("Method A averaged %.2f and method B %.2f, a difference of %.2f.\n", "There %s sufficient evidence at the 5%% level of a real difference."), mean(A), mean(B), mean(A) - mean(B), ifelse(w$p.value <= 0.05, "IS", "is NOT"))) cat("\npooled t-test for comparison: p =", round(t.test(A, B, var.equal = TRUE)$p.value, 5), "\n") cat("Mann-Whitney U (nonparametric): p =", round(wilcox.test(A, B)$p.value, 5), "\n") # ════════════════════════════════════════════════════════════════════════════ # TEST 3 — PAIRED t : reaction time before vs after # ════════════════════════════════════════════════════════════════════════════ d <- pairs$before_ms - pairs$after_ms # positive = improvement cat("\n>>> ASSUMPTION CHECK — the DIFFERENCES (not the originals)\n") cat(sprintf("n pairs = %d, Shapiro-Wilk p = %.4f, outliers = %d\n", length(d), shapiro.test(d)$p.value, length(boxplot.stats(d)$out))) cat(sprintf("correlation between before and after: %.4f -> pairing is worth it\n", cor(pairs$before_ms, pairs$after_ms))) pt_ <- t.test(pairs$before_ms, pairs$after_ms, paired = TRUE, alternative = "greater") pt2 <- t.test(pairs$before_ms, pairs$after_ms, paired = TRUE) # for the two-sided CI d3 <- mean(d) / sd(d) report("TEST 3 — PAIRED t : reaction time before vs after training", "mu_d <= 0", "mu_d > 0 (one-tailed: training is claimed to REDUCE the time)", "t", pt_$statistic, pt_$parameter, pt_$p.value, qt(0.95, pt_$parameter), 0.05, sprintf("Cohen's d = %.4f (%s)", d3, interp_d(d3)), pt2$conf.int, sprintf(paste0("Mean improvement %.2f ms. There %s sufficient evidence at the\n", "5%% level that the training reduces reaction time."), mean(d), ifelse(pt_$p.value <= 0.05, "IS", "is NOT"))) ind <- t.test(pairs$before_ms, pairs$after_ms, alternative = "greater") cat(sprintf("\nTHE COST OF IGNORING THE PAIRING\n")) cat(sprintf(" paired t = %7.4f, p = %.6f\n", pt_$statistic, pt_$p.value)) cat(sprintf(" independent t = %7.4f, p = %.6f <- much weaker\n", ind$statistic, ind$p.value)) cat("Wilcoxon signed-rank (nonparametric): p =", round(wilcox.test(pairs$before_ms, pairs$after_ms, paired = TRUE, alternative = "greater")$p.value, 6), "\n") # ════════════════════════════════════════════════════════════════════════════ # TEST 4 — Z-TEST FOR A PROPORTION # ════════════════════════════════════════════════════════════════════════════ poll_path <- "../04-08-confidence-interval-estimator/data/poll_responses.csv" if (file.exists(poll_path)) { poll <- read.csv(poll_path)$supports x <- sum(poll == "Yes"); n4 <- length(poll); p0 <- 0.50 cat("\n>>> ASSUMPTION CHECK — proportion\n") cat(sprintf("n*p0 = %.0f and n*(1-p0) = %.0f (both >= 5)\n", n4 * p0, n4 * (1 - p0))) se4 <- sqrt(p0 * (1 - p0) / n4) # p0, NOT p-hat ph <- x / n4 z4 <- (ph - p0) / se4 p4 <- 2 * pnorm(-abs(z4)) se_ci <- sqrt(ph * (1 - ph) / n4) # p-hat for the INTERVAL ci4 <- ph + c(-1, 1) * qnorm(0.975) * se_ci report("TEST 4 — Z-TEST FOR A PROPORTION : is support different from 50%?", "p = 0.50", "p != 0.50", "z", z4, NA, p4, c(-1.959964, 1.959964), 0.05, sprintf("p-hat - p0 = %+.4f (%.1f percentage points)", ph - p0, 100 * (ph - p0)), ci4, sprintf(paste0("%d of %d respondents (%.1f%%) said Yes. There %s sufficient\n", "evidence at the 5%% level that support differs from 50%%."), x, n4, 100 * ph, ifelse(p4 <= 0.05, "IS", "is NOT"))) cat("\nNOTE the two different standard errors:\n") cat(sprintf(" TEST uses p0: SE = %.6f\n", se4)) cat(sprintf(" INTERVAL uses p-hat: SE = %.6f\n", se_ci)) } else { cat("\n(Test 4 skipped: poll_responses.csv lives in project 04-08.)\n") } # ── PLOTS ─────────────────────────────────────────────────────────────────── png("testing_plots.png", width = 1100, height = 850) par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1)) boxplot(fill, horizontal = TRUE, col = "#8A5FBF", main = "Test 1: fill volume vs 250 ml target", xlab = "ml") abline(v = 250, col = "#b4122e", lwd = 2) abline(v = mean(fill), col = "#0FA3A3", lwd = 2, lty = 2) boxplot(score ~ method, data = methods, col = c("#5B2A86", "#0FA3A3"), main = "Test 2: method A vs B", ylab = "score") matplot(rbind(pairs$before_ms, pairs$after_ms), type = "b", pch = 19, lty = 1, col = "#5B2A86", xaxt = "n", xlab = "", ylab = "reaction time (ms)", main = "Test 3: each line is one subject") axis(1, at = 1:2, labels = c("Before", "After")) hist(d, breaks = 8, col = "#0B7A7A", border = "white", main = "Test 3: differences (before - after)", xlab = "ms improved") abline(v = 0, col = "#b4122e", lwd = 2) abline(v = mean(d), col = "#0FA3A3", lwd = 2, lty = 2) par(mfrow = c(1, 1)) dev.off() cat("\nWrote testing_plots.png\n")