3.11 Combining studies: meta-analysis

A meta-analysis combines the estimates of several studies into one, giving more weight to the more precise ones. It must also ask whether the studies measure the same effect, and whether the studies we can find are a fair sample of the studies that were done.

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

from statly import style

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


def simulate_studies(
    k: int, effect: float, spread: float, seed: int
) -> tuple[np.ndarray, np.ndarray]:
    """Estimates and standard errors of k studies of varying size; true effects vary by `spread`."""
    local = np.random.default_rng(seed)
    sizes = local.integers(30, 600, k)
    se = 2 / np.sqrt(sizes)  # standard error of a standardised mean difference, roughly
    true = effect + local.normal(0, spread, k)
    return true + local.normal(0, se), se


def pool(estimates: np.ndarray, se: np.ndarray, random: bool) -> tuple[float, float, float]:
    """Inverse-variance pooled estimate and its SE; with `random`, DerSimonian-Laird τ² is added."""
    w = 1 / se**2
    fixed = np.sum(w * estimates) / w.sum()
    q = np.sum(w * (estimates - fixed) ** 2)
    tau2 = max(0.0, (q - (len(se) - 1)) / (w.sum() - np.sum(w**2) / w.sum())) if random else 0.0
    w_star = 1 / (se**2 + tau2)
    return (
        float(np.sum(w_star * estimates) / w_star.sum()),
        float(np.sqrt(1 / w_star.sum())),
        tau2,
    )


def heterogeneity(estimates: np.ndarray, se: np.ndarray) -> tuple[float, float]:
    """Cochran's Q and I² (the share of variation beyond chance)."""
    w = 1 / se**2
    fixed = np.sum(w * estimates) / w.sum()
    q = np.sum(w * (estimates - fixed) ** 2)
    return float(q), float(max(0.0, (q - (len(se) - 1)) / q))

3.11.1 One answer from many studies

Twelve trials of the same treatment report effects of different sizes, some significant and some not. Counting significant results (“vote counting”) is misleading: small studies miss real effects (Chapter 3.10). A meta-analysis instead averages the estimates, weighting each by its precision, one over its squared standard error:

$$\hat\theta = \frac{\sum_i w_i \hat\theta_i}{\sum_i w_i}, \qquad w_i = \frac{1}{\text{SE}_i^2}, \qquad \text{SE}(\hat\theta) = \frac{1}{\sqrt{\sum_i w_i}}$$

The standard picture is the forest plot: one row per study with its estimate and interval, the size of the square showing its weight, and a diamond for the pooled result.

estimates, se = simulate_studies(12, 0.2, 0.0, seed=1)
pooled, pooled_se, _ = pool(estimates, se, random=False)


def forest(
    ax: Axes,
    estimates: np.ndarray,
    se: np.ndarray,
    diamonds: list[tuple[str, float, float, str]],
) -> None:
    """Studies as squares with intervals, and pooled results as diamonds below them."""
    weights = 1 / se**2
    k = len(estimates)
    for i, (est, s, w) in enumerate(zip(estimates, se, weights)):
        ax.plot([est - 1.96 * s, est + 1.96 * s], [i, i], color=style.INK, lw=1)
        ax.scatter(est, i, marker="s", s=200 * w / weights.max(), color=style.BLUE, zorder=3)
    for j, (name, est, s, color) in enumerate(diamonds):
        y = k + 0.8 + j
        ax.fill(
            [est - 1.96 * s, est, est + 1.96 * s, est],
            [y, y - 0.35, y, y + 0.35],
            color=color,
            alpha=0.8,
        )
        ax.text(est + 1.96 * s + 0.03, y, name, va="center", fontsize=8.5)
    ax.axvline(0, color=style.INK, lw=0.8)
    ax.set_yticks(range(k), [f"study {i + 1}" for i in range(k)], fontsize=8)
    ax.invert_yaxis()
    ax.set(xlabel="effect (standardised mean difference)")


fig, ax = plt.subplots(figsize=(8, 4.4))
forest(
    ax,
    estimates,
    se,
    [
        (
            f"pooled {pooled:.2f} ({pooled - 1.96 * pooled_se:.2f} to {pooled + 1.96 * pooled_se:.2f})",
            pooled,
            pooled_se,
            style.CORAL,
        )
    ],
)
plt.show()

Figure 3.11.1: Forest plot of twelve simulated trials of a treatment with a true standardised effect of 0.2. Each square is a study’s estimate with its 95% interval; larger squares carry more weight. Only some of the trials are significant on their own (their intervals exclude 0); the pooled estimate (diamond) is precise and clearly above 0.

3.11.2 Do the studies agree? Heterogeneity and random effects

The calculation above, a fixed-effect meta-analysis, assumes every study estimates the same true effect and differs only by chance. Studies of different populations, doses or follow-up times often estimate genuinely different effects. Cochran’s Q measures how much the estimates scatter beyond what their standard errors predict, and I² expresses that excess as a share:

$$Q = \sum_i w_i (\hat\theta_i - \hat\theta)^2, \qquad I^2 = \frac{Q - (k - 1)}{Q}$$

A random-effects model treats the true effects as drawn from a distribution with spread τ, and adds τ² to each study’s variance before weighting. The pooled estimate is then the average effect, and its interval widens; a prediction interval shows where the effect of a new study is likely to lie.

estimates_h, se_h = simulate_studies(12, 0.2, 0.15, seed=2)
fixed_h, fixed_se_h, _ = pool(estimates_h, se_h, random=False)
random_h, random_se_h, tau2 = pool(estimates_h, se_h, random=True)
q, i2 = heterogeneity(estimates_h, se_h)
predict_half = stats.t.ppf(0.975, 10) * np.sqrt(tau2 + random_se_h**2)

