Есть ли функция R, которая применяет функцию к каждой паре столбцов? - PullRequest
24 голосов
/ 08 марта 2011

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

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 и C ++ все время берет свое):

> 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 

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

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

Ответы [ 4 ]

16 голосов
/ 08 марта 2011

Это не будет быстрее, но вы можете использовать outer для упрощения кода.Для этого требуется векторизованная функция, поэтому здесь я использовал Vectorize для создания векторизованной версии функции для получения корреляции между двумя столбцами.

6 голосов
/ 08 марта 2011

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 ...
6 голосов
/ 08 марта 2011

Я не уверен, что это правильно решит вашу проблему, но взгляните на пакет psych Уильяма Ревелла.corr.test возвращает список матриц с коэффициентами корреляции, числом наблюдений, статистикой t-критерия и значением p.Я знаю, что использую это все время (и 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"
2 голосов
/ 09 марта 2011

Вы можете использовать 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
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...