Calculating a Robust Standard Error by Hand

From the Sandwich Formula to a Number

A heteroskedasticity-robust standard error is not a separate formula detached from the regression data. It is obtained from a covariance matrix built out of the design matrix and the fitted residuals. For a small regression, each matrix multiplication can be carried out by hand: calculate the leverage-adjusted residual weights if required, form a weighted cross-product, apply the two outer factors, and take the square root of the relevant diagonal entry.

The goal here is calculation, not a new diagnostic or a test of a regression coefficient. The model has a linear conditional mean, and the observations are independent. The error terms may have different variances across observations. A robust covariance estimate addresses that unequal variance in estimating coefficient uncertainty; it does not change the fitted least-squares coefficients or make a misspecified mean model correct.

Definition: Let $X$ be the $n \times p$ design matrix, with one row $x_i^\mathsf{T}$ for each observation, and let $\hat e_i$ be the ordinary least-squares residual for observation $i$. For a specified HC estimator, define $\widehat\Omega$ as the diagonal matrix of residual-based weights. The estimated robust covariance matrix is $$ \widehat{\operatorname{Cov}}_{\mathrm{HC}}(\hat\beta) =(X^\mathsf{T}X)^{-1}X^\mathsf{T}\widehat\Omega X(X^\mathsf{T}X)^{-1}. $$ The robust standard error of coefficient $\hat\beta_j$ is the square root of the $j$th diagonal entry of this matrix.

For HC0, the diagonal weight is $\hat e_i^2$. HC1 multiplies every HC0 weight by $n/(n-p)$, where $p$ is the number of columns in $X$, including the intercept when present. HC2 uses $\hat e_i^2/(1-h_{ii})$, and HC3 uses $\hat e_i^2/(1-h_{ii})^2$. Here $h_{ii}$ is the $i$th diagonal entry of the hat matrix $H=X(X^\mathsf{T}X)^{-1}X^\mathsf{T}$. These are alternative estimators; a calculation must use the specific one requested.

The matrix dimensions are a useful first check. The inverse $(X^\mathsf{T}X)^{-1}$ is $p \times p$; the middle matrix $X^\mathsf{T}\widehat\Omega X$ is also $p \times p$. Therefore the covariance estimate is $p \times p$, with one diagonal entry for each fitted coefficient. A negative diagonal entry indicates an arithmetic or setup error: the middle matrix is positive semidefinite when the weights are nonnegative.

Expanding the Middle of the Sandwich

Instead of multiplying several matrices all at once, it is often easier to calculate the middle matrix as a sum. Let $a_i$ be the specified HC weight for observation $i$, and let $x_i^\mathsf{T}$ be its row in the design matrix. Since $\widehat\Omega$ is diagonal,

X^\mathsf{T}\widehat\Omega X =\sum_{i=1}^{n}a_i x_i x_i^\mathsf{T}.

This identity follows by writing the product row by row: the $i$th diagonal entry of $\widehat\Omega$ contributes $a_i x_i x_i^\mathsf{T}$, and adding those contributions gives the full product. For an intercept and one predictor, $x_i^\mathsf{T}=(1,x_i)$, so each contribution is

a_i \begin{pmatrix} 1 & x_i\\ x_i & x_i^2 \end{pmatrix}.

This expansion is the main hand-calculation shortcut. It makes clear that residual size determines the weight, while the design row determines how that weight contributes to the variances and covariance of the coefficients.

Worked Example: A Complete HC2 Covariance Matrix

Suppose a hypothetical four-observation regression has an intercept and one predictor. The predictor values are $-1,0,1,2$, and the supplied least-squares residuals are $1,-1,-1,1$, respectively. Calculate the HC2 covariance matrix and the robust standard error of the slope. The residuals are compatible with least squares because their sum is zero and their predictor-weighted sum is $(-1)(1)+(0)(-1)+(1)(-1)+(2)(1)=0$.

The design matrix and its cross-product are

X= \begin{pmatrix} 1&-1\\ 1&0\\ 1&1\\ 1&2 \end{pmatrix}, \qquad X^\mathsf{T}X= \begin{pmatrix} 4&2\\ 2&6 \end{pmatrix}, \qquad (X^\mathsf{T}X)^{-1} = \begin{pmatrix} 0.3&-0.1\\ -0.1&0.2 \end{pmatrix}.

