Handle multiple testing and significance correctly

In the rest of this track you will not run one statistical test, you will run millions. Every position in a genome, every gene in an expression table, every alignment against a database is its own hypothesis test. The casual statistics you may have picked up (compute a p-value, check it is below 0.05, declare victory) does not just get shaky at that scale, it produces confident nonsense by design. This lesson is about the single most common and most expensive mistake in computational biology, and the honest fixes for it.

Run enough tests and something wins by luck

Start with the intuition, then make it exact. Toss a fair coin ten times and getting ten heads is a genuine surprise, about 1 chance in 1000. Now hand a thousand people ten coins each. Somebody almost certainly gets ten heads, and if you photograph only that person you have manufactured a miracle out of pure chance. Nothing changed about the coins. What changed is how many times you looked.

That is the whole multiple-testing problem, and here is the programmer's version. Picture a test suite of a million assertions where each assertion is a little flaky, failing 5 percent of the time even when the code is perfect. Run the suite once and roughly fifty thousand assertions fail, every time, on flawless code. Cherry-pick one red assertion and file it as a bug and you have invented a defect that does not exist.

The analogy is good, and here is exactly where it cracks, because the crack is the whole point. A flaky test you can simply re-run, and the false failures wash out because each run draws fresh randomness. A multiple-testing false positive does not wash out that way. Your data is fixed. Re-running the same test on the same dataset gives you the identical p-value every time, so re-running proves nothing. The false positive is not noise you can average away by looking again, it is a structural consequence of how many hypotheses you asked at once. The real cures are different: adjust the threshold to account for how many tests you ran, or go get an independent dataset and see if the hit repeats there.

What a p-value actually is, and what it is not

You cannot correct a number you have misdefined, so pin it down. A p-value is this: assume the null hypothesis is true (the null is the boring default, no real effect, no association, this variant does nothing). Given that assumption, the p-value is the probability of seeing a result at least as extreme as the one you actually observed. Small p means your data would be surprising if nothing were going on.

The correction: Bonferroni

If each test carries a chance of a false alarm, and you run many tests, the false alarms add up. Make that exact. Under the null, a test at threshold alpha fires falsely with probability alpha. Run n independent tests and the expected number of false positives is simply n times alpha. Read the worked numbers below before the printed answers, because seeing 50000 fall out of a million tests is the moment this stops being abstract.

# Multiple testing worked on a concrete GWAS-scale run.
# Test one million independent variants, each at the usual alpha = 0.05.

n_tests = 1_000_000
alpha = 0.05

# 1. Expected false positives if NOTHING is real: n times alpha.
expected_false = n_tests * alpha
print("naive expected false hits:", expected_false)     # 50000.0

# 2. Chance of AT LEAST ONE false hit (the family-wise error rate),
#    if every null is true and the tests are independent.
fwer = 1 - (1 - alpha) ** n_tests
print("chance of one or more false hits:", fwer)         # ~1.0, a near certainty

# 3. Bonferroni: to hold the FAMILY-wide error near 0.05, divide alpha by n.
bonferroni_alpha = alpha / n_tests
print("corrected per-test threshold:", bonferroni_alpha) # about 5e-08 (float repr: 5.0000000000000004e-08)

# 4. Re-check the expected false hits at the corrected threshold.
print("corrected expected false hits:", n_tests * bonferroni_alpha)  # 0.05

The Bonferroni correction is that arithmetic run backward. If you want the chance of even one false positive across the whole family of tests, the family-wise error rate, to stay at 0.05, divide your per-test threshold by the number of tests: use alpha divided by n. Then the expected false-positive count becomes n times (alpha over n), which is just alpha again, back to 0.05. It is crude and it makes almost no assumption about the tests, which is exactly why it is trusted: even when tests are correlated, the union bound guarantees Bonferroni still controls the family-wise error at or below alpha.

FDR: a gentler, more powerful question

