Why Sandwich Inference Fails With Few Clusters

Why the Number of Clusters Matters

GEE inference is built from contributions that are independent across clusters. The empirical sandwich estimates variability by combining those contributions, but it must also estimate the regression parameters from the same finite collection of clusters. When that collection is small, fitting can make the observed cluster contributions look less variable than the underlying contributions really are. The resulting covariance estimate may be too small, so Wald intervals can be too narrow and tests can reject too often.

The relevant sample size for this approximation is primarily the number of independent clusters, not the total number of individual observations. Adding many responses within a few clusters may improve information about each cluster, but does not provide the same number of independent contributions as observing many clusters.

Definition: Let $U_j(\beta)$ be the estimating-equation contribution from independent cluster $j$, for $j=1,\ldots,m$, and let $\hat\beta$ solve $\sum_{j=1}^m U_j(\beta)=0$. Define the empirical bread and meat by $$ \hat A=-\sum_{j=1}^m \frac{\partial U_j(\hat\beta)}{\partial\beta^\top}, \qquad \hat B=\sum_{j=1}^m U_j(\hat\beta)U_j(\hat\beta)^\top. $$ The empirical sandwich covariance is $\hat V=\hat A^{-1}\hat B\hat A^{-\top}$, when $\hat A$ is invertible. It estimates the covariance of $\hat\beta$; its square roots on the diagonal are the sandwich standard errors.

What the Large-Cluster Argument Requires

Under standard regularity conditions, the sandwich is justified by a cluster-level central limit theorem and a law of large numbers. For a simple version, suppose the clusters are independent and identically distributed, $E[U_j(\beta_0)]=0$, the score contributions have finite second moments, and $$ A=E\left[-\frac{\partial U_j(\beta_0)}{\partial\beta^\top}\right] \quad\text{and}\quad B=E[U_j(\beta_0)U_j(\beta_0)^\top] $$ exist, with $A$ nonsingular. Also suppose the estimator is consistent, its estimating equation admits a valid first-order Taylor expansion around $\beta_0$, and the empirical bread and meat converge to their respective expectations after division by $m$. Then the cluster-level central limit theorem and the Taylor expansion give $$ \sqrt m(\hat\beta-\beta_0)\xrightarrow{d} N\left(0,A^{-1}BA^{-\top}\right). $$ The empirical sandwich estimates the corresponding covariance at scale $1/m$.

These are large-$m$ statements: they do not promise accurate normal or chi-square approximations for a small fixed $m$. In finite samples, the fitted scores $U_j(\hat\beta)$ are not the same as the unobservable contributions $U_j(\beta_0)$. In particular, fitting forces their sum to be zero. This constraint can reduce the apparent variation in the empirical meat. The size and direction of the resulting error depend on the design and model; downward bias is common, not a universal identity for every GEE problem.

A Finite-Sample Calculation

The intercept-only mean model gives an exact illustration. Regard each cluster as contributing one independent observation $Y_j$, where $Y_j\sim N(\mu,\sigma^2)$, with $\sigma^2>0$. The estimating equation is $U_j(\mu)=Y_j-\mu$. Its solution is $\hat\mu=\bar Y$, and the empirical bread is $m$. Thus the uncorrected empirical sandwich variance is the sum of squared fitted scores divided by $m^2$.

Result: In this model, $$ \hat V_{\mathrm{sand}}=\frac{\sum_{j=1}^m(Y_j-\bar Y)^2}{m^2}, \qquad E[\hat V_{\mathrm{sand}}]=\frac{m-1}{m}\operatorname{Var}(\bar Y). $$ Under $H_0:\mu=0$, the Wald statistic formed from this uncorrected variance has distribution $\sqrt{m/(m-1)}\,t_{m-1}$, not $t_{m-1}$.

