Deriving the GEE Sandwich Covariance

From Cluster Estimating Functions to Covariance

A GEE combines information across clusters, not across individual observations treated as independent. Its uncertainty therefore depends on two distinct features: how strongly the estimating function changes when the parameter changes, and how much the cluster-level estimating functions vary across independent clusters. A first-order Taylor expansion separates these features and leads to the sandwich covariance.

The derivation is a linearization. It gives an exact covariance for the resulting linearized random term under stated moment conditions. It does not, by itself, establish that the covariance of the original estimator is close to that matrix; that stronger claim requires control of the Taylor remainder in second moment.

Definition: For $m$ independent clusters, let $U_i(\beta)$ be the $p$-dimensional estimating function from cluster $i$, and define the summed and normalized estimating functions by $$ U_m(\beta)=\sum_{i=1}^m U_i(\beta), \qquad \Psi_m(\beta)=\frac{1}{m}U_m(\beta). $$ A GEE estimator $\hat\beta$ is a solution of $U_m(\hat\beta)=0$, or equivalently $\Psi_m(\hat\beta)=0$. At the target value $\beta_0$, define the cluster-score covariance $B=\operatorname{Var}(U_i(\beta_0))$ and the sensitivity matrix $A=-E[\partial U_i(\beta_0)/\partial\beta^T]$, when the clusters are identically distributed and these quantities exist. The transpose is denoted by $T$.

For the usual GEE, a cluster contribution has the form $U_i(\beta)=D_i^T V_i^{-1}\{Y_i-\mu_i(\beta)\}$. Here $Y_i$ is the response vector, $\mu_i(\beta)$ its marginal mean vector, $D_i=\partial\mu_i(\beta)/\partial\beta^T$, and $V_i$ is the working covariance matrix. The sensitivity is the negative derivative of the whole estimating function; it is not defined as the inverse working covariance. The cluster-score covariance $B$ is the variance of the entire vector $U_i(\beta_0)$, including any within-cluster dependence in the responses.

The First-Order Linearization

The expansion is made around $\beta_0$, where the mean estimating function is zero. To state the argument precisely, suppose $E[U_i(\beta_0)]=0$, the cluster contributions are independent and identically distributed with finite second moments, and the normalized estimating function is continuously differentiable on a neighborhood of $\beta_0$. Suppose also that $\hat\beta\xrightarrow{p}\beta_0$, that $\sqrt m\,\Psi_m(\beta_0)=O_p(1)$, and that the integrated negative derivative along the segment between $\beta_0$ and $\hat\beta$ converges in probability to a nonsingular matrix $A$. These are explicit conditions for the linearization below; in applications, derivative convergence and consistency must be justified for the particular GEE.

Theorem: Under these conditions, $$ \sqrt m(\hat\beta-\beta_0) = A^{-1}\frac{1}{\sqrt m}\sum_{i=1}^m U_i(\beta_0)+o_p(1). $$ This is a first-order representation of the estimator, not an assertion about its exact finite-sample covariance.

To prove the representation, define the integrated sensitivity $$ \hat A_m = -\int_0^1 \frac{\partial\Psi_m\{\beta_0+t(\hat\beta-\beta_0)\}}{\partial\beta^T}\,dt. $$ By the fundamental theorem of calculus applied along this segment, $$ \Psi_m(\hat\beta) = \Psi_m(\beta_0)-\hat A_m(\hat\beta-\beta_0). $$

Since $\Psi_m(\hat\beta)=0$, this identity gives $\hat\beta-\beta_0=\hat A_m^{-1}\Psi_m(\beta_0)$ whenever $\hat A_m$ is invertible. Convergence $\hat A_m\xrightarrow{p}A$ and nonsingularity of $A$ imply that $\hat A_m$ is invertible with probability tending to one and $\hat A_m^{-1}\xrightarrow{p}A^{-1}$. Multiplying by $\sqrt m$ yields $$ \sqrt m(\hat\beta-\beta_0) = A^{-1}\sqrt m\,\Psi_m(\beta_0) + (\hat A_m^{-1}-A^{-1})\sqrt m\,\Psi_m(\beta_0). $$ The first term is $A^{-1}m^{-1/2}\sum_i U_i(\beta_0)$. The second is $o_p(1)O_p(1)=o_p(1)$, proving the representation. $\square$

