Hypothesis testing

Deep dive

How it works

Parts:

• hypotheses — H0 (no difference) and H1 (one-sided mu_A > mu_B or two-sided !=), written before looking
• test statistic — one number summarising the evidence: t, chi2, a rank sum, a mean difference
• null distribution — law of the statistic if H0 holds: closed form (t, chi-square) or built by shuffling labels
• decision rule — p-value P(T at least as extreme | H0) against alpha
• error budget — alpha (Type I), beta (Type II), a correction when many tests share it
• design — sample size from alpha, power, effect size, variance; the unit of independence

Steps:

1. Fix H0, H1, sidedness, alpha (0.05 is convention), and the smallest effect worth detecting delta.
2. Size: n = 2 (z_a + z_b)^2 sigma^2 / delta^2 per group (z_a = 1.960 two-sided, z_b = 0.842 for 80 % power; sigma from history).
3. Pick the test (table); check that units are independent.
4. Welch: t = (m1 - m2) / sqrt(v1 + v2), vi = si^2/ni, df = (v1+v2)^2 / (v1^2/(n1-1) + v2^2/(n2-1)).
5. Get p from the null law (t with df, or a permutation count).
6. Reject if p <= alpha; with m tests apply Holm or Benjamini-Hochberg.
7. Report effect, CI and p. The interval d +- t(df, 0.975) se holds exactly the delta0 that a two-sided 0.05 test does not reject.

 data ─► statistic T ─► null law of T ──► p = P(|T| >= |t_obs| | H0)
  │                      ▲ (formula | permutation | bootstrap)       │
  └─ units independent?                                           ▼
                                              p <= alpha ? reject H0 : keep
 many tests: sort p ─► Holm (FWER) or BH (FDR) ─► rejection set

Which test

Multiple comparisons. Bonferroni: test each at alpha/m (controls FWER, the chance of any false rejection). Holm (1979): sort p ascending, reject while p_(i) <= alpha/(m - i + 1), stop at the first failure; never weaker than Bonferroni, no extra assumptions. Benjamini-Hochberg (1995): reject the k smallest, k the largest i with p_(i) <= i q/m; controls the expected share of false rejections (FDR) at q for independent tests (positive dependence: Benjamini-Yekutieli 2001).

Code

Welch statistic and df, a seeded permutation test, Benjamini-Hochberg, and checks: a 1-sd shift is detected, 2000 null runs reject within 4 binomial SDs of 0.05, BH matches a hand-worked example. Compiled and run; permutation p-values need no t-CDF.

use std::time::Instant;

