Como obter os resultados de um teste post hoc de Tukey HSD em uma tabela que mostra pares agrupados?

13

Eu adoraria realizar um teste post-hoc de TukeyHSD após meu Anova de duas vias com R, obtendo uma tabela contendo os pares classificados agrupados por diferença significativa. (Desculpe pela redação, ainda sou novo em estatísticas.)

Eu gostaria de ter algo parecido com isto:

insira a descrição da imagem aqui

Então, agrupados com estrelas ou letras.

Qualquer ideia? Testei a função HSD.test()do agricolaepacote, mas parece que não lida com tabelas de mão dupla.

stragu
fonte

Respostas:

22

A agricolae::HSD.testfunçã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%).

insira a descrição da imagem aqui

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

chl
fonte
Muito obrigado por esta resposta exaustiva! Vou tentar esses métodos diferentes assim que tiver alguns minutos. Felicidades!
Stragu 7/07
Eu tentei a função do pacote multcomp, coloquei quando uso a função 'cld ()', recebo o erro 'Erro: sapply (split_names, length) == 2 não é tudo VERDADEIRO' Alguma idéia do porquê?
Stragu
1
@chtfn Parece haver um problema com os rótulos de variáveis. Uma rápida olhada no código-fonte indica que essa mensagem de erro vem da 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.
chl
Eu descobri: eu tinha pontos nos nomes dos meus genótipos e tratamentos, e como qlht () usa um ponto para dividir os nomes dos pares, ficou assustado. Muito obrigado por toda a sua ajuda, chl! :)
stragu 11/07/2012
3
Notei hoje que eu tenho agora para adicionar na fim de obter as mesas, no caso de alguém tenta este e não vê resultado. Provavelmente uma atualização de . console=TRUEHSD.test()agricolae
stragu 18/11/2013
2

Existe uma função chamada TukeyHSDque, 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

smillig
fonte
Obrigado pela sua resposta. Sim, tentei essa função, mas ela me fornece listas brutas de comparações. O que eu gostaria é vê-los agrupados como na imagem da minha pergunta, para ter uma visão clara de qual grupo difere para qual grupo e, eventualmente, adicionar os nomes dos grupos nos meus gráficos (por exemplo: a, ab, abc, bc , c)
stragu