For example, the determinant of $X^\mathsf{T}X$ is $4(6)-2(2)=20$, giving the displayed inverse. The leverages are $h_{ii}=x_i^\mathsf{T}(X^\mathsf{T}X)^{-1}x_i$, where here $x_i^\mathsf{T}=(1,x_i)$. Substitution gives leverages $0.7,0.3,0.3,0.7$. Thus the HC2 weights, $\hat e_i^2/(1-h_{ii})$, are

a_1=\frac{1}{0.3}=\frac{10}{3},\quad a_2=\frac{1}{0.7}=\frac{10}{7},\quad a_3=\frac{1}{0.7}=\frac{10}{7},\quad a_4=\frac{1}{0.3}=\frac{10}{3}.

Using the row-product expansion, the middle matrix is

X^\mathsf{T}\widehat\Omega X = \begin{pmatrix} 200/21&100/21\\ 100/21&380/21 \end{pmatrix}.

For instance, its upper-right entry is the weighted sum of the predictor values: $(-1)(10/3)+(0)(10/7)+(1)(10/7)+(2)(10/3)=100/21$. Multiplying by the inverse on both sides gives

\widehat{\operatorname{Cov}}_{\mathrm{HC2}}(\hat\beta) = \begin{pmatrix} 0.3&-0.1\\ -0.1&0.2 \end{pmatrix} \begin{pmatrix} 200/21&100/21\\ 100/21&380/21 \end{pmatrix} \begin{pmatrix} 0.3&-0.1\\ -0.1&0.2 \end{pmatrix} = \begin{pmatrix} 0.7524&-0.3143\\ -0.3143&0.6286 \end{pmatrix}.

The second diagonal entry is the estimated variance of the slope. Therefore its HC2 standard error is $\sqrt{0.6286}=0.7928$ predictor-response units per predictor unit. The first diagonal entry similarly gives the intercept’s estimated variance, in squared response units. The off-diagonal entries estimate the covariance between the coefficient estimates.

Shortcuts for a Single Coefficient

A full covariance matrix is useful when uncertainty for several coefficients is needed, but a single standard error can sometimes be calculated without writing every entry. If $A=(X^\mathsf{T}X)^{-1}$, the variance estimate for coefficient $j$ is

\widehat{\operatorname{Var}}_{\mathrm{HC}}(\hat\beta_j) =\sum_{i=1}^{n}a_i\left(e_j^\mathsf{T}A x_i\right)^2,

where $e_j$ is the vector with a 1 in position $j$ and zeros elsewhere. To justify the expression, take the $j$th diagonal element of $A(\sum_i a_i x_i x_i^\mathsf{T})A$. Because $A$ is symmetric, the $i$th resulting term is $a_i(e_j^\mathsf{T}A x_i)(x_i^\mathsf{T}A e_j)=a_i(e_j^\mathsf{T}A x_i)^2$. Each term is nonnegative.

For a simple regression with an intercept and predictor values centered so that $\sum_i x_i=0$, the inverse cross-product matrix is diagonal. The slope-variance calculation then reduces to a weighted sum of $x_i^2$, divided by $(\sum_i x_i^2)^2$.

Worked Example: HC0 Slope Standard Error in a Centered Design

In a hypothetical three-observation regression, let the predictor values be $-1,0,1$ and the residuals be $1,-2,1$. These residuals meet the least-squares conditions: they sum to zero and their predictor-weighted sum is $(-1)(1)+(0)(-2)+(1)(1)=0$. Calculate the HC0 slope standard error.

The design cross-product is diagonal: \[ X^\mathsf{T}X= \begin{pmatrix} 3&0\\ 0&2 \end{pmatrix}, \qquad (X^\mathsf{T}X)^{-1}= \begin{pmatrix} 1/3&0\\ 0&1/2 \end{pmatrix}. \] For HC0 the weights are the squared residuals, $1,4,1$. The slope variance is the weighted predictor-square sum divided by $(\sum_i x_i^2)^2$:

\widehat{\operatorname{Var}}_{\mathrm{HC0}}(\hat\beta_1) =\frac{\sum_i \hat e_i^2x_i^2}{(\sum_i x_i^2)^2} =\frac{(1)(1)+(4)(0)+(1)(1)}{(1+0+1)^2} =\frac{2}{4}=0.5.