Bonferroni is strict to a fault. Demanding almost zero chance of even one false positive throws away real signal when you expect many true effects, say thousands of genes really do change between two conditions. The false discovery rate (FDR) asks a softer question: not "did any false positive slip in" but "of the hits I declare, what fraction are false, and can I hold that fraction to, say, 5 percent." The standard method is Benjamini-Hochberg.

# Benjamini-Hochberg FDR control at q = 0.05.
# Sort p-values ascending, compare each to a moving threshold (k / m) * q,
# keep the LARGEST rank k that still clears its bar, accept all up to it.

q = 0.05
pvals = [0.0001, 0.0008, 0.0019, 0.0300, 0.0450, 0.5000]
m = len(pvals)

kept = 0
for k, p in enumerate(pvals, start=1):
    threshold = (k / m) * q
    verdict = "pass" if p <= threshold else "fail"
    if p <= threshold:
        kept = k
    print("rank", k, "p", p, "bar", round(threshold, 4), verdict)

print("BH accepts the", kept, "smallest p-values")   # 4
# Bonferroni here (0.05 / 6 = 0.0083) would accept only 3.

Notice what the moving bar buys you. Rank 4 has p = 0.0300, which Bonferroni's flat 0.0083 bar rejects, but Benjamini-Hochberg accepts it, because it is willing to let a controlled small fraction of its four discoveries be wrong in exchange for finding more of the real ones. FDR trades a little cleanliness for a lot of power. Bonferroni is the right tool when a single false positive is costly (a clinical call). FDR is the right tool for discovery-stage screens where you expect many true effects and can tolerate a known small false-discovery proportion.

E-values: the correction baked into search

You already met this idea in disguise. When you BLAST a query and it reports an e-value, that number is not a p-value, it is the expected number of alignments scoring at least this well by chance in a database of that size (recall from BIO-3.3 that the same alignment gives a bigger e-value against a bigger database). An e-value is a p-value already multiplied by the size of the search space, which makes it a Bonferroni-style correction folded directly into the score. That is why an e-value of 10 is meaningless (you expect ten chance hits that good) while 1e-40 is convincing, and why you read it as an expected count, never as a probability. The multiple-testing correction is not an afterthought bolted onto sequence search, it is built into what the search reports.

Effect size is not statistical significance

One more axis, and it is where large studies mislead even careful people. A p-value depends on two things at once: how big the real effect is, and how much data you collected. The p-value roughly sharpens with the square root of your sample size, so with enough samples any nonzero effect, no matter how trivial, eventually crosses any significance threshold you like. Significance answers "can I distinguish this from zero." It does not answer "is it big enough to care about." That second question is the effect size: the odds ratio, the regression coefficient, the fraction of variance explained.

Why genome-wide studies use 5e-8

Now the threshold you see everywhere in genome-wide association studies (GWAS) falls straight out of the arithmetic above. The naive move is Bonferroni over every common variant, but here the flaky-test analogy cracks a second time, in a way worth naming. Flaky tests are independent. Genetic variants are not: nearby variants are inherited together in blocks (a correlation called linkage disequilibrium), so testing ten million correlated variants is nowhere near ten million independent questions. The number of effectively independent common-variant tests across the human genome works out to roughly one million. Bonferroni on that, 0.05 divided by one million, gives 5e-8 (that is 5 times 10 to the minus 8). The field fixed on that threshold precisely because it corrects for the multiple-testing burden at genome scale, and it holds regardless of exactly how many variants a particular chip happened to type. A GWAS hit at p less than 5e-8 has cleared a bar built to survive a million chances at a fluke.

Where this leaves you

At genome scale, significance stops being a single check and becomes an accounting problem. A p-value means the probability of the data under a world where nothing is happening, never the probability your hit is real. Across a million tests, expect fifty thousand false positives at 0.05, so a lone small p-value is expected noise until it clears a corrected bar. Bonferroni (alpha over n) controls the chance of any false positive, FDR controls the fraction of false discoveries and buys back power, e-values are that correction baked into search, and effect size is a separate axis that significance never measures. Genome-wide 5e-8 is just Bonferroni for a million independent tests. Carry one habit out of this lesson: whenever you see a p-value, ask how many tests it won against, and whether it has replicated.