A agricolae::HSD.test
função faz exatamente isso, mas você precisará informar que está interessado em um termo de interação . Aqui está um exemplo com um conjunto de dados 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)
Isso fornece os resultados mostrados abaixo:
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
Eles correspondem ao que obteríamos com os seguintes 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.
O pacote multcomp também oferece visualização simbólica ('exibição compacta de letras', consulte Algoritmos para exibições compactas de letras: comparação e avaliação para obter mais detalhes) de comparações significativas em pares, embora não as apresente em um formato tabular. No entanto, possui um método de plotagem que permite exibir convenientemente os resultados usando gráficos de caixa. A ordem da apresentação também pode ser alterada (opção decreasing=
) e possui muito mais opções para várias comparações. Há também o pacote multcompView que amplia essas funcionalidades.
Aqui está o mesmo exemplo analisado com 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)
O tratamento que compartilha a mesma letra não é significativamente diferente, no nível escolhido (padrão, 5%).
Aliás, há um novo projeto, atualmente hospedado no R-Forge, que parece promissor: factorplot . Inclui exibições baseadas em linhas e letras, bem como uma visão geral da matriz (através de um gráfico de nível) de todas as comparações aos pares. Um documento de trabalho pode ser encontrado aqui: factorplot: Melhorando a apresentação de contrastes simples nos GLMs
insert_absorb()
qual tenta extrair par de tratamentos. Você pode tentar alterar o separador usado para codificar os níveis do seu termo de interação? Sem um exemplo prático, é difícil dizer o que aconteceu.console=TRUE
HSD.test()
agricolae
Existe uma função chamada
TukeyHSD
que, de acordo com o arquivo de ajuda, calcula um conjunto de intervalos de confiança nas diferenças entre as médias dos níveis de um fator com a probabilidade de cobertura específica da família. Os intervalos são baseados na estatística do intervalo Studentized, no método "Honest Significant Difference" de Tukey. Isso faz o que você quer?http://stat.ethz.ch/R-manual/R-patched/library/stats/html/TukeyHSD.html
fonte