fig, ax = plt.subplots(figsize=(8, 4.8))
forest(
    ax,
    estimates_h,
    se_h,
    [
        (f"fixed effect {fixed_h:.2f}", fixed_h, fixed_se_h, style.GREY),
        (f"random effects {random_h:.2f}", random_h, random_se_h, style.CORAL),
        ("prediction interval", random_h, predict_half / 1.96, style.GOLD),
    ],
)
style.note(
    ax,
    f"Q = {q:.1f} on 11 df, I² = {i2:.0%}, τ = {np.sqrt(tau2):.2f}",
    x=0.55,
    y=0.05,
    color=style.INK,
)
plt.show()

Figure 3.11.2: Twelve simulated studies whose true effects vary around 0.2 with a standard deviation of 0.15. The estimates scatter more than their intervals allow (I² well above 0). The fixed-effect diamond (grey) is too narrow; the random-effects diamond (red) is wider, and the prediction interval for a new study (gold) wider still, reaching below zero.

Is Q of 21 more than chance? If all twelve studies measured the same effect, their estimates would scatter only by their standard errors, and Q would follow a chi-square distribution with k − 1 = 11 degrees of freedom. Simulating studies with exactly these standard errors and no real differences shows it:

def simulated_q() -> float:
    """Q for twelve studies with the same standard errors and no heterogeneity."""
    estimates = 0.2 + rng.normal(0, se_h)
    return heterogeneity(estimates, se_h)[0]


null_q = np.array([simulated_q() for _ in range(5000)])
fig, ax = plt.subplots(figsize=(8, 3.1))
style.null_panel(
    ax, q, samples=null_q, dist=stats.chi2(len(se_h) - 1), tail="right", xlabel="Cochran's Q"
)
plt.show()

Figure 3.11.3: Cochran’s Q for 5,000 simulated sets of twelve studies with the same standard errors as above and one common true effect (grey), with the chi-square curve with 11 degrees of freedom (blue). The observed Q of the heterogeneous studies (dashed) is beyond most of them, but not far beyond: with only twelve studies, the test for heterogeneity has little power, which is why I² and τ are reported whatever its p-value.

When heterogeneity is large, the average effect may describe no real setting well, and the more interesting question is why studies differ: meta-regression relates the study estimates to study characteristics such as dose or population.

3.11.3 Missing studies: publication bias

A meta-analysis can only combine studies it finds, and studies with significant results are more likely to be published. Small studies are significant only if they happen to find a large effect, so the published small studies overstate the effect. The funnel plot shows each study’s estimate against its precision. Without bias, the points form a symmetric funnel, wide at the bottom (small, imprecise studies) and narrow at the top. Publication bias bites a corner out of the funnel.

all_est, all_se = simulate_studies(80, 0.1, 0.0, seed=5)
significant = all_est / all_se > 1.96
published = significant | (np.random.default_rng(6).uniform(size=80) < 0.25)


def funnel(ax: Axes, est: np.ndarray, se: np.ndarray, title: str) -> None:
    """Estimate against standard error (inverted, precise studies on top), with the pooled line."""
    pooled_est, _, _ = pool(est, se, random=False)
    ax.scatter(est, se, s=16, color=style.BLUE, alpha=0.7)
    s_grid = np.linspace(0.001, 0.4, 50)
    ax.plot(pooled_est - 1.96 * s_grid, s_grid, color=style.GREY, ls="--", lw=1)
    ax.plot(pooled_est + 1.96 * s_grid, s_grid, color=style.GREY, ls="--", lw=1)
    ax.axvline(pooled_est, color=style.CORAL, lw=1.5)
    ax.axvline(0.1, color=style.INK, lw=1, ls=":")
    style.note(
        ax,
        f"{len(est)} studies, pooled {pooled_est:.2f}\n(truth 0.10, dotted)",
        x=0.03,
        y=0.97,
        color=style.INK,
    )
    ax.set(
        title=title,
        xlabel="estimated effect",
        ylabel="standard error",
        ylim=(0.4, 0),
        xlim=(-0.8, 1.1),
    )


fig, axes = style.subplots(2, height=3.6, width=10)
funnel(axes[0], all_est, all_se, "All studies")
funnel(axes[1], all_est[published], all_se[published], "Published studies")
fig.tight_layout()
plt.show()

Figure 3.11.4: Funnel plots of 80 simulated studies of a treatment with a true effect of 0.1. Left: all studies. The funnel is symmetric around the truth, and the pooled estimate is close to it. Right: only the studies that were significant, plus a quarter of the others, get published. The small, imprecise studies near the bottom left are missing, the funnel is lopsided, and the pooled estimate of the published studies is biased upwards.

Formal checks (such as Egger’s test, a regression of the standardised estimate on precision) and corrections (such as trim and fill) exist, but are weak with few studies. The real remedy is registering trials before they start and publishing every result; failing that, a funnel plot is a warning, not a proof.

3.11.4 In one picture

• Pool estimates by inverse-variance weights $w_i = 1/\text{SE}_i^2$; show them in a forest plot.
• Q and I² measure disagreement beyond chance; random effects add τ² and widen the interval; report a prediction interval.
• Meta-regression asks why studies differ.
• Publication bias shows as an asymmetric funnel plot; prevention (registration) beats correction.

Next: Bayesian statistics, where the unknowns themselves get probability distributions.