Robust Standard Errors: What the Sandwich Fixes and What It Cannot

อ่านฉบับภาษาไทย (Thai version)

Abstract

Many regression tables report robust standard errors without saying what they repair. A regression has a mean model, for the average outcome, and a variance model, for the scatter around it. When the scatter changes with covariates (heteroskedasticity), the least-squares slope stays unbiased if the mean model is right, but the model-based standard error can mislead. The sandwich estimator rebuilds the variance from each person's squared residual. The change can go either way: the robust standard error is larger when the large residual variance sits at high-leverage covariate values, far from the mean, and can be smaller otherwise. In 2,000 simulated adults whose blood pressure scatters more with age, the robust standard error was 1.13 times the model-based one for age. For an under-40 indicator it was 0.74 times. In repeated studies of 200 simulated adults, robust intervals contained the simulated truth about 95% of the time. This article shows how to compute, choose and report a sandwich standard error, and what it cannot fix.

🔍 Click to enlarge Click anywhere to close

A robust standard error nobody can explain

A research team is drafting a paper from a community cohort of 2,000 adults aged 30 to 80. One question is how systolic blood pressure (SBP), in mmHg, rises with age. The cohort is fictional and its data are simulated.

A scatter plot shows that blood pressure varies far more among older adults than among younger ones. Every regression table in the draft reports a robust standard error. Nobody on the team can say what it repaired, or whether it moved the confidence interval wider or narrower.

Unequal scatter of this kind is called heteroskedasticity: the variance of the outcome around the regression line changes with the covariates. The short answer has two parts. A robust standard error changes the standard error and leaves the estimate alone. Whether it widens or narrows the interval depends on where the large variance sits, and this cohort shows both directions.

Constant and non-constant variance

A straight-line regression writes each person's blood pressure as the mean for their age plus an error:

Here $y_i$ is the SBP of person $i$, $x_i$ their age, $\beta_0$ the intercept, $\beta_1$ the slope and $e_i$ the error. A residual, $\hat e_i$, estimates it: the observed value minus the fitted one. Ordinary least squares (OLS), the usual fitting method, makes the sum of squared residuals as small as possible.

The usual OLS standard error assumes that $\mathrm{Var}(e_i \mid x_i)$, the variance of the error among people with the same age, equals one value, $\sigma^2$, for everyone. Constant variance is called homoskedasticity. When the variance changes with $x_i$, written $\sigma_i^2$, the errors are heteroskedastic.

In the simulated cohort, the random scatter added to each person's SBP has a standard deviation of $4 + 0.008\,(\text{age} - 30)^2$ mmHg. That is 4 mmHg at age 30, 9 at 55 and 24 at 80. Other simulated factors add further variation, so the residuals around a straight line on age are wider than this, most of all at younger ages. A residual plot, of residuals against fitted values or against age, still shows a fan that opens toward the older ages.

What the slope gets right

Heteroskedasticity does not, by itself, bias the OLS slope. The slope is unbiased whenever the mean model is right, meaning the errors average zero at every age, whatever their variances. Unequal variance costs precision instead. OLS gives a noisy 80-year-old the same weight as a quiet 30-year-old, so a method that weights by precision could estimate the slope more tightly.

In the simulated cohort, the fitted slope is 0.66 mmHg per year of age. Simulated blood pressure does not rise along an exact straight line with age, so the target here is the straight-line slope itself, the value a straight-line fit settles on in a very large sample. That simulated truth is 0.67, computed in the large simulated population the cohort was drawn from. The Stata output below shows this fit, called specification A (a specification is one way of setting up the regression), with its model-based standard error.

Stata: specification A with the model-based standard error

* Specification A: SBP on age (the high-leverage ages are the high-variance ages)
regress sbp age

. * Specification A: SBP on age (the high-leverage ages are the high-variance ages)
. regress sbp age

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(1, 1998)      =    951.44
       Model |  189678.648         1  189678.648   Prob > F        =    0.0000
    Residual |  398321.189     1,998  199.359955   R-squared       =    0.3226
-------------+----------------------------------   Adj R-squared   =    0.3222
       Total |  587999.837     1,999  294.146992   Root MSE        =    14.119

------------------------------------------------------------------------------
         sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .6644293   .0215406    30.85   0.000     .6221848    .7066737
       _cons |   88.26957   1.232598    71.61   0.000     85.85226    90.68689
