SLID1 <- read.table("c:/temp/SLID-Ontario.txt",sep = "\t",header=F,skip=1)
names(SLID1)  <- c("age","sex","wages","yearsed")

SLID1$male <- as.numeric(SLID1$sex == "Male")

 # ----------- fit the model on the original scale -----------------------------

SLID1.lm1 <- lm(wages ~ male +age +yearsed, data=SLID1)

summary(SLID1.lm1)

 # automatically generates several useful plots
par(mfrow=c(2,2)) 
plot(SLID1.lm1)
par(mfrow=c(1,1)) 

 # ------- now use HCC "Heterogeneity Consistent Covariance Matrix" estimators -------

library(AER)

 # --- the function vcovHC gives HCC estimators -------------------------------------
 # --- the argument "const" gives the usual LS covariance estimator -----------------
 # --- we extract the diagonal and take a square root to get the standard errors ----

sqrt(diag(vcovHC(SLID1.lm1, type = "const")))


 # --- the argument "HC0" gives the original HCC estimator from White (1980) --------

sqrt(diag(vcovHC(SLID1.lm1, type = "HC0")))


 # --- the argument "HC3" gives the HCC estimator recommended by Long and Ervin (2000) --------

sqrt(diag(vcovHC(SLID1.lm1, type = "HC3")))