Heteroscedasticity robust standard errors with the PLM package

Question

I am trying to learn R after using Stata and I must say that I love it. But now I am having some trouble. I am about to do some multiple regressions with Panel Data so I am using the plm package.

Now I want to have the same results with plm in R as when I use the lm function and Stata when I perform a heteroscedasticity robust and entity fixed regression.

Let's say that I have a panel dataset with the variables Y, ENTITY, TIME, V1.

I get the same standard errors in R with this code

lm.model<-lm(Y ~ V1 + factor(ENTITY), data=data)
coeftest(lm.model, vcov.=vcovHC(lm.model, type="HC1))

as when I perform this regression in Stata

xi: reg Y V1 i.ENTITY, robust

But when I perform this regression with the plm package I get other standard errors

plm.model<-plm(Y ~ V1 , index=C("ENTITY","YEAR"), model="within", effect="individual", data=data)
coeftest(plm.model, vcov.=vcovHC(plm.model, type="HC1))

• Have I missed setting some options?
• Does the plm model use some other kind of estimation and if so how?
• Can I in some way have the same standard errors with plm as in Stata with , robust

Answer 1

By default the plm package does not use the exact same small-sample correction for panel data as Stata. However in version 1.5 of plm (on CRAN) you have an option that will emulate what Stata is doing.

plm.model<-plm(Y ~ V1 , index=C("ENTITY","YEAR"), model="within", 
    effect="individual", data=data)
coeftest(plm.model, vcov.=function(x) vcovHC(x, type="sss"))

This should yield the same clustered by group standard-errors as in Stata (but as mentioned in the comments, without a reproducible example and what results you expect it's harder to answer the question).

For more discussion on this and some benchmarks of R and Stata robust SEs see Fama-MacBeth and Cluster-Robust (by Firm and Time) Standard Errors in R.

See also:

• Clustered standard errors in R using plm (with fixed effects)

Answer 2

Is it possible that your Stata code is different from what you are doing with plm?

plm's "within" option with "individual" effects means a model of the form:

yit = a + Xit*B + eit + ci

What plm does is to demean the coefficients so that ci drops from the equation.

yit_bar = Xit_bar*B + eit_bar

Such that the "bar" suffix means that each variable had its mean subtracted. The mean is calculated over time and that is why the effect is for the individual. You could also have a fixed time effect that would be common to all individuals in which case the effect would be through time as well (that is irrelevant in this case though).

I am not sure what the "xi" command does in STATA, but i think it expands an interaction right ? Then it seems to me that you are trying to use a dummy variable per ENTITY as was highlighted by @richardh.

For your Stata and plm codes to match you must be using the same model.

You have two options:(1) you xtset your data in stata and use the xtreg option with the fe modifier or (2) you use plm with the pooling option and one dummy per ENTITY.

Matching Stata to R:

xtset entity year
xtreg y v1, fe robust 

Matching plm to Stata:

plm(Y ~ V1 + as.factor(ENTITY) , index=C("ENTITY","YEAR"), model="pooling", effect="individual", data=data)

Then use vcovHC with one of the modifiers. Make sure to check this paper that has a nice review of all the mechanics behind the "HC" options and the way they affect the variance covariance matrix.

Hope this helps.