5.4 Cause and effect

A causal effect compares what happened with what would have happened otherwise, and only one of the two is ever seen. Randomisation makes the groups comparable by design; without it, every causal estimate rests on assumptions about how the data came about, which a causal diagram makes explicit.

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes

from statly import style

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


def ols_coef(y: np.ndarray, *columns: np.ndarray) -> np.ndarray:
    """Least squares coefficients (after the intercept) of y on the given columns."""
    return sm.OLS(y, sm.add_constant(np.column_stack(columns))).fit().params[1:]


def plot_estimates(ax: Axes, estimates: dict[str, float], truth: float) -> None:
    """Horizontal bars for several estimates of one effect, with the true value."""
    names = list(estimates)
    values = list(estimates.values())
    colors = [
        style.BLUE if abs(v - truth) < 0.15 * max(abs(truth), 1) else style.GREY for v in values
    ]
    ax.barh(names, values, color=colors, alpha=0.75)
    for i, v in enumerate(values):
        ax.text(max(v, 0), i, f" {v:.2f}", va="center", fontsize=8.5)
    ax.axvline(truth, color=style.CORAL, ls="--", lw=1.5)
    ax.axvline(0, color=style.INK, lw=0.8)
    ax.invert_yaxis()
    ax.set(xlabel="estimated effect (dashed: truth)")


def draw_dag(
    ax: Axes,
    nodes: dict[str, tuple[float, float]],
    edges: list[tuple[str, str]],
    title: str,
    dashed: tuple[str, ...] = (),
) -> None:
    """A small causal diagram: named circles and arrows between them."""
    for name, (x, y) in nodes.items():
        ax.add_patch(
            plt.Circle(
                (x, y),
                0.17,
                color=style.GREY if name in dashed else style.BLUE,
                alpha=0.3 if name in dashed else 0.35,
            )
        )
        ax.text(x, y, name, ha="center", va="center", fontsize=9)
    for a, b in edges:
        (x0, y0), (x1, y1) = nodes[a], nodes[b]
        d = np.hypot(x1 - x0, y1 - y0)
        ux, uy = (x1 - x0) / d, (y1 - y0) / d
        ax.annotate(
            "",
            xy=(x1 - 0.19 * ux, y1 - 0.19 * uy),
            xytext=(x0 + 0.19 * ux, y0 + 0.19 * uy),
            arrowprops=dict(arrowstyle="->", color=style.INK, lw=1.3),
        )
    ax.set(xlim=(-0.3, 2.3), ylim=(-0.4, 1.3), aspect="equal", title=title)
    ax.set_axis_off()

5.4.1 What a causal effect is

Did the treatment lower this patient’s blood pressure? The question compares two worlds: the one where the patient was treated, and the one where, other things equal, they were not. Each person has two potential outcomes, $Y(1)$ if treated and $Y(0)$ if not, and the causal effect for them is $Y(1) - Y(0)$. The trouble is that each person lives in only one world: we see one potential outcome, never both. This is the fundamental problem of causal inference.

What can be estimated is an average, the average treatment effect:

$$\text{ATE} = E[Y(1) - Y(0)] = E[Y(1)] - E[Y(0)]$$

local = np.random.default_rng(5)
untreated = local.normal(150, 12, 8)
treated = untreated - 8 + local.normal(0, 4, 8)
got_treatment = np.array([1, 0, 1, 1, 0, 0, 1, 0], bool)

fig, ax = plt.subplots(figsize=(8, 3.2))
for i in range(8):
    ax.plot([i, i], [untreated[i], treated[i]], color=style.GREY, lw=1)
    ax.scatter(
        i,
        treated[i],
        s=60,
        color=style.CORAL if got_treatment[i] else "none",
        edgecolors=style.CORAL,
        zorder=3,
    )
    ax.scatter(
        i,
        untreated[i],
        s=60,
        color=style.BLUE if not got_treatment[i] else "none",
        edgecolors=style.BLUE,
        zorder=3,
    )