struct Rng(u64); // xorshift64, seeded
impl Rng {
    fn next(&mut self) -> u64 { let mut x = self.0; x ^= x << 13; x ^= x >> 7; x ^= x << 17; self.0 = x; x }
    fn unif(&mut self) -> f64 { ((self.next() >> 11) as f64 + 0.5) / (1u64 << 53) as f64 }
    fn below(&mut self, n: usize) -> usize { (self.next() % n as u64) as usize } // bias ~ n/2^64, negligible
    fn normal(&mut self) -> f64 { (-2.0 * self.unif().ln()).sqrt() * (6.283_185_307_179_586 * self.unif()).cos() }
}
fn mean_var(x: &[f64]) -> (f64, f64) {
    let m = x.iter().sum::<f64>() / x.len() as f64;
    (m, x.iter().map(|v| (v - m).powi(2)).sum::<f64>() / (x.len() - 1) as f64)
}
/// Welch's t statistic and Welch-Satterthwaite degrees of freedom (no equal-variance assumption).
fn welch(a: &[f64], b: &[f64]) -> (f64, f64) {
    let ((ma, va), (mb, vb)) = (mean_var(a), mean_var(b));
    let (sa, sb) = (va / a.len() as f64, vb / b.len() as f64);
    let df = (sa + sb).powi(2) / (sa * sa / (a.len() - 1) as f64 + sb * sb / (b.len() - 1) as f64);
    ((ma - mb) / (sa + sb).sqrt(), df)
}
/// Two-sided permutation test, difference in means: p = (1 + #{|d*| >= |d|}) / (B + 1).
fn perm_test(a: &[f64], b: &[f64], reps: usize, rng: &mut Rng) -> f64 {
    let (na, n) = (a.len(), a.len() + b.len());
    let mut pool: Vec<f64> = a.iter().chain(b).copied().collect();
    let total: f64 = pool.iter().sum();
    let diff = |sum_a: f64| sum_a / na as f64 - (total - sum_a) / (n - na) as f64;
    let obs = diff(a.iter().sum()).abs();
    let mut hits = 0;
    for _ in 0..reps {
        for i in 0..na { let j = i + rng.below(n - i); pool.swap(i, j); } // partial Fisher-Yates: first na = group A
        if diff(pool[..na].iter().sum()).abs() >= obs - 1e-12 { hits += 1; }
    }
    (1 + hits) as f64 / (reps + 1) as f64
}
/// Benjamini-Hochberg step-up: reject the k smallest p-values, k = max{i : p_(i) <= i q / m}.
fn bh(p: &[f64], q: f64) -> Vec<bool> {
    let m = p.len();
    let mut idx: Vec<usize> = (0..m).collect();
    idx.sort_by(|&i, &j| p[i].partial_cmp(&p[j]).unwrap());
    let k = (0..m).rev().find(|&r| p[idx[r]] <= (r + 1) as f64 * q / m as f64).map_or(0, |r| r + 1);
    let mut rej = vec![false; m];
    for &i in &idx[..k] { rej[i] = true; }
    rej
}
fn sample(rng: &mut Rng, n: usize, mu: f64, sd: f64) -> Vec<f64> { (0..n).map(|_| mu + sd * rng.normal()).collect() }

fn main() {
    let mut rng = Rng(0x9E37_79B9_7F4A_7C15);
    // Real shift (1 sd, n = 30 per group): tiny p.
    let (a, b) = (sample(&mut rng, 30, 1.0, 1.0), sample(&mut rng, 30, 0.0, 1.0));
    let (t, df) = welch(&a, &b);
    let p = perm_test(&a, &b, 9_999, &mut rng);
    println!("real shift: t = {t:.2}, df = {df:.1}, permutation p = {p:.4}");
    assert!(p < 0.01 && df > 29.0 && df <= 58.0);
    let (_, df2) = welch(&sample(&mut rng, 10, 0.0, 5.0), &sample(&mut rng, 100, 0.0, 1.0));
    assert!(df2 < 15.0); // noisy small group: df falls toward its n - 1
    // No shift: false-positive rate within 4 binomial SDs of alpha.
    let (runs, reps, alpha) = (2000usize, 999usize, 0.05);
    let t0 = Instant::now();
    let fp = (0..runs).filter(|_| {
        let (x, y) = (sample(&mut rng, 20, 0.0, 1.0), sample(&mut rng, 20, 0.0, 1.0));
        perm_test(&x, &y, reps, &mut rng) <= alpha
    }).count();
    let secs = t0.elapsed().as_secs_f64();
    let rate = fp as f64 / runs as f64;
    let tol = 4.0 * (alpha * (1.0 - alpha) / runs as f64).sqrt();
    println!("null: false-positive rate {rate:.4} (tolerance +-{tol:.4}); {:.0} permutations/s", (runs * reps) as f64 / secs);
    assert!((rate - alpha).abs() < tol);
    // Benjamini-Hochberg at q = 0.05, m = 6: step-up rejects five; Bonferroni (0.05/6) only one.
    let p = [0.04, 0.005, 0.5, 0.03, 0.011, 0.02];
    assert_eq!(bh(&p, 0.05), [true, true, false, true, true, true]);
    assert_eq!(p.iter().filter(|&&v| v <= 0.05 / 6.0).count(), 1);
    println!("ok");
}

