Hypothesis tests and intervals, practically

You have measured a difference between two things — two models, two prompts, two versions of a page. This section is about the question that follows: could a gap that size have turned up by chance alone? It covers what a p-value does and does not say, the several ways it is routinely misread, how to run the standard tests in SciPy, and why a confidence interval is usually the more honest thing to put in the report.

The definition, carefully#

A p-value is the probability of seeing data at least as extreme as yours if there were really no effect. A small p-value means your data would be surprising in a world with no effect. That is all it means.

• It is not the probability that there is no effect.
• It is not the probability your result is a fluke.
• It says nothing about how big or important the effect is. A tiny effect in a huge sample can have a tiny p-value.

The definition has a consequence that makes the whole thing click, and it is easier to see by simulation than to argue. If there is genuinely no effect, the p-value is a uniformly distributed random number between 0 and 1. Not usually large — uniform. So 5 per cent of the time it lands below 0.05, half the time below 0.5, and there is nothing special about any particular value it takes. Run two identical arms against each other twenty thousand times and look.

import math
import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

def two_proportion_test(x1, n1, x2, n2):
    """Successes and trials for each arm. Returns the two-sided p-value."""
    pooled = (x1 + x2) / (n1 + n2)
    se = math.sqrt(pooled * (1 - pooled) * (1 / n1 + 1 / n2))
    z = (x2 / n2 - x1 / n1) / se
    return 2 * (1 - stats.norm.cdf(abs(z)))

null_p = np.array([
    two_proportion_test(rng.binomial(2_000, 0.03), 2_000,
                        rng.binomial(2_000, 0.03), 2_000)
    for _ in range(20_000)
])
print((null_p < 0.05).mean(), (null_p < 0.10).mean(), (null_p < 0.50).mean())
# 0.049 0.09855 0.49675

A worked A/B test#

Take a real-shaped question. The control page converts at 3.0 per cent and the variant at 3.4 per cent, a relative lift of 13 per cent, on 10,000 visitors each. Someone wants to ship it on Friday. The two-proportion test gives p equal to 0.108, so by the usual convention it does not clear the bar — and the more useful output is the interval on the difference, which runs from minus 0.09 to plus 0.89 percentage points. The data is consistent with the variant being slightly worse and with it being a 30 per cent relative improvement. Run it to 40,000 per arm and the same rates give p equal to 0.0013 and an interval of 0.16 to 0.64 points.

def ab(x1, n1, x2, n2):
    p1, p2 = x1 / n1, x2 / n2
    se = math.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    lift = p2 - p1
    return p1, p2, two_proportion_test(x1, n1, x2, n2), lift - 1.96 * se, lift + 1.96 * se

for n in (10_000, 40_000):
    p1, p2, pv, lo, hi = ab(int(0.030 * n), n, int(0.034 * n), n)
    print(n, p1, p2, round(pv, 4), round(lo * 100, 3), round(hi * 100, 3))

# 10000 0.03 0.034 0.108  -0.088  0.888
# 40000 0.03 0.034 0.0013  0.156  0.644

Two common tests#

Proportions are the common case in this field, but two others turn up constantly. Use a t-test when you are comparing averages of a continuous quantity, and pass equal_var=False, because the equal-variance version assumes something you rarely have any reason to believe. Use a chi-square test on a contingency table when you are asking whether two categorical variables are related.

group_a = rng.lognormal(mean=3.0, sigma=0.6, size=400)     # session length, minutes
group_b = rng.lognormal(mean=3.1, sigma=0.6, size=400)
t, p = stats.ttest_ind(group_a, group_b, equal_var=False)
print(round(group_a.mean(), 1), round(group_b.mean(), 1), round(float(t), 2), round(float(p), 4))
# 23.9 26.8 -2.39 0.0169

table = [[120, 880],    # plan A: churned, stayed
         [ 90, 910]]    # plan B
chi2, p2, dof, expected = stats.chi2_contingency(table)
print(round(float(chi2), 2), round(float(p2), 4), dof)
print(np.round(expected, 1).tolist())
# 4.47 0.0344 1
# [[105.0, 895.0], [105.0, 895.0]]

Prefer an interval#

A confidence interval gives a range of effect sizes consistent with the data. 'Plan B churn is 3 points lower, 95% interval 0.2 to 5.8' tells you the direction, the size and the uncertainty at once, which a p-value alone does not.

It also makes the two honest failure modes distinguishable, and a p-value does not. An interval from minus 0.1 to plus 0.9 points says you have not learned much and should keep the experiment running. An interval from minus 0.05 to plus 0.06 points says you have learned something definite: whatever the effect is, it is too small to care about. Both of those come back as not significant, and they call for opposite decisions.

Power#

Power is the chance your study detects an effect of a given size if it exists. Small studies have low power, so they miss real effects, and the effects they do report tend to be exaggerated. Estimate the sample you need before collecting data.

Power is simulable, which means you never have to trust a formula you half remember. Decide the effect you would want to catch, generate data with that effect in it, run your actual test, and count how often it fires. Here the effect is the 3.0 against 3.4 per cent from above, and the question is how many visitors per arm it takes to see it reliably.

for n in (5_000, 20_000, 80_000):
    wins = sum(
        two_proportion_test(rng.binomial(n, 0.030), n, rng.binomial(n, 0.034), n) < 0.05
        for _ in range(4_000)
    )
    print(n, wins / 4_000)
# 5000 0.2125
# 20000 0.61925
# 80000 0.99425

The closed form agrees and is quicker for planning. For two proportions, the sample per arm is 2.8 squared — that constant being 1.96 for the false positive rate plus 0.84 for 80 per cent power — times the sum of the two variances, divided by the squared difference. For these numbers it says 30,353 per arm, and simulating at exactly that n gives 80.2 per cent power.

p1, p2 = 0.030, 0.034
n_needed = (1.96 + 0.84) ** 2 * (p1 * (1 - p1) + p2 * (1 - p2)) / (p2 - p1) ** 2
print(round(n_needed))       # 30353

wins = sum(
    two_proportion_test(rng.binomial(30_353, p1), 30_353, rng.binomial(30_353, p2), 30_353) < 0.05
    for _ in range(4_000)
)
print(wins / 4_000)          # 0.80225

Twenty metrics, one winner#

The most common way a team fools itself is not a bad test, it is many tests. Every dashboard has twenty metrics, and each one gets a 5 per cent chance of a false alarm, so the chance that at least one fires when nothing at all has changed is one minus 0.95 to the twentieth, which is 64 per cent. On average one metric will be significant every time you look.

print(round(1 - 0.95 ** 20, 4))          # 0.6415

flags = rng.random((50_000, 20)) < 0.05   # 20 metrics, nothing happening
print(flags.any(axis=1).mean(), flags.sum(axis=1).mean())
# 0.64172 0.99986