Video Tutorial on Robust Standard Errors
[This article was first published on Coffee and Econometrics in the Morning, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Update: I have included a modified version of this summaryR() command as part of my package tonymisc, which extends mtable() to report robust standard errors. The tonymisc package is available on CRAN through the install.packages() command.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
If you have the right R commands at your disposal, it is simple to correct for heteroskedasticity using the robust correction that is commonly-used among economists. I recorded a video tutorial to describe the simplest (and most flexible) way I know to get R to compute robust standard errors.
The key is to use a “summary-style” command that has an option to correct for heteroskedasticity. The command I like to use is called summaryR(). Here is the script file with the summaryR() command.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## ---------------------------------------------------------------------------------------- ## | |
## Author: John Fox ## | |
## Source: http://r.789695.n4.nabble.com/R-extend-summary-lm-for-hccm-td815004.html ## | |
## Adapted by Tony Cookson. ## | |
## -- Only Change Made: Changed the name of the function (unwisely maybe) ## | |
## to summaryR from summaryHCCM.lm. I also changed the spelling of consistent ## | |
## ---------------------------------------------------------------------------------------- ## | |
summaryR.lm <- function(model, type=c("hc3", "hc0", "hc1", "hc2", "hc4"), ...){ | |
if (!require(car)) stop("Required car package is missing.") | |
type <- match.arg(type) | |
V <- hccm(model, type=type) | |
sumry <- summary(model) | |
table <- coef(sumry) | |
table[,2] <- sqrt(diag(V)) | |
table[,3] <- table[,1]/table[,2] | |
table[,4] <- 2*pt(abs(table[,3]), df.residual(model), lower.tail=FALSE) | |
sumry$coefficients <- table | |
p <- nrow(table) | |
hyp <- cbind(0, diag(p - 1)) | |
sumry$fstatistic[1] <- linearHypothesis(model, hyp,white.adjust=type)[2,"F"] | |
print(sumry) | |
cat("Note: Heteroscedasticity-consistent standard errors using adjustment", type, "\n") | |
} |
I found this function on an R-help discussion board where several people were answering someone’s question about extending the summary.lm() command.
I deserve none of the credit for writing this (credit goes to John Fox), but I consider it my duty to point out how nice this function is. I demonstrate how to use the function in this video
Here are is the script file I used in the video:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## -------------------------------------- ## | |
## R Tutorial on Robust Standard Errors ## | |
## Author: Tony Cookson ## | |
## -------------------------------------- ## | |
## Read Data | |
library(foreign) | |
ed = read.dta("C://R//edudat2.dta") | |
## Obtain regression object | |
educ.lm = lm(log(wage)~gender + educ + pareduc, data=ed) | |
## Summarize Output Assuming Homoskedasticity | |
summary(educ.lm) | |
## Need to load the car library | |
library(car) | |
## If you don't have car installed you'll need to install it | |
## Use this command: | |
## install.packages("car", dependencies = TRUE) | |
## Then, run the library command library(car) | |
## Robust Standard errors | |
summaryR(educ.lm) ## Special Function written by John Fox | |
## Can Specify Type of Heteroskedasticity correction | |
summaryR(educ.lm, type="hc0") ## White SE | |
summaryR(educ.lm, type="hc1") ## Stata's Default | |
summaryR(educ.lm, type="hc2") ## Unbiased under homoskedasticity | |
summaryR(educ.lm, type="hc3") ## Default (conservative) | |
## Test using linearHypothesis | |
linearHypothesis(educ.lm, "educ+pareduc =1", white.adjust="hc3") | |
linearHypothesis(educ.lm, "educ+pareduc =1", white.adjust="hc0") |
Here’s a link to the data set.
To leave a comment for the author, please follow the link and comment on their blog: Coffee and Econometrics in the Morning.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.