Мне часто приходится применять функцию к каждой паре столбцов в матрице данных/матрице и возвращать результаты в матрицу. Теперь я всегда пишу цикл, чтобы сделать это. Например, чтобы создать матрицу, содержащую р-значения корреляций, пишу:
df <- data.frame(x=rnorm(100),y=rnorm(100),z=rnorm(100))
n <- ncol(df)
foo <- matrix(0,n,n)
for ( i in 1:n)
{
for (j in i:n)
{
foo[i,j] <- cor.test(df[,i],df[,j])$p.value
}
}
foo[lower.tri(foo)] <- t(foo)[lower.tri(foo)]
foo
[,1] [,2] [,3]
[1,] 0.0000000 0.7215071 0.5651266
[2,] 0.7215071 0.0000000 0.9019746
[3,] 0.5651266 0.9019746 0.0000000
который работает, но довольно медленный для очень больших матриц. Я могу написать функцию для этого в R (не беспокоясь о времени резания пополам, допуская симметричный результат, как указано выше):
Papply <- function(x,fun)
{
n <- ncol(x)
foo <- matrix(0,n,n)
for ( i in 1:n)
{
for (j in 1:n)
{
foo[i,j] <- fun(x[,i],x[,j])
}
}
return(foo)
}
Или функция с Rcpp:
library("Rcpp")
library("inline")
src <-
'
NumericMatrix x(xR);
Function f(fun);
NumericMatrix y(x.ncol(),x.ncol());
for (int i = 0; i < x.ncol(); i++)
{
for (int j = 0; j < x.ncol(); j++)
{
y(i,j) = as<double>(f(wrap(x(_,i)),wrap(x(_,j))));
}
}
return wrap(y);
'
Papply2 <- cxxfunction(signature(xR="numeric",fun="function"),src,plugin="Rcpp")
Но оба довольно медленные даже на довольно небольшом наборе данных из 100 переменных (я думал, что функция Rcpp будет быстрее, но я думаю, что преобразование между R и С++ все время имеет свои потери):
> system.time(Papply(matrix(rnorm(100*300),300,100),function(x,y)cor.test(x,y)$p.value))
user system elapsed
3.73 0.00 3.73
> system.time(Papply2(matrix(rnorm(100*300),300,100),function(x,y)cor.test(x,y)$p.value))
user system elapsed
3.71 0.02 3.75
Итак, мой вопрос:
plyr, которая делает это? Я искал его, но не смог его найти.Это не будет быстрее, но вы можете использовать outer для упрощения кода. Это требует векторизованной функции, поэтому здесь я использовал Vectorize, чтобы сделать векторизованную версию функции для получения корреляции между двумя столбцами.
df <- data.frame(x=rnorm(100),y=rnorm(100),z=rnorm(100))
n <- ncol(df)
corpij <- function(i,j,data) {cor.test(data[,i],data[,j])$p.value}
corp <- Vectorize(corpij, vectorize.args=list("i","j"))
outer(1:n,1:n,corp,data=df)
Я не уверен, правильно ли это относится к вашей проблеме, но посмотрите на пакет William Revelle psych. corr.test возвращает список матриц с коэффициентами корреляции, # из obs, t-test statistic и p-value. Я знаю, что я использую его все время (и AFAICS вы тоже психолог, поэтому он также может удовлетворить ваши потребности). Пишущие петли - не самый элегантный способ сделать это.
library(psych)
corr.test(mtcars)
( k <- corr.test(mtcars[1:5]) )
Call:corr.test(x = mtcars[1:5])
Correlation matrix
mpg cyl disp hp drat
mpg 1.00 -0.85 -0.85 -0.78 0.68
cyl -0.85 1.00 0.90 0.83 -0.70
disp -0.85 0.90 1.00 0.79 -0.71
hp -0.78 0.83 0.79 1.00 -0.45
drat 0.68 -0.70 -0.71 -0.45 1.00
Sample Size
mpg cyl disp hp drat
mpg 32 32 32 32 32
cyl 32 32 32 32 32
disp 32 32 32 32 32
hp 32 32 32 32 32
drat 32 32 32 32 32
Probability value
mpg cyl disp hp drat
mpg 0 0 0 0.00 0.00
cyl 0 0 0 0.00 0.00
disp 0 0 0 0.00 0.00
hp 0 0 0 0.00 0.01
drat 0 0 0 0.01 0.00
str(k)
List of 5
$ r : num [1:5, 1:5] 1 -0.852 -0.848 -0.776 0.681 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
$ n : num [1:5, 1:5] 32 32 32 32 32 32 32 32 32 32 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
$ t : num [1:5, 1:5] Inf -8.92 -8.75 -6.74 5.1 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
$ p : num [1:5, 1:5] 0.00 6.11e-10 9.38e-10 1.79e-07 1.78e-05 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
.. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
$ Call: language corr.test(x = mtcars[1:5])
- attr(*, "class")= chr [1:2] "psych" "corr.test"
92% времени тратится на cor.test.default и подпрограммы, которые он вызывает, поэтому его безнадежно пытаются получить более быстрые результаты, просто переписывая Papply (кроме экономии от вычисления только тех, что выше или ниже диагонали, предполагая, что ваш функция симметрична в x и y).
> M <- matrix(rnorm(100*300),300,100)
> Rprof(); junk <- Papply(M,function(x,y) cor.test( x, y)$p.value); Rprof(NULL)
> summaryRprof()
$by.self
self.time self.pct total.time total.pct
cor.test.default 4.36 29.54 13.56 91.87
# ... snip ...
Вы можете использовать mapply, но поскольку другие ответы говорят, что он вряд ли будет намного быстрее, поскольку большая часть времени используется cor.test.
matrix(mapply(function(x,y) cor.test(df[,x],df[,y])$p.value,rep(1:3,3),sort(rep(1:3,3))),nrow=3,ncol=3)
Вы могли бы уменьшить объем работы mapply, используя предположение о симметрии и отметив нулевую диагональ, например
v <- mapply(function(x,y) cor.test(df[,x],df[,y])$p.value,rep(1:2,2:1),rev(rep(3:2,2:1)))
m <- matrix(0,nrow=3,ncol=3)
m[lower.tri(m)] <- v
m[upper.tri(m)] <- v