Pairwise comparisons for a regression with sandwich estimates (in R)

The question in short

I run a regression in R and made a boxplot of the response variable with grouping by one of the predictor variables. On this boxplot I'd like to add some information about the statistical model. What information (and how to display it (it is not an issue of programming)) would you suggest me to provide?

The developed question

I have several predictors: two categorical, non-ordinal predictors and one continuous predictor (below coded in R)

set.seed(81)
pred1 = rep(c('Car', 'Bike', 'Train', 'Airplane'), 6)
pred2 = rep(c('High', 'Low', 'Middle'), 8)
pred3 = rnorm(24)
resp = c(rnorm(12, sd = 1), rnorm(12, sd = 5))


resp is the response variable. I ran a regression with sandwich estimates:

require(sandwich)
require(lmtest)

m = aov(resp ~ pred1 + pred2)
coeftest(m, sandwich)

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.49642    0.73911 -0.6716  0.51034
pred1Bike    1.55917    1.16568  1.3376  0.19769
pred1Car     1.23873    1.24080  0.9983  0.33135
pred1Train   2.50882    0.91468  2.7428  0.01338 *
pred2Low     0.11613    1.00540  0.1155  0.90932
pred2Middle  0.51476    0.90924  0.5661  0.57829
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


And I boxplotted the groups for pred1:

require(ggplot2)
ggplot(data.frame(pred1, resp), aes(x=pred1, y=resp)) + geom_boxplot()


On this plot I'd like to add some letters in order to indicate the groups which are statistically similar (p.value < 0.05) as discussed here. Something like this:

ggplot(data.frame(pred1, resp), aes(x=pred1, y=resp)) + geom_boxplot() + annotate('text', x=1:4, y=6, label=c('a','b','a','b'), size = 8, color='red')


My question is:

How can I find these p.values for the pairwise comparisons with my robust regression? I can do what follows where m is a simple aov model:

TukeyHSD(m)


But the following doesn't work:

TukeyHSD(coeftest(m, sandwich))


I might missunderstand what these pairwise comparisons are, and what the results I currently have mean! Please let me know if you feel this! The aim of my question is for me to understand what is the best way to display the results of my statistical model on a boxplot.

Note: the variables pred2 and pred3 are used to withdraw some part of variance that I don't want to be accounted to the effect of pred1 (as pred1, pred2and pred3 are correlated in my case). Therefore, I guess it is better not to run simple pairwise t-test in order to get the p.values I'd like to add at the top of each boxplot.

• Please let me know if my question doesn't make much sense! – Remi.b Dec 9 '13 at 8:28
• This seems confused. You talk about robust regression, but apply analysis of variance. Asking for sandwich standard errors is not everyone's idea of "robust regression". You have box plots and want to enhance them. But boxplots show medians and selected quantiles and are at best indirectly related to analysis of variance. I know that some advise this combination, but I don't understand why. I'd plot raw data in a real problem with any plot showing detail in the data (e.g. stripplots). You do seem to lurch between asking a general question and soliciting R advice. – Nick Cox Dec 9 '13 at 11:03
• @NickCox Thanks for your comment! Using the function aov() or lm() does not change anything as I later read the results with coeftest(). My question is not an issue with R programming. I am not asking for code lines. I'd like to have advice concerning how I should highlight with a graph (boxplot or other) the differences found on for the variable pred1. I thought about a boxplot but didn't know how I could add significant values (not a programming issue!) given that I have n-1 p.values for n groups. – Remi.b Dec 9 '13 at 11:12
• You cite code that "doesn't work". It's best not to complicate your question with such details if they don't matter, and in any case not to post questions about code here as they are off-topic for CV. The one answer to date regards you as asking for code. Wanting both general ideas and specific code is understandable, but you are sending a mixed message. – Nick Cox Dec 9 '13 at 11:23

One solution is actually given as an example in the book on the multcomp package, section 4.6:

Bretz, F., Hothorn, T., & Westfall, P. H. (2011). Multiple comparisons using R. Boca Raton, FL: CRC Press.

One only needs to slightly adapt your code (everything needs to be in one data.frame instead of floating around):

require(multcomp)
require(sandwich)

set.seed(81)
pred3 = rnorm(24)
df <- data.frame(pred1 = rep(c('Car', 'Bike', 'Train', 'Airplane'), 6), pred2 = rep(c('High', 'Low', 'Middle'), 8)
, resp = c(rnorm(12, sd = 1), rnorm(12, sd = 5)))

m <- aov(resp ~ pred1 + pred2, df)

tukey <- glht(m, linfct = mcp(pred1 = "Tukey") , vcov = sandwich)

##          Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = resp ~ pred1 + pred2, data = df)
##
## Linear Hypotheses:
##                       Estimate Std. Error t value Pr(>|t|)
## Bike - Airplane == 0     1.559      1.166    1.34    0.547
## Car - Airplane == 0      1.239      1.241    1.00    0.748
## Train - Airplane == 0    2.509      0.915    2.74    0.058 .
## Car - Bike == 0         -0.320      1.422   -0.23    0.996
## Train - Bike == 0        0.950      1.149    0.83    0.838
## Train - Car == 0         1.270      1.225    1.04    0.726
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## (Adjusted p values reported -- single-step method)


Note that glht as default uses the single-step method as adjustment for alpha-error accumulation. If you want something else, you need to feed it to adjusted()

For example, using the Bonferroni-Holm correction (which I tend to use as I don't understand what single-step actually does):

summary(tukey, test = adjusted("holm"))


If you want no alpha error correction, which I do not recommend, this is also possible:

summary(tukey , test = adjusted("none"))