How to compute two-way clustered standard errors when working with plm package in R?

Question

So far I've tried to compute two-way clustered standard errors for the well-known CRIME4 database popularized by Prof Wooldridge's Introductory Econometrics.

The point is, I'm working with plm package, but I have found no successful way of computing that particular correction. Until now, I have tried with the lmtest and sandwich packages which are my default for these situations, unfortunately it seems like it's not possible to cluster for both county and year, or at least I haven't found a way to do it.

I have nonetheless succeded to do it with the fixest package ( probably my favorite for regression stuff in R).

PD:Right before posting this question, I think I found a way to do it, but results are quite different.

Here's my code and my results.

    library(pacman)
    p_load(tidyverse,rio,AER,plm,estimatr,lmtest,sandwich,stargazer,fixest,lfe,modelsummary,plm,multiwayvcov,wooldridge)

    crime4 = wooldridge::crime4

    cat("using fixest")

    model_fd = feols(data = crime4,
                     clcrmrte ~ d83 +d84 +d85 +d86+ d87 + clprbarr + clprbcon + clprbpri + clavgsen + clpolpc,
                     cluster = c("county"))

    etable(model_fd) # First Differences Model

model_fe = feols(data=crime4,
                 lcrmrte ~ lprbarr +lprbconv + lprbpris +lavgsen + lpolpc|county + year,
                 cluster = c("county","year"))

etable(model_fe) # Fixed Effects model
#====================================================#

cat("using plm")

model_plm = plm(data = crime4,
                formula =lcrmrte ~ lprbarr +lprbconv + lprbpris +lavgsen + lpolpc,
                model = "within",index = c("county","year"),effect = "twoways")

coeftest(model_plm,vcov=vcovDC(model_plm,type = "sss")) 
# This way I was able to get two-way clustered standard errors in plm package, I guess...