5 Randomized Experiments, Regression Adjustment, and Calibration
5.1 Randomized Experiments
For each unit we observe only one of its two potential outcomes, \(Y_i = Y_i(T_i)\), so the unit-level effect \(\tau_i = Y_i(1) - Y_i(0)\) is never observed. This is the fundamental problem of causal inference (Holland 1986). Holland distinguishes two ways around it. The scientific solution measures both outcomes on units that are assumed to be interchangeable, for example identical twins or the same piece of equipment before and after a repair; it rests on a homogeneity assumption that the data cannot check. The statistical solution assigns treatment at random and studies averages. This section develops the statistical solution.
Two kinds of studies differ in who decides the treatment. In a randomized experiment, a known random mechanism decides which units receive \(T_i = 1\); a clinical trial is the standard example. In an observational study, other factors decide \(T_i\): a physician’s judgment, a participant’s choice, an administrative rule. The mechanism is unknown. The rest of this section uses one fact only: in an experiment the assignment mechanism is known, and that knowledge supplies both the identification of the causal effect and the probability distribution used for inference.
Notation. Throughout the design-based parts of this chapter we use the notation of survey sampling. The \(N\) units in the experiment form a finite population \[\mathcal F_N = \bigl\{\bigl(Y_i(1),\, Y_i(0)\bigr) : i = 1, \ldots, N\bigr\},\] treated as fixed numbers. \(N_1\) units are treated and \(N_0 = N - N_1\) are controls, so a sum over \(i = 1, \ldots, N\) is a sum over the population, and a sum over \(\{i : T_i = 1\}\) is a sum over the treated arm. The target is the finite-population average treatment effect \[\bar\tau_N \;=\; \frac1N\sum_{i=1}^N \bigl\{Y_i(1) - Y_i(0)\bigr\} \;=\; \bar Y(1) - \bar Y(0), \qquad \bar Y(t) = \frac1N\sum_{i=1}^N Y_i(t).\]
5.1.1 A Four-Unit Example
A small example shows every idea of this section with numbers that can be checked by hand.
Example 5.1 (The Science Table) Four units have the potential outcomes in Table 5.1. The table is fixed: it describes the four units before any assignment is made, and nothing in it is random. The finite-population ATE is \(\bar\tau_4 = (3 + 3 + 2 + 2)/4 = 2.5\).
We will never see this table. If unit \(i\) is treated we see \(Y_i(1)\), otherwise \(Y_i(0)\); no unit shows both columns, so no \(\tau_i\) is ever observed. The question is how to learn \(\bar\tau_4 = 2.5\) without seeing a single unit-level effect.
| Unit \(i\) | \(Y_i(1)\) | \(Y_i(0)\) | \(\tau_i = Y_i(1) - Y_i(0)\) |
|---|---|---|---|
| 1 | 4 | 1 | 3 |
| 2 | 6 | 3 | 3 |
| 3 | 7 | 5 | 2 |
| 4 | 8 | 6 | 2 |
| Mean | 6.25 | 3.75 | 2.5 |
Suppose we treat \(N_1 = 2\) units and leave \(N_0 = 2\) as controls, and estimate \(\bar\tau_4\) by the difference in means (DIM), \[\hat\tau_{\mathrm{DIM}} \;=\; \bar Y_1 - \bar Y_0 \;=\; \frac1{N_1}\sum_{i:\,T_i=1} Y_i \;-\; \frac1{N_0}\sum_{i:\,T_i=0} Y_i .\] There are \(\binom42 = 6\) ways to choose the treated pair. Table 5.2 lists them. For example, if units 1 and 4 are treated, we observe \(Y_1(1) = 4\), \(Y_4(1) = 8\), \(Y_2(0) = 3\), and \(Y_3(0) = 5\), so \(\hat\tau_{\mathrm{DIM}} = (4 + 8)/2 - (3 + 5)/2 = 2\).
| Case | Treated | Control | \(\bar Y_1 - \bar Y_0\) | \(\hat\tau_{\mathrm{DIM}}\) |
|---|---|---|---|---|
| 1 | 1, 2 | 3, 4 | \(5.0 - 5.5\) | \(-0.5\) |
| 2 | 1, 3 | 2, 4 | \(5.5 - 4.5\) | \(1.0\) |
| 3 | 1, 4 | 2, 3 | \(6.0 - 4.0\) | \(2.0\) |
| 4 | 2, 3 | 1, 4 | \(6.5 - 3.5\) | \(3.0\) |
| 5 | 2, 4 | 1, 3 | \(7.0 - 3.0\) | \(4.0\) |
| 6 | 3, 4 | 1, 2 | \(7.5 - 2.0\) | \(5.5\) |
Only one of the six cases will actually happen, and the estimate we report depends on which one. The spread, from \(-0.5\) to \(5.5\) around the target \(2.5\), is the uncertainty that the design creates.
5.1.2 The Completely Randomized Design and the Design-Based View
Definition 5.1 (Completely Randomized Design) A completely randomized design (CRD) chooses \(N_1\) of the \(N\) units for treatment, giving each of the \(\binom{N}{N_1}\) possible assignment vectors \(\mathbf T = (T_1, \ldots, T_N)\) with \(\sum_i T_i = N_1\) the same probability \(1/\binom{N}{N_1}\). Equivalently, the treated arm is a simple random sample without replacement of size \(N_1\) from the \(N\) units, and the control arm is its complement.
Given \(\mathcal F_N\), the only random quantity is \(\mathbf T\). The distribution of an estimator over the \(\binom{N}{N_1}\) equally likely assignments is its randomization distribution.
Example 5.2 (The Randomization Distribution of the Four-Unit Example) Under the CRD with \(N_1 = 2\), each case in Table 5.2 has probability \(1/6\), so \(\hat\tau_{\mathrm{DIM}}\) takes the values \(-0.5, 1, 2, 3, 4, 5.5\) with probability \(1/6\) each. Its mean is \[\E\bigl(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_4\bigr) = \tfrac16(-0.5 + 1 + 2 + 3 + 4 + 5.5) = \tfrac{15}{6} = 2.5 = \bar\tau_4,\] so the estimator is unbiased over the randomization. The squared deviations from \(2.5\) are \(9, 2.25, 0.25, 0.25, 2.25, 9\), so \[\mathrm{Var}\bigl(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_4\bigr) = \tfrac{23}{6} \approx 3.83 .\]
The design-based view. The example made no assumption about how the potential outcomes are generated. The reference distribution is the sampling distribution produced by repeating the given assignment mechanism on the same units. This is the design-based approach: probability statements refer to the physical randomization, not to a model for the outcomes. Readers who know survey sampling will recognize the logic, in which a known sampling design, not a model, supplies the probability.
Randomization serves two purposes (Fisher 1935). First, it makes the treatment and control groups comparable on average, which is the basis of identification (Section 5.1.3). Second, it gives a “reasoned basis” for inference: the probability distribution used for estimation and testing is created by the design, not assumed (Section 5.1.4 and Section 5.1.7).
5.1.3 Identification under Randomization
The four-unit example showed by enumeration that the difference in means averages to the target. The same fact has a population-level form.
Lemma 5.1 (Identification under Complete Randomization) Assume consistency, \(Y_i = T_iY_i(1) + (1 - T_i)Y_i(0)\), and randomization, \[\bigl\{Y(1),\, Y(0)\bigr\} \indep T .\] Then, for \(t \in \{0, 1\}\), \[\E(Y \mid T = t) \;=\; \E\{Y(t)\}, \tag{5.1}\] so the difference in means targets \(\tau_{\mathrm{ATE}} = \E\{Y(1) - Y(0)\}\).
Proof. By consistency, \(\E(Y \mid T = t) = \E\{Y(t) \mid T = t\}\). By randomization, \(Y(t)\) is independent of \(T\), so \(\E\{Y(t) \mid T = t\} = \E\{Y(t)\}\). \(\square\)
Here \(\E\) refers either to a randomly selected unit of the trial population (the empirical distribution of \(\mathcal F_N\)) or to the two-phase superpopulation model of Section 5.1.6. Randomization makes the potential outcomes independent of assignment unconditionally; no covariate adjustment is needed. This is the unconditional counterpart of ignorability (Definition 4.5). In stratified, blocked, or covariate-adaptive designs, where the assignment probability depends on design variables, the natural statement is conditional on those variables. Throughout this chapter, “randomized experiment” means the CRD unless stated otherwise.
5.1.4 Neyman’s Theorem: Exact Unbiasedness and Design Variance
The four-unit example computed the mean and variance of the randomization distribution by listing all six assignments. Neyman’s theorem gives both in closed form for any \(N\) and \(N_1\) (Neyman 1923). Define the finite-population variances and covariance (divisor \(N - 1\)) \[\begin{aligned} S_t^2 &= \frac1{N-1}\sum_{i=1}^N\bigl\{Y_i(t) - \bar Y(t)\bigr\}^2, \qquad t \in \{0, 1\},\\ S_{01} &= \frac1{N-1}\sum_{i=1}^N\bigl\{Y_i(0) - \bar Y(0)\bigr\}\bigl\{Y_i(1) - \bar Y(1)\bigr\}, \end{aligned}\] and the finite-population variance of the unit-level effects, \(S_\tau^2 = (N-1)^{-1}\sum_{i=1}^N(\tau_i - \bar\tau_N)^2 = S_1^2 + S_0^2 - 2S_{01}\).
Theorem 5.1 (Neyman’s Theorem for the Difference in Means) Under the CRD, with expectations and variances taken over the randomization distribution of \(\mathbf T\) given \(\mathcal F_N\):
- Unbiasedness. \(\E\bigl(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_N\bigr) = \bar\tau_N\).
- Design variance. \[\mathrm{Var}\bigl(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_N\bigr) \;=\; \frac{S_1^2}{N_1} + \frac{S_0^2}{N_0} - \frac{S_\tau^2}{N} \;=\; \Bigl(\frac1{N_1} - \frac1N\Bigr)S_1^2 + \Bigl(\frac1{N_0} - \frac1N\Bigr)S_0^2 + \frac2N\,S_{01}. \tag{5.3}\]
- Conservative variance estimator. Let \(N_1, N_0 \ge 2\) and let \(\hat S_t^2\) be the sample variance of the observed outcomes in arm \(t\). Then \[\hat V \;=\; \frac{\hat S_1^2}{N_1} + \frac{\hat S_0^2}{N_0} \tag{5.4}\] satisfies \(\E(\hat V \mid \mathcal F_N) - \mathrm{Var}(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_N) = S_\tau^2/N \ge 0\). It is unbiased if and only if the unit-level effects are constant.
Example 5.3 (Neyman’s Theorem in the Four-Unit Example) From Table 5.1, \(S_1^2 = 35/12\), \(S_0^2 = 59/12\), and \(S_\tau^2 = 1/3\). With \(N_1 = N_0 = 2\) and \(N = 4\), the design variance is \[\frac{35/12}{2} + \frac{59/12}{2} - \frac{1/3}{4} = \frac{35}{24} + \frac{59}{24} - \frac{1}{12} = \frac{23}{6},\] the variance of the six values found in Example 5.2. Averaging \(\hat V\) over the six assignments gives \(\E(\hat V \mid \mathcal F_4) = 47/12 \approx 3.92\), larger than \(23/6 \approx 3.83\) by \(1/12 = S_\tau^2/N\). The individual values of \(\hat V\), however, range from \(5/4\) (cases 1 and 6) to \(13/2\) (case 4). Conservative here means that \(\hat V\) overestimates the design variance in expectation; it need not exceed it for every realized assignment.
5.1.5 The Difference in Means as a Regression Slope
The difference in means is also an ordinary least squares (OLS) coefficient. Minimizing \(\sum_{i=1}^N(Y_i - \alpha - \tau T_i)^2\) gives \(\hat\alpha = \bar Y_0\) and \(\hat\tau = \bar Y_1 - \bar Y_0 = \hat\tau_{\mathrm{DIM}}\). The regression view connects this section to the covariate-adjusted estimators of Section 5.3 and Section 5.4, but its standard errors need care.
A heteroskedastic regression. In a superpopulation reading, let \(Y_i(t) = \mu_t + e_i(t)\) with \(\E\{e_i(t)\} = 0\) and \(\mathrm{Var}\{e_i(t)\} = \sigma_t^2\); the two errors of a unit may be correlated. Then \[Y_i = \mu_0 + \tau T_i + e_i, \qquad e_i = T_ie_i(1) + (1 - T_i)e_i(0), \qquad \tau = \mu_1 - \mu_0,\] and randomization gives \(\E(e_i \mid T_i) = 0\). But \(\mathrm{Var}(e_i \mid T_i = t) = \sigma_t^2\): the error variance differs between arms. The classical OLS standard error pools the two variances. With \(N_1, N_0 \ge 2\), \[\hat V_{\mathrm{pooled}} - \hat V = \frac{(N - 1)(N_1 - N_0)(\hat S_1^2 - \hat S_0^2)}{(N - 2)N_1N_0},\] so the pooled standard error agrees with Neyman’s estimator for every data set under balanced allocation. With unequal arm sizes the two generally differ, although they also agree when the arm-specific sample variances are equal.
Proposition 5.1 (The Leverage-Corrected Sandwich Variance Equals Neyman’s Estimator) In the regression of \(Y_i\) on \((1, T_i)\), let \(\hat e_i\) be the OLS residuals and \(h_{ii}\) the leverages. The heteroskedasticity-robust (sandwich) variance estimator of \(\hat\tau\) in which each squared residual is divided by \(1 - h_{ii}\) (the “HC2” estimator) equals \(\hat V = \hat S_1^2/N_1 + \hat S_0^2/N_0\), for any \(N_1, N_0 \ge 2\).
Proof. \(\hat\tau = \sum_ia_iY_i\) with \(a_i = 1/N_1\) for treated units and \(a_i = -1/N_0\) for controls, and the robust variance with this residual correction is \(\sum_ia_i^2\hat e_i^2/(1 - h_{ii})\). The fitted value of unit \(i\) is its arm mean, so \(\hat e_i = Y_i - \bar Y_t\) for \(i\) in arm \(t\), and \(h_{ii} = 1/N_t\). For the treated arm, \[\sum_{i:T_i=1}\frac{\hat e_i^2}{N_1^2(1 - 1/N_1)} = \frac{1}{N_1}\cdot\frac{1}{N_1 - 1}\sum_{i:T_i=1}(Y_i - \bar Y_1)^2 = \frac{\hat S_1^2}{N_1},\] and likewise for the control arm. \(\square\)
Section 5.4 adds covariates to this regression. Proposition 5.1 is the simplest instance of a general fact used there: a sandwich variance from the regression estimates the design variance conservatively.
5.1.6 Finite-Population and Superpopulation Targets
Neyman’s theorem conditions on the \(N\) enrolled units and targets \(\bar\tau_N\). Much of the literature instead targets the superpopulation ATE \(\tau_{\mathrm{ATE}} = \E\{Y(1) - Y(0)\}\). Table 5.3 contrasts the two frameworks.
| Finite-population randomization | Superpopulation framework | |
|---|---|---|
| Objects | \(\mathcal F_N\) is fixed | Units are sampled from a population |
| Target | \(\bar\tau_N = \bar Y(1) - \bar Y(0)\) | \(\tau_{\mathrm{ATE}} = \E\{Y(1) - Y(0)\}\) |
| Randomness | Assignment only | Sampling and assignment |
| Scope | The units actually randomized | A wider target population |
| Population link | Not needed | Random sampling or transportability |
Randomization alone justifies design-based inference for \(\bar\tau_N\), the effect among the units randomized (internal validity). Extending the conclusion to a wider population requires random sampling or a transportability argument (external validity). The next remark shows how the two variances relate when the units are themselves a random sample.
The parallel with survey sampling. A completely randomized experiment is simple random sampling without replacement from the finite population of potential outcomes: the treated arm is a sample of the \(Y_i(1)\), and the control arm, its complement, is a sample of the \(Y_i(0)\). This is why Neyman’s variance formula has the form of a survey-sampling variance. Table 5.4 sets out the parallel; it returns in Section 5.4, where the regression estimator of survey sampling becomes a covariate-adjusted estimator of \(\bar\tau_N\).
| Survey sampling | Causal inference | |
|---|---|---|
| Parameter | Descriptive: a finite-population total or mean | Causal: \(\bar\tau_N\), the ATE, the ATT |
| What is missing | \(y_i\) for unsampled units | The unassigned potential outcome of every unit |
| Design-based gold standard | Probability sampling with known inclusion probabilities | Randomized experiment with known assignment probabilities |
| Partial failure of the design | Nonresponse | Noncompliance |
| No randomization | Nonprobability sample (model-based) | Observational study (model-based) |
5.1.7 Fisher’s Randomization Test
Neyman’s theorem uses the randomization to estimate an average effect. Fisher used it to test a hypothesis exactly (Fisher 1935). The hypothesis is the sharp null of no effect for any unit, \[H_0:\; Y_i(1) = Y_i(0) \quad\text{for every } i = 1, \ldots, N.\] Under \(H_0\) each unit’s outcome does not depend on its assignment, so the observed outcome vector \(\mathbf Y\) would have been the same under every assignment. The whole science table is then known, and the randomization distribution of any test statistic can be computed exactly. The two enumerations of this section differ in this respect: Table 5.2 uses the given potential-outcome schedule, so the observed outcomes change with the assignment, while the Fisher enumeration uses the schedule implied by the sharp null, so the outcomes stay fixed.
Example 5.4 (Fisher’s Test in the Four-Unit Example) Suppose case 3 of Table 5.2 happened: units 1 and 4 were treated, and \(\mathbf Y = (4, 3, 5, 8)\), so \(\hat\tau_{\mathrm{DIM}} = 6 - 4 = 2\). Under \(H_0\) the vector \(\mathbf Y\) stays the same under every assignment. Recomputing the difference in means for the six assignments, holding \(\mathbf Y\) fixed, gives \[\text{treated }\{1,2\}: -3,\quad \{1,3\}: -1,\quad \{1,4\}: 2,\quad \{2,3\}: -2,\quad \{2,4\}: 1,\quad \{3,4\}: 3 .\] Two of the six values are at least 2, so the one-sided \(p\)-value is \(2/6 = 1/3\); four have absolute value at least 2, so the two-sided \(p\)-value is \(4/6\). With \(N = 4\) the smallest attainable one-sided \(p\)-value is \(1/6\). For the two-sided test based on \(|\hat\tau_{\mathrm{DIM}}|\) it is \(2/6\), because with balanced allocation complementary assignments give the same absolute value. Either way, this experiment can never reject at the 5% level, however large the effect. The true effects in Table 5.1 are 2 or 3, so \(H_0\) is false; the test simply has no power at \(N = 4\).
Proposition 5.2 (Fisher Randomization Test) Let \(Q = Q(\mathbf T, \mathbf Y)\) be a test statistic, for example \(\hat\tau_{\mathrm{DIM}}\), and let \(\mathbf T^{(1)}, \ldots, \mathbf T^{(M)}\) be all assignments the design allows. Assume that each has the same probability \(1/M\), as under the CRD with \(M = \binom{N}{N_1}\). Define \[p_{\mathrm{FRT}} = \frac1M\sum_{k=1}^M \mathbb{I}\bigl\{Q(\mathbf T^{(k)}, \mathbf Y) \ge Q(\mathbf T, \mathbf Y)\bigr\}.\] Under \(H_0\), \(P(p_{\mathrm{FRT}} \le u) \le u\) for every \(u \in (0, 1)\), in finite samples.
Proof. Under \(H_0\), \(\mathbf Y\) is fixed, so the \(M\) values \(q_k = Q(\mathbf T^{(k)}, \mathbf Y)\) are fixed numbers, and the observed statistic is \(q_K\) for an index \(K\) that is uniform on \(\{1, \ldots, M\}\). Write \(p_k = M^{-1}\#\{j : q_j \ge q_k\}\), so \(p_{\mathrm{FRT}} = p_K\). Since \(p_k\) is nonincreasing in \(q_k\), the set \(A = \{k : p_k \le u\}\) has the form \(\{k : q_k \ge c\}\) for some \(c\). If \(A\) is nonempty, let \(k^\ast \in A\) have the smallest \(q_k\) in \(A\); every member of \(A\) has \(q_k \ge q_{k^\ast}\), so \(|A|/M \le p_{k^\ast} \le u\). Hence \(P(p_{\mathrm{FRT}} \le u) = P(K \in A) = |A|/M \le u\). \(\square\)
For a design with unequal assignment probabilities, the exhaustive tail probability must weight each assignment by its actual probability. Three practical points follow. First, \(M\) grows quickly with \(N\). In practice one draws \(B\) assignments \(\mathbf T^{(1)}, \ldots, \mathbf T^{(B)}\) independently from the actual assignment mechanism, independently of the observed assignment. The fraction of sampled statistics at least as large as the observed one estimates \(p_{\mathrm{FRT}}\), but it is not itself a valid test \(p\)-value for finite \(B\); it can even be zero. A valid, possibly conservative, Monte Carlo \(p\)-value is \[p_{\mathrm{MC}} = \frac{1 + \sum_{b=1}^B \mathbb{I}\bigl\{Q(\mathbf T^{(b)}, \mathbf Y) \ge Q(\mathbf T, \mathbf Y)\bigr\}}{B + 1}:\] under the sharp null the observed statistic and the \(B\) simulated ones are exchangeable, so the rank of the observed statistic among all \(B + 1\) values gives the test (Phipson and Smyth 2010). Second, for a two-sided test use a two-sided statistic, such as \(|\hat\tau_{\mathrm{DIM}}|\). Third, the reference set must match the actual design: a permutation scheme that does not follow the assignment mechanism is not a Fisher randomization test. Fisher’s exact test for a \(2 \times 2\) table is the special case with a binary outcome.
Fisher and Neyman ask different questions. Both frameworks rest on the known assignment mechanism, but they answer different questions (Table 5.5). The two are complementary rather than contradictory.
| Fisher | Neyman | |
|---|---|---|
| Target | The sharp null \(Y_i(1) = Y_i(0)\) for every unit | The average effect \(\bar\tau_N\) |
| Main task | An exact test under the known assignment mechanism | Point estimation and uncertainty quantification |
| Why it works | The sharp null fills in all missing outcomes | The repeated-randomization behavior of \(\hat\tau_{\mathrm{DIM}}\) is analyzable |
| Output | A randomization \(p\)-value | An estimate, a standard error, a confidence interval |
5.2 Ignorability: A Short Review
Section 5.1 showed that randomization makes the potential outcomes independent of the assignment, so the difference in means targets the ATE without any model. In an observational study the treatment is not assigned by the researcher; it is selected by patients, firms, governments, or individuals, often on the basis of characteristics that also affect the outcome. Chapter 4 developed the assumption that replaces randomization in this setting. This section collects the parts of that development that the rest of the chapter uses.
Statistical adjustment can implement ignorability once it is assumed, but it cannot create it. No covariate adjustment, however flexible, closes a back-door path that no observed pretreatment variable intercepts.
5.2.1 Terminology and What Ignorability Requires
Strong ignorability (Definition 4.5; Rosenbaum and Rubin (1983)) combines joint unconfoundedness and overlap: \[\bigl\{Y(0), Y(1)\bigr\} \indep T \mid X \qquad\text{and}\qquad 0 < P(T = 1 \mid X) < 1 \ \text{ almost surely}.\] The weak (pointwise) form, \(Y(t) \indep T \mid X\) for each \(t\) separately, is implied by the joint form and is all that the identification of \(\E\{Y(t)\}\) uses. Neither form identifies the joint law of \(\{Y(0), Y(1)\}\), for example \(P\{Y(1) > Y(0)\}\); that requires cross-world restrictions (Chapter 4). Other names for the same idea are unconfoundedness, selection on observables, no unmeasured confounders, and conditional exchangeability. Authors differ on weak versus joint forms and on whether positivity and consistency are included.
Ignorability is a substantive assumption, not a statistical one. A sufficient graphical condition is an observed pretreatment set \(X\) that blocks every back-door path from \(T\) to \(Y\). This does not require observing every common cause: in the graph \(T \leftarrow U \to W \to Y\) of Chapter 3, the unobserved \(U\) is harmless because the observed \(W\) blocks the only back-door path. The phrase “no unmeasured confounders” should be read in this path-blocking sense.
5.2.2 From Ignorability to Identification
Write \(\mu(t, x) = \E(Y \mid T = t, X = x)\) for the outcome regression.
Theorem 5.2 (Identification of the ATE and the ATT) Assume consistency, weak ignorability \(Y(t) \indep T \mid X\) for \(t \in \{0, 1\}\), overlap, and \(\E|Y(t)| < \infty\). Then \[\tau_{\mathrm{ATE}} = \E\{Y(1) - Y(0)\} = \int \bigl\{\mu(1, x) - \mu(0, x)\bigr\}\, p(x)\, dx, \tag{5.7}\] \[\tau_{\mathrm{ATT}} = \E\{Y(1) - Y(0) \mid T = 1\} = \int \bigl\{\mu(1, x) - \mu(0, x)\bigr\}\, p(x \mid T = 1)\, dx. \tag{5.8}\]
Proof. For the ATE, Proposition 4.2 gives \(\E\{Y(t)\} = \int \mu(t, x)\,p(x)\,dx\) for each \(t\); subtract. For the ATT, iterate expectations within the treated: \(\tau_{\mathrm{ATT}} = \E[\E\{Y(1) - Y(0) \mid T = 1, X\} \mid T = 1]\). Consistency gives \(\E\{Y(1) \mid T = 1, X\} = \mu(1, X)\). For the counterfactual mean, \(\E\{Y(0) \mid T = 1, X\} = \E\{Y(0) \mid X\} = \E\{Y(0) \mid T = 0, X\} = \mu(0, X)\), by ignorability (twice) and consistency; control-arm overlap on the support of \(X \mid T = 1\) makes the last conditional mean well defined. Integrate against \(p(x \mid T = 1)\). \(\square\)
5.2.3 Which Variables to Condition On: The Pretreatment Requirement
The back-door formula conditions on \(X\). Before estimating anything, there is a prior question that is purely graphical: which variables belong in \(X\)? The answer is not “all available variables.”
The practical rule. Classify every candidate covariate as pre-treatment or post-treatment using the causal graph, before running any regression. For the back-door adjustment strategy used here, choose a set of pretreatment covariates that satisfies the criterion; validity is a property of the set jointly, not of each variable separately.
5.3 Estimating the ATE by Outcome Regression
Identification is now settled: under the conditions of Theorem 5.2, \[\tau_{\mathrm{ATE}} = \E\bigl\{\mu(1, X) - \mu(0, X)\bigr\} = \int \bigl\{\mu(1, x) - \mu(0, x)\bigr\}\, p(x)\, dx .\] The right side contains only observable quantities, but it is a population quantity. This section turns it into an estimator and studies its sampling properties. Throughout, the data are \(N\) i.i.d. observations \((Y_i, T_i, X_i)\) from a superpopulation, \(X\) is a set of pretreatment covariates satisfying the back-door criterion, and we write \[P_N f = \frac1N\sum_{i=1}^N f(Y_i, T_i, X_i), \qquad P f = \E\{f(Y, T, X)\}.\] The letter \(N\) is the number of observed units here as in Section 5.1, but the probability model is different: there the units were fixed and the assignment random; here the units are a random sample.
5.3.1 The Plug-In Principle
The identified functional has three unknowns, and each has a natural sample replacement (Table 5.6).
| Unknown | Plug-in replacement | What it costs |
|---|---|---|
| \(p(x)\) | The empirical distribution of \(X_1, \ldots, X_N\): \(\int g(x)\,p(x)\,dx \approx P_N g\) | No model for the distribution of \(X\) |
| \(\mu(1, \cdot)\) | An estimate \(\hat\mu(1, \cdot)\) from the treated units | A model or a smoothing method |
| \(\mu(0, \cdot)\) | An estimate \(\hat\mu(0, \cdot)\) from the control units | A model or a smoothing method |
The plug-in estimator is the difference of two regression estimators, \[\hat\tau_{\mathrm{reg}} = \hat\mu_{1,\mathrm{reg}} - \hat\mu_{0,\mathrm{reg}}, \qquad \hat\mu_{t,\mathrm{reg}} = \frac1N\sum_{i=1}^N \hat\mu(t, X_i), \tag{5.9}\] where \(\hat\mu_{t,\mathrm{reg}}\) is the regression estimator of \(\E\{Y(t)\}\). The averaging step needs no model for the covariate distribution; all assumptions of this section concern \(\hat\mu\).
5.3.2 The Discrete Case: Simpson’s Paradox and Standardization
When \(X\) is categorical, \(\mu(t, x)\) can be estimated by the cell mean \(\bar Y_{t,x}\) of \(Y\) among units with \(T = t\) and \(X = x\), and Equation 5.9 becomes the standardization estimator \[\hat\tau_{\mathrm{STD}} = \sum_{x \in \mathcal X}\bigl(\bar Y_{1,x} - \bar Y_{0,x}\bigr)\,\hat p(x), \qquad \hat p(x) = \frac{N_x}{N}. \tag{5.10}\] This is \(\hat\tau_{\mathrm{reg}}\) with a saturated outcome model, one free parameter per cell, so no functional-form assumption is made. It is also called the g-formula (Robins 1986) and direct standardization.
Example 5.5 (Simpson’s Paradox as a Standardization) Return to the vaccination data of Chapter 1 (Section 1.4.3): \(N = 1{,}000\) people, with \(X = 1\) for the elderly. The vaccinated group contains 360 elderly and 60 young people; the unvaccinated group, 40 elderly and 540 young. The risks of illness are \(\mu(1, \text{elderly}) = 30\%\), \(\mu(0, \text{elderly}) = 50\%\), \(\mu(1, \text{young}) = 5\%\), and \(\mu(0, \text{young}) = 10\%\), so the vaccine lowers the risk in both strata. Three contrasts weight the strata differently.
- Pooled. Each arm uses its own covariate mix, \(P(X = 1 \mid T = 1) = 360/420 = 6/7\) and \(P(X = 1 \mid T = 0) = 40/580 \approx 0.069\): \[100\Bigl(\frac{111}{420} - \frac{74}{580}\Bigr) \approx 26.43\% - 12.76\% \approx +13.7 \text{ percentage points}.\]
- ATE (standardized). Both arms use \(P(X = 1) = 0.40\): \(0.4(-20) + 0.6(-5) = -11.0\) percentage points.
- ATT. Both arms use \(P(X \mid T = 1)\): \(\tfrac67(-20) + \tfrac17(-5) \approx -17.9\) percentage points.
The pooled comparison reverses the sign because the vaccinated are mostly elderly, who have a higher baseline risk. The ATT is more negative than the ATE because the treated are concentrated in the stratum with the larger effect. The reversal itself is a fact about conditional proportions; which comparison is causal is a fact about the graph. The standardized answers are causal only under the maintained assumptions: age is a sufficient adjustment set, with positivity and consistency.
5.3.3 The Outcome Model and the Regression Estimator
When \(X\) is continuous or has many categories, cell means are not available, and \(\mu\) is modeled.
Two linear choices of \(m\) are common.
- Model (a), additive: \(m(t, x; \beta) = \alpha + \tau t + \gamma^\top x\). Then \(m(1, x) - m(0, x) = \tau\) for every \(x\), and \(\hat\tau_{\mathrm{reg}}\) is the OLS coefficient on \(T\). The model imposes the same conditional average effect at every covariate value, so it forces \(\tau_{\mathrm{ATE}} = \tau_{\mathrm{ATT}}\); it does not require identical unit-level effects.
- Model (b), one model per arm: \(m(t, x; \beta) = \alpha_t + \gamma_t^\top x\), fitted separately on each arm. Then \(\hat\tau_{\mathrm{reg}} = (\hat\alpha_1 - \hat\alpha_0) + (\hat\gamma_1 - \hat\gamma_0)^\top \bar X\), with \(\bar X = N^{-1}\sum_iX_i\).
Proposition 5.3 (Model (b): Each Regression Estimator Is an Intercept after Centering) Under model (b), write \(\hat\mu(t, x) = \hat a_t + \hat\gamma_t^\top(x - \bar X)\) with \(\hat a_t = \hat\alpha_t + \hat\gamma_t^\top\bar X\). Then
- \(\hat\mu_{t,\mathrm{reg}} = \hat a_t\), the intercept of the arm-\(t\) regression of \(Y\) on \((1, X - \bar X)\), so \(\hat\tau_{\mathrm{reg}} = \hat a_1 - \hat a_0\);
- in the single regression of \(Y\) on \(\bigl(1,\, T,\, X - \bar X,\, T(X - \bar X)\bigr)\) fitted to all \(N\) units, the coefficient on \(T\) equals \(\hat\tau_{\mathrm{reg}}\).
Proof. (i) Averaging over all \(N\) units, \(N^{-1}\sum_i\hat\mu(t, X_i) = \hat a_t + \hat\gamma_t^\top(\bar X - \bar X) = \hat a_t\). (ii) The single regression has a separate intercept and slope for each arm, so its fitted values on arm \(t\) are those of the arm-\(t\) regression on \((1, X - \bar X)\). Its intercept is \(\hat a_0\), and the coefficient on \(T\) is \(\hat a_1 - \hat a_0\). \(\square\)
Centering at the arm-specific mean \(\bar X_t\) instead of the full-sample mean \(\bar X\) would make the intercept equal \(\bar Y_t\), and no adjustment would take place. Under model (b) the ATT estimator is \(\hat a_1 - \hat a_0 + (\hat\gamma_1 - \hat\gamma_0)^\top(\bar X_1 - \bar X)\), with \(\bar X_1\) the treated-arm covariate mean.
5.3.4 A Worked Example
Example 5.6 (Earnings and a Binary Covariate) Let \(X \in \{0, 1\}\) be a binary covariate and \(Y\) earnings in thousands of dollars. The cell means within each arm are:
| \(X\) | \(N\) (treated) | \(\bar Y\) (treated) | \(N\) (control) | \(\bar Y\) (control) |
|---|---|---|---|---|
| 0 | 30 | 42 | 70 | 35 |
| 1 | 70 | 58 | 30 | 50 |
| All | 100 | 53.2 | 100 | 39.5 |
The unadjusted difference is \(53.2 - 39.5 = 13.7\). But \(X = 1\) units, which have higher baseline earnings, make up 70% of the treated and 30% of the controls, so \(X\) confounds the comparison. With the cell means as \(\hat\mu\), the within-stratum contrasts are \(42 - 35 = 7\) for \(X = 0\) and \(58 - 50 = 8\) for \(X = 1\). In the full sample, half the units have \(X = 1\), so \[\hat\tau_{\mathrm{reg}} = 0.5(7) + 0.5(8) = 7.5, \qquad \hat\tau_{\mathrm{ATT}} = 0.3(7) + 0.7(8) = 7.7 .\] The adjusted estimates are roughly half the unadjusted difference. The ATT is slightly larger because the treated are concentrated in the stratum with the larger contrast.
5.3.5 The Goal, the Obstacle, and the Nuisance Parameter
The goal is statistical inference with \(\hat\tau_{\mathrm{reg}}\): consistency for \(\tau_{\mathrm{ATE}}\), an asymptotic distribution, and a variance estimate, so that a confidence interval can be reported. The obstacle is that \(\hat\tau_{\mathrm{reg}}\) is not a sample mean. It averages fitted values, and the data enter through \(\hat\beta\) inside \(m\), so the law of large numbers and the central limit theorem do not apply to it directly.
The source of the difficulty is a second parameter. The parameter of interest is \(\tau_{\mathrm{ATE}} = \E\{\mu(1, X) - \mu(0, X)\}\). The nuisance parameter is \(\beta\), the coefficient of the outcome model, and more generally the regression function \(\mu\) itself. It is not wanted for its own sake, but \(\tau_{\mathrm{ATE}}\) depends on it, so it must be estimated first, and its estimation error propagates into \(\hat\tau_{\mathrm{reg}}\). The same structure returns in later chapters, with the propensity score (Chapter 6) and with two nuisance functions (Chapters 11 and 12).
For a linear outcome model the propagation can be written exactly. Let \[\mu(t, x) = z(t, x)^\top\beta_0, \qquad z_d(x) = z(1, x) - z(0, x), \qquad Z_i = z(T_i, X_i),\] with \(z\) known, and let \(\hat\beta = (P_NZZ^\top)^{-1}P_NZY\) be the OLS estimator. For model (b), \(z(t, x) = (t,\, 1 - t,\, t\,x^\top,\, (1 - t)\,x^\top)^\top\), \(\beta = (\alpha_1, \alpha_0, \gamma_1^\top, \gamma_0^\top)^\top\), and \(z_d(x) = (1, -1, x^\top, -x^\top)^\top\). For model (a), \(z(t, x) = (1, t, x^\top)^\top\) and \(z_d = (0, 1, 0^\top)^\top\) is constant. Because \(m\) is linear in \(\beta\), the average of the predictions is the prediction at the average: \[\hat\tau_{\mathrm{reg}} = P_N\bigl(z_d^\top\hat\beta\bigr) = (P_Nz_d)^\top\hat\beta, \qquad \tau_{\mathrm{ATE}} = (Pz_d)^\top\beta_0 , \tag{5.11}\] where the second equality uses Theorem 5.2 and the correct specification of the model. For model (b), \(P_Nz_d = (1, -1, \bar X^\top, -\bar X^\top)^\top\). The nuisance enters \(\hat\tau_{\mathrm{reg}}\) through one vector, \(\hat\beta\).
5.3.6 Consistency: Estimation Error and Sampling Error
Subtracting the two displays in Equation 5.11 gives an exact decomposition: \[\hat\tau_{\mathrm{reg}} - \tau_{\mathrm{ATE}} = \underbrace{(P_Nz_d)^\top(\hat\beta - \beta_0)}_{\text{estimation error}} + \underbrace{(P_N - P)\bigl(z_d^\top\beta_0\bigr)}_{\text{sampling error}} . \tag{5.12}\] The estimation error comes from estimating the nuisance parameter; the sampling error comes from averaging over the sampled covariates instead of their population distribution.
Proposition 5.4 (Consistency of the Regression Estimator) Assume the conditions of Theorem 5.2, i.i.d. data, a correctly specified linear model \(\E(Y \mid T, X) = Z^\top\beta_0\), finite second moments of \(Y\) and \(Z\), \(\E\|z_d(X)\| < \infty\), and a nonsingular \(A = \E(ZZ^\top)\). Then \(\hat\tau_{\mathrm{reg}} \pto \tau_{\mathrm{ATE}}\).
Proof. OLS is a continuous function of the averages \(P_NZZ^\top\) and \(P_NZY\), which converge to \(A\) and \(\E(ZY) = A\beta_0\) by the law of large numbers and correct specification; hence \(\hat\beta \pto \beta_0\). Since \(P_Nz_d \pto Pz_d\), the estimation error tends to zero by Slutsky’s theorem. The sampling error is the average of i.i.d. mean-zero terms, so it tends to zero by the law of large numbers. \(\square\)
The two assumptions have two different jobs. Identification (Theorem 5.2) makes \(\E\{\mu(1, X) - \mu(0, X)\}\) the ATE. The outcome model makes \(\hat\tau_{\mathrm{reg}}\) consistent for it, because correct specification is what makes OLS converge to the \(\beta_0\) that defines \(\mu\). Under model (a), \(z_d\) is constant, so the sampling error vanishes and \(\hat\tau_{\mathrm{reg}} = \hat\tau_{\mathrm{OLS}}\).
5.3.7 Asymptotic Normality and Variance Estimation
The decomposition Equation 5.12 reduces the problem to the asymptotics of the nuisance estimator.
Lemma 5.2 (OLS under a Correct Linear Model) Let \(Y = Z^\top\beta_0 + \varepsilon\) with \(\E(\varepsilon \mid Z) = 0\), i.i.d. data, \(\E\|Z\|^4 < \infty\), \(\E\varepsilon^4 < \infty\), and \(A = \E(ZZ^\top)\) nonsingular; \(\mathrm{Var}(\varepsilon \mid Z)\) may depend on \(Z\). With \(B = \E(\varepsilon^2ZZ^\top)\), \[\sqrt N(\hat\beta - \beta_0) = A^{-1}\sqrt N\,P_N(Z\varepsilon) + o_p(1) \;\xrightarrow{d}\; N\bigl(0,\, A^{-1}BA^{-1}\bigr).\]
Proof. Substituting \(Y = Z^\top\beta_0 + \varepsilon\) into \(\hat\beta\) gives the exact identity \(\hat\beta - \beta_0 = (P_NZZ^\top)^{-1}P_N(Z\varepsilon)\). By the law of large numbers, \((P_NZZ^\top)^{-1} \pto A^{-1}\). The terms \(Z_i\varepsilon_i\) are i.i.d. with mean \(\E\{Z\,\E(\varepsilon \mid Z)\} = 0\) and variance \(B\), so \(\sqrt N P_N(Z\varepsilon) \xrightarrow{d} N(0, B)\) by the central limit theorem; Slutsky’s theorem completes the proof. \(\square\)
The matrix \(A^{-1}BA^{-1}\) is the sandwich variance. It is estimated by \(\hat A^{-1}\hat B\hat A^{-1}\) with \(\hat A = P_NZZ^\top\), \(\hat B = P_N\hat\varepsilon^2ZZ^\top\), and \(\hat\varepsilon_i = Y_i - Z_i^\top\hat\beta\). Under homoskedasticity, \(B = \sigma^2A\) and it reduces to the classical \(\sigma^2A^{-1}\).
Theorem 5.3 (Asymptotic Normality of the Regression Estimator) Under the conditions of Proposition 5.4 and Lemma 5.2, assume also \(\E\|z_d(X)\|^2 < \infty\), and let \(c = Pz_d\) and \[\varphi_i = \bigl\{z_d(X_i)^\top\beta_0 - \tau_{\mathrm{ATE}}\bigr\} + c^\top A^{-1}Z_i\varepsilon_i . \tag{5.13}\] Then \(\sqrt N(\hat\tau_{\mathrm{reg}} - \tau_{\mathrm{ATE}}) = \sqrt N\,P_N\varphi + o_p(1) \xrightarrow{d} N\{0, \mathrm{Var}(\varphi)\}\), and \[\mathrm{Var}(\varphi) = \underbrace{\mathrm{Var}\bigl\{z_d(X)^\top\beta_0\bigr\}}_{\text{sampling error}} + \underbrace{c^\top A^{-1}BA^{-1}c}_{\text{estimation error}} . \tag{5.14}\]
Proof. Multiply Equation 5.12 by \(\sqrt N\). By Lemma 5.2, \(\sqrt N(\hat\beta - \beta_0) = O_p(1)\), and \(P_Nz_d = c + o_p(1)\), so the estimation error is \(c^\top A^{-1}\sqrt N P_N(Z\varepsilon) + o_p(1)\). Adding the sampling error gives \(\sqrt N P_N\varphi + o_p(1)\), and the central limit theorem applies to the i.i.d. mean-zero terms \(\varphi_i\), which have finite variance by the moment conditions. For the variance, the correct conditional-mean model gives \(\E(\varepsilon \mid T, X) = 0\), and both \(z_d(X)\) and \(Z\) are functions of \((T, X)\), so the cross term is \[\E\bigl[\{z_d(X)^\top\beta_0 - \tau_{\mathrm{ATE}}\}\,c^\top A^{-1}Z\, \E(\varepsilon \mid T, X)\bigr] = 0 . \qquad \square\]
The moment condition on \(z_d(X)\) is not implied by those on \(Z\) for a general \(z(t, x)\): \(Z\) involves only the received treatment, and a component that matters only under the unreceived treatment can be large where that treatment is rare. For models (a) and (b), \(z_d(X)\) is a subvector of \(Z\) up to sign, or constant, and the condition is automatic.
The function \(\varphi\) is the influence function of \(\hat\tau_{\mathrm{reg}}\): to first order, the estimator behaves like the sample mean of \(\varphi_i\). Its two components are the two errors of Equation 5.12, to first order, and they are uncorrelated because the correct model gives \(\E(\varepsilon \mid T, X) = 0\). The exact finite-sample terms of Equation 5.12 need not be uncorrelated.
Variance estimation. With \(\hat c = P_Nz_d\) and the plug-in quantities above, \[\widehat{\mathrm{Var}}(\hat\tau_{\mathrm{reg}}) = \frac1N\Bigl[P_N\bigl\{z_d(X)^\top\hat\beta - \hat\tau_{\mathrm{reg}}\bigr\}^2 + \hat c^\top\hat A^{-1}\hat B\hat A^{-1}\hat c\Bigr] \tag{5.15}\] is consistent for \(N^{-1}\mathrm{Var}(\varphi)\) by the law of large numbers and continuity, and \(\hat\tau_{\mathrm{reg}} \pm 1.96\,\widehat{\mathrm{SE}}\) is an asymptotic 95% confidence interval. An asymptotically equivalent estimator is \(N^{-2}\sum_i\hat\varphi_i^2\), with \(\hat\varphi_i\) the plug-in influence values. The two are not identical in finite samples: the normal equations do not make the sample cross term between the two components vanish, but it is \(o_p(N^{-1})\).
5.3.8 What Can Go Wrong
Misspecification. If the true \(\mu(t, x)\) is nonlinear and a linear model is fitted, OLS converges to the best linear approximation, not to the \(\beta_0\) that defines \(\mu\), and \(\hat\tau_{\mathrm{reg}}\) absorbs part of the functional-form error. The bias can have either sign.
Extrapolation. Step 2 predicts \(\hat\mu(0, X_i)\) for treated units whose \(X_i\) may lie outside the support of the controls. There the model, not comparable units, supplies the answer. The overlap condition of Chapter 6 formalizes when extrapolation is unavoidable.
Flexible learners and cross-fitting. A more flexible outcome estimator can reduce approximation bias, at the price of more estimation variability. Cross-fitting separates fitting from evaluation, but alone it guarantees neither consistency nor \(\sqrt N\) inference. Chapters 11 and 12 develop estimators whose leading error is less sensitive to errors in the outcome and propensity-score functions.
5.4 Regression Estimation under Complete Randomization
Section 5.3 justified the regression estimator \(\hat\tau_{\mathrm{reg}}\) by a correct outcome model, with the \(N\) units treated as an i.i.d. sample. This section returns to the design-based setting of Section 5.1: the \(N\) units form a fixed finite population \(\mathcal F_N\), \(N_1\) are treated by the CRD, and the only randomness is the assignment. In this setting the difference in means is already exactly unbiased for \(\bar\tau_N\) with no model. If we fit the regression anyway, what happens when the linear model is wrong?
The answer, due to Lin (2013), is that the design, not the model, justifies \(\hat\tau_{\mathrm{reg}}\). It is design consistent whether or not the linear model holds, and asymptotically it cannot be less precise than \(\hat\tau_{\mathrm{DIM}}\). The proof is the survey-sampling argument for the regression estimator, applied to each arm. Table 5.7 places this section beside Section 5.3; it is easy to carry an assumption from one framework into the other by mistake.
| Framework | Target | Source of randomness | Role of the outcome regression |
|---|---|---|---|
| Finite-population CRD (this section) | \(\bar\tau_N\) | Assignment | A working approximation used for precision; correct specification is not needed for the design asymptotics |
| I.i.d. observational data (Section 5.3) | \(\tau_{\mathrm{ATE}}\) | Sampling of \((Y, T, X)\) | Correct specification gives the consistency and variance results proved there |
| I.i.d. units followed by the CRD | \(\tau_{\mathrm{ATE}}\) | Sampling and assignment | Inference must account for assignment uncertainty and for variation in the sampled population |
Example 5.7 (Incentives and Services for College Achievement) Angrist et al. (2009) randomized first-year undergraduates at a Canadian university to support services, financial incentives (awards of $1,000 to $5,000 for meeting a target grade), both, or neither. Following Lin (2013, Table 2), we use the men in two of the four arms, with first-year GPA observed: services and incentives (\(N_1 = 58\)) versus services only (\(N_0 = 99\)). The contrast is the effect of adding incentives when services are offered. The outcome is first-year GPA, and the covariate is high-school GPA, observed for everyone. Randomization makes the two groups comparable on average, but in the realized assignment their high-school GPAs differ a little. Can we use high-school GPA to adjust, and does it help? Section 5.4.6 returns to the data.
5.4.1 Two Regression Estimators and Three Forms
Let \(\mathbf x_i\) be the covariate vector of unit \(i\), including an intercept, and assume each arm’s design matrix has full column rank. Fit OLS separately in the two arms: \[\hat{\mathbf B}_1 = \Bigl(\sum_{i=1}^N T_i\mathbf x_i\mathbf x_i^\top\Bigr)^{-1}\sum_{i=1}^N T_i\mathbf x_iY_i, \qquad \hat{\mathbf B}_0 = \Bigl(\sum_{i=1}^N (1 - T_i)\mathbf x_i\mathbf x_i^\top\Bigr)^{-1}\sum_{i=1}^N (1 - T_i)\mathbf x_iY_i .\] The estimator of \(\bar\tau_N\) is the difference of two regression estimators, \[\hat\tau_{\mathrm{reg}} = \hat\mu_{1,\mathrm{reg}} - \hat\mu_{0,\mathrm{reg}}, \qquad \hat\mu_{t,\mathrm{reg}} = \frac1N\sum_{i=1}^N\mathbf x_i^\top\hat{\mathbf B}_t = \bar{\mathbf x}^\top\hat{\mathbf B}_t, \tag{5.16}\] which is the estimator of Section 5.3 under model (b). Write \(\mathbf x_i = (1, \mathbf u_i^\top)^\top\), so that \(\mathbf u_i\) collects the nonconstant covariates. By Proposition 5.3, \(\hat\tau_{\mathrm{reg}}\) is also the coefficient on \(T\) in the interacted regression of \(Y\) on \(\bigl(1, T, (\mathbf u - \bar{\mathbf u})^\top, T(\mathbf u - \bar{\mathbf u})^\top\bigr)\); the intercept coordinate is omitted from the centered and interaction terms, which would otherwise contain zero columns. With one covariate, \(\mathbf x_i = (1, X_i)^\top\) and \(\hat{\mathbf B}_1 = (\hat\alpha_1, \hat b_1)^\top\), \[\hat\mu_{1,\mathrm{reg}} = \hat\alpha_1 + \bar X\hat b_1 = \hat{\bar Y}_1 + \bigl(\bar X - \hat{\bar X}_1\bigr)\hat b_1,\] where \(\hat{\bar Y}_1\) and \(\hat{\bar X}_1\) are treated-arm means. When the treated arm has a higher covariate mean than the population and the slope is positive, the adjustment lowers the estimated mean.
Lemma 5.3 (Three Forms of the Regression Estimator) Let \(\hat e_i = Y_i - \mathbf x_i^\top\hat{\mathbf B}_1\) be the fitted residuals of the treated units. Then \[\underbrace{\frac1N\sum_{i=1}^N\mathbf x_i^\top\hat{\mathbf B}_1}_{\text{projection form}} = \underbrace{\frac1N\sum_{i=1}^N\mathbf x_i^\top\hat{\mathbf B}_1 + \frac1{N_1}\sum_{i=1}^NT_i\hat e_i}_{\text{projection}\,+\,\text{expansion estimator of the residual mean}} = \underbrace{\hat{\bar Y}_1 + (\bar{\mathbf x} - \hat{\bar{\mathbf x}}_1)^\top\hat{\mathbf B}_1}_{\text{direct}\,+\,\text{augmentation}} . \tag{5.17}\]
Proof. The normal equations of the treated-arm fit are \(\sum_iT_i\mathbf x_i\hat e_i = \mathbf 0\). The first component of \(\mathbf x_i\) is \(1\), so \(\sum_iT_i\hat e_i = 0\) and the added term is zero. Expanding \(\hat e_i\) gives \(N_1^{-1}\sum_iT_i\hat e_i = \hat{\bar Y}_1 - \hat{\bar{\mathbf x}}_1^\top\hat{\mathbf B}_1\), which yields the last form. \(\square\)
The second form, projection plus the expansion (Horvitz–Thompson) estimator of the residual mean, is the generalized regression (GREG) estimator of survey sampling. Its residual term vanishes here because the OLS normal equations include an intercept and the expansion weights are constant within the arm. Under unequal-probability sampling, whether the analogous term vanishes depends on how the regression coefficients and the expansion weights are defined; when it does not, it is the term that makes GREG design consistent. The third form reads the estimator as the direct estimator \(\hat{\bar Y}_1\) corrected by an augmentation term, the covariate imbalance times the slope.
5.4.2 The Survey-Sampling Setup and the Oracle Difference Estimator
Let \(\mathcal F_N = \{(Y_i(1), Y_i(0), \mathbf x_i) : i = 1, \ldots, N\}\) be fixed. Under the CRD, the treated arm is a simple random sample (SRS) of size \(N_1\) without replacement, and the control arm is its complement. Arm 1 is therefore a survey: \(Y_i(1)\) is observed on the SRS, \(\mathbf x_i\) is known for all \(N\) units, and \(\bar{\mathbf x} = N^{-1}\sum_i\mathbf x_i\) is a known population mean. In survey terms, \(\hat\mu_{1,\mathrm{reg}}\) is the regression estimator of \(\bar Y(1)\). All expectations below are over the assignment, given \(\mathcal F_N\), and no model for \(Y_i(t)\) is assumed.
Definition 5.2 (Difference Estimator) For a fixed vector \(\mathbf B\), not estimated from the data, the difference estimator of \(\bar Y(1)\) is \[\hat\mu_{1,\mathrm{diff}}(\mathbf B) = \hat{\bar Y}_1 + (\bar{\mathbf x} - \hat{\bar{\mathbf x}}_1)^\top\mathbf B = \frac1{N_1}\sum_{i=1}^NT_i\bigl\{Y_i(1) - (\mathbf x_i - \bar{\mathbf x})^\top\mathbf B\bigr\}.\]
The difference estimator is an SRS mean of the fixed numbers \(u_i = Y_i(1) - (\mathbf x_i - \bar{\mathbf x})^\top\mathbf B\), whose population mean is \(\bar Y(1)\). It is therefore exactly design unbiased, with \(\mathrm{Var}\{\hat\mu_{1,\mathrm{diff}}(\mathbf B) \mid \mathcal F_N\} = (N_1^{-1} - N^{-1})S_u^2\) by Theorem 5.1. The family contains the plain mean (\(\mathbf B = \mathbf 0\)). Within it, \(S_u^2\) is the variance of \(Y(1)\) around the linear function \((\mathbf x - \bar{\mathbf x})^\top\mathbf B\), which is minimized by the population OLS coefficient (Cochran 1977) \[\mathbf B_1 = \Bigl(\sum_{i=1}^N\mathbf x_i\mathbf x_i^\top\Bigr)^{-1}\sum_{i=1}^N\mathbf x_iY_i(1). \tag{5.18}\] The oracle difference estimator \(\hat\mu_{1,\mathrm{diff}}(\mathbf B_1)\) is the most precise member of the family. It is not available: \(\mathbf B_1\) needs every \(Y_i(1)\). The regression estimator is the same formula with \(\mathbf B_1\) replaced by its treated-arm estimate, \(\hat\mu_{1,\mathrm{reg}} = \hat\mu_{1,\mathrm{diff}}(\hat{\mathbf B}_1)\). The oracle is a device for the proof: the next subsection shows that \(\hat\mu_{1,\mathrm{reg}}\) is asymptotically equivalent to it, so the design properties of the regression estimator are those of an SRS mean of residuals. Define \(\mathbf B_0\) in the same way from the \(Y_i(0)\), and the population residuals \[e_i(t) = Y_i(t) - \mathbf x_i^\top\mathbf B_t, \qquad \sum_{i=1}^Ne_i(t) = 0, \qquad \sum_{i=1}^N\mathbf x_ie_i(t) = \mathbf 0,\] where the last two identities are the population normal equations.
5.4.3 Linearization and Design Consistency
Asymptotic statements in the design-based framework refer to a sequence of finite populations with \(N \to \infty\).
Complete randomization alone does not give stable OLS inverses or the rates used below; the conditions do. We write \(R_N = O_p(a_N)\) if for every \(\varepsilon > 0\) there is \(M < \infty\) with \(P(|R_N| > Ma_N) < \varepsilon\) for all large \(N\), and \(R_N = o_p(a_N)\) if \(R_N/a_N \pto 0\).
Lemma 5.4 (Rates for SRS Means) Under the regularity conditions, \(\hat{\bar{\mathbf x}}_1 - \bar{\mathbf x} = O_p(N^{-1/2})\) and \(\hat{\mathbf B}_1 - \mathbf B_1 = O_p(N^{-1/2})\), and likewise for arm 0.
Proof. For one coordinate \(x\), the treated-arm mean is an SRS mean, so \(\E(\hat{\bar x}_1 \mid \mathcal F_N) = \bar x\) and \(\mathrm{Var}(\hat{\bar x}_1 \mid \mathcal F_N) = (N_1^{-1} - N^{-1})S_x^2 \le S_x^2/N_1\). By Chebyshev’s inequality, \[P\bigl(|\hat{\bar x}_1 - \bar x| > MN^{-1/2} \bigm| \mathcal F_N\bigr) \le \frac{N}{N_1}\cdot\frac{S_x^2}{M^2},\] which is at most \(C/M^2\) for large \(N\) because \(N_1/N \to \rho > 0\) and \(S_x^2\) is bounded; choose \(M\) with \(C/M^2 < \varepsilon\). The same argument applies to every SRS mean of a variable with bounded second moment, in particular to the treated-arm means of \(\mathbf x_i\mathbf x_i^\top\) and \(\mathbf x_iY_i(1)\), whose fourth-moment bounds give bounded variances. \(\hat{\mathbf B}_1\) and \(\mathbf B_1\) are the same smooth function of these means, evaluated at SRS means and at population means respectively, and the limiting Gram matrix is nonsingular, so \(\hat{\mathbf B}_1 - \mathbf B_1 = O_p(N^{-1/2})\). \(\square\)
Theorem 5.4 (Asymptotic Equivalence and Design Consistency) Under the CRD and the regularity conditions, with no model for the potential outcomes, \[\hat\mu_{1,\mathrm{reg}} = \underbrace{\frac1N\sum_{i=1}^N\mathbf x_i^\top\mathbf B_1 + \frac1{N_1}\sum_{i=1}^NT_i\bigl(Y_i - \mathbf x_i^\top\mathbf B_1\bigr)}_{\hat\mu_{1,\mathrm{diff}}(\mathbf B_1)} + \underbrace{(\bar{\mathbf x} - \hat{\bar{\mathbf x}}_1)^\top(\hat{\mathbf B}_1 - \mathbf B_1)}_{O_p(N^{-1})} , \tag{5.19}\] and likewise for arm 0. Consequently \[\hat\tau_{\mathrm{reg}} - \bar\tau_N = \underbrace{\frac1{N_1}\sum_{i=1}^NT_ie_i(1) - \frac1{N_0}\sum_{i=1}^N(1 - T_i)e_i(0)}_{\mathrm{DIME}(e)} + O_p(N^{-1}), \tag{5.20}\] and \(\hat\tau_{\mathrm{reg}} \pto \bar\tau_N\) in the sense that \(\hat\tau_{\mathrm{reg}} - \bar\tau_N \pto 0\).
Proof. By Lemma 5.3, adding the zero term \(N_1^{-1}\sum_iT_i(Y_i - \mathbf x_i^\top\hat{\mathbf B}_1)\) to \(\hat\mu_{1,\mathrm{reg}}\) changes nothing. Replacing \(\hat{\mathbf B}_1\) by \(\mathbf B_1\) in both terms and collecting the difference gives Equation 5.19; the remainder is a product of two \(O_p(N^{-1/2})\) factors by Lemma 5.4. The intercept and the population normal equation give \(N^{-1}\sum_i\mathbf x_i^\top\mathbf B_1 = \bar Y(1)\), and \(Y_i = Y_i(1)\) on the treated arm, so \(\hat\mu_{1,\mathrm{diff}}(\mathbf B_1) = \bar Y(1) + N_1^{-1}\sum_iT_ie_i(1)\). Subtracting the arm-0 version gives Equation 5.20. Each \(e_i(t)\) is fixed with \(\sum_ie_i(t) = 0\), and \(\E(T_i \mid \mathcal F_N) = N_1/N\), so \(\E\{N_1^{-1}\sum_iT_ie_i(1) \mid \mathcal F_N\} = N^{-1}\sum_ie_i(1) = 0\), and its variance is \((N_1^{-1} - N^{-1})S_{e(1)}^2 = O(N^{-1})\). By Chebyshev’s inequality the term tends to zero in probability; the same holds for arm 0. \(\square\)
We write \(\hat\tau_{\mathrm{reg},\ell} = \bar\tau_N + \mathrm{DIME}(e)\) for the linearized estimator. It is the difference of two oracle difference estimators, it is exactly design unbiased, and \(\hat\tau_{\mathrm{reg}} - \hat\tau_{\mathrm{reg},\ell} = O_p(N^{-1}) = o_p(N^{-1/2})\), so the two share their first-order distribution. The model appears nowhere in the argument: \(\mathbf B_1\) is a population projection, and \(Y_i(1)\) need not be linear in \(\mathbf x_i\). Randomization justifies the estimator, not the model.
5.4.4 Design Variance and Conservative Estimation
The linearized estimator is a difference in means of the fixed residuals \(e_i(1)\) and \(e_i(0)\), so Theorem 5.1 applies with \(e(t)\) in place of \(Y(t)\).
Proposition 5.5 (Design Variance of the Regression Estimator) Let \(S_{e(t)}^2\) and \(S_{e(1)-e(0)}^2\) be the finite-population variances (divisor \(N - 1\)) of the residuals and of their difference. Then \[\mathrm{Var}(\hat\tau_{\mathrm{reg},\ell} \mid \mathcal F_N) = \frac{S_{e(1)}^2}{N_1} + \frac{S_{e(0)}^2}{N_0} - \frac{S_{e(1)-e(0)}^2}{N} = \Bigl(\frac1{N_1} - \frac1N\Bigr)S_{e(1)}^2 + \Bigl(\frac1{N_0} - \frac1N\Bigr)S_{e(0)}^2 + \frac2NS_{e(0),e(1)} . \tag{5.21}\]
As in Neyman’s theorem, \(S_{e(1)-e(0)}^2\) is not identified, because no unit shows both residuals. Omitting its negative contribution gives the asymptotically conservative estimator \[\hat V_{\mathrm{reg}} = \frac{s_{\hat e,1}^2}{N_1} + \frac{s_{\hat e,0}^2}{N_0}, \tag{5.22}\] where \(s_{\hat e,t}^2\) is the sample variance of the fitted residuals in arm \(t\). Dropping the covariance term \(2N^{-1}S_{e(0),e(1)}\) of the second form instead would not be conservative, because that term can have either sign: if \(e_i(1) = e_i(0)\) for every unit, it is positive, and dropping it understates the variance. Under the regularity conditions, the sandwich variance of the \(T\) coefficient in the interacted regression also estimates the first-order asymptotic design variance consistently or conservatively, with the same leading \(O(N^{-1})\) term as \(\hat V_{\mathrm{reg}}\) (Lin 2013). Equation Equation 5.21 is exact for the linearized estimator; turning the distributional approximation into a statement about the actual finite-sample variance of \(\hat\tau_{\mathrm{reg}}\) requires moment control of the remainder. The target of this inference is the finite-population \(\bar\tau_N\).
5.4.5 Regression Adjustment Cannot Increase the Asymptotic Variance
The residuals \(e_i(t)\) have smaller spread than the outcomes when \(\mathbf x\) predicts \(Y(t)\). The precise statement uses two facts about simple random sampling.
Lemma 5.5 (Covariances of SRS Means) For fixed population vectors \(u\) and \(v\) with population covariance \(S_{u,v} = (N - 1)^{-1}\sum_i(u_i - \bar u)(v_i - \bar v)\), let \(\hat{\bar u}_t\) denote the arm-\(t\) mean of \(u\). Under the CRD, \[\mathrm{Cov}(\hat{\bar u}_1, \hat{\bar v}_1 \mid \mathcal F_N) = \Bigl(\frac1{N_1} - \frac1N\Bigr)S_{u,v}, \qquad \mathrm{Cov}(\hat{\bar u}_1, \hat{\bar v}_0 \mid \mathcal F_N) = -\frac1NS_{u,v} .\]
Proof. The first identity is the polarization of the variance calculation in the proof of Theorem 5.1, using the assignment covariances Equation 5.5. For the second, \(\sum_i(1 - T_i)v_i = \sum_iv_i - \sum_iT_iv_i\), and \(\sum_iv_i\) is fixed, so \(\mathrm{Cov}(\hat{\bar u}_1, \hat{\bar v}_0) = -(N_1/N_0)\,\mathrm{Cov}(\hat{\bar u}_1, \hat{\bar v}_1) = -(N_1/N_0)(N_1^{-1} - N^{-1})S_{u,v} = -S_{u,v}/N\). \(\square\)
Theorem 5.5 (No Loss of Precision) Under the CRD, \[\mathrm{Var}(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_N) \;\ge\; \mathrm{Var}(\hat\tau_{\mathrm{reg},\ell} \mid \mathcal F_N),\] an exact inequality for the linearized estimator and hence an asymptotic one for \(\hat\tau_{\mathrm{reg}}\) under the regularity conditions. Equality holds if and only if the covariates carry no linear information about \((1 - \rho_N)Y(1) + \rho_NY(0)\), with \(\rho_N = N_1/N\).
Proof. Write \(Y_i(t) = g_i(t) + e_i(t)\) with \(g_i(t) = \mathbf x_i^\top\mathbf B_t\). Since \(\bar e(t) = 0\), \(\bar Y(t) = \bar g(t)\), and \[\hat\tau_{\mathrm{DIM}} - \bar\tau_N = \underbrace{\mathrm{DIME}(e)}_{=\,\hat\tau_{\mathrm{reg},\ell} - \bar\tau_N} + \underbrace{\bigl\{\mathrm{DIME}(g) - \bar g(1) + \bar g(0)\bigr\}}_{\text{mean zero}},\] where \(\mathrm{DIME}(g) = N_1^{-1}\sum_iT_ig_i(1) - N_0^{-1}\sum_i(1 - T_i)g_i(0)\). The population normal equations give \(\sum_ie_i(t)\mathbf x_i = \mathbf 0\) and \(\sum_ie_i(t) = 0\), so \(S_{e(t),\,g(s)} = 0\) for all \(t, s \in \{0, 1\}\). By Lemma 5.5, each of the four covariances between an arm mean of \(e(t)\) and an arm mean of \(g(s)\) is a multiple of \(S_{e(t),g(s)}\) and hence zero, so \(\mathrm{Cov}\{\mathrm{DIME}(e), \mathrm{DIME}(g) \mid \mathcal F_N\} = 0\) and \[\mathrm{Var}(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_N) = \mathrm{Var}\{\mathrm{DIME}(e) \mid \mathcal F_N\} + \mathrm{Var}\{\mathrm{DIME}(g) \mid \mathcal F_N\}.\] The second term is nonnegative. By the second form of Equation 5.3 with \(\rho_N = N_1/N\), it equals \(S^2_{(1-\rho_N)g(1) + \rho_Ng(0)}/\{N\rho_N(1 - \rho_N)\}\), which is zero if and only if \(\mathbf x^\top\{(1 - \rho_N)\mathbf B_1 + \rho_N\mathbf B_0\}\), the population projection of \((1 - \rho_N)Y(1) + \rho_NY(0)\) on \(\mathbf x\), is constant. \(\square\)
The gain is the design variance of \(\mathrm{DIME}(g)\), the part of the outcome explained by \(\mathbf x\). Adjustment buys precision, never validity; it is useless when the covariates do not predict the outcomes.
5.4.6 The Survey-Sampling Precedent and the Example Revisited
Example 5.8 (Incentives and Services, Continued) For the 157 men of Example 5.7, Lin (2013, Table 2) reports:
| Estimator | Estimate | Standard error |
|---|---|---|
| Difference in means \(\hat\tau_{\mathrm{DIM}}\) | \(-0.036\) | \(0.158\) |
| Regression estimator \(\hat\tau_{\mathrm{reg}}\) | \(-0.081\) | \(0.146\) |
The standard errors are uncorrected sandwich standard errors. Both estimators target \(\bar\tau_N\), the effect for these 157 men of offering incentives in addition to services. The estimated standard error falls from 0.158 to 0.146 because the regression estimator’s variance uses residual variances; a smaller estimated standard error in one data set is consistent with, but does not by itself prove, a smaller true variance. The gain is small because high-school GPA predicts first-year GPA only weakly. The two estimates differ by 0.045, well within one standard error: the adjustment corrects a small chance imbalance in high-school GPA.
5.5 Calibration Weighting
Section 5.4 studied \(\hat\tau_{\mathrm{reg}}\) through its coefficients \(\hat{\mathbf B}_t\). This section studies the same estimator through its weights. Each arm’s linear regression estimator, fitted by least squares with an intercept, is a weighted mean of the outcomes in its arm, and the weights have two properties: they reproduce the full-sample covariate mean exactly, and among all weights that do so they are the closest to uniform. This is the calibration view of survey sampling (Deville and Särndal 1992). It explains what regression adjustment does to the covariates, and it is the bridge to stratification (Section 5.6) and to propensity-score weighting (Chapter 6).
The setting is that of Section 5.4: \(\mathbf x_i\) includes an intercept, each arm’s design matrix has full column rank, and each weight vector contains only the units of its own arm. We treat arm 1; arm 0 is identical. Write \[\hat{\mathbf M}_1 = \sum_{i=1}^NT_i\mathbf x_i\mathbf x_i^\top, \qquad \hat{\bar{\mathbf x}}_1 = \frac1{N_1}\sum_{i=1}^NT_i\mathbf x_i, \qquad \bar{\mathbf x} = \frac1N\sum_{i=1}^N\mathbf x_i .\]
5.5.1 The Regression Estimator as a Weighted Sum
Substituting \(\hat{\mathbf B}_1 = \hat{\mathbf M}_1^{-1}\sum_iT_i\mathbf x_iY_i\) into \(\hat\mu_{1,\mathrm{reg}} = \bar{\mathbf x}^\top\hat{\mathbf B}_1\) gives a linear function of the treated outcomes, \[\hat\mu_{1,\mathrm{reg}} = \sum_{i=1}^NT_i\,\omega_{1i}Y_i, \qquad \omega_{1i} = \bar{\mathbf x}^\top\hat{\mathbf M}_1^{-1}\mathbf x_i , \tag{5.23}\] and likewise \(\hat\mu_{0,\mathrm{reg}} = \sum_i(1 - T_i)\,\omega_{0i}Y_i\) with \(\omega_{0i} = \bar{\mathbf x}^\top\{\sum_j(1 - T_j)\mathbf x_j\mathbf x_j^\top\}^{-1}\mathbf x_i\). The weights depend on \(\mathbf x\) and \(T\) only, not on \(Y\), so they can be computed and inspected before any outcome is seen. With \(\mathbf x_i = 1\) alone, they are \(1/N_1\) and \(1/N_0\), and \(\hat\tau_{\mathrm{reg}} = \hat\tau_{\mathrm{DIM}}\).
Proposition 5.6 (The Regression Weights Satisfy the Calibration Constraint) \[\sum_{i=1}^NT_i\,\omega_{1i}\mathbf x_i = \bar{\mathbf x}, \qquad \sum_{i=1}^N(1 - T_i)\,\omega_{0i}\mathbf x_i = \bar{\mathbf x}. \tag{5.24}\]
Proof. \(\sum_iT_i\omega_{1i}\mathbf x_i = \sum_iT_i\mathbf x_i\mathbf x_i^\top\hat{\mathbf M}_1^{-1}\bar{\mathbf x} = \hat{\mathbf M}_1\hat{\mathbf M}_1^{-1}\bar{\mathbf x} = \bar{\mathbf x}\); the control arm is the same. \(\square\)
Each weighted arm reproduces the full-sample covariate mean exactly. The intercept is one of the covariates, so the weights sum to one. Condition Equation 5.24 is called the covariate-balancing or calibration constraint. Balance of these means does not by itself establish exchangeability: it concerns only the covariate functions included in \(\mathbf x\). Many weight vectors satisfy Equation 5.24 when \(N_1 > p\). The rest of this section asks which one the regression estimator uses, and why.
5.5.2 Which Weights? A Working Model
Consider all weighted estimators \(\hat\mu_{\boldsymbol\omega} = \sum_iT_i\omega_iY_i\) of \(\bar Y(1)\) with normalized weights, \(\sum_iT_i\omega_i = 1\); the plain mean has \(\omega_i = 1/N_1\). To compare them, adopt the working model \[Y_i(1) = \mathbf x_i^\top\boldsymbol\beta + \varepsilon_i, \qquad \E(\boldsymbol\varepsilon \mid \mathbf X, \mathbf T) = \mathbf 0, \qquad \mathrm{Var}(\boldsymbol\varepsilon \mid \mathbf X, \mathbf T) = \sigma^2\mathbf I_N ,\] and define the calibration residual \(\mathbf h(\boldsymbol\omega) = \bar{\mathbf x} - \sum_iT_i\omega_i\mathbf x_i\), the target minus the achieved covariate mean. For fixed weights, \[\hat\mu_{\boldsymbol\omega} - \bar Y(1) = -\mathbf h(\boldsymbol\omega)^\top\boldsymbol\beta + \sum_{i=1}^N\Bigl(T_i\omega_i - \frac1N\Bigr)\varepsilon_i ,\] because the target \(\bar Y(1)\) contains the same errors as the estimator. Hence, for normalized weights, \[\mathrm{MSE}\bigl\{\hat\mu_{\boldsymbol\omega} - \bar Y(1) \mid \mathbf X, \mathbf T\bigr\} = \bigl\{\mathbf h(\boldsymbol\omega)^\top\boldsymbol\beta\bigr\}^2 + \sigma^2\Bigl(\sum_{i=1}^NT_i\omega_i^2 - \frac1N\Bigr), \tag{5.25}\] since \(\sum_i(T_i\omega_i - N^{-1})^2 = \sum_iT_i\omega_i^2 - 2N^{-1}\sum_iT_i\omega_i + N^{-1}\). The model supplies the three ingredients of the optimization.
- The constraint. \(\boldsymbol\beta\) is unknown. The estimator is model-unbiased for every \(\boldsymbol\beta\) if and only if \(\mathbf h(\boldsymbol\omega) = \mathbf 0\): the weighted covariate mean must equal the known \(\bar{\mathbf x}\).
- The distance. With the bias removed, minimizing the MSE means minimizing \(\sum_iT_i\omega_i^2\).
- The anchor. On the normalized set, \(\sum_iT_i\omega_i^2 = \sum_iT_i(\omega_i - 1/N_1)^2 + 1/N_1\), so the criterion is the Euclidean distance to the anchor, the initial weights \(1/N_1\).
Only the form of the model is used; \(\boldsymbol\beta\) is never estimated. The model motivates the criterion; the identity between calibration and OLS in the next subsection does not need it to be true, and the design-based results of Section 5.4 do not use it at all.
5.5.3 \(L_2\) Calibration and the Lagrange Multiplier Method
The calibration problem for arm 1 is \[\min_{\boldsymbol\omega}\; D(\boldsymbol\omega) = \frac12\sum_{i=1}^NT_i\Bigl(\omega_i - \frac1{N_1}\Bigr)^2 \quad\text{subject to}\quad \mathbf h(\boldsymbol\omega) = \bar{\mathbf x} - \sum_{i=1}^NT_i\omega_i\mathbf x_i = \mathbf 0, \tag{5.26}\] with Lagrangian \(\mathcal L(\boldsymbol\omega, \boldsymbol\lambda) = D(\boldsymbol\omega) + \boldsymbol\lambda^\top\mathbf h(\boldsymbol\omega)\) and multiplier \(\boldsymbol\lambda \in \mathbb R^p\). The multiplier method turns a problem in \(N_1\) unknowns into one in \(p\) unknowns.
- Minimize in \(\boldsymbol\omega\) for fixed \(\boldsymbol\lambda\). Setting \(\partial\mathcal L/\partial\omega_i = 0\) gives, unit by unit, the linear weight form \[\omega_i^\star(\boldsymbol\lambda) = \frac1{N_1} + \boldsymbol\lambda^\top\mathbf x_i :\] the anchor plus a linear function of \(\mathbf x_i\), with the same \(p\) coefficients for every unit.
- Maximize the dual. Substituting, \(D(\boldsymbol\omega^\star) = \tfrac12\boldsymbol\lambda^\top\hat{\mathbf M}_1\boldsymbol\lambda\) and \(\mathbf h(\boldsymbol\omega^\star) = \bar{\mathbf x} - \hat{\bar{\mathbf x}}_1 - \hat{\mathbf M}_1\boldsymbol\lambda\), so the dual function is the concave quadratic \[\ell(\boldsymbol\lambda) = \mathcal L\{\boldsymbol\omega^\star(\boldsymbol\lambda), \boldsymbol\lambda\} = \boldsymbol\lambda^\top(\bar{\mathbf x} - \hat{\bar{\mathbf x}}_1) - \frac12\boldsymbol\lambda^\top\hat{\mathbf M}_1\boldsymbol\lambda ,\] whose gradient is the calibration residual \(\mathbf h\{\boldsymbol\omega^\star(\boldsymbol\lambda)\}\).
- Return. \(\nabla\ell = \mathbf 0\) gives \(\hat{\boldsymbol\lambda} = \hat{\mathbf M}_1^{-1}(\bar{\mathbf x} - \hat{\bar{\mathbf x}}_1)\) and the regression weights \[\hat\omega_i = \frac1{N_1} + (\bar{\mathbf x} - \hat{\bar{\mathbf x}}_1)^\top\hat{\mathbf M}_1^{-1}\mathbf x_i = \bar{\mathbf x}^\top\hat{\mathbf M}_1^{-1}\mathbf x_i . \tag{5.27}\] The second form holds because the first column of \(\hat{\mathbf M}_1\) is \(\sum_iT_i\mathbf x_i = N_1\hat{\bar{\mathbf x}}_1\), so \(\hat{\mathbf M}_1^{-1}\hat{\bar{\mathbf x}}_1 = \mathbf e_1/N_1\) and \(\hat{\bar{\mathbf x}}_1^\top\hat{\mathbf M}_1^{-1}\mathbf x_i = 1/N_1\).
That the two steps return the solution of Equation 5.26 is proved in Section 5.5.4.
Theorem 5.6 (Calibration Equals Regression (Fuller 1975; Deville and Särndal 1992)) The solution of Equation 5.26 is \(\hat\omega_i = \omega_{1i}\) of Equation 5.23, and for every outcome variable \(Y\), \[\sum_{i=1}^NT_i\,\hat\omega_iY_i = \bar{\mathbf x}^\top\hat{\mathbf M}_1^{-1}\sum_{i=1}^NT_i\mathbf x_iY_i = \bar{\mathbf x}^\top\hat{\mathbf B}_1 = \hat\mu_{1,\mathrm{reg}} .\]
The identity is exact, not asymptotic. The \(N_1\)-dimensional weight problem and the \(p\)-dimensional least-squares problem compute the same number. The weights do not depend on \(Y\), and the coefficient does: calibration separates the work on \(\mathbf x\), done once, from the work on \(Y\), done for each outcome. In an experiment with many outcomes, one set of weights gives the regression estimator for each, and the design consistency of Section 5.4 carries over to the weighted form. Nothing in the problem prevents \(\hat\omega_i < 0\).
5.5.4 KKT Conditions and Strong Duality
The multiplier method is an instance of a general characterization. Write \(n = N_1\) for the size of arm 1, let \(\mathbf X_1\) be the \(n \times p\) matrix of treated covariate rows, and consider the calibration set \[\mathcal C_\omega = \{\boldsymbol\omega \in \mathbb R^n : \mathbf X_1^\top\boldsymbol\omega = \bar{\mathbf x}\}.\] Two subspaces of \(\mathbb R^n\) organize the argument: \[\mathrm{Col}(\mathbf X_1) = \{\mathbf X_1\mathbf a : \mathbf a \in \mathbb R^p\}, \qquad \mathrm{Ker}(\mathbf X_1^\top) = \{\mathbf k \in \mathbb R^n : \mathbf X_1^\top\mathbf k = \mathbf 0\}.\] \(\mathrm{Ker}(\mathbf X_1^\top)\), the null space of \(\mathbf X_1^\top\), contains the directions that leave every weighted \(\mathbf x\)-mean unchanged. \(\mathrm{Col}(\mathbf X_1)\) contains the vectors whose \(i\)th entry is a linear function \(\mathbf x_i^\top\mathbf a\) of the covariates.
Lemma 5.6 (Fundamental Theorem of Linear Algebra) \(\mathrm{Ker}(\mathbf X_1^\top) = \mathrm{Col}(\mathbf X_1)^\perp\), so \(\mathbb R^n = \mathrm{Col}(\mathbf X_1) \oplus \mathrm{Ker}(\mathbf X_1^\top)\) with the two subspaces orthogonal.
Proof. \(\mathbf k \in \mathrm{Ker}(\mathbf X_1^\top)\) if and only if \(\mathbf k^\top\mathbf X_1\mathbf a = (\mathbf X_1^\top\mathbf k)^\top\mathbf a = 0\) for every \(\mathbf a\), that is, if and only if \(\mathbf k \perp \mathrm{Col}(\mathbf X_1)\). Every subspace \(\mathcal S\) of \(\mathbb R^n\) gives \(\mathbb R^n = \mathcal S \oplus \mathcal S^\perp\). \(\square\)
Proposition 5.7 (KKT Conditions for Calibration) Let \(D : \mathbb R^n \to \mathbb R\) be convex and differentiable, and let \(\mathcal C_\omega\) be nonempty. Then \(\hat{\boldsymbol\omega}\) minimizes \(D\) over \(\mathcal C_\omega\) if and only if, for some \(\hat{\boldsymbol\lambda} \in \mathbb R^p\),
feasibility: \(\mathbf X_1^\top\hat{\boldsymbol\omega} = \bar{\mathbf x}\);
stationarity: \(\nabla D(\hat{\boldsymbol\omega}) = \mathbf X_1\hat{\boldsymbol\lambda}\), that is, \(\partial D(\hat{\boldsymbol\omega})/\partial\omega_i = \hat{\boldsymbol\lambda}^\top\mathbf x_i\) for every treated unit.
If \(D\) is strictly convex, the minimizer is unique; if \(\mathbf X_1\) has full column rank, the multiplier is unique.
Proof. Feasible directions. For \(\hat{\boldsymbol\omega} \in \mathcal C_\omega\), the line \(\hat{\boldsymbol\omega} + t\mathbf k\) stays in \(\mathcal C_\omega\) for all \(t \in \mathbb R\) if and only if \(\mathbf X_1^\top\mathbf k = \mathbf 0\). The feasible directions are exactly \(\mathrm{Ker}(\mathbf X_1^\top)\), each in both senses.
Necessity. Let \(\hat{\boldsymbol\omega}\) minimize \(D\) over \(\mathcal C_\omega\) and let \(\mathbf k \in \mathrm{Ker}(\mathbf X_1^\top)\). The function \(g(t) = D(\hat{\boldsymbol\omega} + t\mathbf k)\) is differentiable and minimized at \(t = 0\) over all of \(\mathbb R\), so \(0 = g'(0) = \nabla D(\hat{\boldsymbol\omega})^\top\mathbf k\). This holds for every \(\mathbf k \in \mathrm{Ker}(\mathbf X_1^\top)\), so \(\nabla D(\hat{\boldsymbol\omega}) \in \mathrm{Ker}(\mathbf X_1^\top)^\perp = \mathrm{Col}(\mathbf X_1)\) by Lemma 5.6.
Sufficiency. Let \(\hat{\boldsymbol\omega} \in \mathcal C_\omega\) satisfy (ii), and let \(\boldsymbol\omega \in \mathcal C_\omega\). Then \(\boldsymbol\omega - \hat{\boldsymbol\omega} \in \mathrm{Ker}(\mathbf X_1^\top)\) and \(\nabla D(\hat{\boldsymbol\omega}) \in \mathrm{Col}(\mathbf X_1)\), so \(\nabla D(\hat{\boldsymbol\omega})^\top(\boldsymbol\omega - \hat{\boldsymbol\omega}) = 0\). Convexity gives \(D(\boldsymbol\omega) \ge D(\hat{\boldsymbol\omega}) + \nabla D(\hat{\boldsymbol\omega})^\top(\boldsymbol\omega - \hat{\boldsymbol\omega}) = D(\hat{\boldsymbol\omega})\).
Uniqueness. If \(D\) is strictly convex, two minimizers would give a strictly smaller value at their midpoint, which is feasible. If \(\mathbf X_1\) has full column rank, \(\mathbf X_1\hat{\boldsymbol\lambda} = \mathbf X_1\hat{\boldsymbol\lambda}'\) implies \(\hat{\boldsymbol\lambda} = \hat{\boldsymbol\lambda}'\). \(\square\)
With equality constraints only, the KKT (Karush–Kuhn–Tucker) conditions are (i) and (ii); together they say \(\nabla\mathcal L(\hat{\boldsymbol\omega}, \hat{\boldsymbol\lambda}) = \mathbf 0\), in \(\boldsymbol\omega\) (stationarity) and in \(\boldsymbol\lambda\) (feasibility, \(\nabla_{\boldsymbol\lambda}\mathcal L = \mathbf h\)). For the Euclidean distance, \(\nabla D(\boldsymbol\omega) = \boldsymbol\omega - \boldsymbol\omega^{(0)}\) with anchor \(\boldsymbol\omega^{(0)} = n^{-1}\mathbf 1\), and stationarity is the linear weight form \(\boldsymbol\omega = \boldsymbol\omega^{(0)} + \mathbf X_1\boldsymbol\lambda\) of Section 5.5.3. For a general \(D\) it is the gradient, not the weights, that is linear in \(\mathbf x_i\); a nonlinear \(\nabla D\) gives nonlinear weights.
Proposition 5.8 (Strong Duality) Assume \(D\) is convex and differentiable on \(\mathbb R^n\), the constraints are affine, and the primal problem has a minimizer \(\hat{\boldsymbol\omega}\). With \(\ell(\boldsymbol\lambda) = \inf_{\boldsymbol\omega}\mathcal L(\boldsymbol\omega, \boldsymbol\lambda)\), \[\max_{\boldsymbol\lambda}\ell(\boldsymbol\lambda) = \ell(\hat{\boldsymbol\lambda}) = D(\hat{\boldsymbol\omega}) = \min_{\boldsymbol\omega}\sup_{\boldsymbol\lambda}\mathcal L(\boldsymbol\omega, \boldsymbol\lambda),\] with \(\hat{\boldsymbol\lambda}\) the multiplier of Proposition 5.7.
Proof. By Proposition 5.7, \(\nabla D(\hat{\boldsymbol\omega}) = \mathbf X_1\hat{\boldsymbol\lambda}\), so \(\nabla_{\boldsymbol\omega}\mathcal L(\hat{\boldsymbol\omega}, \hat{\boldsymbol\lambda}) = \mathbf 0\). The function \(\mathcal L(\cdot, \hat{\boldsymbol\lambda}) = D + \hat{\boldsymbol\lambda}^\top\mathbf h\) is convex, so its stationary point \(\hat{\boldsymbol\omega}\) is its global minimizer: step 1 of the multiplier method returns \(\hat{\boldsymbol\omega}\) at \(\hat{\boldsymbol\lambda}\). Since \(\mathbf h(\hat{\boldsymbol\omega}) = \mathbf 0\), \(\ell(\hat{\boldsymbol\lambda}) = \mathcal L(\hat{\boldsymbol\omega}, \hat{\boldsymbol\lambda}) = D(\hat{\boldsymbol\omega})\). For any \(\boldsymbol\lambda\), weak duality gives \(\ell(\boldsymbol\lambda) \le \mathcal L(\hat{\boldsymbol\omega}, \boldsymbol\lambda) = D(\hat{\boldsymbol\omega})\), so \(\hat{\boldsymbol\lambda}\) maximizes \(\ell\). Finally, \(\sup_{\boldsymbol\lambda}\mathcal L(\boldsymbol\omega, \boldsymbol\lambda)\) equals \(D(\boldsymbol\omega)\) on \(\mathcal C_\omega\) and \(+\infty\) off it, so its minimum is \(D(\hat{\boldsymbol\omega})\). \(\square\)
Strong duality justifies steps 2 and 3 of the multiplier method: maximizing the \(p\)-dimensional dual returns the multiplier of the primal solution, and the linear weight form at that multiplier is the calibrated weight vector.
5.5.5 Geometry: Orthogonal Projection
For the Euclidean distance the solution has a geometric reading. \(\mathcal C_\omega\) is an affine set, \(\mathcal C_\omega = \boldsymbol\omega_1 + \mathrm{Ker}(\mathbf X_1^\top)\) for any \(\boldsymbol\omega_1 \in \mathcal C_\omega\), and the level sets of \(D(\boldsymbol\omega) = \tfrac12\|\boldsymbol\omega - \boldsymbol\omega^{(0)}\|^2\) are spheres around the anchor. Growing the spheres until one touches \(\mathcal C_\omega\) gives the solution. At the touching point, moving along \(\mathcal C_\omega\) cannot lower \(D\) to first order, so \(\hat{\boldsymbol\omega} - \boldsymbol\omega^{(0)} \perp \mathrm{Ker}(\mathbf X_1^\top)\), that is, \(\hat{\boldsymbol\omega} - \boldsymbol\omega^{(0)} \in \mathrm{Col}(\mathbf X_1)\): the linear weight form, read as a picture.
Proposition 5.9 (Pythagorean Identity) For every \(\boldsymbol\omega_1 \in \mathcal C_\omega\), \[D(\boldsymbol\omega_1) = D(\hat{\boldsymbol\omega}) + \tfrac12\|\boldsymbol\omega_1 - \hat{\boldsymbol\omega}\|^2 \;\ge\; D(\hat{\boldsymbol\omega}),\] with equality if and only if \(\boldsymbol\omega_1 = \hat{\boldsymbol\omega}\). Hence the \(L_2\) calibration problem has a unique solution.
Proof. \(\boldsymbol\omega_1 - \boldsymbol\omega^{(0)} = (\boldsymbol\omega_1 - \hat{\boldsymbol\omega}) + (\hat{\boldsymbol\omega} - \boldsymbol\omega^{(0)})\). The first vector lies in \(\mathrm{Ker}(\mathbf X_1^\top)\) because both \(\boldsymbol\omega_1\) and \(\hat{\boldsymbol\omega}\) are calibrated; the second lies in \(\mathrm{Col}(\mathbf X_1)\) by the linear weight form. By Lemma 5.6 they are orthogonal, so the cross term vanishes. \(\square\)
The argument proves optimality without multipliers. It is the same right-triangle argument as the analysis-of-variance decomposition in regression.
The two-projection formula. Let \(\mathbf P_1 = \mathbf X_1(\mathbf X_1^\top\mathbf X_1)^{-1}\mathbf X_1^\top\), the hat matrix of the regression on \(\mathbf X_1\) in arm 1. For any \(\boldsymbol\omega_1 \in \mathcal C_\omega\), \(\mathbf X_1^\top(\boldsymbol\omega_1 - \boldsymbol\omega^{(0)}) = \bar{\mathbf x} - \mathbf X_1^\top\boldsymbol\omega^{(0)} = \mathbf X_1^\top(\hat{\boldsymbol\omega} - \boldsymbol\omega^{(0)})\), and \(\hat{\boldsymbol\omega} - \boldsymbol\omega^{(0)} = \mathbf X_1\hat{\boldsymbol\lambda} \in \mathrm{Col}(\mathbf X_1)\), so \(\hat{\boldsymbol\omega} - \boldsymbol\omega^{(0)} = \mathbf P_1(\boldsymbol\omega_1 - \boldsymbol\omega^{(0)})\). Rearranging, \[\hat{\boldsymbol\omega} = (\mathbf I - \mathbf P_1)\boldsymbol\omega^{(0)} + \mathbf P_1\boldsymbol\omega_1 . \tag{5.28}\] The first term is the projection of the anchor onto \(\mathrm{Ker}(\mathbf X_1^\top)\); the second is the projection of any calibrated vector onto \(\mathrm{Col}(\mathbf X_1)\). The anchor matters only through its residual. Here \(\mathbf x_i\) contains an intercept, so the uniform anchor \(\boldsymbol\omega^{(0)} = N_1^{-1}\mathbf 1\) lies in \(\mathrm{Col}(\mathbf X_1)\), the first term vanishes, and \(\hat{\boldsymbol\omega} = \mathbf P_1\boldsymbol\omega_1\) for every calibrated \(\boldsymbol\omega_1\). Since \(\hat{\boldsymbol\omega} \in \mathrm{Col}(\mathbf X_1)\) and every calibrated vector differs from it by an element of \(\mathrm{Ker}(\mathbf X_1^\top)\), \(\|\boldsymbol\omega_1\|^2 = \|\hat{\boldsymbol\omega}\|^2 + \|\boldsymbol\omega_1 - \hat{\boldsymbol\omega}\|^2\): the regression weights are the minimum-norm calibrated weights, which links the geometry back to the variance criterion of Section 5.5.2. In survey sampling the anchor is a vector \(\mathbf d\) of design weights, which is usually not in \(\mathrm{Col}(\mathbf X_1)\). For the unscaled Euclidean distance around \(\mathbf d\), the calibrated weights are \(\hat{\boldsymbol\omega} = \mathbf d + \mathbf X_1(\mathbf X_1^\top\mathbf X_1)^{-1}(\bar{\mathbf x} - \mathbf X_1^\top\mathbf d)\), so \[\hat{\boldsymbol\omega}^\top\mathbf y = \bar{\mathbf x}^\top\hat{\mathbf B}_1 + \mathbf d^\top(\mathbf y - \mathbf X_1\hat{\mathbf B}_1).\] With a general design-weight anchor, the calibrated estimator therefore retains an anchor-weighted residual term. It has the generalized regression form and need not equal the pure OLS prediction \(\bar{\mathbf x}^\top\hat{\mathbf B}_1\).
5.6 Stratification
Stratification (subclassification) divides the units into strata of similar covariate values and compares treated and control units within each stratum. It is the model-free implementation of the back-door formula, and it is also a special case of calibration.
5.6.1 Calibration on Stratum Indicators
Let the strata be \(\mathcal S_k = \{i : X_i \in B_k\}\), \(k = 1, \ldots, K\), with \(N_k\) units in stratum \(k\) and \(N_{tk}\) units of arm \(t\) in stratum \(k\). Take the covariate vector of Section 5.5 to be the stratum indicators, \(\mathbf x_i = (x_{i1}, \ldots, x_{iK})^\top\) with \(x_{ik} = \mathbb{I}(i \in \mathcal S_k)\). This is a full-rank coding with no extra intercept; the indicators sum to one, so the intercept is in their span.
Proposition 5.10 (Stratification Is Calibration on Stratum Indicators) Assume \(N_{tk} \ge 1\) for every \(t\) and \(k\). With stratum indicators as covariates:
- the calibration constraint Equation 5.24 for arm \(t\) says \(\sum_{i:\,T_i=t,\,i\in\mathcal S_k}\omega_{ti} = N_k/N\) for every \(k\): each stratum receives its population share;
- the \(L_2\)-calibrated weights are equal within each arm-by-stratum cell, \(\hat\omega_{ti} = N_k/(N\,N_{tk})\) for \(i \in \mathcal S_k\);
- the resulting regression estimators are \[\hat\mu_{t,\mathrm{reg}} = \sum_{k=1}^K\frac{N_k}{N}\,\bar Y_{t,k}, \qquad \hat\mu_{1,\mathrm{reg}} - \hat\mu_{0,\mathrm{reg}} = \hat\tau_{\mathrm{STRAT}},\] where \(\bar Y_{t,k}\) is the mean of \(Y\) among arm-\(t\) units in stratum \(k\) and \(\hat\tau_{\mathrm{STRAT}}\) is defined in Equation 5.29 below.
Proof. (i) The \(k\)th coordinate of \(\sum_iT_i\omega_{1i}\mathbf x_i\) is the total weight of the treated units in stratum \(k\), and the \(k\)th coordinate of \(\bar{\mathbf x}\) is \(N_k/N\). (ii) \(\hat{\mathbf M}_t = \mathrm{diag}(N_{t1}, \ldots, N_{tK})\), which is invertible because every \(N_{tk} \ge 1\), so by Equation 5.27 \(\hat\omega_{ti} = \bar{\mathbf x}^\top\hat{\mathbf M}_t^{-1}\mathbf x_i = N_k/(NN_{tk})\) for \(i \in \mathcal S_k\). (iii) Summing \(\hat\omega_{ti}Y_i\) over the arm-\(t\) units of stratum \(k\) gives \((N_k/N)\bar Y_{t,k}\). \(\square\)
The full-rank condition \(N_{tk} \ge 1\) is the within-stratum overlap requirement: a stratum with no treated or no control units cannot be calibrated. In a completely randomized experiment the same estimator is the post-stratified difference in means, which is regression on stratum indicators. For a fixed number of strata whose proportions satisfy the regularity conditions of Section 5.4, it is design consistent and has no larger first-order asymptotic variance than \(\hat\tau_{\mathrm{DIM}}\). The finite-sample condition \(N_{tk} \ge 1\) alone is not such a condition, and progressively finer strata, with a growing number of indicators, fall outside that result.
5.6.2 The Stratified Estimator and Its Conditions
Within stratum \(k\), the difference in means is \(\hat\tau_k = \bar Y_{1,k} - \bar Y_{0,k}\), and the stratified estimator averages the within-stratum contrasts with the population shares, \[\hat\tau_{\mathrm{STRAT}} = \sum_{k=1}^K\hat\tau_k\,\frac{N_k}{N}. \tag{5.29}\] For categorical \(X\), with one stratum per category, it is the standardization estimator Equation 5.10. For continuous \(X\), the strata coarsen \(X\), and coarsening does not by itself inherit exchangeability: residual confounding can persist within every finite-width stratum. What justifies the method in an observational study is a combination of conditions:
- exchangeability given the exact \(X\);
- treatment propensity and outcome regressions that vary smoothly within strata;
- stratum widths that shrink as \(N\) grows;
- adequate within-stratum overlap, with growing numbers of treated and control units within strata.
Under these conditions the coarsening bias of \(\hat\tau_k\) vanishes, and \(\hat\tau_k\) is consistent for the stratum-specific effect. With finite widths, coarsening can leave residual confounding within a stratum, so finite-width strata do not generally give unbiased comparisons. Exact categorical strata, or strata within which treatment is randomized, are the exceptions, provided both arms are represented. The last condition matters on its own: a stratum with one unit from each arm has overlap but highly variable cell means.
Cochran (1968) studied one continuous confounder under particular distributional models. With \(K = 5\) strata of equal size, about 90% of the bias is removed; with \(K = 10\), over 95%. This is a benchmark, not a universal fact. More strata are not always better: with too many, empty cells and within-stratum variance dominate, and the mean squared error grows.
5.6.3 Three Implementations of One Formula
Regression adjustment, standardization, and stratification implement the same identification formula Equation 5.7 and rest on the same identification conditions: consistency, conditional exchangeability, and positivity. They differ in how they handle \(X\) at the estimation stage (Table 5.8).
| Method | How it handles \(X\) | Estimation-stage requirement |
|---|---|---|
| Regression / standardization | Estimates \(\mu(t, x)\): parametrically, nonparametrically, or by cells | Sufficient convergence of \(\hat\mu\); correct parametric specification is one route |
| Stratification | Coarsens \(X\) into cells | Vanishing coarsening bias, within-cell overlap, growing treated and control counts within strata |
| Propensity score (Chapter 6) | Reduces \(X\) to the scalar \(\pi(X) = P(T = 1 \mid X)\) | Consistent score estimation or estimator-specific balance conditions, plus overlap |
All three face the curse of dimensionality. With \(p\) binary covariates there are \(2^p\) cells, most of them sparse or empty. Regression interpolates across sparse cells at the price of a model. The propensity score reduces the dimension needed for adjustment to one (Chapter 6), although estimating it from a high-dimensional \(X\) can still be difficult. None of the three repairs unmeasured confounding.
5.7 Lab: Simulation Study of the Outcome Regression Estimator
The identification results of Section 5.2 determine the target before any outcome model is chosen. Estimating it by the regression estimator of Section 5.3 then creates a bias–variance tradeoff: a restrictive working model may be biased, while a flexible estimator pays in variance. Under complete randomization, Theorem 5.4 gives the arm-wise linear regression estimator a design-consistency guarantee even when its working model is wrong; the guarantee disappears when treatment is confounded. This simulation puts the two claims side by side.
Estimators. Both estimators follow the three-step recipe of Section 5.3.3: estimate \(\hat\mu(t, x)\) within each arm, predict both conditional means at every observed \(X_i\), and average, \(\hat\tau_{\mathrm{reg}} = N^{-1}\sum_i\{\hat\mu(1, X_i) - \hat\mu(0, X_i)\}\). They differ only in how \(\hat\mu\) is estimated.
- Linear OR. Within each arm, fit \(Y\) on \((1, X)\) by OLS. The model is misspecified: the true \(\mu(t, x) = t + 3(1 + t)x^2\) has a quadratic term and a treatment–covariate interaction.
- Local-linear OR. Within each arm, local linear regression with a Gaussian kernel and the prescribed bandwidth rule \(h = 0.5\,N^{-1/5}\) (
KernSmooth::locpolyin R), computed on a 401-point grid covering the whole target range \([0, 1]\); predictions at the observed \(X_i\) are obtained by linear interpolation of the fitted curve (approx). Local linear regression removes the leading boundary bias of the Nadaraya–Watson estimator, which matters here because the observational design places treated and control units in different regions of \([0, 1]\).
2024. The Monte Carlo standard errors of the reported means are \(0.0016\), \(0.0015\), \(0.0020\), and \(0.0021\) (rows in order). Generated by lab5_sim.R (R 4.3.3, KernSmooth 2.23-22), in which each replication draws one population \((X, \varepsilon)\) that both designs share.
| Design | Estimator | Mean | Bias | SD | RMSE |
|---|---|---|---|---|---|
| CRD | Linear OR | 1.9964 | \(-0.0036\) | 0.0709 | 0.0710 |
| Local-linear OR | 2.0260 | \(+0.0260\) | 0.0685 | 0.0732 | |
| Observational | Linear OR | 1.9674 | \(-0.0326\) | 0.0895 | 0.0952 |
| Local-linear OR | 2.0273 | \(+0.0273\) | 0.0951 | 0.0989 |
Three lessons stand out (Table 5.9).
1. Under the CRD, the misspecified linear estimator is consistent. Conditional on each simulated finite population, Theorem 5.4 gives \(\hat\tau_{\mathrm{reg}} - \bar\tau_N \pto 0\) under the randomization distribution. Because the units are first drawn i.i.d. from the superpopulation, \(\bar\tau_N \pto 2\), and the two phases together give \(\hat\tau_{\mathrm{reg}} \pto \tau_{\mathrm{ATE}}\). Randomization, not the model, supplies the guarantee.
2. Under the observational design, the same estimator is biased. The estimated absolute bias increases from \(0.0036\) to \(0.0326\). The model and the data-generating process are the same; only the assignment mechanism differs. The two estimates carry different weight. Under the CRD, \(-0.0036\) is about two Monte Carlo standard errors from zero and is consistent with the \(O(N^{-1})\) finite-sample behavior described above. Under the observational design, numerical integration of the population arm-wise fits gives a probability limit of about \(1.9652\), a limiting bias of about \(-0.0348\), which the simulation estimate matches. A single table at \(N = 1{,}000\) illustrates finite-sample performance; consistency follows from the theory, not from one sample size.
3. Flexibility trades bias for variance. Under the observational design the local-linear estimator reduces the absolute bias from \(0.0326\) to \(0.0273\), but its standard deviation is \(0.0951\) against \(0.0895\), and its RMSE is larger (\(0.0989\) against \(0.0952\)). In this finite-sample comparison the parsimonious, misspecified model has the smaller mean squared error.
5.8 Where the Propensity Score Fits
Section 5.6.3 placed regression, standardization, and stratification as three implementations of one formula, and the calibration view of Section 5.5 showed what they share: each reweights the arms so that they look alike on the covariates. Regression does it through calibrated weights, stratification through cells, and the propensity score through a model for treatment. When \(X\) has many components, the propensity score provides an identification-preserving dimension reduction: adjusting for the scalar \(\pi(X)\) suffices for identification even when \(X\) is high-dimensional. The theorem reduces the dimension of the balancing variable, not every statistical difficulty.
5.9 Summary
Randomization. A known assignment mechanism supplies both identification and the probability distribution for inference. Under complete randomization the difference in means is exactly unbiased for the finite-population ATE \(\bar\tau_N\), its design variance is given by Neyman’s formula, and the Neyman variance estimator is conservative because the variance of the unit-level effects is not identified. Fisher’s randomization test gives an exact test of the sharp null. Extending conclusions to a superpopulation needs a sampling or transportability argument.
Ignorability. In an observational study, consistency, weak conditional exchangeability, and positivity identify the ATE and the ATT through the back-door formula. The assumption is substantive, cannot be verified from the observed law, and requires pretreatment adjustment variables.
Outcome regression. The regression estimator \(\hat\tau_{\mathrm{reg}} = \hat\mu_{1,\mathrm{reg}} - \hat\mu_{0,\mathrm{reg}}\) plugs a fitted outcome model and the empirical covariate distribution into the identified functional. Under a correct linear model its error splits into an estimation error, from the nuisance parameter \(\hat\beta\), and a sampling error, from averaging over the sampled covariates. Their first-order influence-function components are uncorrelated, and the influence function gives the asymptotic variance.
Regression under randomization. Under complete randomization the same estimator is asymptotically equivalent to an oracle difference estimator, an SRS mean of fixed residuals. Exact results hold for the oracle (linearized) estimator: its variance is Neyman’s formula with residuals, and it is never less precise than the difference in means. To first order the fitted estimator inherits these properties: it is design consistent without any outcome model and cannot increase the asymptotic variance. The remainder of the linearization is negligible because randomization balances \(\mathbf x\). In an observational study it can be first order, and correct specification of the outcome model is a sufficient, not a necessary, condition for consistency.
Calibration. The arm-wise linear regression estimator studied here, least squares with an intercept, is a weighted mean whose weights reproduce the full-sample covariate mean exactly and, among all such weights, are the closest to uniform. For this estimator, Euclidean calibration equals regression for every outcome; balancing the included covariate means does not by itself establish exchangeability; the KKT conditions characterize optimal calibration weights for any convex distance, and for the Euclidean distance the solution is an orthogonal projection.
Stratification and the propensity score. Stratification is calibration on stratum indicators. Regression, standardization, and stratification implement one identification formula and differ only in how they handle \(X\); the propensity score of Chapter 6 reduces \(X\) to a scalar. None of them repairs unmeasured confounding.
5.10 Problems
1. Randomization and the do-calculus. Let the observational DAG be \(\{U \to T,\; U \to Y,\; X \to T,\; X \to Y,\; T \to Y\}\) with \(U\) unobserved.
- Identify all back-door paths from \(T\) to \(Y\).
- Does \(X\) satisfy the back-door criterion? Does \(U\)? Explain.
- Now suppose treatment is randomized. Draw the randomized-assignment graph and identify all back-door paths. Show that \(p_R(y \mid T = t) = f(y \mid \doop(T = t))\) holds within the randomized law using Rule 2 of the do-calculus. (Hint: construct \(G^R_{\underline{T}}\) from the randomized-assignment graph \(G^R\) and verify \((Y \indep T)_{G^R_{\underline{T}}}\).)
- Under randomization, is covariate adjustment on \(X\) necessary for unbiasedness? Is it ever beneficial? Relate your answer to Section 5.4.
2. Finite-population enumeration. Consider a finite population of \(N = 4\) units with potential-outcome schedule
| Unit \(i\) | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| \(Y_i(0)\) | 0 | 1 | 1 | 2 |
| \(Y_i(1)\) | 1 | 3 | 2 | 4 |
and complete randomization with \(N_1 = N_0 = 2\).
- Compute \(\bar\tau_N\), \(S_1^2\), \(S_0^2\), \(S_{01}\), and \(S_\tau^2\).
- Enumerate the six equally likely assignments and compute \(\hat\tau_{\mathrm{DIM}}\) for each. Verify that \(\E(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_N) = \bar\tau_N\).
- Compute the exact design variance from your enumeration and check it against Equation 5.3.
3. Conservative Neyman variance. Continue with the schedule of the previous problem.
- For each assignment, compute \(\hat V\) of Equation 5.4, and average to obtain \(\E(\hat V \mid \mathcal F_N)\).
- Verify that \(\E(\hat V \mid \mathcal F_N) - \mathrm{Var}(\hat\tau_{\mathrm{DIM}} \mid \mathcal F_N) = S_\tau^2/N\).
- Keeping the \(Y_i(0)\) column, change the \(Y_i(1)\) column so that the gap is exactly zero. Explain in terms of treatment-effect heterogeneity.
4. A Fisher randomization test. Continue with the schedule of Problem 2, and suppose units 2 and 4 are treated.
- Write down the observed outcome vector and compute \(\hat\tau_{\mathrm{DIM}}\).
- Under the sharp null \(Y_i(1) = Y_i(0)\) for every \(i\), compute the difference in means for all six assignments, holding the observed outcomes fixed, and report the one-sided and two-sided randomization \(p\)-values.
- The sharp null is false for this schedule. Explain why the test cannot reject it at the 5% level, and how large \(N\) must be, with \(N_1 = N_0\), before a \(p\)-value strictly below \(0.05\) is attainable, for the one-sided and for the two-sided test.
5. Ignorability and the back-door criterion. Consider the DAG \(\{X \to T,\; X \to Y,\; T \to Y,\; U \to Y\}\) with \(U\) unobserved.
- Does \(\{X\}\) satisfy the back-door criterion? Write the identifying formula for \(\tau_{\mathrm{ATE}}\).
- Add the edge \(U \to T\). Does \(\{X\}\) still satisfy the criterion? What would an adjustment set have to contain?
- In words, what does unconfoundedness given \(X\) require of \(U\) in the data-generating process? Does an unobserved variable that predicts \(Y\) violate it by itself?
6. Standardization. A study of job training (\(T\)) and earnings (\(Y\), in thousands) produces the following cell means:
| \(X\) | \(\hat\mu(1, x)\) | \(\hat\mu(0, x)\) | \(N_x/N\) |
|---|---|---|---|
| 0 | 28 | 22 | 0.4 |
| 1 | 35 | 31 | 0.6 |
- Compute \(\hat\tau_{\mathrm{STD}}\).
- Given \(P(X = 1 \mid T = 1) = 0.8\) and \(P(X = 1 \mid T = 0) = 0.4\), compute the ATT estimate.
- The unadjusted difference in means is \(33.6 - 25.6 = 8\). Compare it with your answers to (a) and (b) and explain the discrepancy.
7. Simpson’s paradox. A hospital reports that patients admitted to intensive care (\(T = 1\)) have a higher mortality rate than those who are not: \(P(Y = 1 \mid T = 1) = 0.30\) and \(P(Y = 1 \mid T = 0) = 0.10\), with 1,000 patients in each group.
- Construct a \(2 \times 2 \times 2\) table, with disease severity \(X \in \{\text{mild}, \text{severe}\}\), that is consistent with the pooled numbers and in which observed mortality is lower among intensive-care patients within both strata.
- Compute the standardized contrast and the contrast standardized to the treated distribution. Under which assumptions (consistency, positivity, and exchangeability given severity) are these the ATE and the ATT?
- Draw the DAG and identify the back-door path that the pooled comparison fails to block.
- A colleague argues that the pooled statistic is the right answer because the hospital treats both kinds of patients. Explain, using potential outcomes, why this is incorrect.
8. The influence function of the regression estimator. Consider model (b) with a scalar covariate, \(\mu(t, x) = \alpha_t + \gamma_tx\), and i.i.d. data satisfying the conditions of Theorem 5.3.
- Show that the sampling component of Equation 5.14 is \((\gamma_1 - \gamma_0)^2\,\mathrm{Var}(X)\). When is it zero?
- Write \(c = Pz_d\) and \(A\) explicitly, and show that the estimation component is the asymptotic variance of \((\hat\alpha_1 - \hat\alpha_0) + (\hat\gamma_1 - \hat\gamma_0)\E(X)\), with \(\E(X)\) treated as known. Here “asymptotic variance” means the variance of the limiting distribution of \(\sqrt N\) times the estimation error; the estimator’s leading variance divides it by \(N\).
- Explain why the sandwich standard error of the coefficient on \(T\) in the centered interacted regression of Proposition 5.3 estimates only the second component.
- Under model (a), show that the sampling component is zero and that Equation 5.15 reduces to the sandwich variance of the coefficient on \(T\).
9. Three forms of the regression estimator. Consider the fully interacted regression of \(Y_i\) on \((1, T_i, X_i - \bar X, T_i(X_i - \bar X))\), with \(\bar X\) the full-sample mean.
- Show that this regression is equivalent to separate OLS fits of \(Y\) on \((1, X)\) within each arm, and that its coefficient on \(T_i\) is \(\hat\mu_{1,\mathrm{reg}} - \hat\mu_{0,\mathrm{reg}}\).
- What does the coefficient on \(T_i\) estimate if the covariate is not centered?
- Verify the three forms of Lemma 5.3 numerically for a data set of your choice, and show that the intercept guarantees that the expansion estimator of the residual mean is zero. (Without an intercept it can still be zero by accident.)
10. Exact unbiasedness versus design consistency. Under complete randomization, \(\hat\tau_{\mathrm{DIM}}\) is exactly unbiased for \(\bar\tau_N\) (Theorem 5.1), while \(\hat\tau_{\mathrm{reg}}\) is only design consistent (Theorem 5.4).
- Identify the step in the proof of Theorem 5.1 that gives exact unbiasedness, and explain why complete randomization does not generally guarantee exact unbiasedness of the fitted regression estimator \(\hat\tau_{\mathrm{reg}}\).
- Locate the \(O_p(N^{-1})\) remainder in the proof of Theorem 5.4, and explain why asymptotic negligibility rather than exact unbiasedness is the relevant conclusion.
- A colleague argues that because \(\hat\tau_{\mathrm{reg}}\) can be biased in finite samples, the difference in means should always be preferred. Evaluate this argument using the finite-sample-bias remark of Section 5.4.5 and Theorem 5.5.
11. The remainder in an observational study. Let \(X \sim \mathrm{Uniform}(0, 1)\), let \(\varepsilon \sim \mathcal N(0, 1)\) be independent of \(X\), and generate \(T \mid X \sim \mathrm{Bern}\{\expit(-2 + 5X)\}\) independently of \(\varepsilon\) given \(X\). Set \(Y(1) = 3X^2 + \varepsilon\). From an i.i.d. sample, fit \(Y\) on \((1, X)\) among the treated units.
- Find \(\mathbf B_1^\ast = \operatorname{plim}\hat{\mathbf B}_1\) as an expression in the moments of \(X\) given \(T = 1\), and compare it with the population projection of \(Y(1)\) on \((1, X)\).
- Show that \(\bar X - \hat{\bar X}_1\) converges to a nonzero constant, and conclude that the remainder of the observational expansion in Section 5.4.3 is of order \(N^{-1/2}\).
- Does the leading term converge to \(\E\{Y(1)\}\)? What changes if \(Y(1) = 1 + 3X + \varepsilon\)?
12. Calibration with an entropy distance. Replace the Euclidean distance in Equation 5.26 by \(D(\boldsymbol\omega) = \sum_{i:\,T_i=1}\{\omega_i\log(\omega_i/\omega_i^{(0)}) - \omega_i + \omega_i^{(0)}\}\), defined for positive weights, with \(\omega_i^{(0)} = 1/N_1\).
- Show that \(D\) is strictly convex on the positive orthant, and that the stationarity condition of Proposition 5.7 gives the exponential-tilting weights \(\omega_i = \omega_i^{(0)}\exp(\boldsymbol\lambda^\top\mathbf x_i)\).
- Proposition 5.7 assumes \(D\) is defined on all of \(\mathbb R^n\). Explain which step of its proof still holds when the minimizer lies in the interior of the positive orthant.
- With one covariate taking the values 0 and 1 in the treated arm, show that no positive weights reproduce a target mean of 2, while the Euclidean problem has a solution. State the general condition for a positive calibrated solution to exist.
13. Stratification as calibration. A binary confounder \(X\) defines two strata. Stratum \(X = 0\) has 200 treated units with mean outcome 10 and 400 controls with mean 8; stratum \(X = 1\) has 300 treated units with mean 15 and 100 controls with mean 12.
- Compute \(\hat\tau_{\mathrm{STRAT}}\) and the unadjusted difference in means, and explain why they differ.
- Compute the \(L_2\)-calibrated weights of Proposition 5.10 for every arm-by-stratum cell, and verify that they reproduce the stratum shares and give \(\hat\tau_{\mathrm{STRAT}}\).
- Explain how forming strata with 15 continuous covariates can produce sparse or empty arm-by-stratum cells.