Sandwich Covariance for Quasi-Likelihood

Why a Sandwich Covariance Is Needed

A quasi-score uses a specified mean and working variance to construct an estimating equation. If the mean model is appropriate but the working variance differs from the observations’ true variance, the equation can still estimate the mean parameters consistently, while the covariance calculated as though the working variance were true can be wrong. The sandwich covariance addresses this distinction: its middle factor estimates the actual variability of the observation-level scores, rather than imposing the working variance on that variability.

The construction is based on three ingredients. The bread measures how quickly the estimating equation changes with the parameter. The meat measures the variation in its individual contributions. The sandwich combines them in the same way that linearizing an estimating equation combines its derivative and its random fluctuation.

Observation-Level Quasi-Scores and the Sandwich

Suppose $Y_1,\ldots,Y_n$ are independent observations, with a model for their means $\mu_i(\beta)$, where $\beta$ is a $p$-dimensional parameter. Write $D_i(\beta)=\partial\mu_i(\beta)/\partial\beta$ as a $p$-dimensional column vector. Let the working variance of $Y_i$ be $\phi v_i(\beta)$, where $v_i(\beta)>0$ and $\phi>0$ is a working dispersion parameter. The observation-level quasi-score contribution and its sum are

\psi_i(\beta)=\frac{D_i(\beta)}{\phi v_i(\beta)}\{Y_i-\mu_i(\beta)\}, \qquad U_n(\beta)=\sum_{i=1}^n\psi_i(\beta).

A quasi-score estimator $\hat\beta$ solves $U_n(\hat\beta)=0$. The contribution $\psi_i$ is a column vector; its outer product $\psi_i\psi_i^T$ is a $p$ by $p$ matrix. This observation-level form matters: the empirical meat is built from these individual vectors, not from the total score alone.

Definition: The empirical sandwich covariance estimator is $$ \widehat{\operatorname{Cov}}_{\mathrm{sand}}(\hat\beta) = \hat H^{-1}\hat M\hat H^{-T}, \qquad \hat H=-\left.\frac{\partial U_n(\beta)}{\partial\beta^T}\right|_{\beta=\hat\beta}, \qquad \hat M=\sum_{i=1}^n\psi_i(\hat\beta)\psi_i(\hat\beta)^T. $$ Here $\hat H$ is the observed bread, $\hat M$ is the empirical meat, and $\hat H^{-T}=(\hat H^{-1})^T$. The estimator requires $\hat H$ to be nonsingular.

The equivalent normalized form uses $\hat A=\hat H/n$ and $\hat B=\hat M/n$:

\widehat{\operatorname{Cov}}_{\mathrm{sand}}(\hat\beta) =\frac{1}{n}\hat A^{-1}\hat B\hat A^{-T}.

These are the same estimator, not two different conventions: the unnormalized bread and meat both grow with $n$, whereas the covariance of $\hat\beta$ shrinks at rate $1/n$.

Worked Example: A Constant Mean With Poisson Working Variance

Suppose four independent counts are $1,2,3,10$, and the mean is a constant $\mu>0$. Use the Poisson working variance $v_i(\mu)=\mu$, with $\phi=1$. Then $\psi_i(\mu)=(Y_i-\mu)/\mu$, and solving $U_4(\mu)=0$ gives

0=\sum_{i=1}^4\frac{Y_i-\hat\mu}{\hat\mu} =\frac{16-4\hat\mu}{\hat\mu}, \qquad \hat\mu=4.

The four fitted score contributions are $-3/4,-1/2,-1/4,3/2$. The bread is obtained by differentiating the total score:

U_4(\mu)=\frac{16}{\mu}-4, \qquad \hat H=-U_4'(4)=\frac{16}{4^2}=1.

The meat is the sum of the four squared contributions:

\hat M=\frac{9}{16}+\frac{4}{16}+\frac{1}{16}+\frac{36}{16} =\frac{50}{16}=\frac{25}{8}.

Thus the empirical sandwich variance is $25/8=3.125$. The variance based directly on the Poisson working model is $\hat\mu/4=1$. The observed counts vary much more than that working variance would suggest, so these two covariance estimates differ substantially. This small-sample calculation illustrates the construction; the sandwich’s consistency is a large-sample claim, not a guarantee of accurate finite-sample calibration.

Large-Sample Basis and Conditions

To state the target precisely, let $\beta_0$ be the parameter at which the estimating equation is centered, so that $E[\psi_i(\beta_0)]=0$ for each $i$. Define the sensitivity and score variability per observation by

