作者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