Asymptotic Normality and Valid Sandwich Inference

Why Cluster-Level Asymptotics Matter

A GEE estimate is obtained by adding estimating-function contributions across clusters. Observations within a cluster may be correlated, so the independent units for large-sample reasoning are the clusters, not the individual observations. Consequently, a large number of observations concentrated in only a few clusters does not by itself justify the usual large-sample normal approximation.

The asymptotic covariance has two distinct ingredients. The sensitivity describes how the expected estimating function changes as the parameter changes; the variability describes the actual variation in the cluster contributions. The sandwich combines both. Its validity does not require the working covariance used to construct the GEE to equal the true covariance, but it does require suitable independence, moment, regularity, and identification conditions.

Cluster Contributions and the Asymptotic Setup

Let $n$ be the number of independent clusters, and let $\beta_0\in\mathbb{R}^p$ be the parameter defined by the marginal mean model. Cluster $i$ contributes the $p$-vector $\psi_{ni}(\beta)$ to the estimating function. In a GEE this is typically a weighted residual contribution such as $D_i^T V_i^{-1}(Y_i-\mu_i)$. The notation allows cluster sizes and covariate patterns to vary with $i$ and $n$.

Definition: The estimating function is $U_n(\beta)=\sum_{i=1}^n\psi_{ni}(\beta)$. The normalized sensitivity and variability matrices at the true parameter are $$ A_n=-\frac{1}{n}\sum_{i=1}^n E\!\left[\frac{\partial\psi_{ni}(\beta_0)}{\partial\beta^T}\right], \qquad B_n=\frac{1}{n}\sum_{i=1}^n \operatorname{Var}\{\psi_{ni}(\beta_0)\}. $$ When their limits exist, write $A_n\to A$ and $B_n\to B$. Here $A$ is $p\times p$, and $B$ is $p\times p$.

A common sufficient setup has independent clusters, $E[\psi_{ni}(\beta_0)]=0$, $A_n\to A$, $B_n\to B$, with $A$ nonsingular. A multivariate Lindeberg condition on the centered cluster contributions ensures $$ \frac{1}{\sqrt n}\sum_{i=1}^n\psi_{ni}(\beta_0)\xrightarrow{d}N(0,B). $$ For example, a sufficient moment condition is a uniform bound on $n^{-1}\sum_i E\|\psi_{ni}(\beta_0)\|^{2+\delta}$ for some $\delta>0$, together with the corresponding vanishing-tail condition; the latter follows from this bound. For varying cluster sizes, these conditions must be checked for the cluster contributions themselves. Bounded cluster sizes and uniformly bounded moments are one familiar route, but bounded sizes are not essential if the required limit and tail conditions hold.

We also require differentiability near $\beta_0$, consistency of the root $\hat\beta$ of $U_n(\beta)=0$, and convergence of the normalized derivative uniformly in a neighborhood shrinking to $\beta_0$: along the line segment from $\beta_0$ to $\hat\beta$, $$ -\frac{1}{n}\frac{\partial U_n(\beta)}{\partial\beta^T}\xrightarrow{p}A. $$ These conditions make the local linear approximation accurate at the scale $n^{-1/2}$. Consistency with independent clusters is the prerequisite that places the root in this local neighborhood.

Asymptotic Normality

Theorem: Suppose the clusters are independent, $E[\psi_{ni}(\beta_0)]=0$, $A_n\to A$ with $A$ nonsingular, and $B_n\to B$. Suppose the cluster-score central limit theorem above holds, the estimating function is differentiable with the stated uniform derivative convergence, and a consistent estimator $\hat\beta$ solves $U_n(\hat\beta)=0$. Then $$ \sqrt n(\hat\beta-\beta_0)\xrightarrow{d}N\!\left(0,A^{-1}BA^{-T}\right). $$

Proof. Apply the multivariate mean-value expansion to $U_n(\hat\beta)$ around $\beta_0$. Equivalently, use the integrated derivative along the segment between these two parameter values. Since $U_n(\hat\beta)=0$, $$ 0=U_n(\beta_0)+\left\{\int_0^1 \frac{\partial U_n\bigl(\beta_0+t(\hat\beta-\beta_0)\bigr)}{\partial\beta^T}\,dt\right\}(\hat\beta-\beta_0). $$ By the derivative condition, the negative of the matrix in braces, divided by $n$, converges in probability to $A$. It is therefore invertible with probability tending to one, and its inverse after normalization converges to $A^{-1}$. Rearranging gives $$ \sqrt n(\hat\beta-\beta_0) = \left\{-\frac{1}{n}\int_0^1 \frac{\partial U_n\bigl(\beta_0+t(\hat\beta-\beta_0)\bigr)}{\partial\beta^T}\,dt\right\}^{-1} \frac{U_n(\beta_0)}{\sqrt n}. $$ The first factor converges in probability to $A^{-1}$; the second converges in distribution to $N(0,B)$. Slutsky’s theorem gives the stated limit, whose covariance is $A^{-1}B(A^{-1})^T=A^{-1}BA^{-T}$. $\square$