A_n=-\frac{1}{n}\sum_{i=1}^n E\left[\frac{\partial\psi_i(\beta_0)}{\partial\beta^T}\right], \qquad B_n=\frac{1}{n}\sum_{i=1}^n E\left[\psi_i(\beta_0)\psi_i(\beta_0)^T\right].

For example, the second expression is the covariance contribution because $E[\psi_i(\beta_0)]=0$. A sufficient set of conditions for the standard sandwich limit is: the observations are independent; $\hat\beta\xrightarrow{p}\beta_0$; $A_n\to A$ and $B_n\to B$, with $A$ nonsingular; the score vectors satisfy a central limit theorem with covariance limit $B$; and the derivatives obey a uniform law of large numbers in a neighborhood of $\beta_0$. For empirical covariance consistency, also require a law of large numbers for the outer products and enough moment and smoothness control that replacing $\beta_0$ by $\hat\beta$ in the scores and derivatives changes the normalized sums by a term converging to zero in probability. These conditions can be established, for instance, using suitable integrable bounds on derivatives and score moments. Merely having a working variance function does not guarantee them.

Theorem: Under these conditions, if $$ \frac{1}{\sqrt n}U_n(\beta_0)\xrightarrow{d}N(0,B), $$ then $$ \sqrt n(\hat\beta-\beta_0)\xrightarrow{d} N\left(0,A^{-1}BA^{-T}\right), \qquad \hat A^{-1}\hat B\hat A^{-T}\xrightarrow{p}A^{-1}BA^{-T}. $$ Consequently, $\hat H^{-1}\hat M\hat H^{-T}$ estimates the asymptotic covariance of $\hat\beta$.

The first conclusion follows by an exact integral expansion, not by assuming that one Jacobian evaluated at a single intermediate point represents the change in a vector-valued function. Put $G_n(\beta)=\partial U_n(\beta)/\partial\beta^T$. Applying the fundamental theorem of calculus componentwise along the line segment from $\beta_0$ to $\hat\beta$ gives

0=U_n(\hat\beta) =U_n(\beta_0)+ \left\{\int_0^1G_n\bigl(\beta_0+t(\hat\beta-\beta_0)\bigr)\,dt\right\} (\hat\beta-\beta_0).

Define the integrated negative derivative

\bar A_n=-\frac{1}{n}\int_0^1G_n\bigl(\beta_0+t(\hat\beta-\beta_0)\bigr)\,dt.

Consistency of $\hat\beta$ and the stated uniform derivative convergence imply $\bar A_n\xrightarrow{p}A$. Since $A$ is nonsingular, $\bar A_n^{-1}\xrightarrow{p}A^{-1}$. Rearranging the exact expansion now yields

\sqrt n(\hat\beta-\beta_0) =\bar A_n^{-1}\frac{U_n(\beta_0)}{\sqrt n} \xrightarrow{d}A^{-1}Z, \qquad Z\sim N(0,B),

by Slutsky’s theorem. The covariance of $A^{-1}Z$ is $A^{-1}BA^{-T}$, proving the distributional claim. For the covariance estimator, the same derivative law gives $\hat A=\hat H/n\xrightarrow{p}A$. The outer-product law and the stipulated control for replacing $\beta_0$ by $\hat\beta$ give $\hat B=\hat M/n\xrightarrow{p}B$. Continuity of matrix inversion at nonsingular $A$ then gives $\hat A^{-1}\hat B\hat A^{-T}\xrightarrow{p}A^{-1}BA^{-T}$, as claimed. $\square$

Working-Variance Covariance and Empirical Meat

For the quasi-score above, suppose the mean model is correct, the observations are independent, and the derivatives and expectations exist. Since $E[Y_i-\mu_i(\beta_0)]=0$, differentiating $\psi_i$ and taking expectations removes terms proportional to the mean residual. The expected negative derivative contribution is therefore

A_n=\frac{1}{n}\sum_{i=1}^n \frac{D_i(\beta_0)D_i(\beta_0)^T}{\phi v_i(\beta_0)}.

By contrast, the meat depends on the true variances $\tau_i^2=\operatorname{Var}(Y_i)$:

B_n=\frac{1}{n}\sum_{i=1}^n \frac{\tau_i^2 D_i(\beta_0)D_i(\beta_0)^T} {\phi^2v_i(\beta_0)^2}.

If the working variance is correct, so that $\tau_i^2=\phi v_i(\beta_0)$, then $B_n=A_n$, and the sandwich reduces to $A_n^{-1}/n$. If it is not correct, that equality need not hold. The empirical meat estimates the score variability without requiring the working variance to equal the true variance. This is the covariance correction; it does not, by itself, repair an incorrectly specified mean model or dependence between observations.

