Imagine that for some model we calculate standard errors clustered on one dimension, e.g. state:

```
using DataFrames, RDatasets, FixedEffectModels
df = dataset("plm", "Cigar")
df[:StateCategorical] = categorical(df[:State])
df[:YearCategorical] = categorical(df[:Year])
reg(df, @model(Sales ~ NDI, fe = StateCategorical + YearCategorical, weights = Pop, vcov = cluster(StateCategorical)))
```

Now I would like to calculate a second set of standard errors, e.g. `vcov = cluster( StateCategorical + YearCategorical)`

.

How can we calculate these without having to run the costly

`reg(df, @model(Sales ~ NDI, fe = StateCategorical + YearCategorical, weights = Pop, vcov = cluster( StateCategorical + YearCategorical)))`

?