5.2 Surveys and sampling

A survey learns about a whole population from a sample. Its precision comes from the sample size; its accuracy comes from how the sample was chosen. A small random sample beats a huge biased one.

Code: a simulated population and style used in this chapter

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes

from statly import style

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

# A simulated country of 200,000 adults in four regions, with weekly hours of exercise.
POP = 200_000
region = rng.choice(4, POP, p=[0.4, 0.3, 0.2, 0.1])
REGION_NAMES = ["north", "south", "east", "west"]
REGION_MEANS = np.array([2.0, 3.5, 5.0, 7.0])
hours = np.maximum(0, REGION_MEANS[region] + rng.normal(0, 1.5, POP))
TRUE_MEAN = hours.mean()


def plot_estimates(ax: Axes, estimates: dict[str, np.ndarray]) -> None:
    """Histograms of repeated survey estimates, with the true population mean."""
    all_values = np.concatenate(list(estimates.values()))
    bins = np.linspace(all_values.min(), all_values.max(), 50)
    for (name, values), color in zip(
        estimates.items(), (style.GREY, style.BLUE, style.GREEN, style.CORAL)
    ):
        ax.hist(
            values, bins=bins, color=color, alpha=0.5, label=f"{name}: sd {values.std():.3f}"
        )
    ax.axvline(TRUE_MEAN, color=style.INK, ls="--", lw=1.2)
    ax.legend(fontsize=8)
    ax.set(xlabel="estimated mean hours (dashed: truth)", ylabel="surveys")

5.2.1 Population, frame and sample

A survey has a target population (all adults in the country), a sampling frame from which the sample is actually drawn (a register, a phone list, the visitors of a website), and the sample of people who end up answering. Errors creep in at each step: people missing from the frame (coverage error), the luck of the draw (sampling error), people who do not answer (nonresponse), and questions that are misunderstood (measurement error). Only the sampling error shrinks as the sample grows.

Code

def target(ax: Axes, bias: float, spread: float, title: str) -> None:
    """Twenty shots around a centre offset by `bias`, with spread `spread`."""
    for radius in (3, 2, 1):
        ax.add_patch(plt.Circle((0, 0), radius, fill=False, color=style.GREY, lw=1))
    shots = rng.normal(0, spread, (20, 2)) + [bias, bias * 0.6]
    ax.scatter(shots[:, 0], shots[:, 1], s=14, color=style.CORAL, zorder=3)
    ax.set(xlim=(-3.3, 3.3), ylim=(-3.3, 3.3), aspect="equal", title=title)
    ax.set_axis_off()


fig, axes = style.subplots(4, height=2.6, width=10.5)
for ax, (bias, spread, title) in zip(
    axes,
    (
        (0, 0.3, "no bias, precise"),
        (0, 1.1, "no bias, imprecise"),
        (1.6, 0.3, "biased, precise"),
        (1.6, 1.1, "biased, imprecise"),
    ),
):
    target(ax, bias, spread, title)
fig.tight_layout()
plt.show()

Figure 5.2.1: Accuracy and precision of a survey, as shots at a target whose centre is the true value. Random sampling error scatters the shots (variance) and shrinks with sample size. Bias moves the whole cloud (a frame that misses part of the population, or a response rate that depends on the answer), and does not shrink at all.

5.2.2 Simple random samples

In a simple random sample, every set of n people is equally likely to be chosen. The sample mean is then unbiased, with standard error

$\text{SE}(\bar y) = \frac{s}{\sqrt{n}} \sqrt{1 - \frac{n}{N}}$

The last factor, the finite population correction, matters only when the sample is a sizeable share of the population N. Remarkably, the precision depends on n, not on N: a random sample of 1,000 people measures a country of 300 million as well as a town of 100,000.

5.2.3 Stratified samples

If the population falls into groups (strata) that differ in the outcome, drawing a separate random sample within each stratum, in proportion to its size, removes the luck of how many people from each stratum happen to be drawn. The estimate is a weighted average of the stratum means, $\bar y_{\text{st}} = \sum_h W_h \bar y_h$ , with $W_h$ the stratum’s share of the population.

5.2.4 Cluster samples

Visiting 1,000 randomly scattered people is expensive. Cluster sampling draws whole groups instead, such as villages, schools or households, and measures everyone (or a sample) in them. It is cheaper per person, but people in the same cluster tend to be alike, so each extra person tells less than an independent one. The loss is measured by the design effect,

$\text{deff} = 1 + (m - 1)\rho$

where m is the number of people per cluster and ρ the intraclass correlation (Chapter 4.7). With 20 people per village and ρ = 0.1, a cluster sample of 1,000 is worth a simple random sample of 1,000 / 2.9 ≈ 345.

Code

village = rng.integers(0, 2000, POP)  # 2,000 villages; add a shared village effect
village_effect = rng.normal(0, 1.2, 2000)
hours_v = np.maximum(0, hours + village_effect[village])
TRUE_MEAN = hours_v.mean()
strata_share = np.bincount(region, minlength=4) / POP
members = [np.flatnonzero(region == h) for h in range(4)]
village_members = [np.flatnonzero(village == v) for v in range(2000)]


def simple_random() -> float:
    """Mean of 400 people drawn at random."""
    return hours_v[rng.choice(POP, 400, replace=False)].mean()


def stratified() -> float:
    """Weighted mean of proportional samples from each region."""
    sizes = np.round(400 * strata_share).astype(int)
    return float(
        sum(
            strata_share[h] * hours_v[rng.choice(members[h], sizes[h], replace=False)].mean()
            for h in range(4)
        )
    )


