Подтвердить что ты не робот

Есть ли функция R, которая применяет функцию к каждой паре столбцов?

Мне часто приходится применять функцию к каждой паре столбцов в матрице данных/матрице и возвращать результаты в матрицу. Теперь я всегда пишу цикл, чтобы сделать это. Например, чтобы создать матрицу, содержащую р-значения корреляций, пишу:

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 

Итак, мой вопрос:

  • Из-за простоты этих функций я предполагаю, что это уже где-то в R. Есть ли функция apply или plyr, которая делает это? Я искал его, но не смог его найти.
  • Если да, то это быстрее?
4b9b3361

Ответ 1

Это не будет быстрее, но вы можете использовать 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)

Ответ 2

Я не уверен, правильно ли это относится к вашей проблеме, но посмотрите на пакет 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"

Ответ 3

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 ...

Ответ 4

Вы можете использовать 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