# Difference in clustering standard errors of coeftest and sum

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

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$$

• Hi, you are right that the demeaned regression and the binary/dummy specification should lead to exactly identical results, which is indeed the case for the point estimate. The difference in t-stats is likely to due how small sample corrections are treated between the two implementations. Have a look at the outstanding standard errors vignette of the fixest package for all the details. May 28, 2022 at 21:53
• @A.Fischer fair point and thanks for the reference to the vignette. I will have a look
– math
May 29, 2022 at 8:24
• Btw, could you please clarify which package the summ function is from? jtools? May 29, 2022 at 8:42

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

tstat(m2, ssc = ssc(adj = TRUE,
fixef.K = "full",
# beertax
#3.049572
tstat(m2, ssc = ssc(adj = TRUE,
fixef.K = "none",

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.