ax.scatter([], [], color=style.CORAL, label="Y(1), if treated")
ax.scatter([], [], color=style.BLUE, label="Y(0), if not")
ax.legend(fontsize=8, loc="upper right")
ax.set_xticks(range(8), [f"person {i + 1}" for i in range(8)], fontsize=8)
style.note(
    ax,
    f"true average effect {np.mean(treated - untreated):.1f}",
    x=0.03,
    y=0.95,
    color=style.INK,
)
ax.set(ylabel="blood pressure")
plt.show()

Figure 5.4.1: Eight people, each with two potential outcomes: blood pressure if treated (red) and if not (blue). The true effect for each is the gap between them, about −8 on average. Each person is observed in one world only (filled), the other stays hidden (open). Comparing the observed treated with the observed untreated recovers the average effect only if the two groups are comparable.

5.4.2 Why randomisation works

If doctors give the treatment mostly to the sickest patients, the treated group starts out with higher blood pressure, and a simple comparison mixes the treatment’s effect with that difference. A confounder is a variable that affects both who gets the treatment and the outcome. Random assignment breaks every such link at once, measured or not: a coin, not the patient’s condition, decides the treatment, so the groups differ only by chance and by the treatment.

n = 2000
severity = rng.normal(0, 1, n)
chosen_by_doctor = rng.uniform(size=n) < 1 / (1 + np.exp(-2 * severity))
chosen_by_coin = rng.uniform(size=n) < 0.5


def blood_pressure(treated: np.ndarray) -> np.ndarray:
    """Outcome: depends strongly on severity, and falls by 8 with treatment."""
    return 150 + 15 * severity - 8 * treated + rng.normal(0, 8, n)


fig, axes = style.subplots(2, height=3.2, width=10)
for ax, assigned, title in (
    (axes[0], chosen_by_doctor, "Doctors choose"),
    (axes[1], chosen_by_coin, "A coin chooses"),
):
    outcome = blood_pressure(assigned)
    for group, color, name in (
        (False, style.BLUE, "untreated"),
        (True, style.CORAL, "treated"),
    ):
        ax.scatter(
            severity[assigned == group],
            outcome[assigned == group],
            s=3,
            alpha=0.3,
            color=color,
            label=name,
        )
    difference = outcome[assigned].mean() - outcome[~assigned].mean()
    style.note(
        ax,
        f"treated − untreated: {difference:+.1f}\n(true effect −8)",
        x=0.03,
        y=0.95,
        color=style.INK,
    )
    ax.legend(fontsize=8, loc="lower right", markerscale=4)
    ax.set(title=title, xlabel="severity", ylabel="blood pressure")
fig.tight_layout()
plt.show()

Figure 5.4.2: Simulated patients whose blood pressure falls by 8 with treatment. Left: in observational data, sicker patients (higher severity) are more likely to be treated, so the treated group has higher blood pressure and the treatment looks harmful. Right: with treatment assigned at random, the groups have the same severity on average and the comparison recovers the true effect.

That is why randomised experiments (Chapter 5.1) are the gold standard. When they are impossible, unethical or too slow, causal claims from observational data need the tools below, and each comes with assumptions that cannot be fully checked from the data.

5.4.3 Causal diagrams: what to adjust for

A causal diagram (a directed acyclic graph, or DAG) draws each variable as a node and each direct cause as an arrow. It answers the question that regression alone cannot: which variables should be adjusted for? Three small patterns cover most cases:

• Fork (confounder): Z → X and Z → Y. Z creates a link between X and Y that is not causal. Adjust for Z.
• Chain (mediator): X → M → Y. M is part of how X works. Adjusting for M removes the effect you want to measure. Do not adjust (for the total effect).
• Collider: X → C ← Y. C is caused by both. X and Y are unrelated, but adjusting for C (or selecting cases by it) creates a link. Do not adjust.

