Форум врачей-аспирантов

Здравствуйте, гость ( Вход | Регистрация )

 
Добавить ответ в эту темуОткрыть тему
> Краскел-Уоллис как регрессия, Все про R
comisora
сообщение 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" предполагается)?
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 27.12.2015 - 21:27
Сообщение #2





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(comisora @ 27.12.2015 - 11:28) *
Всем добрый день.
Возник вопрос касательно порядковых данных и применения одно-двух-повторного анализа данных. В 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. Как он связан с сравнением кумулятивных распределений случайных величин?


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
comisora
сообщение 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)".
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 27.12.2015 - 22:55
Сообщение #4





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(comisora @ 27.12.2015 - 21:54) *
Добрый вечер.
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 .


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
comisora
сообщение 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 нет, хотя они обе выдают одно и тоже. Хочется разобраться ЧЯДНТ.
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 28.12.2015 - 23:59
Сообщение #6





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(comisora @ 27.12.2015 - 23:21) *
kruskal.test(y~x1*x2) не работает, а lrm(y~x1*x2) и orm(y~x1*x2) работают. Для одномерного случая результаты обеих функций неплохо бьются с результатами dunn.test. У lrm контрасты аналогичны попарным сравнениям, а у orm нет, хотя они обе выдают одно и тоже. Хочется разобраться ЧЯДНТ.


Простите, но я все равно не понимаю что вы пытаетесь получить. (и при чем тут dunn.test)


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
comisora
сообщение 29.12.2015 - 01:28
Сообщение #7





Группа: Пользователи
Сообщений: 95
Регистрация: 27.12.2015
Пользователь №: 27815



Цитата(p2004r @ 28.12.2015 - 23:59) *
Простите, но я все равно не понимаю что вы пытаетесь получить. (и при чем тут dunn.test)

dunn.test - это критерий Краскела-Уоллиса с функцией множественных парных сравнений.
Я пытаюсь провести двухфакторный анализ порядковых данных пакетом rms. Для однофакторного анализа хватает критерия Краскела-Уоллиса, а для двухфакторного с взаимодействиями уже нет, для этого нужно использовать регрессию - lrm или orm. Для того, чтобы быть уверенным, что это то, что необходимо, тестирую данные функции на однофакторной модели. В одномерном случае функции lrm и критерий Краскела-Уоллиса дают очень схожие, но не одинаковые результаты при попарном сравнении. orm дает значения коэффициентов регрессии, которые совпадают с таковыми при применении функции lrm, но при попарном сравнении никак не соотносятся с парными сравнениями критерия Данна и контрастов lrm.
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 29.12.2015 - 16:01
Сообщение #8





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(comisora @ 29.12.2015 - 01:28) *
dunn.test - это критерий Краскела-Уоллиса с функцией множественных парных сравнений.
Я пытаюсь провести двухфакторный анализ порядковых данных пакетом rms. Для однофакторного анализа хватает критерия Краскела-Уоллиса, а для двухфакторного с взаимодействиями уже нет, для этого нужно использовать регрессию - lrm или orm. Для того, чтобы быть уверенным, что это то, что необходимо, тестирую данные функции на однофакторной модели. В одномерном случае функции lrm и критерий Краскела-Уоллиса дают очень схожие, но не одинаковые результаты при попарном сравнении. orm дает значения коэффициентов регрессии, которые совпадают с таковыми при применении функции lrm, но при попарном сравнении никак не соотносятся с парными сравнениями критерия Данна и контрастов lrm.


а что мешает просто anova(orm()) позвать?


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
comisora
сообщение 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
#Контрасты всё же различаются. Что-то не так делаю.
#Осталось ко всему этому рассчитать силу эффекта.
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
nokh
сообщение 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
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
comisora
сообщение 12.01.2016 - 11:22
Сообщение #11





Группа: Пользователи
Сообщений: 95
Регистрация: 27.12.2015
Пользователь №: 27815



Цитата(nokh @ 11.01.2016 - 22:33) *
Во-первых, 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
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
comisora
сообщение 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')
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 23.01.2016 - 09:29
Сообщение #13





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(comisora @ 23.01.2016 - 02:25) *
Доброй ночи.
Если кому-то интересно. Написал Харреллу, сказал, что это ошибка http://stats.stackexchange.com/q/191063.
Чтобы с Github'a загрузить фикс, использовать следующий код:
Код
require(rms)
getRs('contrast.s', grepo='rms', dir='R', put='source')


И какую разницу показал расчет с помощью поправленного кода?


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
comisora
сообщение 23.01.2016 - 14:11
Сообщение #14





Группа: Пользователи
Сообщений: 95
Регистрация: 27.12.2015
Пользователь №: 27815



Цитата(p2004r @ 23.01.2016 - 10:29) *
И какую разницу показал расчет с помощью поправленного кода?

Добрый день.
Разницы между orm и lrm нет. В частных случаях результаты функций совпадают с результами wilcox.test() и dunn.test().
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 

Добавить ответ в эту темуОткрыть тему