Sandwich Estimator

Sandwich Estimator for Robust Variance in Health Economic Regression

The sandwich estimator is the formula behind the robust and cluster-robust standard errors that statistical software reports after a regression of patient costs or outcomes. It leaves the fitted coefficients alone and re-estimates only their covariance matrix, replacing the variance the model assumes with the observed spread of each patient's, or each cluster's, contribution to the estimating equations. Health economists use it because patient-level cost data are skewed and heteroskedastic, and because patients in cluster randomised trials and repeated measures studies are not independent.

Bread, meat and the general formula

Most regression estimators used in health economics, including least squares, maximum likelihood and generalised estimating equations, choose the parameters that set a sum of estimating functions to zero. Zeileis describes the large-sample covariance of such an M-estimator as a slice of meat between two slices of bread. The bread is the inverse of the expected derivative of the estimating function and the meat is the variance of the estimating function itself.

$$\widehat{V}(\hat{\theta}) = \hat{A}^{-1}\,\hat{B}\,\hat{A}^{-1}, \qquad \hat{A} = -\sum_{i=1}^{n} \frac{\partial \psi_i(\hat{\theta})}{\partial \theta^{\mathrm{T}}}, \qquad \hat{B} = \sum_{i=1}^{n} \psi_i(\hat{\theta})\,\psi_i(\hat{\theta})^{\mathrm{T}}$$

where $\hat{\theta}$ is the vector of estimated parameters, $n$ is the number of independent units, $\psi_i(\hat{\theta})$ is the estimating function (for maximum likelihood, the score) contribution of unit $i$ at the estimates, $\hat{A}^{-1}$ is the bread and $\hat{B}$ is the meat. Freedman writes the Huber sandwich estimator in this form, with the bread built from the second derivative of the log likelihood, and calls the square roots of its diagonal elements robust or Huber-White standard errors.

When a model fitted by maximum likelihood is correctly specified, the meat and the matrix inside the bread both estimate the Fisher information, so the sandwich collapses to the usual model-based covariance. Zeileis notes that robustness comes from estimating the meat under weaker assumptions, and that the estimating functions must still be correctly specified so that $\hat{\theta}$ is consistent. The result is a variance-covariance matrix whose diagonal gives the robust standard error of each coefficient. The Stata User's Guide records that Huber and White developed the calculation independently and that the sandwich name comes from how the formula looks on the page.

Heteroskedasticity-consistent versions for least squares

For linear regression the bread is the inverse of $X^{\mathrm{T}}X$ and the meat weights each patient's row of regressors by a function of that patient's residual. Zeileis sets out the common variants as different choices of weight.

$$\widehat{V}_{HC}(\hat{\beta}) = (X^{\mathrm{T}}X)^{-1} X^{\mathrm{T}} \hat{\Omega} X (X^{\mathrm{T}}X)^{-1}, \qquad \hat{\Omega} = \operatorname{diag}(\omega_1, \ldots, \omega_n)$$

where $\hat{\beta}$ is the vector of least squares coefficients, $X$ is the $n \times k$ matrix of regressors, $k$ is the number of coefficients and $\omega_i$ is the weight for patient $i$. With residual $\hat{u}_i$ and hat value $h_i$ (the diagonal element of the hat matrix), HC0 sets $\omega_i = \hat{u}_i^2$, HC1 multiplies that by $n/(n-k)$, HC2 divides it by $(1-h_i)$ and HC3 divides it by $(1-h_i)^2$.

Zeileis attributes HC0 to White, on asymptotic arguments, and HC1 to HC3 to MacKinnon and White, for small samples. In their sampling experiments a version based on the jackknife performed best in small samples, and they suggested it may be wise to use it even when no heteroskedasticity has been detected. Zeileis summarises later simulation evidence as favouring HC3 in small samples, because HC3 gives less weight to influential observations.

Software defaults differ. The Stata User's Guide gives $n/(n-k)$ as the finite-sample adjustment for linear regression with the robust option, which is HC1, while the vcovHC function in the R sandwich package uses HC3 by default. A bare report of "robust standard errors" leaves the calculation unclear.

Use in health care cost regressions

Jones, reviewing econometric models for health care costs, notes that individual-level cost data typically show a high degree of heteroskedasticity and that the norm is to use robust standard errors for inference. Poisson regression applied to costs remains consistent when the mean is correctly specified, even if the conditional variance is wrong, and he notes that robust standard errors from the sandwich estimator are then used in place of the maximum likelihood ones. Generalised linear models such as the gamma model share this quasi-likelihood property.

