Few Clusters and Sandwich Inference Failure

Why the Number of Clusters Matters

The cluster-robust sandwich estimates uncertainty by combining contributions from independent clusters. Its large-sample justification therefore relies on having many independent clusters, not simply many observations. A study with 2,000 measurements from 12 clinics still has only 12 independent contributions to the sandwich “meat.” With so few contributions, the estimated meat can be too small, and the usual normal approximation for a test statistic can be poor.

The downward bias has a concrete source: the mean model is fitted using the same data that produce the residual cluster scores. Fitting constrains those scores, especially in clusters with high leverage. The result is that the empirical variability of the fitted scores can understate the variability of the scores that would be seen across repeated samples.

The Cluster Sandwich and Its Finite-Sample Limitation

Definition: For \(K\) independent clusters, let \(U_i(\beta)\) be the estimating-function contribution from cluster \(i\), and let \(\hat\beta\) solve \(\sum_{i=1}^K U_i(\hat\beta)=0\). Define the estimated bread \(\hat A=-\sum_{i=1}^K \partial U_i(\hat\beta)/\partial\beta^\mathsf{T}\) and the empirical meat \(\hat M=\sum_{i=1}^K U_i(\hat\beta)U_i(\hat\beta)^\mathsf{T}\). The ordinary cluster-robust sandwich covariance estimate, often called CR0, is \(\hat V=\hat A^{-1}\hat M\hat A^{-\mathsf{T}}\), when \(\hat A\) is invertible.

The sandwich targets the variability of a sum of independent cluster contributions. Under suitable regularity conditions—including a correctly specified marginal mean, independent clusters, a nonsingular limiting bread, finite moments, and no single cluster dominating the sum—this sum can be approximated by a normal distribution as \(K\) grows. With few clusters, these approximations may be inaccurate even if the estimating equation has a sensible solution.

A useful way to see the problem is to examine ordinary least squares with independent clusters. This is a special case of estimating equations, and its finite-sample calculation makes the role of fitted residuals explicit.

Theorem: Suppose \(y=X\beta+\varepsilon\), where \(X\) is a fixed, full-column-rank \(N\times p\) design matrix, and \(\varepsilon\sim N(0,\sigma^2 I_N)\). Partition the rows and observations into \(K\) clusters, with design block \(X_i\) and residual block \(e_i\). Let \(H=X(X^\mathsf{T}X)^{-1}X^\mathsf{T}\), and let \(H_{ii}\) denote its diagonal block for cluster \(i\). The expected CR0 meat satisfies \[ E[\hat M\mid X]=\sigma^2\sum_{i=1}^K X_i^\mathsf{T}(I-H)_{ii}X_i \preceq \sigma^2 X^\mathsf{T}X. \] Consequently, the expected CR0 covariance is no larger in positive-semidefinite order than the true covariance \(\sigma^2(X^\mathsf{T}X)^{-1}\).

Why Fitting Shrinks the Empirical Meat

Proof. The least-squares residual vector is \(e=(I-H)\varepsilon\). Since \(H\) is a symmetric idempotent projection, \(I-H\) is also symmetric and idempotent. Therefore

E[ee^\mathsf{T}\mid X] =(I-H)\sigma^2 I_N(I-H) =\sigma^2(I-H).

The cluster score is \(X_i^\mathsf{T}e_i\), so its expected outer product is

E[X_i^\mathsf{T}e_i e_i^\mathsf{T}X_i\mid X] =\sigma^2X_i^\mathsf{T}(I-H)_{ii}X_i.

Summing over clusters gives the stated expression. Also, \(X^\mathsf{T}X=\sum_i X_i^\mathsf{T}X_i\), so

\sigma^2X^\mathsf{T}X-E[\hat M\mid X] =\sigma^2\sum_i X_i^\mathsf{T}H_{ii}X_i\succeq 0.

Each \(H_{ii}\) is positive semidefinite because it is a principal block of the positive-semidefinite projection matrix \(H\). The CR0 covariance is \((X^\mathsf{T}X)^{-1}\hat M(X^\mathsf{T}X)^{-1}\), while the true covariance is \(\sigma^2(X^\mathsf{T}X)^{-1}\). Multiplying the meat inequality on both sides by \((X^\mathsf{T}X)^{-1}\) preserves the positive-semidefinite ordering. This proves the result. \(\square\)

