如何在 Kruskal Wallis 事后检验后获得显着性代码答案

【问题标题】：how to get significance codes after Kruskal Wallis post hoc test如何在 Kruskal Wallis 事后检验后获得显着性代码
【发布时间】：2013-06-04 16:33:01
【问题描述】：

在与 Kruskall wallis 检验进行成对比较后，有没有办法获得显着性代码？显着性代码是指分配给总体的字母代码，以指示差异显着的位置。

对于传统的 anova，可以使用 agricolae 库中的 HSD.test 执行此类测试，但对于 anova 的非参数对应物，我找不到任何东西。

一个小玩具示例：

dv  <-  c(runif(100, 5.0, 10))
iv  <-  as.factor( c(rep("I", 10),  rep("II", 10),  rep("III", 10),  rep("IV", 10), rep("V", 10),
                    rep("VI", 10), rep("VII", 10), rep("VIII", 10), rep("IX", 10), rep("X", 10)))

df  <-  data.frame(dv, iv)

# with anova
library(agricolae)
aov.000  <-  aov(dv ~ iv,  data=df)
HSD.test(aov.000, "iv")

# after KW test: 
(kt  <-  kruskal.test(dv ~ iv,  data=df))

library(coin)
library(multcomp)
NDWD <- oneway_test(dv ~ iv, data = df,
        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=1000))

### global p-value
print(pvalue(NDWD))

### sites (I = II) != (III = IV) at alpha = 0.01 (page 244)
print(pvalue(NDWD, method = "single-step"))

【问题讨论】：

可能重复stackoverflow.com/questions/2478272/…

标签： r statistics

【解决方案1】：

因为它可能对其他人有用，以下似乎有效（使用multcompView 库）：

library(multcompView)
mat  <-  data.frame(print(pvalue(NDWD, method = "single-step")))
(a   <-  c(mat[, 1]));  names(a)  <-  rownames(mat)
multcompLetters(a)

以下方法也可以使用：

test  <-  pairwise.wilcox.test(dv, iv, p.adj="bonferroni", exact=FALSE)
# test  <-  pairwise.wilcox.test(et.ef, s.t, p.adj="holm", exact=FALSE)

library(multcompView)
test$p.value
library(reshape)
(a <- melt(test$p.value))
a.cc  <-  na.omit(a)
a.pvals  <-  a.cc[, 3]
names(a.pvals)  <-  paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)

【讨论】：

太棒了！总有一天这对我来说会很方便。
运行mat <- data.frame(print(pvalue(NDWD, method = "single-step")))后出现以下错误Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) : cannot coerce class ""ci"" to a data.frame

【解决方案2】：

您还可以使用 rcompanion 包中的 cldList 函数（请参阅https://rcompanion.org/rcompanion/d_06.html）。示例：

k_test <- k_test$res

library(rcompanion)

cldList(comparison = k_test$Comparison,
    p.value    = PT$P.adj,
    threshold  = 0.05)


Error: No significant differences.

我将它与 Dunn post-hoc 结合使用，效果很好。

【讨论】：

【解决方案3】：

您至少可以使用 multicomp 包以图形方式完成：

dv  <-  c(runif(100, 5.0, 10))
iv  <-  as.factor( c(rep("I", 10),  rep("II", 10),  rep("III", 10),  rep("IV", 10), rep("V", 10),
                rep("VI", 10), rep("VII", 10), rep("VIII", 10), rep("IX", 10), rep("X", 10)))
df  <-  data.frame(dv, iv)
anova_results  <-  aov(dv ~ iv,  data=df)
library(multcomp)
tuk <- glht(anova_results, linfct = mcp(iv = "Tukey"))
summary(tuk)          # standard display
tuk.cld <- cld(tuk)   # letter-based display
opar <- par(mai=c(1,1,1.5,1))
plot(tuk.cld)
par(opar)

当然给定你随机生成的数据，结果图不是很有趣，但会给你分组-

这是我的情节之一，使用相同的方法：

最后，如果你不想要图形，你可以深入包中，很容易找到存储分组信息的字符串，以便在其他地方使用。

【讨论】：

我不想使用方差分析。这个例子确实可以用 anaova 分析（一切都是正态分布的，方差相等），但是对于我正在处理的真实数据，我需要使用非参数测试。
是的，很抱歉没有发现您想要进行非参数测试的事实（我已经有一段时间没有进行 KW 测试了）。截至目前，我不知道如何在非参数对比测试上绘制类似的图。但是，您可能知道，特别是对于单向比较等简单设计，您可能会尝试对秩执行方差分析和对比度测试。如果有的话，您可以将您的初始 KW 测试与排名上的方差分析进行比较，以查看结果有何不同，如果没有，则应用排名上的方差分析以便能够执行事后对比测试。

【解决方案4】：

如果您想要 Kruskal 测试的紧凑型字母显示，同一个库 agricolae 似乎允许它与函数 kruskal 一起使用。使用您自己的数据：

library(agricolae)
kruskal(df$dv, df$iv, group=TRUE, p.adj="bonferroni")$groups
####    trt  means M
#### 1  VI    59.2 a
#### 2  VII   57.0 a
#### 3  IX    56.4 a
#### 4  II    55.0 a
#### ...

（好吧，在这个例子中，这些组不被认为是不同的......）

【讨论】：