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

Solution 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)

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

Comments

[Joris Meys]

this is something you better ask at crossvalidated.com, they'll be able to help you more. And it would be nice to have some reproducible code while you're at it, together with the expected outcome. This often clears a problem up quite faster.

[Richard Herron]

I don't know stata, but it looks like your stata regression is a pooled linear model of Y = a0 + a1*V1 + a2*ENTITY + epsilon with robust het se, which is what you're doing with lm, so the results match. In the plm model you're doing an FE regression Y = a0 + a1*V1 + ui + epsilon, where ui is the FE for each "individual", which by index you've specified to be ENTITY. So I think your stata and R results match in the first case because you're doing a pooled panel with entity as an ind var in both cases. But I don't know stata.

[Helix123]

plm's type="sss" does not replicate Stata's small sample adjustment exactly because Stata has an intercept in its FE model and thus, uses a different amount of coefficients in the adjustment. plm does not have this intercept and, thus, uses "amount of coefficients - 1" relative to Stata in the adjustment. For small samples, this gives slightly different results but will go unnoticed in large samples.

[landroni]

@Helix123 Indeed. I've discussed this with the authors of plm and overall it is not clear what is the theoretical justification for Stata to include the intercept in these computations. While an intercept can be computed and plm does it as well if requested (within_intercept(), possibly still in SVN), it seems to be somewhat of an artificial concept in the context of within models. For now it didn't seem wise to include a "numerically exact" replicate of Stata clustered SEs... This said, Gretl seems to replicate these computations exactly à la Stata.