Robust Standard Errors for GEE Coefficients

Why GEE Standard Errors Use Clusters

In generalized estimating equations, responses within a cluster may be correlated, while distinct clusters are treated as independent. The working covariance describes the within-cluster dependence used to estimate the regression coefficients. If that working covariance is not the true covariance, a standard error calculated as though it were correct can misrepresent coefficient uncertainty. The empirical sandwich covariance instead uses the observed variation in the estimating contributions from independent clusters.

The key unit in this calculation is therefore the cluster, not an individual response considered in isolation. Each cluster contributes a vector to the coefficient estimating equation. The empirical sandwich measures how variable those vectors are across clusters and combines that variation with the sensitivity of the estimating equation to changes in the coefficients.

Cluster Contributions and the Sandwich Formula

Suppose there are \(m\) independent clusters. In cluster \(i\), let \(Y_i\) be the \(n_i\)-vector of responses, let \(\mu_i(\beta)\) be its marginal mean vector, and let \(D_i=\partial\mu_i(\beta)/\partial\beta^T\), an \(n_i\)-by-\(p\) matrix. Write \(V_i\) for the working covariance matrix, which may depend on fitted mean and correlation parameters. The GEE contribution from cluster \(i\) is

U_i(\beta)=D_i^T V_i^{-1}\{Y_i-\mu_i(\beta)\}.

The coefficient estimate \(\hat\beta\) solves \(\sum_{i=1}^m U_i(\hat\beta)=0\). After fitting, define the sensitivity, often called the bread, and the empirical variability, often called the meat, by

\widehat H=\sum_{i=1}^m D_i^T V_i^{-1}D_i, \qquad \widehat B=\sum_{i=1}^m U_i(\hat\beta)U_i(\hat\beta)^T.

Here every quantity on the right is evaluated at the fitted parameters. In particular, the residual vector in \(U_i(\hat\beta)\) is \(Y_i-\mu_i(\hat\beta)\), and the contribution is calculated for the whole cluster. Provided \(\widehat H\) is nonsingular, the empirical sandwich covariance estimator is

\widehat{\operatorname{Cov}}(\hat\beta) = \widehat H^{-1}\widehat B\,\widehat H^{-1}.

The estimated robust standard error for coefficient \(\hat\beta_j\) is the square root of the \(j\)th diagonal entry of this matrix. The off-diagonal entries estimate covariances between coefficient estimates. The formula is called a sandwich because the meat \(\widehat B\) is enclosed by two bread factors \(\widehat H^{-1}\).

Definition: The cluster-level estimating-function contribution is \(U_i(\hat\beta)\), the \(p\)-vector contributed by cluster \(i\) to the fitted GEE. The empirical sandwich covariance uses the sum of the outer products of these vectors as its meat. The contributions are vectors, so each outer product \(U_iU_i^T\) is a \(p\)-by-\(p\) matrix.

Worked Calculations

Worked Example: A Scalar Intercept Estimate

Suppose four independent clusters each contribute one scalar response to an intercept-only GEE with working variance one. The responses are 2, 4, 4, and 10. The estimating contribution is \(U_i(\beta)=Y_i-\beta\), so the fitted intercept is their average:

\hat\beta=\frac{2+4+4+10}{4}=5.

The fitted contributions are \(2-5=-3\), \(4-5=-1\), \(4-5=-1\), and \(10-5=5\). They sum to zero, as required by the estimating equation. The sensitivity is \(\widehat H=4\), and the meat is

\widehat B=(-3)^2+(-1)^2+(-1)^2+5^2=9+1+1+25=36.

Thus

\widehat{\operatorname{Var}}(\hat\beta)=\frac{1}{4}(36)\frac{1}{4} =\frac{36}{16}=2.25, \qquad \widehat{\operatorname{SE}}(\hat\beta)=\sqrt{2.25}=1.5.

For comparison, the model-based variance in this simple working model is \(\widehat H^{-1}=1/4\), with standard error \(0.5\). The difference illustrates the role of the empirical meat: it reflects the observed spread of cluster contributions rather than imposing the working variance as the full source of uncertainty.

Worked Example: Two Coefficients and a Matrix Meat

Consider four independent clusters, each with one response and covariate \(x_i\), in a linear mean model \(\mu_i=\beta_0+\beta_1x_i\) with working variance one. Let the covariates be \((-1,-1,1,1)\), and let the responses be \((0,2,3,7)\). The fitted intercept is the grand mean, 3. The fitted slope is the difference between the means at \(x=1\) and \(x=-1\), divided by 2:

