¿Cómo obtener los resultados de una prueba post-hoc Tukey HSD en una tabla que muestra pares agrupados?

13

Me encantaría realizar una prueba post-hoc de TukeyHSD después de mi Anova bidireccional con R, obteniendo una tabla que contenga los pares ordenados agrupados por diferencia significativa. (Perdón por la redacción, todavía soy nuevo con las estadísticas).

Me gustaría tener algo como esto:

ingrese la descripción de la imagen aquí

Entonces, agrupados con estrellas o letras.

¿Alguna idea? Probé la función HSD.test()desde el agricolaepaquete, pero parece que no maneja tablas de dos vías.

Stragu
fuente

Respuestas:

22

La agricolae::HSD.testfunción hace exactamente eso, pero deberá hacerle saber que está interesado en un término de interacción . Aquí hay un ejemplo con un conjunto de datos Stata:

library(foreign)
yield <- read.dta("http://www.stata-press.com/data/r12/yield.dta")
tx <- with(yield, interaction(fertilizer, irrigation))
amod <- aov(yield ~ tx, data=yield)
library(agricolae)
HSD.test(amod, "tx", group=TRUE)

Esto da los resultados que se muestran a continuación:

Groups, Treatments and means
a        2.1     51.17547 
ab       4.1     50.7529 
abc      3.1     47.36229 
 bcd     1.1     45.81229 
  cd     5.1     44.55313 
   de    4.0     41.81757 
    ef   2.0     38.79482 
    ef   1.0     36.91257 
     f   3.0     36.34383 
     f   5.0     35.69507 

Coinciden con lo que obtendríamos con los siguientes comandos:

. webuse yield
. regress yield fertilizer##irrigation
. pwcompare fertilizer#irrigation, group mcompare(tukey)

-------------------------------------------------------
                      |                           Tukey
                      |     Margin   Std. Err.   Groups
----------------------+--------------------------------
fertilizer#irrigation |
                 1 0  |   36.91257   1.116571    AB    
                 1 1  |   45.81229   1.116571      CDE 
                 2 0  |   38.79482   1.116571    AB    
                 2 1  |   51.17547   1.116571         F
                 3 0  |   36.34383   1.116571    A     
                 3 1  |   47.36229   1.116571       DEF
                 4 0  |   41.81757   1.116571     BC   
                 4 1  |    50.7529   1.116571        EF
                 5 0  |   35.69507   1.116571    A     
                 5 1  |   44.55313   1.116571      CD  
-------------------------------------------------------
Note: Margins sharing a letter in the group label are
      not significantly different at the 5% level.

El paquete multicomp también ofrece visualización simbólica ('pantallas de letras compactas', consulte Algoritmos para pantallas de letras compactas: comparación y evaluación para más detalles) de comparaciones significativas por pares, aunque no las presenta en formato tabular. Sin embargo, tiene un método de trazado que permite mostrar convenientemente los resultados usando diagramas de caja. El orden de presentación también puede modificarse (opción decreasing=), y tiene muchas más opciones para comparaciones múltiples. También está el paquete multcompView que extiende esas funcionalidades.

Aquí está el mismo ejemplo analizado con glht:

library(multcomp)
tuk <- glht(amod, linfct = mcp(tx = "Tukey"))
summary(tuk)          # standard display
tuk.cld <- cld(tuk)   # letter-based display
opar <- par(mai=c(1,1,1.5,1))
plot(tuk.cld)
par(opar)

El tratamiento que comparte la misma letra no es significativamente diferente, en el nivel elegido (predeterminado, 5%).

ingrese la descripción de la imagen aquí

Por cierto, hay un nuevo proyecto, actualmente alojado en R-Forge, que parece prometedor: factorplot . Incluye pantallas basadas en líneas y letras, así como una descripción general de la matriz (a través de un gráfico de nivel) de todas las comparaciones por pares. Puede encontrar un documento de trabajo aquí: factorplot: Mejora de la presentación de contrastes simples en GLM

chl
fuente
¡Muchas gracias por esta respuesta exhaustiva! Probaré esos métodos diferentes tan pronto como tenga unos minutos. ¡Salud!
stragu
Probé la función de paquete multicomp. Cuando uso la función 'cld ()', aparece el error 'Error: sapply (split_names, length) == 2 no es todo VERDADERO' ¿Alguna idea de por qué?
stragu
1
@chtfn Parece que hay un problema con las etiquetas variables. Un vistazo rápido al código fuente indica que este mensaje de error proviene del insert_absorb()que intenta extraer un par de tratamientos. ¿Quizás puede intentar cambiar el separador que utilizó para codificar los niveles de su término de interacción? Sin un ejemplo de trabajo, es difícil saber qué sucedió.
chl
Lo descubrí: tenía puntos en los nombres de mis genotipos y tratamientos, y como qlht () usa un punto para dividir los nombres de los pares, se asustó. Muchas gracias por toda su ayuda, chl! :)
stragu
3
Noté hoy que ahora tengo que añadir console=TRUEen HSD.test()el fin de obtener las tablas, en caso de que alguien trata de esto y no ve ningún resultado. Probablemente una actualización de agricolae.
stragu
2

Hay una función llamada TukeyHSDque, de acuerdo con el archivo de ayuda, calcula un conjunto de intervalos de confianza sobre las diferencias entre las medias de los niveles de un factor con la probabilidad de cobertura especificada por la familia. Los intervalos se basan en la estadística de rango Studentizado, el método de "Diferencia significativa honesta" de Tukey. ¿Hace esto lo que quieres?

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/TukeyHSD.html

smillig
fuente
Gracias por su respuesta. Sí, probé esta función, pero me da listas crudas de comparaciones. Lo que me gustaría es verlos agrupados como en la imagen de mi pregunta, para tener una visión clara de qué grupo difiere en cada grupo, y eventualmente agregar los nombres de grupo en mis gráficos (por ejemplo: a, ab, abc, bc , c)
stragu