This is the second of two lectures on the linear model. It builds on the simple linear model, OLS, unbiasedness and the zero conditional mean assumption, standard errors and the \(t\)-test.
The first half covers estimation and testing: multiple regression, control variables, \(t\)- and \(F\)-tests, dummy variables and interactions, and reading statistical output. It produces the estimate \(\hat{\beta}_j\).
The second half covers inference when the classical assumptions fail: heteroskedasticity-robust and cluster-robust standard errors. It determines whether the standard error printed next to that estimate is correct.
The corresponding reading is Wooldridge, Chapters 3 and 4, plus Chapter 8 for heteroskedasticity. For clustering, Wooldridge Chapter 14 or Angrist and Pischke, Mostly Harmless Econometrics, Chapter 8.
What should you be able to do?
By the end of this lecture, you should be able to:
Interpret a coefficient in a multiple regression, and say what “holding constant” means mechanically.
Decide whether a candidate control belongs in the regression, and why.
Carry out and interpret \(t\)- and \(F\)-tests, and say how the two relate.
Interpret dummy variables, interaction terms, and sets of category dummies.
Explain what heteroskedasticity-robust standard errors correct, and what they do not.
Decide when to cluster, at what level, and predict how much difference it will make.
Guiding questions
We return to the simulated world of Lecture 2, where the true return to education is 1.2 and where leaving ability out of the regression produced an estimate of roughly 1.85.
First: which controls remove that bias, and which controls increase it?
Second: is the standard error reported next to the estimate correct?
Prediction
Write down by what factor you expect a default standard error to be wrong when the errors are correlated within groups.
We measure that factor by simulation later in the lecture and compare it with your prediction.
Controlling for what we can measure
More than one regressor
The simple regression had no way of holding other factors constant, so the coefficient on \(x\) absorbed every unobserved factor correlated with it. Multiple regression includes those factors as additional regressors.
Definition: the multiple linear regression model
\[
y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_k x_k + u
\]
\(\beta_j\) is the effect of a one-unit change in \(x_j\) on \(y\), holding all the other explanatory variables constant.
The final clause carries the interpretation: the coefficient isolates the contribution of one variable while holding the other included variables fixed, which approximates ceteris paribus without an experiment.
Estimating the model
The estimation principle is unchanged: choose \(\hat{\beta}_0, \dots, \hat{\beta}_k\) to minimize the sum of squared residuals. Writing the solution in closed form requires matrix algebra, but the software requires only one extra term in the formula.
R: feols(y ~ x1 + x2, data = df)
Python: pf.feols("y ~ x1 + x2", data = df) after import pyfixest as pf
The interpretation of the output is the subject of the rest of this lecture. Which variables belong on the right-hand side is the subject of the next three slides.
A good control is a confounder
Our goal is usually the causal effect of one variable \(x_1\). The other regressors exist to hold constant the factors that would otherwise contaminate it.
A confounder is a pre-treatment variable correlated with both \(x_1\) and \(y\). Leave it out and \(\hat{\beta}_1\) absorbs its effect, which is Lecture 2’s omitted variable bias.
Example
In estimating the effect of education on wages, innate ability is the classic confounder: it raises schooling and it raises earnings.
Controlling for a proxy such as a test score stops the return to education from being biased upwards.
The validity of a causal estimate depends on which controls are chosen. Including the wrong ones can produce more bias than including none at all, which is the subject of the next slide.
A bad control is a mediator or a collider
Two kinds of variable induce bias when you control for them.
A mediator lies on the causal pathway from \(x_1\) to \(y\). Controlling for it closes one of the channels through which \(x_1\) acts, so part of the effect being estimated is held constant.
A collider is caused by both \(x_1\) and \(y\). Controlling for it creates an association between them even when none exists.
Example: a collider
Does a degree (\(x_1\)) raise job performance (\(y\))? A firm hires candidates on the basis of both, so being hired (\(z\)) is caused by the degree and by performance.
Among the hired, the two can appear negatively related: a candidate without a degree needed high performance to be hired, while a candidate with a degree was hired at average performance.
Restricting the sample to those who were hired conditions on hiring. Including hiring as a control also conditions on this collider and can introduce bias.
Collider bias in a scatter plot
In this simulation, education and potential performance are independent. Hiring requires their standardized scores to sum to more than 1.
Among hires, less education must be offset by higher performance. Conditioning on hiring creates the negative association.
The three cases as causal diagrams
Arrows point from cause to effect. In each diagram the question is whether the third variable should be included in the regression.
The rule of thumb: control for causes of the treatment, and not for consequences of it.
Adding controls one at a time
The simulated world of Lecture 2 set \(\beta_0 = 1.5\) and \(\beta_1 = 1.2\). We keep that return to education and add two variables: experience, and innate ability, which raises both schooling and wages.
One extra year of education therefore still raises wages by 1.2. The table reports the coefficient on educ as controls are added.
Column (1) gives 1.854, well above the true 1.2. Ability is omitted, so the education coefficient absorbs the wage premium that belongs to ability.
Column (2) adds experience, and the coefficient moves only from 1.854 to 1.849. Experience is a valid regressor, but it is uncorrelated with education, so omitting it caused no bias.
Column (3) adds ability, and the coefficient falls to 1.212, close to the true value.
Selecting controls
Bias is removed by the control that is correlated with both the treatment and the outcome. Adding regressors unrelated to the variable of interest raises the \(R^2\) slightly and leaves the estimate unchanged, while adding a mediator or a collider increases the bias.
The bias in column (1) is upward because ability raises both schooling and pay. Its sign is the product of those two signs, as Lecture 2’s formula states.
The variance of a coefficient
Adding regressors also changes the variance of the estimate, in either direction.
Variance of the OLS estimator with several regressors
\(SST_j = \sum_i (x_{ij} - \bar{x}_j)^2\) is the total variation in \(x_j\).
\(R_j^2\) is the \(R^2\) from regressing \(x_j\) on all the other regressors.
The term \(R_j^2\) measures multicollinearity: the share of the variation in \(x_j\) that the other regressors already explain.
Multicollinearity and the cost of a control
As \(R_j^2\) approaches 1, the denominator approaches zero and the variance rises without bound. The no perfect collinearity assumption states that \(R_j^2 \neq 1\).
This creates a trade-off. A control that is strongly correlated with the treatment is often the control needed to remove bias, and it is also the control that costs the most precision.
The single \(\hat{\sigma}^2\) in the numerator assumes that the error variance is constant and that the errors are uncorrelated across observations. The second half of this lecture examines what happens when each of those assumptions fails.
Testing
Is the number distinguishable from zero?
Our estimates come from a sample, so they are subject to sampling variability. The true population relationship may be zero even when the estimated \(\hat{\beta}_1\) is not.
Two tests address this, and they answer different questions:
The \(t\)-test asks about one coefficient at a time.
The \(F\)-test asks about several coefficients jointly, or about the model as a whole.
The \(t\)-test
The \(t\)-test handles a hypothesis about a single coefficient, most often that the variable has no effect: \(H_0: \beta_j = 0\) against \(H_A: \beta_j \neq 0\).
where \(se(\hat{\beta}_j)\) is the square root of the variance on the previous slide.
It measures how many standard errors the estimate lies from the hypothesized value. The statistic is the one from Lecture 2, computed with the multiple-regression standard error.
The \(F\)-test
The \(F\)-test handles multiple coefficients at once. Its most common use is testing the overall significance of a model.
\(H_0\): none of the explanatory variables affect \(y\), so \(\beta_1 = \beta_2 = \dots = \beta_k = 0\). The model has no explanatory power.
\(H_A\): at least one coefficient is non-zero.
More generally, the \(F\)-test compares a restricted model, in which some coefficients are set to zero under the null, against an unrestricted model in which they are free. A large \(F\)-statistic means the unrestricted model explains substantially more variation, and as usual we read the \(p\)-value.
Where the \(F\)-statistic comes from
Definition: the \(F\) statistic
\[
F = \frac{(SSR_{r} - SSR_{u}) / q}{SSR_{u} / (n - k - 1)}
\]
\(SSR_{r}\): sum of squared residuals from the restricted model, with the null imposed.
\(SSR_{u}\): sum of squared residuals from the unrestricted model.
\(q\): the number of restrictions. Count the equals signs in \(H_0\): for \(H_0: \beta_1 = \beta_2 = 0\), \(q = 2\).
\(n\): observations; \(k\): regressors in the unrestricted model, excluding the intercept.
The numerator is how much worse the fit gets when you impose the null, per restriction. The denominator scales that by how much unexplained variation there was to begin with.
The \(F\) distribution
The distribution takes two degrees-of-freedom parameters, built from \(q\), \(n\), and \(k\). With few observations relative to parameters, the right tail is heavier, so a larger \(F\) is required to reject the null at a given significance level.
An \(F\)-test by hand
Take the fully specified model and ask whether experience and ability jointly matter, so \(H_0: \beta_{exper} = \beta_{ability} = 0\) and \(q = 2\).
We build the statistic from the formula rather than from a built-in test. It requires two regressions and their sums of squared residuals.
m_restr <-feols(wage ~ educ, data = dat) # H0 imposedm_unrestr <-feols(wage ~ educ + exper + ability, data = dat) # H0 relaxedssr_r <-sum(resid(m_restr)^2)ssr_u <-sum(resid(m_unrestr)^2)q <-2# two restrictionsk <-3# regressors in the unrestricted model, excluding the interceptf_stat <- ((ssr_r - ssr_u) / q) / (ssr_u / (nobs(m_unrestr) - k -1))c(SSR_restricted = ssr_r, SSR_unrestricted = ssr_u, F = f_stat)
SSR_restricted SSR_unrestricted F
24484.2994 17882.5404 368.4351
Code
# fixest will of course do all of this for us:wald(m_unrestr, "exper|ability")
Wald test, H0: joint nullity of exper and ability
stat = 368.4, p-value < 2.2e-16, on 2 and 1,996 DoF, VCOV: IID.
Code
import pyfixest as pffrom scipy.stats import f as f_distm_restr = pf.feols("wage ~ educ", data=r.dat) # H0 imposedm_unrestr = pf.feols("wage ~ educ + exper + ability", data=r.dat) # H0 relaxedssr_r = (m_restr.resid() **2).sum()ssr_u = (m_unrestr.resid() **2).sum()q, k, N =2, 3, r.dat.shape[0]f_stat = ((ssr_r - ssr_u) / q) / (ssr_u / (N - k -1))print(f"SSR_restricted {ssr_r:.1f}, SSR_unrestricted {ssr_u:.1f}, F {f_stat:.1f}")
SSR_restricted 24484.3, SSR_unrestricted 17882.5, F 368.4
Code
print(f"p-value {f_dist.sf(f_stat, q, N - k -1):.3g}")
p-value 6.52e-137
The SSR rises by 37% when we drop the two variables, giving \(F \approx 368\). We reject, and wald() confirms the hand computation exactly.
The \(F\)-test with a single restriction
With one restriction the \(F\)-test reproduces the \(t\)-test. If \(q = 1\) then \(F = t^2\), and the two report the same \(p\)-value. The code checks this on the estimated model:
The printed \(F\)-statistic assumes homoskedasticity
The \(F\)-statistic reported in standard regression output is built from sums of squared residuals, in the same way as the default standard errors. It therefore assumes a constant error variance.
Once the classical assumptions fail, which is the subject of the second half of this lecture, that statistic is no longer valid and a robust version is required. wald() reports it once a robust vcov has been set.
Comparing the two tests
A group of variables can be jointly significant even when no single one of them is individually significant, so both tests are needed.
\(t\)-test
\(F\)-test
Scope
One coefficient
Two or more coefficients
Typical use
Is this variable significant?
Is this group jointly significant, or the model as a whole?
Null hypothesis
\(H_0: \beta_j = 0\)
\(H_0: \beta_1 = \beta_2 = \dots = 0\)
Statistic
\(t = \hat{\beta}_j / se(\hat{\beta}_j)\)
Compares restricted and unrestricted sums of squares
Example question
Does education affect wages, holding other factors constant?
Do exper and exper^2 jointly affect wages?
The standard case for the \(F\)-test is two strongly collinear regressors. Neither is individually distinguishable from zero, while the two are jointly significant.
Dummies and interactions
A dummy shifts the intercept
Lecture 2 introduced a dummy variable as the only regressor. We now add it to a model that also contains a continuous regressor:
Two lines with a different intercept and the same slope. So \(\beta_0\) is the intercept for the reference group, \(\beta_1\) the wage difference between women and men with the same education, and \(\beta_2\) the return to education, which this model forces to be equal for both.
Parallel slopes
The two lines are separated vertically by exactly \(\beta_1\), at every level of education.
Interactions let the slope differ too
The parallel-slopes model assumes the return to education is identical for men and women. What if it is not?
To let the slope differ, add the product of the dummy and the continuous variable:
So women have intercept \((\beta_0 + \beta_1)\) and slope \((\beta_2 + \beta_3)\).
Reading an interaction
The presence of an interaction changes what the “main effects” mean, which is where most mistakes happen.
\(\beta_2\) is the effect of education for the reference group only, not for everyone.
\(\beta_1\) is the wage difference between women and men when \(Educ = 0\). If zero education never occurs in the data, this number is not interpretable on its own.
\(\beta_3\), the interaction coefficient, is the difference in slopes. It says by how much the return to education differs for women, and it is usually the coefficient of interest.
So the effect of education is \(\beta_2\) for men and \(\beta_2 + \beta_3\) for women.
Different slopes, visually
The lines now converge, so the size of the gap depends on where you stand on the horizontal axis.
The effect depends on another variable
With an interaction, the effect of a variable is a function of another variable rather than a single coefficient.
Two derived quantities
The marginal effect of education for women: \[
\frac{\partial E[wage \mid F = 1, Educ]}{\partial Educ} = \beta_2 + \beta_3.
\]
The wage gap between women and men: \[
\begin{aligned}
E[wage \mid F=1] - E[wage \mid F=0]
&= \big[(\beta_0 + \beta_1) + (\beta_2 + \beta_3)Educ\big] - \big[\beta_0 + \beta_2 Educ\big] \\
&= \beta_1 + \beta_3 Educ.
\end{aligned}
\]
The gap is not constant. It depends on the level of education.
This determines which test to use. To test whether the return to education differs, test \(H_0: \beta_3 = 0\) with a \(t\)-test. To test whether gender has any effect at all, the two lines must coincide, so test \(H_0: \beta_1 = \beta_3 = 0\) jointly with an \(F\)-test.
More than two categories
A categorical variable with \(k\) categories, such as region, cannot simply be coded 1, 2, 3, 4. That would impose an ordering and equal spacing that do not exist.
Instead, create \(k - 1\) binary dummies. Each observation takes the value 1 on its own category’s dummy and 0 on all the others.
The dummy variable trap
One category must always be omitted, and it becomes the baseline or reference.
Including a dummy for every category produces perfect multicollinearity: knowing that an observation is not in the North, South, or East implies that it is in the West. The model cannot separate perfectly redundant regressors, so the software drops one or returns an error.
With West as the baseline: \[
y_i = \beta_0 + \beta_1 \text{North}_i + \beta_2 \text{South}_i + \beta_3 \text{East}_i + \dots + u_i.
\]
Interpreting a set of dummies
The intercept \(\beta_0\) is now the average outcome for the baseline category, with all other regressors at zero.
Each dummy coefficient is the average difference between that category and the baseline, holding other variables constant. It is not that category’s average.
Statistical significance follows the same logic: a significant \(p\)-value on the North dummy says North differs from West, not that it differs from South or East.
Example: predicting wages
With “High School” as the baseline, suppose \[
\widehat{\text{Wage}}_i = 35{,}000 + 15{,}000 \cdot \text{Bachelors}_i + 30{,}000 \cdot \text{Masters}_i.
\]
High School earns €35,000; Bachelors earn €15,000 more, so €50,000; Masters earn €30,000 more, so €65,000.
Changing the baseline
The choice of baseline changes the individual coefficients and their interpretation. It does not change the explanatory power of the categorical variable as a whole.
For a set of dummies representing one categorical predictor, the \(F\)-test evaluates the joint null that all of them are zero: does education have any explanatory power?
Whichever category is dropped, the model uses the same information and explains the same variation, so the \(R^2\) and the \(F\)-statistic are unchanged. The conclusion about whether education matters does not depend on the reference level.
The individual \(t\)-tests do depend on it. Each tests one comparison against the baseline, so changing the baseline changes which comparisons are being tested.
Three categories, visually
Each dummy coefficient is a vertical distance from the dashed baseline, and nothing else.
Reading regression output
The full regression output
We can now read all parts of the standard output. The model is the fully specified wage regression:
Coefficients.Estimate gives \(\hat{\beta}_0\) through \(\hat{\beta}_3\). Std. Error is \(se(\hat{\beta}_j)\), measuring sampling uncertainty. t value is their ratio, and Pr(>|t|) is the \(p\)-value of the two-sided test against zero.
Goodness of fit. The residual standard error is the SER from Lecture 2, here 2.99, close to the \(\sigma = 3\) we built into the data. The \(R^2\) is 0.64, so the model accounts for about 64% of the variation in wages, against 51% for the simple regression on education alone.
The rise in \(R^2\) reports improved fit. The coefficient on education also changed by about a third, and that change, rather than the fit, is what bears on the credibility of the estimate.
Interpretation with controls
Every coefficient has an economic interpretation, whatever its statistical significance. The standard error, \(t\)-statistic, confidence interval, and \(p\)-value quantify sampling uncertainty; they do not decide whether a number means anything.
Education:\(\hat{\beta}_1 \approx\) 1.21. Holding experience and ability constant, one more year of education is associated with €1.21 per hour more.
Experience:\(\hat{\beta}_2 \approx\) 0.19. Holding education and ability constant, one more year of experience is associated with €0.19 per hour more.
Ability:\(\hat{\beta}_3 \approx\) 1.98. A one standard deviation higher score is associated with €1.98 per hour more.
Comparison with the true values
Because the data were simulated, the estimates can be compared with the parameters that generated them: 1.2, 0.20, and 2.0. All three estimates lie within about one standard error of the truth.
That comparison is unavailable in real data, so the assumptions behind those standard errors carry the weight of the inference. They are the subject of the second half of the lecture.
Heteroskedasticity-robust standard errors
Where the variance comes from
The first half of the lecture took the standard errors in the output as given. We now examine the assumptions they rest on.
The variance of the OLS estimator, conditional on the regressors, is
Evaluating the second factor requires assumptions about the errors.
The two default assumptions
The default output in R and Python imposes two restrictions on the errors:
Errors are uncorrelated across observations: \(\operatorname{Cov}(u_i, u_j) = 0\) for \(i \neq j\).
Homoskedasticity: the variance is constant, \(\operatorname{Var}(u_i) = \sigma^2\).
This section relaxes assumption 2, and the next section relaxes assumption 1.
Notation in the second half
The remainder of the lecture uses the single-regressor formulas and writes \(\operatorname{Var}(\hat{\beta}_1)\) for the conditional variance \(\operatorname{Var}(\hat{\beta}_1 \mid X)\). The results carry over to the multivariate case.
The variance under homoskedasticity
Under homoskedasticity the variance reduces to the formula from Lecture 2:
The assumption is convenient, and in economic data it usually does not hold. Heteroskedasticity means that the error variance is not constant across observations. Ignoring it invalidates the inference while leaving the estimates unbiased.
Illustration: constant and non-constant spread
The two panels share the same true line. Only the spread of the observations around it differs.
OLS fits the same line in both cases and is unbiased in both. What differs is the precision of that line, and therefore the standard error that should be reported.
The variance under heteroskedasticity
Formally, heteroskedasticity means \(\operatorname{Var}(u_i) = \sigma_i^2\): the variance carries an \(i\) subscript, so it can differ for every observation.
The food expenditure example illustrates this. Low-income households have limited budgets, so their spending varies little. High-income households have discretion over how much they spend on food, so their spending varies more. The error variance rises with income.
With non-constant variance, the true variance of the OLS estimator is
This expression is exact whenever the \(\sigma_i^2\) differ, and the classical formula is the special case obtained by setting \(\sigma_i^2 = \sigma^2\) for every observation.
Estimating the variance from the residuals
The individual error variances \(\sigma_i^2\) are unknown.
White (1980) proposed replacing each \(\sigma_i^2\) by the squared OLS residual \(\hat{u}_i^2\) of that observation. A single squared residual estimates a single variance imprecisely, but the sum over observations is a consistent estimate of the sum that appears in the variance formula.
Its square root is the heteroskedasticity-robust standard error, also called the White or HC standard error. Because it relies on that summation argument, it is valid in large samples.
Robust standard errors in software
The robust variance estimator must be requested explicitly:
Heteroskedasticity is common in cross-sectional data, and the current convention is to report HC-robust standard errors by default.
The cost of using them when the errors are in fact homoskedastic is small in large samples. The cost of omitting them when they are needed is invalid inference.
What robust standard errors do not fix
HC-robust standard errors correct for heteroskedasticity and for nothing else. They do not correct for serial correlation, and they do not correct for omitted variable bias.
For serial correlation in time series, HAC or Newey-West standard errors are the usual choice. When the dependence follows a group structure, clustered standard errors are appropriate, which is the subject of the next section.
Three common misconceptions
“Robust standard errors fix my estimates.” They do not. OLS remains unbiased and consistent under heteroskedasticity; what it loses is efficiency. Robust standard errors correct the inference and leave the estimator unchanged. Recovering efficiency requires a different method, such as weighted least squares.
“Robust standard errors are always larger.” They are larger when observations with unusual \(x\) values have large residuals, which is common, but they can also come out smaller. Robust means valid under weaker assumptions, and does not imply conservative.
“Changing vcov changed my results.” The variance estimator cannot change \(\hat{\beta}\). The point estimates are identical for every choice of variance estimator, and only the standard errors change. A coefficient that moves indicates that the specification changed as well.
Clustered standard errors
Correlated errors within groups
We now turn to the other default assumption, \(\operatorname{Cov}(u_i, u_j) = 0\). In many economic datasets it does not hold: students within schools, individuals within states, firms observed over time.
Unobserved factors at the group level induce correlation among the errors within that group. Data with this structure are called clustered data.
The next slides measure, by simulation, how far the default standard error is from the true sampling standard deviation when that correlation is ignored.
The simulation design
We build 50 schools of 30 students each, so \(N = 1500\). Each student’s error has two components, \(u_{gi} = a_g + e_{gi}\): a shock \(a_g\)shared by everyone in school \(g\), plus individual noise \(e_{gi}\).
The share of error variance carried by the shared piece is
\[
\rho_u = \frac{\operatorname{Var}(a_g)}{\operatorname{Var}(a_g) + \operatorname{Var}(e_{gi})} = \operatorname{Corr}(u_{gi}, u_{gj}), \quad i \neq j \text{ in the same school},
\]
so \(\rho_u\) is exactly the intra-cluster correlation of the errors. We set \(\rho_u = 0.5\).
A policy is assigned by coin flip to whole schools, and its true effect is exactly zero. A 5% test applied to a true null should reject about 5% of the time. We draw 1000 datasets and record the rejection rate.
What the decomposition means
Each dot is one student’s error, grouped by school. The red bar marks the school’s mean error, which estimates the shared shock \(a_g\).
The spread of the red bars around zero reflects \(\operatorname{Var}(a_g)\); the spread of the dots around their own bar reflects \(\operatorname{Var}(e_{gi})\). With \(\rho_u = 0.5\) the two are equal. One student’s error therefore predicts part of a schoolmate’s error, because both contain the same \(a_g\).
Estimating \(\rho_u\) in data
The errors are unobserved, so we use the OLS residuals. Regressing the residuals on school dummies fits each school’s mean residual. The \(R^2\) of that regression is the share of residual variance that lies between schools, which estimates \(\rho_u\). Regressing \(x\) on the school dummies in the same way estimates \(\rho_x\).
In this dataset the estimate is \(\hat{\rho}_u =\) 0.505, close to the true 0.5. The estimate \(\hat{\rho}_x\) equals 1, because the policy is constant within each school.
The simulation in code
One replication draws a new dataset, runs feols once, and returns the estimate together with all three standard errors. run_sim repeats that 1000 times and records how often each standard error rejects.
Code
set.seed(2026)# One replication of a world in which the null hypothesis is TRUE.one_replication <-function(cluster_level_x) { G <-50; Ng <-30; rho_u <-0.5 g <-rep(1:G, each = Ng)# Each school gets a shared shock a_g on top of individual noise e_gi. u <-rep(rnorm(G, 0, sqrt(rho_u)), each = Ng) +rnorm(G * Ng, 0, sqrt(1- rho_u))# The policy either varies only across schools, or freely within them. x <-if (cluster_level_x) rep(rbinom(G, 1, 0.5), each = Ng) elserbinom(G * Ng, 1, 0.5) y <-0* x + u # the TRUE effect of x is exactly zero d <-data.frame(y = y, x = x, g = g) m <-feols(y ~ x, data = d, vcov ="iid")c(b =coef(m)["x"],classical =se(m)["x"],hc1 =se(summary(m, vcov ="hc1"))["x"],cr1 =se(summary(m, vcov =~g))["x"])}run_sim <-function(cluster_level_x, reps =1000) { out <-t(replicate(reps, one_replication(cluster_level_x)))colnames(out) <-c("b", "classical", "hc1", "cr1")# Software uses t with G-1 degrees of freedom for clustered tests, N-2 otherwise. cv_n <-qt(0.975, 1500-2); cv_g <-qt(0.975, 50-1)list(reject =c(classical =mean(abs(out[, "b"] / out[, "classical"]) > cv_n),hc1 =mean(abs(out[, "b"] / out[, "hc1"]) > cv_n),cr1 =mean(abs(out[, "b"] / out[, "cr1"]) > cv_g)),mean_se =colMeans(out[, c("classical", "hc1", "cr1")]),true_sd =sd(out[, "b"]) )}sim_cluster <-run_sim(TRUE)sim_individual <-run_sim(FALSE)
Results: policy assigned by school
Policy assigned at the school level, true effect zero, nominal significance level 5%:
Standard error used
Rejects a true null
Average SE reported
True sampling sd of \(\hat{\beta}_1\)
Classical (the default)
64.4%
0.052
0.213
HC1 robust
64.6%
0.052
0.213
CR1 clustered
7.5%
0.203
0.213
The default standard error is about 4.1 times too small, and a test that should reject 5% of the time rejects 64% of the time. Compare that factor with the one you predicted at the start.
Heteroskedasticity-robust standard errors do not improve on this, 64.6% against 64.4%. HC-robust errors correct the diagonal of the variance matrix, while the correlation here sits off the diagonal.
The clustered standard error comes close to 5% without reaching it. The reason is discussed at the end of the section.
Results: policy assigned within schools
Now assign the policy within schools, by flipping a coin for each student separately. Everything else, including \(\rho_u = 0.5\), is unchanged.
Standard error used
Rejects a true null
Average SE reported
True sampling sd of \(\hat{\beta}_1\)
Classical (the default)
5.4%
0.051
0.052
HC1 robust
5.4%
0.051
0.052
CR1 clustered
5.9%
0.051
0.052
Every standard error is now correct, and clustering changes nothing.
When clustering matters
Correlated errors alone leave the standard error unaffected. The default standard error breaks down only when correlated errors are combined with a regressor that is also correlated within clusters.
The question to ask is at what level the treatment varies, rather than whether the errors are correlated.
Notation for clustered data
Let \(g\) index the cluster, such as a school, and \(i\) the individual within it, so an observation is \(gi\):
\[
y_{gi} = \beta_0 + \beta_1 x_{gi} + u_{gi}.
\]
The defining feature of clustered data is the error structure:
\(\operatorname{Cov}(u_{gi}, u_{gj}) \neq 0\) for \(i \neq j\): correlation within a cluster.
\(\operatorname{Cov}(u_{gi}, u_{g'j}) = 0\) for \(g \neq g'\): no correlation between clusters.
The standard OLS formula sets all of those within-cluster covariance terms to zero. When they are on average positive, which is the common case, the variance of \(\hat{\beta}_1\) is systematically underestimated and the standard errors are too small.
What each estimator assumes
Each of the three estimators is a statement about which entries of the error covariance matrix \(E[uu']\) may be non-zero. Shaded cells are free; white cells are assumed to be zero.
Reading left to right, each estimator relaxes an assumption rather than applying a correction. Clustering allows the off-diagonal blocks to be non-zero.
How the clustered estimator is computed
The clustered estimator is computed in four steps.1 Run OLS to obtain the residuals; group them by cluster; for each cluster compute a term capturing the variances and covariances within it; sum those terms across clusters.
The difference from the HC estimator is the order of summing and squaring.
HC squares each observation’s contribution, then sums \(N\) of them.
The clustered estimator sums the contributions within a cluster first, squares that cluster total, then sums \(G\) of them.
Squaring after summing is what preserves the cross-products \(\hat{u}_{gi}\hat{u}_{gj}\) that the HC estimator sets to zero.
This treats each cluster as a single larger observation. When observations within a cluster carry similar information, the cluster contributes less independent information than the same number of independent observations would, so the resulting standard errors are typically larger.
where \(\rho_u\) is the intra-cluster correlation of the errors, \(\rho_x\) that of the regressor, and \(\bar{N}_g\) the average cluster size.
The factor equals 1, so clustering changes nothing, whenever any one of the three terms is absent: no shared shocks (\(\rho_u \approx 0\)), a regressor that varies randomly within clusters (\(\rho_x \approx 0\)), or small clusters.
The factor is largest when all three terms are large, which describes many policy evaluations: a state-level law observed on many individuals per state.
In our simulation \(\rho_u = 0.5\), \(\rho_x = 1\), and \(N_g = 30\), giving \(\sqrt{1 + 0.5 \times 29} = 3.94\). The simulation delivered 4.13.
Why the ratio grows with cluster size
Ask what happens to the average error in a school as we enrol more students:
Only the individual noise averages away. The school-wide shock \(a_g\) never does, because every student carries the same\(a_g\). A school of 1000 students is barely more informative about that school than a school of 30.
The classical formula does not account for this. It treats each additional student as an independent draw and shrinks the standard error at rate \(1/\sqrt{N}\), while precision is governed by the number of schools.
The gap therefore widens as \(N_g\) grows, and \(\sqrt{1 + \rho_u \rho_x (N_g - 1)}\) measures it: each additional student adds \(N_g - 1\) positive covariance terms that the classical formula sets to zero.
The Moulton factor, visually
A small intra-cluster correlation is sufficient. At \(\rho_u = 0.05\), which is not visible in a scatter plot, the standard error doubles once clusters contain 60 observations.
An interactive illustration
Treatment varies
Standard error
Mean
Rejects
Classical
—
—
HC1 robust
—
—
CR1 clustered
—
—
True sd of β̂1
—
5% ideal
Moulton factor √(1 + ρuρx(Ng−1)) =
—
300 datasets in which the true effect is zero. Bars are the actual spread of
β̂1; curves are what each standard error claims that spread is.
Set ρu to zero, or let treatment vary within clusters, and all three curves
collapse onto the bars.
When to cluster
The decision follows from the research design. Cluster when the treatment is assigned at the cluster level, or when the sample was drawn in clusters, as in multi-stage surveys.
Panel data satisfy both conditions, since each individual is a cluster of their own repeated observations. That is why this material directly precedes the panel lectures.
The presence of a group variable in the data is not itself a reason to cluster. As the second simulation showed, a regressor that varies randomly within groups leaves the standard errors unaffected.
Failing to cluster when the design requires it biases the standard errors downwards, sometimes by a factor of four or more, which produces overstated precision and spurious findings.
At what level should you cluster?
Cluster at the level at which the unobserved components are shared, and at which treatment was assigned.
With students in classrooms within schools, and a policy rolled out by school, cluster on school. This allows for correlation within classrooms and between classrooms in the same school.
Clustering at too fine a level, here the classroom, leaves the between-classroom correlation unaccounted for, so the standard errors remain too small.
The constraint in the other direction is that the chosen level must contain enough clusters, which is the subject of the next slide.
Clustering in practice
Clustering does not change \(\hat{\beta}\). Only the standard errors, \(t\)-statistics, and \(p\)-values change. Coefficients that move when the vcov changes indicate that something else in the specification changed.
Cluster-robust standard errors require many clusters rather than many observations. They are consistent as \(G \to \infty\), not as \(N \to \infty\). A dataset with 500,000 observations in 12 states has 12 clusters, and its clustered standard errors remain biased downwards. This is why the simulation rejected 7.5% of the time rather than exactly 5%, even with \(G = 50\).
With fewer than roughly 40 clusters, use the wild cluster bootstrap. In R, fwildclusterboot::boottest(); in Python, pf.feols(...).wildboottest(reps = 9999, param = "x", cluster = "cluster_variable").
Clustering in software
Like HC-robust standard errors, clustered standard errors must be requested explicitly:
feols changes its default once fixed effects are added. Writing feols(y ~ x | firm, data = df) clusters on firm automatically, without that being requested.
Read the vcov line that fixest prints under the coefficient table rather than assuming the default. Setting vcov explicitly makes the choice visible.
Summary
What did we learn?
Multiple regression holds the other measured factors constant.
A valid control is a confounder; a mediator or a collider adds bias. Selecting a sample on a variable is equivalent to controlling for it.
Bias is removed by the right control rather than by more controls.
Controls cost precision through \(R_j^2\), and the most useful control is often the most expensive.
The \(t\)-test concerns one coefficient and the \(F\)-test several. With one restriction, \(F = t^2\).
An interaction makes the effect a function of another variable.
Default standard errors assume constant variance and uncorrelated errors.
HC-robust errors relax the first assumption and cluster-robust errors the second. Neither changes \(\hat{\beta}\) or removes bias.
Exit questions
Answer these questions without referring to the preceding slides:
Adding experience changed the education coefficient by almost nothing. Why?
Why is controlling for a collider harmful when the collider is correlated with both variables?
When are a \(t\)-test and an \(F\)-test guaranteed to give the same \(p\)-value?
In an interacted model, what does the coefficient on the dummy alone mean?
Our simulation had \(\rho_u = 0.5\) throughout, yet clustering mattered in one design and not the other. What differed?
You have 500,000 observations in 12 states, and the policy varies by state. What is your effective sample size for inference?
Answers to the exit questions
Why experience did nothing. Bias is \(\beta_2 \delta_1\), and experience is uncorrelated with education, so \(\delta_1 \approx 0\). A valid control that is orthogonal to the treatment cannot have caused bias.
Why a collider harms. It is a consequence of both variables rather than a cause of either. Conditioning on a common effect creates an association between its causes where none existed.
When \(t\) and \(F\) agree. When there is exactly one restriction, \(q = 1\), in which case \(F = t^2\) identically.
The dummy in an interacted model. The group difference evaluated at zero on the interacted variable. If zero is outside the data, it is not interpretable on its own.
What differed between the designs.\(\rho_x\). With treatment assigned to whole schools \(\rho_x = 1\); assigned within schools, \(\rho_x \approx 0\), and the Moulton factor collapses to 1.
Effective sample size. 12. Cluster-robust inference is consistent in the number of clusters, so those clustered standard errors are still biased downwards, and this is a case for the wild cluster bootstrap.
Next lecture
Controls work only for confounders that can be measured. Ability was observable here only because the data were simulated.
Lecture 4 introduces fixed effects, which control for unmeasured confounders provided that they do not change over time.
Appendix: cluster-robust SEs
The setting and the assumption
This appendix derives the cluster-robust formula used in the lecture. It adds no new economics, and shows where the formula comes from and why the summation sits inside the square.
The setting. There are \(G\) clusters, indexed by \(g\), and cluster \(g\) holds \(N_g\) observations indexed by \(i\). As in the main text we keep a single regressor:
The assumption. Errors may be correlated in any way within a cluster, but not at all across clusters:
\[
\operatorname{Cov}(u_{gi}, u_{gj}) \ \text{ unrestricted}, \qquad \operatorname{Cov}(u_{gi}, u_{g'j}) = 0 \ \text{ for } g \neq g'.
\]
The variance in bread-and-meat form
Every variance formula in this section has the same bread and meat structure: a factor built only from the \(x\)’s, multiplying a term built from the errors.
The derivation has three steps:
Write \(\operatorname{Var}(\hat{\beta}_1)\) as bread times meat.
Split the meat cluster by cluster, which is where the covariance terms remain.
Replace the unknown errors \(u_{gi}\) by the observed residuals \(\hat{u}_{gi}\).
Colours used on the slides that follow
\(\color{#0066a1}{\textbf{Blue}}\): the bread, \(\left(\sum (x_{gi} - \bar{x})^2\right)^{-2}\). Identical in the classical, HC and clustered estimators.
\(\color{#2e7d32}{\textbf{Green}}\): the within-cluster sum, the “sum inside” step that creates clustering.
\(\color{#c43c39}{\textbf{Red}}\): the within-cluster covariance terms, precisely what the classical and HC estimators set to zero.
Expressing the variance
We express \(\operatorname{Var}(\hat{\beta}_1)\) in a form that recognises the cluster structure:
The blue factor is the bread, and it is now fixed for the rest of the derivation. In the final line we only re-indexed: summing over all \(N\) observations is the same as summing within each cluster and then over clusters.
The meat, cluster by cluster
Restating the result, with the meat still to be evaluated:
The black line is the only part the HC estimator retains. The red line is the entire difference, and it is what a classical or HC standard error sets to zero.
From errors to residuals
We cannot compute the true variance, because the errors \(u_{gi}\) are unknown. We estimate it from the OLS residuals \(\hat{u}_{gi} = y_{gi} - \hat{y}_{gi}\).
Define a score for each observation, \(s_{gi} = (x_{gi} - \bar{x})\hat{u}_{gi}\). It is built from the observed residual rather than from the unobserved error. Let the cluster total be
This green object is the sample counterpart of the green term on the previous slides.
Squaring after summing
The single-cluster variance can be written as an expected squared sum, \(E\left[ \left( \sum_{i} (x_{gi} - \bar{x})u_{gi} \right)^2 \right]\).1
A natural estimator is the squared cluster total, \(\color{#2e7d32}{(s_g)}^2\). Squaring after summing is the step that keeps the red covariance terms: expanding \((s_{g1} + \dots + s_{gN_g})^2\) produces every cross-product \(s_{gi}s_{gj}\).
Summing over clusters estimates the meat in Equation 1:
The blue bread is the same in both. The only difference is where the summation sits relative to the square.
HC squares each observation separately and sums \(N\) of them, which assumes independence. The clustered estimator sums the green term inside first, squares that cluster total, and sums \(G\) of them. Expanding the square recovers every \(\color{#c43c39}{\text{red covariance term}}\) without estimating any of them individually.
