r - Clustered standard errors different in plm vs lfe -

when run cluster standard error panel specification plm , lfe results differ @ second significant figure. know why differ in calculation of se's?

set.seed(572015) library(lfe) library(plm) library(lmtest) # clustering example x <- c(sapply(sample(1:20), rep, times = 1000)) + rnorm(20*1000, sd = 1) y <- 5 + 10*x + rnorm(20*1000, sd = 10) + c(sapply(rnorm(20, sd = 10), rep, times = 1000)) facx <- factor(sapply(1:20, rep, times = 1000)) mydata <- data.frame(y=y,x=x,facx=facx, state=rep(1:1000, 20)) model <- plm(y ~ x, data = mydata, index = c("facx", "state"), effect = "individual", model = "within") plmtest <- coeftest(model,vcov=vcovhc(model,type = "hc1", cluster="group")) lfetest <- summary(felm(y ~ x | facx | 0 | facx)) data.frame(lfeclusterse=lfetest$coefficients[2],        plmclusterse=plmtest[2])  lfeclusterse plmclusterse 1   0.06746538   0.06572588 

the difference in degrees-of-freedom adjustment. usual first guess when looking differences in supposedly similar standard errors (see e.g., different robust standard errors of logit regression in stata , r). here, problem can illustrated when comparing results (1) plm+vcovhc, (2) felm, (3) lm+cluster.vcov (from package multiwayvcov).

first, refit models:

m1 <- plm(y ~ x, data = mydata, index = c("facx", "state"),   effect = "individual", model = "within") m2 <- felm(y ~ x | facx | 0 | facx, data = mydata) m3 <- lm(y ~ facx + x, data = mydata) 

all lead same coefficient estimates. m3 fixed effects explicitly reported while not m1 , m2. hence, m3 last coefficient extracted tail(..., 1).

all.equal(coef(m1), coef(m2)) ## [1] true all.equal(coef(m1), tail(coef(m3), 1)) ## [1] true 

the non-robust standard errors agree.

se <- function(object) tail(sqrt(diag(object)), 1) se(vcov(m1)) ##          x  ## 0.07002696  se(vcov(m2)) ##          x  ## 0.07002696  se(vcov(m3)) ##          x  ## 0.07002696  

and when comparing clustered standard errors can show felm uses degrees-of-freedom correction while plm not:

se(vcovhc(m1)) ##          x  ## 0.06572423  m2$cse ##          x  ## 0.06746538  se(cluster.vcov(m3, mydata$facx)) ##          x  ## 0.06746538  se(cluster.vcov(m3, mydata$facx, df_correction = false)) ##          x  ## 0.06572423