3.5 Many groups and tables of counts

With more than two groups, or with counts in a table, the test compares two kinds of variation: the variation the groups explain, and the variation left over. If the first is much larger than chance allows, the groups differ.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from matplotlib.axes import Axes
from scipy import stats
from sklearn import datasets
from statsmodels.stats.multicomp import pairwise_tukeyhsd

from statly import style

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

wine = datasets.load_wine(as_frame=True).frame  # 178 Italian wines from three cultivars
CULTIVARS = ["cultivar 1", "cultivar 2", "cultivar 3"]
CULTIVAR_COLORS = [style.BLUE, style.GREEN, style.CORAL]
fair = sm.datasets.fair.load_pandas().data  # 6,366 married women, survey from 1974
health = sm.datasets.randhie.load_pandas().data  # 20,190 people in the RAND health experiment
anes = sm.datasets.anes96.load_pandas().data  # 944 voters, US election survey 1996

3.5.1 More than two groups: analysis of variance

The alcohol content of 178 wines from three grape cultivars. Do the cultivars differ? We could run a t-test for every pair, but with many groups the pairs multiply, and so do the false alarms (Chapter 3.10). One test should look at all groups at once.

Analysis of variance (ANOVA) splits the variation into two parts. Each wine’s distance from the overall average $\bar x$ is the sum of two distances: from the overall average to its group’s average $\bar x_g$ (between groups), and from the group’s average to the wine (within the group). Squared and added up, the parts add up exactly:

$$\underbrace{\sum_{i} (x_i - \bar x)^2}_{\text{total}} = \underbrace{\sum_{g} n_g (\bar x_g - \bar x)^2}_{\text{between}} + \underbrace{\sum_{g} \sum_{i \in g} (x_i - \bar x_g)^2}_{\text{within}}$$

some = pd.concat(
    [wine[wine.target == g].sample(10, random_state=1) for g in range(3)]
).reset_index(drop=True)
grand = some.alcohol.mean()
group_means = some.groupby("target").alcohol.transform("mean")


def plot_part(ax: Axes, start: pd.Series | float, end: pd.Series, title: str) -> None:
    """Vertical lines from start to end for every wine, coloured by cultivar."""
    for i, row in some.iterrows():
        color = CULTIVAR_COLORS[int(row.target)]
        top = start if isinstance(start, float) else start[i]
        ax.plot([i, i], [top, end[i]], color=color, lw=1.4)
        ax.plot(i, row.alcohol, "o", ms=4, color=color)
    for g in range(3):
        ax.plot(
            [g * 10 - 0.4, g * 10 + 9.4],
            [group_means[g * 10]] * 2,
            color=CULTIVAR_COLORS[g],
            lw=2.5,
            alpha=0.7,
        )
    ax.axhline(grand, color=style.INK, ls="--", lw=1)
    ax.set(title=title, xticks=[], xlabel="wines, by cultivar")


def plot_sums(ax: Axes) -> None:
    """Between and within sums of squares stacked, next to the total, for all wines."""
    grand_all = wine.alcohol.mean()
    means = wine.groupby("target").alcohol.transform("mean")
    between = ((means - grand_all) ** 2).sum()
    within = ((wine.alcohol - means) ** 2).sum()
    ax.bar(0, between, color=style.GOLD, alpha=0.7, label="between")
    ax.bar(0, within, bottom=between, color=style.GREY, alpha=0.7, label="within")
    ax.bar(1, between + within, color=style.BLUE, alpha=0.5, label="total")
    ax.set_xticks([0, 1], ["between\n+ within", "total"])
    ax.legend(fontsize=7.5, loc="upper left")
    ax.set(title="All 178", ylabel="sum of squares")


fig = plt.figure(figsize=(11, 3.2))
grid = fig.add_gridspec(1, 4, width_ratios=[3, 3, 3, 1.6])
axes = [fig.add_subplot(grid[i]) for i in range(4)]
plot_part(axes[0], grand, some.alcohol, "Total")
plot_part(axes[1], grand, group_means, "Between")
plot_part(axes[2], group_means, some.alcohol, "Within")
plot_sums(axes[3])
axes[0].set_ylabel("alcohol (%)")
for ax in axes[1:3]:
    ax.sharey(axes[0])
