I have a random effects model about stock returns for different companies across several quarters. To handle clustered standard errors, I would like to use vcovHC and cluster by both company and quarter. I found several answers that would help to cluster a regular linear model by both dimensions but random effects would not be present, so asking for your help to include that.
Would the code below work for that or what would be your approach?
clustered_se <- vcovHC(random_model, type = "HC1", cluster = "group", group = "time")
coeftest(random_model, clustered_se)