Bread, Meat, and Model-Based Covariance

What the Sandwich Pieces Measure

A GEE estimate is built from cluster-level estimating functions. Its covariance estimate has two different inputs: how sensitively the estimating equation changes as the regression parameters change, and how much the observed cluster contributions vary. Separating these inputs clarifies what the sandwich estimator is estimating—and why it need not equal the covariance obtained by treating the working covariance as correct.

Let clusters be indexed by $i=1,\ldots,m$. For cluster $i$, let $Y_i$ be its response vector, $\mu_i(\beta)$ its marginal mean vector, $D_i=\partial\mu_i(\beta)/\partial\beta^T$, and $V_i$ its working covariance matrix. The cluster estimating function is $U_i(\beta)=D_i^T V_i^{-1}\{Y_i-\mu_i(\beta)\}$. The GEE estimate $\hat\beta$ solves $\sum_i U_i(\hat\beta)=0$.

Definition: In the usual GEE covariance calculation, the sensitivity matrix (the bread) and the empirical cluster-score covariance (the meat) are $$ A(\beta)=\sum_{i=1}^m D_i^T V_i^{-1}D_i, \qquad B(\beta)=\sum_{i=1}^m U_i(\beta)U_i(\beta)^T. $$ When $A(\hat\beta)$ is nonsingular, the sandwich covariance estimate is $A(\hat\beta)^{-1}B(\hat\beta)A(\hat\beta)^{-T}$. In the standard GEE form $A$ is symmetric, so this is $A(\hat\beta)^{-1}B(\hat\beta)A(\hat\beta)^{-1}$. The model-based covariance estimate is $A(\hat\beta)^{-1}$.

The bread describes sensitivity: a large $A$, in a given parameter direction, means that the estimating equation changes substantially for a small change in that parameter. Its inverse therefore translates a given amount of score variation into less parameter variation. The meat records the outer products of whole-cluster estimating functions. It measures variation in the contributions that the independent clusters make to the estimating equation. Each $U_i$ is a vector, so its outer product retains variation in each parameter direction and covariation between directions.

The units help explain the construction. The meat has the scale of score covariance; multiplying on the left and right by inverse sensitivities expresses that variability on the parameter scale. It is a “sandwich” because the meat sits between two slices of bread.

Why the Sandwich Propagates Score Variation

The covariance calculation can be understood by linearizing the estimating equation near the target parameter $\beta_0$. Suppose the mean model is correctly specified, the clusters are independent, and the estimating functions have mean zero at $\beta_0$. Suppose also that the estimating equation is differentiable and its relevant sensitivity is nonsingular. The first-order approximation is $\hat\beta-\beta_0\approx A(\beta_0)^{-1}\sum_i U_i(\beta_0)$. This approximation explains the covariance formula; it does not assert that the sandwich estimate is an exact finite-sample covariance for $\hat\beta$.

Result: Under that first-order linearization, if the cluster estimating functions are independent with finite covariance matrices, then $$ \operatorname{Cov}(\hat\beta-\beta_0) \ \approx\ A(\beta_0)^{-1} \left\{\sum_{i=1}^m \operatorname{Cov}\bigl(U_i(\beta_0)\bigr)\right\} A(\beta_0)^{-T}. $$ The empirical meat $B(\hat\beta)$ estimates the middle, score-variability component.

To see the covariance propagation explicitly, write $S=\sum_i U_i(\beta_0)$. For fixed matrices $C$, the covariance identity $\operatorname{Cov}(CS)=C\operatorname{Cov}(S)C^T$ follows directly from centering: $CS-E[CS]=C\{S-E[S]\}$. Applying this with $C=A(\beta_0)^{-1}$ gives the two outer factors. Independence gives $\operatorname{Cov}(S)=\sum_i\operatorname{Cov}(U_i)$, since each cross-covariance between distinct clusters is zero. These two identities yield the displayed expression. Replacing the unknown score covariances with their empirical cluster outer products gives the sandwich estimate. $\square$

For a sample covariance of cluster scores, the outer products are not centered by subtracting the sample average score. At the target parameter the estimating functions have mean zero, and in the fitted sample their total is zero by the estimating equation. The conventional GEE meat is the sum of the cluster outer products evaluated at $\hat\beta$. It is not a covariance matrix computed from individual observations as if they were independent.