fig.tight_layout()
plt.show()

Figure 3.5.1: The split behind ANOVA, for 30 of the wines (10 per cultivar). Left: each wine’s distance to the overall average (dashed): the total variation. Middle: the part between groups, from the overall average to each group’s average (thick lines). Right: the part within groups, from each group’s average to each wine. The bars on the far right show the three sums of squares for all 178 wines: between + within = total.

The F statistic compares the two parts, each divided by its degrees of freedom (k − 1 for k groups, N − k for N wines):

$$F = \frac{\text{between} / (k - 1)}{\text{within} / (N - k)}$$

If the cultivars did not differ, the group averages would differ only by chance, and F would be around 1. It follows the F distribution with k − 1 and N − k degrees of freedom.

def plot_by_group(ax: Axes, values: pd.Series, title: str) -> float:
    """Dots by cultivar with group averages; return the F statistic."""
    groups = [values[wine.target == k].to_numpy() for k in range(3)]
    for k, (group, color) in enumerate(zip(groups, CULTIVAR_COLORS)):
        ax.scatter(k + rng.uniform(-0.15, 0.15, len(group)), group, s=9, color=color, alpha=0.5)
        ax.plot([k - 0.3, k + 0.3], [group.mean()] * 2, color=color, lw=2.5)
    ax.axhline(values.mean(), color=style.INK, ls="--", lw=1)
    f = stats.f_oneway(*groups).statistic
    label = f"F = {f:.0f}" if f >= 10 else f"F = {f:.1f}"
    style.note(ax, label, x=0.97, y=0.1, ha="right", color=style.INK)
    ax.set_xticks(range(3), ["1", "2", "3"])
    ax.set(title=title, xlabel="cultivar", ylabel="alcohol (%)")
    return f


def closer(values: pd.Series, factor: float = 1 / 6) -> pd.Series:
    """The same spread within groups, with the group averages pulled towards the overall one."""
    means = values.groupby(wine.target).transform("mean")
    return values.mean() + factor * (means - values.mean()) + (values - means)


def plot_f(ax: Axes, observed: list[tuple[float, str]], dfs: tuple[int, int]) -> None:
    """The F distribution under no difference, with observed F values marked."""
    x = np.linspace(0.01, 8, 400)
    ax.plot(x, stats.f.pdf(x, *dfs), color=style.BLUE)
    for f, label in observed:
        shown = min(f, 7.8)
        ax.axvline(shown, color=style.CORAL, lw=1.4)
        ax.text(
            shown - 0.15, 0.6, label, rotation=90, ha="right", fontsize=9, color=style.CORAL
        )
    ax.set(title="If cultivars did not differ", xlabel="F", ylabel="how common", yticks=[])


fig, axes = style.subplots(3, height=3.1, width=10.5)
f_real = plot_by_group(axes[0], wine.alcohol, "The wines")
f_closer = plot_by_group(axes[1], closer(wine.alcohol), "Closer averages, same noise")
plot_f(
    axes[2], [(f_real, f"real: F = {f_real:.0f} →"), (f_closer, "closer")], (2, len(wine) - 3)
)
fig.tight_layout()
plt.show()

Figure 3.5.2: Left: alcohol content of the wines by cultivar, with group averages (lines) and the overall average (dashed). Middle: the same wines with the group averages moved closer together, made up for comparison. Right: the F distribution for no difference, with the F of both; the real one is far beyond the curve.

The noise inside the groups is the same in the left two panels; only the gaps between the averages differ. Far apart, they stand out from the noise (F = 135); close together, they blur into it. ANOVA judges the gaps between averages against the spread inside the groups, just like the t-test, which is ANOVA’s special case for two groups ($F = t^2$). The share of the total that lies between groups, between / total, is called $\eta^2$ (“eta squared”): the ANOVA’s R² (Chapter 4.1). Here it is 0.61.

