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%).
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