Robust and conventional standard errors can be close. In Jones's illustrative analysis of 2,955 older people with positive medical expenditure, tests showed strong heteroskedasticity in the linear cost model, yet Huber-White robust estimates did little to change the standard errors.

Clusters, repeated measures and generalised estimating equations

When patients are grouped within practices or hospitals, or observed repeatedly over time, the meat is built from cluster totals rather than individual contributions. Freedman describes this as grouping the score terms by cluster before forming the outer product, with the two outer factors unchanged. The Stata User's Guide notes that the clustered form requires observations to be independent only across clusters.

$$\widehat{V}_{CR}(\hat{\beta}) = (X^{\mathrm{T}}X)^{-1} \left( \sum_{g=1}^{G} X_g^{\mathrm{T}} \hat{u}_g \hat{u}_g^{\mathrm{T}} X_g \right) (X^{\mathrm{T}}X)^{-1}$$

where $G$ is the number of clusters, $X_g$ holds the regressor rows for cluster $g$ and $\hat{u}_g$ is the vector of least squares residuals for that cluster. Cameron and Miller stress that each cluster's own residual outer product is a very poor estimate of that cluster's error covariance, and that the estimator works only through averaging across clusters. Its accuracy therefore depends on the number of clusters, not the number of patients.

Cameron and Miller recommend at least a small-sample scaling, such as the factor $\frac{G}{G-1} \times \frac{N-1}{N-K}$ that Stata applies, where $N$ is the number of observations and $K$ the number of regressors, together with critical values from a t distribution with $G-1$ degrees of freedom. Even then, Wald tests generally over-reject with few clusters, and a test with nominal size 0.05 may have a true size of 0.10. Cameron and Miller give no fixed threshold, noting that "few" may mean fewer than 20 clusters in some balanced designs, fewer than 50 in others and more when cluster sizes are unbalanced.

Mancl and DeRouen call this the robust covariance estimator of generalised estimating equations and proposed a bias-corrected version. In their simulations with 10 to 40 subjects (the clusters), tests from the robust and jackknife estimators could reach two to three times the nominal size, while the corrected estimator stayed close to it.

Evidence from cluster randomised cost-effectiveness analyses

Gomes and colleagues simulated cost-effectiveness analyses of cluster randomised trials and compared the coverage of 95% confidence intervals for incremental net benefit. In the base case, with 20 clusters per arm, seemingly unrelated regression that ignored the cluster effect gave coverage of 0.89, against 0.940 with a robust standard error and 0.933 for a generalised estimating equation with a robust standard error. With 3 clusters per arm, coverage for the two robust methods fell to 0.856 and 0.841.

The authors concluded that 8 to 15 clusters per arm are needed for the asymptotic assumptions to hold, more when cluster sizes are highly imbalanced, and did not recommend either robust method with few clusters. Multilevel (hierarchical) models and a two-stage bootstrap with shrinkage correction performed well throughout.

Worked example with an illustrative two-arm cost regression

The numbers below are illustrative. A small trial compares a new service (T, 4 patients) with usual care (U, 8 patients). Costs in arm T are GBP 1,000, GBP 1,500, GBP 2,500 and GBP 7,000. In arm U they are GBP 800, GBP 900, GBP 1,000, GBP 1,000, GBP 1,100, GBP 1,100, GBP 1,200 and GBP 1,300. Cost is regressed on an intercept and a treatment indicator.

1. Coefficients. Mean cost is GBP 3,000 in arm T and GBP 1,050 in arm U, so the intercept is 1,050 and the treatment coefficient, the incremental cost, is 3,000 minus 1,050, or 1,950. These point estimates stay the same under every variance option below.

2. Residuals. In arm T the residuals are -2,000, -1,500, -500 and 4,000, whose squares sum to 22,500,000. In arm U they lie between -250 and 250, and their squares sum to 180,000.

3. Bread and meat. With 12 patients, 4 of them treated, $X^{\mathrm{T}}X$ has rows (12, 4) and (4, 4), determinant 12 × 4 minus 4 × 4, or 32, so the bread has rows (0.125, -0.125) and (-0.125, 0.375). Summing $\hat{u}_i^2 x_i x_i^{\mathrm{T}}$, where $x_i$ is patient $i$'s column of regressors (1 and the treatment indicator), gives an HC0 meat with rows (22,680,000, 22,500,000) and (22,500,000, 22,500,000). For this design the treatment element of the sandwich reduces to a simple form.

