Цитата(DrgLena @ 12.02.2013 - 13:47)

Да, относительно референтной группы интересно получить результат, но у меня не сходится ручной расчет по этому коду при сравнении конкретно табл 2хk. Пример из учебника есть тут, на стр 175-176
http://books.google.com.ua/books?id=EpEfWv...ple&f=false Я получаю руками 0,604 как в книге и 0,310 если поменять местами группы, а пакет выдает 0,578 и 0,473. Или иранский студент, автор кода, что то не так делает, но результат не сходится. Проверить можно только в SAS , но профессиональный статистик, работающий под прикрытием TL и который нам поведал про RIDI о результатх умалчивает. Я надеюсь, что Вы, уважаемы p2012 сможете свой код написать с понятным вводом таблицы сопряженности, может я не верно ввожу таблицу, не указываю, что она cross.
что делает иранский студент легко увидеть набрав ?ridit.raw
Код
## The function is currently defined as
function (x, g, ref = NULL)
{
x = as.numeric(x)
x = as.vector(x)
g = as.factor(g)
levels = levels(g)
levels(g) = 1:length(levels)
g = as.vector(g)
g = as.character(g)
code = is.numeric(ref)
ref = as.vector(ref)
ref = as.character(ref)
if (length(ref) > 1) {
x = c(x, ref)
g = c(g, rep(".ref", length(ref)))
levels = c(".ref", levels)
}
crosstab = t(as.matrix(table(x, g)))
rownames(crosstab) = levels
refindex = NULL
if (length(ref) == 1) {
if (!code)
refindex = which(levels == ref)
if (code && ref >= 1 && ref <= nrow(crosstab))
refindex = as.numeric(ref)
}
else if (length(ref) > 1)
refindex = which(levels == ".ref")
if (length(refindex) != 0)
refrow = crosstab[refindex, ]
else refrow = apply(crosstab, 2, sum)
if (length(refindex) == 0)
msg = paste("Reference: Total of all groups", sep = "")
else msg = paste("Reference: Group = ", refindex, ", Label = ",
levels[refindex], sep = "")
nref = sum(refrow)
ridit = 0.5 * refrow[1]/nref
for (i in 2:length(refrow)) {
iridit = (sum(refrow[1:i - 1]) + 0.5 * refrow[i])/nref
ridit = c(ridit, iridit)
}
n = apply(crosstab, 1, sum)
meanRidit = c()
for (i in 1:nrow(crosstab)) {
itable = crosstab[i, ]
meanRidit = c(meanRidit, sum(ridit * itable)/n[i])
}
n0 = sum(n)
rbar0 = sum(n * meanRidit)/n0
t = apply(crosstab, 2, sum)
f = 1 - (sum(t * (t - 1) * (t + 1)))/(n0 * (n0 - 1) * (n0 +
1))
teststatistic = (12 * n0 * sum(n * (meanRidit - rbar0)^2))/((n0 +
1) * f)
testdf = nrow(crosstab) - 1
pvalue = pchisq(q = teststatistic, df = testdf, lower.tail = FALSE)
if (length(ref) == 0)
ref = NULL
names(meanRidit) = rownames(crosstab)
output = list(MeanRidit = meanRidit, TestStatistic = teststatistic,
df = testdf, Sig = pvalue, x = x, g = g, ref = ref, crosstab = crosstab,
msg = msg)
class(output) <- c("ridit", class(output))
output
}
если посмотреть в сравнении с немодифицированным тестом, то сумма сходится

Код
> str(ridit(as.table(matrix(c(1,3,7,46,0,0,9,58,0,2,22,43,0,2,22,43,1,6,14,46), nrow = 5, byrow = TRUE)),1))
List of 9
$ MeanRidit : Named num [1:5] 0.535 0.573 0.462 0.462 0.473
..- attr(*, "names")=8322456 [1:5] "A" "B" "C" "D" ...
$ TestStatistic: num 13.2
$ df : num 4
$ Sig : num 0.0104
$ x : num [1:325] 1 2 2 2 3 3 3 3 3 3 ...
$ g :8322456 [1:325] "1" "1" "1" "1" ...
$ ref : NULL
$ crosstab : 'table' int [1:5, 1:4] 1 0 0 0 1 3 0 2 2 6 ...
..- attr(*, "dimnames")=List of 2
.. ..$ g:8322456 [1:5] "A" "B" "C" "D" ...
.. ..$ x:8322456 [1:4] "1" "2" "3" "4"
$ msg :8322456 "Reference: Total of all groups"
- attr(*, "class")=8322456 [1:2] "ridit" "list"
> ridit(as.table(matrix(c(1,3,7,46,0,0,9,58,0,2,22,43,0,2,22,43,1,6,14,46), nrow = 5, byrow = TRUE)),1)$x
[1] 1 2 2 2 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[38] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
[75] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[112] 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[149] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[186] 4 4 4 4 4 4 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4
[223] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1
[260] 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[297] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
> ridit(as.table(matrix(c(1,3,7,46,0,0,9,58,0,2,22,43,0,2,22,43,1,6,14,46), nrow = 5, byrow = TRUE)),1)$g
[1] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[19] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[37] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[55] "1" "1" "1" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
[73] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
[91] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
[109] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "3" "3"
[127] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
[145] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
[163] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
[181] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "4" "4" "4" "4" "4" "4" "4"
[199] "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4"
[217] "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4"
[235] "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4"
[253] "4" "4" "4" "4" "4" "4" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5"
[271] "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5"
[289] "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5"
[307] "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5"
[325] "5"
> kruskal.test(ridit(as.table(matrix(c(1,3,7,46,0,0,9,58,0,2,22,43,0,2,22,43,1,6,14,46), nrow = 5, byrow = TRUE)),1)$x, as.factor(ridit(as.table(matrix(c(1,3,7,46,0,0,9,58,0,2,22,43,0,2,22,43,1,6,14,46), nrow = 5, byrow = TRUE)),1)$g))
Kruskal-Wallis rank sum test
data: ridit(as.table(matrix(c(1, 3, 7, 46, 0, 0, 9, 58, 0, 2, 22, 43, and as.factor(ridit(as.table(matrix(c(1, 3, 7, 46, 0, 0, 9, 58, 0, 0, 2, 22, 43, 1, 6, 14, 46), nrow = 5, byrow = TRUE)), 1)$x and 2, 22, 43, 0, 2, 22, 43, 1, 6, 14, 46), nrow = 5, byrow = TRUE)), ridit(as.table(matrix(c(1, 3, 7, 46, 0, 0, 9, 58, 0, 2, 22, 43, and 1)$g)
Kruskal-Wallis chi-squared = 13.1813, df = 4, p-value = 0.01042
> ridit(as.table(matrix(c(1,3,7,46,0,0,9,58,0,2,22,43,0,2,22,43,1,6,14,46), nrow = 5, byrow = TRUE)),1)
Ridit Analysis:
Group Label Mean Ridit
----- ----- ----------
1 A 0.5351
2 B 0.5729
3 C 0.4621
4 D 0.4621
5 E 0.4731
Reference: Total of all groups
chi-squared = 13.1813, df = 4, p-value = 0.01042
учебник пока не смотрел...