Difference in clustering standard errors of coeftest and sum

Question

Suppose I have a panel data and would like to look at time fixed effects, i.e. effects constant across "state" but varying over time.

My understanding was always, that the estimated coefficient using either

1. de-meaned estimation
2. binary variable method

yield the same coefficient and standard errors. However, I get different standard errors while same coefficient values. Where does this difference in standard errors (p-value) come frome:

library(AER)

data("Fatalities")
Fatalities$fatal_rate <- Fatalities$fatal / Fatalities$pop * 10000


m1 <- plm(fatal_rate ~ beertax, 
          data = Fatalities,
          index = c("state","year"), 
          model = "within",
          effect=c("time"))
coeftest(m1,vcov = vcovHC, type = "HC1", cluster = "group")

m2 <- lm(fatal_rate ~ beertax + year -1, data = Fatalities)
summ(m2, robust="HC1", cluster="state", digits=4)

the results look like this

coeftest(m1,vcov = vcovHC, type = "HC1", cluster = "group")

t test of coefficients:

        Estimate Std. Error t value Pr(>|t|)   
beertax  0.36634    0.11904  3.0774 0.002264 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

while for sum

summ(m2, robust="HC1", cluster="state", digits=4)
MODEL INFO:
Observations: 336
Dependent Variable: fatal_rate
Type: OLS linear regression 

MODEL FIT:
F(8,328) = 588.7231, p = 0.0000
R² = 0.9349
Adj. R² = 0.9333 

Standard errors: Cluster-robust, type = HC1
---------------------------------------------------
                   Est.     S.E.    t val.        p
-------------- -------- -------- --------- --------
beertax          0.3663   0.1214    3.0176   0.0027
year1982         1.8948   0.1413   13.4080   0.0000
year1983         1.8128   0.1304   13.8995   0.0000
year1984         1.8231   0.1219   14.9580   0.0000
year1985         1.7843   0.1176   15.1747   0.0000
year1986         1.8787   0.1181   15.9042   0.0000
year1987         1.8793   0.1156   16.2614   0.0000
year1988         1.8938   0.1127   16.8073   0.0000
---------------------------------------------------

as we can see the coefficient for the variable of interest beertaxis the same, $0.3663$. However, the $t$-values are different: $3.0774$ vs $3.0176$

Answer 1

For cluster robust standard errors, we assume that the model of interest is

$$ y_g = X_g \beta + u_g $$

for all clusters $g = 1, ..., G$.

In both packages, (CRV1-) Cluster robust standard errors are computed as

$$ Var(\beta) = m \times (X'X)^{-1} \Omega(X'X)^{-1} $$

with $\Omega = \sum_{g=1}^{G} X_g u_g u_g' X_g'$.

The t-statistics above differ in how m is chosen. The 'standard way' is to set

$$ m = \frac{N-1}{N-k} \times \frac{G}{G-1} $$

where $N$ is the number of observations, $k$ is the number of estimated parameters and $G$ the number of clusters. For fixed effects that are swept out, the question obviously arises if $k$ should include fixed effects parameters or not.

You can investigate the effect of different choices of small sample corrections via the fixest package:

library(fixest)
m2 <- feols(
  fatal_rate ~ beertax |year, 
  data = Fatalities, 
  cluster = ~state
)
tstat(m2, ssc = ssc(adj = TRUE, 
                   fixef.K = "full"))["beertax"]
# beertax 
3.017638 
tstat(m2, ssc = ssc(adj = TRUE, 
                    fixef.K = "none"))["beertax"]
#  beertax 
3.049668 
tstat(m2, ssc = ssc(adj = TRUE, 
                    fixef.K = "nested"))["beertax"]
# beertax 
3.017638 

# include cluster.adj?
tstat(m2, ssc = ssc(adj = TRUE, 
                    fixef.K = "full", 
                    cluster.adj = FALSE))["beertax"]
# beertax 
#3.049572 
tstat(m2, ssc = ssc(adj = TRUE, 
                    fixef.K = "none", 
                    cluster.adj = FALSE))["beertax"]
# beertax 
#3.081941 
tstat(m2, ssc = ssc(adj = TRUE, 
                    fixef.K = "nested", 
                    cluster.adj = FALSE))["beertax"]
# beertax 
#3.049572 

I have managed to match the t-stats provided by summ() with adj = TRUE, fixef.K="full" and cluster.adj = TRUE (3.017). In consequence, summ by default includes the fixed effects when calculating $k$. I have, unfortunately, not managed to match the results produced by plm.