25

多くの場合、データフレーム/マトリックス内の列の各ペアに関数を適用し、結果をマトリックスで返す必要があります。今、私はいつもこれを行うためにループを書いています。たとえば、相関の 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 のどこかにあると思いますplyr。これを行う適用または関数はありますか? 私はそれを探しましたが、見つけることができませんでした。
  2. もしそうなら、それはより速いですか?
4

4 に答える 4

18

高速ではありませんがouter、コードを簡素化するために使用できます。ベクトル化された関数が必要なので、ここVectorizeでは関数のベクトル化バージョンを作成して 2 つの列間の相関関係を取得しました。

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)
于 2011-03-08T14:20:50.753 に答える
6

これがあなたの問題に適切な方法で対処するかどうかはわかりませんが、William Revelle のpsychパッケージを見てください。corr.test相関係数、obs の数、t 検定統計量、および p 値を含む行列のリストを返します。私はいつもそれを使用していることを知っています(そしてAFAICS、あなたは心理学者でもあるので、あなたのニーズにも合うかもしれません). ループを書くことは、これを行うための最も洗練された方法ではありません。

> library(psych)
> ( 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"
于 2011-03-08T14:04:17.500 に答える
6

時間の 92% は とそれが呼び出すルーチンに費やされcor.test.defaultているため、単純に書き直してより高速な結果を得ようとしても絶望的です(関数が と で対称であるとPapply仮定して、対角線の上または下のみを計算することによる節約を除いて)。 xy

> 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 ...
于 2011-03-08T14:09:25.017 に答える
2

を使用できます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
于 2011-03-09T11:24:10.540 に答える