It's been quite a while. I hope you've found the answer to your question by now. However, I thought I'd post this in case anyone else faces a similar problem.

The core issue is that sandwich::vcovCL and plm::vcovHC use different finite-sample adjustments by default.

Denote the raw sandwich covariance matrix as $\widehat{\text{Var}}^{\text{raw}}$.

The default behavior of plm::vcovHC(..., type = "HC0", cluster = "group") is to return exactly this raw estimator with no additional scaling factor:

$$ \widehat{\text{Var}}^{\text{plmHC0}} = \widehat{\text{Var}}^{\text{raw}} $$

By default, sandwich::vcovCL applies the cadjust = TRUE argument to adjust for the number of clusters $G$:

$$ \widehat{\text{Var}}^{\text{sandwich}} = \frac{G}{G-1} \cdot \widehat{\text{Var}}^{\text{raw}} $$

Therefore, comparing plm::vcovHC(type="HC0") (raw) with sandwich::vcovCL(..., type="HC0", cadjust=TRUE) (adjusted) will produce different standard errors. The sandwich variances will be larger by a factor of $\frac{G}{G-1}$ (and thus the standard errors will be larger by a factor of $\sqrt{\frac{G}{G-1}}$).

To obtain the raw (unadjusted) standard errors with sandwich::vcovCL, you need to set cadjust = FALSE. If you prefer to use plm but want the default sandwich adjustment, you must manually multiply the raw covariance matrix by $\frac{G}{G-1}$.

Here is a complete reproducible example:

# Load required packages
library(plm)
library(sandwich)
library(lmtest)

# Load data
data("Grunfeld", package = "plm")

# --- Estimation ---

# OLS estimation with explicit dummies
est_lm <- lm(
  inv ~ value + capital + as.factor(firm) + as.factor(year),
  data = Grunfeld
)

# Two-way fixed effects plm estimation
est_plm <- plm(
  inv ~ value + capital,
  data = Grunfeld,
  model = "within",
  effect = "twoways"
)

# --- Comparison 1: Raw (unadjusted) standard errors ---
# These two must be identical
coeftest(est_lm, vcov = vcovCL(est_lm, cluster = ~ firm, type = "HC0", cadjust = FALSE))[2:3, ]
coeftest(est_plm, vcov = vcovHC(est_plm, type = "HC0", cluster = "group"))

# --- Comparison 2: Default sandwich adjustment (G/(G-1)) ---
# Default sandwich
coeftest(est_lm, vcov = vcovCL(est_lm, cluster = ~ firm, type = "HC0", cadjust = TRUE))[2:3, ]

# Manually replicate this in plm
G <- length(unique(Grunfeld$firm))
vcov_plm_corrected <- vcovHC(est_plm, type = "HC0", cluster = "group") * (G / (G - 1))
coeftest(est_plm, vcov = vcov_plm_corrected)