Why GEE Errors Are Easy to Make
GEE calculations combine a mean model, a working covariance, and independent clusters. Errors often arise when one of these pieces is given the wrong role: an observation is mistaken for an independent cluster, matrix factors are rearranged as if they were numbers, or a working correlation is interpreted as a model for outcomes conditional on one another. Interpretation can also drift from the marginal mean to a subject-specific effect.
The most reliable checks are structural. Track the unit of independence, verify each matrix dimension, and state what expectation the mean model describes. These checks do not require a particular working-correlation structure.
Set Up the Estimating Equation
Let \(Y_i\) be the \(m_i\)-vector of responses in cluster \(i\), for \(i=1,\ldots,n\). Clusters are independent; responses within a cluster may be dependent. Conditional on the covariates, the specified mean vector is \(\mu_i(\beta)=E[Y_i\mid X_i]\). Define \(D_i=\partial\mu_i(\beta)/\partial\beta^T\), a matrix with \(m_i\) rows and \(p\) columns, and let \(V_i\) be the working covariance matrix, with \(m_i\) rows and columns.
The bread is \(p\times p\), as is the meat, so the covariance estimate is \(p\times p\). In this standard formulation \(H_n\) is symmetric, and \(H_n^{-T}=H_n^{-1}\). Keeping the transpose in the general expression helps prevent an order error and makes the covariance structure explicit.
Cluster Scores, Not Observation Scores, Form the Meat
The sandwich meat estimates the variability of the entire estimating equation. Since independent units are clusters, its empirical contributions are the cluster vectors \(U_i\), not separate observation-level pieces treated as independent. To see the issue, write a cluster score as a sum of observation contributions, \(U_i=\sum_j q_{ij}\). Its outer product is
The terms with \(j\ne k\) are within-cluster cross-products. Replacing the cluster outer product with \(\sum_jq_{ij}q_{ij}^T\) drops those terms. That replacement is justified only under additional within-cluster independence conditions; it is not the robust GEE calculation when responses within a cluster may be dependent.
Worked Example: The Missing Within-Cluster Cross-Products
Suppose one scalar-parameter cluster has two observation-level score pieces \(q_{i1}=1\) and \(q_{i2}=1\). The cluster score is \(U_i=1+1=2\), so its contribution to the robust meat is \(U_i^2=4\).
If the pieces are incorrectly treated as independent units, their squared contributions sum to \(1^2+1^2=2\). The missing amount is \(4-2=2\), exactly the two cross-products \(q_{i1}q_{i2}+q_{i2}q_{i1}=1+1\). The cluster calculation retains the covariance information that the observation-level calculation discards.
This arithmetic illustrates the unit-of-independence rule; it does not assert that one cluster suffices to estimate a population variance. In practice, the meat sums \(U_iU_i^T\) over the independent clusters.
Why the Sandwich Order Cannot Be Reversed
For a matrix-valued sandwich, factors generally cannot be rearranged. The robust covariance has one inverse bread factor on each side of the meat. The transpose on the right is important in the general form: it guarantees a symmetric result when the meat is symmetric.
Worked Example: A Reversed Factor Produces a Nonsymmetric Result
Take \(H_n=\begin{pmatrix}1&0\\0&2\end{pmatrix}\) and \(J_n=\begin{pmatrix}1&1\\1&1\end{pmatrix}\). Both are symmetric, and \(H_n^{-1}=\begin{pmatrix}1&0\\0&1/2\end{pmatrix}\).
The correctly ordered sandwich is
If the left inverse is incorrectly applied twice, the result is
The off-diagonal entries differ, so the second matrix is not symmetric and cannot be a covariance matrix. This example shows why “the same factors are present” is not enough: their positions matter.
Deriving the Cluster Sandwich
The robust form follows by linearizing the estimating equation around the true parameter \(\beta_0\). The derivation makes the role of independent clusters explicit. Assume the mean model is correct, clusters are independent, the estimating equation has a consistent root \(\hat\beta\), and the usual differentiability and moment conditions for a Taylor expansion and a central limit theorem hold. Also assume \(H_n/n\) converges to a nonsingular matrix \(H\), \(n^{-1}\sum_i\operatorname{Var}(U_i(\beta_0)\mid X_i)\) converges to \(J\), and no single cluster dominates the sum. These conditions support the expansion and limit below.
Use a Taylor expansion of \(\sum_iU_i(\hat\beta)=0\) about \(\beta_0\). The derivative is approximately \(-H_n\), so \(0=\sum_iU_i(\beta_0)-H_n(\hat\beta-\beta_0)+o_p(\sqrt n)\).
Rearranging gives \(\hat\beta-\beta_0=H_n^{-1}\sum_iU_i(\beta_0)+o_p(n^{-1/2})\), because \(H_n\) is of order \(n\) and the score sum is of order \(\sqrt n\).
Conditional on the covariates, independence across clusters gives \(\operatorname{Var}(\sum_iU_i)=\sum_i\operatorname{Var}(U_i)\). The correct-mean condition gives \(E[U_i\mid X_i]=0\), so no nonzero conditional means add extra terms.
For a matrix \(A\), \(\operatorname{Var}(AZ)=A\operatorname{Var}(Z)A^T\). Setting \(A=H_n^{-1}\) gives the two bread factors on opposite sides of the score-sum variance.
Finally, replacing each unknown cluster score by its fitted value gives \(J_n=\sum_iU_i(\hat\beta)U_i(\hat\beta)^T\). This substitution estimates the cluster-level variability without requiring the working covariance to equal the true covariance. The result uses independent clusters; it does not justify treating responses within a cluster as independent.
A Working Correlation Is a Weighting Device
The working covariance \(V_i\) specifies how residual components are weighted in the estimating equation. It need not be the true covariance. In particular, specifying a working correlation does not by itself specify a conditional distribution of one response given another, nor does it turn the GEE mean into such a conditional mean.
Worked Example: Residual Weighting Is Not a Conditional Prediction
For a two-response Gaussian cluster, suppose the working covariance is \(V_i=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\) and the residual vector at a parameter value is \(r_i=Y_i-\mu_i=\begin{pmatrix}2\\-1\end{pmatrix}\). The inverse working covariance is
Thus \(V_i^{-1}r_i=\begin{pmatrix}10/3\\-8/3\end{pmatrix}\). The score contribution is \(D_i^T\begin{pmatrix}10/3\\-8/3\end{pmatrix}\). This is a transformed residual vector used in estimating \(\beta\). The calculation does not say that the second response has conditional mean \(\mu_{i2}\) given the first response, or provide a conditional distribution for it.
The working covariance affects the estimating weights. The mean specification remains the stated marginal mean, and the validity of robust inference relies on the appropriate cluster-level conditions rather than on interpreting those weights as an outcome model.
Marginal Effects Are Not Subject-Specific Effects
In GEE, \(\mu_{ij}=E[Y_{ij}\mid X_i]\) describes a response mean conditional on covariates, not on the other responses in the cluster. A coefficient is interpreted through this specified mean and link. “Marginal” distinguishes the mean model from a model that conditions on a latent cluster effect; it does not mean that covariates are ignored or averaged over an unspecified covariate distribution.
Worked Example: Interpreting a Logistic GEE Coefficient
Suppose a binary response has the marginal mean model \(\operatorname{logit}(\mu_{ij})=\beta_0+\beta_1x_{ij}\), with \(\beta_0=\operatorname{logit}(0.4)\) and \(\beta_1=0.7\). For \(x=0\), the fitted marginal probability is \(0.4\). For \(x=1\), the log odds increase by \(0.7\), so the marginal odds ratio is \(\exp(0.7)\approx2.014\).
At \(x=1\), the fitted probability is
The coefficient describes the covariate association in the specified marginal mean. It is not automatically a log odds ratio conditional on a latent subject-specific effect, and it is not an effect conditional on the cluster’s other observed responses. Those are different estimands.
Common Mistakes and Checks
| Tempting but incorrect move | Correction |
|---|---|
| Sum outer products of observation scores as if every response were independent | First sum score pieces within each cluster; then form one outer product per independent cluster. |
| Put both inverse bread factors on the same side of the meat | Use \(H_n^{-1}J_nH_n^{-T}\); verify that the result is symmetric. |
| Read a working correlation as a fitted conditional outcome distribution | Treat \(V_i\) as a working residual-weighting matrix unless a separate conditional model has been specified. |
| Describe a marginal GEE coefficient as subject-specific or conditional on other responses | Interpret it through the covariate-conditional marginal mean model and its link. |
A useful final audit is to check that \(U_i\) has \(p\) entries, \(U_iU_i^T\) is \(p\times p\), the meat sums over independent clusters, and the sandwich has one inverse bread on each side. Then state the mean being modeled before translating a coefficient into words.
Check Your Understanding
Use the cluster unit, matrix order, and mean-model interpretation to answer each question.
- If a cluster score is \(q_{i1}+q_{i2}\), which terms are omitted by replacing its outer product with \(q_{i1}q_{i1}^T+q_{i2}q_{i2}^T\)?
- Why should a sandwich covariance estimate be symmetric, and what matrix-order check can detect a reversed-factor error?
- In a GEE, what role does the working covariance play in the estimating equation?
- For a logistic GEE coefficient, explain what is marginal and what the coefficient is not automatically describing.
- Which units must be independent for the cluster-level sandwich derivation to apply?