Краскел-Уоллис как регрессия, Все про R Опции

[comisora]

Всем добрый день.
Возник вопрос касательно порядковых данных и применения одно-двух-повторного анализа данных. В 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" предполагается)?

[comisora]

Добрый день.
Ничего не мешает запустить. Выдаёт хи, 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
#Контрасты всё же различаются. Что-то не так делаю.
#Осталось ко всему этому рассчитать силу эффекта.