16 votos

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

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

Me gustaría tener algo así:

enter image description here

Así que, agrupados con estrellas o letras.

¿Alguna idea? He probado la función HSD.test() de la agricolae paquete, pero parece que no maneja mesas de doble sentido.

25voto

DavLink Puntos 101

El agricolae::HSD.test hace exactamente eso, pero tendrá que hacerle saber que está interesado en un término de interacción . He aquí un ejemplo con un conjunto de datos de 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 multcomp también ofrece una visualización simbólica ("visualizaciones de letras compactas", véase Algoritmos para la visualización de letras compactas: Comparación y evaluación para más detalles) de las comparaciones significativas por pares, aunque no las presenta en formato tabular. Sin embargo, dispone de un método de trazado que permite mostrar cómodamente los resultados mediante boxplots. El orden de presentación también puede modificarse (opción decreasing= ), y tiene muchas más opciones para las comparaciones múltiples. También está el multcompView que amplía 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)

Los tratamientos que comparten la misma letra no son significativamente diferentes, al nivel elegido (por defecto, 5%).

enter image description here

Por cierto, hay un nuevo proyecto, actualmente alojado en R-Forge, que parece prometedor: factorplot . Incluye visualizaciones basadas en líneas y letras, así como una visión general de la matriz (mediante un gráfico de niveles) de todas las comparaciones por pares. Se puede encontrar un documento de trabajo aquí: factorplot: mejora de la presentación de los contrastes simples en los GLM

0 votos

Muchas gracias por esta respuesta tan exhaustiva. Probaré esos diferentes métodos en cuanto tenga unos minutos. ¡Saludos!

0 votos

He probado la función del paquete multcomp, pero cuando utilizo la función 'cld()' me sale el error 'Error: sapply(split_names, length) == 2 is not all TRUE' ¿Alguna idea de por qué?

1 votos

@chtfn Parece que hay un problema con las etiquetas de las variables. Un rápido vistazo al código fuente indica que este mensaje de error proviene de insert_absorb() que intenta extraer un par de tratamientos. ¿Tal vez pueda 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é ha pasado.

3voto

Jeff Davis Puntos 1999

Hay una función llamada TukeyHSD que, según 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 familia. Los intervalos se basan en el estadístico de rango estudiado, el método de "diferencia significativa honesta" de Tukey. ¿Hace esto lo que usted quiere?

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

0 votos

Gracias por su respuesta. Sí, he probado esta función, pero me da listas crudas de comparaciones. Lo que me gustaría es verlas agrupadas como en la imagen de mi pregunta, para tener una visión clara de qué grupo difiere de qué grupo, y eventualmente añadir los nombres de los grupos en mis gráficos (por ejemplo: a, ab, abc, bc, c)

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X