------------------------------------------------------------------------------

What the model-based standard error assumes

The variance assumption enters through the standard error. Write $\bar x$ for the mean age and $S_{xx} = \sum_i (x_i - \bar x)^2$ for the spread of ages around it. A hat marks an estimate. The model-based variance of the slope is then

Here $s^2 = \sum_i \hat e_i^2 / (n - 2)$ pools the squared residuals of all $n$ people into one variance. The formula is exact only when everyone shares that variance. When the variance changes with age, $s^2$ is the wrong average, because the slope leans hardest on people far from the mean age. The standard error, the confidence interval and the P value can then be wrong in either direction.

Looking for unequal variance

Two plots usually settle the question. Plot the residuals against the fitted values and look for a fan or a funnel. Then plot them against each covariate, here age, or compare their standard deviation across age bands.

A formal test can back up what the plots show. The Breusch-Pagan test asks whether the covariates predict the squared residuals; the version run here does not assume normal errors. In the simulated cohort it gives a chi-square statistic of 184.24 on 1 degree of freedom.

Treat the test as supporting evidence only. In a large sample it can flag differences too small to matter, and in a small study it can miss ones that do. The paper that brought the robust standard error into wide use for linear regression also proposed a broader direct test [1].

Stata: the Breusch-Pagan test after the model-based fit

estat hettest age, iid

Assumption: i.i.d. error terms
Variable: age

H0: Constant variance

    chi2(1) = 184.24
Prob > chi2 = 0.0000

The sandwich estimator

