Run and analyze an A/B test, and report it honestly

Sessions b1-b3 got you a trustworthy test: a clean estimand, enough power, and a randomization you can defend. Now the data has landed and someone wants a number by end of day. The failure mode is not the arithmetic - statsmodels does that in one line. The failure mode is the story around it: quoting a lift with no interval, confusing "statistically significant" with "worth shipping", and running an eye over five metrics until one crosses 0.05. Today we do the whole readout the way a staff data scientist would - the test, the interval, the effect size, and the multiple-comparisons correction that keeps you honest.

Conversion is a binary outcome per user - bought or did not - so each arm is a proportion. To ask whether variant 3.6% beats control 3.2%, we compare two proportions. The z-test pools the two arms into one rate under the null (they are equal), builds a standard error from that pooled rate, and measures how many standard errors apart the observed difference sits. Big z, small p. But the p-value alone never ships a change - you also need the interval on the difference, because that is what tells you the plausible size of the effect.

Every test at α = 0.05 has a 5% chance of a false positive when nothing is going on. Test one metric, that risk is 5%. Test 20 metrics, you expect about one to cross 0.05 by pure luck - and the more you slice (per device, per geo, per product line), the surer you are to "find" something. This is the multiple-comparisons problem, and it is why dashboards full of green stars are so often noise. Two corrections tame it: Bonferroni tests each metric at α/m, which controls the family-wise error rate but is conservative; Benjamini-Hochberg controls the false discovery rate and is more powerful when you have many metrics.

The test came back: control 3.2% (1024 of 32000), variant 3.6% (1152 of 32000). We run the z-test and, in the same breath, the confidence interval on the difference - because a p-value without an interval is half a readout.

import numpy as np
from statsmodels.stats.proportion import proportions_ztest, confint_proportions_2indep

# Lumen hero-first product-page test, primary metric = checkout conversion
n_ctrl, n_var       = 32_000, 32_000
conv_ctrl, conv_var = 1024,   1152          # 3.2% vs 3.6%

count = np.array([conv_var, conv_ctrl])     # successes, variant first
nobs  = np.array([n_var,   n_ctrl])         # trials

z, p = proportions_ztest(count, nobs, alternative='two-sided')
print(f"z = {z:.3f}   p = {p:.4f}")         # z ≈ 2.79   p ≈ 0.005

low, high = confint_proportions_2indep(
    conv_var, n_var, conv_ctrl, n_ctrl,
    compare='diff', alpha=0.05)

print(f"diff 95% CI = [{low:.4f}, {high:.4f}]")
# ≈ [0.0012, 0.0068]  ->  [0.12pp, 0.68pp], excludes 0

A single "the variant won" is not a report. We compute both flavours of lift and put an interval on the relative lift, then say out loud where statistical significance and practical significance diverge.

p_ctrl, p_var = conv_ctrl / n_ctrl, conv_var / n_var   # 0.032, 0.036

abs_lift = p_var - p_ctrl                # 0.004  -> 0.4pp
rel_lift = abs_lift / p_ctrl             # 0.125  -> +12.5%
print(f"absolute lift = {abs_lift*100:.2f}pp")
print(f"relative lift = {rel_lift*100:.1f}%")

rel_low, rel_high = low / p_ctrl, high / p_ctrl
print(f"relative lift 95% CI = [{rel_low*100:.1f}%, {rel_high*100:.1f}%]")
# ≈ [3.7%, 21.3%]  ->  point estimate +12.5%, but the true lift
# could plausibly be as small as ~4% or as large as ~21%

The product team did not just watch conversion. They eyeballed AOV, add-to-cart, refund rate, and page latency too - and refund rate "improved" at p = 0.048. We correct for the five simultaneous tests with multipletests and watch that naive win disappear.

from statsmodels.stats.multitest import multipletests

metrics = ["conversion", "add_to_cart", "AOV", "refund_rate", "latency"]
pvals   = [0.005,        0.009,         0.028,  0.048,         0.55]
m = len(pvals)

naive = [p < 0.05 for p in pvals]
print("naive winners (α=0.05):", [x for x, w in zip(metrics, naive) if w])
# ['conversion', 'add_to_cart', 'AOV', 'refund_rate'] - 4 "wins"

bonf = multipletests(pvals, alpha=0.05, method='bonferroni')[0]
bh   = multipletests(pvals, alpha=0.05, method='fdr_bh')[0]

print("Bonferroni survivors:", [x for x, s in zip(metrics, bonf) if s])
# ['conversion', 'add_to_cart']   (threshold 0.05/5 = 0.01)
print("Benjamini-Hochberg :",     [x for x, s in zip(metrics, bh)   if s])
# ['conversion', 'add_to_cart', 'AOV']   (more powerful - keeps AOV)

• Bootstrap the CI to cross-check. Resample each arm with replacement 10,000 times, recompute the difference in rates each time, and take the 2.5th / 97.5th percentiles. Confirm the bootstrap interval lands close to the analytic confint_proportions_2indep result - two independent methods agreeing is your sanity check.
• Re-run the sweep with more metrics. Add device and geo splits until you have ~15 secondary metrics, feed random-ish p-values, and watch how many "winners" Bonferroni vs BH keep. Feel the power gap widen as m grows.
• Flip absolute and relative. Recompute the lift assuming the baseline were 8% instead of 3.2%. Same 0.4pp absolute lift, but the relative lift shrinks to +5% - see why baseline choice changes the headline.
• Sequential teaser. You analyzed this test once, at the pre-committed N of ~32,000/arm. What if a PM had peeked on day 3 and called it early? That is the peeking problem - fixed-horizon p-values are only valid at the planned N, and repeated looks can push Type I error past 30%. We fix it properly in b6 (sequential testing and peeking).

Analyzing a two-arm test is the working core of every A/B course. We teach ~80% of the practical readout here - the test, the interval, effect size, and multiple comparisons; the fuller inference theory and platform-specific stats engines live in the sources below.