作者DrRd (就這樣吧)
站內Statistics
標題[程式] Steel-Dwass procedure in R program
時間Wed Jun 8 13:04:21 2011
[軟體程式類別]:
R2.13
[程式問題]:
Post Hoc for nonparamatric test
[軟體熟悉度]:
中(3個月到1年)
[問題敘述]:
我所處理的資料有三組，而資料分佈都是正偏態，所以我用Kruskal-Wallis做無母數比較
但是用Bonferroni做事後檢定的話因為太保守了，一些想要有差異的組別卻沒有差異
我在網路上看到有人是用Steel-Dwass 這個方法來做post hoc
但是算法我有找到兩個，兩個算出來的結果卻不一樣
大家是否可以給個意見，那一個才是正確的？
[程式範例]:
程式一：從coin package裡面來的(說明pdf:
http://ppt.cc/c(Ix )
### Nemenyi-Damico-Wolfe-Dunn test (joint ranking)
### Hollander & Wolfe (1999), page 244
### (where Steel-Dwass results are given)
if (require("multcomp")) {
NDWD <- oneway_test(length ~ site, data = YOY,
ytrafo = function(data) trafo(data, numeric_trafo = rank),
xtrafo = function(data) trafo(data, factor_trafo = function(x)
model.matrix(~x - 1) %*% t(contrMat(table(x), "Tukey"))),
teststat = "max", distribution = approximate(B = 90000))
### global p-value
print(pvalue(NDWD))
### sites (I = II) != (III = IV) at alpha = 0.01 (page 244)
print(pvalue(NDWD, method = "single-step"))
}
程式二：在網路上找到的，跟書上對照看起來是直接寫公式(
http://ppt.cc/1_An )
並且這個應該是單尾的
Steel.Dwass <- function(data, group)
{
OK <- complete.cases(data, group)
data <- data[OK]
group <- group[OK]
n.i <- table(group)
ng <- length(n.i)
t <- combn(ng, 2, function(ij) {
i <- ij[1]
j <- ij[2]
r <- rank(c(data[group == i], data[group == j]))
R <- sum(r[1:n.i[i]])
N <- n.i[i]+n.i[j]
E <- n.i[i]*(N+1)/2
V <- n.i[i]*n.i[j]/(N*(N-1))*(sum(r^2)-N*(N+1)^2/4)
return(abs(R-E)/sqrt(V))
})
p <- ptukey(t*sqrt(2), ng, Inf, lower.tail=FALSE)
result <- cbind(t, p)
rownames(result) <- combn(ng, 2, paste, collapse=":")
return(result)
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 203.68.96.125