This theorem is an exact result under its stated linear-model assumptions, not a universal theorem that every sandwich estimate is downward biased. In more general GEE settings, the direction and size of finite-sample bias can depend on the mean model, working covariance, and cluster structure. The result does establish a common mechanism: estimating the mean from the data reduces the apparent variability of the fitted cluster scores.

Worked Example: Four Clusters and an Intercept

Suppose four independent cluster-level observations satisfy \(Y_i=\beta+\varepsilon_i\), with \(\varepsilon_i\) independent and distributed \(N(0,\sigma^2)\). The estimate is \(\hat\beta=\bar Y\), and the CR0 estimate of its variance is \(\sum_i(Y_i-\bar Y)^2/4^2\).

The residuals sum to zero, so only three independent residual directions remain. More precisely, \(\sum_i(Y_i-\bar Y)^2/\sigma^2\) has a chi-squared distribution with \(4-1=3\) degrees of freedom, whose expectation is 3. Thus

E[\hat V_{\mathrm{CR0}}] =\frac{3\sigma^2}{16}, \qquad \operatorname{Var}(\bar Y)=\frac{\sigma^2}{4} =\frac{4\sigma^2}{16}.

The expected CR0 variance is \(3/4\) of the true variance. For a numerical check, take observed values \(1,2,3,4\). Their mean is \(2.5\), and the residuals are \(-1.5,-0.5,0.5,1.5\), with squared sum \(5\). CR0 is \(5/16=0.3125\). The usual sample-variance estimate of the variance of the mean is \((5/3)/4=5/12\), larger by the factor \(4/3\). The observed comparison illustrates the calculation; the expectation argument above establishes the bias under the model.

Leverage and Unequal Cluster Sizes

Leverage describes how strongly a cluster’s observations influence the fitted mean parameters. In the linear-model theorem, the block \(H_{ii}\) measures the cluster’s share of the projection onto the fitted model. When that block is large in a direction important for estimating a coefficient, fitting the model can leave little residual variation in that direction. The corresponding cluster score then contributes less to the empirical meat than its unsampled variability would warrant.

Worked Example: High Leverage in a Three-Cluster Regression

Consider three independent Gaussian observations, one per cluster, in a regression with an intercept and predictor values \(-1,0,1\). The design rows are \((1,-1),(1,0),(1,1)\). The diagonal leverages are \(h_i=1/3+x_i^2/2\), giving \(5/6,1/3,5/6\). For one-observation clusters, the theorem gives the expected meat as \(\sigma^2\sum_i(1-h_i)x_i x_i^\mathsf{T}\).

The first and third clusters each have weight \(1-h_i=1/6\); the middle cluster has weight \(2/3\). Adding their weighted outer products gives

E[\hat M\mid X] =\sigma^2 \begin{pmatrix} 1 & 0\\ 0 & 1/3 \end{pmatrix}.

By contrast, the target meat is \(\sigma^2X^\mathsf{T}X=\sigma^2\begin{pmatrix}3&0\\0&2\end{pmatrix}\). The expected meat for the slope is only one-sixth of its target value. Here the two extreme predictor values have leverage \(5/6\), and the slope depends heavily on those observations. This example is deliberately small, but it shows why a high-leverage cluster can make a sandwich estimate particularly unstable.

Unequal cluster sizes create a related concern. A large cluster may contribute much more to the fitted estimating equation than a small one. If a few clusters dominate, the nominal count \(K\) can overstate how many independent contributions effectively determine the estimate. The following calculation makes this visible even when the individual outcomes are independent.

Worked Example: Unequal Cluster Sizes Under an Independence Working Model

Suppose \(K=4\) clusters contain \(m_1=1,m_2=1,m_3=1,m_4=9\) independent observations each, all distributed as \(N(\mu,1)\). Use the common-mean estimating equation with working independence. There are \(N=12\) observations, \(\hat\mu\) is their overall mean, and the cluster score is \(U_i(\hat\mu)=\sum_{j=1}^{m_i}(Y_{ij}-\hat\mu)\). The bread is \(N\), so CR0 is \(\sum_i U_i(\hat\mu)^2/N^2\).

