3.4 Comparing two groups

A difference between two groups means something only next to the noise around it. Every two-group test divides the signal by the noise; they differ in how they measure each.

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)

ccard = sm.datasets.ccard.load_pandas().data  # credit card spending of 72 cardholders
owners = ccard.AVGEXP[ccard.OWNRENT == 1].to_numpy()  # average monthly spending, home owners
renters = ccard.AVGEXP[ccard.OWNRENT == 0].to_numpy()  # and renters

fertility = sm.datasets.fertility.load_pandas().data  # births per woman, by country and year
countries = fertility.loc[
    ~fertility["Country Name"].str.contains(
        "income|members|small states"
    ),  # drop region totals
    ["Country Name", "1990", "2010"],
].dropna()
before, after = countries["1990"].to_numpy(), countries["2010"].to_numpy()


def plot_groups(ax: Axes, groups: list[tuple[str, np.ndarray, str]]) -> None:
    """Every value as a dot, spread sideways a little, with each group's average."""
    for i, (name, values, color) in enumerate(groups):
        ax.scatter(
            i + rng.uniform(-0.15, 0.15, len(values)), values, s=14, color=color, alpha=0.6
        )
        ax.plot([i - 0.28, i + 0.28], [values.mean()] * 2, color=color, lw=2.5)
    ax.set_xticks(range(len(groups)), [name for name, *_ in groups])
    ax.set_xlim(-0.6, len(groups) - 0.4)


def plot_t(ax: Axes, t: float, df: float, title: str = "If there were no difference") -> None:
    """The t distribution under the null hypothesis, with both tails beyond the observed t shaded."""
    limit = max(4.5, abs(t) + 1)
    x = np.linspace(-limit, limit, 400)
    curve = stats.t.pdf(x, df)
    ax.plot(x, curve, color=style.BLUE)
    for tail in (x >= abs(t), x <= -abs(t)):
        ax.fill_between(x[tail], curve[tail], color=style.CORAL, alpha=0.5, lw=0)
    ax.axvline(t, color=style.CORAL, lw=1.6)
    p = 2 * stats.t.sf(abs(t), df)
    style.note(
        ax, f"t = {t:.2f}\np = {p:.3f}" if p >= 0.001 else f"t = {t:.2f}\np < 0.001", x=0.03
    )
    ax.set(title=title, xlabel="t", ylabel="how common", yticks=[])

3.4.1 One group against a fixed value

The simplest comparison is one group against a number. About 2.1 births per woman keep a population stable in the long run, the replacement level. Across 194 countries in 2010, is the average above it?

The one-sample t-test divides the distance between the sample mean and the fixed value μ0\mu_0 by the standard error of the mean:

t=x‾−μ0s/n t = \frac{\bar x - \mu_0}{s / \sqrt{n}}

With no real difference, t follows a t distribution with n − 1 degrees of freedom.

Code 
one_sample = stats.ttest_1samp(after, 2.1)


def plot_values(ax: Axes) -> None:
    """The countries' values, the fixed value and the mean."""
    style.hist(ax, after, bins=30, color=style.GREEN)
    ax.axvline(2.1, color=style.INK, ls="--", lw=1.4)
    ax.axvline(after.mean(), color=style.CORAL, lw=1.8)
    style.note(ax, f"mean {after.mean():.2f}\nmedian {np.median(after):.2f}", x=0.6)
    ax.set(title="Countries in 2010", xlabel="births per woman", ylabel="countries")


fig, axes = style.subplots(2, height=3.1, width=9.5)
plot_values(axes[0])
plot_t(axes[1], one_sample.statistic, one_sample.df, "If the average were 2.1")
fig.tight_layout()
plt.show()
Figure 3.4.1: Left: births per woman in 194 countries in 2010, with the replacement level of 2.1 (dashed) and the countries’ average (red). Right: the observed t against the t distribution it would follow if the true average were 2.1.

The average is 2.92, t = 7.6: far above replacement, as a world average. But note the shape: most countries are near 2, with a long tail of countries with 5 or more. The median, 2.39, is much closer to 2.1. A t-test answers a question about the mean; whether that is the right question depends on what you want to know.

3.4.2 Two groups: the t-test