ANOVA assumes roughly normal values within each group and similar spreads; Welch’s version relaxes the second, as for the t-test. The rank version, for skewed data, is the Kruskal-Wallis test, which runs ANOVA on the ranks. Its statistic H measures how far the groups’ average ranks lie from the overall average rank; when the groups do not differ, it follows a chi-square distribution with (groups − 1) degrees of freedom, and shuffling the labels shows the same:

alcohol_ranks = stats.rankdata(wine.alcohol)
kruskal = stats.kruskal(*[wine.alcohol[wine.target == c] for c in range(3)])


def shuffled_h() -> float:
    """Kruskal-Wallis H after shuffling the cultivar labels."""
    labels = rng.permutation(wine.target.to_numpy())
    return float(stats.kruskal(*[wine.alcohol[labels == c] for c in range(3)]).statistic)


null_h = np.array([shuffled_h() for _ in range(5000)])
fig, axes = style.subplots(2, height=3.2, width=10.5)
for c, color in enumerate(CULTIVAR_COLORS):
    ranks_c = alcohol_ranks[wine.target == c]
    axes[0].scatter(
        c + rng.uniform(-0.15, 0.15, len(ranks_c)), ranks_c, s=10, color=color, alpha=0.6
    )
    axes[0].plot([c - 0.3, c + 0.3], [ranks_c.mean()] * 2, color=color, lw=2.5)
axes[0].axhline(alcohol_ranks.mean(), color=style.INK, ls="--", lw=1)
axes[0].set_xticks(range(3), CULTIVARS)
axes[0].set(ylabel="rank of alcohol (1 = lowest)", title="Ranks")
style.null_panel(
    axes[1],
    kruskal.statistic,
    samples=null_h,
    dist=stats.chi2(2),
    tail="right",
    xlabel="Kruskal-Wallis H",
)
axes[1].set_title("H when the labels are shuffled")
fig.tight_layout()
plt.show()

Figure 3.5.3: The Kruskal-Wallis test for alcohol by cultivar. Left: the wines’ alcohol as ranks, 1 to 178, with each cultivar’s average rank (lines) and the overall average rank (dashed). Right: H for 5,000 shufflings of the cultivar labels (grey) and the chi-square curve with 2 degrees of freedom (blue). The observed H, far to the right, gives p < 10⁻²⁰.

3.5.2 Which groups differ?

A small p-value says that some cultivars differ, not which. A post-hoc comparison of all pairs, corrected for the number of pairs, answers that:

tukey = pairwise_tukeyhsd(wine.alcohol, wine.target.map(dict(enumerate(CULTIVARS))))
table = pd.DataFrame(tukey.summary().data[1:], columns=tukey.summary().data[0])


def plot_pairs(ax: Axes) -> None:
    """One interval per pair of cultivars, against zero."""
    for i, row in table.iterrows():
        ax.plot([row.lower, row.upper], [i, i], color=style.BLUE, lw=2.5)
        ax.plot(row.meandiff, i, "o", ms=7, color=style.BLUE)
    ax.axvline(0, color=style.INK, lw=1)
    labels = [f"{row.group2[-1]} − {row.group1[-1]}" for _, row in table.iterrows()]
    ax.set_yticks(range(len(table)), [f"cultivar {label}" for label in labels])
    ax.invert_yaxis()
    ax.set(xlabel="difference in average alcohol (percentage points)")


fig, ax = plt.subplots(figsize=(8, 2.4))
plot_pairs(ax)
plt.show()

Figure 3.5.4: Every pair of cultivars: the difference in average alcohol with a 95% interval, widened for having three comparisons (Tukey’s method). All three intervals leave out zero, so all three pairs differ.

Tukey’s method widens every interval just enough that the chance of any of them missing is 5%. Simpler corrections, such as Bonferroni (Chapter 3.10), work for any set of comparisons but are wider. If only a few comparisons were planned before seeing the data, test just those.

3.5.3 Two factors at once