def clustered() -> float:
    """Mean of 20 people in each of 20 random villages."""
    chosen = rng.choice(2000, 20, replace=False)
    return float(
        np.mean(
            [hours_v[rng.choice(village_members[v], 20, replace=False)].mean() for v in chosen]
        )
    )


fig, ax = plt.subplots(figsize=(8, 3.3))
plot_estimates(
    ax,
    {
        "simple random": np.array([simple_random() for _ in range(2000)]),
        "stratified by region": np.array([stratified() for _ in range(2000)]),
        "20 villages × 20": np.array([clustered() for _ in range(2000)]),
    },
)
plt.show()

Figure 5.2.2: Estimating average weekly exercise in the simulated population (true mean: dashed line), with 2,000 repeated surveys of 400 people under three designs. Stratifying by region (blue) is more precise than simple random sampling (grey), because regions differ a lot. Sampling 20 villages of 20 people (green), where people in a village resemble each other, is much less precise.

5.2.5 Weights

When some groups are deliberately sampled more heavily, for instance to have enough people from a small region, each person stands for a different number of people in the population. Every respondent gets a weight, one over their probability of selection, and estimates use weighted averages:

$\bar y_w = \frac{\sum_i w_i y_i}{\sum_i w_i}, \qquad w_i = \frac{1}{\pi_i}$

Weights are also used to repair nonresponse: if young people answered less often, they are weighted up so that the sample matches the population’s known age distribution (post-stratification, or raking when several variables are matched at once).

Code

equal = np.array([100, 100, 100, 100])


def equal_allocation(weighted: bool) -> float:
    """Mean from 100 people per region, with or without design weights."""
    values = [hours_v[rng.choice(members[h], equal[h], replace=False)] for h in range(4)]
    if not weighted:
        return float(np.mean(np.concatenate(values)))
    return float(sum(strata_share[h] * values[h].mean() for h in range(4)))


fig, axes = style.subplots(2, height=3.2, width=10.5)
x = np.arange(4)
axes[0].bar(x - 0.2, strata_share, width=0.4, color=style.BLUE, alpha=0.7, label="population")
axes[0].bar(x + 0.2, equal / 400, width=0.4, color=style.GREY, alpha=0.7, label="sample")
axes[0].set_xticks(x, REGION_NAMES)
axes[0].yaxis.set_major_formatter(style.PERCENT)
axes[0].legend(fontsize=8)
axes[0].set(ylabel="share")
plot_estimates(
    axes[1],
    {
        "unweighted": np.array([equal_allocation(False) for _ in range(2000)]),
        "weighted": np.array([equal_allocation(True) for _ in range(2000)]),
    },
)
fig.tight_layout()
plt.show()

Figure 5.2.3: A survey of 400 that takes 100 people from each region, although the west holds only 10% of the population and exercises most. Left: share of the population and of the sample by region. Right: 2,000 repeated surveys. The unweighted mean (grey) is biased upwards, because the west is overrepresented; weighting each person by their region’s share over its share of the sample (blue) removes the bias.

5.2.6 Big biased samples

Weights can only fix what they can see. If the chance of answering depends on the answer itself, people who exercise more being more willing to answer a survey about exercise, no weighting by region or age removes the bias, and a larger sample only makes the wrong answer more precise. Online opinion polls, app data and convenience samples of any size suffer from this.

Code

willingness = 1 / (1 + np.exp(-(hours_v - 4) / 1.5))
SIZES = np.array([100, 300, 1000, 3000, 10_000, 30_000, 100_000])


def typical_error(n: int, self_selected: bool, repeats: int = 30) -> float:
    """Median absolute error of the sample mean."""
    p = willingness / willingness.sum() if self_selected else None
    errors = [
        abs(hours_v[rng.choice(POP, n, replace=False, p=p)].mean() - TRUE_MEAN)
        for _ in range(repeats)
    ]
    return float(np.median(errors))


fig, ax = plt.subplots(figsize=(8, 3.3))
ax.plot(
    SIZES,
    [typical_error(n, False) for n in SIZES],
    "o-",
    color=style.GREY,
    lw=2,
    label="random sample",
)
ax.plot(
    SIZES,
    [typical_error(n, True) for n in SIZES],
    "o-",
    color=style.CORAL,
    lw=2,
    label="self-selected sample",
)
ax.set(
    xscale="log",
    yscale="log",
    xlabel="respondents (log scale)",
    ylabel="typical error (hours, log scale)",
)
ax.legend(fontsize=8)
plt.show()

Figure 5.2.4: Error of the estimated mean exercise, against the number of respondents (log scale). Grey: simple random samples; the error falls steadily with size. Red: a self-selected sample in which the chance of responding grows with the hours exercised. Its error stops falling at the size of the bias: 100,000 self-selected respondents are worse than 100 random ones.

5.2.7 In one picture

Errors come from the frame, the draw, nonresponse and measurement; only the draw shrinks with n.
A simple random sample gives $\text{SE} = s/\sqrt{n}$ ; precision depends on n, not on the population size.
Stratifying on what predicts the outcome gains precision; clustering saves money and loses precision by the design effect $1 + (m-1)\rho$ .
Unequal selection needs weights $1/\pi_i$ ; weights can also correct known nonresponse.
If response depends on the answer itself, a bigger sample does not help.

Next: what to do when some of the values are missing.