To verify the expectation, write $S=\sum_{j=1}^m(Y_j-\bar Y)^2$. The usual sample variance is $S/(m-1)$, which is unbiased for $\sigma^2$, so $E[S]=(m-1)\sigma^2$. Therefore $E[\hat V_{\mathrm{sand}}]=(m-1)\sigma^2/m^2$, whereas $\operatorname{Var}(\bar Y)=\sigma^2/m$. Their ratio is $(m-1)/m$.

For the distribution statement, under the null, $Z=\sqrt m\,\bar Y/\sigma$ is standard normal and $Q=S/\sigma^2$ has a chi-square distribution with $m-1$ degrees of freedom; $Z$ and $Q$ are independent. The raw sandwich Wald statistic is $$ W=\frac{\bar Y}{\sqrt{S/m^2}} =\frac{m\bar Y}{\sqrt S} =\sqrt{\frac{m}{m-1}}\, \frac{Z}{\sqrt{Q/(m-1)}}. $$ The final ratio is $t_{m-1}$ by the definition of the Student $t$ distribution. This proves the stated result. $\square$

The factor matters: the finite-sample variance expectation is smaller by $(m-1)/m$, and the raw standard error is smaller by a factor of $\sqrt{(m-1)/m}$ relative to a variance estimate with the corresponding degrees-of-freedom adjustment. This calculation is not a claim that every GEE sandwich has this exact bias factor. It is a transparent example of why “robust” does not mean “automatically well calibrated with few clusters.”

Worked Example: Five Cluster-Level Outcomes

Suppose five independent clusters contribute outcomes $2,3,4,5,6$, and the target is their common mean. The fitted mean is $$ \bar Y=(2+3+4+5+6)/5=4. $$

The fitted cluster scores are $-2,-1,0,1,2$, so their squared sum is $4+1+0+1+4=10$. The bread is $5$, and the raw empirical sandwich variance and standard error are $$ \hat V_{\mathrm{sand}}=10/5^2=0.4, \qquad \widehat{\operatorname{SE}}_{\mathrm{sand}}=\sqrt{0.4}\approx0.632. $$

The sample variance of the five outcomes is $10/(5-1)=2.5$. Its square root divided by $\sqrt5$ is $\sqrt{2.5/5}=\sqrt{0.5}\approx0.707$, which is larger than the raw sandwich standard error. This arithmetic illustrates the finite-sample difference between the empirical meat as defined above and a variance using a degrees-of-freedom adjustment. It does not establish the unknown true variance from these five observations.

Consequences for Wald Tests

A scalar Wald test divides an estimate minus its null value by its estimated standard error and compares the result with a standard normal reference distribution. For a $q$-dimensional contrast $C\beta=r$, the usual statistic is $$ W=(C\hat\beta-r)^\top(C\hat V C^\top)^{-1}(C\hat\beta-r), $$ which is compared with a chi-square distribution with $q$ degrees of freedom under the large-cluster approximation. If $\hat V$ is systematically too small, the statistic tends to be too large. Even if the covariance estimate is not biased downward, a small number of independent terms can make its sampling distribution far from the approximation used for calibration.

Worked Example: Null Rejection With Six Clusters

Return to the normal intercept model with $m=6$ and $H_0:\mu=0$. The uncorrected empirical sandwich gives $$ W=\frac{\bar Y}{\sqrt{\sum_{j=1}^{6}(Y_j-\bar Y)^2/36}} =\sqrt{\frac65}\,t_5 $$ in distribution under the null. A two-sided test using the standard-normal cutoff $1.96$ rejects when $|W|>1.96$. Equivalently, it rejects when $$ |t_5|>1.96\sqrt{\frac56}\approx1.789. $$

Consequently the exact null rejection probability is $2P(t_5>1.96\sqrt{5/6})$, approximately $0.13$, rather than $0.05$. The statistic is more dispersed than a standard normal because it is the $t_5$ ratio multiplied by $\sqrt{6/5}$. In particular, this calculation uses the raw sandwich variance $S/36$ throughout; using $S/[6(5)]$ would define a different variance estimate and a different statistic. The discrepancy is consequential even in this unusually simple model.