n = 5000
z = rng.normal(size=n)
x_fork = z + rng.normal(size=n)
y_fork = 1.0 * x_fork + 2 * z + rng.normal(size=n)
x_chain = rng.normal(size=n)
m_chain = x_chain + rng.normal(size=n)
y_chain = 1.0 * m_chain + rng.normal(size=n)  # all of X's effect runs through M
x_coll = rng.normal(size=n)
y_coll = rng.normal(size=n)  # no effect of X on Y at all
c_coll = x_coll + y_coll + rng.normal(size=n)
CASES = [
    (
        "Fork: adjust",
        {"Z": (1, 1), "X": (0, 0), "Y": (2, 0)},
        [("Z", "X"), ("Z", "Y"), ("X", "Y")],
        x_fork,
        y_fork,
        z,
        "Z",
        1,
    ),
    (
        "Chain: do not adjust",
        {"X": (0, 0), "M": (1, 1), "Y": (2, 0)},
        [("X", "M"), ("M", "Y")],
        x_chain,
        y_chain,
        m_chain,
        "M",
        1,
    ),
    (
        "Collider: do not adjust",
        {"X": (0, 1), "Y": (2, 1), "C": (1, 0)},
        [("X", "C"), ("Y", "C")],
        x_coll,
        y_coll,
        c_coll,
        "C",
        0,
    ),
]

fig, axes = plt.subplots(2, 3, figsize=(10.5, 5), gridspec_kw={"height_ratios": [1, 1]})
for col, (title, nodes, edges, x, y, third, name, truth) in enumerate(CASES):
    draw_dag(axes[0, col], nodes, edges, title)
    estimates = {
        "without " + name: ols_coef(y, x)[0],
        "adjusting for " + name: ols_coef(y, x, third)[0],
    }
    plot_estimates(axes[1, col], estimates, truth)
    axes[1, col].set(xlim=(-0.8, 2.3))
fig.tight_layout()
plt.show()

Figure 5.4.3: Top: the three patterns. Bottom: the estimated effect of X on Y from regressions with and without the third variable, on simulated data (true effect of X on Y: 1 in the fork and the chain, 0 in the collider). Adjusting is right for the fork, and wrong for the chain and the collider.

The collider is the least intuitive and among the most common traps. If a hospital only admits patients who have disease A or disease B (or both), A and B look negatively related among its patients even if they are unrelated in the population. Studying only successful companies, published studies, or people who answered a survey can create links the same way (selection bias).

The backdoor rule generalises the three patterns: to estimate the effect of X on Y, adjust for a set of variables that blocks every path from X to Y that starts with an arrow into X, without adjusting for anything caused by X.

5.4.4 Adjusting for measured confounders: propensity scores

If all confounders are measured, there are several ways to compare like with like. Regression adds them as predictors. The propensity score, the probability of treatment given the confounders, $e(x) = P(\text{treated} \mid x)$, offers two more:

• Matching: pair each treated case with an untreated case of similar propensity, and compare within pairs.
• Inverse probability weighting: weight each treated case by $1/e(x)$ and each untreated case by $1/(1 - e(x))$. Rarely treated cases that were treated count more, so the weighted groups resemble the whole population. The estimate compares the two weighted averages:

$$\widehat{\text{ATE}} = \frac{\sum_i T_i Y_i / e(x_i)}{\sum_i T_i / e(x_i)} - \frac{\sum_i (1 - T_i) Y_i / (1 - e(x_i))}{\sum_i (1 - T_i) / (1 - e(x_i))}$$

n = 3000
age = rng.normal(60, 10, n)
severity = rng.normal(0, 1, n)
logit = -0.5 + 0.06 * (age - 60) + 1.2 * severity
treated = rng.uniform(size=n) < 1 / (1 + np.exp(-logit))
outcome = 140 + 0.5 * (age - 60) + 12 * severity - 8 * treated + rng.normal(0, 8, n)
confounders = np.column_stack([age, severity])
propensity = sm.Logit(treated.astype(float), sm.add_constant(confounders)).fit(disp=0).predict()


def matched_difference() -> float:
    """Mean difference between each treated case and the untreated case nearest in propensity."""
    controls = np.flatnonzero(~treated)
    order = np.argsort(propensity[controls])
    sorted_scores = propensity[controls][order]
    gaps = []
    for i in np.flatnonzero(treated):
        j = np.clip(np.searchsorted(sorted_scores, propensity[i]), 1, len(order) - 1)
        nearest = (
            order[j - 1]
            if abs(sorted_scores[j - 1] - propensity[i]) < abs(sorted_scores[j] - propensity[i])
            else order[j]
        )
        gaps.append(outcome[i] - outcome[controls[nearest]])
    return float(np.mean(gaps))


