Recall that if heteroskedasticity is present in our data sample, the OLS estimator will still be unbiased and consistent, but it will **not** be efficient. Specifically, estimated standard errors will be biased, a problem we cannot solve with a larger sample size.

To correct for this bias, it may make sense to adjust your estimated standard errors. Two popular ways to tackle this are to use:

- “Robust” standard errors (a.k.a. White’s Standard Errors, Huber–White standard errors, Eicker–White or Eicker–Huber–White)
- Clustered Standard Errors

In practice, heteroskedasticity-robust and clustered standard errors are usually larger than standard errors from regular OLS — however, this is not always the case. For further detail on when robust standard errors are smaller than OLS standard errors, see Jorn-Steffen Pische’s response on Mostly Harmless Econometrics’ Q&A blog.

## “Robust” standard errors

The following example will use the ` CRIME3.dta `

Because one of this blog’s main goals is to translate STATA results in R, first we will look at the ` robust `

command in STATA. For backup on the calculation of heteroskedasticity-robust standard errors, see the following link: http://www.stata.com/support/faqs/stat/cluster.html. The formulation is as follows:

where number of observations, and the number of regressors (including the intercept). This returns a Variance-covariance (VCV) matrix where the diagonal elements are the estimated heteroskedasticity-robust coefficient variances — the ones of interest. Estimated coefficient standard errors are the square root of these diagonal elements.

#### STATA:

reg cmrdrte cexec cunem if year==93, robust

#### R:

The following bit of code was written by Dr. Ott Toomet (mentioned in the Dataninja blog). I added a degrees of freedom adjustment so that the results mirror STATA’s ` robust `

command results.

## Heteroskedasticity-robust standard error calculation. summaryw <- function(model) { s <- summary(model) X <- model.matrix(model) u2 <- residuals(model)^2 XDX <- 0 ## Here one needs to calculate X'DX. But due to the fact that ## D is huge (NxN), it is better to do it with a cycle. for(i in 1:nrow(X)) { XDX <- XDX + u2[i]*X[i,]%*%t(X[i,]) } # inverse(X'X) XX1 <- solve(t(X)%*%X) # Variance calculation (Bread x meat x Bread) varcovar <- XX1 %*% XDX %*% XX1 # degrees of freedom adjustment dfc <- sqrt(nrow(X))/sqrt(nrow(X)-ncol(X)) # Standard errors of the coefficient estimates are the # square roots of the diagonal elements stdh <- dfc*sqrt(diag(varcovar)) t <- model$coefficients/stdh p <- 2*pnorm(-abs(t)) results <- cbind(model$coefficients, stdh, t, p) dimnames(results) <- dimnames(s$coefficients) results }

To use the function written above, simply replace ` summary() `

with `summaryw() `

to look at your regression results — like this:

require(foreign) mrdr = read.dta(file="/Users/kevingoulding/data/MURDER.dta") # run regression reg4 = lm(cmrdrte ~ cexec + cunem, data = subset(mrdr,year == 93)) # see heteroskedasticity-robust standard errors summaryw(reg4)

These results should match the STATA output exactly.

## Heteroskedasticity-robust LM Test

It may also be important to calculate heteroskedasticity-robust restrictions on your model (e.g. an F-test).

## Clustered Standard Errors

Let’s say that you want to relax your homoskedasticity assumption, and account for the fact that there might be a bunch of covariance structures that vary by a certain characteristic – a “cluster” – but are homoskedastic within each cluster. Similar to heteroskedasticity-robust standard errors, you want to allow more flexibility in your variance-covariance (VCV) matrix. The result is clustered standard errors, a.k.a. cluster-robust.

#### STATA:

use wr-nevermar.dta reg nevermar impdum, cluster(state)

#### R:

In R, you first must run a function here called ` cl() `

written by Mahmood Ara in Stockholm University – the backup can be found here.

cl <- function(dat,fm, cluster){ attach(dat, warn.conflicts = F) library(sandwich) M <- length(unique(cluster)) N <- length(cluster) K <- fm$rank dfc <- (M/(M-1))*((N-1)/(N-K)) uj <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum)); vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N) coeftest(fm, vcovCL) }

After running the code above, you can run your regression with clustered standard errors as follows:

nmar = read.dta("http://www.montana.edu/econ/cstoddard/562/wr-nevermar.dta") # Run a plain linear regression regt = lm(nevermar ~ impdum, data = nmar) # apply the 'cl' function by choosing a variable to cluster on. # here, we are clustering on state. cl(nmar, regt, nmar$state)