Often the groups are defined by two things. In the RAND health study, people rated their health (excellent, good, fair, poor) and said whether a physical limitation held them back. Both probably affect doctor visits. Two-way ANOVA separates the effect of each factor, and asks a third question: does the effect of one depend on the other? That is an interaction (Chapter 4.2 has the regression view).

health["rating"] = pd.Categorical(
    np.select(
        [health.hlthp == 1, health.hlthf == 1, health.hlthg == 1],
        ["poor", "fair", "good"],
        "excellent",
    ),
    categories=["excellent", "good", "fair", "poor"],
    ordered=True,
)
health["limited"] = np.where(health.physlm > 0, "limited", "not limited")
health["log_visits"] = np.log1p(health.mdvis)
two_way = smf.ols("log_visits ~ C(rating) * C(limited)", health).fit()
anova_table = sm.stats.anova_lm(two_way, typ=2)


def plot_interaction(ax: Axes) -> None:
    """Group means with intervals, one line per limitation group."""
    for group, color, shift in (
        ("not limited", style.BLUE, -0.06),
        ("limited", style.CORAL, 0.06),
    ):
        part = health[health.limited == group].groupby("rating", observed=True).log_visits
        means, se = part.mean(), part.std() / np.sqrt(part.count())
        x = np.arange(4) + shift
        ax.errorbar(
            x, means, yerr=1.96 * se, fmt="o-", color=color, capsize=3, lw=1.6, label=group
        )
    ax.set_xticks(range(4), ["excellent", "good", "fair", "poor"])
    ax.legend(fontsize=8)
    ax.set(xlabel="self-rated health", ylabel="average log(1 + visits)")


fig, ax = plt.subplots(figsize=(8, 3.3))
plot_interaction(ax)
plt.show()

Figure 3.5.5: Average doctor visits (on a log scale: log of 1 + visits) by self-rated health, for people with (red) and without (blue) a physical limitation, with 95% intervals. Worse health and a limitation both go with more visits, and the lines are not parallel: the gap from a limitation grows as health gets worse. That is an interaction.

shown = anova_table.rename(
    index={
        "C(rating)": "health rating",
        "C(limited)": "limitation",
        "C(rating):C(limited)": "interaction",
        "Residual": "within groups",
    }
)
shown.columns = ["sum of squares", "degrees of freedom", "F", "p-value"]
shown["p-value"] = [
    "< 0.001" if p < 0.001 else f"{p:.3f}" if p == p else "" for p in shown["p-value"]
]
shown["F"] = [f"{f:.1f}" if f == f else "" for f in shown["F"]]
print(shown.round(1).to_markdown(colalign=["left"] + ["right"] * 4))

Table 3.5.1: Two-way ANOVA for log doctor visits: each row tests one source of variation.

	sum of squares	degrees of freedom	F	p-value
health rating	42.9	3	20.9	< 0.001
limitation	193.2	1	282.4	< 0.001
interaction	18.4	3	9.0	< 0.001
within groups	13808.3	20182

All three rows are clearly significant. With 20,190 people, even small effects are, so read the figure for size: the limitation adds about 0.2 on this log scale for people in excellent health, and twice that for people in poor health.

3.5.4 Tables of counts: the chi-square test

Some data are just counts in a table. In a 1974 survey, 6,366 married women said how religious they were and whether they had had an affair. Does the share with an affair depend on how religious a woman is?

If the two did not depend on each other, every row of the table would have the same share of affairs, the overall 32%. That gives the expected count in every cell: row total × column total / grand total. The chi-square statistic adds up how far each observed count O lies from its expected count E, squared and scaled:

$$\chi^2 = \sum_{\text{cells}} \frac{(O - E)^2}{E}$$

With no relation, it follows the chi-square distribution of Chapter 2.1, with (rows − 1) × (columns − 1) degrees of freedom.

had_affair = fair.affairs > 0
religious = fair.religious.astype(int)
observed = pd.crosstab(religious, had_affair)
expected = np.outer(observed.sum(axis=1), observed.sum(axis=0)) / observed.to_numpy().sum()
chi2 = stats.chi2_contingency(observed)
residuals = (observed.to_numpy() - expected) / np.sqrt(expected)


