Deriving Sandwich Covariance for Quasi-Estimators

Why Two Matrices Determine Uncertainty

An estimating equation identifies a parameter by requiring a sum of contributions to equal zero. For a quasi-estimator, those contributions are built from the specified mean and variance relationship; they need not be likelihood scores. The estimator’s uncertainty depends on two distinct features: how sharply the equation changes when the parameter changes, and how variable the equation is at the true parameter.

The first feature is called sensitivity; the second is score variability (or variability, for a general estimating equation). Their combination gives the sandwich covariance. The derivation is a local linearization of the estimating equation. In multiple dimensions, that linearization must use an integral of Jacobians along the segment between the true parameter and the estimator: a single intermediate Jacobian does not, in general, give an exact multivariate Taylor formula.

Estimating Equations and the Two Matrices

Definition: Let $\beta\in\mathbb{R}^p$ be a parameter and let $U_n(\beta)=\sum_{i=1}^n\psi_{ni}(\beta)$ be an estimating function, where each contribution $\psi_{ni}(\beta)\in\mathbb{R}^p$. An estimator $\hat\beta$ solves $U_n(\hat\beta)=0$. At the true parameter $\beta_0$, the sensitivity matrix is the probability limit of the negative normalized Jacobian, when that limit exists. The score-variability matrix is the limiting normalized covariance of $U_n(\beta_0)$, when that limit exists.

More precisely, for independent contributions with $E[\psi_{ni}(\beta_0)]=0$, define $B_n=n^{-1}\operatorname{Var}(U_n(\beta_0))=n^{-1}\sum_{i=1}^n\operatorname{Var}(\psi_{ni}(\beta_0))$. If $B_n$ converges to $B$, then $B$ is the score-variability matrix. For the sensitivity, the derivation below uses the normalized Jacobian averaged along the segment from $\beta_0$ to $\hat\beta$; the required condition is that its negative converges in probability to a nonsingular matrix $A$.

In a GLM quasi-score, a typical contribution has the form $\psi_i(\beta)=x_i\mu_i'(\beta)\{Y_i-\mu_i(\beta)\}/\{\phi V(\mu_i(\beta))\}$, where $\mu_i'(\beta)$ denotes the derivative of the mean with respect to the linear predictor, $V$ is the specified variance function, and $\phi$ is the dispersion parameter. The covariance derivation applies to this equation without requiring a full probability likelihood.

The Sandwich Derivation

Theorem: Suppose $p$ is fixed, $\beta_0$ is an interior point of an open parameter set, and $U_n$ is continuously differentiable on a neighborhood containing the segment from $\beta_0$ to a consistent root $\hat\beta$. Suppose $$ -\frac{1}{n}\int_0^1 \frac{\partial U_n\{\beta_0+t(\hat\beta-\beta_0)\}}{\partial\beta^\mathsf{T}}\,dt \xrightarrow{p} A, $$ where $A$ is nonsingular. Also suppose $n^{-1/2}U_n(\beta_0)\xrightarrow{d}N(0,B)$, with $B$ a finite covariance matrix. Then $$ \sqrt{n}(\hat\beta-\beta_0) =A^{-1}\frac{U_n(\beta_0)}{\sqrt{n}}+o_p(1) \xrightarrow{d}N\!\left(0,A^{-1}BA^{-\mathsf{T}}\right). $$

Proof. Write $\Delta_n=\hat\beta-\beta_0$. For each component of $U_n$, apply the one-dimensional fundamental theorem of calculus to the function $t\mapsto U_n(\beta_0+t\Delta_n)$. The chain rule gives the exact vector identity

U_n(\hat\beta)-U_n(\beta_0) = \left[ \int_0^1 \frac{\partial U_n(\beta_0+t\Delta_n)}{\partial\beta^\mathsf{T}}\,dt \right]\Delta_n.

This formula uses an averaged Jacobian matrix, not a Jacobian at one unspecified intermediate point. Since $U_n(\hat\beta)=0$, define

C_n=-\frac{1}{n}\int_0^1 \frac{\partial U_n(\beta_0+t\Delta_n)}{\partial\beta^\mathsf{T}}\,dt.

The exact identity becomes $0=U_n(\beta_0)-nC_n\Delta_n$, hence

