Cor.test() принимает в качестве аргументов векторы x и y, но у меня есть целая матрица данных, которую я хочу проверить, попарно. Cor() воспринимает эту матрицу как аргумент просто отлично, и я надеюсь найти способ сделать то же самое для Cor.test().
Общим советом других людей, по-видимому, является использование cor.prob():
https://stat.ethz.ch/pipermail/r-help/2001-November/016201.html
Но эти значения p не совпадают с теми, которые генерируются Cor.test()!!! Cor.test() также лучше подходит для обработки парного удаления (у меня довольно много отсутствующих данных в моем наборе данных), чем cor.prob().
Есть ли у кого-нибудь альтернативы cor.prob()? Если решение включает вложенные для циклов, пусть будет так (я уже достаточно для R, даже если это будет проблематично для меня).
corr.test в пакете psych предназначен для этого:
library("psych")
data(sat.act)
corr.test(sat.act)
Как отмечено в комментариях, чтобы реплицировать значения p из базовой функции cor.test() по всей матрице, вам необходимо отключить настройку p-значений для нескольких сравнений (по умолчанию используется метод Холма регулировки):
corr.test(sat.act, adjust = "none")
[Но будьте осторожны при интерпретации этих результатов!]
Если вы строго следуете за pvalues в матричном формате от cor.test, здесь решение бесстыдно украдено у Vincent (LINK):
cor.test.p <- function(x){
FUN <- function(x, y) cor.test(x, y)[["p.value"]]
z <- outer(
colnames(x),
colnames(x),
Vectorize(function(i,j) FUN(x[,i], x[,j]))
)
dimnames(z) <- list(colnames(x), colnames(x))
z
}
cor.test.p(mtcars)
Примечание: Tommy также обеспечивает более быстрое решение, хотя и менее легкое для внедрения. Ох и нет для циклов:)
Изменить У меня есть функция v_outer в моем пакете qdapTools, которая делает эту задачу довольно простой:
library(qdapTools)
(out <- v_outer(mtcars, function(x, y) cor.test(x, y)[["p.value"]]))
print(out, digits=4) # for more digits
Вероятно, самый простой способ - использовать rcorr() из Hmisc. Он будет принимать только матрицу, поэтому используйте rcorr(as.matrix(x)), если ваши данные находятся в data.frame. Он вернет вам список: 1) матрица r попарно, 2) матрица попарно n, 3) матрица значений p для r. Он автоматически игнорирует отсутствующие данные.
В идеале, функция такого типа должна также принимать данные. Также выводить доверительные интервалы в соответствии с Новая статистика.
Принятое решение (функция corr.test в пакете psych) работает, но для больших матриц очень медленно. Я работал с матрицей экспрессии генов (~ 20000 на ~ 1000), коррелировал с матрицей чувствительности к лекарственным средствам (~ 1000 на ~ 500), и мне пришлось остановить ее, потому что она велась навсегда.
Я взял код из пакета psych и использовал функцию cor() напрямую и получил гораздо лучшие результаты:
# find (pairwise complete) correlation matrix between two matrices x and y
# compare to corr.test(x, y, adjust = "none")
n <- t(!is.na(x)) %*% (!is.na(y)) # same as count.pairwise(x,y) from psych package
r <- cor(x, y, use = "pairwise.complete.obs") # MUCH MUCH faster than corr.test()
cor2pvalue = function(r, n) {
t <- (r*sqrt(n-2))/sqrt(1-r^2)
p <- 2*(1 - pt(abs(t),(n-2)))
se <- sqrt((1-r*r)/(n-2))
out <- list(r, n, t, p, se)
names(out) <- c("r", "n", "t", "p", "se")
return(out)
}
# get a list with matrices of correlation, pvalues, standard error, etc.
result = cor2pvalue(r,n)
Даже с двумя матрицами размером 100 x 200 разница была ошеломляющей. Второе или два против 45 секунд.
> system.time(test_func(x,y))
user system elapsed
0.308 2.452 0.130
> system.time(corr.test(x,y, adjust = "none"))
user system elapsed
45.004 3.276 45.814