R - Самый чистый способ запустить статистический тест на каждую перестановку нескольких групп населения

Question

У меня есть три популяции, хранящиеся как отдельные векторы.Мне нужно запустить статистический тест (Уилкоксон, если это имеет значение) для каждой пары из этих трех групп населения.

Я хочу ввести три вектора в некоторый блок кода и получить в качестве вывода вектор из 6 p-значений(одно значение p является результатом одного теста и является двойным).

У меня есть метод, который работает, но я новичок в 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-тестов, пропорциональных тестов и тестов Вилкоксона; Вот пример использования 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, чтобы правильно настроить множественные сравнения.

По вашим комментариям вас интересуют тесты, в которых важен порядок; это действительно трудно представить (и это не так для теста Уилкоксона, о котором вы упоминали), но все же ...

Это функция 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

Answer 2

Вот пример одного подхода, который использует combn(), у которого есть аргумент функции, который можно использовать для простого применения wilcox.test() ко всем комбинациям переменных.

set.seed(234)

# Create dummy data
df <- data.frame(replicate(3, sample(1:5, 100, replace = TRUE)))

# Apply wilcox.test to all combinations of variables in data frame.
res <- combn(names(df), 2, function(x) list(data = c(paste(x[1], x[2])), p = wilcox.test(x = df[[x[1]]], y =  df[[x[2]]])$p.value), simplify = FALSE)

# Bind results
do.call(rbind, res) 

     data    p         
[1,] "X1 X2" 0.45282   
[2,] "X1 X3" 0.06095539
[3,] "X2 X3" 0.3162251