Clustered/Robust standard errors in R

Question

I am estimating an OLS regression without fixed effects and an OLS regression with fixed effects in R Studio. I have read, that it is common to use robust standard errors, when estimating a simple OLS model and to use clustered standard errors when estimating a fixed-effects model...however, I am not sure If I did it correctly! I thank you for your advice in advance!

my data is paneldatafinal, and I have a panel structure with IND and year.

This is the OLS model with robust standard error:

# MODEL EPS MARKET BASED:

model_mb_basic_ols <- lm(diff_log_VATFP_I ~ log_Impen_Mil + RuDintensity_vapc +
                           dummy_var*EPS_MKT_growth_MA3 + log_PRDK_growth_rate + GAP, 
                         data = paneldatafinal)


# Compute robust standard errors:
robust_se <- vcovHC(model_mb_basic, type = "HC1")

# Display the results:
summary_result <- coeftest(model_mb_basic, vcov = robust_se)

# Print summary result with robust standard errors:
print(summary_result)

This is the Fixed effect model:

## Model with clustered Standard Errors & fixed effects EPS TOTAL:

model_basicfe <- plm(diff_log_VATFP_I ~ log_Impen_Mil + RuDintensity_vapc + 
                       dummy_var*EPS_growth_MA3 + log_PRDK_growth_rate + GAP,
                    data = paneldatafinal, model = "within", index = c("IND", "year"))

# Calculate clustered standard errors:
vcov_clustered <- vcovHC(model_basicfe, type = "HC1", cluster = "group")

# Apply the clustered standard errors to the model:
coeftest(model_basicfe, vcov = vcov_clustered)

# Use stargazer to present the results:
stargazer(model_basicfe, type = "text", se = list(sqrt(diag(vcov_clustered))), header = FALSE)

Answer 1

Yes, the coding is correct. Note that there is no need for lmtest::coeftest here (it is a great package/function anyways!), but you can use plm's built-in summary method on a plm object. This has the advantage that you will also get adjusted F-test results. Further, for random effect models summary will automatically use z-tests for coefficients and Chi-sq test for the Wald test (whereas for coeftest you would set argument df = Inf manually).

You can also input vcovHC as a function with parameters. See/copied from ?plm::vcovHC:

library(plm)
data("Produc", package = "plm")
zz <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp,
          data = Produc, model = "within")
## as function input to plm's summary method (with and without additional arguments):
summary(zz, vcov = vcovHC)
#> Oneway (individual) effect Within Model
#> 
#> Note: Coefficient variance-covariance matrix supplied: vcovHC
#> 
#> Call:
#> plm(formula = log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, 
#>     data = Produc, model = "within")
#> 
#> Balanced Panel: n = 48, T = 17, N = 816
#> 
#> Residuals:
#>      Min.   1st Qu.    Median   3rd Qu.      Max. 
#> -0.120456 -0.023741 -0.002041  0.018144  0.174718 
#> 
#> Coefficients:
#>             Estimate Std. Error t-value  Pr(>|t|)    
#> log(pcap) -0.0261497  0.0603262 -0.4335   0.66480    
#> log(pc)    0.2920069  0.0617425  4.7294 2.681e-06 ***
#> log(emp)   0.7681595  0.0816652  9.4062 < 2.2e-16 ***
#> unemp     -0.0052977  0.0024958 -2.1226   0.03411 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Total Sum of Squares:    18.941
#> Residual Sum of Squares: 1.1112
#> R-Squared:      0.94134
#> Adj. R-Squared: 0.93742
#> F-statistic: 406.02 on 4 and 47 DF, p-value: < 2.22e-16
summary(zz, vcov = function(x) vcovHC(x, method="arellano", type="HC1"))
#> Oneway (individual) effect Within Model
#> 
#> Note: Coefficient variance-covariance matrix supplied: function(x) vcovHC(x, method = "arellano", type = "HC1")
#> 
#> Call:
#> plm(formula = log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, 
#>     data = Produc, model = "within")
#> 
#> Balanced Panel: n = 48, T = 17, N = 816
#> 
#> Residuals:
#>      Min.   1st Qu.    Median   3rd Qu.      Max. 
#> -0.120456 -0.023741 -0.002041  0.018144  0.174718 
#> 
#> Coefficients:
#>             Estimate Std. Error t-value  Pr(>|t|)    
#> log(pcap) -0.0261497  0.0604746 -0.4324   0.66557    
#> log(pc)    0.2920069  0.0618944  4.7178 2.834e-06 ***
#> log(emp)   0.7681595  0.0818661  9.3831 < 2.2e-16 ***
#> unemp     -0.0052977  0.0025020 -2.1174   0.03455 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Total Sum of Squares:    18.941
#> Residual Sum of Squares: 1.1112
#> R-Squared:      0.94134
#> Adj. R-Squared: 0.93742
#> F-statistic: 404.03 on 4 and 47 DF, p-value: < 2.22e-1