def plot_shares(ax: Axes) -> None:
    """Share with an affair per level, with intervals."""
    share = observed[True] / observed.sum(axis=1)
    half = 1.96 * np.sqrt(share * (1 - share) / observed.sum(axis=1))
    ax.bar(
        share.index, share, color=style.BLUE, alpha=0.6, yerr=half, capsize=4, ecolor=style.INK
    )
    ax.axhline(had_affair.mean(), color=style.INK, ls="--", lw=1)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set_xticks(share.index)
    ax.set(title="Share with an affair", xlabel="how religious", ylabel="share of women")


def plot_expected(ax: Axes) -> None:
    """Observed and expected numbers of women with an affair, per level."""
    levels = observed.index.to_numpy()
    ax.bar(
        levels - 0.2, observed[True], width=0.4, color=style.BLUE, alpha=0.6, label="observed"
    )
    ax.bar(
        levels + 0.2, expected[:, 1], width=0.4, color=style.GREY, alpha=0.6, label="expected"
    )
    ax.legend(fontsize=8)
    style.note(
        ax,
        f"χ² = {chi2.statistic:.0f}, df = 3\np = {chi2.pvalue:.0e}",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.set_ylim(0, observed[True].max() * 1.35)
    ax.set_xticks(levels)
    ax.set(title="Observed and expected", xlabel="how religious", ylabel="women with an affair")


def plot_residuals(ax: Axes) -> None:
    """Standardised residuals of every cell."""
    limit = abs(residuals).max()
    ax.imshow(residuals.T, cmap="RdBu_r", vmin=-limit, vmax=limit)
    for i in range(residuals.shape[0]):
        for j in range(2):
            ink = "white" if abs(residuals[i, j]) > 0.6 * limit else "black"
            ax.text(
                i,
                j,
                f"{residuals[i, j]:+.1f}",
                ha="center",
                va="center",
                fontsize=10,
                color=ink,
            )
    ax.set_xticks(range(4), ["1", "2", "3", "4"])
    ax.set_yticks([0, 1], ["no affair", "affair"])
    ax.grid(False)
    ax.set(title="Residuals", xlabel="how religious")


fig, axes = style.subplots(3, height=3.1, width=11)
plot_shares(axes[0])
plot_expected(axes[1])
plot_residuals(axes[2])
fig.tight_layout()
plt.show()

Figure 3.5.6: Left: the share who had an affair, by how religious the women said they were (1 = not, 4 = strongly), with 95% intervals; the dashed line is the overall share. Middle: observed counts of affairs against the counts expected if religion made no difference. Right: each cell’s standardised residual, (O − E)/√E: red cells have more than expected, blue fewer.

How far is “far beyond”? If religion and affairs were unrelated, the religion answers could be shuffled among the women without changing anything. Each shuffled table has its own chi-square statistic, and their distribution is the yardstick:

had = had_affair.to_numpy()
religion_values = religious.to_numpy()


def shuffled_chi2() -> float:
    """Chi-square statistic after shuffling the religion answers."""
    table = pd.crosstab(rng.permutation(religion_values), had)
    return float(stats.chi2_contingency(table, correction=False).statistic)


null_chi2 = np.array([shuffled_chi2() for _ in range(2000)])
fig, ax = plt.subplots(figsize=(8, 3.1))
style.null_panel(
    ax,
    chi2.statistic,
    samples=null_chi2,
    dist=stats.chi2(3),
    tail="right",
    xlabel="chi-square statistic",
)
plt.show()

Figure 3.5.7: The chi-square statistic for 2,000 tables in which the women’s religion answers were shuffled (grey), and the chi-square curve with (4 − 1) × (2 − 1) = 3 degrees of freedom (blue), which the shuffles follow. Shuffled tables rarely exceed 10. The observed statistic (dashed) is far to the right of all of them.

The statistic is far beyond anything chance gives: the less religious women reported affairs more often. The residuals show where the table departs from independence: most of all, too few affairs among the strongly religious. A common effect size is Cramér’s V, $\sqrt{\chi^2 / (N \cdot (\min(\text{rows}, \text{columns}) - 1))}$, from 0 (no relation) to 1; here it is 0.13, a weak relation that is very clear only because N is large.

Cautions: the test says the rows differ, not why; religion may stand for age, region or something else. And the chi-square approximation needs expected counts of about 5 or more in every cell. For small tables with few counts, use Fisher’s exact test (Chapter 3.4).

3.5.5 Does a distribution fit? The goodness-of-fit test

The same statistic compares observed counts with any set of expected counts, not just those of independence. Benford’s law says that in many collections of numbers that span several orders of magnitude, the first digit is 1 about 30% of the time and 9 only about 5%: $P(d) = \log_{10}(1 + 1/d)$. Do the populations of the voters’ home towns in the 1996 survey follow it?

populations = anes.popul[anes.popul > 0]
first_digit = (populations / 10 ** np.floor(np.log10(populations))).astype(
    int
)  # first significant digit
counts = first_digit.value_counts().sort_index().reindex(range(1, 10), fill_value=0)
benford = np.log10(1 + 1 / np.arange(1, 10)) * len(populations)
fit = stats.chisquare(counts.to_numpy(), benford)


def plot_benford(ax: Axes) -> None:
    """Observed first-digit counts against Benford's expected counts."""
    ax.bar(counts.index, counts, color=style.BLUE, alpha=0.6, label="observed")
    ax.plot(counts.index, benford, "o-", color=style.CORAL, ms=5, label="Benford's law")
    style.note(
        ax,
        f"χ² = {fit.statistic:.0f}, df = 8, p = {fit.pvalue:.0e}",
        x=0.45,
        y=0.8,
        color=style.INK,
    )
    ax.legend(fontsize=8)
    ax.set(xlabel="first digit of the town's population", ylabel="voters", xticks=range(1, 10))


def simulated_fit() -> float:
    """Goodness-of-fit statistic for 716 first digits drawn from Benford's law itself."""
    draws = rng.multinomial(len(populations), benford / benford.sum())
    return float(((draws - benford) ** 2 / benford).sum())


null_fit = np.array([simulated_fit() for _ in range(10_000)])
fig, axes = style.subplots(2, height=3.2, width=10.5)
plot_benford(axes[0])
style.null_panel(
    axes[1],
    fit.statistic,
    samples=null_fit,
    dist=stats.chi2(8),
    tail="right",
    xlabel="chi-square statistic",
)
axes[1].set_title("If the law held")
fig.tight_layout()
plt.show()

Figure 3.5.8: Left: first digits of the home-town populations of 716 voters (bars) against Benford’s law (red). The falling pattern is there, but there are too many 3s and 9s and too few 4s and 8s. Right: the goodness-of-fit statistic for 10,000 simulated samples of 716 digits drawn from Benford’s law itself (grey), with the chi-square curve with 8 degrees of freedom (blue). Real samples from the law rarely go beyond 20; the voters’ digits give 38.

With k categories and nothing estimated from the data, the statistic has k − 1 degrees of freedom: 8 here. If parameters were fitted first (say the mean of a Poisson), subtract one more degree of freedom per parameter. The survey records populations in thousands, rounded, which may explain the misfit; tests like this one are used to spot doctored accounts, whose first digits are often too evenly spread.

3.5.6 In one picture

Question	Test	Without the normal shape
do several group averages differ?	one-way ANOVA (F test), Welch’s ANOVA	Kruskal-Wallis
which pairs differ?	Tukey’s comparisons	pairwise rank tests, corrected
two factors and their interaction	two-way ANOVA
are two counted variables related?	chi-square test of independence	Fisher’s exact test (small tables)
do counts fit a given distribution?	chi-square goodness-of-fit test

• ANOVA splits the total variation into between and within; F compares them.
• An interaction means one factor’s effect depends on the other: lines that are not parallel.
• The chi-square test compares observed counts with expected ones: $\sum (O - E)^2/E$.
• Report effect sizes: $\eta^2$ for ANOVA, Cramér’s V for tables. With big samples, everything is significant.

Next: correlation, and testing it.