Module 2: Standard Errors & Clustering

How Sure Are We, Really?

Author

Austin Tucker

Published

July 23, 2026

1 The Sampling Distribution

Every regression table you will ever read has two columns that matter: the coefficient and its standard error. Most students learn to read the first column carefully and treat the second as decoration. This module is about the second column — because when the standard error is wrong, everything downstream (t-statistics, p-values, confidence intervals, policy conclusions) is wrong with it.

1.1 \(\hat{\beta}\) is a random variable

Suppose the population model for child hemoglobin (in g/dL) is:

\[ \text{hgb}_i = \beta_0 + \beta_1 \, \text{visits}_i + \beta_2 \, \text{age}_i + u_i \]

where visits counts home visits from a community health worker (CHW). In our simulated population of \(N = 10{,}000\) children, the true effect of each visit is \(\beta_1 = 0.6\) g/dL.

Here is the key conceptual move. When you draw a sample of \(n = 200\) children and run OLS, you get one number, say \(\hat{\beta}_1 = 0.57\). But that number depends on which 200 children happened to land in your sample. A different draw gives a different estimate. The estimator \(\hat{\beta}_1\) is therefore a random variable, and it has a distribution across hypothetical repeated samples — the sampling distribution.

What a standard error is

The standard error of \(\hat{\beta}_1\) is the standard deviation of the sampling distribution of \(\hat{\beta}_1\). Nothing more, nothing less.

In real research you get exactly one sample, so you never see this distribution. The reported SE is an estimate of its spread, computed from your single sample using a formula. Whether that formula is right is the entire subject of this module.

In a simulation we can do what real research cannot: draw the sample again and again, and watch the sampling distribution materialize.

1.2 Exercise 1: Watch the sampling distribution materialize

Draw 500 random samples of \(n = 200\) children each. In each sample, fit the regression of hemoglobin on visits and age, and store the coefficient on visits. Then plot the histogram of the 500 estimates — and note that the standard deviation of this histogram is the standard error.

Solution

set.seed(20260201)

beta1hat <- numeric(500)

for (i in 1:500) {
  samp <- pop[sample(nrow(pop), 200, replace = FALSE), ]
  mod  <- lm(hgb ~ visits + age_months, data = samp)
  beta1hat[i] <- coef(mod)["visits"]
}

cat("Mean of the 500 estimates:", round(mean(beta1hat), 3), "\n")
cat("True beta_1:               0.600\n")
cat("SD of the 500 estimates:  ", round(sd(beta1hat), 4), "\n")
cat("  ^ this SD *is* the standard error of beta_1-hat at n = 200\n")

ggplot(data.frame(estimate = beta1hat), aes(x = estimate)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 30, fill = slate_blue, color = "white", alpha = 0.8) +
  geom_vline(xintercept = mean(beta1hat), color = terracotta,
             linewidth = 1, linetype = "solid") +
  geom_vline(xintercept = 0.6, color = sage,
             linewidth = 1, linetype = "dashed") +
  labs(
    title = "Sampling Distribution of the Visits Coefficient",
    subtitle = "Terracotta = mean of estimates | Sage dashed = true value (0.6)",
    x = "Estimated coefficient on visits",
    y = "Density"
  ) +
  theme_house()

Interpretation question

The histogram is centered on 0.6, so the estimator is unbiased. If a single study reported \(\hat{\beta}_1 = 0.52\) with SE \(= 0.05\), where would that study sit in this histogram? Is it “wrong”?

Discussion: That study sits about 1.6 standard errors below the truth — entirely unremarkable. Roughly a third of well-conducted studies land more than one SE from the true value, by design of probability itself. This is why we never interpret a single point estimate as “the” effect: the honest statement is the estimate plus its uncertainty. Every real study is one draw from a histogram it can never see. The standard error is our estimate of how wide that invisible histogram is.

2 The Conventional SE Formula

2.1 From simulation to formula

We got the SE in Exercise 1 by brute force: 500 samples, take the SD. Real research gets one sample. The classical solution is a formula. Under the Gauss-Markov assumptions — in particular homoskedasticity, \(\text{Var}(u_i \mid X) = \sigma^2\) for all \(i\) — the variance of the OLS coefficient vector is:

\[ \text{Var}(\hat{\beta} \mid X) = \sigma^2 (X'X)^{-1} \]

We estimate \(\sigma^2\) with the residual variance \(\hat{\sigma}^2 = \frac{\sum_i \hat{u}_i^2}{n - k}\) (where \(k\) counts the regressors including the intercept), and the standard errors are the square roots of the diagonal of \(\hat{\sigma}^2 (X'X)^{-1}\). This is exactly what summary(lm()) prints.

For a single coefficient the formula has an illuminating scalar form:

\[ \text{Var}(\hat{\beta}_j) = \frac{\sigma^2}{\text{SST}_j \, (1 - R_j^2)} \]

where \(\text{SST}_j = \sum_i (x_{ij} - \bar{x}_j)^2\) and \(R_j^2\) is the \(R^2\) from regressing \(x_j\) on the other regressors.

What shrinks a standard error

Read the scalar formula like a clinician reads a differential. The SE shrinks when:

\(n\) grows — \(\text{SST}_j\) accumulates with every observation; SEs shrink at rate \(1/\sqrt{n}\). Quadrupling the sample halves the SE.
\(x_j\) varies more — if every child got exactly 3 visits, \(\text{SST}_j = 0\) and the effect of visits is not estimable at all. Variation in exposure is statistical oxygen.
\(\sigma^2\) is smaller — less unexplained noise in the outcome. Adding good predictors of hgb (like age) reduces \(\hat{\sigma}^2\).
\(x_j\) is less collinear with other regressors — \(R_j^2\) near 1 inflates the variance.

2.2 Exercise 2: Compute the conventional SE by hand

Draw one sample of 200 children and fit the model. Then rebuild the standard errors from raw matrix algebra — \(\hat{\sigma}^2 (X'X)^{-1}\) — and confirm they match summary(lm) to the last decimal.

Solution

set.seed(20260202)
samp <- pop[sample(nrow(pop), 200, replace = FALSE), ]

mod <- lm(hgb ~ visits + age_months, data = samp)

# Build the pieces of sigma2-hat * (X'X)^-1 by hand
X <- model.matrix(mod)         # 200 x 3: intercept, visits, age_months
u <- residuals(mod)
n <- nrow(X)
k <- ncol(X)

sigma2_hat <- sum(u^2) / (n - k)
V_conv     <- sigma2_hat * solve(t(X) %*% X)
se_hand    <- sqrt(diag(V_conv))

cat("Hand-computed SEs:\n")
print(round(se_hand, 5))
cat("\nsummary(lm) SEs:\n")
print(round(summary(mod)$coefficients[, "Std. Error"], 5))

Interpretation question

Your hand-computed SE on visits should match summary(lm) exactly. Compare it also to the SD of the histogram from Exercise 1. Why are those two numbers close but not identical?

Discussion: They estimate the same quantity from different directions. The Exercise 1 number is a Monte Carlo estimate of the true SE (limited by 500 replications); the formula number is a single-sample estimate (limited by having one \(\hat{\sigma}^2\) from one draw). Both converge to the same truth — provided the homoskedasticity assumption baked into the formula actually holds. When it does not, the formula quietly reports the wrong number while the simulation keeps telling the truth. That divergence is where we go next.

2.3 Exercise 3: What does quadrupling the sample buy you?

The \(1/\sqrt{n}\) law is the most budget-relevant fact in this course: precision is bought at a quadratically increasing price. Verify it by simulation — compare the sampling SD of \(\hat{\beta}_1\) at \(n = 200\) versus \(n = 800\) across 200 replications each.

Solution

set.seed(20260203)

b_small <- numeric(200)
b_big   <- numeric(200)

for (i in 1:200) {
  samp_s <- pop[sample(nrow(pop), 200, replace = FALSE), ]
  samp_b <- pop[sample(nrow(pop), 800, replace = FALSE), ]
  b_small[i] <- coef(lm(hgb ~ visits + age_months, data = samp_s))["visits"]
  b_big[i]   <- coef(lm(hgb ~ visits + age_months, data = samp_b))["visits"]
}

se_small <- sd(b_small)
se_big   <- sd(b_big)

cat("Empirical SE at n = 200:", round(se_small, 4), "\n")
cat("Empirical SE at n = 800:", round(se_big, 4), "\n")
cat("Ratio (should be near 2):", round(se_small / se_big, 2), "\n")

Interpretation question

Your study budget allows either \(n = 800\) children with a cheap hemoglobin test (adds noise to hgb) or \(n = 200\) with a gold-standard test. Which considerations from the scalar variance formula are in play?

Discussion: Both levers appear in \(\text{Var}(\hat{\beta}_j) = \sigma^2 / [\text{SST}_j (1 - R_j^2)]\). Quadrupling \(n\) quadruples \(\text{SST}_j\), halving the SE. The noisy test inflates \(\sigma^2\), pushing the other way. Which wins depends on how noisy the cheap test is — a question you can and should answer with exactly this kind of pre-study simulation. Field epidemiologists call this power calculation; it is nothing but the SE formula run in reverse.

3 Heteroskedasticity

3.1 When the error variance depends on \(x\)

The conventional formula leans on one load-bearing assumption: \(\text{Var}(u_i \mid X) = \sigma^2\), a single constant for everyone. Health data violate this constantly:

Health spending vs. income: poor households cluster tightly near subsistence-level spending; rich households range from frugal to catastrophic. The variance of spending fans out with income.
Length of stay vs. severity: mildly ill patients all go home in 1–2 days; severely ill patients range from 3 days to 3 months.
Any grouped average: if your outcome is a facility-level mean, facilities with fewer patients have noisier means.

When \(\text{Var}(u_i \mid x_i) = \sigma_i^2\) varies with \(i\), we have heteroskedasticity.

What heteroskedasticity does and does not do

\(\hat{\beta}\) remains unbiased and consistent. The point estimates are fine.
The conventional SEs are wrong — they can be too small or too large, so t-statistics and p-values are unreliable.
The fix does not change the coefficients at all; it changes only the uncertainty attached to them.

Let’s diagnose it visually first, in the household spending data (het_data), where true spending is \(400 + 18 \cdot \text{income} + \varepsilon\) with \(\text{sd}(\varepsilon) = 0.05 \cdot \text{income}^2\): the noise grows with the square of income.

The funnel shape is the classic signature: residual spread grows with income. No formal test needed — your eyes are a consistent estimator here.

3.2 The sandwich: robust SEs by hand

The fix, due to White (1980), replaces the single \(\sigma^2\) with the squared residuals themselves. In matrix form:

\[ \widehat{\text{Var}}_{\text{HC0}}(\hat{\beta}) = \underbrace{(X'X)^{-1}}_{\text{bread}} \underbrace{\left( \sum_{i=1}^{n} \hat{u}_i^2 \, x_i x_i' \right)}_{\text{meat}} \underbrace{(X'X)^{-1}}_{\text{bread}} \]

Bread, meat, bread — hence “sandwich estimator.” Each observation contributes its own squared residual \(\hat{u}_i^2\) to the meat, so observations in the high-variance region correctly count for more uncertainty. The HC1 variant multiplies by the small-sample correction \(\frac{n}{n-k}\) (this is what Stata’s , robust reports).

Notice the payoff of the sandwich structure: under homoskedasticity, \(\hat{u}_i^2 \approx \sigma^2\) for everyone, the meat collapses to \(\sigma^2 X'X\), and the whole sandwich reduces to the conventional \(\sigma^2 (X'X)^{-1}\). Robust SEs are not a different philosophy — they are the same formula with one fewer assumption.

3.3 Exercise 4: Compute robust SEs by hand

Fit the spending regression, then compute HC0 and HC1 robust standard errors from the sandwich formula and compare them with the conventional SEs. No packages — three lines of matrix algebra.

Solution

mod_het <- lm(spend ~ income, data = het_data)

X <- model.matrix(mod_het)
u <- residuals(mod_het)
n <- nrow(X)
k <- ncol(X)

# Bread: (X'X)^-1
bread <- solve(t(X) %*% X)

# Meat: sum of u_i^2 * x_i x_i'  (note: X * u^2 scales each row i by u_i^2)
meat <- t(X) %*% (X * u^2)

# Sandwich
V_hc0 <- bread %*% meat %*% bread
V_hc1 <- (n / (n - k)) * V_hc0

cat("Conventional SEs:\n")
print(round(summary(mod_het)$coefficients[, "Std. Error"], 3))
cat("\nRobust HC0 SEs:\n")
print(round(sqrt(diag(V_hc0)), 3))
cat("\nRobust HC1 SEs (Stata's ', robust'):\n")
print(round(sqrt(diag(V_hc1)), 3))

Interpretation question

The robust SE on income should come out about 20–25% larger than the conventional SE. Why larger here — and could robust SEs ever come out smaller?

Discussion: Here the high-variance observations (rich households) are also the high-leverage observations (far from mean income), so the conventional formula — which averages the noise evenly over everyone — understates the uncertainty about the slope. Robust SEs can come out smaller when the variance is highest near the center of the \(x\) distribution and low in the tails. The direction is an empirical question; that is exactly why we let the data speak through \(\hat{u}_i^2\) rather than assuming a single \(\sigma^2\). And note your coefficients did not move — only the honesty of the uncertainty did.

4 Clustering

4.1 The big one for field research

Heteroskedasticity is a misdemeanor; clustering is a felony. Nearly every health policy study has structure like this:

Children within health facilities — our running example
Patients within physicians within hospitals
Villagers within villages in a cluster-randomized trial
Residents within states in a difference-in-differences design

And, critically, the treatment is assigned at the group level: the facility adopts the CHW program, the village gets the intervention, the state passes the law. Every child in a treated facility has the same value of treat.

Children in the same facility share things the model does not capture: the same catchment-area diet, the same water source, the same nurse who is either meticulous or burned out. Formally, the error has a shared component:

\[ \text{hgb}_{ig} = \beta_0 + \beta_1 \, \text{treat}_g + \underbrace{f_g + e_{ig}}_{u_{ig}}, \qquad f_g \sim N(0, \sigma_f^2), \quad e_{ig} \sim N(0, \sigma_e^2) \]

Two children in the same facility \(g\) share \(f_g\), so their errors are correlated:

\[ \rho = \text{Corr}(u_{ig}, u_{jg}) = \frac{\sigma_f^2}{\sigma_f^2 + \sigma_e^2} \]

This \(\rho\) is the intra-cluster correlation (ICC). Our clustered population (fac_data) has 40 facilities of 50 children each, with \(\sigma_f = 0.6\) and \(\sigma_e = 1.0\).

Look at the large points: entire facilities sit above or below the pack, dragging all 50 of their children with them. That vertical scatter of facility means is \(f_g\) made visible.

4.2 Why independence fails — and what it costs

The conventional formula credits you with \(n\) independent pieces of information. But the 50 children in a facility partially echo one another: once you know the facility’s shared shock, each additional child adds less news. The design effect quantifies the loss. For a regressor constant within clusters of size \(m\):

\[ \text{DEFF} = 1 + (m - 1)\,\rho, \qquad n_{\text{eff}} = \frac{n}{\text{DEFF}} \]

The brutal arithmetic of clustering

Even a small ICC is devastating when clusters are large, because \(\rho\) is multiplied by \(m - 1\). With \(\rho = 0.05\) and \(m = 100\), DEFF \(= 5.95\): your 4,000 patients carry the information of 672. For facility-level treatments, the effective sample size is governed by the number of clusters, not the number of individuals. Forty facilities is closer to \(n = 40\) than to \(n = 2{,}000\).

4.3 Exercise 5: ICC, design effect, effective sample size

Compute the ICC implied by our DGP (\(\sigma_f = 0.6\), \(\sigma_e = 1.0\)), the design effect for clusters of \(m = 50\), and the effective sample size of our \(n = 2{,}000\) children.

Solution

sigma_f <- 0.6   # SD of the facility-level shock
sigma_e <- 1.0   # SD of the child-level error
m       <- 50    # children per facility
n       <- 2000  # total children

# Intra-cluster correlation
rho <- sigma_f^2 / (sigma_f^2 + sigma_e^2)

# Design effect
deff <- 1 + (m - 1) * rho

# Effective sample size
n_eff <- n / deff

cat("ICC (rho):            ", round(rho, 3), "\n")
cat("Design effect (DEFF): ", round(deff, 2), "\n")
cat("Nominal sample size:  ", n, "children\n")
cat("Effective sample size:", round(n_eff), "children\n")
cat("Number of facilities:  40\n")

Interpretation question

You should find \(\rho \approx 0.26\), DEFF \(\approx 14\), and an effective sample size of roughly 143 — barely more than the 40 facilities. If you could afford 2,000 more child observations, would you rather add 50 facilities of 40 children or 40 more children to each existing facility?

Discussion: Add facilities, almost always. Enrolling more children per facility mostly re-measures the same facility shocks \(f_g\) — information you already have. New facilities contribute new independent draws of \(f_g\), which is the binding constraint on precision when treatment varies at the facility level. This is why cluster-randomized trials are powered by the number of clusters, and why a trial with 6 enormous hospitals can be less informative than one with 30 small clinics.

4.4 The cluster-robust sandwich

The repair generalizes the White sandwich: instead of letting each observation contribute its own squared residual, we let each cluster contribute the outer product of its summed score. With clusters \(g = 1, \dots, G\):

\[ \widehat{\text{Var}}_{\text{CR}}(\hat{\beta}) = (X'X)^{-1} \left( \sum_{g=1}^{G} X_g' \hat{u}_g \, \hat{u}_g' X_g \right) (X'X)^{-1} \]

where \(X_g\) and \(\hat{u}_g\) stack the rows for cluster \(g\). Summing residuals within a cluster before squaring is what captures the within-cluster correlation: if all 50 residuals in a facility are pushed up together by \(f_g\), their sum is large, and the meat records it. The CR1 small-sample correction multiplies by \(\frac{G}{G-1} \cdot \frac{n-1}{n-k}\) (this is Stata’s , cluster()).

4.5 Exercise 6: Cluster-robust SEs by hand

Fit the treatment-effect regression on the clustered data and build the cluster-robust SE from scratch. The one new tool: rowsum(X * u, group = fac_id) computes \(X_g'\hat{u}_g\) for every facility at once, returning a \(G \times k\) matrix of cluster score sums.

Solution

mod_fac <- lm(hgb_c ~ treat, data = fac_data)

X <- model.matrix(mod_fac)
u <- residuals(mod_fac)
n <- nrow(X)
k <- ncol(X)
G <- length(unique(fac_data$fac_id))

# Bread
bread <- solve(t(X) %*% X)

# Meat: sum over clusters of (X_g' u_g)(X_g' u_g)'
Xu   <- rowsum(X * u, group = fac_data$fac_id)   # G x k matrix
meat <- t(Xu) %*% Xu

# CR1 small-sample correction
adj <- (G / (G - 1)) * ((n - 1) / (n - k))

V_cr   <- adj * bread %*% meat %*% bread
se_cr  <- sqrt(diag(V_cr))

cat("Coefficient on treat:  ", round(coef(mod_fac)["treat"], 3), "\n\n")
cat("Conventional SE:       ",
    round(summary(mod_fac)$coefficients["treat", "Std. Error"], 3), "\n")
cat("Cluster-robust SE:     ", round(se_cr["treat"], 3), "\n")
cat("Ratio (cluster/conv):  ",
    round(se_cr["treat"] /
          summary(mod_fac)$coefficients["treat", "Std. Error"], 2), "\n")

Interpretation question

The cluster-robust SE should be roughly 3–4 times the conventional SE — close to \(\sqrt{\text{DEFF}} = \sqrt{14} \approx 3.7\) from Exercise 5. The treatment effect that looked overwhelmingly significant now hangs near the threshold. Which analysis would you defend to a journal referee, and why?

Discussion: The clustered one, without hesitation. The conventional SE claims 2,000 independent observations; Exercise 5 showed the design carries about 143 observations’ worth of information about a facility-level treatment. The sandwich did not “destroy” significance — it revealed that the conventional analysis was manufacturing certainty out of correlated noise. A t-statistic near 2 with honest SEs beats a t-statistic near 8 with dishonest ones. Referees know the DEFF arithmetic too.

5 The Over-Rejection Simulation

5.1 The centerpiece: how often do we cry wolf?

Here is the single most persuasive demonstration in applied econometrics, in the spirit of Bertrand, Duflo & Mullainathan (2004). We assign a placebo treatment — the true effect is exactly zero — at the facility level, and ask: how often does each SE lead us to declare a significant effect at \(\alpha = 0.05\)?

A valid test rejects a true null 5% of the time. That 5% is a promise. Let’s audit it.

The plan, repeated 500 times:

Generate a fresh clustered population: 40 facilities \(\times\) 25 children, facility shocks \(f_g \sim N(0, 0.6)\), child errors \(e \sim N(0, 1)\).
Randomly assign a placebo “program” to 20 of the 40 facilities. The DGP for hemoglobin does not include the treatment: the true effect is 0.
Regress hgb on treat; test \(H_0: \beta_1 = 0\) at \(\alpha = 0.05\) twice — once with conventional SEs, once with hand-coded cluster-robust SEs (the cluster_se() function from setup, identical to your Exercise 6 code).

5.2 Exercise 7: Run the placebo audit

Solution

set.seed(20260204)

n_sim   <- 500
n_fac_s <- 40
m_s     <- 25
fid     <- rep(1:n_fac_s, each = m_s)

reject_conv <- logical(n_sim)
reject_clus <- logical(n_sim)

for (s in 1:n_sim) {
  # Fresh clustered data; note: NO treatment effect in the DGP
  f_g   <- rnorm(n_fac_s, mean = 0, sd = 0.6)
  e     <- rnorm(n_fac_s * m_s, mean = 0, sd = 1.0)
  trt   <- as.numeric(fid %in% sample(1:n_fac_s, 20))  # placebo: 20 facilities
  hgb_s <- 10.2 + 0 * trt + f_g[fid] + e               # true effect = 0

  mod_s   <- lm(hgb_s ~ trt)
  b       <- coef(mod_s)["trt"]
  se_conv <- summary(mod_s)$coefficients["trt", "Std. Error"]
  se_clus <- cluster_se(mod_s, fid)["trt"]

  reject_conv[s] <- abs(b / se_conv) > 1.96
  reject_clus[s] <- abs(b / se_clus) > 1.96
}

rates <- data.frame(
  se_type = factor(c("Conventional", "Cluster-robust"),
                   levels = c("Conventional", "Cluster-robust")),
  rate    = c(mean(reject_conv), mean(reject_clus))
)

cat("Rejection rate of a TRUE null at alpha = 0.05:\n")
cat("  Conventional SEs:  ", round(mean(reject_conv), 3), "\n")
cat("  Cluster-robust SEs:", round(mean(reject_clus), 3), "\n")
cat("  Nominal (promised):  0.050\n")

ggplot(rates, aes(x = se_type, y = rate, fill = se_type)) +
  geom_col(width = 0.55, alpha = 0.9) +
  geom_hline(yintercept = 0.05, color = charcoal,
             linetype = "dashed", linewidth = 0.8) +
  geom_text(aes(label = sprintf("%.1f%%", 100 * rate)),
            vjust = -0.5, color = charcoal, size = 4.5) +
  scale_fill_manual(values = c("Conventional" = terracotta,
                               "Cluster-robust" = sage)) +
  scale_y_continuous(limits = c(0, 0.65)) +
  labs(
    title = "False-Positive Rates Under a Placebo Facility-Level Treatment",
    subtitle = "True effect = 0 | Dashed line = the promised 5% | 500 simulations",
    x = NULL,
    y = "Share of simulations rejecting the (true) null"
  ) +
  theme_house() +
  theme(legend.position = "none")

Interpretation question

You should see the conventional SEs rejecting a true null roughly 40–55% of the time — closer to a coin flip than to the promised 5% — while the cluster-robust SEs sit near 5–8%. Imagine 100 published evaluations of facility-level programs that truly do nothing, all analyzed with conventional SEs. What does the literature look like?

Discussion: Roughly half of those 100 null programs get published as “statistically significant” successes — a literature full of confident findings about interventions that do nothing. This is not hypothetical: Bertrand, Duflo & Mullainathan (2004) showed that placebo state-level “laws” in difference-in-differences designs were declared significant up to 45% of the time with unclustered SEs. The point estimates in all these studies are unbiased; it is purely the uncertainty that is misstated. The cluster-robust rate of 5–8% is slightly above nominal because 40 clusters is adequate but not luxurious — with fewer clusters it degrades further, which brings us to the practical guidance.

6 Practical Guidance

6.1 Cluster at the level of treatment assignment

The default rule, and the one referees will hold you to:

The clustering rule of thumb

Cluster your standard errors at the level at which treatment is assigned (or, if treatment is assigned below the level of correlated shocks, at the level of the shocks).

Table 1: Where to cluster

Design	Treatment varies at	Cluster at
CHW program adopted by facilities	Facility	Facility
Cluster-randomized trial in villages	Village	Village
State insurance expansion (DiD)	State	State
Individually randomized drug trial	Patient	Usually none needed

Clustering at too fine a level (e.g., patient, when the program varies by facility) recreates the over-rejection you just simulated. Clustering coarser than needed costs some power but rarely lies. When in doubt, cluster coarser — and say so.

6.2 The few-clusters caveat

Cluster-robust SEs are justified asymptotically in the number of clusters \(G\), not the number of observations. With \(G\) below roughly 40, the CR1 sandwich is itself too small and over-rejection creeps back — our simulation’s 5–8% instead of 5% was a mild preview. With \(G\) of 10 states or 12 hospitals, it can be severe. The standard remedy is the wild cluster bootstrap (Cameron, Gelman-style discussions aside, see Cameron, Gelbach & Miller 2008), which resamples cluster-level residuals with random sign flips. We will not implement it here — just recognize the situation: few clusters + cluster-level treatment = do not trust the plain sandwich, reach for the wild bootstrap.

6.3 Reading SEs in published tables

A field guide for the tables you will read in health policy journals:

Find the footnote first. “Robust standard errors in parentheses” vs. “standard errors clustered by facility (40 clusters)” are entirely different claims. If a facility-level intervention reports only heteroskedasticity-robust SEs, mentally multiply the SEs by \(\sqrt{\text{DEFF}}\) — significance often evaporates.
Count the clusters, not the \(n\). A table boasting \(n = 60{,}000\) patients across 8 hospital systems has 8 effective units for any system-level policy variable.
Stars are not effect sizes. A precisely estimated tiny effect and a noisily estimated large one can carry the same stars. Read the coefficient, the SE, and the implied confidence interval — in clinical units (here, g/dL of hemoglobin) — before reading the asterisks.
Be suspicious of t-statistics above 8 on group-level treatments with individual-level data. You now know exactly how those are manufactured.

6.4 Summary

Table 2: Choosing a standard error

Situation	Conventional SE	Correct SE	Coefficients affected?
Homoskedastic, independent errors	Correct	Same	—
Heteroskedastic errors	Wrong (either direction)	HC1 robust sandwich	No
Clustered errors, group-level treatment	Far too small	Cluster-robust (CR1) sandwich	No
Clustered errors, \(G \lesssim 40\)	Far too small	Wild cluster bootstrap	No

What’s next?

You now know what a standard error is (the SD of the sampling distribution) and how to keep it honest under heteroskedasticity and clustering. Module 3 (Inference & Asymptotics) takes the next step: why the ratio \(\hat{\beta}/\text{SE}\) is compared to 1.96 in the first place, what the Central Limit Theorem does and does not guarantee in finite samples, and how confidence intervals and hypothesis tests inherit — or fail to inherit — the honesty of the standard errors underneath them.

References

Angrist, Joshua D., and Jorn-Steffen Pischke. 2009. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton University Press. (Chapter 8: “Nonstandard Standard Error Issues.”)
Bertrand, Marianne, Esther Duflo, and Sendhil Mullainathan. 2004. “How Much Should We Trust Differences-in-Differences Estimates?” Quarterly Journal of Economics 119 (1): 249–275.
Cameron, A. Colin, Jonah B. Gelbach, and Douglas L. Miller. 2008. “Bootstrap-Based Improvements for Inference with Clustered Errors.” Review of Economics and Statistics 90 (3): 414–427.
Cameron, A. Colin, and Douglas L. Miller. 2015. “A Practitioner’s Guide to Cluster-Robust Inference.” Journal of Human Resources 50 (2): 317–372.
Abadie, Alberto, Susan Athey, Guido W. Imbens, and Jeffrey M. Wooldridge. 2023. “When Should You Adjust Standard Errors for Clustering?” Quarterly Journal of Economics 138 (1): 1–35.
White, Halbert. 1980. “A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity.” Econometrica 48 (4): 817–838.
Wooldridge, Jeffrey M. 2020. Introductory Econometrics: A Modern Approach. 7th ed. Cengage Learning.