Seventy-two credit card holders: 27 own their home and spent on average 362 a month on the card, 45 rent and spent 203. Is the difference real, or could it be chance?

The two-sample t-test answers with a ratio: the difference between the averages, divided by its standard error. The standard errors of the two means combine like the sides of a right triangle (Chapter 2.1):

t=x‾1−x‾2s12/n1+s22/n2 t = \frac{\bar x_1 - \bar x_2}{\sqrt{s_1^2 / n_1 + s_2^2 / n_2}}

A large t means the difference stands out from the noise. With no real difference, t follows a t distribution, which gives the p-value.

Code 
welch = stats.ttest_ind(owners, renters, equal_var=False)
student = stats.ttest_ind(owners, renters, equal_var=True)

fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_groups(axes[0], [("renters", renters, style.BLUE), ("owners", owners, style.CORAL)])
axes[0].set(title="Spending", ylabel="average monthly card spending")
plot_t(axes[1], welch.statistic, welch.df)
fig.tight_layout()
plt.show()
Figure 3.4.2: Left: monthly spending of every cardholder, with each group’s average (line). Right: the t distribution for no difference, with the observed t (red line); the red tails hold the p-value.

t = 2.1 and p = 0.04: a difference this large would be rare if owners and renters spent the same. The version above is Welch’s t-test, which lets the two groups have different spreads. Its degrees of freedom come from a formula (Welch and Satterthwaite’s) and need not be whole: here about 54. The older Student’s t-test assumes equal spreads and pools them into one estimate,

sp2=(n1−1)s12+(n2−1)s22n1+n2−2,t=x‾1−x‾2sp1/n1+1/n2 s_p^2 = \frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 - 2}, \qquad t = \frac{\bar x_1 - \bar x_2}{s_p\sqrt{1/n_1 + 1/n_2}}

with n1+n2−2n_1 + n_2 - 2 degrees of freedom. Here both give t = 2.1. When spreads and group sizes differ, Student’s test can be badly wrong, while Welch’s is never much worse, so use Welch’s by default.

3.4.3 How big is the difference?

A p-value mixes the size of an effect with the size of the sample. To describe the size alone, divide the difference by the spread of the values (not of the averages). That is Cohen’s d:

d=x‾1−x‾2sp d = \frac{\bar x_1 - \bar x_2}{s_p}

how many standard deviations apart the groups are.

Code 
x = np.linspace(-4, 5, 400)


def plot_effect(ax: Axes, d: float, label: str) -> None:
    """Two normal curves d standard deviations apart, with their overlap."""
    a, b = stats.norm.pdf(x), stats.norm.pdf(x, d)
    ax.fill_between(x, a, color=style.BLUE, alpha=0.3, lw=0)
    ax.fill_between(x, b, color=style.CORAL, alpha=0.3, lw=0)
    ax.plot(x, a, color=style.BLUE, lw=1.5)
    ax.plot(x, b, color=style.CORAL, lw=1.5)
    above = stats.norm.cdf(d / np.sqrt(2))  # chance a random member of B beats one of A
    style.note(ax, f"B beats A\n{above:.0%} of the time", x=0.66, color=style.INK)
    ax.set(title=f"{label}: d = {d}", xlabel="standard deviations", yticks=[])


fig, axes = style.subplots(3, height=2.6, width=10.5)
for ax, (d, label) in zip(axes, ((0.2, "Small"), (0.5, "Medium"), (0.8, "Large"))):
    plot_effect(ax, d, label)
axes[0].set_ylabel("how common")
fig.tight_layout()
plt.show()
Figure 3.4.3: What effect sizes look like: two normal groups, A (blue) and B (red), d standard deviations apart. Even a large effect (d = 0.8) leaves the groups overlapping heavily.

The owners and renters are about half a standard deviation apart, a medium effect. The labels small, medium and large are rough guides from psychology; in your own field, compare with effects people care about. And note how much the groups overlap even at d = 0.8: a difference between averages says little about any one person.

3.4.4 When the data are skewed: the Mann-Whitney test