Covariance of the Linearized Term

The random term on the right side of the representation is a deterministic linear transformation of the summed cluster scores. Its covariance follows directly from independence. This calculation is the source of the two sides of the sandwich: sensitivity on the outside and cluster-score variability in the middle.

Theorem: Let $U_1(\beta_0),\ldots,U_m(\beta_0)$ be independent, identically distributed, mean-zero $p$-vectors with finite covariance matrix $B$, and let $A$ be a fixed nonsingular matrix. Define the deterministic-sensitivity linearized term $$ \tilde\beta = \beta_0+A^{-1}\frac{1}{m}\sum_{i=1}^m U_i(\beta_0). $$ Then its covariance is exactly $$ \operatorname{Var}(\tilde\beta) = \frac{1}{m}A^{-1}BA^{-T}, $$ where $A^{-T}=(A^{-1})^T$.

For the proof, the mean of $\tilde\beta$ is $\beta_0$, because each cluster score has mean zero. Independence makes the cross-covariances between distinct cluster scores zero, so $$ \operatorname{Var}\left(\frac{1}{m}\sum_{i=1}^m U_i(\beta_0)\right) = \frac{1}{m^2}\sum_{i=1}^m \operatorname{Var}(U_i(\beta_0)) = \frac{B}{m}. $$ For any fixed matrix $C$, $\operatorname{Var}(CZ)=C\operatorname{Var}(Z)C^T$. Taking $C=A^{-1}$ gives the stated formula. $\square$

Thus $A^{-1}BA^{-T}/m$ is the exact covariance of this deterministic-sensitivity linearized term. The Taylor representation alone does not make it the exact covariance, or even necessarily a covariance approximation, for $\hat\beta$. A remainder that is small in probability can still have a large second moment on rare events.

Worked Example: A Scalar Cluster Estimating Equation

Suppose each independent cluster contributes $U_i(\beta)=Y_i-\beta$, where $E[Y_i]=\beta_0=2$ and $\operatorname{Var}(Y_i)=3$. The estimating equation is $\sum_i(Y_i-\beta)=0$.

The derivative of each contribution with respect to $\beta$ is $-1$, so $A=-E[-1]=1$. The cluster-score variance at the target is $B=\operatorname{Var}(Y_i-2)=3$. The linearized term is $$ \tilde\beta = 2+\frac{1}{m}\sum_{i=1}^m(Y_i-2) = \frac{1}{m}\sum_{i=1}^mY_i. $$ Its covariance, here a scalar variance, is exactly $A^{-1}BA^{-T}/m=3/m$. In this example the estimating equation is linear, so the root itself equals the displayed linearized term; its variance is exactly $3/m$, not merely a first-order approximation.

Worked Example: Repeated Binary Outcomes

Consider independent clusters, each containing two binary responses $Y_{i1},Y_{i2}$, with common marginal mean $p_0=0.5$. Suppose each response has variance $0.25$, and their within-cluster correlation is $0.4$. Use an intercept estimating equation with working covariance $V=0.25R$, where $$ R= \begin{pmatrix} 1&0.2\\ 0.2&1 \end{pmatrix}. $$ The working correlation is $0.2$, which differs from the actual correlation $0.4$.

For the scalar parameter $p$, $D=(1,1)^T$, and the cluster contribution is $U_i(p)=D^TV^{-1}(Y_i-pD)$. Since $R^{-1}D=(1/1.2)D$, this becomes $$ U_i(p)=\frac{4}{1.2}\{Y_{i1}+Y_{i2}-2p\}. $$ Differentiating gives $A=-E[\partial U_i(p_0)/\partial p]=8/1.2=6.6667$.

The actual covariance of the two responses is $0.4(0.25)=0.1$. Therefore $$ \operatorname{Var}(Y_{i1}+Y_{i2}) = 0.25+0.25+2(0.1) = 0.7. $$ It follows that $B=\operatorname{Var}(U_i(p_0))=(4/1.2)^2(0.7)=7.7778$. The sandwich variance of the linearized estimator is consequently $$ \frac{1}{m}A^{-2}B = \frac{1}{m}\left(\frac{1.2}{8}\right)^2(7.7778) = \frac{0.175}{m}. $$ As a check, the mean of the two responses in one cluster has variance $0.7/4=0.175$, matching the sandwich calculation. The working correlation affects the sensitivity and score variance separately, but those factors combine to give this result.

