Why do I get different heteroscedasticity robust standard errors in R when using the plm package?

Question

I am now writing my bachelors thesis and I have come across some difficulties. I am about to do some panel regressions with time and entity fixed effects and I would therefore like to use the plm package. But when I do add fixed effects and want to have heteroscedasticity robust standard errors they seem to be incorrect.

Does anyone know why the HC standard errors differ?

Here is my code

# Load data
load(file="panel")
attach(panel)

# Load packages
library(lmtest)
library(plm)


# Create two models. The lm.model is a linear model and as the
# LAND variable is a factor variable representing countries
# (Land = Country in Swedish) this model will have entity fixed
# effects. In the plm.model the plm package is used and
# individual effects and within model is turned on (which is
# the same as entity fixed effects)
lm.model<-lm(NETTOSPARANDE ~ EURO + LAND, data=panel)
plm.model<-plm(NETTOSPARANDE ~ EURO, index=c("LAND","ÅR"), effect="individual", model="within", data=panel)

# When looking at the coefficents without heteroscadisity robust
# standard errors they are identical. They do also have the same
# value in stata.
coeftest(lm.model)[1:2,]
coeftest(plm.model)

# But when looking at the coefficents using heteroscadisity
# robust standard errors the lm.model and the plm.model produces
# different standard errors.
coeftest(lm.model, vcov.=vcovHC(lm.model, method="white2", type="HC1"))[1:2,]
coeftest(plm.model, vcov.=vcovHC(plm.model, method="white2", type="HC1"))

If you want to test the data it can be found here (the panel file) [1]: https://sourceforge.net/projects/emumoralhazard/files/ R-data

And here is my output

1> # Load data
1> load(file="panel")

1> attach(panel)

1> # Load packages
1> library(lmtest)
Loading required package: zoo

1> library(plm)
Loading required package: kinship
Loading required package: survival
Loading required package: splines
Loading required package: nlme
Loading required package: lattice
[1] "kinship is loaded"
Loading required package: Formula
Loading required package: MASS
Loading required package: sandwich

1> # Create two models. The lm.model is a linear model and as the
1> # LAND variabel is a factor variable representing countries
1> # (Land = Country in swedish) this model will have entity fixed
1> # effects. In the plm.model the plm package is used and
1> # individual effects and within model is turned on (which is
1> # the same as entity fixed effects)
1> lm.model<-lm(NETTOSPARANDE ~ EURO + LAND, data=panel)

1> plm.model<-plm(NETTOSPARANDE ~ EURO, index=c("LAND","ÅR"), effect="individual", model="within", data=panel)

1> # When looking at the coefficients without heteroscedasticity robust
1> # standard errors they are identical. They do also have the same
1> # value in Stata.
1> coeftest(lm.model)[1:2,]
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -3.731024  0.7731778 -4.825570 1.726921e-06
EURO1        2.187170  0.4076720  5.365024 1.112984e-07

1> coeftest(plm.model)

t test of coefficients:

      Estimate Std. Error t value  Pr(>|t|)    
EURO1  2.18717    0.40767   5.365 1.113e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


1> # But when looking at the coefficients using heteroscedasticity 
1> # robust standard errors the lm.model and the plm.model produces
1> # different standard errors.
1> coeftest(lm.model, vcov.=vcovHC(lm.model, method="white2", type="HC1"))[1:2,]
             Estimate Std. Error    t value     Pr(>|t|)
(Intercept) -3.731024  0.3551280 -10.506138 5.102122e-24
EURO1        2.187170  0.3386029   6.459395 2.009894e-10

1> coeftest(plm.model, vcov.=vcovHC(plm.model, method="white2", type="HC1"))

t test of coefficients:

      Estimate Std. Error t value  Pr(>|t|)    
EURO1  2.18717    0.33849  6.4615 1.983e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Answer 1

Though numerically both of your model coefficients estimates coincide, you are actually fitting two different models: least square dummy variables and fixed effects. They differ on their assumptions, so the robust standard errors are calculated differently.

Function vcovHC is a wrapper function, different functions are actually used in case of lm model and plm. In the first case the function vcovHC comes from the sandwich package, in the second from the plm package.

The difference in results arises because vcovHC for lm object does not exploit panel data structure. Comparing

plm:::vcovHC.plm

and

vcovHC.default

reveals that for panel data there is additional argument cluster = c("group", "time").

Look into Wooldridge book for more explanation why least squares dummy variables and fixed effects regressions differ. It has whole chapter dedicated for that.

Answer 2

@Skolnick - After an admittedly quick look at your problem, it seems that plm() truncates and/or rounds results to less digits than lm(). I think this is the most likely cause of the differences in your heteroscedasticity-robust standard error estimates. Explore summary(plm.model) and summary(lm.model) for differences in residual quantiles.