Address unequal variance between groups before applying contrasts for a linear model? (r)

Question

My Goal: I have an ordinal factor variable (5 levels) to which I would like to apply contrasts to test for a linear trend. However, the factor groups have heterogeneity of variance.

What I've done: Upon recommendation, I applied the contrasts to the independent variable and used lmRob() from robust pckg to create a robust linear model that tests the trend.

# assign the codes for a linear contrast of 5 groups, save as object
contrast5 <- contr.poly(5)

# set contrast property of sf1 to contain the weights
contrasts(SCI$sf1) <- contrast5

# fit and save a robust model (exhaustive instead of subsampling)
robmod.sf1 <- lmRob(ICECAP_A ~ sf1, data = SCI, nrep = Exhaustive)

summary.lmRob(robmod.sf1)

My problem: I have since been reading that robust regression is more suited to address outliers, and not heterogeneity of variance. (bottom of https://stats.idre.ucla.edu/r/dae/robust-regression/_ ) This UCLA page (among others) suggests the sandwich package to get heteroskedastic-consistent (HC) standard errors (such as in https://thestatsgeek.com/2014/02/14/the-robust-sandwich-variance-estimator-for-linear-regression-using-r/ ).

But these examples use a series of functions/calls to generate output that gives you the HC std errors that could be used to calculate confidence intervals, t-values, p-values etc.

My thinking is that if I use vcovHC(), I could get the HC std errors, but the HC std errors would not have been 'applied' to/a property of the model, so I couldn't pass the model (with the new HC std errors) through a function to apply the contrasts that I ultimately want. I don't want to conflate two separate concepts, but surely if a function addresses/down-weights outliers, that should at least somewhat address unequal variances as well?

Can anyone confirm if my reasoning is sound (and thus remain with lmRob()? Or suggest how I could just correct my standard errors and still apply the contrasts?

UPDATE: Following Noah's advice, I've rewritten my code as:

# Fit a regular linear model
sf1_lm <- lm(ICECAP_A ~ sf1, data = SCI)
summary(sf1_lm)

# Create contrast matrix specifying coefs for a linear trend w 5 vars
cont5 <- matrix(c(-2, -1, 0, 1, 2), 1)

# Run model, specify contrast set & SEs robust to heteroscedasticity
sf1_hc  <- glht(sf1_lm, linfct = cont5, vcov. = vcovHC )
summary(sf1_hc)

This output produced results (significant) that were consistent with my expectation based on a visual depiction of group means as well as both the regular lm() and lmRob() results. Yay!

New follow-up question: I checked this method again on a variable that tested significant for unequal variances (Levene's), where the means of 2 of the 3 factor levels were very similar, and the CIs were (visually) not hugely different. Using lm() produced a (barely) significant linear trend (p=0.0485) while lmRob() (which I know now was suboptimal) was insignificant (p=0.402).

... BUT when I tested with the method under "UPDATE", the p-value for the linear trend became tiny (p<2e-16). I am not surprised by the p-value fluctuating a bit so that it becomes 'significant' again, but a change of this magnitude has me worried. Should I be? Could it be that the standard errors were 'over adjusted' by vcov. = vcovHC ?

Further concern: To double-check this that the code works properly, I re-levelled the groups so that the group with the lower mean is in the second position, while the 2 groups with the similar means are positioned 1st and 3rd, respectively. This should give a curved pattern, not a linear trend. BUT when I re-ran the analysis, it still gave me (p<2e-16) which it should not. This indicates to me that something is off.

After re-levelling:

# Fit a regular linear model
educ_lm <- lm(ICECAP_A ~ edu_relevel, data = SCI)

# Create contrast matrix specifying coefs for a linear trend w 3 vars
cont3 <- matrix(c( -1, 0, +1), 1)

# Run model, specify contrast set & SEs robust to heteroscedasticity
educ_hc  <- glht( educ_lm , linfct = cont3 , vcov. = vcovHC )
summary( educ_hc )

SOLUTION? I believe that the follow-up problem was due to an incomplete argument. Need to add linfct = mcp(), otherwise it might only compares the 1st (lowest) and 3rd (highest) groups, which had a difference in means. Now, this test produces a non-significant p-value (as expected). And, when I remove the vcov. argument, the p-value matches that produced through another method (applying contrasts to the variable then running through ANOVA).

educ_hc  <- glht( educ_lm , linfct = mcp(educ_col = cont3) , vcov. = vcovHC )
summary( educ_hc )

Answer 1

You are correct in that robust regression is not designed to address heteroscedasticity but rather is designed to reduce the impact of outliers. Note that these are two separate issues: outliers increase the variance of estimated coefficients, while heteroscedasticity invalidates the use of the common estimator of the standard error of the estimated coefficients. Based on your description, your primary concern with with the latter. Robust standard errors were designed specifically to solve this problem.

To perform contrasts, you can use the glht function of the multcomp package. In the documentation for glht, you'll see that the description for ... is

additional arguments to function modelparm in all glht methods.

If we look at modelparm, we see that one of its arguments is vcov.:

an accessor function for the covariance matrix of the model parameters. Alternatively, the covariance matrix directly.

vcovHC() from the sandwich package is "an accessor function for the covariance matrix of the model parameters". So simply supplying vcov. = vcovHC to glht() will produce standard errors that are robust to heteroscedasticity for the standard errors of your contrasts.