# A Simple Formula for Clustered Standard Errors for a Univariate Regression with a Dummy Regressor

Suppose I have a simple cross-sectional regression like:

$$y_i = \alpha + \beta \cdot D_i + \varepsilon_i,$$

where $D_i$ is a binary treatment flag and $i$ indexes individuals. $\beta$ is the parameter of interest. We can estimates $$SE(\hat \beta)=\frac{\hat \sigma}{\sqrt{\sum_i^N(D_i-\bar D)^2}}.$$

Now to complicate things, let's say that we have $G$ groups/blocks and that everyone in the same group receives an identical treatment assignment. To adjust for this, we typically cluster standard errors at group-level (assuming there are enough groups for the asymptotics to kick in), which is usually accomplished (switching to matrix/vector notation below) with $$Var(\hat \beta)=(X'X)^{-1}\left(\sum_{g=1}^{G} X_g' \hat \varepsilon_g'\hat\varepsilon_g X_g\right)(X'X)^{-1}.$$

Does there exist an non-matrix equivalent of the formula for the clustered SEs case like there is in iid case with the simple univariate regression?

• Intuitively you should arrive at some sort of a heteroskedatstic t-test. – StasK Aug 7 '18 at 13:15

OK, let's parse this term by term.

$$\mathbf{x}_{gi} = (1, D_{gi})$$

$$(X'X) = \sum_g \sum_i (1, D_{gi})' (1, D_{gi}) = n \left( \begin{array}{cc} 1 & p \\ p & p \end{array} \right)$$

where $$p=\sum_{gi} D_{gi} / n$$ = proportion of 1s in the sample.

$$(X'X)^{-1} = \frac 1n \left( \begin{array}{cc} p & -p \\ -p & 1 \end{array} \right)$$

Within cluster sum of scores: $$\mathbf{X}_g \mathbf{\epsilon}_g = \sum_i (1, D_{gi})' \epsilon_{gi} = n_g (\bar \epsilon_g, \bar\epsilon_g^I), \quad \bar \epsilon_g = \frac1{n_g} \sum_i \epsilon_{gi}, \bar\epsilon_g^I = \frac1{n_g} \sum_i D_{gi} \epsilon_{gi} = \frac{\sum_{i: D_{gi}=1} \epsilon_{gi}}{\sum_{i: D_{gi}=1} 1}$$

Within cluster outer product of scores: $$\mathbf{X}_g' \mathbf{\epsilon}_g' \mathbf{\epsilon}_g \mathbf{X}_g = n_g^2 (\bar \epsilon_g, \bar\epsilon_g^I)' (\bar \epsilon_g, \bar\epsilon_g^I) = n_g^2 \left( \begin{array}{cc} v_g^A & v_g^{AI} \\ v_g^{AI} & v_g^I \end{array} \right)$$

where $$A$$ stands for "all in cluster", $$I$$ stands for "those in cluster with $$D_{gi}=1$$, and $$AI$$ is a mix:

$$v_g^A = \frac1{n_g^2} \sum_i \epsilon_{gi}^2 + \frac1{n_g^2} \sum_{i\neq j} \epsilon_{gi} \epsilon_{gj}, \, v_g^I = \frac1{n_g^2} \sum_{i: D_{gi}=1} \epsilon_{gi}^2 + \frac1{n_g^2} \sum_{\begin{array}{c}i\neq j \\ D_{gi}=1 \\ D_{gj}=1 \end{array}} \epsilon_{gi} \epsilon_{gj},$$

$$v_g^{AI} = v_g^I + \frac1{n_g^2} \sum_{\begin{array}{c}D_{gi}=0 \\ D_{gj}=1 \end{array}} \epsilon_{gi} \epsilon_{gj}$$

And then you sum these monsters over clusters.

$$\sum_g \mathbf{X}_g' \mathbf{\epsilon}_g' \mathbf{\epsilon}_g \mathbf{X}_g= \sum_g n_g^2 \left( \begin{array}{cc} v_g^A & v_g^{AI} \\ v_g^{AI} & v_g^I \end{array} \right)$$

Frankly, I am not seeing how this would improve the visual understanding of the issue. You could start making some simplifying assumptions about $$\mathbb{V}\epsilon_{gi}$$ when $$D_{gi}=1$$ vs when $$D_{gi}=0$$, and about within-cluster correlations, and about proportions of the $$D$$ within clusters (the whole cluster treated, or equal in all clusters, or a Beta of some kind arcoss clusters), and maybe you can derive heteroskedastic ANOVA variance decompositions, or improved effective degrees of freedom, or a nice saddlepoint approximation :). That's how the heteroskedasticity corrections like HC2 were derived: assume homoskedasticity, and make sure that the heteroskedasticty-robust estimator is, at least, unbiased in that case.

Update: I really missed the point in the OP that everybody in the cluster has the same $$D_{gi}$$. Let's see how this would simplify things. I don't need a superscript over the sums of squared residuals, because the sums are always over the whole cluster, $$v_g \equiv v_g^A$$.

$$\mathbf{X}_g' \mathbf{\epsilon}_g' \mathbf{\epsilon}_g \mathbf{X}_g = n_g^2 (\bar \epsilon_g, \bar\epsilon_g^I)' (\bar \epsilon_g, \bar\epsilon_g^I) = n_g^2 \left( \begin{array}{cc} v_g & v_g \\ v_g & v_g \end{array} \right), \quad D_g = 1$$

$$\mathbf{X}_g' \mathbf{\epsilon}_g' \mathbf{\epsilon}_g \mathbf{X}_g = n_g^2 (\bar \epsilon_g, \bar\epsilon_g^I)' (\bar \epsilon_g, \bar\epsilon_g^I) = n_g^2 \left( \begin{array}{cc} v_g & 0 \\ 0 & 0 \end{array} \right), \quad D_g = 0$$

$$\sum_g \mathbf{X}_g' \mathbf{\epsilon}_g' \mathbf{\epsilon}_g \mathbf{X}_g = \left( \begin{array}{cc} \sum_g n_g^2 v_g & \sum_g^I n_g^2 v_g \\ \sum_g^I n_g^2 v_g & \sum_g^I n_g^2 v_g \end{array} \right) = \left( \begin{array}{cc} \sum_g^I n_g^2 v_g + \sum_g^O n_g^2 v_g & \sum_g^I n_g^2 v_g \\ \sum_g^I n_g^2 v_g & \sum_g^I n_g^2 v_g \end{array} \right)$$

where $$\sum_g^I \equiv \sum_{g: D_g = 1}$$ are summations only over the "treated" clusters where $$D_g=1$$; and $$\sum_g^O \equiv \sum_{g: D_g = 0}$$ are summations only over the "comparison" clusters where $$D_g=0$$.

Considering that this is simple enough, we can multiply by $$(X'X)^{-1}$$ derived above:

$$\frac1n \left( \begin{array}{cc} p & -p \\ -p & 1 \end{array} \right) \times \left( \begin{array}{cc} \sum_g^I n_g^2 v_g + \sum_g^O n_g^2 v_g & \sum_g^I n_g^2 v_g \\ \sum_g^I n_g^2 v_g & \sum_g^I n_g^2 v_g \end{array} \right) =$$

$$= \frac1n \left( \begin{array}{cc} p \sum_g^O n_g^2 v_g & 0 \\ (1-p) \sum_g^I n_g^2 v_g - p \sum_g^O n_g^2 v_g & (1-p) \sum_g^I n_g^2 v_g \end{array} \right);$$

The whole thing = $$\frac1n \left( \begin{array}{cc} p \sum_g^O n_g^2 v_g & 0 \\ (1-p) \sum_g^I n_g^2 v_g - p \sum_g^O n_g^2 v_g & (1-p) \sum_g^I n_g^2 v_g \end{array} \right) \times \frac1n \left( \begin{array}{cc} p & -p \\ -p & 1 \end{array} \right) =$$

$$= \frac1{n^2} \left( \begin{array}{cc} p^2 \sum_g^O n_g^2 v_g & - p^2 \sum_g^O n_g^2 v_g \\ - p^2 \sum_g^O n_g^2 v_g & p^2 \sum_g^O n_g^2 v_g + (1-p)^2 \sum_g^I n_g^2 v_g \end{array} \right)$$

For variance of residuals $$\mathbb{V}\epsilon_{gi} = \sigma^2_g$$ and the intraclass correlation $$\frac{1}{\sigma^2_g} \mathbb{E} \epsilon_{gi} \epsilon_{gj} = \rho_g$$, we can have $$\mathbb{E} v_g = \frac1{n_g^2} \sum_i \mathbb{E} \epsilon_{gi}^2 + \frac1{n_g^2} \sum_{i\neq j} \mathbb{E} \epsilon_{gi} \epsilon_{gj} = \frac1{n_g} \sigma^2_g + \frac{n_g - 1}{n_g} \rho_g \sigma^2_g = \frac1{n_g} \sigma^2_g \left[ 1 + (n_g-1) \rho_g \right]$$

The square brackets is the design effect for cluster samples that survey statisticians have known since 1950s... definitely since Kish (1965) book. It describes by how much a naive variance $$\sigma^2/n_g$$ would have to go up because of the cluster effects.

This is getting more useful, as you can explore how this stuff behaves depending on the relative size of clusters, relative variances of clusters, etc. Intra-class correlations are notoriously difficult to estimate, so keeping them cluster specific $$\rho_g$$ makes little sense, you would be better off with some common $$\rho_g \equiv \rho$$. I have provided a brief overview of what I had known in 2014 about ICCs as a part of this chapter. Note also that the estimator should be presumed biased; Korn and Graubard (1999) sampling statistics book has a full set of exercises (2.3-6 through 2.3-11) showing that biases can go each and every way depending on the interplay of whether the big clusters have variances larger or smaller than those in small clusters.

This paper by Cameron & Miller has a little bit of that in it.

If I understand you correctly, you're asking for a closed form solution, not in matrix notation? I think this would be:

$$V_{clu}(\hat \beta)=\frac{\sum_i\sum_j\tilde{D}_i\tilde{D}_jE[u_iu_j]1(i, j \text{ in same cluster)}}{\left(\sum_i\tilde{D_i}^2\right)^2}$$ where, $\tilde{D}_i=D_i-\bar{D}$ and $1(i, j \text{ in same cluster)}$ is an indicator function that takes the value 1 iff observation $i$ and $j$ are in the same cluster.

The corresponding estimator: $$V_{clu}(\hat \beta)=\frac{\sum_i\sum_j\tilde{D}_i\tilde{D}_j\hat{u_i}\hat{u_j}1(i, j \text{ in same cluster)}}{\left(\sum_i\tilde{D_i}^2\right)^2}$$

This is still not as "simple" as the one you put up there for the homescedastic case. But I guess that's simply because the problem is more complex, and probably that's the reaons matrix notation is more popular here.