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 )