treated_weight, control_weight = treated / propensity, ~treated / (1 - propensity)
weighted = np.average(outcome, weights=treated_weight) - np.average(
    outcome, weights=control_weight
)
estimates = {
    "naive difference": outcome[treated].mean() - outcome[~treated].mean(),
    "regression": ols_coef(outcome, treated.astype(float), age, severity)[0],
    "matching": matched_difference(),
    "weighting": weighted,
}

fig, axes = style.subplots(2, height=3.2, width=10)
bins = np.linspace(0, 1, 30)
axes[0].hist(propensity[~treated], bins=bins, color=style.BLUE, alpha=0.5, label="untreated")
axes[0].hist(propensity[treated], bins=bins, color=style.CORAL, alpha=0.5, label="treated")
axes[0].legend(fontsize=8)
axes[0].set(xlabel="propensity score", ylabel="patients")
plot_estimates(axes[1], estimates, -8)
fig.tight_layout()
plt.show()

Figure 5.4.4: Simulated observational data with two measured confounders (age and severity) that raise both the chance of treatment and the outcome; the true effect is −8. Left: the estimated propensity scores of treated and untreated patients. They overlap, which is what makes a comparison possible; where only one group exists, no method can compare. Right: the naive difference is badly biased; regression, nearest-neighbour matching on the propensity score, and inverse probability weighting all come close to the truth.

All of these rest on the same untestable assumption: no unmeasured confounding. If an important confounder was not measured, none of them can fix it. A sensitivity analysis asks how strong an unmeasured confounder would have to be to explain the effect away.

5.4.5 Before and after, with a comparison: difference in differences

When a change hits one group and not another, comparing the two groups’ changes removes everything that differs between the groups but stays constant over time. The difference in differences estimate is

$$(\bar Y_{\text{treated, after}} - \bar Y_{\text{treated, before}}) - (\bar Y_{\text{control, after}} - \bar Y_{\text{control, before}})$$

The key assumption is parallel trends: without the change, both groups would have moved in parallel. Several periods before the change let you check that the trends were parallel then.

months = np.arange(1, 25)
after = months > 12
trend = 0.6 * months
region_b = 40 + trend + rng.normal(0, 1, 24)
region_a = 50 + trend + 5 * after + rng.normal(0, 1, 24)
did = (region_a[after].mean() - region_a[~after].mean()) - (
    region_b[after].mean() - region_b[~after].mean()
)

fig, ax = plt.subplots(figsize=(8, 3.4))
ax.plot(months, region_a, "o-", color=style.CORAL, ms=4, label="region A (policy)")
ax.plot(months, region_b, "o-", color=style.BLUE, ms=4, label="region B (no policy)")
counterfactual = region_b + (region_a[~after].mean() - region_b[~after].mean())
ax.plot(
    months[after],
    counterfactual[after],
    "--",
    color=style.CORAL,
    lw=1.5,
    label="A without the policy",
)
ax.axvline(12.5, color=style.INK, lw=1)
style.note(
    ax, f"difference in differences: {did:.1f} (true 5)", x=0.03, y=0.95, color=style.INK
)
ax.legend(fontsize=8, loc="lower right")
ax.set(xlabel="month", ylabel="sales")
plt.show()

Figure 5.4.5: Monthly sales in two regions (simulated); a policy starts in region A after month 12 and raises sales by 5. The regions differ in level and both drift upwards, so neither ‘A after minus A before’ nor ‘A minus B after’ gives the effect. The dashed line is A’s path without the policy, projected from B’s change: the difference in differences recovers the effect. The parallel lines before month 12 support the assumption.

5.4.6 Natural experiments: instrumental variables

Sometimes a variable shifts the treatment but has no other route to the outcome: an instrument. Distance to the nearest specialist hospital affects whether a patient gets a specialist procedure, but, one hopes, not their health otherwise. Since the instrument is as good as random, the part of the treatment it explains is free of confounding. Two-stage least squares uses exactly that part: first predict the treatment from the instrument, then regress the outcome on the prediction. With one instrument it amounts to a ratio of two slopes:

$$\hat\beta_{\text{IV}} = \frac{\text{effect of } Z \text{ on } Y}{\text{effect of } Z \text{ on } X}$$

