R - 对多个群体的每个排列进行统计测试的最简洁方法

Question

我将三个群体存储为单独的向量。我需要对这三个群体中的每一对进行统计检验(wilcoxon，如果重要的话)。

我想将三个向量输入到某个代码块中，并作为输出获得一个包含 6 个 p 值的向量(一个 p 值是一次测试的结果，是一个 double 值)。

我有一个可行的方法，但我是 R 的新手，根据我一直在阅读的内容，我觉得应该有更好的方法，可能涉及将向量存储为数据框并使用向量化来编写此代码.

这是我的代码:

library(arrangements)

runAllTests <- function(pop1,pop2,pop3) {
    populations <- list(pop1=pop1,pop2=pop2,pop3=pop3)
    colLabels <- c("pop1", "pop2", "pop3")

    #This line makes a data frame where each column is a pair of labels
    perms <- data.frame(t(permutations(colLabels,2)))

    pvals <- vector()

    #This for loop gets each column of that data frame
    for (pair in perms[,]) {
        pair <- as.vector(pair)
        p1 <- as.numeric(unlist(populations[pair[1]]))
        p2 <- as.numeric(unlist(populations[pair[2]]))

        pvals <- append(pvals, wilcox.test(p1, p2,alternative=c("less"))$p.value)
    }

    return(pvals)
}

编写此代码的 R 更合适的方法是什么？

注意:生成总体并将它们相互比较是一件很常见的事情(并且编码起来很棘手)，我认为这个问题将适用于比我自己更多的人。

编辑:我忘记了我的实际人口规模不同。这意味着我不能用向量制作数据框(据我所知)。不过我可以列出向量。我已经用一个有效的版本更新了我的代码。

Answer 1

是的，这确实很常见；确实如此普遍，以至于 R 有一个内置函数正好适用于这种情况:pairwise.table。

p <- list(pop1, pop2, pop3)
pairwise.table(function(i, j) { 
   wilcox.test(p[[i]], p[[j]])$p.value 
}, 1:3)

t检验、比例检验、Wilcoxon检验也有特定版本；这是一个使用 pairwise.wilcox.test 的示例。

p <- list(pop1, pop2, pop3)
d <- data.frame(x=unlist(p), g=rep(seq_along(p), sapply(p, length)))
with(d, pairwise.wilcox.test(x, g))

此外，请确保查看 p.adjust.method 参数以正确调整多重比较。

根据您的评论，您对顺序很重要的测试感兴趣；这真的很难想象(并且对于您提到的 Wilcoxon 测试而言并非如此)但仍然...

这是 pairwise.table 函数，经过编辑以在两个方向上进行测试。

pairwise.table.all <- function (compare.levels, level.names, p.adjust.method) {
  ix <- setNames(seq_along(level.names), level.names)
  pp <- outer(ix, ix, function(ivec, jvec) 
    sapply(seq_along(ivec), function(k) {
             i <- ivec[k]; j <- jvec[k]
             if (i != j) compare.levels(i, j) else NA }))
  pp[] <- p.adjust(pp[], p.adjust.method)
  pp
}

这是 pairwise.wilcox.test 的一个版本，它使用上述函数，也运行在向量列表上，而不是长格式的数据帧上。

pairwise.lazerbeam.test <- function(dat, p.adjust.method=p.adjust.methods) {
  p.adjust.method <- match.arg(p.adjust.method)
  level.names <- if(!is.null(names(dat))) names(dat) else seq_along(dat)
  PVAL <- pairwise.table.all(function(i, j) { 
    wilcox.test(dat[[i]], dat[[j]])$p.value 
  }, level.names, p.adjust.method = p.adjust.method)
  ans <- list(method = "Lazerbeam's special method", 
              data.name = paste(level.names, collapse=", "), 
              p.value = PVAL, p.adjust.method = p.adjust.method)
  class(ans) <- "pairwise.htest"
  ans
}

整理前后的输出如下所示:

> p <- list(a=1:5, b=2:8, c=10:16)
> out <- pairwise.lazerbeam.test(p)

> out

    Pairwise comparisons using Lazerbeams special method 

data:  a, b, c 

  a      b      c     
a -      0.2821 0.0101
b 0.2821 -      0.0035
c 0.0101 0.0035 -     

P value adjustment method: holm 

> pairwise.lazerbeam.test(p) %>% broom::tidy()
# A tibble: 6 x 3
  group1 group2 p.value
  <chr>  <chr>    <dbl>
1 b      a      0.282  
2 c      a      0.0101 
3 a      b      0.282  
4 c      b      0.00350
5 a      c      0.0101 
6 b      c      0.00350