Тест Крускала-Уоллиса с подробными сведениями о попарных сравнениях

Стандартная статистика:: модуль kruskal.test позволяет вычислить тест kruskal-wallis на наборе данных:

>>> data(diamonds)
>>> kruskal.test(price~carat, data=diamonds)

Kruskal-Wallis rank sum test

data:  price by carat by color 
Kruskal-Wallis chi-squared = 50570.15, df = 272, p-value < 2.2e-16

Это правильно, это дает мне вероятность того, что все группы в данных имеют одинаковое среднее значение.

Тем не менее, я хотел бы иметь детали для каждого сравнения пар, например, если алмазы цветов D и E имеют одинаковую среднюю цену, как и некоторые другие программы (SPSS), когда вы запрашиваете тест Kruskal.

Я нашел kruskalmc из пакета pgirmess, который позволяет мне делать то, что я хочу сделать:

> kruskalmc(diamonds$price, diamonds$color)
Multiple comparison test after Kruskal-Wallis 
p.value: 0.05 
Comparisons
      obs.dif critical.dif difference
D-E  571.7459     747.4962      FALSE
D-F 2237.4309     751.5684       TRUE
D-G 2643.1778     726.9854       TRUE
D-H 4539.4392     774.4809       TRUE
D-I 6002.6286     862.0150       TRUE
D-J 8077.2871    1061.7451       TRUE
E-F 2809.1767     680.4144       TRUE
E-G 3214.9237     653.1587       TRUE
E-H 5111.1851     705.6410       TRUE
E-I 6574.3744     800.7362       TRUE
E-J 8649.0330    1012.6260       TRUE
F-G  405.7470     657.8152      FALSE
F-H 2302.0083     709.9533       TRUE
F-I 3765.1977     804.5390       TRUE
F-J 5839.8562    1015.6357       TRUE
G-H 1896.2614     683.8760       TRUE
G-I 3359.4507     781.6237       TRUE
G-J 5434.1093     997.5813       TRUE
H-I 1463.1894     825.9834       TRUE
H-J 3537.8479    1032.7058       TRUE
I-J 2074.6585    1099.8776       TRUE

Однако этот пакет допускает только одну категориальную переменную (например, я не могу изучить цены, сгруппированные по цвету и карат, как я могу сделать с kruskal.test), и я ничего не знаю о пакете pgirmess, поддерживается ли оно или нет, или если оно проверено.

Можете ли вы порекомендовать мне пакет для выполнения теста Kruskal-Wallis, который возвращает данные для каждого сравнения? Как бы вы справились с этой проблемой?

Ответ 1

Еще один подход, помимо kruskal:: agricolae, упомянутый Marek, - это тест Nemenyi-Damico-Wolfe-Dunn, реализованный на странице справки oneway_test в пакете монет, который использует multcomp. Используя настройку hadley и уменьшая значение B = для приближенной() функции, чтобы она заканчивалась за конечное время:

#updated translation of help page implementation of NDWD
NDWD <- 
    independence_test(dv ~ iv, data = sum_codings1, distribution = approximate(B = 10000), 
                          ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo), 
                          xtrafo = mcp_trafo(iv = "Tukey"))


    ### global p-value
    print(pvalue(NDWD))

    ### sites (I = II) != (III = IV) at alpha = 0.01 (page 244)
    print(pvalue(NDWD, method = "single-step"))

Более стабильные результаты для этого более крупного набора данных могут потребовать увеличения значения B и увеличения терпения пользователя.

Январь: 2012: Недавно была опубликована публикация в R-help, в которой были получены неожиданные результаты из этого метода, поэтому я отправил это письмо сопровождающему. Марк Диффорд сказал, что он подтвердил проблемы и предложил альтернативные тесты с пакетом nparcomp: https://stat.ethz.ch/pipermail/r-help/2012-January/300100.html

На той же неделе была еще пара предложений по rhelp для пост-hoc-контрастов с испытаниями KW: kruskalmc, предложенный Марио Гарридо Эскудеро и rms::polr, за которым следует rms::contrasts, предложенный Фрэнком Харрелом https://stat.ethz.ch/pipermail/r-help/2012-January/300329.html

Ноябрь 2015: Согласитесь с toto_tico, что код страницы монеты был изменен за прошедшие годы. Страница справки ?independence_test теперь предлагает многовариантный тест KW, а справочная страница ?oneway_test заменяет предыдущую реализацию кодом выше, чем функция independence_test.

Ответ 2

Вы можете использовать пакет PMCMR. Здесь есть больше информации об этом.

Spelling_Grades <- c(90,87,89,90,75,88,97,99,78,85,72,76,77,79,70)
Methods <- c("A","A","A","A","B","B","B","B","B","B","C","C","C","C","C")
kruskalmc(Spelling_Grades~Methods)

#This method doesn't accept characters that why I've changed the methods to integer
Methods <- c(1,1,1,1,2,2,2,2,2,2,3,3,3,3,3)
posthoc.kruskal.nemenyi.test(Spelling_Grades~Methods)

Два вышеприведенных метода дают одинаковые результаты.

Ответ 3

К сожалению, я не знаю о такой функции. Если его уже нет, было бы интересной задачей построить функцию, которая возвращает матрицу со всеми сравнениями попарного лечения. Контраст считается существенным, если выполняется неравенство

alt text http://www.statsdirect.com/help/image/stat0199_wmf.gif

где T - статистическая оценка Крускала-Уоллиса для k выборок, S ^ 2 - знаменатель статистики T, N - общее число (все ni), а Ri - сумма рангов (из всех сэмплов) для i-го образца, t - квантиль из распределения Стьюдента t на Nk степенях свободы.

Я знаю, что не очень помог:)
Я также жду ответа.

Ответ 4

Я бы подумал, что вы сможете сделать следующее:

data(diamonds, package = "ggplot2")

library(coin)
library(multcomp)

kt <- kruskal_test(price ~ clarity, data = diamonds)
glht(kt, mcp(clarity = "Tukey"))

Но похоже, что multcomp не поддерживает объекты coin (пока?).