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

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

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(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. Если да, то это быстрее?

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

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

    размер блока данных в HDFS, почему 64 МБ?

    Как разрешить конфликты слияния в Git

    Приложение, использующее профиль bluetooth SPP, не работает после обновления с Android 4.2 до Android 4.3

    Переменные оболочки, находящиеся внутри цикла, не видны вне его

    Сколько виртуальной памяти?

    Есть ли способ заставить Outlook 2003 отправлять электронное письмо в определенное время, а не сразу?

    Мне нужен быстрый синтаксический анализатор времени выполнения

    Где находится предварительный просмотр макета Android Studio?

    Все JFrame замерзают, когда JavaMail

    Запуск активности из виджета

    RuntimeException: ваш контент должен иметь ListView, чей атрибут id – ‘android.R.id.list’

    Выполнять команды с использованием sudo на удаленном сервере после входа в PuTTY через командный файл

    Можно ли расшифровывать hashи md5?

    Максимальный размер метода в Java 7 и 8

    Как получить часть файла после строки, которая соответствует выражению grep? (первое совпадение)

    Давайте будем гением компьютера.