Mean Misspecification and Robust GEE Inference

What the GEE Sandwich Can—and Cannot—Protect

A generalized estimating equation (GEE) uses a working covariance to combine responses within a cluster. If the marginal mean model is correct, a wrong working correlation does not, by itself, make the mean estimating equation biased. The sandwich covariance is designed to account for the resulting discrepancy between the estimating equation’s sensitivity and the actual variability of its cluster-level scores.

A different problem arises when the marginal mean model is wrong. The estimating equation may still have a well-defined solution, and the sandwich may estimate that solution’s sampling variability. But the solution is then generally a pseudo-true coefficient, not a parameter that makes the specified mean model correct. Robust standard errors do not remove that distinction or repair the model’s interpretation.

Estimating Equations and Their Target

Definition: For independent clusters $i=1,\ldots,n$, let $Y_i$ be the response vector and $X_i$ the covariates. A GEE has cluster score $g_i(\beta)=D_i(\beta)^T W_i(\beta)^{-1}\{Y_i-\mu_i(\beta)\}$, where $\mu_i(\beta)$ is the specified marginal mean vector, $D_i(\beta)=\partial\mu_i(\beta)/\partial\beta^T$, and $W_i(\beta)$ is a positive-definite working covariance matrix. The estimator $\hat\beta$ solves $\sum_{i=1}^n g_i(\beta)=0$. The pseudo-true coefficient $\beta^*$ is a root of $E[g_i(\beta)]=0$. If the mean model is correct at $\beta_0$, then typically $E[g_i(\beta_0)\mid X_i]=0$, so $\beta_0$ is a root; if the mean model is wrong, $\beta^*$ instead describes the root selected by the chosen estimating equation.

The definition assumes independent, identically distributed clusters. Within-cluster responses need not be independent. For clarity, take cluster sizes to be bounded, the score to have finite second moments near $\beta^*$, and the expected score to have a unique root there. Let $A=-E[\partial g_i(\beta^*)/\partial\beta^T]$ and $B=\operatorname{Var}\{g_i(\beta^*)\}$. Inference requires $A$ to be nonsingular; a nondegenerate asymptotic distribution additionally requires positive variance in the contrasts being considered.

The “meat” $B$ is the covariance of the whole cluster score, not a sum of marginal variances that assumes the working covariance is true. Under mean misspecification, score means conditional on covariates can vary across clusters. Their variation contributes to $B$, even though the unconditional score has mean zero at $\beta^*$.

Why the Sandwich Estimates Variability Around the Root

Theorem: Suppose the clusters are iid; $E[g_i(\beta^*)]=0$; $A$ is nonsingular; $B$ is finite; and $\hat\beta\xrightarrow{p}\beta^*$. Assume the score is differentiable near $\beta^*$, its sample derivative converges in probability to its expectation, and the iid central limit theorem applies to $g_i(\beta^*)$. Then $$ \sqrt{n}(\hat\beta-\beta^*)\xrightarrow{d}N\!\left(0,A^{-1}BA^{-T}\right). $$

Proof. Write $\bar g_n(\beta)=n^{-1}\sum_i g_i(\beta)$. Since $\bar g_n(\hat\beta)=0$, a Taylor expansion about $\beta^*$ gives $$ 0=\bar g_n(\beta^*)+\left\{\frac{1}{n}\sum_{i=1}^n \frac{\partial g_i(\tilde\beta)}{\partial\beta^T}\right\}(\hat\beta-\beta^*), $$ where $\tilde\beta$ lies between $\hat\beta$ and $\beta^*$ (coordinatewise, or in the integral form of the expansion). Consistency and the assumed derivative convergence make the matrix in braces converge in probability to $-A$. Thus $$ \sqrt{n}(\hat\beta-\beta^*)= \left\{-\frac{1}{n}\sum_{i=1}^n \frac{\partial g_i(\tilde\beta)}{\partial\beta^T}\right\}^{-1} \sqrt{n}\bar g_n(\beta^*). $$ The iid central limit theorem gives $\sqrt n\,\bar g_n(\beta^*)\xrightarrow{d}N(0,B)$. The matrix factor converges in probability to $A^{-1}$; Slutsky’s theorem yields the stated limit. $\square$