The spending data have long right tails: a few cardholders spent far more than the rest. Those few move the averages, and the t-test leans on the averages. The Mann-Whitney test (also called the Wilcoxon rank-sum test) ignores the sizes and asks only about the order. Take every pair of one owner and one renter, 27 × 45 = 1,215 pairs, and count the pairs in which the owner spent more. That count is the statistic U.

Code 
mann_whitney = stats.mannwhitneyu(owners, renters)
wins = np.sort(owners)[:, None] > np.sort(renters)[None, :]
everyone = np.concatenate([renters, owners])
ranks = stats.rankdata(everyone)


def plot_pairs(ax: Axes) -> None:
    """The grid of owner-renter pairs, red where the owner spent more."""
    ax.imshow(
        wins,
        cmap=plt.matplotlib.colors.ListedColormap([style.BLUE, style.CORAL]),
        alpha=0.6,
        aspect="auto",
        origin="lower",
    )
    style.note(
        ax,
        f"owner higher in {wins.mean():.0%}\nof {wins.size:,} pairs",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.set(
        title="All pairs",
        xlabel="renters, lowest to highest",
        ylabel="owners, lowest to highest",
        xticks=[],
        yticks=[],
    )
    ax.grid(False)


fig, axes = style.subplots(2, height=3.3, width=9.5)
plot_pairs(axes[0])
plot_groups(
    axes[1],
    [
        ("renters", ranks[: len(renters)], style.BLUE),
        ("owners", ranks[len(renters) :], style.CORAL),
    ],
)
axes[1].set(title="Ranks", ylabel="rank (1 = lowest spending)")
style.note(
    axes[1],
    f"t-test: p = {welch.pvalue:.3f}\nMann-Whitney: p = {mann_whitney.pvalue:.3f}",
    x=0.03,
    y=0.98,
)
fig.tight_layout()
plt.show()
Figure 3.4.4: Left: every owner (rows) against every renter (columns), both sorted by spending. A cell is red when the owner spent more. The red share, 70%, is the chance that a random owner outspends a random renter; with no difference it would be about 50%. Right: the same data as ranks, 1 for the lowest spender of all 72.

Where does the p-value come from? If owning or renting had nothing to do with spending, the labels “owner” and “renter” could be shuffled among the 72 cardholders without changing anything. Shuffling them thousands of times and recomputing U each time shows how U varies when nothing is going on:

Code 
n_pairs = len(owners) * len(renters)
observed_u = float(wins.sum())


def shuffled_u() -> float:
    """U after shuffling the owner/renter labels among all cardholders."""
    order = rng.permutation(everyone)
    shuffled_owners, shuffled_renters = order[len(renters) :], order[: len(renters)]
    return float((shuffled_owners[:, None] > shuffled_renters[None, :]).sum())


null_u = np.array([shuffled_u() for _ in range(10_000)])
spread = np.sqrt(len(owners) * len(renters) * (len(everyone) + 1) / 12)
fig, ax = plt.subplots(figsize=(8, 3.2))
style.null_panel(
    ax,
    observed_u,
    samples=null_u,
    dist=stats.norm(n_pairs / 2, spread),
    tail="both",
    centre=n_pairs / 2,
    xlabel="U: pairs in which the owner spent more",
)
plt.show()
Figure 3.4.5: The Mann-Whitney statistic U (pairs in which the owner spent more) for 10,000 random shufflings of the owner and renter labels (grey bars), and the normal curve that approximates it (blue), centred on half of the 1,215 pairs. The observed U lies far out; the red areas in both tails, as far from the centre as the observed value or further, are the two-sided p-value.

U/(n1n2)U / (n_1 n_2) is the probability of superiority: here 0.70, an effect size that needs no normal shape. The test finds the difference more clearly than the t-test, p = 0.004 against 0.04, because the extreme spenders no longer add noise. Taking logs first, as in Chapter 2.1, helps the t-test in a similar way.

The price of ranks: the test no longer speaks about averages in money, only about which group tends to be higher. When the data are close to normal, the t-test is a little more powerful; when they have long tails, the rank test often wins by a lot.

3.4.5 Comparing whole distributions

The t-test compares means and the Mann-Whitney test compares which group tends to be higher. Two groups can also differ in spread or shape with the same centre. The two-sample Kolmogorov-Smirnov test compares the whole cumulative curves (Chapter 2.3) and uses the largest vertical gap between them:

Code 
ks = stats.ks_2samp(owners, renters)


def plot_staircases(ax: Axes) -> None:
    """Both groups' cumulative staircases, with the largest gap."""
    for values, color, name in (
        (renters, style.BLUE, "renters"),
        (owners, style.CORAL, "owners"),
    ):
        v = np.sort(values)
        ax.step(
            np.r_[0, v],
            np.r_[0, np.arange(1, len(v) + 1) / len(v)],
            where="post",
            color=color,
            lw=1.8,
            label=name,
        )
    at = ks.statistic_location
    low, high = (
        np.searchsorted(np.sort(g), at, side="right") / len(g) for g in (owners, renters)
    )
    ax.annotate(
        "",
        xy=(at, low),
        xytext=(at, high),
        arrowprops=dict(arrowstyle="<->", color=style.INK, lw=1.5),
    )
    style.note(
        ax, f"D = {ks.statistic:.2f}, p = {ks.pvalue:.3f}", x=0.45, y=0.3, color=style.INK
    )
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8, loc="lower right")
    ax.set(xlabel="average monthly card spending", ylabel="share spending less", xlim=(0, 1000))