The order of the matrices matters: the sensitivity acts on both sides of the variability matrix. If the working covariance is correct, the two matrices may simplify under appropriate model assumptions. In general, they need not agree, and the sandwich retains both.

Worked Example: A Scalar Estimator with Unequal Cluster Variability

Suppose the independent cluster contributions are $\psi_{ni}(\beta)=Z_i-\beta$, where $Z_i\sim N(2,1)$ for even $i$ and $Z_i\sim N(2,3)$ for odd $i$. The estimating equation gives $\hat\beta=n^{-1}\sum_i Z_i$. For even $n$, calculate the asymptotic variance.

The true value is $\beta_0=2$, because each cluster has $E[Z_i]=2$. The derivative of each contribution is $-1$, so $A_n=1$ and $A=1$. Half the cluster variances are 1 and half are 3, giving $$ B_n=\frac{1}{n}\left(\frac n2\cdot1+\frac n2\cdot3\right)=2. $$ The normal limit is therefore $$ \sqrt n(\hat\beta-2)\xrightarrow{d}N(0,2). $$ Thus the approximate variance of $\hat\beta$ is $2/n$, not the variance obtained by pretending all clusters had variance 1. Independent clusters with different contributions are allowed; their average variability determines $B$.

Estimating the Limiting Covariance

The limiting matrices are unknown, so inference uses their empirical counterparts. Define the normalized empirical sensitivity and cluster “meat” by $$ \hat A_n=-\frac{1}{n}\sum_{i=1}^n\frac{\partial\psi_{ni}(\hat\beta)}{\partial\beta^T}, \qquad \hat B_n=\frac{1}{n}\sum_{i=1}^n\psi_{ni}(\hat\beta)\psi_{ni}(\hat\beta)^T. $$ The cluster-robust sandwich estimator of the covariance of $\hat\beta$ is $$ \widehat{\operatorname{Var}}(\hat\beta) =\frac{1}{n}\hat A_n^{-1}\hat B_n\hat A_n^{-T}. $$ Each outer product uses one whole cluster contribution. It does not treat the observations within a cluster as separate independent score contributions.

Theorem: In addition to the asymptotic-normality conditions, suppose $\hat A_n\xrightarrow{p}A$ and $\hat B_n\xrightarrow{p}B$, with $A$ nonsingular. Then $$ n\,\widehat{\operatorname{Var}}(\hat\beta)\xrightarrow{p}A^{-1}BA^{-T}. $$

Justification. The assumed convergence of $\hat A_n$ to a nonsingular matrix implies $\hat A_n^{-1}\xrightarrow{p}A^{-1}$. Also, $\hat B_n\xrightarrow{p}B$. Matrix multiplication is continuous, so the continuous mapping theorem gives $$ \hat A_n^{-1}\hat B_n\hat A_n^{-T}\xrightarrow{p}A^{-1}BA^{-T}. $$ The definition of the sandwich estimator then yields the claim. $\square$

These consistency conditions can be checked directly. For example, a law of large numbers for independent cluster outer products gives $n^{-1}\sum_i\psi_{ni}(\beta_0)\psi_{ni}(\beta_0)^T\xrightarrow{p}B$ under suitable moment and tail conditions. To replace $\beta_0$ by $\hat\beta$, one needs the average change in these outer products to converge to zero in probability. Smoothness and appropriate moment bounds provide a common sufficient condition. The analogous derivative law of large numbers, together with uniform control near $\beta_0$, yields $\hat A_n\xrightarrow{p}A$. Independence alone is not enough: the relevant averages must also stabilize, and no small set of clusters may dominate them.

Worked Example: Checking a Two-Parameter Sandwich Calculation

Suppose the limiting sensitivity and variability matrices are $$ A=\begin{pmatrix}1&0\\-1/2&1\end{pmatrix}, \qquad B=\begin{pmatrix}1&0\\0&2\end{pmatrix}. $$ Compute the limiting covariance matrix.

