# R | Robust standard errors in panel regression clustered at level != Group Fixed Effects

I want to run a regression on a panel data set in R, where robust standard errors are clustered at a level that is not equal to the level of fixed effects. That is, I have a firm-year panel and I want to inlcude Industry and Year Fixed Effects, but cluster the (robust) standard errors at the firm-level.

Until now, I only had regressions where the group fixed effects were also the level of clustering. Hence, I used the following plm based procedure in R:

library(sandwich)
library(lmtest)
table <- plm(Dependent ~ Independent,
data=panel_1, index=c("FIRM","YEAR"), model="within", effect = "twoways")
robust_se <- sqrt(diag(vcovHC(table, type = "HC1", cluster="group")))
coeftest(table, vcov = vcovHC(table, type = "HC1", cluster="group"))


where I used the robust_se variable for stargazer output and the coeftest for console output.

However, now I can't use the cluster="group" option, since the index is defined at the industry level:

table <- plm(Dependent ~ Independent,
data=panel_1, index=c("INDUSTRY","YEAR"), model="within", effect = "twoways")


Do you have ideas how to solve this in R?

It is possible to cluster standard errors at a level other than the "group" or "time" indices of the plm model (e.g. industry). You would need to use the vcovCR function in the clubSandwich R package. Here is some example code I used: vcovCR(mypanelmodel, PanelData\$Owner_ID,type="CR2").
for reference, the lfe package provides multiway standard error clustering, for example:
felm(y ~x1 + x2 | fixedeffect1 + fixedeffect2 | Ivvariable | clustervar1 + clustervar2, data = data)