Как получить результаты специального теста Tukey HSD в таблице, показывающей сгруппированные пары?

13

Я хотел бы выполнить специальный тест TukeyHSD после моей двусторонней Anova с R, получив таблицу, содержащую отсортированные пары, сгруппированные по значительным различиям. (Извините за формулировку, я все еще новичок со статистикой.)

Я хотел бы иметь что-то вроде этого:

введите описание изображения здесь

Итак, сгруппированы по звездам или буквам.

Любая идея? Я протестировал функцию HSD.test()из agricolaeпакета, но, похоже, она не обрабатывает двусторонние таблицы.

stragu
источник

Ответы:

22

agricolae::HSD.testФункция делает именно это, но вам нужно будет дать ему знать , что вы заинтересованы в перспективе взаимодействия . Вот пример с набором данных 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)

Это дает результаты, показанные ниже:

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 

Они соответствуют тому, что мы получили бы с помощью следующих команд:

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

Пакет multcomp также предлагает символическую визуализацию («компактные буквенные дисплеи», см. « Алгоритмы компактного буквенного отображения: сравнение и оценка» для более подробной информации) значимых парных сравнений, хотя он не представляет их в табличном формате. Тем не менее, он имеет метод построения графиков, который позволяет удобно отображать результаты с помощью коробочных диаграмм. Порядок представления также может быть изменен (опция decreasing=), и он имеет гораздо больше возможностей для нескольких сравнений. Существует также пакет multcompView , который расширяет эти функции.

Вот тот же пример, проанализированный с помощью 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)

Обработка с использованием одной и той же буквы существенно не отличается на выбранном уровне (по умолчанию 5%).

введите описание изображения здесь

Кстати, есть новый проект, в настоящее время размещенный на R-Forge, который выглядит многообещающе: factorplot . Он включает в себя линейные и буквенные дисплеи, а также обзор матрицы (через график уровня) всех парных сравнений. Рабочий документ можно найти здесь: factorplot: Улучшение представления простых контрастов в GLM

хл
источник
Большое спасибо за этот исчерпывающий ответ! Я попробую эти разные методы, как только получу несколько минут. Ура!
Страгу
Я попробовал функцию пакета multcomp, поставил, когда я использую функцию 'cld ()', я получаю ошибку 'Ошибка: sapply (split_names, length) == 2 не все ИСТИНА' Есть идеи, почему?
Страгу
1
@chtfn Кажется, проблема с метками переменных. Беглый взгляд на исходный код показывает, что это сообщение об ошибке, из insert_absorb()которого пытается извлечь пару обработок. Возможно, вы можете попытаться изменить разделитель, который вы использовали для кодирования уровней вашего термина взаимодействия? Без рабочего примера трудно сказать, что случилось.
chl
Я понял это: у меня были точки в именах моих генотипов и обработок, и, поскольку qlht () использует точку для разделения имен пар, она взбесилась. Большое спасибо за вашу помощь, чел! :)
Страгу
3
Я заметил сегодня , что теперь я должен добавить console=TRUEв HSD.test()для того , чтобы получить таблицы, в случае , если кто - то пытается это и не видит никакого результата. Вероятно, обновление agricolae.
Страгу
2

Есть функция, которая называется TukeyHSD, согласно файлу справки, вычисляет набор доверительных интервалов на основе различий между средними уровнями фактора с заданной вероятностью охвата всей семьей. Интервалы основаны на статистике диапазона Studentized, методе «честного значительного различия» Тьюки. Делает ли это то, что вы хотите?

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

smillig
источник
Благодарю за ваш ответ. Да, я пробовал эту функцию, но она дает мне сырые списки сравнений. Я хотел бы видеть их сгруппированными, как показано на изображении в моем вопросе, чтобы иметь четкое представление о том, какая группа отличается от какой группы, и в конечном итоге добавить имена групп на мои графики (например: a, ab, abc, bc , в)
страгу