def shuffled_d() -> float:
    """The largest gap between the two staircases after shuffling the group labels."""
    order = rng.permutation(everyone)
    return float(stats.ks_2samp(order[len(renters) :], order[: len(renters)]).statistic)


null_d = np.array([shuffled_d() for _ in range(5000)])
fig, axes = style.subplots(2, height=3.3, width=10.5)
plot_staircases(axes[0])
style.null_panel(axes[1], ks.statistic, samples=null_d, tail="right", xlabel="largest gap D")
axes[1].set_title("D when the labels are shuffled")
fig.tight_layout()
plt.show()
Figure 3.4.6: Left: the cumulative curves of spending for renters (blue) and owners (red). The arrow marks the largest vertical gap, D = 0.44: at a spending of about 230 a month, 78% of renters but only 33% of owners spent less. Right: D for 5,000 shufflings of the owner and renter labels. Shuffled groups also differ by chance, and D is never zero, but a gap as large as 0.44 is rare: the share of shuffles at least as large is the p-value. Only large gaps count as evidence, so the test is one-sided.

It is the most general of the three, and for that reason the least powerful for a plain shift in the centre.

3.4.6 Paired data

Sometimes the two groups are the same units measured twice: patients before and after treatment, countries in two years. Then compare each unit with itself. The differences between units, often large, drop out.

Code 
change = after - before


def plot_lines(ax: Axes) -> None:
    """One line per country, from 1990 to 2010."""
    for b, a in zip(before, after):
        ax.plot([0, 1], [b, a], color=style.CORAL if a > b else style.BLUE, lw=0.7, alpha=0.4)
    ax.set_xticks([0, 1], ["1990", "2010"])
    ax.set(title="Each country", ylabel="births per woman", xlim=(-0.2, 1.2))


def plot_changes(ax: Axes) -> None:
    """The change per country, with the average change and its interval."""
    style.hist(ax, change, bins=30, color=style.GREEN)
    ax.axvline(0, color=style.INK, lw=1)
    mean = change.mean()
    half = stats.t.ppf(0.975, len(change) - 1) * change.std(ddof=1) / np.sqrt(len(change))
    style.note(
        ax,
        f"average change {mean:.2f}\n95% interval {mean - half:.2f} to {mean + half:.2f}",
        x=0.03,
    )
    ax.set(title="Change per country", xlabel="change in births per woman", ylabel="countries")


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_lines(axes[0])
plot_changes(axes[1])
fig.tight_layout()
plt.show()
Figure 3.4.7: Births per woman in 194 countries, 1990 and 2010. Left: one line per country; most go down. Right: the change for each country. Countries differ a lot from each other, but the changes are much more alike, and that makes the paired test far sharper.

Births per woman fell by about one child in twenty years. A paired t-test is simply the one-sample t-test of Section 3.4.1 on the changes, against 0. The standard deviation of the changes (0.87) is much smaller than the spread between countries (1.9 in 1990), so the paired test gives t = −17, against t = −6 if we wrongly treated the two years as separate groups. Pairing removes the noise of countries being different.

