How can I get the entire covariance matrix from coeftest() of a panel model

I'm running panel regressions using the plm package like this:

library("plm")
Data <- data.frame(id = c(rep("a",50), rep("b", 50)),
               y = rnorm(100),
               x = c(rnorm(50), rnorm(50, sd = 5)),
               z = c(rnorm(50), rnorm(50, sd = 3)))

panelmodel <- plm(y ~ x + z, 
                  data = Data, 
                  effect = "individual", 
                  model = "within", 
                  index = "id")

I want to adjust my standard errors for clustering which I can do using the lmtest package like:

library("lmtest")
coeftest(panelmodel, 
     vcov = vcovHC(panelmodel, 
                   type = "sss"), # I need this exact type of standard errors
     cluster = "id")

But I need the cluster adjusted covariance matrix that was used to calculate the standard errors here. The non-adjusted covariance matrix is used above and it is this:

vcovHC(panelmodel, type = "sss")

But this is not yet adjusted. How can I get the covariance matrix that coeftest has calculated to return the adjusted standard errors? I need that for further calculations with different covariances.

The multiwayvcov package unfortunately does not work with plm objects so that one is out.

Solution

See documentation: ?plm::vcovHC.plm: vcovHC(panelmodel, type = "sss") is already what you are looking for, it is the same as vcovHC(panelmodel, type = "sss", cluster = "group").

The argument cluster you are passing to coeftest() (and not to vcovHC()) is not evaluated as it is an unknown argument of the function and the function it is passed to (for your example, at least). You can check by comparing the results of these three commands, they are all the same:

coeftest(panelmodel, vcov. = vcovHC(panelmodel, type = "sss"))
coeftest(panelmodel, vcov. = vcovHC(panelmodel, type = "sss"), cluster = "id")
coeftest(panelmodel, vcov. = vcovHC(panelmodel, type = "sss"), cluster = "nonsense")