Sandwich Inference for a Working GLS Estimator

Why Working GLS Needs Robust Inference

GLS uses a specified covariance structure to choose how observations are weighted. When that structure is correct, its model-based covariance formula can provide valid standard errors. But a working covariance matrix can be useful for estimation even when it is not the true covariance: under suitable conditions, the coefficient estimator can remain consistent while its conventional standard errors are wrong. Sandwich inference addresses that mismatch by estimating uncertainty from the variability of the estimator’s score contributions rather than relying on the working covariance as if it were correct.

This distinction matters in weighted regression. If the weights describe the relative precision of observations imperfectly, the fitted coefficients need not change much, but the reported uncertainty can. Robust standard errors are not automatically larger; they can be either larger or smaller than model-based standard errors.

Definition: Let $y$ be an $n \times 1$ outcome vector, $X$ an $n \times k$ full-column-rank design matrix, and $u=y-X\beta$ the regression error. Let $W$ be a specified positive-definite working precision matrix. The working GLS estimator is $$ \hat{\beta}_{W}=(X'WX)^{-1}X'Wy. $$ The bread is the matrix $X'WX$, the negative derivative of the estimating equation with respect to $\beta$. The meat is the covariance of that equation’s score contributions. The sandwich covariance places the meat between two inverse bread matrices.

Assumptions and the Sandwich Formula

For the basic result, take $X$ and $W$ as fixed, assume $X'WX$ is nonsingular, and suppose the conditional mean is correctly specified: $E[u\mid X]=0$. The true conditional covariance of $u$ is $\Sigma$, which need not equal the working covariance $W^{-1}$. For a sampling interpretation, observations or independent clusters must satisfy suitable regularity conditions, and no small number of units should dominate the estimator. These conditions support consistency and asymptotic normality; they do not make the sandwich an exact finite-sample distributional result.

The estimator solves the estimating equation $U(\beta)=X'W(y-X\beta)=0$. Since $E[u\mid X]=0$, its expectation at the true $\beta$ is zero. Its deviation from the true coefficient is exactly $$ \hat{\beta}_{W}-\beta=(X'WX)^{-1}X'Wu. $$ Consequently, conditional on $X$, its covariance is