The resulting HC0 slope standard error is $\sqrt{0.5}=0.7071$. The middle observation contributes no slope variance in this centered design because its predictor value is zero, even though its residual is relatively large.

Leverage Adjustments Change the Weights

The same data and design can produce different robust standard errors under different HC choices. For the three-observation example, the leverage is $h_{ii}=1/3+x_i^2/2$: it is $5/6$ at $x_i=-1$ and $x_i=1$, and $1/3$ at $x_i=0$. HC3 divides each squared residual by $(1-h_{ii})^2$. The endpoint weights are therefore $1/(1/6)^2=36$, and the middle weight is $4/(2/3)^2=9$. Since only the endpoints have nonzero $x_i^2$, the HC3 slope variance is

\widehat{\operatorname{Var}}_{\mathrm{HC3}}(\hat\beta_1) =\frac{36(1)+9(0)+36(1)}{(1+0+1)^2} =\frac{72}{4}=18, \qquad SE_{\mathrm{HC3}}(\hat\beta_1)=\sqrt{18}=4.243.

This is much larger than the HC0 result for these deliberately small hypothetical data. The endpoint residuals are being enlarged substantially because those observations have high leverage. This calculation illustrates the effect of the correction; it does not imply that HC3 must always exceed HC0 for every coefficient and data set.

Worked Example: Applying the HC1 Multiplier

Use the same three-observation design and residuals to calculate the HC1 slope standard error. There are $n=3$ observations and $p=2$ coefficients, an intercept and a slope. HC1 multiplies every HC0 weight by $n/(n-p)=3$.

Because this common multiplier also multiplies the entire covariance matrix, the HC1 slope variance is three times the HC0 slope variance just calculated:

\widehat{\operatorname{Var}}_{\mathrm{HC1}}(\hat\beta_1) =\frac{3}{3-2}(0.5)=1.5, \qquad SE_{\mathrm{HC1}}(\hat\beta_1)=\sqrt{1.5}=1.225.

The multiplier applies to the covariance estimate, not directly to the standard error. Equivalently, the HC1 standard error is the HC0 standard error multiplied by $\sqrt{3}$, not by $3$.

A Reliable Hand-Calculation Sequence

Set up the design.
Write the rows of $X$, identify the number of coefficients $p$, and calculate $(X^\mathsf{T}X)^{-1}$.

Choose the requested HC weights.
Start from the supplied residuals and, for HC2 or HC3, calculate the required leverages before adjusting the squared residuals.

Form the middle matrix.
Add the weighted outer products $a_i x_i x_i^\mathsf{T}$, checking the entries and units as you go.

Complete the sandwich and extract the standard error.
Multiply by $(X^\mathsf{T}X)^{-1}$ on both sides, select the coefficient’s diagonal entry, and take its nonnegative square root.

Common Mistakes and Interpretation

Using the wrong HC adjustment. HC0 uses squared residuals; HC2 and HC3 additionally use leverage. Check the specified estimator before forming any weights.
Taking the square root too early. First calculate the covariance matrix or the selected variance. The standard error is the square root of a coefficient variance, not the square root of each residual weight.
Multiplying HC1 standard errors by its full correction factor. HC1 scales the variance by $n/(n-p)$, so it scales the standard error by the square root of that factor.
Ignoring units or the model’s limits. A slope standard error has the slope’s units, while an intercept standard error has response units. A robust standard error changes uncertainty estimation; it does not correct an incorrect conditional-mean form or establish causation.

In practice, calculations are usually performed by statistical software, but the hand procedure makes the result auditable. It also helps diagnose discrepancies: different results may arise from different HC choices, different design columns, or a leverage calculation error. Robust standard errors quantify estimated sampling variability under the fitted mean model and the chosen covariance estimator; they are not exact guarantees for a small sample.

Check Your Understanding

Use the formulas and calculation sequence above to answer these questions.

For a design matrix with an intercept and two predictors, what are the dimensions of the robust covariance matrix?
With predictor value $x_i=2$, residual $\hat e_i=-3$, and leverage $h_{ii}=0.5$, what is this observation’s HC2 weight?
Why does an HC1 multiplier of $n/(n-p)$ multiply the variance by that factor but the standard error by its square root?
In the centered three-observation example, why does the middle observation contribute zero to the HC0 slope variance?
Does a robust standard error change the least-squares coefficient estimate? Explain briefly.