A practical sandwich estimator uses $$ \widehat{\operatorname{Var}}(\hat\beta)=\frac{1}{n}\hat A^{-1}\hat B\hat A^{-T}, \qquad \hat A=-\frac{1}{n}\sum_{i=1}^n \frac{\partial g_i(\hat\beta)}{\partial\beta^T}, \qquad \hat B=\frac{1}{n}\sum_{i=1}^n g_i(\hat\beta)g_i(\hat\beta)^T. $$ Under the theorem’s conditions and the corresponding laws of large numbers, $\hat A\xrightarrow{p}A$ and $\hat B\xrightarrow{p}B$, so this estimates the covariance in the limit for $\hat\beta$. If the working covariance contains estimated nuisance parameters, the same conclusion requires conditions ensuring that their estimation does not change this first-order limit; the theorem as stated treats the score as specified.

Worked Example: Wrong Working Correlation, Correct Mean

Suppose each cluster contains two responses with true means $E[Y_{i1}]=E[Y_{i2}]=\mu$, variances $\sigma^2$, and correlation $\rho=0.6$. Clusters are independent. Fit an intercept-only GEE using the working covariance $W_i=\sigma^2 I$, which assumes zero within-cluster correlation.

With $\beta=\mu$, the score is $g_i(\mu)=\{(Y_{i1}-\mu)+(Y_{i2}-\mu)\}/\sigma^2$. Its expectation is zero at the true mean regardless of $\rho$. The sensitivity and score variance are $$ A=-E\left[\frac{\partial g_i}{\partial\mu}\right] =\frac{2}{\sigma^2}, \qquad B=\operatorname{Var}\{g_i(\mu)\} =\frac{2\sigma^2+2\rho\sigma^2}{\sigma^4} =\frac{2(1+\rho)}{\sigma^2}. $$ Therefore the asymptotic variance of $\hat\mu$ is $$ \frac{1}{n}A^{-1}BA^{-T} =\frac{\sigma^2(1+\rho)}{2n} =\frac{0.8\sigma^2}{n}. $$ The working-independence model-based variance is $\sigma^2/(2n)=0.5\sigma^2/n$, which understates uncertainty here. The sandwich uses the observed cluster-score variation and, under the stated independent-cluster conditions, estimates the correct variance even though the working correlation is wrong.

Mean Misspecification Changes the Target

When the mean model is wrong, the theorem can still apply: it gives the variability of $\hat\beta$ around $\beta^*$. But $\beta^*$ depends on the estimating equation, including its working weights. It may summarize a weighted compromise among covariate-specific means rather than a true coefficient describing those means. Consequently, “robust” refers to the covariance calculation under the stated sampling assumptions; it does not mean robust to every choice of mean model.

Worked Example: An Omitted Group Effect

Let clusters consist of one observation each. A group indicator $G$ is sampled independently with probabilities $P(G=0)=P(G=1)=1/2$, and $Y\mid G=0\sim\operatorname{Bernoulli}(0.1)$, $Y\mid G=1\sim\operatorname{Bernoulli}(0.9)$. Fit a logistic intercept-only model $p=\operatorname{expit}(\beta)$, omitting $G$.

The score is $g(\beta)=Y-p$. Its expected value is zero when $p=E[Y]=(0.1+0.9)/2=0.5$, so $\beta^*=\operatorname{logit}(0.5)=0$. The fitted model predicts probability $0.5$ in both groups, despite their true probabilities being $0.1$ and $0.9$.

At the root, $A=p(1-p)=0.25$. Because group membership is random in this iid mixture, the marginal distribution of $Y$ is Bernoulli with probability $0.5$, and hence $B=\operatorname{Var}(Y)=0.25$. The sandwich variance is therefore $$ \frac{1}{n}\frac{B}{A^2} =\frac{1}{n}\frac{0.25}{0.25^2} =\frac{4}{n}. $$ This calculation includes the variation in conditional means across groups. It does not make the intercept a group-specific effect: it estimates the pooled marginal mean on the logit scale under this sampling mixture.

If instead the design fixes exactly equal numbers in the two groups, the conditional sampling variance of the pooled response mean is $0.09/n$, because the average within-group Bernoulli variance is $\{0.1(0.9)+0.9(0.1)\}/2=0.09$. The corresponding logit variance is $0.09/(0.25^2n)=1.44/n$. That is a different, fixed-allocation design calculation. The iid pooled sandwich above includes between-group score variation and is not an estimator of this within-stratum variance; the design must be reflected in the variance calculation.

Worked Example: A Misspecified Logistic Slope

Let independent observations have $X\in\{-1,0,1\}$, each value with probability $1/3$. Conditional on $X=-1,0,1$, let $Y$ be Bernoulli with success probabilities $0.2,0.8,0.2$, respectively. Fit the logistic mean model $\operatorname{expit}(\beta_0+\beta_1X)$ with independent working variance.

