Causality: selection bias, randomisation and omitted variables

Prediction asks “what will happen?” Econometrics asks “what will happen if I change something?” Does health insurance make people healthier? Does this execution algorithm actually reduce slippage, or do we just use it on easier orders? Answering needs more than a good fit.

The selection-bias problem

Let $D_i = 1$ if person $i$ has insurance and $0$ otherwise. Each person has two potential outcomes: $Y_{1i}$ with insurance and $Y_{0i}$ without. We only ever see one. Comparing group averages gives

The selection-bias term asks: would the insured have been healthier even without insurance? Usually yes: richer, more educated people are both more likely to buy insurance and healthier anyway. The naive comparison mixes the two.

Randomise

If a coin decides who is treated, treatment can’t be related to anything else about people. By the law of large numbers both groups have the same average $Y_0$, so the selection-bias term goes to zero and the simple difference in means is the causal effect. After randomising, “check for balance”: compare the groups’ averages on age, income and so on; they should match.

Statistics for comparing two groups

• Sample mean $\bar Y$ is unbiased for $E[Y]$, with sampling variance $\sigma_Y^2/n$ and standard error $\widehat{\operatorname{SE}}(\bar Y) = S(Y)/\sqrt n$.
• $t$-statistic for “the mean is $\mu$”: $t = (\bar Y - \mu)/\widehat{\operatorname{SE}}(\bar Y)$; $|t| > 2$ is the usual bar, and $\bar Y \pm 2\,\widehat{\operatorname{SE}}$ is an approximate 95% interval.
• Two groups: $\widehat{\operatorname{SE}}(\bar Y^1 - \bar Y^0) = S(Y)\sqrt{\tfrac{1}{n_1} + \tfrac{1}{n_0}}$ (pooled), and $t = (\bar Y^1 - \bar Y^0)/\widehat{\operatorname{SE}}$.

A large $t$ needs both a real effect and a small standard error. A small $t$ may mean there is no effect, or just not enough data to see it.

Regression as a substitute for randomisation

Often you can’t randomise. Regression controls for observed differences:

where $A_i$ holds control variables (age, income, …). The key assumption is the conditional independence assumption: once you compare people with the same controls, treatment is as good as random. Continuous variables enter as they are; categories enter as one-hot (dummy) variables. Regressing on $\log Y$ makes $\beta$ read as an approximate percentage effect.

Omitted-variable bias

What if you left out an important control? Compare the short regression (without the control $A$) and the long one (with it):

where $\gamma$ is $A$’s effect on $Y$ in the long regression and $\pi_1$ is the slope from regressing $A$ on the treatment. Bias needs both: the omitted variable must move with the treatment ($\pi_1 \ne 0$) and matter for the outcome ($\gamma \ne 0$). Even without data on $A$, you can reason about the sign of the bias.

Regression facts worth knowing

• Regression gives the best linear approximation to the conditional expectation function $E[Y\mid X]$, which is itself the minimum-mean-squared-error predictor: $E[Y\mid X] = \arg\min_{m} E\big[(Y - m(X))^2\big]$.
• Residuals have mean zero and are uncorrelated with every input and with the fitted values.
• ANOVA identity: $\operatorname{Var}(Y) = \operatorname{Var}\big(E[Y\mid X]\big) + E\big[\operatorname{Var}(Y\mid X)\big]$.
• If the error variance changes with $X$ (heteroskedasticity, which any nonlinear CEF or a 0/1 outcome produces), use robust standard errors. They are usually close to the conventional ones.
• Sampling weights (inverse sampling probabilities) fix non-random samples; heteroskedasticity alone is not a reason to use weighted least squares.
• A saturated model (a dummy for every combination of discrete inputs) fits the CEF exactly, at the cost of noisy coefficients.

In code

import numpy as np
rng = np.random.default_rng(11)

n = 100_000
health = rng.standard_normal(n)                        # unobserved
chooses = rng.random(n) < 1 / (1 + np.exp(-2.2 * health))
y0 = 70 + 7 * health + 3 * rng.standard_normal(n)
effect = 2.0

naive = (y0[chooses] + effect).mean() - y0[~chooses].mean()
print("people choose:", round(naive, 2))              # far above 2: selection bias

coin = rng.random(n) < 0.5
print("coin flip:", round((y0[coin] + effect).mean() - y0[~coin].mean(), 2))   # ≈ 2