When Few-Cluster Problems Are More Severe

The count of clusters is a first diagnostic, not a complete one. A nominally moderate cluster count can still provide weak evidence for a particular coefficient if only a few clusters meaningfully inform it. The following features deserve special attention:

High leverage. A few clusters may have unusual covariate patterns or large influence on the fitted coefficient. The fitted scores from those clusters can dominate the empirical meat.
Unequal cluster sizes. Large clusters can contribute much more information and influence than small ones. The effective number of clusters for a coefficient may be well below the raw count.
Few clusters in an important group. A treatment comparison is especially fragile when treatment is assigned at cluster level and only a small number of clusters receive one treatment.
Strong within-cluster dependence. Many responses within a cluster do not become independent replications. Treating the total number of responses as the basis for a large-sample approximation overstates the available independent information.
Many fitted parameters relative to clusters. Estimating several coefficients from a small number of independent score contributions leaves less information for estimating their variability.

These features do not imply that GEE itself is unsuitable. They mean that the usual large-cluster covariance and reference distribution need a credibility assessment rather than automatic acceptance. The working correlation and mean model remain important, but choosing a plausible working correlation cannot manufacture additional independent clusters.

Worked Example: A Cluster-Level Simulation Check

Suppose a proposed analysis has six independent clusters, and the scientific null is a zero mean. To assess the calibration of the raw sandwich Wald test in the simple normal setting above, simulate $R=100{,}000$ data sets. For each replicate $r$, generate six independent $Y_{rj}\sim N(0,1)$, calculate the mean and the raw sandwich standard error, and record whether $|W_r|>1.96$.

The estimate of the rejection probability is the proportion of replicates that reject. Its Monte Carlo standard error is approximately $\sqrt{\hat p(1-\hat p)/R}$. Near $\hat p=0.13$, this is $\sqrt{0.13(0.87)/100{,}000}\approx0.0011$. Thus a simulation should typically report a rejection rate near $0.13$, with Monte Carlo variation of roughly one tenth of a percentage point. The exact benchmark derived above makes this check independently verifiable.

For an actual GEE design, simulate from a scientifically plausible model that preserves the planned cluster sizes, covariates, treatment allocation, and within-cluster dependence. Refit the GEE and calculate the intended sandwich Wald test in every replicate. Under a null scenario, the empirical rejection rate estimates the test's type I error for that design; under alternatives, the rejection rate estimates power. Report the number of clusters and the Monte Carlo uncertainty, and examine whether a small subset of clusters drives the result. A simulation only assesses the scenarios encoded in it, so it complements rather than replaces design judgment.

Common Mistake

Wrong reasoning: “There are hundreds of patient records, so the sandwich standard error must be reliable, even though they come from only a few clinics.” This counts observations as if they were independent contributions to the sandwich approximation.

Correction: The asymptotic argument aggregates independent cluster contributions. If patients within a clinic are correlated, the clinic—not each patient—is the independent unit for this argument. Report the number of clusters, inspect cluster sizes and influence, and use design-specific simulation when the count is small or the allocation is unbalanced. The raw empirical sandwich is not a guarantee of valid small-sample inference.

Key takeaway: Sandwich inference is robust to some forms of covariance misspecification under large-cluster conditions; it is not automatically reliable with few clusters. Check which clusters identify the contrast, distinguish the raw empirical sandwich from any adjusted version, and evaluate the actual test's null calibration when the asymptotic approximation is doubtful.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

Why does the number of individual observations within a cluster not replace the number of independent clusters in the sandwich asymptotic argument?
For $m$ independent normal cluster-level outcomes, derive the expectation of the raw empirical sandwich variance for their mean.
Under the normal null model with six clusters, what is the distribution of the raw sandwich Wald statistic, and why is it not $t_5$?
Name two features of a GEE design that can make its effective information for a coefficient much smaller than the raw cluster count suggests.
What does a null simulation's empirical rejection rate assess, and why should its Monte Carlo uncertainty be reported?