5.1 Experiments

In an experiment we decide who gets what, and deciding by chance makes the groups comparable. Good design then makes the comparison as precise as possible: compare like with like, measure what varies, and plan the size before starting.

Code: data and style used in this chapter 
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats

from statly import style

style.use()
rng = np.random.default_rng(1)


def spread_of_estimates(ax: Axes, estimates: dict[str, np.ndarray], truth: float) -> None:
    """Histograms of an effect estimate under several designs, with their standard deviations."""
    bins = np.linspace(
        min(e.min() for e in estimates.values()), max(e.max() for e in estimates.values()), 40
    )
    for (name, values), color in zip(estimates.items(), (style.GREY, style.BLUE, style.GREEN)):
        ax.hist(
            values, bins=bins, color=color, alpha=0.5, label=f"{name}: sd {values.std():.2f}"
        )
    ax.axvline(truth, color=style.CORAL, ls="--", lw=1.5)
    ax.legend(fontsize=8)
    ax.set(xlabel="estimated effect (dashed: truth)", ylabel="repeated experiments")

5.1.1 Random assignment balances everything

In an experiment the researcher assigns the treatment. If the assignment is random, nothing about the units, measured or unmeasured, can influence who is treated, so the groups differ only by chance (Chapter 5.4). Chance imbalances still happen, and they are larger in small experiments:

Code 
SIZES = (20, 100, 1000)


def age_gaps(n: int, repeats: int = 2000) -> np.ndarray:
    """Difference in mean age between two randomly assigned halves of n people."""
    gaps = np.empty(repeats)
    for r in range(repeats):
        ages = rng.normal(50, 15, n)
        treated = rng.permutation(n) < n // 2
        gaps[r] = ages[treated].mean() - ages[~treated].mean()
    return gaps


fig, ax = plt.subplots(figsize=(8, 3.2))
bins = np.linspace(-20, 20, 81)
for n, color in zip(SIZES, (style.CORAL, style.GOLD, style.BLUE)):
    ax.hist(age_gaps(n), bins=bins, color=color, alpha=0.5, label=f"{n} people")
ax.legend(fontsize=8)
ax.set(xlabel="difference in average age between groups (years)", ylabel="experiments")
plt.show()
Figure 5.1.1: Random assignment of people to two equal groups, repeated 2,000 times for experiments of three sizes. Shown is the difference in average age between the groups (people’s ages have standard deviation 15 years). Randomisation balances age on average in every size; in a trial of 20 people, a difference of 10 years is not rare; in a trial of 1,000, it never happens. The same holds for every other characteristic, including ones nobody measured.

How to randomise matters less than that it is done properly: by a computer or a table of random numbers, decided before anyone sees the units, and hidden from those who enrol them, so that no one can steer the sicker patients or the better plots to one group (allocation concealment).

5.1.2 Controls, placebos and blinding

Without a control group, a before-and-after comparison credits the treatment with everything else that changed: the passage of time, the season, and regression to the mean (Chapter 4.1). Patients often enrol in a trial when their symptoms are at their worst, and would improve anyway.

Code 
n = 4000
own_average = rng.normal(145, 12, n)
screening = own_average + rng.normal(0, 10, n)
enrolled = screening > 160
later = own_average + rng.normal(0, 10, n)  # the treatment changes nothing
arm = rng.permutation(enrolled.sum()) < enrolled.sum() // 2

fig, axes = style.subplots(2, height=3.2, width=10)
for i in np.flatnonzero(enrolled)[:150]:
    axes[0].plot([0, 1], [screening[i], later[i]], color=style.GREY, lw=0.4, alpha=0.6)
axes[0].plot(
    [0, 1], [screening[enrolled].mean(), later[enrolled].mean()], color=style.CORAL, lw=3
)
style.note(
    axes[0],
    f"average change {later[enrolled].mean() - screening[enrolled].mean():+.1f}",
    x=0.4,
    y=0.95,
    color=style.CORAL,
)
axes[0].set_xticks([0, 1], ["screening", "after treatment"])
axes[0].set(ylabel="blood pressure", title="No control group")
for group, color, name in ((arm, style.CORAL, "treatment"), (~arm, style.BLUE, "placebo")):
    before, after = screening[enrolled][group], later[enrolled][group]
    axes[1].plot(
        [0, 1],
        [before.mean(), after.mean()],
        "o-",
        color=color,
        lw=2.5,
        label=f"{name}: {after.mean() - before.mean():+.1f}",
    )