\sqrt{n}\Delta_n=C_n^{-1}\frac{U_n(\beta_0)}{\sqrt{n}}.

Because $C_n\xrightarrow{p}A$ and $A$ is nonsingular, $C_n$ is invertible with probability tending to one and $C_n^{-1}\xrightarrow{p}A^{-1}$. Slutsky’s theorem gives the stated limiting normal distribution. Its covariance is $A^{-1}BA^{-\mathsf{T}}$: applying the linear map $A^{-1}$ to a vector with covariance $B$ transforms covariance to $A^{-1}B(A^{-1})^\mathsf{T}$. This establishes both the linear representation and the sandwich covariance. $\square$

For independent, mean-zero contributions, the covariance of their sum is the sum of their covariances, which explains the definition of $B_n$. The central limit assumption in the theorem is a separate requirement; in applications it is usually justified by a suitable central limit theorem. The theorem does not assert that a covariance formula is valid without the stated consistency, Jacobian convergence, nonsingularity, and distributional limit conditions.

Worked Examples

Worked Example: A Weighted Scalar Estimating Equation

Suppose there are $n$ independent observations, with $n$ even. Half have distribution $Y_i\sim N(\beta_0,1)$, and half have distribution $Y_i\sim N(\beta_0,4)$. Use weight $w_i=1$ for the first group and $w_i=1/2$ for the second, and define $\psi_i(\beta)=w_i(Y_i-\beta)$. Find the sandwich variance.

The estimating equation is $\sum_i w_i(Y_i-\hat\beta)=0$. Its derivative with respect to $\beta$ is $-\sum_iw_i$, so the normalized negative derivative is $A_n=n^{-1}\sum_iw_i=(1+1/2)/2=3/4$. This is constant in $\beta$, so $A=3/4$.

Each contribution has mean zero at $\beta_0$. Its variance is $w_i^2\operatorname{Var}(Y_i)$. Thus $B_n=n^{-1}\sum_iw_i^2\operatorname{Var}(Y_i)=\{1^2(1)+(1/2)^2(4)\}/2=1$, so $B=1$. The sandwich variance per observation is

A^{-1}BA^{-1} =\left(\frac{4}{3}\right)(1)\left(\frac{4}{3}\right) =\frac{16}{9}.

Therefore the asymptotic variance of $\hat\beta$ is $16/(9n)$. In fact, here this is also its exact variance: the weighted numerator has variance $n$, and the squared denominator is $(3n/4)^2$, giving $n/(9n^2/16)=16/(9n)$. The example separates sensitivity from variability: knowing the derivative alone does not determine the estimator’s uncertainty.

Worked Example: Multiplying the Matrices

Consider a two-dimensional estimating equation whose sensitivity and variability limits are $A=\begin{pmatrix}2&1\\1&2\end{pmatrix}$ and $B=\begin{pmatrix}5&1\\1&2\end{pmatrix}$. For example, these matrices arise from contributions $\psi_i(\beta)=Z_i-A(\beta-\beta_0)$, where the independent random vectors $Z_i\sim N_2(0,B)$. Compute the sandwich covariance.

The determinant of $A$ is $2(2)-1(1)=3$, so $A^{-1}=\frac13\begin{pmatrix}2&-1\\-1&2\end{pmatrix}$. Both matrices are symmetric, so $A^{-\mathsf{T}}=A^{-1}$. Multiplying first on the left gives

\begin{pmatrix}2&-1\\-1&2\end{pmatrix} \begin{pmatrix}5&1\\1&2\end{pmatrix} = \begin{pmatrix}9&0\\-3&3\end{pmatrix}.

Multiplying by the remaining factor and including the two factors of $1/3$ gives

A^{-1}BA^{-\mathsf{T}} =\frac{1}{9} \begin{pmatrix}9&0\\-3&3\end{pmatrix} \begin{pmatrix}2&-1\\-1&2\end{pmatrix} =\frac{1}{9} \begin{pmatrix}18&-9\\-9&9\end{pmatrix} = \begin{pmatrix}2&-1\\-1&1\end{pmatrix}.

This matrix is a valid covariance matrix: its diagonal entries are positive and its determinant is $2(1)-(-1)^2=1>0$. The off-diagonal entries are negative, indicating negative asymptotic covariance between the two estimator components in this example.