Worked Example: Grouped Binomial Responses With Extra Variation

Consider two independent groups, each containing $m=10$ trials, with observed success counts $2$ and $8$. Fit a common success probability $p$ using the logistic intercept $\beta$, where $p=e^\beta/(1+e^\beta)$, and use the binomial working variance. The group score contributions are $\psi_i=y_i-mp$. Their sum is zero at $\hat p=(2+8)/20=1/2$, so $\hat\beta=0$.

The derivative of each score with respect to $\beta$ is $-mp(1-p)$. Thus the total bread at the estimate is

\hat H=2(10)(1/2)(1/2)=5.

The two fitted score contributions are $2-5=-3$ and $8-5=3$, so $\hat M=(-3)^2+3^2=18$. The empirical sandwich variance for $\hat\beta$ is $18/5^2=0.72$. The working binomial covariance is $\hat H^{-1}=1/5=0.20$. The two calculations differ because the empirical group-level score contributions are more variable than the working binomial calculation represents. With only two groups, neither estimate should be treated as a reliable large-sample approximation; the arithmetic shows precisely which quantities the two methods use.

Worked Example: A Two-Parameter Linear Mean Model

Take three independent observations at $x_i=0,1,2$, with responses $1,2,5$. Fit $\mu_i=\beta_0+\beta_1x_i$ with working variance $1$. Let $z_i=(1,x_i)^T$. The estimating equation is $\sum_i z_i(Y_i-z_i^T\beta)=0$, so

\hat\beta=(X^TX)^{-1}X^TY = \begin{pmatrix}3&3\\3&5\end{pmatrix}^{-1} \begin{pmatrix}8\\12\end{pmatrix} = \begin{pmatrix}2/3\\2\end{pmatrix}.

The fitted residuals are $1/3,-2/3,1/3$. With unit working variance, $\hat H=X^TX$, and

\hat H^{-1}=\frac{1}{6} \begin{pmatrix}5&-3\\-3&3\end{pmatrix}, \qquad \hat M=\sum_{i=1}^3 (Y_i-z_i^T\hat\beta)^2z_iz_i^T = \begin{pmatrix}2/3&2/3\\2/3&8/9\end{pmatrix}.

Multiplying the matrices gives the empirical sandwich covariance:

\hat H^{-1}\hat M\hat H^{-1} = \begin{pmatrix} 7/54&-1/18\\ -1/18&1/18 \end{pmatrix}.

For comparison, the covariance calculated directly from the unit working variance is $\hat H^{-1}$, whose entries are $5/6,-1/2,1/2$. The difference illustrates why the empirical meat is a separate part of the calculation: even when the bread is the same, the covariance depends on how the observed score contributions vary. This three-observation example is an algebraic illustration, not a setting in which large-sample accuracy should be expected.

Common Pitfalls and Practical Interpretation

Using the working variance as the meat. That substitution assumes the working variance correctly describes the score variability. The empirical meat instead sums the outer products of fitted observation-level scores.
Using the total score as one meat contribution. The relevant quantity is $\sum_i\psi_i(\hat\beta)\psi_i(\hat\beta)^T$, not $U_n(\hat\beta)U_n(\hat\beta)^T$. At an exact solution the latter is zero.
Forgetting the transpose on the right bread. For a general estimating equation the covariance is $H^{-1}MH^{-T}$, not necessarily $H^{-1}MH^{-1}$.
Assuming robustness means protection from every error. The sandwich addresses covariance misspecification under suitable conditions. It does not make a misspecified mean target correct, nor does the independent-observation formula account for unmodelled dependence.
Interpreting large-sample justification as a small-sample guarantee. With few independent observations, the empirical meat can be unstable. The sandwich formula alone does not provide a finite-sample correction.

Key takeaway: Form one quasi-score vector per independent observation, estimate sensitivity with the negative derivative of the total score, and estimate score variability with the sum of individual outer products. Their sandwich estimates the large-sample covariance without requiring the working variance to be the true variance.

Check Your Understanding

Answer these questions using the definitions and derivations above.

For an estimating equation with $p$ parameters, what dimensions do the bread, meat, and sandwich covariance have?
Why does the empirical meat use observation-level score contributions rather than the total score at the fitted parameter?
Under what condition does the expected quasi-score meat equal the sensitivity for the mean-and-variance construction?
In the integral expansion, why is it valid to use an integrated Jacobian rather than assert the existence of one intermediate point that gives an exact matrix equality?
Give one form of model misspecification that the sandwich covariance does not correct.