Break the relationship on purpose, a thousand times, and see where reality lands.
The JonStats course opens its treatment with a story worth keeping intact. Coin A was flipped 12 times, each result written on a red card; Coin B 8 times, onto blue cards. The coins are gone — the cards are all you have. Red shows 9 heads of 12 ($0.750$); blue shows 3 of 8 ($0.375$). The difference is $0.375$. Were the two coins actually different, or is a gap that size just what chance does with twenty flips?
The permutation move: if the coins were identical (the null hypothesis), then the colour of each card is an irrelevant label — any assignment of these twenty outcomes to colours is equally plausible. So shuffle the outcomes across the colours (without replacement — every card keeps existing exactly once), recompute the difference, repeat. That builds the distribution of differences a boring world produces. Then see where the real $0.375$ falls in it.
The p-value before it had a name. Run the full thousand: 106 shuffles produce a difference of $0.375$ or more, so $p = 0.106$ — not a table lookup, a count. The exact answer (Fisher's test, which enumerates every possible shuffle rather than sampling them) is $0.113$. The classical $\chi^2$ test on this data, incidentally, comes with a warning: half its cells hold fewer than five observations. The shuffle needs no such apology — small samples are where it is most at home.
Nothing about shuffling is specific to coins. To test whether exercise predicts resting heart rate at all, shuffle the heart rates across the exercise values — deliberately destroying any real X–Y link while keeping both variables' own distributions intact — and refit the line each time. The slopes of those broken datasets form the null distribution; the real slope either sits among them or it doesn't.
Not one in a thousand. The full run's most extreme null slope falls well short of the observed $-1.271$: zero of 1,000 shuffles get there, so $p < 0.001$. Note the null distribution's spread ($\text{sd} \approx 0.337$) — the shuffle has rebuilt, from scratch, roughly the standard error that the curvature route derived analytically ($0.29$), which is exactly what it should do under the null.
Counting extremes: for the cards we counted one tail (the question was “is Coin A better?”); for the slope we counted both tails ($|\beta_1^*| \ge |\hat\beta_1|$ — a relationship in either direction would have been news). Decide which question you are asking before you count. And a historical footnote from JonStats: the $\chi^2$ and $t$ distributions were, in effect, labour-saving devices for exactly these shuffles — algebra once, instead of a thousand permutations by hand. Silicon has made the original method cheap again.
Real runs (native random streams, so digits differ slightly from the seeded interactives):
# The cards (after the JonStats base-R example)
outcomes <- c(1,1,1,0,1,1,0,1,1,1,0,1, 0,1,0,0,1,0,1,0)
set.seed(42)
nulls <- replicate(1000, {
s <- sample(outcomes) # shuffle = sample w/o replacement
mean(s[1:12]) - mean(s[13:20])
})
mean(nulls >= 0.375) # one-sided p, by counting
# The regression slope
set.seed(42)
perm_slopes <- replicate(1000, coef(lm(sample(hr$y) ~ hr$x))[2])
sum(abs(perm_slopes) >= 1.2711) # two-sided count
# cards one-sided p: 0.109
# slope extreme count: 0 of 1000 (p < 0.001)
cards = np.array([1,1,1,0,1,1,0,1,1,1,0,1, 0,1,0,0,1,0,1,0])
rng = np.random.default_rng(42)
nulls = []
for _ in range(1000):
s = rng.permutation(cards) # shuffle = draw w/o replacement
nulls.append(s[:12].mean() - s[12:].mean())
print(np.mean(np.array(nulls) >= 0.375)) # one-sided p, by counting
rng = np.random.default_rng(42)
perm = [np.polyfit(x, rng.permutation(y), 1)[0] for _ in range(1000)]
print(int(np.sum(np.abs(perm) >= 1.2711))) # two-sided count
# cards one-sided p: 0.131
# slope extreme count: 0 of 1000 (p < 0.001)
infer-package version of the
same test.