Worked Example: A Quasi-Poisson Intercept

Let $Y_1,\ldots,Y_n$ be independent negative-binomial counts with mean $\mu=4$ and size $r=4$, using the parameterization $\operatorname{Var}(Y_i)=\mu+\mu^2/r$. The variance is $4+16/4=8$. Fit an intercept-only log-mean equation with working dispersion $\phi=1$, so $\mu(\beta)=e^\beta$ and $\beta_0=\log 4$.

For the quasi-Poisson variance function $V(\mu)=\mu$, the quasi-score contribution is $\psi_i(\beta)=\mu(\beta)\{Y_i-\mu(\beta)\}/\{\phi V(\mu(\beta))\}=Y_i-e^\beta$. Its derivative is $-e^\beta$, so the sensitivity per observation at $\beta_0$ is $A=4$. At $\beta_0$, the contribution is $Y_i-4$, whose variance is $8$, so $B=8$.

The sandwich covariance for $\sqrt{n}(\hat\beta-\beta_0)$ is therefore

A^{-1}BA^{-1} =\frac{1}{4}(8)\frac{1}{4} =\frac{1}{2}.

Thus the estimator’s asymptotic variance is $1/(2n)$. For comparison, using the working Poisson variability $B=4$ would give $4/(16n)=1/(4n)$, only half as large. The actual count variance, not just the equation’s derivative, matters for uncertainty.

Why Rescaling the Equation Does Not Change the Sandwich

An estimating equation can be multiplied by a fixed nonsingular matrix without changing its roots. The sandwich covariance is likewise invariant under this transformation. To verify this, let $C$ be a nonsingular $p\times p$ matrix and replace each contribution by $C\psi_{ni}(\beta)$. The new sensitivity and variability are $A_*=CA$ and $B_*=CBC^\mathsf{T}$. Therefore

A_*^{-1}B_*A_*^{-\mathsf{T}} =(CA)^{-1}(CBC^\mathsf{T})(CA)^{-\mathsf{T}} =A^{-1}C^{-1}CBC^\mathsf{T}C^{-\mathsf{T}}A^{-\mathsf{T}} =A^{-1}BA^{-\mathsf{T}}.

The cancellation confirms that the sandwich describes the estimator’s uncertainty, not an arbitrary scale or coordinate choice for the equation. This result assumes the multiplier is fixed and nonsingular; it does not claim invariance under every data-dependent transformation.

Common Mistake: Using One Intermediate Jacobian

A tempting multivariate argument writes the change in $U_n$ as its Jacobian at a single intermediate parameter, multiplied by the parameter change. That is not generally an exact vector identity when $p>1$: the scalar mean value theorem does not extend in that form to vector-valued functions. The correction is the integral Taylor identity proved above, which averages the Jacobian along the whole segment. The limit of the negative normalized integral Jacobian is the sensitivity matrix used in the sandwich derivation.

A second related mistake is to report $A^{-1}$ as the covariance. Sensitivity measures the response of the equation to parameter changes; it does not measure the random fluctuation of the equation. The variability matrix $B$ must also be included. Only in special settings where the matrices satisfy additional identities can the sandwich simplify.

Key takeaway: Under consistency, convergence of the averaged negative Jacobian to a nonsingular $A$, and a central limit theorem for the estimating function with covariance limit $B$, the quasi-estimator has asymptotic covariance $A^{-1}BA^{-\mathsf{T}}/n$. Sensitivity determines how equation fluctuations translate into parameter fluctuations; score variability determines the size and structure of those fluctuations.

Check Your Understanding

Use the definitions and derivation above to answer each question.

Why does a vector-valued Taylor expansion use an integral of Jacobians rather than necessarily one Jacobian at an intermediate point?
For independent, mean-zero estimating-function contributions, how is the covariance of their sum related to the individual covariance matrices?
If $A$ is nonsingular and $B$ is a covariance matrix, write the asymptotic covariance of $\sqrt{n}(\hat\beta-\beta_0)$.
In the weighted scalar example, what would the sandwich variance per observation be if $B$ were incorrectly set to $4$ while $A=3/4$?
Why does multiplying an estimating equation by a fixed nonsingular matrix leave its sandwich covariance unchanged?