Краскел-Уоллис как регрессия, Все про R |
Здравствуйте, гость ( Вход | Регистрация )
Краскел-Уоллис как регрессия, Все про R |
27.12.2015 - 11:28
Сообщение
#1
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
Всем добрый день.
Возник вопрос касательно порядковых данных и применения одно-двух-повторного анализа данных. В R есть пакет "dunn.test" с помощью которого можно провести одномерный вариант с попарными сравнениями (dunn.test(x,g, "BY") - обычный KW. Есть пакет rms с функцией orm, частным случаем которой должен быть KW для одного вектора с порядковыми данными с разделением по подгруппам. Результат функции он выдает в виде регрессии с коэффициентами, однако возникло желание привести результат к "привычному" виду. Возникла мысль использовать конструкцию contrasts(Glm(M), list(g="a"), list(g="b")) потом "a" и "c" и "b" и "c". Результат схож (но не идентичен) с результатом dunn.test, однако уверенности в том, что так можно делать нет. 1. Насколько правомерно так тестировать контрасты для данных, измеренных в порядковой шкале? 2. Можно ли таким способом тестировать двух-трех-многовходовые порядковые данные для вычленения парных различий (использование поправки "BY" предполагается)? |
|
27.12.2015 - 21:27
Сообщение
#2
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Всем добрый день. Возник вопрос касательно порядковых данных и применения одно-двух-повторного анализа данных. В R есть пакет "dunn.test" с помощью которого можно провести одномерный вариант с попарными сравнениями (dunn.test(x,g, "BY") - обычный KW. Есть пакет rms с функцией orm, частным случаем которой должен быть KW для одного вектора с порядковыми данными с разделением по подгруппам. Результат функции он выдает в виде регрессии с коэффициентами, однако возникло желание привести результат к "привычному" виду. Возникла мысль использовать конструкцию contrasts(Glm(M), list(g="a"), list(g="b")) потом "a" и "c" и "b" и "c". Результат схож (но не идентичен) с результатом dunn.test, однако уверенности в том, что так можно делать нет. 1. Насколько правомерно так тестировать контрасты для данных, измеренных в порядковой шкале? 2. Можно ли таким способом тестировать двух-трех-многовходовые порядковые данные для вычленения парных различий (использование поправки "BY" предполагается)? 1. А что такое "одно-двух-повторного анализа данных"? 2. Как он связан с сравнением кумулятивных распределений случайных величин? |
|
27.12.2015 - 21:54
Сообщение
#3
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
Добрый вечер.
1. Однофакторный анализ (y~x). Соответственно можно думать о (y~x1*x2) и тп. 2. http://r.789695.n4.nabble.com/kruskal-wall...-td4288008.html - "The Kruskal-Wallis test is a special case of the proportional odds ordinal logistic model. You can get any contrast you want by testing regression coefficients. In a couple of weeks the rms package's contrast function will allow for individual confidence intervals of effects that together have a 0.05 type I error, by using the multcomp package (called automatically from contrast.rms)". |
|
27.12.2015 - 22:55
Сообщение
#4
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Добрый вечер. 1. Однофакторный анализ (y~x). Соответственно можно думать о (y~x1*x2) и тп. 2. http://r.789695.n4.nabble.com/kruskal-wall...-td4288008.html - "The Kruskal-Wallis test is a special case of the proportional odds ordinal logistic model. You can get any contrast you want by testing regression coefficients. In a couple of weeks the rms package's contrast function will allow for individual confidence intervals of effects that together have a 0.05 type I error, by using the multcomp package (called automatically from contrast.rms)". Тогда надо использовать просто ?kruskal.test . |
|
27.12.2015 - 23:21
Сообщение
#5
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
kruskal.test(y~x1*x2) не работает, а lrm(y~x1*x2) и orm(y~x1*x2) работают. Для одномерного случая результаты обеих функций неплохо бьются с результатами dunn.test. У lrm контрасты аналогичны попарным сравнениям, а у orm нет, хотя они обе выдают одно и тоже. Хочется разобраться ЧЯДНТ.
|
|
28.12.2015 - 23:59
Сообщение
#6
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
kruskal.test(y~x1*x2) не работает, а lrm(y~x1*x2) и orm(y~x1*x2) работают. Для одномерного случая результаты обеих функций неплохо бьются с результатами dunn.test. У lrm контрасты аналогичны попарным сравнениям, а у orm нет, хотя они обе выдают одно и тоже. Хочется разобраться ЧЯДНТ. Простите, но я все равно не понимаю что вы пытаетесь получить. (и при чем тут dunn.test) |
|
29.12.2015 - 01:28
Сообщение
#7
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
Простите, но я все равно не понимаю что вы пытаетесь получить. (и при чем тут dunn.test) dunn.test - это критерий Краскела-Уоллиса с функцией множественных парных сравнений. Я пытаюсь провести двухфакторный анализ порядковых данных пакетом rms. Для однофакторного анализа хватает критерия Краскела-Уоллиса, а для двухфакторного с взаимодействиями уже нет, для этого нужно использовать регрессию - lrm или orm. Для того, чтобы быть уверенным, что это то, что необходимо, тестирую данные функции на однофакторной модели. В одномерном случае функции lrm и критерий Краскела-Уоллиса дают очень схожие, но не одинаковые результаты при попарном сравнении. orm дает значения коэффициентов регрессии, которые совпадают с таковыми при применении функции lrm, но при попарном сравнении никак не соотносятся с парными сравнениями критерия Данна и контрастов lrm. |
|
29.12.2015 - 16:01
Сообщение
#8
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
dunn.test - это критерий Краскела-Уоллиса с функцией множественных парных сравнений. Я пытаюсь провести двухфакторный анализ порядковых данных пакетом rms. Для однофакторного анализа хватает критерия Краскела-Уоллиса, а для двухфакторного с взаимодействиями уже нет, для этого нужно использовать регрессию - lrm или orm. Для того, чтобы быть уверенным, что это то, что необходимо, тестирую данные функции на однофакторной модели. В одномерном случае функции lrm и критерий Краскела-Уоллиса дают очень схожие, но не одинаковые результаты при попарном сравнении. orm дает значения коэффициентов регрессии, которые совпадают с таковыми при применении функции lrm, но при попарном сравнении никак не соотносятся с парными сравнениями критерия Данна и контрастов lrm. а что мешает просто anova(orm()) позвать? |
|
11.01.2016 - 19:08
Сообщение
#9
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
Добрый день.
Ничего не мешает запустить. Выдаёт хи, df и p, как и должно выдавать. Лучше, наверное, сформулировать задачку. Цитата #Пусть x - показатель, измеренный в порядковой шкале от 0 до 10. Есть 3 группы q,w,e. Задача: выяснить, есть ли различия в распределении между группами; если различия есть, то между какими группами именно; рассчитать силу эффекта, доверительный интервал силы эффекта для общего теста и парных сравнений.
#Дополнение к задаче: учесть влияние фактора f. > x1<-c(1,2,1,2,3,4,0,4,2,1) > x2<-c(5,4,6,2,5,4,7,5,6,8) > x3<-c(7,8,6,9,6,8,7,7,9,10) > x<-c(x1,x2,x3) > f <- factor(rbinom(n = 30,size = 1,prob = 0.5), labels = c("no","yes")) > g <- gl(n = 3,k = 10,length = 30,labels = c("q","w","e")) > t<-data.frame(g,f,x) > attach(t) #Проверяем различие распределений в подгруппах q,w,e. > kruskal.test(x~g) #Kruskal-Wallis rank sum test #data: x by g #Kruskal-Wallis chi-squared = 21.674, df = 2, p-value = 1.966e-05 #Подгруппы различны. Чтобы понять за счет чего, можно провести парные сравнения. > dunn.test(x,g,"BY",kw = FALSE) # Comparison of x by g # (Benjamini-Yekuteili) #Col Mean-| #Row Mean | q w #---------+---------------------- # w | -2.516784 # | 0.0163* # | # e | -4.650301 -2.133517 # | 0.0000* 0.0301 # Теперь это все это нужно сделать в более общем виде. > m<-orm(x~g) > anova(m) # Wald Statistics Response: x # # Factor Chi-Square d.f. P # g 21.1 2 <.0001 # TOTAL 21.1 2 <.0001 > ctr<-contrast(m, list(g=c("q","q","w")),list(g=c("w","e","e"))) > ctr # Contrast S.E. Lower Upper Z Pr(>|z|) #1 1 -1.8661910 1.300574 -4.415270 0.6828876 -1.43 0.1513 #2 2 -4.9520552 1.566679 -8.022690 -1.8814200 -3.16 0.0016 #3* 3 -0.8692409 1.054338 -2.935705 1.1972231 -0.82 0.4097 # #Redundant contrasts are denoted by * # #Confidence intervals are 0.95 individual intervals > p.adjust(ctr$Pvalue,"BY",3) # 1 2 3 #0.416119667 0.008651706 0.751094958 #Результат контрастов отличается от того, что было вызвано функцией dunn.test(). Может не та команда? Попробуем lrm. > m<-lrm(x~g) > anova(m) # Wald Statistics Response: x # # Factor Chi-Square d.f. P # g 21.1 2 <.0001 # TOTAL 21.1 2 <.0001 > ctr<-contrast(m, list(g=c("q","q","w")),list(g=c("w","e","e"))) > ctr # Contrast S.E. Lower Upper Z Pr(>|z|) #1 1 -4.082814 1.300574 -6.631893 -1.533736 -3.14 0.0017 #2 2 -7.168679 1.566679 -10.239314 -4.098043 -4.58 0.0000 #3* 3 -3.085864 1.054338 -5.152328 -1.019400 -2.93 0.0034 # #Redundant contrasts are denoted by * # #Confidence intervals are 0.95 individual intervals > p.adjust(ctr$Pvalue,"BY",3) # 1 2 3 #4.658137e-03 2.610279e-05 6.278054e-03 #Результат более схож с результатами dunn.test(). #Проведём двухфакторный анализ порядковых данных. > kruskal.test(x~g*f) #Error in kruskal.test.formula(x ~ g * f) : # 'formula'должна быть в виде отклик ~ группа > anova(lrm(x~g*f)) # Wald Statistics Response: x # # Factor Chi-Square d.f. P # g (Factor+Higher Order Factors) 21.70 4 0.0002 # All Interactions 0.76 2 0.6845 # f (Factor+Higher Order Factors) 1.29 3 0.7311 # All Interactions 0.76 2 0.6845 # g * f (Factor+Higher Order Factors) 0.76 2 0.6845 # TOTAL > m<-lrm(x~g*f) > contrast(m, list(g="q", f="no"), list(g="w",f="no")) # f Contrast S.E. Lower Upper Z Pr(>|z|) #1 no -4.734745 1.621549 -7.912923 -1.556566 -2.92 0.0035 # #Confidence intervals are 0.95 individual intervals > m<-orm(x~g*f) > contrast(m, list(g="q", f="no"), list(g="w",f="no")) # f Contrast S.E. Lower Upper Z Pr(>|z|) #1 no -2.963172 1.621549 -6.141351 0.2150061 -1.83 0.0676 # #Confidence intervals are 0.95 individual intervals #Контрасты всё же различаются. Что-то не так делаю. #Осталось ко всему этому рассчитать силу эффекта. |
|
11.01.2016 - 22:33
Сообщение
#10
|
|
Группа: Пользователи Сообщений: 1202 Регистрация: 13.01.2008 Из: Челябинск Пользователь №: 4704 |
Добавлю свои 5 копеек.
Во-первых, dunn.test - это не "критерий Краскела - Уоллиса с функцией множественных парных сравнений". Это - критерий Данна. Он был предложен в Оливом Данном в 1964 г., однако распространение получил только после 1988 г, когда без указания авторства был описан в ставшей популярной монографии [Siegel, S., & Castellan, N. J. (1988). Nonparametric statistics for the behavioral sciences (2nd ed.) New York: McGraw-Hill] и потому в некоторых работах необоснованно называется именами её авторов как Siegel-Castellan test. Это - не аналог контрастов, это - прямой ранговый аналог критерия Холма - Шидака (Holm-Šidak multiple t-test), относящийся к процедурам пошагового спуска. Удерживает ошибку семейства гипотез (Familywise Error Rate, FWER) на заданном уровне. Эффективный, но отчасти консервативный метод. Используется для непараметрических множественных сравнений в пакетах Statistica и Graphpad Prizm. Таким образом, результаты линейных контрастов на рангах и результаты теста Данна и не должны быть идентичными - это разные методы. Более того, линейные контрасты по Шеффе (и какой-то их ранговый аналог, который в R получаете вы) относятся к запланированным сравнениям, тогда как критерий Данна - к множественным апостериорным (post-hoc) сравнениям. При числе групп k запланированных сравнений может быть только k-1, тогда как апостериорных - куда больше: 0,5k*(k-1). Поэтому мощность запланированных сравнений всегда выше и p должны быть меньше. Во-вторых, помимо критерия Данна для пост-хок сравнений после Краскела-Уоллиса есть ещё куча критериев: Неменьи - прямой ранговый аналог множественных сравнений по Шеффе (не путать с контрастами по Шеффе), часто некорректно называемый Schaich-Hamerle test, Стила - Двасса - прямой ранговый аналог Tukey's HSD, Коновера - Инмана - прямой ранговый аналог Fisher's LSD др. Поэтому полагаю, что результат, который вы получаете в контрастах на рангах, также имеет право на существование. Это будет подход в духе Коновера, который вообще рекомендовал всегда перед использованием параметрических тестов предварительно ранжировать значения, а расчёты проводить с рангами - типа для увеличения робастности (Conover W.J., Iman R.L. Rank Transformations as a Bridge Between Parametric and Nonparametric Statistics // The American Statistician. 1981. V. 35. N 3. P. 124-129.). Сообщение отредактировал nokh - 11.01.2016 - 23:08 |
|
12.01.2016 - 11:22
Сообщение
#11
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
Во-первых, dunn.test - это не "критерий Краскела - Уоллиса с функцией множественных парных сравнений". Это - критерий Данна. Однако https://cran.r-project.org/web/packages/dun.../dunn.test.pdf: dunn.test (x, g=NA, method=p.adjustment.methods, kw=TRUE, label=TRUE, wrap=FALSE, table=TRUE, list=FALSE, rmc=FALSE, alpha=0.05). Идентичными они быть, конечно, не должны. Но резко различаться не должны тоже. Мне пока не понятно, чем именно отличаются olm() и rlm(). Сообщение отредактировал comisora - 12.01.2016 - 11:24 |
|
23.01.2016 - 02:25
Сообщение
#12
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
Доброй ночи.
Если кому-то интересно. Написал Харреллу, сказал, что это ошибка http://stats.stackexchange.com/q/191063. Чтобы с Github'a загрузить фикс, использовать следующий код: Код require(rms)
getRs('contrast.s', grepo='rms', dir='R', put='source') |
|
23.01.2016 - 09:29
Сообщение
#13
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Доброй ночи. Если кому-то интересно. Написал Харреллу, сказал, что это ошибка http://stats.stackexchange.com/q/191063. Чтобы с Github'a загрузить фикс, использовать следующий код: Код require(rms) getRs('contrast.s', grepo='rms', dir='R', put='source') И какую разницу показал расчет с помощью поправленного кода? |
|
23.01.2016 - 14:11
Сообщение
#14
|
|
Группа: Пользователи Сообщений: 95 Регистрация: 27.12.2015 Пользователь №: 27815 |
|
|