Estimating the Sandwich Covariance

Why Estimate Both Parts of the Sandwich?

The asymptotic distribution of an estimator defined by an estimating equation depends on two features: how the equation changes as the parameter moves, and how much its contributions vary across observations. The first determines how strongly the equation responds to parameter changes; the second describes its sampling fluctuation. A covariance estimate that replaces only one of these features can be misleading, especially when the estimating equation is not a likelihood score under a correctly specified model.

This tutorial develops plug-in estimates of both features. Throughout, let \(Z_1,\ldots,Z_n\) be independent and identically distributed observations, let \(\theta\in\mathbb{R}^p\), and let \(\psi(Z,\theta)\in\mathbb{R}^p\) be an estimating function. The estimator \(\hat\theta\) solves \(\sum_{i=1}^n\psi(Z_i,\hat\theta)=0\), or solves it up to an error negligible for its first-order asymptotic distribution.

Definition: Suppose \(E[\psi(Z,\theta_0)]=0\), and the derivative matrix exists and is integrable. Define the sensitivity matrix \(A=E[\partial\psi(Z,\theta_0)/\partial\theta']\) and the score-variability matrix \(B=E[\psi(Z,\theta_0)\psi(Z,\theta_0)']\). Since the estimating function has mean zero at \(\theta_0\), \(B=\operatorname{Var}(\psi(Z,\theta_0))\). If \(A\) is nonsingular, the sandwich covariance for the first-order approximation to \(\hat\theta\) is \(n^{-1}A^{-1}BA^{-T}\), where \(A^{-T}=(A^{-1})'\).

The derivative convention here uses the derivative of \(\psi\), not its negative. Thus \(A\) may have negative entries; that sign is retained in the formula, where it cancels appropriately. The asymptotic covariance of \(\sqrt{n}(\hat\theta-\theta_0)\) is \(A^{-1}BA^{-T}\), while the covariance scale for \(\hat\theta\) itself includes the factor \(1/n\).

The Empirical Sandwich Estimate

The natural plug-in procedure evaluates each observation’s derivative and estimating-function contribution at the fitted parameter. Define

\hat A=\frac{1}{n}\sum_{i=1}^n \frac{\partial\psi(Z_i,\hat\theta)}{\partial\theta'}, \qquad \hat B=\frac{1}{n}\sum_{i=1}^n \psi(Z_i,\hat\theta)\psi(Z_i,\hat\theta)'.

When \(\hat A\) is nonsingular, the estimated covariance matrix is

\widehat{\operatorname{Cov}}(\hat\theta) =\frac{1}{n}\hat A^{-1}\hat B\hat A^{-T}.

The product order matters for vector parameters: the rightmost factor is the transpose of the inverse on the left. This construction is often called the empirical sandwich estimator. The middle matrix uses observation-level estimating-function contributions, not a single aggregate contribution.

A common alternative uses centered contributions: \(\hat B_c=n^{-1}\sum_i(\psi(Z_i,\hat\theta)-\bar\psi)(\psi(Z_i,\hat\theta)-\bar\psi)'\), where \(\bar\psi=n^{-1}\sum_i\psi(Z_i,\hat\theta)\). For an exact root, \(\bar\psi=0\), so \(\hat B_c=\hat B\). For an approximate root, centering removes the empirical mean; under the conditions below, this distinction vanishes asymptotically if \(\bar\psi\xrightarrow{p}0\).

Why the Plug-In Matrices Are Consistent

A plug-in calculation is justified by consistency of \(\hat\theta\), laws of large numbers, and continuity of matrix inversion. Here is a precise sufficient set of conditions. Assume \(\hat\theta\xrightarrow{p}\theta_0\). In a neighborhood of \(\theta_0\), suppose \(\psi(z,\theta)\) is continuously differentiable in \(\theta\), and its derivative and the entries of \(\psi(z,\theta)\psi(z,\theta)'\) satisfy integrable envelope conditions that ensure uniform laws of large numbers on that neighborhood. Assume \(E[\psi(Z,\theta_0)]=0\), \(E\|\psi(Z,\theta_0)\|^2<\infty\), and \(A\) is nonsingular. These conditions imply the required continuity of the population expectations at \(\theta_0\).

Theorem: Under these conditions, \(\hat A\xrightarrow{p}A\) and \(\hat B\xrightarrow{p}B\). Consequently, with probability tending to one \(\hat A\) is invertible, and \(\hat A^{-1}\hat B\hat A^{-T}\xrightarrow{p}A^{-1}BA^{-T}\). Thus the plug-in covariance estimate consistently estimates the asymptotic covariance divided by \(n\).

To see the first convergence, the uniform law of large numbers controls the empirical derivative uniformly near \(\theta_0\). Since \(\hat\theta\xrightarrow{p}\theta_0\), with probability tending to one the fitted value lies in that neighborhood, and therefore

\left\|\hat A-A\right\| \leq \sup_{\theta\text{ near }\theta_0} \left\| \frac{1}{n}\sum_{i=1}^n \frac{\partial\psi(Z_i,\theta)}{\partial\theta'} - E\left[\frac{\partial\psi(Z,\theta)}{\partial\theta'}\right] \right\| + \left\| E\left[\frac{\partial\psi(Z,\hat\theta)}{\partial\theta'}\right]-A \right\|.

The first term converges in probability to zero by the uniform law; the second does so by consistency and continuity. Hence \(\hat A\xrightarrow{p}A\). Apply the same argument to the matrix-valued function \(\psi(z,\theta)\psi(z,\theta)'\), using its assumed integrable envelope. This gives \(\hat B\xrightarrow{p}E[\psi(Z,\theta_0)\psi(Z,\theta_0)']=B\).

Nonsingularity of \(A\) and continuity of the inverse map imply that \(\hat A\) is nonsingular with probability tending to one and \(\hat A^{-1}\xrightarrow{p}A^{-1}\). Transposition is continuous as well. The continuous mapping theorem and the product rule for convergence in probability now give the asserted sandwich convergence. Multiplying by \(1/n\) yields a consistent estimate of the covariance on the scale of \(\hat\theta\). \(\square\)

Worked Calculations

Worked Example: A Scalar Mean Equation

Suppose the estimating function is \(\psi(Z,\theta)=Z-\theta\), and the observed values are \(2,4,6\). The estimating equation is \( (2-\theta)+(4-\theta)+(6-\theta)=0\), so \(\hat\theta=4\).

The derivative is \(\partial\psi/\partial\theta=-1\) for every observation. Thus \(\hat A=(-1-1-1)/3=-1\). The three fitted contributions are \(-2,0,2\), so

\hat B=\frac{(-2)^2+0^2+2^2}{3}=\frac{8}{3}.

The estimated variance and standard error are therefore

\widehat{\operatorname{Cov}}(\hat\theta) =\frac{1}{3}(-1)^{-1}\left(\frac{8}{3}\right)(-1)^{-1} =\frac{8}{9}, \qquad \widehat{\operatorname{SE}}(\hat\theta)=\sqrt{\frac{8}{9}}\approx0.943.

The two negative scalar inverse factors multiply to a positive value. This plug-in estimate uses denominator \(n\) in \(\hat B\); it is an asymptotic covariance estimate, not a claim of finite-sample unbiasedness.

Worked Example: A Weighted Estimating Equation

Consider \(\psi((x,y),\theta)=x^2(y-\theta)\) and the four observations \((x,y)=(1,2),(1,4),(2,5),(2,8)\). Solving the empirical equation gives

\hat\theta =\frac{\sum_i x_i^2y_i}{\sum_i x_i^2} =\frac{2+4+20+32}{1+1+4+4} =\frac{58}{10}=5.8.

The fitted contributions are \(-3.8,-1.8,4(5-5.8)=-3.2,\) and \(4(8-5.8)=8.8\). They sum to zero, as required. The derivative of the \(i\)th contribution is \(-x_i^2\), so \(\hat A=-(1+1+4+4)/4=-2.5\). The score-variability estimate is

\hat B =\frac{(-3.8)^2+(-1.8)^2+(-3.2)^2+(8.8)^2}{4} =\frac{14.44+3.24+10.24+77.44}{4} =\frac{105.36}{4}=26.34.

Hence

\widehat{\operatorname{Cov}}(\hat\theta) =\frac{26.34}{4(-2.5)^2} =\frac{26.34}{25}=1.0536, \qquad \widehat{\operatorname{SE}}(\hat\theta)=\sqrt{1.0536}\approx1.026.

This example illustrates why the observed variability of the estimating contributions belongs in the middle of the sandwich: the derivative alone does not determine the estimator’s uncertainty.

Worked Example: A Two-Parameter Wald Region

Let \(Z\) be a two-dimensional observation and use the mean equation \(\psi(Z,\theta)=Z-\theta\), with observed vectors \((0,0),(2,0),(0,2),(2,2)\). Their mean is \(\hat\theta=(1,1)'\). The fitted contributions are \((-1,-1)',(1,-1)',(-1,1)',(1,1)'\). Since \(\hat A=-I_2\), and averaging the four outer products gives \(\hat B=I_2\), the covariance estimate is

\widehat{\operatorname{Cov}}(\hat\theta) =\frac{1}{4}(-I_2)^{-1}I_2(-I_2)^{-T} =\frac{1}{4}I_2.

If the estimator satisfies the usual asymptotic normality conditions and the covariance estimate is consistent, a \(95\%\) asymptotic Wald confidence region for \(\theta\) is

\left\{\theta: (\hat\theta-\theta)' \widehat{\operatorname{Cov}}(\hat\theta)^{-1} (\hat\theta-\theta) \leq \chi^2_{2,0.95} \right\}.

Here the inverse covariance estimate is \(4I_2\), and \(\chi^2_{2,0.95}\approx5.991\). Thus the region is \(4[(1-\theta_1)^2+(1-\theta_2)^2]\leq5.991\), or a disk centered at \((1,1)'\) with squared radius \(5.991/4\). This is an asymptotic confidence region, not an exact finite-sample guarantee.

Wald Inference and Common Pitfalls

More generally, suppose \(\sqrt n(\hat\theta-\theta_0)\xrightarrow{d}N(0,V)\), where \(V=A^{-1}BA^{-T}\) is positive definite, and \(\hat V=\hat A^{-1}\hat B\hat A^{-T}\xrightarrow{p}V\). Slutsky’s theorem gives the Wald statistic

n(\hat\theta-\theta_0)'\hat V^{-1}(\hat\theta-\theta_0) \xrightarrow{d}\chi^2_p.

Indeed, the standardized vector converges to \(N(0,I_p)\), and the sum of squares of its \(p\) coordinates has a \(\chi^2_p\) distribution. A confidence region with asymptotic coverage \(1-\alpha\) replaces the right side of the inequality in the example with \(\chi^2_{p,1-\alpha}\). For a single component, the estimated standard error is the square root of the corresponding diagonal entry of \(\hat V/n\).

Several errors are especially easy to make:

Dropping the factor \(1/n\). The matrix \(\hat A^{-1}\hat B\hat A^{-T}\) estimates the covariance of the \(\sqrt n\)-scaled estimator. The covariance estimate for \(\hat\theta\) includes division by \(n\).
Using the wrong matrix order. For vectors, the sandwich is \(A^{-1}BA^{-T}\), not \(A^{-1}BA^{-1}\) unless the relevant matrices happen to be symmetric.
Replacing score variability with the derivative. The sensitivity and variability are distinct quantities. They coincide in particular likelihood settings under additional information identities, but they need not coincide for a general estimating equation.
Changing derivative signs without checking the formula. Defining sensitivity as the negative derivative is also common, but the convention and sandwich formula must be changed together. With the derivative convention used here, both inverse factors use \(A\).
Treating a Wald approximation as exact. The chi-square calibration is justified asymptotically under normality and consistency conditions; it does not automatically describe small-sample coverage.

The practical value of the sandwich estimate is that it uses the observed estimating contributions to quantify variability while separately accounting for the equation’s sensitivity to parameter changes. Reporting both the covariance matrix and the standard errors derived from its diagonal makes the uncertainty calculation auditable.

Check Your Understanding

Work through these questions to check the construction and interpretation of the sandwich estimate.

For a vector estimating function, why does the covariance formula contain \(A^{-T}\) on the right?
If \(\hat A\) estimates the derivative matrix and \(\hat B\) estimates score variability, which matrix estimates the covariance of \(\sqrt n(\hat\theta-\theta_0)\), and where does the sample-size factor enter for \(\hat\theta\)?
When does the centered empirical score-variability estimate equal the uncentered estimate exactly?
In the weighted-equation example, verify that the four fitted contributions sum to zero and that their squared values sum to \(105.36\).
State the asymptotic conditions needed to interpret a Wald confidence region based on the sandwich covariance estimate.