0

My question is on the permutation of correlation coefficients.

        A<-data.frame(A1=c(1,2,3,4,5),B1=c(6,7,8,9,10),C1=c(11,12,13,14,15 ))

        B<-data.frame(A2=c(6,7,7,10,11),B2=c(2,1,3,8,11),C2=c(1,5,16,7,8))

          cor(A,B)

          #           A2        B2       C2
          # A1 0.9481224 0.9190183 0.459588
          # B1 0.9481224 0.9190183 0.459588
          # C1 0.9481224 0.9190183 0.459588

I obtained this correlation and then wanted to perform permutation tests to check if the correlation still holds.

I did the permutation as follows:

              A<-as.vector(t(A))
              B<-as.vector(t(B))

     corperm <- function(A,B,1000) {
     # n is the number of permutations
     # x and y are the vectors to correlate
    obs = abs(cor(A,B))
    tmp = sapply(1:n,function(z) {abs(cor(sample(A,replace=TRUE),B))})
   return(1-sum(obs>tmp)/n)
     }

The result was 

         [1] 0.645

and using "cor.test"

cor.test(A,B)

Pearson's product-moment correlation

data:  A and B
t = 0.4753, df = 13, p-value = 0.6425
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4089539  0.6026075
sample estimates:
cor 
0.1306868

How could I draw a plot or a histogram to show the actual correlation and the permuted correlation value from the permuted data ???

4

2 に答える 2

2

まず第一に、このように正確に行うことはできません...

> corperm = function(A,B,1000) {
Error: unexpected numeric constant in "corperm = function(A,B,1000"

3 番目の引数には名前がありませんが、名前が必要です。おそらくあなたが意味した

> corperm <- function(A, B, n=1000) {
# etc

次に、何を達成したいのかを考える必要があります。最初に、3 つの変数を持つ 2 つのデータ セットがあり、次にそれらを 2 つのベクトルに折りたたんで、並べ替えられたベクトル間の相関を計算します。なぜそれが理にかなっているのですか？置換されたデータ セットの構造は、元のデータ セットと同じである必要があります。

obs = abs(cor(A,B))
tmp = sapply(1:n,function(z) {abs(cor(sample(A,replace=TRUE),B))})
return(1-sum(obs>tmp)/n)

ここで replace=TRUE を使用するのはなぜですか? これは、ブートストラップ CI が必要な場合に理にかなっていますが、(a) パッケージ ブートからのブートなど、専用の関数を使用する方がよいでしょう。また、(B) B でも同じことを行う必要があります。つまり、サンプル (B、置換 = TRUE)。

置換テストの場合、置換なしでサンプリングし、A と B の両方に対して行うか、A のみに対して行うかに違いはありません。

そして、ヒストグラムを取得する方法は？hist(tmp) は並べ替えられた値のヒストグラムを描画し、obs は観察された相関の絶対値です。

HTHAB

（編集）

corperm <- function(x, y, N=1000, plot=FALSE){
    reps <- replicate(N, cor(sample(x), y))
    obs <- cor(x,y)
    p <- mean(reps > obs) # shortcut for sum(reps > obs)/N
    if(plot){
        hist(reps)
        abline(v=obs, col="red")
        }
     p
     }

これで、変数の単一のペアでこれを使用できます。

corperm(A[,1], B[,1])

すべてのペアに適用するには、forまたはを使用しますmapply。 の方が理解しやすいので、可能なすべてのペアを取得するためforに使用することを主張しません。mapply

res <- matrix(NA, nrow=NCOL(A), ncol=NCOL(B))
for(iii in 1:3) for(jjj in 1:3) res[iii,jjj] <- corperm(A[,iii], B[,jjj], plot=FALSE)
rownames(res)<-names(A)
colnames(res) <- names(B)
print(res)

すべてのヒストグラムを作成するには、上記の plot=TRUE を使用します。

于 2013-09-24T11:30:06.370 に答える
0

2 つのバリアントの相関分析に順列検定を行う意味はあまりないと思います。このcor.test()関数は、順列検定と同じ効果を持つ "p.value" を提供するためです。

于 2015-12-26T13:13:12.087 に答える