At $\beta_1=0$, the fitted probability is constant. The intercept equation requires that probability to equal $E[Y]=(0.2+0.8+0.2)/3=0.4$, so $\beta_0=\operatorname{logit}(0.4)$. The slope equation is $E[X(Y-0.4)]=0$: the contributions at $X=-1$ and $X=1$ cancel, and the $X=0$ contribution is zero. Thus $\beta^*=(\operatorname{logit}(0.4),0)^T$. The zero fitted slope conceals the nonconstant, nonmonotone true probabilities.

For the logistic score $g=(1,X)^T(Y-p)$, the sensitivity matrix at this root is diagonal. Since $p(1-p)=0.4(0.6)=0.24$ and $E[X^2]=2/3$, $$ A_{11}=0.24,\qquad A_{22}=0.24\left(\frac{2}{3}\right)=0.16,\qquad A_{12}=A_{21}=0. $$ The off-diagonal terms vanish because $E[X]=0$. For the slope score, the meat is $$ B_{22}=E[X^2(Y-0.4)^2] =\frac{1}{3}\left\{(0.2-0.4)^2+0+(0.2-0.4)^2\right\} +\frac{1}{3}\{0.16+0+0.16\} =\frac{2}{75}+\frac{8}{75} =\frac{2}{15}. $$ The first term is variation of the conditional means around the fitted mean; the second is average conditional Bernoulli variance. Hence the slope’s asymptotic variance is $$ (A^{-1}BA^{-T})_{22} =\frac{2/15}{(0.16)^2} =\frac{125}{24}. $$ Omitting the conditional-mean term would give an incorrect meat and standard error. Even the correctly calculated sandwich only quantifies uncertainty around the pseudo-true zero slope; it does not make that slope a useful description of the true pattern.

Wald Tests and Practical Interpretation

Theorem: Under the preceding asymptotic conditions, let $C$ be a fixed $r\times p$ matrix of rank $r$. If $C A^{-1}BA^{-T}C^T$ is positive definite, then the sandwich Wald statistic for $H_0:C\beta^*=c$, $$ W=\{C\hat\beta-c\}^T \{C\widehat{\operatorname{Var}}(\hat\beta)C^T\}^{-1} \{C\hat\beta-c\}, $$ converges in distribution to $\chi^2_r$ when $H_0$ holds.

Indeed, the preceding limit theorem and multiplication by $C$ give a centered normal limit for $C\hat\beta-c$, with covariance $C A^{-1}BA^{-T}C^T/n$. Consistency of the sandwich estimate and the continuous mapping theorem give the stated chi-square limit. When the mean model is misspecified, the null concerns $\beta^*$, not automatically a scientifically desired coefficient in the true conditional or marginal response pattern. A small p-value can reject a restriction on the estimating-equation target; a large p-value does not establish that the mean model is correct.

For interpretation, specify what population, covariate distribution, and working equation define the reported coefficient. Then assess whether the specified mean shape is scientifically defensible and whether conclusions change under plausible alternatives. These checks address questions the sandwich does not answer.

Common Mistakes

“Robust” means the whole model is safe. The sandwich protects the variance calculation against working-covariance misspecification under its assumptions; it does not correct a wrong mean structure.
Using only conditional variances for the meat. For iid clusters, $B=\operatorname{Var}\{g_i(\beta^*)\}$. Variation of $E[g_i(\beta^*)\mid X_i]$ across covariates contributes to this variance. The logistic example shows the omitted term explicitly.
Interpreting a pseudo-true coefficient as a true effect. A root of the expected estimating equation can be a weighted compromise. Its interpretation depends on the model and estimating weights.
Ignoring the sampling design. An iid-mixture sandwich need not estimate variance under fixed group counts or other non-iid designs. The variance calculation must match how clusters were sampled.

Key takeaway: With independent clusters and regularity conditions, the sandwich estimates the sampling covariance around the root of the GEE. Correct marginal means make that root a scientifically specified mean parameter; misspecified means generally leave a pseudo-true target whose interpretation and tests require care.

Check Your Understanding

Use the assumptions and examples above to answer these questions.

In the theorem, what are $A$ and $B$, and why is $B$ defined from the complete cluster score?
Under what condition does a wrong working correlation leave the mean estimating equation unbiased at the true coefficient?
In the omitted-group example, why is $B=0.25$ for iid random group membership rather than the average conditional variance $0.09$?
What does a sandwich Wald test assess when the marginal mean model is misspecified?
Why can a zero fitted logistic slope in the final worked example fail to describe the true response pattern?