n = 5000
u = rng.normal(size=n)
instrument = rng.normal(size=n)
x_iv = 0.8 * instrument + u + rng.normal(size=n)
y_iv = 2 * x_iv - 3 * u + rng.normal(size=n)
first_stage = sm.OLS(x_iv, sm.add_constant(instrument)).fit().fittedvalues
iv_estimate = ols_coef(y_iv, first_stage)[0]

fig, axes = style.subplots(2, height=3, width=10)
draw_dag(
    axes[0],
    {"Z": (0, 0), "X": (1, 0), "Y": (2, 0), "U": (1.5, 1)},
    [("Z", "X"), ("X", "Y"), ("U", "X"), ("U", "Y")],
    "Instrument Z",
    dashed=("U",),
)
plot_estimates(
    axes[1],
    {"regression of Y on X": ols_coef(y_iv, x_iv)[0], "two-stage least squares": iv_estimate},
    2,
)
fig.tight_layout()
plt.show()

Figure 5.4.6: Left: the diagram. An unmeasured confounder U (grey) drives both treatment X and outcome Y, so regressing Y on X is biased. The instrument Z moves X and has no other path to Y. Right: on simulated data with a true effect of 2, ordinary regression is far off; two-stage least squares, using Z, recovers the effect, with a wider uncertainty.

Instruments must satisfy two conditions: they must really move the treatment (a weak instrument gives wild estimates), and they must affect the outcome only through the treatment, the exclusion restriction, which can be argued but never tested. Randomised experiments with imperfect compliance are the cleanest case: the random assignment is the instrument.

5.4.7 Thresholds: regression discontinuity

When treatment is given to everyone above a cut-off on some score (a scholarship above a test mark, a drug above a blood value), people just above and just below the cut-off are nearly identical except for the treatment. Fitting the outcome on each side of the threshold and measuring the jump at the cut-off estimates the effect there: regression discontinuity.

score = rng.uniform(20, 100, 1500)
tutored = score >= 60
grade = 30 + 0.4 * score + 4 * tutored + rng.normal(0, 4, 1500)
WINDOW = 15


def side_line(side: np.ndarray) -> np.ndarray:
    """Intercept and slope of grade on (score − 60) within the window, on one side."""
    near = side & (np.abs(score - 60) < WINDOW)
    return sm.OLS(grade[near], sm.add_constant(score[near] - 60)).fit().params


below, above = side_line(~tutored), side_line(tutored)
fig, ax = plt.subplots(figsize=(8, 3.4))
ax.scatter(score, grade, s=4, c=np.where(tutored, style.CORAL, style.BLUE), alpha=0.35)
for params, lo, hi in ((below, 60 - WINDOW, 60), (above, 60, 60 + WINDOW)):
    xs = np.array([lo, hi])
    ax.plot(xs, params[0] + params[1] * (xs - 60), color=style.INK, lw=2.5)
ax.axvline(60, color=style.INK, ls="--", lw=1)
style.note(
    ax, f"jump at 60: {above[0] - below[0]:.1f} (true 4)", x=0.03, y=0.95, color=style.INK
)
ax.set(xlabel="entrance test score", ylabel="later grade")
plt.show()

Figure 5.4.7: Simulated students: those scoring 60 or more in an entrance test get extra tuition, which raises a later grade by 4 points. Grades rise with the test score anyway. Lines fitted separately on each side, within 15 points of the cut-off, show the jump at 60, the estimated effect for students near the threshold.

The estimate is local: it says what tuition does for students near 60, not for those at 30 or 95. The design fails if people can push themselves across the threshold, so check that the number of cases just below and just above the cut-off is similar.

5.4.8 In one picture

• A causal effect compares potential outcomes $Y(1) - Y(0)$; only one is ever observed.
• Randomisation makes groups comparable in everything, measured or not.
• Causal diagrams decide what to adjust for: adjust for confounders (forks), not for mediators (chains) or colliders.
• With all confounders measured: regression, matching or weighting on the propensity score; check overlap; no method fixes unmeasured confounding.
• Difference in differences needs parallel trends; instrumental variables need an instrument that acts only through the treatment; regression discontinuity needs a sharp, unmanipulated threshold.

Next: Part 6, supervised learning, where the goal shifts from understanding an outcome to predicting it.