Two rank-free alternatives for pairs:

  • The sign test looks only at the direction of each change. 176 of the 194 countries went down. If ups and downs were equally likely, that would be like 176 heads in 194 tosses: the binomial test of Chapter 3.1 gives p < 10⁻³⁰.
  • The Wilcoxon signed-rank test ranks the sizes of the changes, ignoring their sign, and then adds up the ranks of the increases. Small sum, mostly decreases. It uses more information than the sign test and, like Mann-Whitney, does not care about extreme values.

Both have a simple picture of “no change”: each country’s change would be equally likely to go up or down. So flip a fair coin for the sign of every change, recompute the statistic, and repeat:

Code 
down = int((change < 0).sum())
signed_ranks = stats.rankdata(np.abs(change))
observed_w = float(signed_ranks[change > 0].sum())
null_w = np.array(
    [signed_ranks[rng.uniform(size=len(change)) < 0.5].sum() for _ in range(10_000)]
)
n_changes = len(change)
w_mean = n_changes * (n_changes + 1) / 4
w_sd = np.sqrt(n_changes * (n_changes + 1) * (2 * n_changes + 1) / 24)

fig, axes = style.subplots(2, height=3.2, width=10.5)
counts = np.arange(n_changes + 1)
probability = stats.binom.pmf(counts, n_changes, 0.5)
axes[0].bar(counts, probability, color=style.BLUE, alpha=0.7, width=1)
axes[0].axvline(down, color=style.CORAL, ls="--", lw=1.5)
style.note(axes[0], f"observed {down} of {n_changes}\np < 10⁻³⁰", x=0.97, y=0.95, ha="right")
axes[0].set(
    title="Sign test",
    xlabel="countries going down, if up and down were equally likely",
    ylabel="probability",
    yticks=[],
)
style.null_panel(
    axes[1],
    observed_w,
    samples=null_w,
    dist=stats.norm(w_mean, w_sd),
    tail="both",
    centre=w_mean,
    xlabel="sum of the ranks of the increases",
)
axes[1].set_title("Wilcoxon signed-rank test")
fig.tight_layout()
plt.show()
Figure 3.4.8: The null distributions of the two paired tests for the 194 countries. Left: the sign test. If ups and downs were equally likely, the number of countries going down would follow the binomial distribution (bars); the observed 176 lies so far to the right that the p-value is smaller than 10⁻³⁰. Right: the Wilcoxon signed-rank statistic, the sum of the ranks of the increases, for 10,000 sets of random signs (grey), with its normal approximation (blue). The observed sum, from only 18 increases that are also small ones, sits at the far left.

3.4.7 Two shares

When the outcome is yes or no, each group has a share. In the teaching experiment of Chapter 3.1, 8 of 14 students improved with the new method, 3 of 18 without. The two-proportion z-test divides the difference in shares by its standard error, using the pooled share p̂\hat p of both groups together (what the share would be if the groups did not differ):

z=p̂1−p̂2p̂(1−p̂)(1/n1+1/n2) z = \frac{\hat p_1 - \hat p_2}{\sqrt{\hat p(1 - \hat p)\,(1/n_1 + 1/n_2)}}

With small counts, the normal curve is a poor approximation. Fisher’s exact test avoids it. It fixes the margins of the 2 × 2 table, 11 improvers among 32 students, 14 of them taught the new way, and asks: if the method did nothing, how many of the 11 improvers would land in the new group by chance? That count follows a known distribution (the hypergeometric), and the p-value is the probability of tables at least as lopsided as the one observed.

Code 
table = np.array([[8, 6], [3, 15]])  # rows: new method, usual; columns: improved, not
fisher = stats.fisher_exact(table)
z_p = stats.chi2_contingency(table, correction=False).pvalue
k = np.arange(0, 12)
probability = stats.hypergeom.pmf(k, 32, 11, 14)
as_unlikely = probability <= probability[8] + 1e-12