The robust standard error, also called the sandwich or heteroskedasticity-consistent standard error, replaces the one pooled variance with each person's own squared residual [1, 2]. In matrix form, let $X$ be the design matrix, with one row per person, a column of 1s for the intercept and one column per covariate. Let $X'$ be its transpose, with rows and columns swapped, and $(X'X)^{-1}$ the matrix inverse of $X'X$, the matrix version of dividing by it.

Here $\hat\beta$ holds all the estimated coefficients, and $\mathrm{diag}(\hat e_i^2)$ is a diagonal matrix of the squared residuals. The two outer factors, $(X'X)^{-1}$, are the bread, and the middle factor is the meat, which gives the sandwich estimator its name. If every squared residual is replaced by the pooled $s^2$, the meat becomes $s^2 X'X$ and the formula collapses to the model-based $s^2 (X'X)^{-1}$. For a single slope it reduces to

Each squared residual counts in proportion to $(x_i - \bar x)^2$, the squared distance of that person's age from the mean age. The worked example below shows where this comes from.

HC0, HC1, HC2 and HC3

The formula above is HC0, the original version. Squared residuals tend to understate the variances they stand in for, most of all for high-leverage people, so later versions inflate them [3]. The leverage $h_{ii}$ of person $i$ measures how far their covariate values sit from everyone else's. With one covariate it is $1/n + (x_i - \bar x)^2 / S_{xx}$, and a high-leverage person pulls the fitted line toward themselves.

• HC1 multiplies HC0 by $n/(n - k)$, where $k$ is the number of coefficients; Stata's vce(robust) uses it.
• HC2 divides each $\hat e_i^2$ by $1 - h_{ii}$.
• HC3 divides each $\hat e_i^2$ by $(1 - h_{ii})^2$, inflating high-leverage residuals most; Stata's vce(hc3) uses it.

In large samples the four versions agree. In small samples HC0 and HC1 tend to be too small, and a simulation study recommends HC3 when the sample is small [4]. In R, the sandwich package computes all four, and its vcovHC function uses HC3 unless told otherwise [5]. To reproduce Stata's vce(robust) in R, ask for type = "HC1".

Which way the standard error moves

The direction of the change is not fixed. The robust standard error is larger when the large residual variance sits at high-leverage values of $x$, and can be smaller otherwise. The simulated cohort shows both, with the same 2,000 adults and the same blood pressures in two specifications, summarised in the table.

Specification A: blood pressure on age, robust larger

Ages run evenly from 30 to 80, so the high-leverage ages are the youngest and the oldest. The oldest are also the adults whose blood pressure scatters most. So the sandwich weights large squared residuals heavily, and the robust standard error rises from 0.0215 to 0.0243 mmHg per year. That is 1.13 times the model-based value, against 1.15 in the large simulated population.

The interval widens accordingly, while the slope itself does not move. HC3 gives the same 0.0243 to four decimals, because with 2,000 adults no single person has much leverage. R gives the same numbers.

Stata: specification A with HC1, the robust standard error of vce(robust)

regress sbp age, vce(robust)

. regress sbp age, vce(robust)

Linear regression                               Number of obs     =      2,000
                                                F(1, 1998)        =     750.63
                                                Prob > F          =     0.0000
                                                R-squared         =     0.3226
                                                Root MSE          =     14.119

------------------------------------------------------------------------------
             |               Robust
         sbp | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .6644293   .0242514    27.40   0.000     .6168686    .7119899
       _cons |   88.26957   1.203478    73.35   0.000     85.90937    90.62978
------------------------------------------------------------------------------

Stata: specification A with HC3

regress sbp age, vce(hc3)

. regress sbp age, vce(hc3)

Linear regression                               Number of obs     =      2,000
                                                F(1, 1998)        =     749.16
                                                Prob > F          =     0.0000
                                                R-squared         =     0.3226
                                                Root MSE          =     14.119

------------------------------------------------------------------------------
             |             Robust HC3
         sbp | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .6644293   .0242751    27.37   0.000     .6168221    .7120364
       _cons |   88.26957   1.204652    73.27   0.000     85.90707    90.63208
------------------------------------------------------------------------------

R: specification A with HC3 from the sandwich package

print(lmtest::coeftest(fa, vcov. = sandwich::vcovHC(fa, type = "HC3")))

> print(lmtest::coeftest(fa, vcov. = sandwich::vcovHC(fa,
+     type = "HC3")))

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)
(Intercept) 88.269574   1.204652  73.274 < 2.2e-16 ***
age          0.664429   0.024275  27.371 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Specification B: an under-40 indicator, robust smaller

Specification B replaces age with an indicator that is 1 for adults under 40 and 0 otherwise; 387 of the 2,000 adults have it. This dichotomises age only to show the opposite direction of the change; a real analysis would keep age continuous. The coefficient is now a difference in mean SBP, under-40s minus the rest: -17.65 mmHg, against a simulated truth of -17.86.

Now the high-leverage people are the 387 under-40s, the smaller group, and they are also the quiet one. The large residual variance sits among the older adults, the bigger, low-leverage group. The model-based standard error pools the two groups' variances and applies the result to both, which overstates the uncertainty of the small, quiet group and understates that of the large, noisy one. The first error wins, because the mean of the 387 under-40s carries most of the uncertainty in the difference.

The robust standard error falls from 0.8871 to 0.6558 mmHg, 0.74 times the model-based value, against 0.75 in the large simulated population. This time the robust interval is the narrower one. In R, fb <- lm(sbp ~ under40, data = d) gives the same standard errors to four decimals, as square roots of the variances from vcov(fb) for the model-based one and from sandwich::vcovHC(fb, type = "HC1"), which matches Stata's vce(robust), or type = "HC3" for the robust ones.

Stata: specification B with the model-based standard error

* Specification B: SBP under 40 versus 40 and over (the smaller group is the low-variance group)
regress sbp under40

. * Specification B: SBP under 40 versus 40 and over (the smaller group is the low-variance group)
. regress sbp under40

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(1, 1998)      =    395.85
       Model |  97232.6389         1  97232.6389   Prob > F        =    0.0000
    Residual |  490767.198     1,998  245.629228   R-squared       =    0.1654
-------------+----------------------------------   Adj R-squared   =    0.1649
       Total |  587999.837     1,999  294.146992   Root MSE        =    15.673

------------------------------------------------------------------------------
         sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     under40 |  -17.65015     .88712   -19.90   0.000    -19.38993   -15.91037
       _cons |   128.4365    .390232   329.13   0.000     127.6711    129.2018
------------------------------------------------------------------------------

Stata: specification B with HC1

regress sbp under40, vce(robust)

. regress sbp under40, vce(robust)

Linear regression                               Number of obs     =      2,000
                                                F(1, 1998)        =     724.34
                                                Prob > F          =     0.0000
                                                R-squared         =     0.1654
                                                Root MSE          =     15.673

------------------------------------------------------------------------------
             |               Robust
         sbp | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     under40 |  -17.65015   .6558071   -26.91   0.000    -18.93629   -16.36401
       _cons |   128.4365   .4172295   307.83   0.000     127.6182    129.2547
------------------------------------------------------------------------------

Stata: specification B with HC3

regress sbp under40, vce(hc3)

. regress sbp under40, vce(hc3)

Linear regression                               Number of obs     =      2,000
                                                F(1, 1998)        =     722.47
                                                Prob > F          =     0.0000
                                                R-squared         =     0.1654
                                                Root MSE          =     15.673

------------------------------------------------------------------------------
             |             Robust HC3
         sbp | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     under40 |  -17.65015   .6566548   -26.88   0.000    -18.93795   -16.36235
       _cons |   128.4365   .4172796   307.79   0.000     127.6181    129.2548
------------------------------------------------------------------------------

Coverage: how often the intervals contain the truth

A standard error is judged by the intervals it builds. The coverage of a 95% confidence interval is the share of repeated studies whose interval contains the value being estimated. It should be close to 95%. To check it, the simulation drew 5,000 new studies of 200 adults each from the cohort's equations and fitted both specifications in Stata.

The Monte Carlo standard error measures how much a coverage figure would vary if the whole simulation were run again, because only 5,000 studies were drawn. A coverage $p$ estimated from $m$ studies has a Monte Carlo standard error of $\sqrt{p(1 - p)/m}$, which is 0.0031 at exactly 95%.

In specification A, model-based intervals contained the simulated truth 91.7% of the time, because they were too narrow; HC1 and HC3 reached 95.0% and 95.1%. In specification B, model-based intervals contained it 99.0% of the time, because they were too wide; HC1 and HC3 reached 94.7% and 95.1%. All four robust coverages came close to 95%, though with only 200 people per study a robust interval can still fall a little short of it [4].

Alternatives: model the variance instead

A robust standard error patches the uncertainty and leaves the estimate as it was. When the variance pattern is understood, modelling it can give a more precise estimate as well. Weighted least squares (WLS) weights each person by the inverse of their variance, so quiet younger adults count more and noisy older ones less. Generalized least squares (GLS) extends this to a full variance and correlation structure; WLS is its special case without correlation.

A transformation, such as the logarithm of the outcome, can steady a variance that grows with the mean. It also changes the scale of the effect, and so the question being answered. A generalized linear model can instead use a family whose variance function matches the data. The variance function, $V(\mu)$, ties the variance to the mean $\mu$: a gamma family, for example, suits costs whose spread grows with their mean.

The link and family part of this series shows how the family sets the standard error. A modelled variance can itself be wrong, so it may be paired with a robust standard error as insurance.

The same sandwich in other models

The sandwich is used wherever the variance part of a model is known to be approximate. Modified Poisson regression estimates a risk ratio for a binary outcome, and the sandwich repairs the standard error that its Poisson variance gets wrong. Generalized estimating equations (GEE) model average outcomes over repeated visits, and a sandwich summed within each patient protects the standard error if the assumed correlation between visits is wrong. The Andersen-Gill model, a regression for the rate of repeated hospital admissions, clusters the sandwich on patient because one patient's admissions are correlated.

What the sandwich cannot fix

A sandwich standard error repairs the standard error when the variance assumption is wrong. It does not repair a wrong mean model, and a biased coefficient keeps its bias. When blood pressure rises along a curve and the model forces a straight line, the robust standard error is honest only about that straight-line slope [6]. The simulated cohort is such a case: both of its standard errors describe the straight-line slope, not the curve.

A large gap between the robust and model-based standard errors is better read as a prompt to check the mean model than as a fix [7]. The sandwich also needs enough data. With few observations it tends to be too small, which is why HC3 is preferred in small samples [4]. The clustered form used in GEE and Andersen-Gill needs many clusters, not just many rows, and with few it can be too small even when the model is right.

What to do in your own analysis

• Plot the residuals against the fitted values and against each main covariate before reading any standard error.
• Check the mean model first, including the shape of each continuous covariate.
• If the variance may change across covariates, consider prespecifying a robust standard error, HC3 when the sample is small, and name the version used.
• When the model-based and robust standard errors differ noticeably, consider reporting both in a supplement and saying which way the interval moved.
• For clustered data, consider clustering on the unit that was sampled or allocated, such as the patient or the clinic, and check that there are enough clusters.

Related in the wiki: [[risk-regression-models-epidemiology]] [[repeated-measures-modeling-guide]]