\hat\beta_0=3, \qquad \hat\beta_1=\frac{(3+7)/2-(0+2)/2}{2}=2.

The fitted values are \(1,1,5,5\), so the residuals are \(-1,1,-2,2\). The derivative row for observation \(i\) is \((1,x_i)\). Multiplying each row by its residual gives the four cluster contributions:

U_1=\begin{pmatrix}-1\\1\end{pmatrix},\quad U_2=\begin{pmatrix}1\\-1\end{pmatrix},\quad U_3=\begin{pmatrix}-2\\-2\end{pmatrix},\quad U_4=\begin{pmatrix}2\\2\end{pmatrix}.

The sensitivity and meat are

\widehat H= \begin{pmatrix}4&0\\0&4\end{pmatrix}, \qquad \widehat B= \sum_{i=1}^4 U_iU_i^T = \begin{pmatrix}10&6\\6&10\end{pmatrix}.

For example, the first two outer products sum to \(\begin{pmatrix}2&-2\\-2&2\end{pmatrix}\), and the last two sum to \(\begin{pmatrix}8&8\\8&8\end{pmatrix}\). Since \(\widehat H^{-1}\) is one quarter of the identity matrix,

\widehat{\operatorname{Cov}}(\hat\beta) = \frac{1}{16} \begin{pmatrix}10&6\\6&10\end{pmatrix} = \begin{pmatrix}0.625&0.375\\0.375&0.625\end{pmatrix}.

Both robust standard errors are \(\sqrt{0.625}\approx0.791\). The estimated covariance between the intercept and slope is \(0.375\). Keeping the off-diagonal entries matters when assessing combinations of coefficients, even though each individual standard error uses only a diagonal entry.

Worked Example: How Uneven Cluster Contributions Affect the Meat

Consider an intercept-only GEE with six independent clusters, each containing two responses, and working covariance equal to the two-dimensional identity matrix. Then \(\widehat H=12\). Suppose the fitted cluster contributions are \((-3,-1,0,0,1,3)\); they sum to zero. The meat and robust standard error are

\widehat B=9+1+0+0+1+9=20, \qquad \widehat{\operatorname{SE}}(\hat\beta) =\sqrt{\frac{20}{12^2}} =\frac{\sqrt{20}}{12} \approx0.373.

For a second data set with the same six clusters and the same sensitivity, suppose the contributions are \((-1,-1,0,0,1,1)\). They also sum to zero, but now

\widehat B=1+1+0+0+1+1=4, \qquad \widehat{\operatorname{SE}}(\hat\beta) =\sqrt{\frac{4}{12^2}} =\frac{1}{6} \approx0.167.

The number of clusters and the sensitivity are identical in these two calculations. The larger standard error in the first comes from more variable cluster contributions. In general, a few large contributions can dominate the empirical meat, so it is useful to inspect their sizes as well as the total cluster count.

Why the Sandwich Has This Form

The sandwich formula follows from approximating the estimating equation near the true coefficient. The next result states the conditions needed for that argument. It is an asymptotic result: the number of independent clusters \(m\) increases, while the number of coefficients \(p\) stays fixed.

Theorem: Suppose clusters are independent and identically distributed, the true coefficient is \(\beta_0\), and \(E[U_i(\beta_0)]=0\). Assume \(U_i\) has finite second moments, \(m^{-1/2}\sum_{i=1}^m U_i(\beta_0)\xrightarrow{d}N(0,B)\), and \(m^{-1}\widehat H\xrightarrow{p}H\), where \(H\) and \(B\) are finite and \(H\) is nonsingular. Assume also that \(\hat\beta\xrightarrow{p}\beta_0\), the estimating equation is differentiable near \(\beta_0\), and its summed derivative is locally well approximated by \(-\widehat H\). If a working-correlation parameter is estimated, assume its plug-in effect on the first-order expansion is negligible. Then \(\sqrt m(\hat\beta-\beta_0)\xrightarrow{d}N(0,H^{-1}BH^{-1})\).

Expand the estimating equation.
Write \(0=\sum_i U_i(\hat\beta)\). A first-order Taylor expansion around \(\beta_0\) gives \(0=\sum_i U_i(\beta_0)-\widehat H(\hat\beta-\beta_0)+r_m\), where the remainder \(r_m\) is negligible at the scale of the summed estimating contributions under the stated local approximation.

Isolate the coefficient error.
Rearranging and multiplying by \(\sqrt m\) yields \(\sqrt m(\hat\beta-\beta_0)=(\widehat H/m)^{-1}\{m^{-1/2}\sum_i U_i(\beta_0)\}+o_p(1)\).