axes[1].set_xticks([0, 1], ["screening", "later"])
axes[1].legend(fontsize=8)
axes[1].set(ylabel="average blood pressure", title="Randomised control group")
fig.tight_layout()
plt.show()
Figure 5.1.2: Simulated patients whose blood pressure varies from day to day around their own average. Those measuring above 160 on the screening day are enrolled and given a treatment that does nothing. Left: before and after, every enrolled patient’s measurement (grey lines) and the average (red). Blood pressure falls by about 10, with no treatment effect at all. Right: with a randomised control group given a placebo, both groups fall alike, and the difference between them is about zero.

A placebo makes the control group’s experience as similar as possible to the treated group’s, so that expectations act on both. Blinding keeps participants, and if possible those who treat and measure them, unaware of who got what, because knowing changes behaviour and judgement. A double-blind trial hides it from both.

5.1.3 Blocking: compare like with like

Units often differ in ways that affect the outcome a lot: plots of land in fertility, patients in severity, machines in age. Randomisation keeps such differences from biasing the comparison, but they still add noise. Blocking groups similar units together (neighbouring plots, patients of similar severity) and randomises within each block, so that the treatment is compared within blocks, where the units are alike. The extreme case is the paired design (Chapter 3.4), and the crossover design in which each person receives both treatments in random order.

Code 
fertility = np.linspace(0, 20, 40)


def field_experiment(blocked: bool) -> tuple[float, np.ndarray]:
    """One experiment: the estimated effect and which plots were fertilised."""
    if blocked:
        fertilised = np.ravel([rng.permutation([True, False]) for _ in range(20)])
    else:
        fertilised = rng.permutation(40) < 20
    yield_ = 50 + fertility + 2 * fertilised + rng.normal(0, 2, 40)
    if blocked:
        pairs = yield_.reshape(20, 2)
        sign = np.where(fertilised.reshape(20, 2)[:, 0], 1, -1)
        return float(np.mean(sign * (pairs[:, 0] - pairs[:, 1]))), fertilised
    return float(yield_[fertilised].mean() - yield_[~fertilised].mean()), fertilised


fig = plt.figure(figsize=(10.5, 3.3))
grid = fig.add_gridspec(1, 2, width_ratios=[1.1, 1])
ax = fig.add_subplot(grid[0])
for row, blocked in ((1, False), (0, True)):
    _, fertilised = field_experiment(blocked)
    for i in range(40):
        shade = plt.cm.Greens(0.2 + 0.6 * i / 40)
        ax.add_patch(plt.Rectangle((i, row), 0.9, 0.7, color=shade))
        if fertilised[i]:
            ax.scatter(i + 0.45, row + 0.35, s=14, color=style.CORAL, zorder=3)
    ax.text(
        -1,
        row + 0.35,
        "blocks of 2" if blocked else "completely\nrandom",
        ha="right",
        va="center",
        fontsize=8.5,
    )
ax.set(xlim=(-9, 40.5), ylim=(-0.3, 2))
ax.set_axis_off()
ax.set_title("fertility rises to the right; dots: fertilised", fontsize=9)
spread_of_estimates(
    fig.add_subplot(grid[1]),
    {
        "completely random": np.array([field_experiment(False)[0] for _ in range(5000)]),
        "blocked": np.array([field_experiment(True)[0] for _ in range(5000)]),
    },
    2,
)
fig.tight_layout()
plt.show()
Figure 5.1.3: A field of 40 plots with a strong fertility gradient from left to right (shading), testing a fertiliser that adds 2 units of yield. Left: completely random assignment (top) against blocks of two neighbouring plots, one of each treatment per block (bottom; filled: fertilised). Right: the estimated effect in 5,000 repetitions of each design. Both are unbiased, but blocking, which compares neighbours, is far more precise.

“Block what you can, randomise what you cannot.”

5.1.4 Adjusting for a baseline

