Generalized Least Squares and Clustered Dependence

Lecture 7

Author
Affiliation

Henrique Veras

PIMES/UFPE

0.1 Where we are

Where this lecture sits

DGPIdentificationEstimationAsymptoticsInference
Computation

Lectures 4–5 already made heteroskedasticity the default: \(\E[\beps\beps'\mid\bX]=\bD\) is a general positive-definite matrix, OLS is consistent, and the robust (sandwich) standard errors deliver valid inference with no homoskedasticity assumption. So the classical “violation of A4” story is behind us. Two questions remain — and they are what this lecture is about.

The reordering of the course did the heavy lifting here. What used to be “the heteroskedasticity lecture” — deriving the sandwich variance, White’s estimator, showing that \(s^2(\bX'\bX)^{-1}\) is wrong — now lives in Lectures 4 and 5, as the general case. Building on that, this lecture asks the two things robust OLS does not settle: whether we can estimate \(\bbeta\) more efficiently, and what to do when errors are correlated across observations.

Two questions robust OLS does not answer

The two questions have very different answers. Efficiency leads to generalized least squares — a better estimator if we know (or can model) \(\bD\), but one that trades robustness for that efficiency. Dependence leads to cluster-robust inference — not a better estimator at all, but a different variance, because the design of the sample changed.

1 Part I · The generalized model and the cost to OLS

1.1 The generalized regression model

One assumption, dropped long ago

The model is the familiar linear one, with the error covariance left general: \[\bY=\bX\bbeta+\beps,\qquad \E[\beps\mid\bX]=\bzero,\qquad \E[\beps\beps'\mid\bX]=\bD,\] with \(\bD\) symmetric positive-definite (an \(n\times n\) matrix). Homoskedasticity, \(\bD=\sigma^2\bI\), is the special case — not the starting point.

Definition 1 (The two leading structures) Heteroskedasticity (this lecture’s focus): errors uncorrelated across observations but with different variances — \(\bD=\operatorname{diag}(\sigma_1^2,\dots,\sigma_n^2)\), i.e. \(\E[\varepsilon_i^2\mid\mathbf{x}_i]=\sigma_i^2(\mathbf{x}_i)\). Dependence: \(\bD\) has non-zero off-diagonal entries — serial correlation (time series) or within-cluster correlation (Part V).

OLS still works — it just is not the best

Nothing about consistency or valid inference breaks under a general \(\bD\):

This is the exact sense in which the lecture is about efficiency, not validity. OLS gives you the right answer with honest uncertainty. The open question is whether another linear unbiased estimator has a smaller variance — and Lecture 5’s “asymptotic efficiency” slide already told us where to look: the estimator that attains the efficiency bound is generalized least squares.

2 Part II · Generalized Least Squares

2.1 The cost of ignoring D

Every observation weighted equally

A general \(\bD\) can break A4 in two ways — heteroskedasticity (the diagonal of \(\bD\) varies across observations) or autocorrelation (non-zero off-diagonals), or both. OLS stays unbiased and consistent either way, but it pays for ignoring \(\bD\) in precision. Compare the variance OLS assumes with the true one: \[\underbrace{\sigma^2(\bX'\bX)^{-1}}_{\text{under A4}} \qquad\text{vs.}\qquad \underbrace{(\bX'\bX)^{-1}\bX'\bD\bX(\bX'\bX)^{-1}}_{\text{true, general }\bD}.\]

Following Greene: the conventional variance is built from \(\sum_i\mathbf{x}_i\mathbf{x}_i'\) — weight \(1\) on every \(\mathbf{x}_i\mathbf{x}_i'\) — whereas the true variance weights the \(i\)-th term by \(\sigma_i^2\) (through \(\bX'\bD\bX\)). The two agree only when all \(\sigma_i^2\) are equal, i.e. A4. Otherwise OLS spreads its attention equally over observations that deserve unequal trust, and a re-weighted estimator must do better. That estimator is GLS; the weight that “should” appear, \(1/\sigma_i^2\), emerges on its own.

2.2 GLS with known D

Transforming to constant variance

If the trouble is unequal (and correlated) variances, the cure is to rescale the observations so the variances become equal and the correlations vanish — then the classical model, and OLS’s optimality, are back. Since \(\bD\) is symmetric positive-definite, it has a spectral decomposition \(\bD=\mathbf{C}\boldsymbol{\Lambda}\mathbf{C}'\), with \(\mathbf{C}\) orthogonal (\(\mathbf{C}'\mathbf{C}=\bI\), columns the eigenvectors) and \(\boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1,\dots,\lambda_n)\), \(\lambda_i>0\). Following Greene, build \(\mathbf{P}=\boldsymbol{\Lambda}^{-1/2}\mathbf{C}'\), which satisfies \[\mathbf{P}\,\bD\,\mathbf{P}'=\bI,\qquad\text{equivalently}\qquad \bD^{-1}=\mathbf{P}'\mathbf{P}\] — the matrix analogue of “dividing by the standard deviation”, now along the principal axes of \(\bD\).

On the Board

Pre-multiply the model by \(\mathbf{P}\): with \(\bY_*=\mathbf{P}\bY\), \(\bX_*=\mathbf{P}\bX\), \(\beps_*=\mathbf{P}\beps\), \[\bY_*=\bX_*\bbeta+\beps_*,\qquad \E[\beps_*\beps_*'\mid\bX]=\mathbf{P}\,\bD\,\mathbf{P}'=\bI.\] Two moves: \(\mathbf{C}'\) rotates the errors into the eigenbasis (decorrelating them — \(\E[\mathbf{C}'\beps\,\beps'\mathbf{C}\mid\bX]=\boldsymbol{\Lambda}\)), then \(\boldsymbol{\Lambda}^{-1/2}\) rescales each axis by \(1/\sqrt{\lambda_i}\) (equalizing the variances to \(1\)). The transformed errors are homoskedastic and uncorrelated — the classical model holds for \((\bY_*,\bX_*)\), so OLS on the transformed data is efficient. That is GLS.

The GLS estimator

Theorem 1 (Generalized least squares) OLS on the transformed model is \[\hat{\bbeta}_{\text{GLS}} =(\bX_*'\bX_*)^{-1}\bX_*'\bY_* =(\bX'\bD^{-1}\bX)^{-1}\bX'\bD^{-1}\bY.\] Among linear unbiased estimators it has the smallest variance — the generalized Gauss–Markov (Aitken) theorem — with \[\Var[\hat{\bbeta}_{\text{GLS}}\mid\bX]=(\bX'\bD^{-1}\bX)^{-1}.\] OLS is the special case \(\bD=\sigma^2\bI\) (weight matrix \(\bI\) instead of \(\bD^{-1}\)).

GLS re-weights: it downweights high-variance observations and upweights precise ones, which is exactly the information OLS ignores by treating every observation alike. The efficiency gain is the payoff Lecture 5 promised — \(\Var[\hat{\bbeta}_{\text{GLS}}]\preceq\Var[\bb]\) in the psd order, whenever \(\bD\neq\sigma^2\bI\). Under i.i.d. sampling the result is stronger still (Hansen’s modern Gauss–Markov): \((\bX'\bD^{-1}\bX)^{-1}\) is the smallest variance among all unbiased estimators — not only linear ones — a semiparametric efficiency bound.

Weighted least squares: the heteroskedastic case

Corollary 1 (Weighted least squares) When the violation is pure heteroskedasticity, \(\bD=\operatorname{diag}(\sigma_1^2,\dots,\sigma_n^2)\), the transformation is just dividing each observation by its own \(\sigma_i\) (\(\mathbf{P}=\operatorname{diag}(1/\sigma_1,\dots,1/\sigma_n)\)). GLS then takes the familiar weighted form \[\hat{\bbeta}_{\text{WLS}} =\Big(\textstyle\sum_i \tfrac{1}{\sigma_i^2}\,\mathbf{x}_i\mathbf{x}_i'\Big)^{-1}\Big(\textstyle\sum_i \tfrac{1}{\sigma_i^2}\,\mathbf{x}_i y_i\Big) \ -\ \text{exactly the weighting OLS was missing.}\]

3 Part III · Feasible GLS

3.1 When D is unknown

The skedastic regression

GLS needs \(\bD\), which we never know. But an unrestricted \(\bD\) has \(n\) variances — too many to estimate from \(n\) observations. Model the variance with a few parameters:

Definition 2 (The skedastic regression) Specify \(\sigma_i^2=\sigma^2(\mathbf{z}_i,\boldsymbol{\alpha})\) for a small parameter vector \(\boldsymbol{\alpha}\) and observed \(\mathbf{z}_i\) (often a subset/transform of \(\mathbf{x}_i\)) — for example the linear \(\sigma_i^2=\alpha_0+\mathbf{z}_i'\boldsymbol{\alpha}_1\) or the multiplicative \(\sigma_i^2=\exp(\alpha_0+\mathbf{z}_i'\boldsymbol{\alpha}_1)\). Estimate \(\boldsymbol{\alpha}\) by regressing the squared OLS residuals \(\hat e_i^2\) on \(\mathbf{z}_i\) — the skedastic regression.

The squared residual \(\hat e_i^2\) is a noisy but consistent proxy for \(\sigma_i^2\) (OLS is consistent, so \(\hat e_i\) tracks \(\varepsilon_i\)). Regressing \(\hat e_i^2\) on \(\mathbf{z}_i\) recovers how the variance moves with \(\mathbf{z}_i\). Replacing the true \(\sigma_i^2\) by the fitted \(\hat\sigma_i^2\) is what makes GLS feasible.

The FGLS estimator

Theorem 2 (Feasible GLS) Form \(\widehat{\bD}=\operatorname{diag}(\hat\sigma_1^2,\dots,\hat\sigma_n^2)\) from the skedastic regression and plug it in: \[\hat{\bbeta}_{\text{FGLS}}=(\bX'\widehat{\bD}^{-1}\bX)^{-1}\bX'\widehat{\bD}^{-1}\bY.\] If the skedastic model is correctly specified, estimating \(\boldsymbol{\alpha}\) is asymptotically free: \(\hat{\bbeta}_{\text{FGLS}}\) has the same asymptotic distribution as infeasible GLS.

On the Board

Why “asymptotically free”: the sampling error in \(\hat{\boldsymbol{\alpha}}\) enters \(\hat{\bbeta}_{\text{FGLS}}\) only through \(\widehat{\bD}\pto\bD\), and a \(\sqrt{n}\)-consistent \(\hat{\boldsymbol{\alpha}}\) perturbs \(\hat{\bbeta}\) by a term that vanishes faster than \(1/\sqrt{n}\). So FGLS and GLS share the limit \(\Normal(\bzero,(\bX'\bD^{-1}\bX/n)^{-1})\).

3.2 Why not always FGLS?

The efficiency–robustness trade-off

If FGLS is efficient and asymptotically free, why is OLS still the default? Three reasons — this is the heart of the lecture.

OLS is a robust estimator. The deepest reason is the third. If the equation is a projection (a best linear approximation) rather than a true CEF, OLS still estimates that projection consistently — but weighting by \(\widehat{\bD}^{-1}\) changes the target, so FGLS estimates something else. Efficiency and robustness pull in opposite directions: FGLS buys a smaller variance by betting the conditional mean is exactly right. The standard advice is OLS with robust (or cluster-robust) standard errors — the safe default — with FGLS reserved for cases where the variance structure is genuinely known.

This resolves the guiding question. Robust inference already handles heteroskedasticity for validity; FGLS offers efficiency, but conditionally and at a real cost. So the practical hierarchy is: OLS + robust by default; WLS/FGLS when you have a credible, verified model for \(\bD\) (e.g. known grouping, survey weights, a physical variance law).

4 Part IV · Testing for heteroskedasticity

4.1 The skedastic regression, tested

White and Breusch–Pagan are one test

Theorem 3 (Testing constant variance) Homoskedasticity is \(H_0:\boldsymbol{\alpha}_1=\bzero\) in the skedastic regression \(\hat e_i^2=\alpha_0+\mathbf{z}_i'\boldsymbol{\alpha}_1+\nu_i\) — the variance does not depend on \(\mathbf{z}_i\). Test it with the joint Wald/\(F\) on \(\boldsymbol{\alpha}_1\) (equivalently \(nR^2\dto\chi^2_q\), \(q=\dim\boldsymbol{\alpha}_1\)). The two classic choices differ only in \(\mathbf{z}_i\): White uses the squares and cross-products of \(\mathbf{x}_i\); Breusch–Pagan uses a general \(\mathbf{z}_i\). Under independence they are essentially the same test.

Both read off a single auxiliary regression of squared residuals. White’s choice of \(\mathbf{z}_i\) makes the test omnibus (it also picks up misspecification); Breusch–Pagan targets a direction you specify. The statistic is the joint significance of that regression — nothing new relative to Part II of Lecture 6.

What the test is — and is not — for

Under the Hood—Do not misuse the test

A heteroskedasticity test answers one scientific question: does \(\sigma^2\) depend on the regressors? It is not a tool to choose between OLS and FGLS, nor between classical and robust standard errors — “hypothesis tests are not designed for these purposes.” You use robust errors regardless of the test. Run the test only if whether \(\sigma^2(\mathbf{x})\) varies is itself of economic interest.

5 Part V · Clustered dependence

5.1 When observations come in groups

Clusters: dependence by design

Definition 3 (Clustered sampling) The sample is \(C\) groups (clusters) — firms in industries, pupils in schools, people in villages: \[y_{i,c}=\mathbf{x}_{i,c}'\bbeta+\varepsilon_{i,c}.\] Errors are correlated within a cluster but independent across them — this is dependence, a trait of the sampling design, not heteroskedasticity (not a statement about each observation’s own variance). Stacking by cluster, \(\bD=\E[\beps\beps'\mid\bX]\) is block-diagonal: each cluster its own block \(\bD_c=\E[\beps_c\beps_c'\mid\bX_c]\), zeros between blocks. OLS stays unbiased when \(\E[\beps_c\mid\bX_c]=\bzero\) (within-cluster interactions must be in the model — a pupil’s outcome must not depend on classmates’ regressors).

Where does \(\rho\) come from? Write the error as a cluster component plus an idiosyncratic one, \[\varepsilon_{i,c}=u_c+\nu_{i,c},\qquad \Var(u_c)=\sigma_u^2,\quad \Var(\nu_{i,c})=\sigma_\nu^2,\] so \(\sigma^2=\sigma_u^2+\sigma_\nu^2\), \(\Cov(\varepsilon_{i,c},\varepsilon_{j,c})=\sigma_u^2\), and the intra-cluster (intraclass) correlation is \(\rho=\sigma_u^2/(\sigma_u^2+\sigma_\nu^2)\). The within-cluster covariance is then exchangeable, \[\bD_c=\sigma^2\big[(1-\rho)\bI_N+\rho\,\mathbf{J}_N\big],\qquad \mathbf{J}_N=\mathbf{1}_N\mathbf{1}_N',\] with one large eigenvalue \(\sigma^2[1+\rho(N-1)]\) along the common mode \(\mathbf{1}_N\) and \(N-1\) small ones \(\sigma^2(1-\rho)\) orthogonal to it. The cluster sum lives in that common mode, so \[\Var(S_c)=\mathbf{1}_N'\bD_c\mathbf{1}_N=N\sigma^2\big[1+\rho(N-1)\big]=N^2\sigma_u^2+N\sigma_\nu^2.\] That is why a group-level regressor is where the damage concentrates — it, too, lives in the common mode.

How bad can it get? For a general regressor the inflation is governed by the product of the error and regressor intra-cluster correlations — the design effect (Moulton) \[\tau=1+(N-1)\,\rho_x\,\rho_\varepsilon\qquad\Longrightarrow\qquad \text{s.e.}\times\sqrt{\tau},\] with \(\rho_x=1\) for a cluster-constant regressor (worst case) and \(\rho_x\approx 0\) for an individual one (almost no inflation). For a group-level regressor, \(\sqrt{1+(N-1)\rho_\varepsilon}\) is:

\(N\) \(\rho_\varepsilon=0.05\) \(0.10\) \(0.25\)
10 1.20 1.38 1.80
30 1.57 1.97 2.87
50 1.86 2.43 3.64
100 2.44 3.30 5.07

At \(\rho_\varepsilon=0.25,\ N\approx 48\) the correct standard error is about three to four times the conventional one — ignoring clustering does not nudge the errors, it can shrink them several-fold. Equivalently, the design effect makes the effective sample size \(n/[1+(N-1)\rho]\approx C\) (the number of clusters) for a group-level regressor, which is why the relevant asymptotics is \(C\to\infty\), not \(n\to\infty\) (next slides).

A map of \(\bD\): block-diagonal

Three clusters (here of size two) — correlated within a block (\(\blacksquare\)), independent across blocks (\(0\)): \[\bD=\left[\begin{array}{cc|cc|cc} \blacksquare&\blacksquare&0&0&0&0\\ \blacksquare&\blacksquare&0&0&0&0\\\hline 0&0&\blacksquare&\blacksquare&0&0\\ 0&0&\blacksquare&\blacksquare&0&0\\\hline 0&0&0&0&\blacksquare&\blacksquare\\ 0&0&0&0&\blacksquare&\blacksquare \end{array}\right] \qquad \bD=\begin{bmatrix}\bD_1&&\\&\bD_2&\\&&\bD_3\end{bmatrix}\]

The picture is the whole argument. Within a block the errors share the cluster’s common shocks, so the entries are non-zero; across blocks the sampling makes them independent, so those entries are exactly zero. White’s estimator forms \(\sum_i \hat e_i^2\,\mathbf{x}_i\mathbf{x}_i'\) — a sum over the diagonal alone — and so throws away precisely the \(\blacksquare\) off-diagonal terms that clustering creates. The fix (next slide) sums over blocks: \(\sum_c \bX_c'\hat{\beps}_c\hat{\beps}_c'\bX_c\).

The cluster-robust variance

Theorem 4 (Cluster-robust covariance) With \(\bX_c,\beps_c\) the regressors and errors of cluster \(c\), OLS is unchanged, and \[\Var[\bb\mid\bX]=(\bX'\bX)^{-1}\Big[\textstyle\sum_{c=1}^{C}\bX_c'\,\bD_c\,\bX_c\Big](\bX'\bX)^{-1}, \qquad \bD_c=\E[\beps_c\beps_c'\mid\bX_c].\] It is estimated by the cluster sandwich, summing the cluster score vectors \(\bX_c'\hat{\be}_c\): \[\widehat{\Asyvar}[\bb]=(\bX'\bX)^{-1}\Big[\textstyle\sum_{c=1}^{C}(\bX_c'\hat{\be}_c)(\bX_c'\hat{\be}_c)'\Big](\bX'\bX)^{-1}.\] White’s estimator is the special case of one observation per cluster (\(\bD_c\) scalar).

The construction mirrors White’s exactly, but the “meat” now sums outer products of cluster totals \(\bX_c'\hat{\be}_c\) instead of individual \(\hat e_i^2\mathbf{x}_i\mathbf{x}_i'\) — that is what captures the within-cluster covariances. Everything else — the bread \((\bX'\bX)^{-1}\), the sandwich shape — is identical.

Inside a block. With \(\beps_c=(\varepsilon_{1,c},\dots,\varepsilon_{n_c,c})'\), the block \(\bD_c=\E[\beps_c\beps_c'\mid\bX_c]\) is the \(n_c\times n_c\) matrix \[\bD_c=\begin{bmatrix}\sigma_{1,c}^2 & \sigma_{12,c} & \cdots\\[2pt] \sigma_{12,c} & \sigma_{2,c}^2 & \cdots\\[2pt] \vdots & & \ddots\end{bmatrix}, \qquad \sigma_{i,c}^2=\E[\varepsilon_{i,c}^2\mid\bX_c],\quad \sigma_{ij,c}=\E[\varepsilon_{i,c}\varepsilon_{j,c}\mid\bX_c].\] The diagonal holds the conditional variances; the off-diagonal \(\sigma_{ij,c}\) are the within-cluster covariances — zero under heteroskedasticity alone, non-zero under clustering.

Why it works without estimating \(\bD_c\). A block has \(n_c(n_c+1)/2\) distinct entries, and each covariance would be “estimated” from a single product \(\hat e_{i,c}\hat e_{j,c}\) — far too noisy. We never form a \(\widehat{\bD}_c\). Instead, with the true errors, \[\E\!\big[(\bX_c'\beps_c)(\bX_c'\beps_c)'\mid\bX_c\big]=\bX_c'\bD_c\bX_c,\] so the single outer product \(\mathbf{g}_c\mathbf{g}_c'\) (with \(\mathbf{g}_c=\bX_c'\hat{\be}_c\)) is approximately unbiased for the block \(\bX_c'\bD_c\bX_c\). It is a rank-one, hopeless estimate of \(\bD_c\) entry by entry — but that is irrelevant: what we need is the sum \(\sum_c\bX_c'\bD_c\bX_c\), and \(\sum_c\mathbf{g}_c\mathbf{g}_c'\) is consistent for it as \(C\to\infty\). Consistency comes from averaging over many independent clusters, not from pinning down any single \(\bD_c\).

Where the \(1/C\) in that statement comes from. The estimator itself is the sum \(\sum_c\mathbf{g}_c\mathbf{g}_c'\) — but a sum of \(C\) positive matrices grows with \(C\), so the law of large numbers cannot act on it directly. It acts on the average: both \(\tfrac1C\sum_c\mathbf{g}_c\mathbf{g}_c'\pto\E[\mathbf{g}_c\mathbf{g}_c']\) and \(\tfrac1C\bX'\bX\pto\E[\bX_c'\bX_c]\) converge. Since \(\bX'\bX\) also grows like \(C\), each bread carries a hidden \(\tfrac1C\): \((\bX'\bX)^{-1}=\tfrac1C\big(\tfrac1C\bX'\bX\big)^{-1}\). Factoring the sandwich, the three factors of \(C\) collapse to a net \(1/C\) — the correct order of \(\Var[\bb]\) — so the normalization was already inside the two breads; it is not added to the estimator. “Consistent as \(C\to\infty\)” is thus a statement about the averaged meat \(\tfrac1C\sum_c\mathbf{g}_c\mathbf{g}_c'\), and the object with a finite limit is \(C\,\widehat{\Asyvar}[\bb]\).

The catch: asymptotics in the number of clusters

Under the Hood—Few clusters break the sandwich

The cluster-robust variance is consistent as the number of clusters \(C\to\infty\) — not as \(n\to\infty\). Stata’s finite-sample correction multiplies it by \(a_n=\dfrac{n-1}{n-k}\cdot\dfrac{C}{C-1}\) (the \(C/(C-1)\) helps when \(C\) is small; \(C=n\) recovers the i.i.d. factor). But no correction rescues small \(C\): \(C\) is the effective sample size — \(C=50\) clusters is like heteroskedasticity-robust inference with \(n=50\) observations, and unequal cluster sizes make it worse. Treat it as a small-sample problem, not a fixed reference distribution.

The reliable finite-sample tools are resampling, in Lecture 8. The wild cluster bootstrap multiplies each cluster’s residual vector by a single \(\pm1\) (Rademacher) weight — one draw per cluster, preserving within-cluster dependence — then re-computes the statistic (the restricted version, under \(H_0\), is recommended). The delete-cluster jackknife drops one whole cluster at a time; both are formula-free alternatives to the sandwich.

5.2 At what level to cluster?

The cluster level is a choice

The cluster level is a choice, and it is a genuine trade-off — not “always cluster more”.

The trade-off, made precise. The chosen level partitions the data into blocks \(b\), and the cluster-robust meat keeps within-block covariances and zeros the rest, \(\hat{\mathbf M}=\sum_b(\bX_b'\hat{\be}_b)(\bX_b'\hat{\be}_b)'\). Write \(\mathbf M=\sum_c\bX_c'\bD_c\bX_c\) for the meat at the true dependence level.

Clustering too fine (sub-clusters \(s\) of a true cluster) keeps only the within-sub-block terms and drops the within-true-cluster covariances across sub-blocks: \[\E[\hat{\mathbf M}_{\text{fine}}]=\mathbf M-\underbrace{\sum_c\sum_{s\neq s'\in c}\bX_s'\bD_{ss'}\bX_{s'}}_{\succeq\,\bzero\ \text{under positive within-correlation}}\ \prec\ \mathbf M,\] a downward bias — standard errors too small, spurious significance: the Moulton problem reintroduced.

Clustering too coarse (super-clusters that merge independent true clusters) adds pairs whose true covariance is zero, so in expectation it is unbiased, \[\E[\hat{\mathbf M}_{\text{coarse}}]=\mathbf M+\underbrace{\{\text{cross-true-cluster terms}\}}_{\E[\cdot]=\bzero}=\mathbf M,\] but it is built from fewer, larger, more unequal blocks — small \(C\), hence a noisy estimator and the finite-sample breakdown of the previous slide.

The trade-off is thus in the variance estimator, not in \(\bb\): finer buys precision at the risk of bias, coarser is unbiased but imprecise. This is why a changing standard error is not a decision rule — the change may be bias reduction (the fine level was dropping real covariances) or sampling noise (the coarse level has few clusters), and the magnitude alone cannot tell them apart. The correct level is a judgement about the design — sampling, treatment assignment — not something the data select for you.

5.3 Clustering in practice

Classic studies and their clustering level

The cluster level is an empirical judgement — here is how canonical studies made it.

Study Setting Cluster level Lesson
Moulton (1986, 1990) Aggregate regressors on micro data Aggregate group Conventional s.e. understated several-fold
Bertrand, Duflo & Mullainathan (2004) DiD, state\(\times\)year panel State Serial correlation inflates false significance
Angrist & Lavy (1999) Class size (Maimonides’ rule) Class / school Group-level instrument → cluster there
Petersen (2009) Finance panels Firm and year (two-way) Many published s.e. were wrong
Cameron, Gelbach & Miller (2008) Few clusters — Wild cluster bootstrap
Cameron, Gelbach & Miller (2011) Trade, networks Multiway Clustering on several dimensions
Abadie–Athey–Imbens–Wooldridge (2023) When to cluster Design level Clustering is about design, not just correlation

The common thread: cluster where the design induces dependence — the unit of sampling or of treatment assignment — and allow the broadest plausible correlation (in panels, serial correlation within the unit).

A menu of standard errors

Same coefficient, same regression — the standard error depends entirely on the dependence you allow for (\(\hat\beta=0.050\) throughout; numbers illustrative):

Standard error s.e. \(t\) verdict (5%)
i.i.d. (conventional) 0.010 5.0 reject
robust (White) 0.012 4.2 reject
cluster by state 0.028 1.8 fail to reject
two-way (state, year) 0.031 1.6 fail to reject

The point estimate never moves — only the inference does.

This is the format of Petersen’s (2009) finance table: reporting the same slope under i.i.d., White, one-way and two-way clustering exposed how many published standard errors were too small. In practice, report the clustered s.e. as the default, and note how far it exceeds the conventional one (the empirical Moulton factor). If significance hinges on the clustering choice, that fragility is itself a finding.

5.4 When robust is not enough

The limits of a robust variance

The Four Questions—What robust standard errors do and do not fix
  • They fix the variance — the sandwich (White or cluster) gives valid standard errors under heteroskedasticity or clustering.
  • They do not fix inconsistency — if \(\E[\mathbf{x}\varepsilon]\neq\bzero\) (endogeneity, Lecture 10), no variance formula saves you; the point estimate is wrong.
  • They do not fix dependence you failed to model — the White errors are wrong under clustering; you must match the sandwich to the actual dependence structure.

5.5 Summary

What to take away

Back to top