Apply the central limit theorem.
The scaled sum converges to a normal vector with covariance \(B\), and \(\widehat H/m\) converges in probability to \(H\). Slutsky's theorem gives the stated covariance \(H^{-1}BH^{-1}\).

For completeness, the covariance transformation in the last step is \(H^{-1}B(H^{-1})^T\). Here \(H\) is a limit of the symmetric sensitivity matrices \(m^{-1}\widehat H\), so \(H^{-1}\) is symmetric and this equals \(H^{-1}BH^{-1}\). This derivation explains why the meat must be built from whole-cluster contributions: the central limit theorem applies to the independent cluster-level vectors \(U_i\), not to response-level pieces treated as independent.

Result: Under the theorem's conditions, if \(m^{-1}\widehat B\xrightarrow{p}B\) and \(m^{-1}\widehat H\xrightarrow{p}H\), then the empirical sandwich consistently estimates the asymptotic covariance of \(\hat\beta\). In particular, its \(j\)th diagonal entry estimates the asymptotic variance of \(\hat\beta_j\).

To see this consistency directly, substitute the scaled matrices into the estimator:

\widehat H^{-1}\widehat B\,\widehat H^{-1} = \frac{1}{m} \left(\frac{\widehat H}{m}\right)^{-1} \left(\frac{\widehat B}{m}\right) \left(\frac{\widehat H}{m}\right)^{-1}.

By the assumed probability limits and continuity of matrix inversion at a nonsingular matrix, the right side is \(m^{-1}\{H^{-1}BH^{-1}+o_p(1)\}\). The quantity in braces estimates the covariance in the limiting distribution of \(\sqrt m(\hat\beta-\beta_0)\); the factor \(1/m\) gives the corresponding scale for \(\hat\beta\). Since taking the square root is continuous for a positive variance, the square root of a diagonal entry estimates its standard error. \(\square\)

What Cluster Count and Variability Tell You

The empirical meat is a sum of \(m\) outer products. Its reliability depends on having enough independent clusters for that sum to represent the population variability of the contributions. A large number of responses inside only a few clusters does not provide the same independent information about between-cluster variability as a large number of clusters. Cluster sizes can affect the contribution vectors, but they do not change which units are independent in this calculation.

The number of clusters alone is not the whole story. If a small number of clusters contribute much larger vectors than the rest, their outer products can dominate \(\widehat B\). Conversely, contributions that are similar across clusters may yield a comparatively stable empirical meat. This is why it is useful to examine the cluster contributions and consider whether a few clusters drive the estimated covariance, while remembering that such inspection does not itself establish that the standard errors are accurate.

The sandwich is robust to misspecification of the working covariance under the regularity conditions behind the estimating-equation argument; it is not a guarantee of accurate finite-sample uncertainty in every data set. The standard formula shown here is the unadjusted empirical sandwich. It should not be silently mixed with an additional finite-sample multiplier: any adjustment changes the estimator and requires its own justification.

Common Mistakes

Forming the meat from individual residuals as if responses were independent. The GEE contribution is \(D_i^TV_i^{-1}\{Y_i-\mu_i\}\) for a whole cluster. Use its outer product, not a collection of response-level terms that discards the cluster structure.
Using the working covariance as the robust meat. The working covariance enters the contribution and the sensitivity. The empirical meat is instead \(\sum_i U_iU_i^T\).
Taking a square root of the entire covariance matrix. A coefficient's standard error is the square root of its own diagonal variance entry. Off-diagonal entries are covariances, not standard errors.
Assuming many observations mean many independent contributions. For this empirical covariance, the independent units are clusters. The number of response rows and the number of clusters need not be equal.

Key takeaway: Calculate one fitted estimating contribution per independent cluster, sum their outer products for the empirical meat, and sandwich that meat between inverse sensitivity matrices. The resulting diagonal square roots are the robust coefficient standard errors, whose reliability depends on the number and variability of independent cluster contributions.

Check Your Understanding

Use the formulas and distinctions in this tutorial to answer the following.

For a two-coefficient GEE, what are the dimensions of \(U_i\), \(U_iU_i^T\), and the empirical meat?
In an intercept-only GEE with working variance one and three cluster responses 1, 2, and 6, calculate the fitted intercept, the fitted contributions, the meat, and the robust standard error.
Why does the empirical sandwich use one contribution per cluster rather than treating every response as an independent contribution?
If two fits have the same sensitivity and cluster count but one has a larger empirical meat, what happens to its robust covariance?
What assumptions about the number of independent clusters and the sensitivity matrix support the asymptotic sandwich argument?