Comparaciones de pares para una regresión con sándwich de estima (en R)

Question

La pregunta en corto

Puedo ejecutar una regresión robusta en R e hizo un boxplot con sólo uno de los predictores de la variable de agrupación de la variable de respuesta. En este boxplot me gustaría añadir algo de información sobre el modelo estadístico. Qué información (y cómo mostrar (no es un problema de programación)) me sugieren?

---

El desarrollado pregunta

Tengo varios predictores: dos no ordinal predictor y un continuo predictor (abajo codificado en 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 es la variable de respuesta. Me encontré con una robusta de regresión con sándwich de estimaciones:

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

Y yo boxploted los grupos para pred1:

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

En esta trama me gustaría añadir algunas letras para indicar los grupos que son estadísticamente similares (p.valor < 0.05) como se describe aquí. Algo como esto:

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')

Mi pregunta es:

¿Cómo puedo encontrar estos p.los valores de las comparaciones por pares con mi regresión robusta? Puedo hacer lo que sigue donde m es un simple aov modelo:

TukeyHSD(m)

Pero el siguiente no funciona:

TukeyHSD(coeftest(m, sandwich))

Yo podría missunderstand ¿cuáles son estas las comparaciones por pares y cuáles son los resultados que tengo actualmente significa! Por favor, hágamelo saber si usted se siente esto! El objetivo de mi pregunta es para mí entender lo que es la mejor manera de mostrar los resultados de mi modelo estadístico en un boxplot.

Nota: las variables pred2 y pred3 se utiliza para retirar parte de la varianza que no quiero ser tenido en cuenta para el efecto de la pred1 (como pred1, pred2y pred3 están correlacionadas en mi caso). Por lo tanto, supongo que es mejor no correr simple pares de t-test, con el fin de obtener el p.los valores que me gustaría añadir en la parte superior de cada boxplot.

Answer 1

Una solución es en realidad un ejemplo en el libro de la multcomp paquete, sección 4.6:

Bretz, F., Hothorn, T., & Westfall, P. H. (2011). Comparaciones múltiples mediante R. Boca Raton, FL: CRC Press.

Sólo se necesita ligeramente adaptar su código (todo lo que necesita para estar en uno de los datos.marco en lugar de flotar a su alrededor):

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)

summary(tukey, test = adjusted())

##          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)

Tenga en cuenta que glht como predeterminado utiliza el solo paso como el método de ajuste para la alfa-error de acumulación. Si usted desea algo más, usted necesita para alimentar a adjusted()

Por ejemplo, utilizando la corrección de Bonferroni-Holm corrección (que yo tiendo a usar como no entiendo lo de un solo paso, en realidad no):

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

Si no desea una alfa de corrección de error, que yo no recomiendo, esto también es posible:

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