Does using robust standard errors change effect size? I ran some linear regressions in R using lm, with an interaction term (cat x cat) as the predictor of interest (and also incorporating a covariate). To calculate the effect sizes of the interaction as a Cohen's d, I followed Jake Westfall's method as explained at http://jakewestfall.org/blog/index.php/2015/05/27/follow-up-what-about-uris-2n-rule/.
Specifically, in R, I ran:
RSS <- c(crossprod(model$residuals))
MSE <- RSS / length(model$residuals)
RMSE <- sqrt(MSE)
[code omitted to store the four cell means]
d=((A1-B1) - (A2-B2))/(2*RMSE)

But then I determined my data had some heteroscedasticity issues and redid my models with robust standard errors.
Specifically, I used the sandwich package to adjust my t- and p-values for reporting with:
coeftest(model, vcov. = vcovHC(model, type='HC3'))

And I reported the robust SEs with emmeans:
modelEM <- emmeans(model, ~MayoFile*TreatmentMDX,
                  vcov = sandwich::vcovHC(model,type='HC3'))

All fine and dandy, but now I'm stuck with a question I ought to know the answer to: Does the effect size change too now I'm using robust SEs? My gut says "No, the effect size is the same, what's changed is whether or not it's significant" - but I'm not entirely confident in that conclusion.
Thanks for reading!

There are some subtleties here, depending on why robust standard errors were estimated and what is meant by the effect size.
Robust standard errors for regression coefficients do not change estimates of expected mean square in linear regression models. McNeish et al, Psychological Methods 22:114-140 (2017), in discussing cluster-robust standard errors (CR-SEs), say explicitly:

CR-SEs can output model $R^2$ and effect size measures that are identical to what would be obtained through OLS because quantities used in these calculations (sum of squares, expected mean squares) are unaffected by the statistical correction to the standard error estimates and the computational formulas are equivalent to a single-level model.

That's straightforward when the CR-SEs are used to account for correlation structures and there is no heteroscedasticity. You might think about the robust correction to coefficient standard errors as changing the effective number of degrees of freedom while using the same mean-square error. For an effect size like Cohen's d, it's the residual standard deviation that's important, not the standard errors of the regression coefficient estimates. Even if the latter are changed with CR-SEs, the former isn't.
If there is heteroscedasticity (i.e., different within-group standard deviations), as this question posits, then does it make sense to report a single effect size for an interaction based on some single standard deviation estimate? I suppose that you could extend formulas for pooled standard deviation estimates under heteroscedasticity, which derive from formulas for variances of sums. But is that what's most useful to your audience?
It might make more sense to choose a different "effect size" estimate for the interaction term, for example comparing the magnitude of the difference in predictions with the interaction versus what you might predict based on A and B alone without the interaction. To my mind, at least, that practical significance of how much the interaction improves the model is more informative than comparing the interaction to a measure of residual standard deviation.
