4. Experiments & Linear Regression

Causal Data Science for Business Analytics

Christoph Ihl

Hamburg University of Technology

Monday, 4. May 2026

Randomized Experiments

Randomization is the only design that buys you independence from every covariate — observed and unobserved.

Open the session by establishing why we are spending an entire 90-minute lecture on the simplest causal-inference design. Randomization is the only research design that creates statistical independence between the treatment indicator and every potential confounder, observed or not. Every later session in the course (matching, IV, DiD, synthetic controls) is in some sense a workaround for situations where randomization is infeasible. Students should leave this section with three lenses on why randomization works: a do-calculus lens (slides 2–3), a potential-outcomes bias decomposition (slides 4–6), and the exchangeability statement that closes both arguments (slide 7).

Randomization makes the treatment indicator independent of all covariates

Under a randomized assignment protocol, each unit’s treatment value is drawn from a known distribution that does not depend on any of the unit’s other characteristics:

\[P(D_i = d) = \frac{1}{k} \quad\text{for } k \text{ equally-sized treatment groups.}\]

Definition: Covariate Balance

We have covariate balance if the distribution of covariates \(\mathbf{X}\) is the same across treatment groups:

\[P(\mathbf{X} \mid D = d) \stackrel{d}{=} P(\mathbf{X} \mid D = d') \quad \text{for all } d, d'.\]

Equivalently:

\[\mathbf{X} \perp\!\!\!\perp D.\]

Randomization implies covariate balance — across all covariates, including unobserved ones. This is the central reason RCTs dominate observational designs.

The defining feature of an RCT is that the treatment is drawn from a coin flip (or a multi-armed analogue), so by construction it is independent of every covariate, observed or not. The “unobserved” qualifier is the crucial one: matching, regression adjustment, and IPW can only balance covariates we have measured; randomization balances them all in expectation. Students often think of RCTs as “fairness” or “ethics” devices — push them toward the statistical role: it is the design that gives us \(X \perp D\) for free.

Perspective 1a: Backdoor adjustment collapses to a marginal under randomization

Let \(\mathbf{X}\) be a sufficient adjustment set including all (un-)observed confounders. Backdoor adjustment gives:

\[P(Y \mid do(D = d)) = \sum_x P(Y \mid D = d,\, \mathbf{X} = \mathbf{x}) P(\mathbf{X} = \mathbf{x}).\]

Multiply by \(\frac{P(D = d \mid \mathbf{X})}{P(D = d \mid \mathbf{X})}\) to recover the joint distribution in the numerator:

\[P(Y \mid do(D = d)) = \sum_x \frac{P(Y \mid D = d, \mathbf{X}) \, P(D = d \mid \mathbf{X}) \, P(\mathbf{X})}{P(D = d \mid \mathbf{X})} = \sum_x \frac{P(Y, D, \mathbf{X})}{P(D = d \mid \mathbf{X})}.\]

Apply randomization \(\mathbf{X} \perp\!\!\!\perp D\), so \(P(D \mid \mathbf{X}) = P(D)\):

\[P(Y \mid do(D = d)) = \sum_x \frac{P(Y, D, \mathbf{X})}{P(D)} = \sum_x P(Y, \mathbf{X} \mid D = d) = P(Y \mid D = d). \quad\blacksquare\]

Intervention distribution collapses to the conditional — association is causation under randomization.

This is the do-calculus payoff slide. We proved the backdoor formula in Session 03; this slide shows the special case where randomization simplifies it to a single observable quantity. The algebraic moves are: (1) start from the backdoor formula; (2) introduce a 1 = P(D|X)/P(D|X) factor to assemble the joint; (3) randomization implies P(D|X) = P(D) so we can pull it out of the sum; (4) the marginalization collapses. Students should see that the do-operator’s purpose was to define the causal target — randomization makes the observational distribution equal to the causal target, no adjustment needed.

Perspective 1b: Randomization blocks every backdoor path — empty adjustment set suffices

Graphically: randomization replaces every incoming arrow into \(D\) with a single edge from a coin-flip node \(d\). All backdoor paths are mutilated, so the empty set \(\emptyset\) satisfies the backdoor criterion.

Pre-intervention world

\(X\) confounds \(D \to Y\): backdoor path \(D \leftarrow X \to Y\) is open.

Post-intervention world (do-mutilation)

The arrow \(X \to D\) is deleted: every backdoor path is blocked.

Two graphs side by side make the do-mutilation lemma concrete. In the pre-intervention world, X confounds D and Y. In the post-intervention world, the intervention replaces the natural mechanism for D with a coin-flip node, deleting all incoming arrows into D. The backdoor path D <- X -> Y becomes D … X -> Y (the leftmost edge gone), which is no longer a path because there is no edge connecting D back to X. So the empty conditioning set satisfies the backdoor criterion — formally what Perspective 1a proved algebraically. Students should connect this to S03’s do-mutilation derivation: randomization is just an experimental implementation of do(D=d).

Perspective 2: Decomposing the associational ATE

The causal ATE is defined on potential outcomes: \(ATE_{cs} = \mathbb{E}[Y_i(1) - Y_i(0)]\).

Let \(p = P(D_i = 1)\). Decompose by conditioning on \(D\):

\[ATE_{cs} = p \underbrace{\mathbb{E}[Y_i(1) - Y_i(0) \mid D_i = 1]}_{\text{ATT}} + (1-p) \underbrace{\mathbb{E}[Y_i(1) - Y_i(0) \mid D_i = 0]}_{\text{ATU}}\]

Expanding the conditional expectations and naming the four group means:

\[ \begin{aligned} ATE_{cs} &= p\, \mathbb{E}[Y_i(1) \mid D_i = 1] - p\, \mathbb{E}[Y_i(0) \mid D_i = 1] + (1-p)\, \mathbb{E}[Y_i(1) \mid D_i = 0] - (1-p)\, \mathbb{E}[Y_i(0) \mid D_i = 0] \\ &= p\, \mu_{11} - p\, \mu_{01} + (1-p)\, \mu_{10} - (1-p)\, \mu_{00}, \quad \text{where } \mu_{yd} = \mathbb{E}[Y_i(y) \mid D_i = d]. \end{aligned} \]

The associational ATE uses observed quantities only:

\[ATE_{as} = \mathbb{E}[Y_i \mid D_i = 1] - \mathbb{E}[Y_i \mid D_i = 0] = \mu_{11} - \mu_{00}.\]

The second equality uses consistency: when \(D_i = d\), the observed \(Y_i\) equals the potential outcome \(Y_i(d)\), so \(\mathbb{E}[Y_i \mid D_i = d] = \mathbb{E}[Y_i(d) \mid D_i = d]\).

Set up the bias decomposition by carefully naming the four counterfactual group means. The notation mu_{yd} reads “the mean of potential outcome Y(y) among units with treatment D=d”. Two of these four (mu_{11} and mu_{00}) are observable in any sample; the other two (mu_{01} and mu_{10}) require knowing what would have happened under the counterfactual treatment. The associational ATE only uses the two observable means. The next slide will show this differs from the causal ATE by exactly two bias terms — confounding bias and heterogeneity bias. The full algebraic add-and-subtract derivation is in the session notes; on the slide we just state the setup.

Perspective 2: Decomposing the associational ATE

Working through the algebra (see session notes for derivation), the associational ATE decomposes as:

\[ \begin{aligned} ATE_{as} \;=\;\; & \underbrace{\mathbb{E}[Y_i(1) - Y_i(0)]}_{ATE_{cs}} \\ \;+\;\; & \underbrace{\mathbb{E}[Y_i(0) \mid D_i = 1] - \mathbb{E}[Y_i(0) \mid D_i = 0]}_{\color{#00C1D4}{\text{Confounding Bias}}} \\ \;+\;\; & \underbrace{(1 - \mathbb{E}[D_i]) \big[\underbrace{\mathbb{E}[Y_i(1) \mid D_i = 1] - \mathbb{E}[Y_i(0) \mid D_i = 1]}_{\text{ATT}} - \underbrace{(\mathbb{E}[Y_i(1) \mid D_i = 0] - \mathbb{E}[Y_i(0) \mid D_i = 0])}_{\text{ATU}}\big]}_{\color{#FF7E15}{\text{Heterogeneity Bias}}} \end{aligned} \]

This is the gold equation of the session — we will use it again in Exercise 1 (Christmas / socks). Two distinct biases separate the observed mean difference from the true causal ATE:

• Confounding bias — treated and untreated differ in their baseline potential outcome \(Y(0)\).
• Heterogeneity bias — effect itself differs across treated and untreated subpopulations (ATT \(\neq\) ATU).

This is the load-bearing identity of the entire session. The associational ATE — what a naive mean comparison gives you — equals the causal ATE plus two extra terms. Confounding bias arises when the baseline (untreated) potential outcome differs across the treated and untreated subpopulations: in observational data, healthier people often choose the treatment, so mu_01 > mu_00 and the mean difference is biased downward (or upward, depending on direction). Heterogeneity bias is more subtle: even if the baselines were equal, if the treatment effect itself differs across subpopulations (ATT != ATU), and if the populations have different sizes ((1-E[D]) != 0), then the associational ATE picks up a weighted contrast of those effects. The next slide will show that under exchangeability — i.e., under randomization — both terms are exactly zero. Note the two terms are TIE-coloured inside MathJax via \color{} — \<mark\> only renders in HTML, not inside display math.

Perspective 2: Exchangeability zeroes both bias terms

Exchangeability under randomization: \(\{Y(1), Y(0)\} \perp\!\!\!\perp D\), which implies:

• Before treatment: \(\mathbb{E}[Y_i(0) \mid D_i = 0] = \mathbb{E}[Y_i(0) \mid D_i = 1] = \mathbb{E}[Y_i(0)]\).
• After treatment: \(\mathbb{E}[Y_i(1) \mid D_i = 1] = \mathbb{E}[Y_i(1) \mid D_i = 0] = \mathbb{E}[Y_i(1)]\).

\[ \begin{aligned} ATE_{as} \;=\;\; & \underbrace{\mathbb{E}[Y_i(1) - Y_i(0)]}_{ATE_{cs}} \\ \;+\;\; & \underbrace{\mathbb{E}[Y_i(0) \mid D_i = 1] - \mathbb{E}[Y_i(0) \mid D_i = 0]}_{\color{#00C1D4}{\text{Confounding Bias} \;=\; 0}} \\ \;+\;\; & \underbrace{(1 - \mathbb{E}[D_i]) \big[\mathbb{E}[Y_i(1) \mid D_i = 1] - \mathbb{E}[Y_i(0) \mid D_i = 1] - (\mathbb{E}[Y_i(1) \mid D_i = 0] - \mathbb{E}[Y_i(0) \mid D_i = 0])\big]}_{\color{#FF7E15}{\text{Heterogeneity Bias} \;=\; 0}} \end{aligned} \]

\[\boxed{ATE_{as} \;=\; ATE_{cs} \quad \text{under randomization.}}\]

Closing the bias-decomposition arc. Randomization implies the strong form of exchangeability — the joint potential-outcome vector is independent of treatment, so the conditional expectations of Y(0) and Y(1) given D do not depend on D at all. Both bias terms in the previous slide vanish, and the associational ATE equals the causal ATE. This is the same conclusion as Perspective 1 — but reached entirely through potential outcomes rather than through the do-calculus. Two perspectives, one conclusion: the associational mean difference is the causal mean difference, no adjustment needed. The boxed equation is what the rest of the session will operationalize via OLS and statistical inference.

Blocked randomization: enforce balance within strata

Simple randomization gives balance in expectation — but in any one realisation, treated and control groups can drift apart on important pre-treatment covariates (especially in small samples).

Block (stratified) randomization: partition units into strata \(S_1, \ldots, S_K\) on a pre-treatment covariate, then randomize within each stratum. Balance on \(S\) is now guaranteed exactly, not just in expectation.

Stratified ATE estimator

Within each block \(k\), compute the local mean difference \(\hat\tau_k = \bar Y_{1,k} - \bar Y_{0,k}\). Combine using stratum weights \(\pi_k = n_k / n\):

\[\hat\tau_{\text{strat}} \;=\; \sum_{k=1}^{K} \pi_k\, \hat\tau_k.\]

Variance is reduced when within-stratum outcome variance is smaller than the unconditional outcome variance — which is exactly the case when \(S\) is predictive of \(Y\).

When to block: prior-period outcome (\(Y_{\text{pre}}\)), demographic strata, geographic regions, or any pre-treatment variable known to be strongly predictive of \(Y\). Avoid blocking on too many strata (sparse cells, lost variance reduction).

Blocking is the schedule.yml topic that the legacy deck completely omitted. The intuition: simple randomization fails to balance in any single realisation; stratified randomization fixes balance on whichever covariate you condition on. The estimator is a precision-weighted combination of within-stratum effects. Variance reduction is the main reason to block — it’s the design-stage version of the variance reduction we’ll see at the analysis stage via covariate-adjusted regression. The trade-off is that with too many strata you get sparse cells (some strata might have only treated or only controls), so practitioners typically block on 1-3 categorical covariates with predictable strong effects on Y. This matters in tech A/B tests: blocking on prior-period spend or app version is standard practice. The next slide shows a worked example.

Worked example: blocking on prior-period spend in a marketing experiment

Simulated experiment with \(n = 1{,}000\) users, pre-period spend \(X\) predicts outcome spend \(Y\), and we compare simple vs. blocked random assignment.

• R
• Python

library(randomizr)
library(dplyr)

set.seed(42)
n <- 1000
X <- rnorm(n, mean = 50, sd = 15)
block_id <- cut(X, breaks = quantile(X, probs = 0:4/4),
                labels = paste0("Q", 1:4), include.lowest = TRUE)

# Simple random assignment vs blocked random assignment
D_simple  <- complete_ra(N = n, prob = 0.5)
D_blocked <- block_ra(blocks = block_id, prob = 0.5)

# True effect = 5, plus pre-period spend predicts Y strongly
Y <- 0.8 * X + 5 * D_blocked + rnorm(n, sd = 5)

# SE of the ATE under each design (HC1-robust, matches the rest of the deck)
library(estimatr)

se_simple  <- lm_robust(Y ~ D_simple,             se_type = "HC1")$std.error["D_simple"]
se_blocked <- lm_robust(Y ~ D_blocked + block_id, se_type = "HC1")$std.error["D_blocked"]

c(SE_simple = round(unname(se_simple), 3), SE_blocked = round(unname(se_blocked), 3))

 SE_simple SE_blocked 
     0.834      0.434 

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(42)
n = 1000
X = rng.normal(loc=50, scale=15, size=n)
block_id = pd.qcut(X, q=4, labels=["Q1", "Q2", "Q3", "Q4"])

# Simple random assignment
D_simple = rng.choice([0, 1], size=n)

# Blocked random assignment: shuffle within each quartile
D_blocked = np.zeros(n, dtype=int)
for q in block_id.unique():
    idx = np.where(block_id == q)[0]
    half = len(idx) // 2
    perm = rng.permutation(idx)
    D_blocked[perm[:half]] = 1

# True effect = 5
Y = 0.8 * X + 5 * D_blocked + rng.normal(scale=5, size=n)

df = pd.DataFrame({"Y": Y, "D_blocked": D_blocked,
                   "D_simple": D_simple, "block": block_id})

# SE of the ATE under each design (HC1-robust, matches the rest of the deck)
se_simple  = smf.ols("Y ~ D_simple",              data=df).fit(cov_type="HC1").bse["D_simple"]
se_blocked = smf.ols("Y ~ D_blocked + C(block)",  data=df).fit(cov_type="HC1").bse["D_blocked"]

print({"SE_simple": round(se_simple, 3), "SE_blocked": round(se_blocked, 3)})

{'SE_simple': np.float64(0.82), 'SE_blocked': np.float64(0.428)}

The blocked design typically yields a substantially smaller standard error when \(X\) is predictive of \(Y\) — variance reduction by design.

Walk through the simulation results. With X strongly predictive of Y (slope 0.8), blocking on quartiles of X reduces residual variance and hence the SE of the ATE estimator. R uses randomizr (the canonical experimental-design package); Python rolls the same logic with numpy.random.permutation per stratum. Both panels produce ballpark-equivalent SE comparisons. Students should walk away with: blocking is a design-stage variance reduction technique that pairs with analysis-stage covariate adjustment (covered in the Including covariates slide of the Extensions stack). In practice, modern A/B test platforms like Optimizely and GrowthBook block on user attributes by default.

Effect Identification

From a population identity to an estimator: mean comparison and linear regression are the same thing.

Transition into Identification stack. So far we have shown — three ways — that under randomization the associational mean difference equals the causal ATE. Identification is the formal name for that step: expressing the unobservable causal target in terms of population-level observable quantities. The estimator (the sample-level recipe) is the natural plug-in of the population identity. We then show that mean-comparison and linear regression of Y on D produce the same estimator under a binary treatment — which is why OLS will be the workhorse method for the rest of the session.

Identification = expressing the unobservable ATE in observable quantities

Identification is the step of expressing the inherently unobservable ATE in terms of observable population quantities — by appeal to randomization or further assumptions that remove confounding and heterogeneity bias.

Under randomization of \(D\), the mean comparison identifies the ATE:

\[\tau_{ATE} = \mathbb{E}[Y_i(1) - Y_i(0)] \;=\; \mathbb{E}[Y_i \mid D_i = 1] - \mathbb{E}[Y_i \mid D_i = 0].\]

The corresponding sample estimator is the difference of two sample means:

\[\hat\tau_{ATE} \;=\; \frac{1}{\sum_i D_i} \sum_i D_i Y_i \;-\; \frac{1}{\sum_i (1 - D_i)} \sum_i (1 - D_i) Y_i.\]

Make the population-vs-sample distinction explicit. The first equation is an identification result: it tells us which population quantity we are after. The second is the natural estimator: a sample analog of the population identity. Under randomization, the difference in sample means is unbiased, consistent, and asymptotically normal for the ATE — properties we will spell out in Stack 3. Students should see the parallel: in S03 we identified causal effects as adjustment-formula sums; under randomization the adjustment formula collapses to a simple mean difference.

Job Corps RCT: 9,240 disadvantaged 16–24 year olds, weekly earnings in year 4

Job Corps is a US Department of Labor education program targeting disadvantaged 16–24 year olds. The RCT randomized program access between November 1994 and February 1996; outcome is weekly earnings four years post-assignment.

• R
• Python

library(readr)
library(here)

jc <- read_csv(here("materials/data/jc.csv"), show_col_types = FALSE)

# Mean earnings, by assignment group
mean(jc$earny4[jc$assignment == 0])  # control

[1] 197.9258

mean(jc$earny4[jc$assignment == 1])  # treated

[1] 213.981

# ATE
mean(jc$earny4[jc$assignment == 1]) - mean(jc$earny4[jc$assignment == 0])

[1] 16.05513

import pandas as pd

jc = pd.read_csv("../../materials/data/jc.csv")

# Mean earnings, by assignment group
print(jc.loc[jc.assignment == 0, "earny4"].mean())  # control

197.92583695570693

print(jc.loc[jc.assignment == 1, "earny4"].mean())  # treated

213.9809642758061

# ATE
print(jc.loc[jc.assignment == 1, "earny4"].mean()
      - jc.loc[jc.assignment == 0, "earny4"].mean())

16.05512732009916

Sources: Schochet, Burghardt & Glazerman (2001); Schochet, Burghardt & McConnell (2008). Data: 1994–96 US, distributed via the causalweight R package.

Job Corps remains the canonical RCT example in causal-inference textbooks because it has all the features: large sample (n ~ 9,240), binary random assignment, well-measured outcomes, and a non-trivial effect size. The naive mean comparison gives roughly $16/week additional earnings four years out — a positive and statistically meaningful effect. We’ll re-estimate this via OLS and HC-robust SEs at the end of Stack 3. Note the disclaimer: the data is 1994–96, so the effect-size magnitudes are not directly applicable to today’s labor market. The estimator and methodology, however, are unchanged.

CATE is identified the same way under conditional independence

A randomized experiment also delivers conditional independence for any pre-treatment covariate \(X\):

\[\{Y(1), Y(0)\} \perp\!\!\!\perp D \mid X.\]

Hence the Conditional Average Treatment Effect (CATE) is identified:

\[\tau_{CATE}(x) \;=\; \mathbb{E}[Y_i(1) - Y_i(0) \mid X_i = x] \;=\; \mathbb{E}[Y_i \mid D_i = 1, X_i = x] - \mathbb{E}[Y_i \mid D_i = 0, X_i = x].\]

X must be pre-treatment

The conditioning set \(X\) must be measured before treatment assignment — formally, we condition on \(X\), not \(X(D)\). Conditioning on a post-treatment variable can both block part of the causal effect and open collider paths through unobservables (back to Session 02’s collider lesson — and previewed again in the Pre- vs. post-treatment covariates slide).

Estimators of CATEs (causal forests, meta-learners) are the topic of Session 09: Treatment Effect Heterogeneity.

CATE generalises ATE to subgroup-level effects. The randomization buys conditional exchangeability “for free” — for any pre-treatment X. The identifying formula is just the mean comparison, applied within each level of X. In practice this is what stratified analysis looks like: estimate the effect within each demographic cell, then aggregate (or report cell-by-cell). The hard pre-treatment caveat is critical and we’ll revisit it in the Pre- vs. post-treatment covariates slide: conditioning on something the treatment caused (e.g., engagement metrics measured after assignment) introduces both mediator-blocking and collider-opening bias. The CATE estimation problem proper — when X is high-dimensional and the function tau(x) is unknown — is Session 09’s territory.

From mean comparison to OLS — derivation

Start from the observed-outcome equation for each unit:

\[Y_i = Y_i(1) \cdot D_i + Y_i(0) \cdot (1 - D_i) = Y_i(0) + (Y_i(1) - Y_i(0))\, D_i.\]

Take the conditional expectation given \(D_i\) in the population:

\[\mathbb{E}[Y_i \mid D_i] = \mathbb{E}[Y_i(0) \mid D_i] + \big(\mathbb{E}[Y_i(1) \mid D_i] - \mathbb{E}[Y_i(0) \mid D_i]\big) \, D_i.\]

Apply exchangeability — under randomization, \(\{Y(1), Y(0)\} \perp\!\!\!\perp D\), so the conditional expectations do not depend on \(D\):

\[\mathbb{E}[Y_i(0) \mid D_i] = \mathbb{E}[Y_i \mid D_i = 0], \qquad \mathbb{E}[Y_i(1) \mid D_i] = \mathbb{E}[Y_i \mid D_i = 1].\]

Substituting:

\[\mathbb{E}[Y_i \mid D_i] \;=\; \underbrace{\mathbb{E}[Y_i \mid D_i = 0]}_{\beta_0} \;+\; \underbrace{\big(\mathbb{E}[Y_i \mid D_i = 1] - \mathbb{E}[Y_i \mid D_i = 0]\big)}_{\beta_1} \cdot D_i.\]

This is the population linear-regression model with \(\beta_1 = ATE\).

Walk students through the algebra carefully. Step 1 uses the consistency assumption (Y_i = Y_i(D_i)). Step 2 takes conditional expectations. Step 3 — the explicitly added step per the architect’s improve-note — invokes exchangeability to drop the conditioning on D from the right-hand-side terms. The result is the population OLS model E[Y|D] = beta_0 + beta_1 D with beta_1 exactly equal to the ATE. This is the key bridge from the potential-outcomes language to the regression workhorse: under randomization, OLS of Y on D is not just a statistical tool, it is the natural estimator of the causal effect. The next slide visualizes this with a 2-group scatter.

Residuals visualised: \(\beta_0\), \(\beta_1\), \(\epsilon_i\) on a 2-group scatter

The residual is the gap between the observed outcome and the conditional mean: \(\epsilon_i \;=\; Y_i - \underbrace{(\beta_0 + \beta_1 D_i)}_{\mathbb{E}[Y_i \mid D_i]}\).

• R
• Python

library(ggplot2)
library(dplyr)

set.seed(42)
n <- 50
df <- bind_rows(
  data.frame(D = 0, Y = rnorm(n, mean = 2, sd = 1)),
  data.frame(D = 1, Y = rnorm(n, mean = 4, sd = 1))
)

means <- df |> group_by(D) |> summarise(Ybar = mean(Y), .groups = "drop")
mu0 <- means$Ybar[means$D == 0]
mu1 <- means$Ybar[means$D == 1]

ggplot(df, aes(x = D, y = Y)) +
  geom_point(shape = 1, alpha = 0.7) +
  geom_segment(aes(x = 0, y = mu0, xend = 1, yend = mu1),
               colour = "#005e73", linewidth = 0.6) +
  geom_hline(yintercept = mu0, linetype = "dotted") +
  geom_hline(yintercept = mu1, linetype = "dotted") +
  geom_segment(aes(x = 0.99, y = mu1, xend = 0.99, yend = max(df$Y)),
               colour = "#FF7E15", linewidth = 0.8,
               arrow = arrow(length = unit(0.18, "cm"))) +
  annotate("text", x = 0.94, y = 4.8,
           label = "epsilon[i]", parse = TRUE,
           colour = "#FF7E15", size = 5) +
  geom_segment(aes(x = 0.5, y = mu0, xend = 0.5, yend = mu1),
               colour = "#00C1D4", linewidth = 0.8,
               arrow = arrow(length = unit(0.18, "cm"), ends = "both")) +
  annotate("text", x = 0.55, y = 3,
           label = "beta[1]", parse = TRUE,
           colour = "#00C1D4", size = 5) +
  geom_segment(aes(x = 0.02, y = 0, xend = 0.02, yend = mu0),
               colour = "#00C1D4", linewidth = 0.8,
               arrow = arrow(length = unit(0.18, "cm"))) +
  annotate("text", x = 0.07, y = 1,
           label = "beta[0]", parse = TRUE,
           colour = "#00C1D4", size = 5) +
  scale_x_continuous(breaks = c(0, 1), limits = c(-0.05, 1.05)) +
  scale_y_continuous(breaks = 0:6, limits = c(0, 6)) +
  labs(x = "D", y = "Y")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
n = 50
df = pd.DataFrame({
    "D": np.r_[np.zeros(n), np.ones(n)],
    "Y": np.r_[rng.normal(2, 1, n), rng.normal(4, 1, n)],
})
mu0 = df.loc[df.D == 0, "Y"].mean()
mu1 = df.loc[df.D == 1, "Y"].mean()

fig, ax = plt.subplots(figsize=(7, 4))
_ = ax.scatter(df["D"], df["Y"], facecolors="none", edgecolors="black", alpha=0.7)
_ = ax.plot([0, 1], [mu0, mu1], color="#005e73", linewidth=1.4)
_ = ax.axhline(mu0, linestyle=":", color="grey")
_ = ax.axhline(mu1, linestyle=":", color="grey")
_ = ax.annotate("", xy=(0.99, df["Y"].max()), xytext=(0.99, mu1),
                arrowprops=dict(arrowstyle="->", color="#FF7E15", lw=1.5))
_ = ax.text(0.94, 4.8, r"$\epsilon_i$", color="#FF7E15", fontsize=14)
_ = ax.annotate("", xy=(0.5, mu1), xytext=(0.5, mu0),
                arrowprops=dict(arrowstyle="<->", color="#00C1D4", lw=1.5))
_ = ax.text(0.55, 3, r"$\beta_1$", color="#00C1D4", fontsize=14)
_ = ax.annotate("", xy=(0.02, mu0), xytext=(0.02, 0),
                arrowprops=dict(arrowstyle="->", color="#00C1D4", lw=1.5))
_ = ax.text(0.07, 1, r"$\beta_0$", color="#00C1D4", fontsize=14)
_ = ax.set(xticks=[0, 1], xlim=(-0.05, 1.05), ylim=(0, 6),
           xlabel="D", ylabel="Y")
plt.tight_layout()
plt.show()

The 2-group scatter visualizes the regression objects geometrically. beta_0 is the conditional mean in the control group; beta_0 + beta_1 is the conditional mean in the treated group; beta_1 is the vertical gap between the two horizontal mean lines. epsilon_i is the residual — the vertical distance from each observation to the relevant group mean. With binary D, OLS literally fits two horizontal lines, one per group; the slope beta_1 is the difference between the two heights. This is the geometric statement of “OLS of Y on D = mean comparison”.

OLS = minimising the sum of squared residuals

The sample estimator is the OLS minimiser:

\[\hat\beta_0,\, \hat\beta_1 \;=\; \arg\min_{\beta_0^*,\, \beta_1^*} \sum_{i=1}^n \underbrace{(Y_i - \beta_0^* - \beta_1^* D_i)^2}_{\hat\epsilon_i^2}.\]

Solving the first-order conditions yields the closed-form estimators:

\[\hat\beta_1 \;=\; \frac{\widehat{\text{Cov}}(Y_i, D_i)}{\widehat{\text{Var}}(D_i)},\]

\[\widehat{\text{Cov}}(Y_i, D_i) = \tfrac{1}{n-1} \sum_{i=1}^n \left(Y_i - \tfrac{1}{n}\sum_j Y_j\right)\left(D_i - \tfrac{1}{n}\sum_j D_j\right),\]

\[\widehat{\text{Var}}(D_i) = \tfrac{1}{n-1} \sum_{i=1}^n \left(D_i - \tfrac{1}{n}\sum_j D_j\right)^2, \qquad \hat\beta_0 \;=\; \bar{Y} - \hat\beta_1 \bar{D}.\]

The closed-form solution. Beta_1-hat equals the sample covariance of Y and D divided by the sample variance of D. With binary D this simplifies algebraically to the difference of group means — a useful sanity-check exercise. Beta_0-hat is then chosen so the regression line passes through the sample mean point (D-bar, Y-bar). The legacy slide had a sign error on beta_0-hat (a leading minus sign on Y-bar); the correct formula is Y-bar minus beta_1-hat times D-bar. Students don’t need to memorise these — modern OLS software handles the closed form — but they should recognise that linear regression is doing pure linear algebra.

Statistical Inference

A point estimate is not a finding — uncertainty quantification is what turns a number into a claim.

Transition into Inference. Identification gave us a recipe for the population ATE; estimation gave us a sample analog. Statistical inference is the third leg of the stool: how confident can we be that the sample estimate is close to the population truth? We’ll cover three flavors: (1) the BLUE properties — what OLS guarantees in large samples; (2) Neyman vs Fisher inference — two philosophies for quantifying uncertainty; (3) the practical machinery — robust SEs, t-tests, confidence intervals, and the bootstrap. The output of Inference is what you actually report in the paper or the dashboard.

Under random assignment, OLS estimators are unbiased, consistent, asymptotically normal

Under standard regularity conditions (Wooldridge 2010), the OLS estimators of \(\beta_0, \beta_1\) have three desirable properties:

Unbiased

On average across hypothetical samples, equal to the true parameter: \(\mathbb{E}[\hat\beta_1] = \beta_1, \;\; \mathbb{E}[\hat\beta_0] = \beta_0\).

Consistent

Converges in probability to the truth as \(n \to \infty\): \(\hat\beta_1 \xrightarrow{p} \beta_1, \;\; \hat\beta_0 \xrightarrow{p} \beta_0\).

Asymptotically normally distributed

Approximately normal in large samples — the basis for \(t\)-tests and confidence intervals:

\[\sqrt{n}(\hat\beta_1 - \beta_1) \xrightarrow{d} \mathcal{N}(0,\, \sigma^2_{\beta_1}), \qquad \sqrt{n}(\hat\beta_0 - \beta_0) \xrightarrow{d} \mathcal{N}(0,\, \sigma^2_{\beta_0}).\]

Proofs: Wooldridge (2010), Econometric Analysis of Cross Section and Panel Data, Ch. 4. Skipped in lecture.

Three properties of OLS that justify everything we do downstream. Unbiasedness is the small-sample property — the expected value of the estimator equals the parameter. Consistency is the large-sample property — the estimator converges to the parameter as n grows. Asymptotic normality is the distributional property — for large enough n, the sampling distribution of beta_1-hat is approximately normal, which lets us construct t-tests and confidence intervals using normal-distribution theory. Proofs require some moment-condition algebra (E[Xu] = 0 etc.) — we skip them; Wooldridge has the canonical treatment.

In an RCT, who actually randomizes? The coin flip — not nature drawing from an infinite population.

Plant the wedge before showing the slide. Standard CLT inference assumes the sample is one i.i.d. draw from an infinite superpopulation; in an RCT, the only thing that is actually random is the treatment-assignment vector among the specific units in front of us. This is the conceptual flip that motivates design-based inference. Two traditions answer it: Neyman averages over hypothetical re-randomizations and takes n large; Fisher enumerates re-randomizations and stays exact. Different stances on the same randomness; different sweet spots.

Neyman vs. Fisher inference: large-sample asymptotics vs. exact randomization tests

Both fix the units in front of you and treat the assignment vector \(\mathbf{D}\) as the only randomness — they disagree on what to do with it.

Neyman (1923) — large-sample

Asks “how big is the average effect, and how precisely?” — tests the weak null \(\mathbb{E}[Y(1) - Y(0)] = 0\). Treats the potential outcomes as fixed; randomization is the source of variability.

\[\sqrt{n}(\hat\tau - \tau) \xrightarrow{d} \mathcal{N}(0,\, V_N).\]

Output: point estimate \(\hat\tau\) + asymptotic SE + CI. Variance \(V_N\) is conservatively estimated under treatment-effect heterogeneity (Neyman’s variance bound).

Pros: closed-form CIs, scales to large \(n\), dominates textbooks. Cons: only valid asymptotically; conservative in small samples.

Fisher (1925) — exact

Asks “could this difference be a coincidence under no effect for anyone?” — tests the sharp null \(H_0: Y_i(1) = Y_i(0)\) for every unit \(i\).

Under \(H_0\), observed outcomes are fixed; the test-statistic distribution comes from re-shuffling treatment labels across all \(\binom{n}{n_1}\) permutations.

Output: an exact finite-sample p-value \[p = \frac{1}{|\mathcal{P}|} \sum_{\pi \in \mathcal{P}} \mathbf{1}\{|\hat\tau_\pi| \geq |\hat\tau_{\text{obs}}|\}.\]

Pros: no asymptotic approximation; works in tiny samples. Cons: tests sharp null (not average null); compute-intensive.

The two inference traditions for randomized experiments. Neyman’s view fixes the potential outcomes and treats randomization as the source of sampling variation; the resulting variance estimator is asymptotic and conservatively biased upward when treatment effects are heterogeneous. Fisher’s view tests the sharp null that the treatment had zero effect on every single unit — under that null, swapping treatment labels around shouldn’t change the test statistic distribution, which gives an exact p-value via the randomization distribution itself. The Fisher approach is increasingly the default in modern tech-company A/B testing because it gives valid p-values even with n in the dozens, which Neyman-style asymptotics struggle with. The trade-off: Fisher tests a stronger null (no effect for any unit, rather than no average effect), and it is computationally more expensive (need to enumerate or sample many permutations). For most lecture purposes, Neyman is the working framework, and we use it for the rest of this stack.

From here on we ride with Neyman — but Fisher’s permutation logic returns next session as the engine behind A/B-test p-values.

Resolve the slide before moving on. Neyman is our working framework for the rest of this lecture (closed-form SEs, t-tests, CIs) — that’s what the next slide builds. But flag the forward link: Fisher’s randomization-test idea is exactly what powers permutation-based p-values in modern A/B-test platforms, which we revisit in Session 5 alongside CUPED and stratified A/B testing. The framing students should leave with: two legitimate stances on the same randomness, picked by sample size and computational budget, not by ideology.

Inference goals: standard errors, hypothesis tests, confidence intervals

A point estimate \(\hat\beta_1\) has zero information about its own precision. We need three primitives:

• Standard errors — \(se(\hat\beta_1)\) measures the sample-to-sample variability of the estimator.

• Hypothesis tests — assess whether the data are inconsistent with a specific null value (typically \(\beta_1 = 0\)).

• Confidence intervals — report a range of plausible values for \(\beta_1\) at a chosen confidence level.

These three are not three separate calculations: they are three views on the same sampling-distribution object. Once you have the SE, you have everything else by simple arithmetic.

A point estimate without a measure of precision is statistically meaningless — you don’t know whether 0.3 is “exactly 0.3” or “0.3 plus or minus 5”. The three inference primitives are different ways of summarising the same underlying object: the sampling distribution of the estimator. SEs report its spread directly. Hypothesis tests ask “is the value zero (or some other reference) compatible with our data?”. CIs invert hypothesis tests to report a range of compatible parameter values. Once you have the SE you can compute t-statistics, p-values, and CI endpoints with grade-school arithmetic.

Heteroskedasticity-robust standard error of \(\hat\beta_1\)

The asymptotic variance of \(\hat\beta_1\) in the population (heteroskedasticity-robust / sandwich form):

\[\text{Var}(\hat\beta_1) \;=\; \frac{\mathbb{E}\!\left[\epsilon^2 (D_i - \mathbb{E}[D_i])^2\right]}{n \cdot (\text{Var}(D_i))^2}.\]

Plug-in sample analog using residuals \(\hat\epsilon_i = Y_i - \hat\beta_0 - \hat\beta_1 D_i\):

\[\widehat{\text{Var}}(\hat\beta_1) \;=\; \frac{\sum_i \hat\epsilon_i^2 (D_i - \bar D)^2}{\big[\sum_i (D_i - \bar D)^2\big]^2}, \qquad se(\hat\beta_1) = \sqrt{\widehat{\text{Var}}(\hat\beta_1)}.\]

This is the HC0 sandwich estimator. R’s sandwich::vcovHC(fit, type = "HC1") and Python’s cov_type="HC1" apply a small-sample correction factor of \(n/(n-k)\).

HC vs. classical SE

The expression above is the heteroskedasticity-robust (HC0/HC1) sandwich form. The classical “homoskedastic” SE simplifies to \(\sigma_\epsilon^2 / (n \cdot \text{Var}(D))\) and assumes \(\text{Var}(\epsilon \mid D) = \sigma_\epsilon^2\) — a constraint we rarely want to impose. Default to HC.

The legacy deck called this the “asymptotic variance” but the formula given is already the heteroskedasticity-robust sandwich form — the classical homoskedastic version would be much simpler. We have re-titled the slide and added the callout to flag the distinction. In practice: always use HC SEs (HC1 for finite-sample correction) — they are valid under heteroskedasticity AND under homoskedasticity, while classical SEs are valid only under homoskedasticity. The R package sandwich (Zeileis), Python’s statsmodels (cov_type=‘HC1’), and Stata’s , robust all do the same thing. The t-distribution rather than normal is needed because the SE is itself estimated; with df = n - 2 (binary D, two coefficients estimated), and for any reasonable n, t_{n-2} is essentially N(0,1).

Reading the t-distribution: 5% two-sided vs. 5% one-sided critical values

The standardised statistic with known \(sd\) is normal; with the estimated \(se\) it is \(t\)-distributed with \(df = n - 2\):

\[z_1 = \frac{\hat\beta_1 - \beta_1}{sd(\hat\beta_1)} \xrightarrow{d} \mathcal{N}(0, 1) \quad \text{vs.} \quad t_1 = \frac{\hat\beta_1 - \beta_1}{se(\hat\beta_1)} \sim t_{n-2}.\]

Two ggplots side by side make the two-sided / one-sided distinction visceral. The blue shaded area is the rejection region — when our observed t-statistic falls there under the null, we reject. For two-sided tests at alpha = 5%, the rejection region is split: 2.5% in each tail, with critical value +/- 1.96. For one-sided tests at alpha = 5%, all 5% goes in one tail, with critical value 1.65 (or -1.65 for the lower-tail version). For df = 998 (effectively n = 1000), the t-distribution is indistinguishable from N(0,1) — but we still write t-distribution because the SE is estimated, not known.

For the curious: the t-distribution is defined as \(t_{df} = Z / \sqrt{V/df}\) with \(Z \sim N(0,1)\) and \(V \sim \chi^2(df)\), \(Z \perp V\). Our \(t_1 = (\hat\beta_1 - \beta_1)/se\) matches this template: the numerator is a \(z\)-statistic, and dividing the unknown \(\sigma\) out via the estimated \(\hat\sigma\) introduces the \(\sqrt{V/df}\) in the denominator (with \(V = (n-2)\hat\sigma^2/\sigma^2 \sim \chi^2(n-2)\) under classical assumptions). Wooldridge ch. 4 has the full derivation. We state-and-use rather than derive in lecture.

Hypothesis testing — null/alternative + significance level

Step 1. State the null \(H_0\) and alternative \(H_1\):

• Two-sided (undirected): \(H_0\): \(\beta_1 = 0\) vs. \(H_1\): \(\beta_1 \neq 0\).
• One-sided (directed, positive): \(H_0\): \(\beta_1 \leq 0\) vs. \(H_1\): \(\beta_1 > 0\).
• One-sided (directed, negative): \(H_0\): \(\beta_1 \geq 0\) vs. \(H_1\): \(\beta_1 < 0\).

Step 2. Set the significance level \(\alpha\):

• The maximally accepted type I error probability — the chance of rejecting \(H_0\) when it is in fact true.
• \(\alpha = 0.05\) means at most a 5% rate of false rejections in repeated samples.
• Conventional alternatives: \(\alpha = 0.01\) (1%), \(\alpha = 0.10\) (10%).

Two-sided is the workhorse: we reject H0 if the data are far from zero in either direction. One-sided is appropriate only if directional theory rules out one tail before seeing the data — which is rare and usually a red flag for HARKing if applied post hoc. Significance level alpha is the false-rejection rate we are willing to tolerate; alpha = 0.05 is the convention but it has no statistical privilege over alpha = 0.01 or 0.10. For applied A/B testing, alpha = 0.05 is typical; for high-stakes scientific findings, lower alphas are common.

Hypothesis testing — critical value, t-statistic, p-value

Step 3. Compute the critical value \(c\):

• Quantile of the \(t\)-distribution at \(\alpha\), given \(H_0\).
• Two-sided, \(\alpha = 0.05\), \(df = 1000-2\): \(c = t_{\alpha/2 = 0.025,\, 998} \approx 1.96\).
• One-sided, \(\alpha = 0.05\), \(df = 998\): \(c = t_{\alpha = 0.05,\, 998} \approx 1.65\).

Step 4. Compute the test statistic and p-value:

• Test statistic: \(t_1 = \dfrac{\hat\beta_1 - 0}{se(\hat\beta_1)}\).
• p-value: probability of a test statistic as extreme as \(t_1\) under \(H_0\):
  • Two-sided: \(p = P(|A| \geq |t_1|)\) for \(A \sim t_{n-2}\).
  • One-sided: \(p = P(A \geq t_1)\) or \(p = P(A \leq t_1)\).

Decision rules:

• Two-sided: reject \(H_0\) if \(|t_1| > c\), equivalently \(p < \alpha\).
• One-sided: reject \(H_0\) if \(t_1 > c\) (or \(t_1 < -c\)), equivalently \(p < \alpha\).

Steps 3 and 4 turn the abstract H0/H1 setup into a concrete decision. The critical value c is the threshold beyond which we reject; the test statistic is the standardised distance of our estimate from the null. Reject if the test statistic exceeds the critical value — equivalently, reject if the p-value falls below the significance level. The p-value answers: “if the null were true, how likely would we be to see data at least as extreme as ours?” Small p means our data is unlikely under the null, so we reject. Note: small p does NOT mean the null is unlikely; it means our data is unlikely given the null. Misinterpreting this is the classic “prosecutor’s fallacy”.

Confidence intervals: \(\hat\beta_1 \pm c \cdot se(\hat\beta_1)\)

A \((1 - \alpha)\) confidence interval is the range of values \(\beta_1^*\) for which the corresponding null \(H_0: \beta_1 = \beta_1^*\) would not be rejected at level \(\alpha\).

Frequentist interpretation: across hypothetical repeated samples, \((1 - \alpha)\) of the constructed CIs would cover the true \(\beta_1\).

Two-sided CI:

\[CI \;=\; \hat\beta_1 \;\pm\; c \cdot se(\hat\beta_1).\]

Worked example: \(\hat\beta_1 = 0.30\), \(se = 0.05\), \(c = 1.96\) → \(CI = 0.30 \pm 1.96 \cdot 0.05 = [0.202,\; 0.398]\).

One-sided CI:

\[CI \;=\; \hat\beta_1 \;-\; c \cdot se(\hat\beta_1) \quad\text{(upper open)}.\]

Worked example with \(c = 1.65\): \(CI = [0.218,\; \infty)\).

CI is the dual of a hypothesis test. The two-sided 95% CI is the set of all null values that we would NOT reject at the 5% level given the observed data. The frequentist interpretation is subtle: the 95% probability statement is about the procedure (how often it generates intervals that cover the truth), not about any one realised interval (which either covers or doesn’t, with no probability statement attached). One-sided CIs are unbounded on one end and useful when only one direction of the parameter is substantively interesting.

Common misinterpretations of \(\alpha\) and CIs — and where the intuitive reading actually lives

The frequentist machinery is precise but not intuitive. Students (and reviewers, and journalists) routinely misread it.

What people want it to mean

• “\(p = 0.03\) means there’s a 3% chance the null is true.”
• “There’s a 95% chance the true \(\beta_1\) lies in this CI.”
• “\(\alpha = 0.05\) is the probability I’m wrong.”
• “A non-significant result means \(H_0\) is true.”

These are all probability statements about hypotheses or parameters given data, \(P(H_0 \mid \text{data})\) or \(P(\beta_1 \in CI \mid \text{data})\).

What it actually means

• \(p\): probability of data this extreme given \(H_0\) true, \(P(\text{data} \mid H_0)\).
• 95% CI: across hypothetical re-samples, 95% of the constructed intervals cover \(\beta_1\). Any single realised interval either covers or doesn’t.
• \(\alpha\): long-run type I error rate of the procedure.
• Non-significant ≠ “no effect” — only “no detectable effect at this \(n\) and \(\alpha\).”

All are statements about the procedure, not the parameter.

The intuitive readings are legitimate — under Bayesian inference. With a prior on \(\beta_1\), the posterior \(P(\beta_1 \mid \text{data})\) supports direct statements like “there’s a 95% chance \(\beta_1 \in [0.20,\, 0.40]\)”. That’s the credible interval — coming up in Session 5: online A/B testing.

A slide of warnings and a forward link. The four misreadings on the left are the canonical ones — all of them treat \(p\), CI, or \(\alpha\) as making a probability statement about the parameter or hypothesis given the data, which is exactly what frequentist inference cannot do. The right column gives the correct, procedure-level reading. The fragment closes the loop: those intuitive readings are not wrong-headed, they are just Bayesian — under a prior on \(\beta_1\), the posterior gives you exactly the kind of “95% chance the parameter is in the interval” statement everyone wants. Forward link to S05 where Bayesian A/B testing (and the credible interval) gets a proper treatment alongside CUPED + stratified A/B testing. Pedagogical purpose: send students into the next session predisposed to see Bayesian inference as fixing a real problem, not as a competing ideology.

R-squared in the single-treatment regression equals the squared correlation

General: percentage of variance in \(Y\) explained by all variables in the model:

\[R^2 \;=\; 1 - \frac{SSR}{SST}, \qquad SSR = \sum_i (Y_i - \hat Y_i)^2, \qquad SST = \sum_i (Y_i - \bar Y)^2.\]

Specific case — only a single binary treatment \(D\) in the model:

\[R^2 \;=\; r^2 \;=\; \left( \frac{\sum_i (Y_i - \bar Y)(D_i - \bar D)}{\sqrt{\sum_i (Y_i - \bar Y)^2} \cdot \sqrt{\sum_i (D_i - \bar D)^2}} \right)^{\!2}.\]

In a causal-inference context, \(R^2\) is largely diagnostic — high \(R^2\) does not mean the regression is causally identified, and low \(R^2\) does not preclude precise ATE estimation. We rarely report it as a primary result.

R-squared has the standard “fraction of variance explained” interpretation. The specific case for binary D collapses to the squared Pearson correlation between Y and D. Important interpretive caveat: in a causal-inference course, R-squared is largely diagnostic. A randomized treatment effect can be small relative to total Y variance (low R-squared) yet perfectly identified and statistically significant. Conversely, a regression with high R-squared can still be confounded if it includes the wrong covariates. Don’t optimise for R-squared in causal work — optimise for unbiased identification.

Job Corps via OLS with HC-robust SEs

Estimate the ATE of Job Corps assignment on year-4 weekly earnings using OLS with heteroskedasticity-robust standard errors.

• R
• Python

library(readr)
library(here)
library(sandwich)
library(modelsummary)

jc <- read_csv(here("materials/data/jc.csv"),
               show_col_types = FALSE)
ols <- lm(earny4 ~ assignment, data = jc)

modelsummary(ols,
  vcov = sandwich::vcovHC,
  estimate = "{estimate} (se = {std.error}, t = {statistic}){stars}",
  statistic = "p = {p.value}, CI = [{conf.low}, {conf.high}]",
  stars = TRUE,
  gof_map = c("r.squared", "nobs"),
  output = "markdown")

import pandas as pd
import statsmodels.formula.api as smf

jc = pd.read_csv("../../materials/data/jc.csv")
ols = smf.ols("earny4 ~ assignment",
              data=jc).fit(cov_type="HC1")
print(ols.summary().tables[1])

(Intercept)	197.926 (se = 3.073, t = 64.416)***
	p = <0.001, CI = [191.903, 203.949]
assignment	16.055 (se = 4.074, t = 3.941)***
	p = <0.001, CI = [8.069, 24.041]
R2	0.002
Num.Obs.	9240

Significance: + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001.

This is the standard reporting workflow for a causal RCT analysis. R uses lm() + sandwich::vcovHC for HC SEs, with modelsummary for a publication-ready table. Python uses statsmodels’ formula API with cov_type=‘HC1’ for the same SE adjustment. Both produce numerically identical (within rounding) point estimates and SEs. The reported coefficient is roughly $16/week — a small but statistically distinguishable effect on a sample of 9,240 individuals. The 95% CI excludes zero, so we reject the null of no effect at alpha = 5%. This is the analysis we’d write up in a paper.

Bootstrap: resample to estimate \(se(\hat\beta_1)\)

The sample stands in for the population: resampling with replacement mimics fresh draws, and the resulting variability proxies the true sampling variability.

The procedure (Huber 2023):

1. Draw \(B\) bootstrap samples of size \(n\) with replacement from the original sample.
2. Estimate \(\hat\beta_1^{(b)}\) in each.
3. The bootstrap SE is the standard deviation across the \(B\) replicates:

\[se(\hat\beta_1) \;=\; \sqrt{\tfrac{1}{B - 1} \sum_{b=1}^{B} (\hat\beta_1^{(b)} - \bar{\hat\beta_1})^2}.\]

• R
• Python

library(readr); library(here); library(boot)
jc <- read_csv(here("materials/data/jc.csv"),
               show_col_types = FALSE)
bs <- function(data, idx) {
  coef(lm(earny4 ~ assignment, data = data[idx, ]))[2]
}
set.seed(1)
res <- boot(data = jc, statistic = bs, R = 1999)
c(est = res$t0, boot_se = sd(res$t))

est.assignment        boot_se 
      16.05513        3.95481 

import pandas as pd, numpy as np, statsmodels.formula.api as smf
jc = pd.read_csv("../../materials/data/jc.csv")
rng = np.random.default_rng(1)
def bs(df, idx):
    return smf.ols("earny4 ~ assignment",
                   data=df.iloc[idx]).fit().params["assignment"]
B = 1999
boots = np.array([
    bs(jc, rng.integers(0, len(jc), len(jc)))
    for _ in range(B)
])
print({"est": boots.mean().round(2), "boot_se": boots.std(ddof=1).round(3)})

Bootstrap SE ≈ HC-robust SE in well-behaved samples.

When bootstrap > HC-robust: small \(n\), heavy tails, clustered/dependent data, or non-smooth statistics (medians, quantiles, ratios). Otherwise HC1 is faster and equally honest.

The bootstrap is non-parametric inference: instead of relying on asymptotic normal theory or a specific functional form for the variance, we resample from the data itself and use the empirical sampling distribution. The trick is “the empirical distribution is the population” — by resampling with replacement, we mimic drawing fresh samples from the unknown population. R uses boot::boot with B = 1999 replicates; Python rolls the same logic with numpy. For OLS with HC-robust SEs, the bootstrap SE typically matches the analytic HC SE within a few percent — they’re both estimating the spread of the same sampling distribution. The bootstrap is most useful when (a) the analytic SE is unavailable or known to be biased, and (b) for non-standard functionals (medians, quantiles, complex test statistics).

Extensions

OLS extends naturally to multivalued and continuous treatments — and to covariate-adjusted variance reduction.

Closing stack. So far D has been binary. We now extend the OLS framework to (a) multivalued treatments via dummy encoding, (b) continuous treatments via marginal effects, and (c) covariate adjustment for variance reduction. The covariate slot is the most subtle — under randomization, covariates are NOT for de-confounding (that’s S06’s job), they are for variance reduction. We close with a warning: post-treatment covariates do the opposite of variance reduction; they re-introduce bias.

Multivalued treatments: dummy-encode each level

Treatment \(D\) can take discrete values, ordered or unordered:

• Ordered: \(D \in \{0, 1, 2, \ldots\}\) (e.g., 0 / 1 / 2 weeks of training).
• Unordered: \(D \in \{A, B, C, \ldots\}\) (e.g., no training, IT, management).

The independence assumption generalises to all treatment levels:

\[\{Y_i(0), Y_i(1), Y_i(2), \ldots, Y_i(J)\} \perp\!\!\!\perp D_i.\]

Encode each non-baseline level as a binary dummy \(D_{ij} = \mathbf{1}\{D_i = j\}\) and run OLS:

\[\mathbb{E}[Y_i \mid D_i] \;=\; \underbrace{\beta_0}_{\mathbb{E}[Y \mid D = 0]} + \underbrace{\beta_1}_{\Delta \text{ vs } 0 \text{ at } 1} D_{i1} + \underbrace{\beta_2}_{\Delta \text{ vs } 0 \text{ at } 2} D_{i2} + \ldots + \underbrace{\beta_J}_{\Delta \text{ vs } 0 \text{ at } J} D_{iJ}.\]

Each \(\beta_j\) is the ATE of level \(j\) vs. the baseline level \(0\).

Multivalued treatments are common in practice: dose-response studies, multi-armed A/B tests, factorial designs. The encoding is straightforward — pick a baseline level (typically 0 or “control”), create J dummy variables for the other levels, and read off J ATEs from the regression coefficients. Each beta_j is the ATE comparing level j to the baseline level 0. Interpretation note: by default the baseline is the omitted level, so a 5-armed experiment has 4 ATEs in the regression output. To get an ATE between two non-baseline levels, take the difference of the corresponding coefficients (or re-parameterize with a different reference).

Multivalued treatments example

804 university students randomized into three mutually exclusive arms: control (default survey: professional-plans questions first, then personal-info questions), information (shown a graph of median monthly Swiss wages before being asked their own expectation), order (default sequence reversed — personal questions first, professional plans second — pure framing, no extra information). Outcome wexpect2 = own wage expectation 3 years after graduation, in 500 CHF brackets (Fernandes, Huber, and Vaccaro 2020).

• R
• Python

library(readr); library(here)
library(sandwich); library(modelsummary)
wexpect <- read_csv(here("materials/data/wexpect.csv"),
                    show_col_types = FALSE)
ols <- lm(wexpect2 ~ treatmentinformation + treatmentorder,
          data = wexpect)
modelsummary(ols,
  vcov = sandwich::vcovHC,
  estimate = "{estimate} (se = {std.error}, t = {statistic}){stars}",
  statistic = "p = {p.value}, CI = [{conf.low}, {conf.high}]",
  stars = TRUE,
  gof_map = c("r.squared", "nobs"),
  output = "markdown")

import pandas as pd, statsmodels.formula.api as smf
wexpect = pd.read_csv("../../materials/data/wexpect.csv")
ols = smf.ols(
    "wexpect2 ~ treatmentinformation + treatmentorder",
    data=wexpect).fit(cov_type="HC1")
print(ols.summary().tables[1])

(Intercept)	9.408 (se = 0.159, t = 59.268)***
	p = <0.001, CI = [9.096, 9.719]
treatmentinformation	0.345 (se = 0.243, t = 1.421)
	p = 0.156, CI = [-0.132, 0.822]
treatmentorder	-0.173 (se = 0.234, t = -0.741)
	p = 0.459, CI = [-0.633, 0.286]
R2	0.006
Num.Obs.	804

Significance: + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001.

Single 3-level treatment encoded as 2 dummies, with control as baseline. The information arm is informational (real wage data); the order arm is psychological (which question is asked first). Both have plausibly causal effects on stated wage expectations. The two-coefficient OLS reads off both ATEs vs. control simultaneously under randomization. Interactions don’t apply here because the arms are mutually exclusive (no subject is both info AND order); if instead this were a 2×2 factorial you’d add treatmentinformation:treatmentorder.

Continuous treatments: ATE = marginal effect \(\partial \mathbb{E}[Y \mid D] / \partial D\)

Treatment \(D\) takes (continuum-) many values respecting cardinality: marketing spend in EUR, years of education, dose in mg.

Independence assumption for every value of the continuous treatment:

\[Y_i(d) \perp\!\!\!\perp D_i \qquad \text{for all } d.\]

Two strategies:

1. Discretise by binning into narrow brackets — \(D_{i0} = \mathbf{1}\{D_i \leq 100\}\), \(D_{i1} = \mathbf{1}\{100 < D_i \leq 200\}\), …; reduce to multivalued case.
2. Include \(D_i\) continuously in the regression: \(\mathbb{E}[Y_i \mid D_i] = \beta_0 + \beta_1 D_i\). Then the ATE = marginal effect:

\[ATE \;=\; \frac{\partial \mathbb{E}[Y_i \mid D_i]}{\partial D_i} \;=\; \beta_1.\]

Continuous treatments require a small generalisation of the independence assumption: it must hold for every potential outcome Y(d), not just two. Under that assumption, the ATE becomes a marginal effect — the partial derivative of the conditional mean with respect to D. With a linear specification, this is just the slope coefficient beta_1. The interpretation is “the average effect on Y of a one-unit increase in D”, evaluated at any point in the support (linearity makes it constant). For a non-linear specification, the marginal effect varies with D — see the next slide.

Continuous treatments — non-linearity via quadratic specification

Linearity assumption: \(\partial \mathbb{E}[Y_i(D = t')] / \partial D_i \;=\; \partial \mathbb{E}[Y_i(D = t)] / \partial D_i\) for all \(t \neq t'\) in the support of \(D_i\) — i.e., the ATE is constant in \(D\).

Quadratic specification: \(\mathbb{E}[Y_i \mid D_i] = \beta_0 + \beta_1 D_i + \beta_2 D_i^2 \;\;\Rightarrow\;\; ATE(D_i) = \partial \mathbb{E}[Y_i \mid D_i] / \partial D_i = \beta_1 + 2 \beta_2 D_i\).

For maximum flexibility: higher-order polynomials, splines, or kernel regression.

Linearity is a strong assumption — it says the slope of E[Y|D] is constant across the support of D. The two scatterplots show the visual difference. Left: linear DGP, the line fits well and the slope is constant. Right: non-linear DGP, the loess curve bends, so the marginal effect at D=1 differs from the marginal effect at D=4. To capture this with regression, add a quadratic term — the marginal effect is then beta_1 + 2beta_2D, which varies with D. Higher-order polynomials, splines, and kernel regression are all flexible alternatives. Session 08 (Double Machine Learning) returns to flexible nuisance estimation in a more general framework.

Continuous treatments example: newspaper ad spend → sales (kernel regression)

Kernel regression estimates \(\widehat{\mathbb{E}}[Y \mid D = d]\) as a locally weighted average of nearby \(Y_i\) — observations close to \(d\) get high weight, distant ones get low weight; the bandwidth controls “close”.

200 brand-week observations, marketing-spend dataset (datarium::marketing). Outcome: sales; treatment: newspaper-ad spend (continuous, USD).

• R
• Python

library(readr); library(here); library(np)
options(np.messages = FALSE)
mkt <- read_csv(here("materials/data/marketing.csv"),
                show_col_types = FALSE)
bw <- npregbw(sales ~ newspaper, data = mkt)
plot(bw, plot.errors.method = "asymptotic")

import pandas as pd, numpy as np, matplotlib.pyplot as plt
from statsmodels.nonparametric.kernel_regression import KernelReg
mkt = pd.read_csv("../../materials/data/marketing.csv")
D = mkt["newspaper"].values
Y = mkt["sales"].values
kr = KernelReg(endog=Y, exog=D, var_type="c")
grid = np.linspace(D.min(), D.max(), 100)
yhat, _ = kr.fit(grid)
fig, ax = plt.subplots(figsize=(6, 3.5))
ax.scatter(D, Y, alpha=0.4, color="#005e73")
ax.plot(grid, yhat, color="#FF7E15", linewidth=1.8)
ax.set_xlabel("newspaper"); ax.set_ylabel("sales")
plt.tight_layout(); plt.show()

Kernel regression is one way to get a flexible (non-parametric) estimate of E[Y|D]. R uses np::npregbw, which automatically chooses bandwidth by cross-validation. Python uses statsmodels’ KernelReg with a ‘c’ (continuous) variable type; for production work scikit-learn’s KernelRidge is also a reasonable choice. Both produce the same flexible curve. The marginal effect at any point D = d is just the slope of the curve there — visibly larger near D = 0 and smaller at higher spend, suggesting diminishing returns to newspaper advertising.

Including covariates reduces variance — under randomization, NOT bias correction

Under randomization, \(D \perp X\), so adding pre-treatment covariates does not change the population value of \(\hat\beta_1\) — it cannot fix bias that does not exist.

What it does change is precision: the residual variance shrinks, and so does \(se(\hat\beta_1)\).

Augment the regression with covariates:

\[Y_i \;=\; \underbrace{\hat\beta_0 + \hat\beta_1 D_i + \hat\beta_{X_1} X_{i1} + \ldots + \hat\beta_{X_K} X_{iK}}_{\widehat{\mathbb{E}}[Y_i \mid D_i, X_i]} + \hat\epsilon_i.\]

The covariates explain part of the variation in \(Y\) that was previously residual:

\[R^2 \uparrow, \qquad \frac{\widehat{\text{Var}}(\hat\epsilon_i)}{\widehat{\text{Var}}(Y_i)} \downarrow, \qquad se(\hat\beta_1) \downarrow.\]

Under non-random assignment, covariates serve a different purpose

When \(D\) is not randomized, including covariates is for bias correction (de-confounding), not just variance reduction — Session 06 spends a full lecture on that. Don’t confuse the two roles.

This is the load-bearing legacy point — preserve the framing. Covariates do TWO different jobs depending on the assignment design. Under randomization (this session): covariates do NOT change the expected value of the ATE estimator (treatment is already independent of X), but they DO reduce its variance by absorbing residual variation. This is precision improvement, not bias correction. Under non-random assignment (Session 06): covariates are used to BLOCK backdoor paths and de-confound the estimator. Same regression syntax, fundamentally different statistical roles. Confusing the two is a common conceptual error.

Lin (2013) covariate-adjusted estimator: interact \(D\) with mean-centred \(X\)

Freedman (2008): naively adding covariates additively, \(Y \sim D + X\), can introduce small-sample bias when treatment-effect heterogeneity is present and \(X\) is unbalanced.

Why: \(Y \sim D + X\) assumes one effect \(\beta_1\) for everyone. With true heterogeneity \(\tau(X)\), OLS instead recovers a \(\widehat{\text{Var}}(D \mid X)\)-weighted average of conditional ATEs — and finite-sample imbalance in \(X\) tilts those weights off the marginal density that defines the ATE.

Lin (2013) proposes the canonical fix — interact \(D\) with mean-centred covariates:

\[Y_i \;=\; \beta_0 + \beta_1 D_i + \boldsymbol{\gamma}^\top (X_i - \bar X) + \boldsymbol{\delta}^\top D_i \cdot (X_i - \bar X) + \epsilon_i.\]

The coefficient \(\hat\beta_1\) is the Lin estimator of the ATE.

We’ll use the Lin estimator in the coding lab and compare its SE to the unadjusted mean comparison and to the additive lm(Y ~ D + X).

The Lin estimator is the modern best-practice default for covariate adjustment in randomized experiments. The story: Freedman (2008) showed that mechanically adding covariates to an OLS of Y on D can actually HURT under treatment-effect heterogeneity. Lin (2013) showed that interacting D with mean-centred X repairs this — the resulting estimator is asymptotically at least as efficient as the unadjusted mean comparison and is robust to misspecification of the outcome model. The mean-centring is critical: it ensures that the coefficient on D recovers the ATE, not the conditional effect at X = 0. In R, estimatr::lm_lin does this in one call. In Python, the formula API supports the centred-interaction syntax directly. Recommended workflow: always default to Lin when you have pre-treatment covariates available; never use lm(Y ~ D + X) without centering and interacting unless you really mean to model the conditional effect.

Pre-treatment vs. post-treatment covariates: only the former are safe

Pre-treatment \(X\) — safe

\(X\) is exogenous to \(D\). Adding \(X\) to the regression reduces variance without changing the expectation of \(\hat\beta_1\).

Post-treatment \(X\) — unsafe

Conditioning on \(X\) here opens a collider path through \(U\) — bridge to S03’s collider lesson.

Two harms from conditioning on post-treatment \(X\): (1) mediator-blocking — closes part of \(D \to X \to Y\); (2) collider-opening — opens \(D \to X \leftarrow U \to Y\), breaking randomization.

Two side-by-side DAGs make the pre/post distinction visceral. Left: X is a parent of Y but unrelated to D — a typical pre-treatment covariate (think: prior-period spend, age, demographic). Adding X to the regression reduces residual variance without affecting the expectation of beta_1-hat. Right: X is downstream of D and shares a hidden parent U with Y — a post-treatment “covariate” (think: a mid-experiment metric like in-app engagement). Conditioning on X here does TWO bad things at once: (1) it blocks part of the causal pathway D -> X -> Y; (2) it opens a collider path D -> X <- U -> Y, contaminating the estimator with confounding from the unobserved U. The collider story is the punchline from Session 02 returning in a randomized-experiment context. Bottom line: never condition on anything measured after the treatment unless you have a very specific identification argument for doing so.

Including-covariates example: coffee-leaflet RCT

522 university students, randomly assigned to receive (or not) an environmental-impact leaflet about coffee. Outcome: awareness about waste production (awarewaste); covariates: mother’s education (mumedu) and student sex.

• R
• Python

library(readr); library(here)
library(sandwich); library(modelsummary)
cl <- read_csv(here("materials/data/coffeeleaflet.csv"),
               show_col_types = FALSE)
ols <- lm(awarewaste ~ treatment + mumedu + sex, data = cl)
modelsummary(ols,
  vcov = sandwich::vcovHC,
  estimate = "{estimate} (se = {std.error}, t = {statistic}){stars}",
  statistic = "p = {p.value}, CI = [{conf.low}, {conf.high}]",
  stars = TRUE,
  gof_map = c("r.squared", "nobs"),
  output = "markdown")

import pandas as pd, statsmodels.formula.api as smf
cl = pd.read_csv("../../materials/data/coffeeleaflet.csv")
ols = smf.ols(
    "awarewaste ~ treatment + mumedu + sex",
    data=cl).fit(cov_type="HC1")
print(ols.summary().tables[1])

(Intercept)	1.187 (se = 0.249, t = 4.770)***
	p = <0.001, CI = [0.698, 1.676]
treatment	0.332 (se = 0.096, t = 3.449)***
	p = <0.001, CI = [0.143, 0.521]
mumedu	0.272 (se = 0.090, t = 3.007)**
	p = 0.003, CI = [0.094, 0.450]
sex	0.137 (se = 0.100, t = 1.370)
	p = 0.171, CI = [-0.060, 0.333]
R2	0.046
Num.Obs.	466

Significance: + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001.

The coefficient on treatment is the ATE of the leaflet on waste-awareness, with mumedu and sex absorbing residual variation in \(Y\) — variance reduction, not bias correction (since assignment was randomized).

The coffeeleaflet study is a small-but-clean RCT. Treatment is a randomly distributed leaflet about environmental impacts of coffee production; outcome is a self-reported waste-awareness scale; the two covariates (mother’s education and sex) are pre-treatment demographics. The legacy code used attach() — bad practice that masks the data frame in the global environment. We’ve rewritten using the data argument to lm() and the formula API in Python. Reading the regression: the ATE of treatment is the coefficient on treatment, with HC-robust SE; the mumedu and sex coefficients are absorbed-variance covariates and shouldn’t be interpreted as ATEs of those variables (they’re not randomized).

References

Fernandes, Ana, Martin Huber, and Giannina Vaccaro. 2020. “Gender Differences in Wage Expectations.” https://arxiv.org/abs/2003.11496.

Huber, Martin. 2023. Causal Analysis: Impact Evaluation and Causal Machine Learning with Applications in R. Cambridge, MA: MIT Press.

Wooldridge, Jeffrey M. 2010. Econometric Analysis of Cross Section and Panel Data. 2nd ed. Cambridge, MA: MIT Press.

Thank you for your attention!
	www.startupengineer.io/team/ihl/ christoph-ihl christophihl