Worked Example: A Two-Parameter Sandwich

Suppose a cluster estimating function is locally linear: $$ U_i(\beta)=Z_i-A(\beta-\beta_0), \qquad A= \begin{pmatrix} 2&1\\ 0&1 \end{pmatrix}. $$ Let $Z_i$ take the four values $(1,0)^T,(-1,0)^T,(0,2)^T,(0,-2)^T$, each with probability $1/4$. Then $E[Z_i]=0$, and direct calculation gives $$ B=\operatorname{Var}(Z_i) = \begin{pmatrix} 1/2&0\\ 0&2 \end{pmatrix}. $$

The sensitivity is $A$, because $-\partial U_i/\partial\beta^T=A$. Its inverse and the middle product are $$ A^{-1} = \begin{pmatrix} 1/2&-1/2\\ 0&1 \end{pmatrix}, \qquad A^{-1}B = \begin{pmatrix} 1/4&-1\\ 0&2 \end{pmatrix}. $$ Multiplying by $A^{-T}$ gives $$ A^{-1}BA^{-T} = \begin{pmatrix} 5/8&-1\\ -1&2 \end{pmatrix}. $$ Thus the covariance of the linearized term is $m^{-1}\begin{pmatrix}5/8&-1\\-1&2\end{pmatrix}$. The matrix is symmetric, and its determinant is $(5/8)(2)-1=1/4>0$, consistent with a positive-definite covariance matrix.

When the Formula Does—and Does Not—Approximate Estimator Covariance

The first-order representation has a remainder $R_m$, defined by $$ \hat\beta=\tilde\beta+R_m, \qquad \tilde\beta=\beta_0+A^{-1}\frac{1}{m}\sum_{i=1}^mU_i(\beta_0). $$ The representation proved above yields $R_m=o_p(m^{-1/2})$ under its conditions. That statement concerns probability, not second moments. It does not rule out a remainder that is very large on events whose probabilities shrink rapidly.

A sufficient additional condition for covariance approximation is $E[\|R_m\|^2]=o(1/m)$, along with the finite covariance condition for the linearized term. To see why, write centered variables using $X^\circ=X-E[X]$. Then $\operatorname{Var}(R_m)$ has entries bounded in magnitude by $E[\|R_m\|^2]=o(1/m)$. Each cross-covariance between $\tilde\beta$ and $R_m$ is bounded by Cauchy–Schwarz by a constant times $$ \sqrt{E[\|\tilde\beta-E[\tilde\beta]\|^2]E[\|R_m\|^2]} = O(m^{-1/2})o(m^{-1/2}) = o(1/m). $$ Expanding the covariance of $\tilde\beta+R_m$ therefore shows that its difference from $\operatorname{Var}(\tilde\beta)$ is $o(1/m)$. Without such second-moment control, that covariance conclusion does not follow from the Taylor expansion in probability.

Common Pitfall

Wrong reasoning: “The Taylor remainder is $o_p(m^{-1/2})$, so the estimator’s covariance must be $A^{-1}BA^{-T}/m$.” This treats convergence in probability as if it guaranteed convergence of second moments.
Correction: The sandwich is exactly the covariance of the deterministic-sensitivity linearized term under the stated independence and finite-variance conditions. Approximating the estimator’s covariance additionally requires suitable second-moment control of the remainder, such as $E[\|R_m\|^2]=o(1/m)$.

The distinction matters in practice: the sandwich derivation identifies the two sources of variability, but an estimator-level covariance claim needs both a valid linearization and adequate control of its errors. The cluster, rather than each response within it, is the independent unit in the middle covariance calculation.

Check Your Understanding

Use the definitions and derivations above to check the following points.

For an estimating function $U_i(\beta)$, write the definitions of $A$ and $B$, and state what each measures.
Why does independence across clusters make the covariance of $m^{-1}\sum_i U_i(\beta_0)$ equal to $B/m$?
In the repeated-binary example, what actual response covariance was used to calculate $B$, and why was it not the working covariance?
What is the exact scope of the formula $A^{-1}BA^{-T}/m$ under the conditions given?
Why does $R_m=o_p(m^{-1/2})$ alone not establish that the estimator covariance is approximated by the linearized covariance?