Measured on this Mac (Apple M4 Pro, one core, release): real shift t = 4.79, df = 57.3, p = 0.0001 (the floor 1/(B+1), B = 9,999); null false-positive rate 0.054 over 2000 runs (B = 999, n = 20 per group); about 21-24 million permutations/s at n = 40 (three runs).

Complexity and constants

• t, z: one pass, O(1) state. Rank tests, BH, Holm: a sort.
• Permutation: O(B n); Monte Carlo error sqrt(p(1-p)/B), about 0.007 at p = 0.05, B = 999; p cannot go below 1/(B+1).
• Halving the detectable effect quadruples n; (z_a + z_b)^2 is 7.85 at 80 % power, 10.5 at 90 %.

Efficiency

• Cost: Welch is one memory-bound pass; a permutation costs about 43 ns at n = 40 (measured above, one core); the real cost is the data collected, n scales with sigma^2 / delta^2.
• Levers:
  • variance reduction (CUPED, paired or stratified design) — variance times 1 - rho^2 (Deng et al. 2013) — a correlated covariate exists
  • closed-form p instead of resampling — B-fold less compute — null law known
  • sequential designs — fewer expected samples at equal error rates (SPRT) — data arrive over time
• Naive → standard → state of the art:
  • test until p < 0.05 · 14.5 % false positives at 5 looks (measured)
  • pre-registered fixed-n Welch or rank test, Holm / BH · alpha held, n from power
  • variance-reduced or group-sequential design · fewer units at the same alpha and power

Implementation notes

• Default to Welch; pooled t gains little unless variances and sizes are equal.
• Permutation p-value: (1 + hits)/(B + 1), never hits/B (p = 0 breaks alpha).
• Seed resampling tests and store seed and B with the result.
• Compare floats with a tolerance; a different summation order creates false ties.

Pitfalls

• Reading p as P(H0 | data), "chance alone", or effect size: ASA principles 2 and 5.
• Peeking and optional stopping: five looks give about 14.5 % false positives (measured in A/B testing); use sequential methods.
• p-hacking: many metrics, subgroups or outlier rules, reporting the one that works; fix the plan or correct for all (ASA principle 4).
• Non-independent units (sessions per user, patients per clinic): analyse at cluster level or the standard error is too small.
• Statistical vs practical significance: huge n makes trivial effects significant; report the CI.
• Not rejecting is not "no effect"; check power or CI width. Choosing sidedness after seeing the sign doubles alpha.
• Permuting means with unequal variances and sizes breaks exchangeability and over-rejects; use a studentized statistic or bootstrap.

Variants and alternatives

In production

• Clinical trials: pre-registered endpoints, alpha spent across interim looks, multiplicity control.
• Manufacturing: two-sample process comparison (NIST/SEMATECH handbook ch. 7).
• Experiment platforms: Welch-style tests per metric (layering); genomics: per-gene tests with BH q-values.
• Software performance gates: branch versus main, corrected across benchmarks.
• Model evaluation: paired tests on overlapping folds need care (time-series CV).

Implementations

Sources

• B. L. Welch, The generalization of "Student's" problem when several different population variances are involved, Biometrika 34(1-2), 1947. paper
• S. Holm, A simple sequentially rejective multiple test procedure, Scandinavian Journal of Statistics 6(2), 65-70, 1979.
• Y. Benjamini, Y. Hochberg, Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing, JRSS B 57(1), 289-300, 1995. paper
• Y. Benjamini, D. Yekutieli, The control of the false discovery rate in multiple testing under dependency, Annals of Statistics 29(4), 1165-1188, 2001. paper
• R. L. Wasserstein, N. A. Lazar, The ASA Statement on p-Values: Context, Process, and Purpose, The American Statistician 70(2), 129-133, 2016. paper · ASA release with the six principles
• NIST/SEMATECH, e-Handbook of Statistical Methods, section 7.3.1 "Do two processes have the same mean?". handbook
• A. Deng, Y. Xu, R. Kohavi, T. Walker, Improving the Sensitivity of Online Controlled Experiments by Utilizing Pre-Experiment Data, WSDM 2013. paper