A function that reports White Standard Errors

Update: Since writing this post, I found a natural extension written by John Fox that does what the function in this post does and I have extended it in my CRAN-contributed library tonymisc. In my package, this modified and extended function is called summaryR().

One of the frustrating things about using R for econometrics is that R has some cumbersome functionality for computing the standard (in econometrics, anyway) White standard errors. As one of my favorite references on R learning points out, there are two ways to retrieve these precious standard errors.

1. You could install and load the sandwich library. Then, apply the function vcovHC to your lmHC0. That will give you the White covariance matrix. Next, take the square root of the diagonals, and you have retrieved your standard errors.
2. Install and load the car library. Then, apply the function hccm to your lm object, specifying the option hc0. Again, this only gives the covariance matrix, which you have to transform to obtain standard errors.

Both options produce the exact same answers (I tried it myself), but both are clunky if you want to use the White standard errors as a routine check for whether it is important. To mitigate this issue (at least for myself), I wrote a function that does the work for you.

The function is called summarize.r(), and it works just like the summary command. In fact, my function calls summary.lm, and then modifies the contents that are printed using method 1 from above.

To run my function, (1) download and install the sandwich library, (2) run the code that defines the function (see below), and (3) after fitting your regression using lm(), you can use the function by calling summarize.r(yourregression.lm).

Without further ado, here's the code you should run.

summarize.r = function(reg.lm, ...){
z = summary.lm(reg.lm)
W.VCOV = vcovHC(reg.lm,type="HC0")
sig.sq = ((summary(reg.lm)$sigma)^2)
coef = z$coefficients[,1]


White.SE <- sqrt(diag(W.VCOV)) t.robust <- coef/White.SE ## Element-by-element division df <- reg.lm$df.residual p.val.robust <- 2*(1-pt(abs(t.robust),df)) z$cov.unscaled = W.VCOV/sig.sq z$coefficients[,2] = White.SE z$coefficients[,3] = t.robust z$coefficients[,4] = p.val.robust z }

And, here's an example of putting the code to use. First, note what the summary output gives (cutting the call, residuals, and significance codes). Next, see the output from my summarize.r command.

>  summary(mathcomp.lm)

Coefficients:
       Estimate Std. Error t value Pr(>|t|)
(Intercept)   643.59       2.06 312.388  <> summarize.r(mathcomp.lm)

Coefficients:
       Estimate Std. Error t value Pr(>|t|)
(Intercept)  643.587      2.053 313.495  <>