For each cluster, the fitted residual sum has variance \(m_i-m_i^2/N\). This follows by writing \(U_i(\hat\mu)=\sum_j(Y_{ij}-\mu)-m_i(\hat\mu-\mu)\) and using the covariance induced by the shared overall mean. Consequently,

E[\hat V_{\mathrm{CR0}}] =\frac{1}{N^2}\sum_i\left(m_i-\frac{m_i^2}{N}\right) =\frac{1}{N}\left(1-\sum_i\frac{m_i^2}{N^2}\right).

Here \(\sum_i m_i^2=1+1+1+81=84\), so the expected CR0 variance is \((1/12)(1-84/144)=5/144\). The true variance of \(\hat\mu\) is \(1/12=12/144\). The ratio is \(5/12\), substantially below one. The largest cluster holds \(9/12\) of all observations; treating the four clusters as four equally informative contributions would hide this imbalance.

A useful descriptive quantity is \(K_{\mathrm{eff}}=1/\sum_i w_i^2\), where \(w_i=m_i/N\). In this example, \(K_{\mathrm{eff}}=144/84=12/7\), far below the nominal \(K=4\). This is a diagnostic for imbalance, not a universal replacement for \(K\) in an inferential formula.

When Asymptotic Inference Is Most at Risk

The examples point to several warning signs. They do not provide a sharp threshold at which inference becomes valid: reliability depends on the design and the distribution of cluster contributions as well as on the count.

Few independent clusters. The asymptotic approximation concerns \(K\), not the total number of individual records. When \(K\) is small, the empirical meat is built from few terms, and normal critical values can be unreliable.
High cluster leverage. A cluster that strongly determines one or more fitted coefficients can have unusually constrained residual scores. A small number of influential clusters is especially concerning.
Unequal cluster sizes or influence. If a few large clusters dominate the estimating equation, the effective amount of independent information can be much smaller than the nominal cluster count.
Many fitted parameters relative to clusters. Fitting more mean parameters uses more of the available cluster-level information. When the model is rich relative to \(K\), residual cluster scores have little room to reflect their full variability.
Fragile working-model estimates. Estimating dispersion or correlation parameters from few clusters adds uncertainty. A plug-in sandwich may not account well for that uncertainty in a small sample.

A common mistake is to argue, “There are many observations, so the robust standard error must be reliable.” That reasoning counts dependent observations as though they were independent replications. The correction is to assess the number and influence of the independent clusters, the leverage of the fitted model, and the balance of cluster sizes. A large total sample does not repair a small number of independent units.

Corrections and Practical Interpretation

One familiar adjustment multiplies CR0 by a scalar such as \(K/(K-1)\), or by a factor that also accounts for the number \(p\) of fitted parameters. In the four-cluster intercept example, multiplying by \(4/3\) removes the exact expectation bias. But that result is specific to that simple setting. A scalar correction does not generally remove leverage effects, address severe cluster-size imbalance, or make a normal reference distribution accurate. Likewise, changing the reference distribution alone does not restore information absent from the data.

A defensible analysis should report the number of independent clusters and examine whether a few clusters dominate the fit. When the count is small, describe the sandwich inference as potentially fragile, and consider whether the study design or a suitable small-sample method supports stronger conclusions. Do not present a conventional robust standard error as automatically trustworthy merely because software labels it “robust.”

Key takeaway: Cluster-robust inference is an approximation across independent clusters. Fitting the mean can shrink residual cluster scores, and few clusters, high leverage, or extreme size imbalance can make the ordinary sandwich variance too small. The direction is not guaranteed in every GEE, but the risk must be assessed at the cluster level.

Check Your Understanding

Use the calculations and diagnostics in this tutorial to answer the following.

A study has 900 repeated measurements from 10 clinics. Which count is central to the large-sample justification for the cluster sandwich, and why?
In the Gaussian linear-model theorem, what feature of the fitted residuals causes the expected meat to differ from the target meat?
In the three-cluster regression example, why is the slope especially affected by the high-leverage clusters?
For cluster sizes \(2,2,2,6\), calculate \(K_{\mathrm{eff}}=1/\sum_i w_i^2\), where \(w_i=m_i/\sum_jm_j\). How does it compare with the nominal cluster count?
Why does multiplying CR0 by a degrees-of-freedom factor not guarantee reliable inference for an arbitrary GEE with few clusters?