def plot_table(ax: Axes) -> None:
    """The counts in a 2 × 2 grid."""
    for i, row in enumerate(table):
        for j, count in enumerate(row):
            ax.add_patch(plt.Rectangle((j, 1 - i), 1, 1, fill=False, ec=style.INK))
            ax.text(j + 0.5, 1.5 - i, str(count), ha="center", va="center", fontsize=18)
    ax.set_xticks([0.5, 1.5], ["improved", "did not"])
    ax.set_yticks([1.5, 0.5], ["new method", "usual"])
    ax.set(xlim=(0, 2), ylim=(0, 2), aspect="equal")
    ax.grid(False)
    ax.tick_params(length=0)


def plot_hypergeometric(ax: Axes) -> None:
    """The distribution of improvers in the new-method group, if the method did nothing."""
    ax.bar(k, probability, color=np.where(as_unlikely, style.CORAL, style.BLUE), alpha=0.7)
    style.note(
        ax,
        f"Fisher: p = {fisher.pvalue:.3f}\nz-test: p = {z_p:.3f}\nodds ratio {fisher.statistic:.1f}",
        x=0.03,
        color=style.INK,
    )
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set(
        title="If the method did nothing",
        xlabel="improvers among the 14 taught the new way",
        ylabel="probability",
    )


fig = plt.figure(figsize=(9.5, 3.2))
grid = fig.add_gridspec(1, 2, width_ratios=[1, 2])
plot_table(fig.add_subplot(grid[0]))
plot_hypergeometric(fig.add_subplot(grid[1]))
fig.tight_layout()
plt.show()
Figure 3.4.9: Left: the 2 × 2 table of the teaching experiment. Right: if the method did nothing, the number of the 11 improvers who would be among the 14 taught the new way (bars). The observed 8, and anything at least as unlikely (red), add up to Fisher’s p-value.

The z-test (which is the same as the chi-square test of Chapter 3.5 for a 2 × 2 table) gives 0.017, Fisher’s exact test 0.027. With counts this small, trust Fisher’s. The odds ratio, (8/6) / (3/15) = 6.7, compares the odds of improving in the two groups; it returns in logistic regression (Chapter 4.4).

3.4.8 Paired yes/no: McNemar’s test

Yes/no outcomes can be paired too: the same patients tested with two methods, the same voters asked before and after a debate. Only the pairs that disagree carry information about a difference. McNemar’s test compares the two kinds of disagreement.

Code 
paired = np.array([[60, 10], [24, 106]])  # rows: old test positive/negative; columns: new test
discordant = paired[0, 1] + paired[1, 0]
mcnemar = stats.binomtest(int(paired[1, 0]), int(discordant), 0.5)


def plot_paired_table(ax: Axes) -> None:
    """The paired table, with the disagreeing cells highlighted."""
    for i in range(2):
        for j in range(2):
            disagree = i != j
            ax.add_patch(
                plt.Rectangle(
                    (j, 1 - i),
                    1,
                    1,
                    color=style.GOLD if disagree else "none",
                    alpha=0.25,
                    ec=style.INK,
                )
            )
            ax.text(j + 0.5, 1.5 - i, str(paired[i, j]), ha="center", va="center", fontsize=18)
    ax.set_xticks([0.5, 1.5], ["new +", "new −"])
    ax.set_yticks([1.5, 0.5], ["old +", "old −"])
    ax.set(xlim=(0, 2), ylim=(0, 2), aspect="equal")
    ax.grid(False)
    ax.tick_params(length=0)


def plot_split(ax: Axes) -> None:
    """How the disagreeing pairs would split if the tests were equally sensitive."""
    k = np.arange(discordant + 1)
    probability = stats.binom.pmf(k, discordant, 0.5)
    extreme = abs(k - discordant / 2) >= abs(paired[1, 0] - discordant / 2)
    ax.bar(k, probability, color=np.where(extreme, style.CORAL, style.BLUE), alpha=0.7)
    style.note(ax, f"McNemar (exact): p = {mcnemar.pvalue:.3f}", x=0.03, color=style.INK)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set(
        title="If both tests were equal",
        xlabel="pairs flagged only by the new test, of 34",
        ylabel="probability",
    )