\operatorname{Var}(\hat{\beta}_{W}\mid X) =(X'WX)^{-1}X'W\Sigma WX(X'WX)^{-1}.

This expression is the sandwich form: the inverse bread $(X'WX)^{-1}$ appears on each side, and the middle matrix $X'W\Sigma WX$ is the covariance of the estimating equation. If the working covariance is correct, $\Sigma=W^{-1}$, the middle reduces to $X'WX$, giving the familiar model-based covariance $(X'WX)^{-1}$. Under misspecification, that simplification is not justified.

Score Contributions and Robust Estimation

When $W$ is diagonal, write its diagonal entries as $w_i>0$, and let $x_i'$ be row $i$ of $X$. The estimating equation is a sum of observation-level contributions:

U(\beta)=\sum_{i=1}^{n}\psi_i(\beta), \qquad \psi_i(\beta)=w_i x_i(y_i-x_i'\beta).

Thus the score contribution for observation $i$ is the vector $w_i x_i u_i$. The bread is $A=X'WX=\sum_i w_i x_ix_i'$. If observations are independent but their variances can differ, the covariance of the sum is estimated by summing the outer products of the fitted score contributions. This yields the heteroskedasticity-robust estimator

\widehat{\operatorname{Var}}_{\mathrm{rob}}(\hat{\beta}_{W}) =(X'WX)^{-1} \left[ \sum_{i=1}^{n}w_i^2x_ix_i'\hat{u}_i^2 \right] (X'WX)^{-1}, \qquad \hat{u}_i=y_i-x_i'\hat{\beta}_{W}.

The squared weight is important: each score contribution contains one factor of $w_i$, so its outer product contains $w_i^2$. The residuals are used in place of unobserved errors. This estimator is often called the HC0 sandwich for diagonal-weight regression; small-sample adjustments are available, but do not remove the need for appropriate assumptions.

With correlated observations within independent clusters, the independent contributions are cluster scores, not individual observation scores. If $W$ is block diagonal by cluster, define $X_g$, $y_g$, and $W_g$ for cluster $g$. The cluster score is $\psi_g=X_g'W_g(y_g-X_g\beta)$, and the robust meat is $\sum_g\hat{\psi}_g\hat{\psi}_g'$. Treating correlated observations as independent rows can understate or otherwise misstate uncertainty. If dependence also crosses the chosen clusters, this cluster formula does not by itself account for it.

Worked Example: Deriving the Sandwich Covariance

Start with fixed $X$ and $W$, a nonsingular $X'WX$, and mean-zero errors with true covariance $\Sigma$, which may differ from $W^{-1}$. The estimating equation gives

\hat{\beta}_{W}-\beta=(X'WX)^{-1}X'Wu.

For any fixed matrix $C$, $\operatorname{Var}(Cu\mid X)=C\Sigma C'$. Here $C=(X'WX)^{-1}X'W$. Since $W$ is symmetric, $C' = WX(X'WX)^{-1}$, so

\begin{aligned} \operatorname{Var}(\hat{\beta}_{W}\mid X) &=(X'WX)^{-1}X'W\Sigma WX(X'WX)^{-1}. \end{aligned}

If $\Sigma=W^{-1}$, then $W\Sigma W=W$, and the middle becomes $X'WX$. The expression reduces to $(X'WX)^{-1}$. If instead $\Sigma\ne W^{-1}$, this cancellation does not occur. The two inverse bread matrices remain, while the true covariance determines the middle. This establishes why using the conventional GLS formula under misspecification can give incorrect standard errors.

Calculating Robust Standard Errors for Weighted Regression

For diagonal weights, the sandwich can be computed directly from the weighted fit and its residuals. The following hypothetical data are deliberately small to make the matrix arithmetic visible; with only four observations, asymptotic robust inference should be treated cautiously.

Worked Example: Heteroskedasticity-Robust Weighted Regression

Suppose a weighted regression with an intercept and one predictor uses the following hypothetical observations. The weights are positive working precision weights; the data-generating error variances are not assumed to equal their reciprocals.

$i$	$x_i$	$y_i$	$w_i$
1	0	1	1
2	1	4	2
3	2	3	1
4	3	5	2

Use $x_i'=(1,x_i)$. The weighted cross-product matrix is

X'WX= \begin{pmatrix} \sum_i w_i & \sum_i w_ix_i\\ \sum_i w_ix_i & \sum_i w_ix_i^2 \end{pmatrix} = \begin{pmatrix} 6&10\\ 10&24 \end{pmatrix}.

Also, $X'Wy=(22,44)'$, so the weighted estimate is

\hat{\beta}_{W} = \begin{pmatrix}6&10\\10&24\end{pmatrix}^{-1} \begin{pmatrix}22\\44\end{pmatrix} = \begin{pmatrix}2\\1\end{pmatrix}, \qquad (X'WX)^{-1}=\frac{1}{44} \begin{pmatrix}24&-10\\-10&6\end{pmatrix}.

The fitted values are $2,3,4,5$, giving residuals $-1,1,-1,0$. The score contributions $w_ix_i\hat{u}_i$, listed as column vectors, are $(-1,0)'$, $(2,2)'$, $(-1,-2)'$, and $(0,0)'$. Therefore the meat is

\sum_i w_i^2x_ix_i'\hat{u}_i^2 = \begin{pmatrix}1&0\\0&0\end{pmatrix} + \begin{pmatrix}4&4\\4&4\end{pmatrix} + \begin{pmatrix}1&2\\2&4\end{pmatrix} = \begin{pmatrix}6&6\\6&8\end{pmatrix}.

Substituting the bread and meat gives

\widehat{\operatorname{Var}}_{\mathrm{rob}}(\hat{\beta}_{W}) =\frac{1}{44^2} \begin{pmatrix}24&-10\\-10&6\end{pmatrix} \begin{pmatrix}6&6\\6&8\end{pmatrix} \begin{pmatrix}24&-10\\-10&6\end{pmatrix} = \begin{pmatrix} 0.7107&-0.2355\\ -0.2355&0.08678 \end{pmatrix}.

The robust standard errors are the square roots of the diagonal entries: approximately $0.8431$ for the intercept and $0.2946$ for the slope, with the slope measured in outcome units per unit of $x$. For comparison, treating the working covariance as correct gives $(X'WX)^{-1}$, with standard errors $\sqrt{24/44}=0.7385$ and $\sqrt{6/44}=0.3693$. Here the robust intercept standard error is larger and the robust slope standard error is smaller. The comparison illustrates that robust standard errors are not guaranteed to increase.

Cluster-Level Scores

The same bread-meat-bread structure applies when independent units are clusters, but each cluster contributes one score vector. This is needed when observations within a cluster may be correlated. The calculation below uses a simple intercept-only fit to make the aggregation explicit.

Worked Example: Cluster-Robust Intercept Standard Error

Suppose six hypothetical outcomes are arranged in three independent clusters of two observations each: $(1,3)$, $(2,4)$, and $(4,6)$. Use identity working precision, so the GLS estimate is the pooled mean. Assume the three clusters are independent, while allowing dependence within a cluster.

The estimate is $\hat{\beta}_{W}=20/6=3.333$. The bread is $X'WX=6$. For an intercept-only model, the score for cluster $g$ is the sum of its two residuals. The three cluster scores are $4-2(20/6)=-2.667$, $6-2(20/6)=-0.667$, and $10-2(20/6)=3.333$. Their sum is zero, as required by the fitted intercept equation.

The cluster meat is the sum of their squared scores:

(-2.667)^2+(-0.667)^2+(3.333)^2 \approx 18.67.

Thus the uncorrected cluster sandwich variance is $(168/9)/(6^2)\approx 0.5185$, and its standard error is $\sqrt{0.5185}\approx 0.7201$ outcome units. This is the cluster-level HC0 calculation. With only three clusters it is a weak basis for conventional large-sample inference; in practice, cluster-robust methods generally require a substantial number of independent clusters, and small-sample corrections or other methods may be appropriate.

Conditions, Interpretation, and Common Pitfalls

Robust covariance estimation does not repair every problem with a regression. The zero conditional mean assumption remains important for interpreting the coefficient estimator as targeting $\beta$; the sandwich addresses covariance misspecification, not omitted-variable bias or endogeneity. Likewise, a robust covariance formula does not create independence: the meat must reflect the actual independent sampling units. For diagonal-weight HC inference those units are observations; for cluster inference they are clusters.

The standard errors summarize estimated sampling variability of the coefficient estimates under the stated sampling assumptions. A confidence interval or test built from them still relies on the relevant large-sample approximation and should be interpreted accordingly. In particular, a numerically computed sandwich standard error does not ensure that a normal approximation is reliable in a tiny sample or with very few clusters.

Wrong: “The weights were specified, so the model-based GLS standard errors must be valid.” Correction: They have that justification only when the working covariance is correct (or under other conditions that validate the formula). Under misspecification, use a suitable sandwich estimate.
Wrong: “The robust meat is always $\sum_i x_ix_i'\hat{u}_i^2$.” Correction: With diagonal precision weights, the observation score includes $w_i$, so the meat includes $w_i^2$. With within-cluster dependence, aggregate score contributions at the cluster level instead.
Wrong: “Robust standard errors must be larger.” Correction: The bread and meat jointly determine each variance; as the weighted-regression example shows, some robust standard errors can be larger and others smaller than their model-based counterparts.
Wrong: “Robust standard errors make the coefficient estimate unbiased.” Correction: The sandwich changes the estimated covariance, not the fitted coefficients or the assumptions needed for their consistency.

Key takeaway: A working GLS estimator has sandwich covariance when its estimating equation is mean-zero and the sampling units are handled correctly. Identify the score contributions first, form their covariance as the meat, and place it between inverse bread matrices. Use observation scores for independent rows and cluster scores when dependence is contained within independent clusters.

Check Your Understanding

Try these questions using the definitions and formulas above.

For a diagonal working precision matrix, what is the score contribution of observation $i$, and why does its squared weight appear in the robust meat?
Under what condition does the sandwich covariance reduce to the conventional model-based GLS covariance?
Why can the robust standard error for a slope be smaller than its conventional GLS standard error?
If observations may be correlated within schools but schools are independent, what should count as a score contribution for robust inference?
Does using a sandwich covariance formula address a failure of the conditional mean assumption? Explain.

\(i\)	\(x_i\)	\(y_i\)	\(w_i\)
1	0	1	1
2	1	4	2
3	2	3	1
4	3	5	2