When the Model-Based Covariance Appears

The model-based expression uses the same bread on both sides, but omits an empirical meat. This omission has a precise justification if the working covariance is the actual conditional covariance of each response vector. Treat covariates and design quantities as fixed in the following calculation. At the true parameter, assume $E[Y_i]=\mu_i$, $\operatorname{Cov}(Y_i)=V_i$, and that $D_i$ and $V_i$ are nonrandom given those quantities.

Result: Under these assumptions, $$ \operatorname{Cov}\bigl(U_i(\beta_0)\bigr)=D_i^T V_i^{-1}D_i. $$ If clusters are independent, the covariance of the total estimating function is therefore $A(\beta_0)$. Consequently, the first-order covariance in the preceding result reduces to $A(\beta_0)^{-1}$.

Proof. Since $E[Y_i-\mu_i]=0$, the score contribution has mean zero. For a fixed matrix $C$, $\operatorname{Cov}(CZ)=C\operatorname{Cov}(Z)C^T$; this follows by writing $CZ-E[CZ]=C\{Z-E[Z]\}$ in the definition of covariance. Here $C=D_i^TV_i^{-1}$, and $V_i$ is a covariance matrix, so it is symmetric and $V_i^{-T}=V_i^{-1}$. Thus

\begin{aligned} \operatorname{Cov}(U_i) &=D_i^T V_i^{-1}\operatorname{Cov}(Y_i)V_i^{-1}D_i\\ &=D_i^T V_i^{-1}V_iV_i^{-1}D_i\\ &=D_i^T V_i^{-1}D_i. \end{aligned}

Summing over independent clusters gives $\operatorname{Cov}(\sum_i U_i)=\sum_iD_i^TV_i^{-1}D_i=A$. Substitution into the first-order covariance expression gives $A^{-1}AA^{-T}=A^{-1}$, because $A$ is symmetric. This is the model-based covariance. $\square$

The argument identifies the key comparison. The model-based covariance substitutes the sensitivity $A$ for the score covariance in the middle of the sandwich. The robust estimate instead uses the empirical $B$. If the working covariance is correct, their population targets agree in this calculation; if the observed cluster scores vary differently from that implication, the sandwich estimate reflects that difference.

Worked Comparisons

Worked Example: A Scalar Robust Variance Larger Than the Model-Based Variance

Suppose there is one regression parameter and, at the fitted GEE, the sensitivity is $A=20$. Suppose four cluster scores are $5,-5,0,0$. Their sum is zero, as required at a solution of the estimating equation, and the empirical meat is $B=5^2+(-5)^2+0^2+0^2=50$.

The model-based variance estimate is $\hat{\operatorname{Var}}_{\mathrm{MB}}(\hat\beta)=1/A=1/20=0.05$. The robust variance estimate is $\hat{\operatorname{Var}}_{\mathrm{rob}}(\hat\beta)=B/A^2=50/20^2=50/400=0.125$. The corresponding standard errors are $\sqrt{0.05}\approx0.224$ and $\sqrt{0.125}\approx0.354$.

Here the observed score variability is larger than the model-based covariance would imply: the meat is 50, whereas the bread-based implication for the total score covariance is 20. The sandwich variance is consequently larger. This comparison is about the score variation and sensitivity at the fitted model, not a universal rule that robust standard errors must be larger.

Worked Example: Comparing Two-Parameter Covariance Matrices

Consider the fitted matrices $$ A=\begin{pmatrix}10&0\\0&5\end{pmatrix}, \qquad B=\begin{pmatrix}40&10\\10&25\end{pmatrix}. $$ The meat is a valid positive-definite matrix: its leading diagonal entry is 40 and its determinant is $40(25)-10^2=900>0$. It can, for example, be formed from four cluster score vectors $v,-v,w,-w$, where $v=(2\sqrt{5},\sqrt{5}/2)^T$ and $w=(0,3\sqrt{5}/2)^T$. Their outer products sum to $B$.

