About sim_slopes function to conduct clustered robust

Hello, help! I'm wondering how to use the sim_slopes function to conduct **clustered robust** plotting based on the `gvkey `of enterprises. I've tried using `cluster`, `v.cov`, and `v.cov.args`, but none of them worked:

_stata code: reg y c.x##c.mod i.fyear i.industry, cluster(gvkey)_

```
fit1 <- lm(y ~ x * mod + factor(industry) + factor(fyear), data = data)
cov_matrix <- vcovCL(fit1, cluster = data$gvkey)
cov_matrix2 = coeftest(fit1, vcov = cov_matrix)

sim_slopes(fit1, pred = x, modx = mod, jnplot = TRUE, cluster = "gvkey", data = merged_data)
# There was no error message, but the plotting result was the same as when not adding cluster = "gvkey". I guess it didn't take effect.

sim_slopes(fit1, pred = x, modx = mod, jnplot = TRUE, v.cov = sandwich::vcovCL(fit1, cluster = data$gvkey), v.cov.args = fit1)
# Error information: vcovargs[[which(sapply(vcovargs, function(x) length(x[[1]]) == : attempt to select less than one element in OneIndex

sim_slopes(fit1,
pred = x,
modx = mod,
jnplot = TRUE,
v.cov = sandwich::vcovCL(fit1, cluster = data$gvkey), # Specify the function for dynamically calculating the covariance
v.cov.args = fit1) # Specify the parameters for vcovCL
# Error information: vcovargs[[which(sapply(vcovargs, function(x) length(x[[1]]) == : attempt to select less than one element in OneIndex
```

Start by reproducing the reported calls to sim_slopes with fit1, cluster, v.cov, and v.cov.args using the clustered covariance examples in the issue. Compare the plotting result with and without cluster = "gvkey" and investigate the reported v.cov.args error. Done means the clustered robust calculation either works as documented or the limitation and correct usage are clearly established.

Written by the indexing model from the issue text.