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

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
attach(panel)

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> attach(panel)

1> library(lmtest)

1> library(plm)

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


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.

• Thank you for the answear. But is the first model pooled? I have added the LAND term wich is the entity dummy variable (as it's a facor variable. And as I can se the parameter estimates (EURO) are identical? – Skolnick Jan 10 '11 at 22:27
• @Skolnick, I've missed that at first, so I changed my answer accordingly – mpiktas Jan 10 '11 at 22:31
• I have tried to get reach of that book but the library has a copy that's not available for the moment. But I can perhaps ask which model that has the most accurate standard errors when HC-robust, LSDV or within? – Skolnick Jan 19 '11 at 9:32

@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.