$$\widehat{V}_{HC0}(\hat{\beta}_1) = \frac{\sum_{i \in T} \hat{u}_i^2}{n_T^2} + \frac{\sum_{i \in U} \hat{u}_i^2}{n_U^2} = \frac{22{,}500{,}000}{16} + \frac{180{,}000}{64} = 1{,}409{,}062.5$$

where $\hat{\beta}_1$ is the treatment coefficient, $n_T = 4$ and $n_U = 8$ are the arm sizes and the sums run over each arm's residuals. The HC0 standard error is the square root, GBP 1,187.0.

4. Comparison of versions. The conventional estimate pools the residual variance, 22,680,000 / 10 = 2,268,000, and multiplies it by 0.375, giving 850,500. Every patient in arm T has hat value 0.25 and every patient in arm U has hat value 0.125, which makes HC2 equal to the unequal-variance (Welch) standard error.

5. Interpretation. The conventional formula borrows the small spread of the larger usual-care arm and applies it to the smaller, far more variable treatment arm, so it understates uncertainty about the incremental cost. The sandwich standard errors are larger than the conventional one by a factor between 1.29 for HC0 and 1.72 for HC3. With 12 patients and one extreme cost the versions disagree widely, so the example shows the mechanics only and supports no firm conclusion about the cost difference.

Assumptions, choices and neighbouring methods

The estimator relaxes assumptions about the variance and, in its clustered form, about independence, but it still needs a correctly specified mean or estimating function and enough independent units to average over. Cameron and Miller advise clustering at a level broad enough to capture correlation whenever both regressors and errors might be correlated within the cluster, for example at practice level when treatment is assigned by practice.

A multilevel model writes the within-cluster correlation into the model itself and a bootstrap standard error obtains precision by resampling, whereas the sandwich corrects the variance of an existing estimator. Its covariance matrix feeds the Wald test for one coefficient or several.

Misreadings and limits

Most problems arise when the estimator is credited with more than it does. Others arise when the version used goes unreported.

1. Treating robust standard errors as a cure for a wrong model. Freedman argues that if the model is nearly correct the usual standard errors are nearly correct too, and that if it is seriously wrong the sandwich may help on the variance side while the parameters being estimated are likely to be meaningless. The algorithm also ignores bias.
2. Expecting robust standard errors always to be larger. Cameron and Miller note that White heteroskedasticity-robust standard errors are in practice sometimes larger and sometimes smaller than the default.
3. Leaving the method unreported. An analysis should state the variance version, the cluster variable and the number of clusters.

Sources

• Cameron AC, Miller DL. A practitioner's guide to cluster-robust inference. Journal of Human Resources. 2015;50(2):317-372. https://doi.org/10.3368/jhr.50.2.317
• Freedman DA. On the so-called "Huber sandwich estimator" and "robust standard errors". The American Statistician. 2006;60(4):299-302. https://doi.org/10.1198/000313006X152207
• Gomes M, Ng ES, Grieve R, Nixon R, Carpenter J, Thompson SG. Developing appropriate methods for cost-effectiveness analysis of cluster randomized trials. Medical Decision Making. 2012;32(2):350-361. https://doi.org/10.1177/0272989X11418372
• Jones AM. Models for health care. HEDG Working Paper 10/01. University of York; January 2010. https://www.york.ac.uk/media/economics/documents/herc/wp/10_01.pdf
• MacKinnon JG, White H. Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties. Journal of Econometrics. 1985;29(3):305-325. https://doi.org/10.1016/0304-4076(85)90158-7
• Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics. 2001;57(1):126-134. https://doi.org/10.1111/j.0006-341X.2001.00126.x
• StataCorp. Stata User's Guide, Release 19. Chapter 20, section 20.22, Obtaining robust variance estimates. Accessed 28 October 2026. https://www.stata.com/manuals/u20.pdf
• Zeileis A. Econometric computing with HC and HAC covariance matrix estimators. Journal of Statistical Software. 2004;11(10). CRAN vignette version. https://cran.r-project.org/web/packages/sandwich/vignettes/sandwich.pdf
• Zeileis A. Object-oriented computation of sandwich estimators. Journal of Statistical Software. 2006;16(9). CRAN vignette version. https://cran.r-project.org/web/packages/sandwich/vignettes/sandwich-OOP.pdf