The inverse bread is $A^{-1}=\begin{pmatrix}1/10&0\\0&1/5\end{pmatrix}$. Therefore the model-based covariance is $$ \hat{\operatorname{Var}}_{\mathrm{MB}}(\hat\beta) =\begin{pmatrix}0.1&0\\0&0.2\end{pmatrix}. $$ For the sandwich, premultiplying and postmultiplying $B$ by $A^{-1}$ scales its first row and column by $1/10$, and its second row and column by $1/5$: $$ \hat{\operatorname{Var}}_{\mathrm{rob}}(\hat\beta) =\begin{pmatrix} 40/100&10/50\\ 10/50&25/25 \end{pmatrix} =\begin{pmatrix}0.4&0.2\\0.2&1\end{pmatrix}. $$

The robust variances are larger in both directions in this example, and the robust estimate also has a nonzero covariance between the two coefficient estimates. Its estimated correlation is $0.2/\sqrt{0.4(1)}\approx0.316$; the model-based matrix gives zero correlation. Thus the meat can alter both marginal uncertainty and the estimated relationship between parameter errors.

Worked Example: Cluster Scores Include Within-Cluster Cross-Products

Suppose two clusters each contain two observations, and their individual contributions to a scalar estimating function are $1,2$ in the first cluster and $-1,-2$ in the second. The cluster scores are the sums $U_1=1+2=3$ and $U_2=-1-2=-3$; their total is zero. The correct cluster-level meat is $B=3^2+(-3)^2=18$.

Expanding a cluster's squared score shows why its observations cannot simply be treated as independent contributions: $(1+2)^2=1^2+2^2+2(1)(2)=9$. The cross-product $2(1)(2)=4$ is included in the cluster score outer product. The second cluster similarly contributes $(-1-2)^2=1+4+4=9$. Squaring the four individual contributions separately would instead give $1+4+1+4=10$, omitting both within-cluster cross-products.

If the fitted sensitivity is $A=6$, the model-based variance is $1/6$, while the cluster-robust variance is $18/6^2=18/36=1/2$. The calculation shows why the independent units for the meat are clusters: it preserves the joint variation of responses within each cluster.

Interpreting the Difference and Avoiding a Common Mistake

The sandwich covariance is a plug-in estimate of uncertainty based on the fitted sensitivity and empirical cluster-score variation. The model-based covariance is the corresponding first-order calculation under the working covariance assumption. They are two covariance calculations for the same GEE estimate, not two different coefficient estimates. If the working covariance is correct under the assumptions of the derivation, the two calculations agree at the level of their covariance target. In a finite sample their numerical estimates need not match.

A common mistake is to assume that robust covariance must be larger because it is called “robust,” or that a smaller robust estimate proves the working covariance is correct. Neither conclusion follows. In one parameter, the comparison is $B/A^2$ against $1/A$, so the robust variance is larger exactly when $B>A$, equal when $B=A$, and smaller when $B<A$. With multiple parameters, the comparison is directional: for a vector $c$, compare $c^TA^{-1}BA^{-T}c$ with $c^TA^{-1}c$. There is no general requirement that one covariance matrix exceed the other in every direction.

Do not replace cluster scores by observation-level residuals. The outer product of a cluster sum includes within-cluster cross-products.
Do not interpret the bread alone as empirical variability. It describes sensitivity; the meat supplies the score variability.
Do not treat a numerical difference as a definitive diagnosis. The estimates can differ because of finite-sample score variation as well as a mismatch between the working covariance and actual response variation.

Key takeaway: The bread translates score changes into parameter changes; the meat records how cluster scores vary. The model-based covariance uses the bread as the score covariance implied by a correct working covariance, while the sandwich replaces that implication with the empirical cluster-score outer products.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

In the GEE sandwich estimator, what does the sensitivity matrix measure, and what does the empirical meat measure?
If $A=12$ and the scalar empirical meat is $B=9$, calculate both covariance estimates and state which is larger.
Why does the model-based covariance reduce to $A^{-1}$ when the working covariance equals the true response covariance?
For a cluster with two scalar contributions $a$ and $b$, write the cluster contribution to the meat and identify the within-cluster cross-product.
Can a robust covariance estimate be smaller than its model-based counterpart? Explain using the one-parameter comparison.