Measuring the outcome before treatment as well as after, and including the baseline as a predictor in the analysis (an analysis of covariance, or ANCOVA), removes much of the person-to-person variation from the comparison. It is free precision, as long as the baseline is measured before randomisation.

Code 
def one_trial(n: int = 60) -> tuple[float, float, float]:
    """Effect estimates from follow-up only, from change scores, and from ANCOVA."""
    baseline = rng.normal(50, 10, n)
    treated = rng.permutation(n) < n // 2
    follow_up = 15 + 0.7 * baseline + 3 * treated + rng.normal(0, 10 * np.sqrt(1 - 0.7**2), n)
    change = follow_up - baseline
    ancova = (
        sm.OLS(follow_up, sm.add_constant(np.column_stack([treated, baseline]))).fit().params[1]
    )
    return (
        follow_up[treated].mean() - follow_up[~treated].mean(),
        change[treated].mean() - change[~treated].mean(),
        ancova,
    )


trials = np.array([one_trial() for _ in range(5000)])
fig, ax = plt.subplots(figsize=(8, 3.2))
spread_of_estimates(
    ax,
    {
        "follow-up only": trials[:, 0],
        "change from baseline": trials[:, 1],
        "ANCOVA": trials[:, 2],
    },
    3,
)
plt.show()
Figure 5.1.4: A trial of 60 patients where the follow-up score correlates 0.7 with the baseline score and the treatment adds 3 points. Estimated effects in 5,000 simulated trials, from three analyses: comparing follow-up scores, comparing changes from baseline, and regression of follow-up on treatment and baseline (ANCOVA). All are unbiased; ANCOVA is the most precise.

5.1.5 Several factors at once: factorial designs

Testing two factors one at a time needs two experiments. A factorial design crosses them: every combination of the levels is tried. A 2 × 2 design answers both questions with the same units, each effect being estimated from all of them, and it is the only way to see whether the factors interact (Chapter 3.5): whether the effect of one depends on the level of the other.

Code 
CELL = 25
water = np.repeat([0, 0, 1, 1], CELL)
fertiliser = np.repeat([0, 1, 0, 1], CELL)
growth = 10 + 3 * water + 2 * fertiliser + 3 * water * fertiliser + rng.normal(0, 2.5, 4 * CELL)
table = sm.OLS(
    growth, sm.add_constant(np.column_stack([water, fertiliser, water * fertiliser]))
).fit()

fig, axes = style.subplots(2, height=3.2, width=10)
for f, color, name in ((0, style.GREY, "no fertiliser"), (1, style.GREEN, "fertiliser")):
    means = [growth[(water == w) & (fertiliser == f)].mean() for w in (0, 1)]
    errors = [growth[(water == w) & (fertiliser == f)].std() / np.sqrt(CELL) for w in (0, 1)]
    axes[0].errorbar(
        [0, 1], means, yerr=errors, marker="o", color=color, lw=2, capsize=4, label=name
    )
axes[0].set_xticks([0, 1], ["low water", "high water"])
axes[0].legend(fontsize=8)
style.note(
    axes[0],
    f"interaction: {table.params[3]:+.1f} (p = {table.pvalues[3]:.3f})",
    x=0.03,
    y=0.95,
    color=style.INK,
)
axes[0].set(ylabel="growth", title="Cell means")
for i, (w, f) in enumerate(((0, 0), (0, 1), (1, 0), (1, 1))):
    axes[1].add_patch(
        plt.Rectangle(
            (w, f), 0.95, 0.95, color=style.GREEN if f else style.GREY, alpha=0.25 + 0.3 * w
        )
    )
    axes[1].text(
        w + 0.47,
        f + 0.47,
        f"{CELL} units\nwater {'high' if w else 'low'}\n{'fertiliser' if f else 'none'}",
        ha="center",
        va="center",
        fontsize=8.5,
    )
axes[1].set(xlim=(-0.1, 2.05), ylim=(-0.1, 2.05), title="The 2 × 2 design")
axes[1].set_axis_off()
fig.tight_layout()
plt.show()
Figure 5.1.5: A 2 × 2 factorial experiment (simulated, 25 units per cell): watering (low, high) crossed with fertiliser (none, added). Left: the four cell means with their standard errors. Fertiliser helps more when water is plentiful: the lines are not parallel, an interaction. Right: the design. Each main effect compares two whole halves of the experiment, so every unit contributes to every estimate.

