19. Building a p-value by hand — the permutation test · Dr. Seok's Tumour-Infiltrating Lymphocytes, First Steps

In lesson 14 we got z and then a p-value from 2 × norm.sf(|z|). That one line hides a large assumption — that z follows a standard normal distribution.

The assumption holds when the sample is large enough and conditions are right. But how do we know ours is?

By building it ourselves.

Stomach was p = 0.87 in lesson 14 — a place with nothing in it. Good for watching a situation where the null is true.

Python Runs in your browser

import til
import numpy as np

p = til.patients().dropna(subset=["OS", "OS_time"])
g = p[p.tumor == "STAD"]                   # stomach - the type with no result in lesson 14
hi, lo = til.split(g)
real_z, formula_p = til.logrank(hi.OS_time.values, hi.OS.values,
                                lo.OS_time.values, lo.OS.values)
print("%d stomach cancer patients" % len(g))
print("real z = %.3f   p from the formula = %.4f" % (real_z, formula_p))

Mean 0.01, sd 0.95, tails 4.5%. Close enough to theory (0, 1, 5%) — we have built ourselves the evidence that the normal approximation is fine here.

Python Runs in your browser

import til
import numpy as np

p = til.patients().dropna(subset=["OS", "OS_time"])
g = p[p.tumor == "STAD"]
t, e = g.OS_time.values, g.OS.values
n_hi = len(g) // 2

rng = np.random.default_rng(0)
null = []
for _ in range(1000):                       # takes a few seconds
    i = rng.permutation(len(g))
    a, b = i[:n_hi], i[n_hi:]
    null.append(til.logrank(t[a], e[a], t[b], e[b])[0])

null = np.array(null)
print("%s z values obtained by shuffling" % format(len(null), ","))
print("  mean    %+.3f   (theory 0)" % null.mean())
print("  sd       %.3f   (theory 1)" % null.std())
print("  fraction with |z| > 1.96: %.1f%%   (theory 5.0%%)" % (100*(np.abs(null) > 1.96).mean()))

Press "shuffle". The violet line is the real z; the red bars are results at least that extreme.

Loading…

The melanoma z never once appears in 1,000 shuffles. Permutation p < 0.001.

Python Runs in your browser

import til
import numpy as np
import matplotlib.pyplot as plt

p = til.patients().dropna(subset=["OS", "OS_time"])
g = p[p.tumor == "SKCM"]                  # this time, a type that did have a result
hi, lo = til.split(g)
real_z, formula_p = til.logrank(hi.OS_time.values, hi.OS.values,
                                lo.OS_time.values, lo.OS.values)

t, e = g.OS_time.values, g.OS.values
rng = np.random.default_rng(1)
null = np.array([til.logrank(t[i[:len(hi)]], e[i[:len(hi)]],
                             t[i[len(hi):]], e[i[len(hi):]])[0]
                 for i in (rng.permutation(len(g)) for _ in range(1000))])

perm_p = (np.abs(null) >= abs(real_z)).mean()
plt.figure(figsize=(8, 3.6))
plt.hist(null, bins=40, color="#3b82f6")
plt.axvline(real_z, color="#8b5cf6", lw=2.5)
plt.xlabel("z under the null"); plt.ylabel("count")
plt.tight_layout()

print("real z = %.3f" % real_z)
print("shuffles at least this extreme: %d / 1000" % int((np.abs(null) >= abs(real_z)).sum()))
print("permutation p = %.4f   (formula p = %.2e)" % (perm_p, formula_p))

Here is where a permutation test really earns its keep — you can shuffle an entire procedure, whatever it is.

Remember the cutoff hacking of lesson 15: try seventeen cut points and keep the best p. There is no formula for that procedure's p-value. But a permutation gives you one —

shuffle the labels
on the shuffled data, try all seventeen cut points and keep the smallest p
repeat 1,000 times to get the distribution of "best p"
see how extreme our own "best p" is in that distribution

That gives you an honest p-value that includes the searching. It takes a while to compute, but the principle is exactly what we just did.

🤔 Check yourself

What world does shuffling the labels create in a permutation test?

🤔 Check yourself

What is the strongest reason to use a permutation test?