#list #r
#Список #r
Вопрос:
Я использую следующую функцию для вычисления t-stat для данных во фрейме данных (x):
wilcox.test.all.genes<-function(x,s1,s2) {
x1<-x[s1]
x2<-x[s2]
x1<-as.numeric(x1)
x2<-as.numeric(x2)
wilcox.out<-wilcox.test(x1,x2,exact=F,alternative="two.sided",correct=T)
out<-as.numeric(wilcox.out$statistic)
return(out)
}
Мне нужно написать цикл for, который будет повторяться определенное количество раз. Для каждой итерации столбцы необходимо перетасовывать, выполнять вышеуказанную функцию и сохранять максимальное значение t-stat в списке.
Я знаю, что могу использовать sample() функцию для перетасовки столбцов фрейма данных и max() функцию для определения максимального значения t-stat, но я не могу понять, как соединить их вместе, чтобы получить работоспособный код.
Комментарии:
1. Приведенный ниже ответ работает, но мне все еще немного неясно, что вы на самом деле пытаетесь сделать. Сравнение вашей реальной тестовой статистики с результатами распределения, полученными в результате этих перестановок, даст вам эмпирическое значение p. Однако вы упомянули о максимальном значении, которое обычно используется для исправления множественных сравнений, если вы выполняете много подобных тестов. Однако меня смущает то, что ваш x является вектором, и вы выполняете только один тест. В любом случае, если вы немного проясните, я буду рад соответствующим образом отредактировать свой ответ.
2. Фрейм данных, который я использую, содержит данные от двух классов индивидуумов, я разделяю данные на основе класса. Функция представляет собой тест Уилкокса на данных из двух классов. После завершения первоначального теста я хочу перетасовать столбцы фрейма данных, снова выполнить функцию, получить максимальную статистику теста из новых данных и сохранить это число в списке. Я хочу выполнить перетасовку 500 раз, в результате чего будет получен список из 500 максимальных тестовых статистических данных из перетасованных данных. Я буду использовать этот список тестовой статистики для определения 95% тестовой статистики для сравнения с исходным набором данных. Спасибо
3. Я понимаю. Это то, что я подумал. Было бы полезно, если бы вы предоставили примерные данные, поскольку у нас не было возможности узнать, что у вас есть фрейм данных. Кроме того, возникает небольшая путаница, потому что для выполнения подобного тестирования перестановок вы переставляете метки классов (т. Е. строки ), а не столбцы. Смотрите мой отредактированный ответ ниже для вашего решения. Удачи!
Ответ №1:
Вы пытаетесь сгенерировать эмпирические p-значения с поправкой на множественные сравнения, которые вы проводите из-за множества столбцов в ваших данных. Сначала давайте смоделируем пример набора данных:
# Simulate data
n.row = 100
n.col = 10
set.seed(12345)
group = factor(sample(2, n.row, replace=T))
data = data.frame(matrix(rnorm(n.row*n.col), nrow=n.row))
Вычислите критерий Уилкоксона для каждого столбца, но мы будем повторять это много раз, переставляя принадлежность наблюдений к классу. Это дает нам эмпирическое нулевое распределение этой тестовой статистики.
# Re-calculate columnwise test statisitics many times while permuting class labels
perms = replicate(500, apply(data[sample(nrow(data)), ], 2, function(x) wilcox.test(x[group==1], x[group==2], exact=F, alternative="two.sided", correct=T)$stat))
Вычислите нулевое распределение максимальной тестовой статистики путем свертывания по нескольким сравнениям.
# For each permuted replication, calculate the max test statistic across the multiple comparisons
perms.max = apply(perms, 2, max)
Теперь, просто отсортировав результаты, мы можем определить критическое значение p = 0.05.
# Identify critical value
crit = sort(perms.max)[round((1-0.05)*length(perms.max))]
Мы также можем построить наше распределение вместе с критическим значением.
# Plot
dev.new(width=4, height=4)
hist(perms.max)
abline(v=crit, col='red')
Наконец, сравнение реальной тестовой статистики с этим распределением даст вам эмпирическое значение p, скорректированное для множественных сравнений путем регулирования ошибки по семейству до p<0,05. Например, давайте представим, что реальная тестовая статистика была 1600. Затем мы могли бы вычислить p-значение следующим образом:
> length(which(perms.max>1600))/length(perms.max)
[1] 0.074