fig = plt.figure(figsize=(9.5, 3.2))
grid = fig.add_gridspec(1, 2, width_ratios=[1, 2])
plot_paired_table(fig.add_subplot(grid[0]))
plot_split(fig.add_subplot(grid[1]))
fig.tight_layout()
plt.show()
Figure 3.4.10: Simulated: 200 patients, each tested with an old and a new test. Left: the paired table. Most patients get the same result from both tests and say nothing about which is more sensitive. Right: of the 34 patients where the tests disagree, 24 were flagged only by the new test and 10 only by the old one. If the tests were equally likely to flag, the split would be like a coin toss (bars).

The 166 patients on whom the tests agree are left out entirely. It is the sign test of Section 3.4.6 again, on the disagreeing pairs.

3.4.9 Showing that there is no difference

A large p-value does not show that two groups are the same (Chapter 3.1). To show equivalence, decide first how large a difference would matter, the margin, and then show that the confidence interval lies entirely inside it. This is the two one-sided tests (TOST) procedure: one test that the difference is above −margin, one that it is below +margin.

Code 
STUDIES = [
    ("A: equivalent", 0.8, -2.5, 4.1),
    ("B: inconclusive", 1.5, -4.5, 7.5),
    ("C: different, but trivially", 2.4, 0.8, 4.0),
]
MARGIN = 5


def plot_equivalence(ax: Axes) -> None:
    """Intervals against zero and against the equivalence margin."""
    ax.axvspan(-MARGIN, MARGIN, color=style.GREEN, alpha=0.1, lw=0)
    ax.axvline(0, color=style.INK, lw=1)
    for i, (name, estimate, low, high) in enumerate(STUDIES):
        inside = -MARGIN < low and high < MARGIN
        color = style.GREEN if inside else style.GREY
        ax.plot([low, high], [i, i], color=color, lw=2.5)
        ax.plot(estimate, i, "o", ms=7, color=color)
    ax.set_yticks(range(len(STUDIES)), [s[0] for s in STUDIES])
    ax.invert_yaxis()
    ax.set(xlabel="difference from the original (units)", xlim=(-8, 9))


fig, ax = plt.subplots(figsize=(8, 2.6))
plot_equivalence(ax)
plt.show()
Figure 3.4.11: Three imagined studies of a new generic drug against the original, with a margin of ±5 units (shaded) set in advance. A: the 90% interval lies inside the margin: equivalent. B: not significantly different, but the interval reaches beyond the margin: inconclusive. C: significantly different, yet inside the margin: different, but by too little to matter.

Two one-sided tests at 5% correspond to a 90% interval, which is why the figure uses 90%.

3.4.10 Which test when

Situation Test Without the normal shape
one group against a fixed value one-sample t-test Wilcoxon signed-rank, sign test
two separate groups, means Welch’s t-test Mann-Whitney
two separate groups, whole distributions Kolmogorov-Smirnov
the same units twice paired t-test Wilcoxon signed-rank, sign test
one share against a fixed value z-test for a share binomial test (exact)
two shares, separate groups two-proportion z-test Fisher’s exact test
two shares, paired McNemar’s test
showing no difference two one-sided tests (TOST)

And check the data first. Plot both groups. If there are long tails, use ranks or logs. If the units are paired, use the pairing. Chapter 3.13 puts all the tests of this part on one map.

3.4.11 In one picture

  • Two-group tests divide a difference by its standard error: signal over noise.
  • One sample: t=(x‾−μ0)/(s/n)t = (\bar x - \mu_0)/(s/\sqrt n); two samples: t=(x‾1−x‾2)/s12/n1+s22/n2t = (\bar x_1 - \bar x_2)/\sqrt{s_1^2/n_1 + s_2^2/n_2}.
  • Use Welch’s t-test by default; it allows different spreads.
  • Report an effect size: Cohen’s d, or the probability of superiority.
  • With long tails, compare ranks (Mann-Whitney); for shapes, Kolmogorov-Smirnov.
  • With paired data, test the changes: paired t, signed-rank, sign test; for yes/no, McNemar.
  • For shares, use the z-test with large counts and Fisher’s exact test with small ones.
  • To show no difference, show the interval lies inside a margin you set in advance.

Next: more than two groups, and tables of counts.