With many factors, a full factorial grows quickly (2k2^k cells for k two-level factors). Fractional factorial designs run a carefully chosen fraction of the cells, giving up some high-order interactions that are usually small.

5.1.6 A/B tests

Online, experiments are cheap: visitors are randomly shown version A or B of a page, and the share who buy, click or stay is compared. The statistics are the two-share comparison of Chapter 3.4, and the size needed follows from the power calculation of Chapter 3.10. For two shares pAp_A and pBp_B, each group needs about

n≈(z1−α/2+z1−β)2[pA(1−pA)+pB(1−pB)](pA−pB)2 n \approx \frac{(z_{1-\alpha/2} + z_{1-\beta})^2 \, [p_A(1-p_A) + p_B(1-p_B)]}{(p_A - p_B)^2}

visitors. Small effects on small rates need enormous samples:

Code 
def visitors_needed(
    base: float, lift: np.ndarray, power: float = 0.8, alpha: float = 0.05
) -> np.ndarray:
    """Per-group sample size for a two-share comparison."""
    p_b = base * (1 + lift)
    z = stats.norm.ppf(1 - alpha / 2) + stats.norm.ppf(power)
    return z**2 * (base * (1 - base) + p_b * (1 - p_b)) / (p_b - base) ** 2


lifts = np.linspace(0.02, 0.3, 100)
days = np.arange(1, 31)
daily = 2000
conv_a = np.cumsum(rng.binomial(daily, 0.05, 30))
conv_b = np.cumsum(rng.binomial(daily, 0.055, 30))
seen = days * daily
diff = conv_b / seen - conv_a / seen
se = np.sqrt(
    (conv_a / seen) * (1 - conv_a / seen) / seen + (conv_b / seen) * (1 - conv_b / seen) / seen
)

fig, axes = style.subplots(2, height=3.2, width=10.5)
for base, color in ((0.02, style.CORAL), (0.05, style.GOLD), (0.2, style.BLUE)):
    axes[0].plot(
        lifts, visitors_needed(base, lifts), color=color, lw=2, label=f"base rate {base:.0%}"
    )
axes[0].set(
    yscale="log", xlabel="relative improvement", ylabel="visitors per group (log scale)"
)
axes[0].xaxis.set_major_formatter(style.PERCENT)
axes[0].legend(fontsize=8)
axes[1].fill_between(days, diff - 1.96 * se, diff + 1.96 * se, color=style.BLUE, alpha=0.25)
axes[1].plot(days, diff, color=style.BLUE, lw=2)
axes[1].axhline(0, color=style.INK, lw=0.8)
axes[1].axhline(0.005, color=style.CORAL, ls="--", lw=1)
axes[1].yaxis.set_major_formatter(style.PERCENT)
axes[1].set(xlabel="day", ylabel="B − A conversion (dashed: truth)")
fig.tight_layout()
plt.show()
Figure 5.1.6: Left: visitors needed per group to detect an improvement in a conversion rate with 80% power at the 5% level, for base rates of 2%, 5% and 20%, as a function of the relative improvement. Lifting a 2% rate by 5% (to 2.1%) needs about 315,000 visitors per group. Right: a simulated A/B test, B truly 10% better than A’s 5%. The difference in conversion and its 95% interval, day by day. The interval narrows as visitors accumulate; stopping the first time it excludes zero, by peeking every day, inflates false positives (Chapter 3.10).

Practical traps of A/B tests: deciding the size and the metric in advance; checking that the two groups really got the planned shares of visitors (a sample ratio mismatch signals a broken assignment); novelty effects that fade after a few days; and running many variants or metrics at once (Chapter 3.10).

5.1.7 In one picture

  • Random assignment balances measured and unmeasured characteristics on average; chance imbalance shrinks with size.
  • A control group separates the treatment from time, placebo effects and regression to the mean; blinding keeps expectations out.
  • Blocking and pairing compare like with like; baseline adjustment (ANCOVA) adds precision for free.
  • Factorial designs test several factors with the same units and reveal interactions.
  • A/B tests are experiments: fix size and metric in advance, check the assignment, do not peek.

Next: surveys, where the question is not what a treatment does but what a population looks like.