两个以上样本中位数差异的假设检验

机器算法验证 r 假设检验 多重比较 意思是 中位数
2022-02-17 02:12:13

问题

三组人的测试分数在 R 中保存为单独的向量。

set.seed(1)
group1 <- rnorm(100, mean = 75, sd = 10)
group2 <- rnorm(100, mean = 85, sd = 10)
group3 <- rnorm(100, mean = 95, sd = 10)

我想知道这些组之间的中位数是否存在显着差异。我知道我可以使用 Wilcoxon 测试来测试第 1 组与第 2 组，就像这样。

wilcox.test(group1, group2)

但是，这一次只比较两组，我想同时比较所有三组。我想要一个在 0.05 显着性水平上产生 ap 值的统计检验。有人可以帮忙吗？

编辑 #1 - Mood 的中值测试

按照用户休眠的建议答案，我尝试了 Mood 的中值测试。

median.test <- function(x, y){
    z <- c(x, y)
    g <- rep(1:2, c(length(x), length(y)))
    m <- median(z)
    fisher.test(z < m, g)$p.value
}

median.test(group1, group2)

但是，这种方法允许我一次测试仅两组的中位数之间的显着差异。我不确定如何使用它同时比较所有三个的中位数。

编辑 #2 - Kruskal-Wallis 测试

用户 dmartin 的建议答案似乎或多或少是我需要的，并允许我同时测试所有三个组。

kruskal.test(list(group1, group2, group3))

编辑#3

用户 Greg Snow 在他的回答中很有帮助地指出，Kruskal-Wallis 检验是合适的，只要它做出严格的假设，使其也是一种手段检验。

4个回答

首先，Wilcoxon 检验（或 Mann-Whitney 检验）不是中位数检验（除非您做出非常严格的假设，使其成为均值检验）。对于超过 2 组的比较，Wilcoxon 检验可能会导致一些自相矛盾的结果（参见Efron 的骰子）。

由于 Wilcoxon 检验只是置换检验的一个特例，而且您对中位数特别感兴趣，因此我建议对中位数进行置换检验。

首先选择差异的度量，例如 3 个中位数中的最大值减去 3 个中位数中的最小值（或 3 个中位数的方差，或 MAD 等）。

现在计算原始数据的统计数据。

将所有数据集中在一组中，然后将值随机分成 3 组

与原始大小相同并计算相同的统计数据。

重复多次（如 9998）

比较来自真实数据的统计信息与测试的所有统计信息的分布情况。

Mood 中位数检验是一种非参数检验，用于检验两个或多个总体的中位数是否相等。请参阅此处了解问题的 R 部分。另请参阅此处的相关问题。也从这里：

Mood 的中位数测试是最容易手动完成的测试：计算出（所有数据的）总体中位数，并计算每组中有多少值高于和低于中位数。如果各组大致相同，则观察值应在每组的总体中位数上下约 50-50... 低于中位数和高于中位数的计数...形成一个双向表，其中然后使用卡方检验进行分析。Mood 的中位数测试很像推广到两组或更多组的符号测试。

编辑： 对于三个组，您可以考虑我链接到​​的 R 代码的简单概括：

median.test2 <- function(x, y, z) {
  a <- c(x, y, z)
  g <- rep(1:3, c(length(x), length(y), length(z)))
  m <- median(a)
  fisher.test(a < m, g)$p.value
}

也可以使用 Kruskal-Wallis 检验，因为它是非参数方差分析。此外，它通常被认为比Mood 的中值测试更强大。它可以在 R 中使用stats 包中的kruskal.test函数在 R 中实现。

为了响应您的编辑，解释 KW 类似于单向方差分析。显着的 p 值对应于拒绝所有三个均值相等的空值。您必须使用后续测试（再次，就像 ANOVA）来回答有关特定组的问题。这通常遵循您可能遇到的特定研究问题。仅通过查看模拟参数，如果您进行后续测试，所有三个组应该彼此之间存在显着差异（因为它们都是 1 SD，N = 100）。

我知道这已经很晚了，但我也找不到适合 Mood 中值测试的包，所以我自己在 R 中创建了一个似乎可以解决问题的函数。 

#Mood's median test for a data frame with one column containing data (d),
#and another containing a factor/grouping variable (f)

moods.median = function(d,f) {

    #make a new matrix data frame
    m = cbind(f,d)
    colnames(m) = c("group", "value")


    #get the names of the factors/groups
    facs = unique(f)

    #count the number of factors/groups
    factorN = length(unique(f))


    #Make a 2 by K table that will be saved to the global environment by using "<<-":
    #2 rows (number of values > overall median & number of values <= overall median)
    #K-many columns for each level of the factor
    MoodsMedianTable <<- matrix(NA, nrow = 2, ncol = factorN)

    rownames(MoodsMedianTable) <<- c("> overall median", "<= overall median")
    colnames(MoodsMedianTable) <<- c(facs[1:factorN])
    colnames(MoodsMedianTable) <<- paste("Factor: ",colnames(MoodsMedianTable))


    #get the overall median
    overallmedian = median(d)



    #put the following into the 2 by K table:
    for(j in 1:factorN){ #for each factor level

        g = facs[j] #assign a temporary "group name"


        #count the number of observations in the factor that are greater than
        #the overall median and save it to the table
        MoodsMedianTable[1,j] <<- sum(m[,2][ which(m[,1]==g)] > overallmedian)


        #count the number of observations in the factor that are less than
        # or equal to the overall median and save it to the table
        MoodsMedianTable[2,j] <<- sum(m[,2][ which(m[,1]==g)] <= overallmedian)

    }


    #percent of cells with expected values less than 5
    percLT5 = ((sum(chisq.test(MoodsMedianTable)$expected < 5)) /
        (length(chisq.test(MoodsMedianTable)$expected)))


    #if >20% of cells have expected values less than 5
    #then give chi-squared stat, df, and Fisher's exact p.value
    if (percLT5 > 0.2) {
        return(list(
            "Chi-squared" = chisq.test(MoodsMedianTable)$statistic,
            "df" = chisq.test(MoodsMedianTable)$parameter,
            "Fisher's exact p.value" = fisher.test(MoodsMedianTable)$p.value))

    }


    #if <= 20% of cells have expected values less than 5
    #then give chi-squared stat, df, and chi-squared p.value
    if (percLT5 <= 0.2) {
        return(list(
            "Chi-squared" = chisq.test(MoodsMedianTable)$statistic,
            "df" = chisq.test(MoodsMedianTable)$parameter,
            "Chi-squared p.value" = chisq.test(MoodsMedianTable)$p.value))

    }

}

对于 OP 的问题，您将首先运行它以创建一个新的数据框，以使用匹配的“组”变量保存三个组向量中的值。

require(reshape2)
df = cbind(group1, group2, group3)
df = melt(df)
colnames(df) = c("observation", "group", "value")

并运行 Mood 中位数测试的函数moods.median(df$value, df$group)

其它你可能感兴趣的问题
上一篇了解 QR 分解 下一篇吨t-测试与zz-测试？