The inverse sensitivity and its transpose are $$ A^{-1}=\begin{pmatrix}1&0\\1/2&1\end{pmatrix}, \qquad A^{-T}=\begin{pmatrix}1&1/2\\0&1\end{pmatrix}. $$ First multiply $A^{-1}B$, then multiply by $A^{-T}$: $$ A^{-1}B= \begin{pmatrix}1&0\\1/2&2\end{pmatrix}, \qquad A^{-1}BA^{-T} = \begin{pmatrix}1&0\\1/2&2\end{pmatrix} \begin{pmatrix}1&1/2\\0&1\end{pmatrix} = \begin{pmatrix}1&1/2\\1/2&9/4\end{pmatrix}. $$ In the lower-right entry, the product is $(1/2)(1/2)+2(1)=1/4+2=9/4=2.25$. The resulting matrix is symmetric, as a covariance matrix must be, and its determinant is $9/4-1/4=2>0$, so it is positive definite. The off-diagonal entry is $1/2$, not $-1/2$.

Worked Example: Binary Outcomes Within Independent Pairs

There are $n$ independent clusters, each containing two binary outcomes $Y_{i1},Y_{i2}$. Each has marginal success probability $p=1/2$, and within-cluster correlation $\rho=0.4$. Use an independence-working GEE for the common logit intercept $\beta$, so $p=\{1+\exp(-\beta)\}^{-1}$, and find the limiting variance of $\sqrt n(\hat\beta-\beta_0)$.

For a pair, the score contribution is $$ \psi_i(\beta)=Y_{i1}+Y_{i2}-2p. $$ At $\beta_0=0$, $p=1/2$ and $dp/d\beta=p(1-p)=1/4$, so the derivative of the score is $-2(1/4)=-1/2$. Hence $A=1/2$. Each outcome has variance $1/4$, and their covariance is $\rho\sqrt{(1/4)(1/4)}=0.4/4=0.1$. Therefore $$ B=\operatorname{Var}(Y_{i1}+Y_{i2}) =\frac14+\frac14+2(0.1)=0.7. $$ The sandwich variance is $$ A^{-1}BA^{-T}=\frac{0.7}{(0.5)^2}=2.8. $$ For a valid joint distribution with these margins and correlation, take $P(Y_{i1}=1,Y_{i2}=1)=P(Y_{i1}=0,Y_{i2}=0)=0.35$ and each discordant-pair probability to be $0.15$. These probabilities sum to 1, give marginal success probability $0.5$, and give covariance $0.35-0.25=0.1$. The calculation thus corresponds to an actual clustered binary model. Treating the two outcomes as independent would instead use $B=0.5$ and produce variance $0.5/(0.5)^2=2$, understating the true asymptotic variance here.

Inference and a Common Pitfall

Once the sandwich is consistent and the normal approximation is adequate, a coefficient can be divided by its sandwich standard error to form a Wald statistic. For a scalar component $\beta_j$, the estimated standard error is the square root of the $j$th diagonal entry of $\widehat{\operatorname{Var}}(\hat\beta)$. This supports large-sample confidence intervals and tests, provided the number of independent clusters is sufficiently informative for the approximation.

Common mistake: counting observations instead of clusters. The incorrect reasoning is: “There are many observations, so the independent-observation central limit theorem applies to their individual residual contributions.” That fails when observations within a cluster are dependent: the sum of contributions in a cluster is the independent unit, and its variance includes within-cluster covariance. The correction is to form one score vector per cluster, apply the limit argument to those independent vectors, and use their outer products in the empirical meat. The asymptotic rate is governed by $n$, the number of independent clusters, under the conditions stated above.

Key takeaway: GEE asymptotic normality combines a cluster-score central limit theorem with local linearization of the estimating equation. Valid sandwich inference additionally requires empirical sensitivity and cluster-score variability to converge to their population limits.

Check Your Understanding

Answer each question by referring to the cluster-level quantities in the tutorial.

What roles do the sensitivity matrix $A$ and variability matrix $B$ play in the limiting covariance?
Why must the outer product in the empirical sandwich be formed from an entire cluster contribution?
In the binary-pair example, what would change in the limiting variance if the within-pair covariance were zero?
State two conditions, beyond independent clusters, that support the cluster-score central limit theorem.
Why does consistency of $\hat A_n$ to a nonsingular matrix imply consistency of its inverse?