Fixed effect regression with clustered standard errors

I'm currently working on a fixed effects regression in r.

I have a dataset consisting of the following variables: total compensation of the CEO of a firm (TOTAL_COMP), the firm code (GVKEY), the fiscal year (FISCAL_YEAR), a number of firm characteristics (like assets (AT)) and the industry of the firm (SIC).

Although the data is proprietary it looks like:

e.g.

TOTAL_COMP <- c(100, 200, 50, 150, 300, 200, 150, 75)

GVKEY <- c("103", "103", "103", "103", "104", "104", "104", "104")

AT <- c(1000, 1100, 1200, 1300, 2000, 2100, 2200, 2300)

FISCAL_YEAR <- c(2015, 2016, 2017, 2018, 2015, 2016, 2017, 2018)

SIC <- c(78, 78, 78, 78, 80, 80, 80, 80)

I used the plm function and ran a fixed effects regression with industry (SIC) fixed effects.

model <- plm(TOTAL_COMP ~ AT, index = "SIC", model = "within", data = COMBINED_DATA)

I subsequently want to cluster the standard errors at the firm (GVKEY) and not the industry (SIC) level. I am, however, not sure which function to use. I tried:

coeftest(model, vcov = vcovHC(regression, type = "HC0", cluster = "group"))

but am afraid that this clusters the standard errors at the industry (SIC) instead of the firm (GVKEY) level. How do I compute the standard errors at the firm (GVKEY) level?

(packages used: plm and lmtest)

Solution

You can set SIC as dummy variable and cluster the standard errors at the firm level:

# The fixed effects model
model <- lm(TOTAL_COMP ~ AT + factor(SIC), data = COMBINED_DATA)

# The fixed effects model with cluster settings
library(estimatr)
Clu_robust <- lm_robust(TOTAL_COMP ~ AT + factor(SIC), 
                        data = COMBINED_DATA,
                        clusters = GVKEY, 
                        se_type = 'stata')