16. Test twelve times and — correcting for multiplicity · Dr. Seok's Tumour-Infiltrating Lymphocytes, First Steps

In lesson 14 we tested each cancer type. That is 12 tests. Three of them came in under p < 0.05. How should that be read?

Get a feel for it first. With no relationship at all, how many of 12 tests come out significant?

Even with no relationship at all, 3 or more of 12 come out significant about 2% of the time.

Python Runs in your browser

import numpy as np

rng = np.random.default_rng(0)
# repeat "12 tests on two completely unrelated values" 2,000 times
counts = []
for _ in range(2000):
    ps = []
    for _ in range(12):
        a, b = rng.normal(size=60), rng.normal(size=60)
        from scipy import stats
        ps.append(stats.mannwhitneyu(a, b).pvalue)
    counts.append(sum(p < 0.05 for p in ps))

counts = np.array(counts)
print("out of 12 tests, number significant - mean %.2f (theory 0.60)" % counts.mean())
print()
for k in range(5):
    print("  probability of exactly %d: %5.1f%%" % (k, 100*(counts == k).mean()))
print("  probability of 3 or more:  %5.1f%%" % (100*(counts >= 3).mean()))

2%. Low, but not zero. And real studies often run not 12 tests but hundreds. Test 20,000 genes one by one? A thousand of them become significant for free.

Hence correction. Two methods are in wide use.

3 before correction → 2 after BH and 2 after Bonferroni. Breast cancer drops out.

Python Runs in your browser

import til
import numpy as np

p = til.patients()
rows = []
for t, g in p.dropna(subset=["OS", "OS_time"]).groupby("tumor"):
    r = til.surv_test(g)
    if r:
        rows.append([t, r["z"], r["p"]])
rows.sort(key=lambda x: x[2])
ps = np.array([r[2] for r in rows])
k = len(ps)

bonf = np.minimum(ps * k, 1)
bh   = til.bh(ps)

print("%d tests" % k)
print("%-6s%7s%11s%11s%11s" % ("type", "z", "p", "Bonf", "BH(q)"))
for (t, z, pv), b, q in zip(rows, bonf, bh):
    star = "  *" if q < 0.05 else ""
    print("%-6s%7.2f%11.2e%11.2e%11.2e%s" % (t, z, pv, b, q, star))
print()
print("significant uncorrected: %d" % int((ps < .05).sum()))
print("significant after Bonferroni: %d" % int((bonf < .05).sum()))
print("significant after BH: %d" % int((np.array(bh) < .05).sum()))

Switch the correction method. The same data reaches a different conclusion.

Loading…

To summarise —

SKCM — survives any correction (BH q = 1.4×10⁻⁴). A solid result.
CESC — survives BH (q = 0.024) and squeaks past Bonferroni (0.047).
BRCA — p = 0.017 uncorrected, but BH q = 0.067.

The BRCA result is not thereby false. It means "as one of twelve, it is not persuasive enough". Confirming it needs another cohort.

🤔 Check yourself

With no relationship at all, how many of 12 tests come out at p < 0.05 on average?

🤔 Check yourself

Why is BH less stingy than Bonferroni?