R 大量の変数間の相関行列を求め、有意で相関の高いペアのみ抽出する

目的

Rの標準で入っているライブラリlibrary(stats)の関数cor.testを使い、変数間の相関行列を求めます。
Rでは、変数のデータフレームを与えて一気に変数間の相関行列を出力させることが可能です。ですが変数が膨大にある場合、相関行列が出力されてもその読みとりは困難です。そこで、相関が強いと判定されたペアのみデータとして出力するスクリプトを作りました。

想定

「相関が強い」の閾値は、絶対値で0.4としています。
母相関係数の無相関の検定の有意水準は5％です。

アウトプットイメージ

こんな感じになります（結果はダミーです）。
i列とj列：変数の列番号
cor: 相関係数
p_val : p値

データ準備

1. 全ての列が量的な変数のデータ。
   (factorの変数があれば、あらかじめ除去します）

2. 変数ラベルリスト
   データとは別に変数ラベルリストを用意。
   下記のvar_no は連番
   ファイル名:labelfile.txt
   

Rのサンプルコード

corr_pair.R

#-----------------------------------------------------
# cor.testの結果、相関係数絶対値0.4以上のペアのみ出力する
#-----------------------------------------------------
#相関行列を計算するデータを読み込む。
df<-read.table("testdata.dat",header=T,sep="\t")

#変数ラベルを読み込む
var_list<-read.table("labelfile.txt",header=T)

#相関係数を計算する。
#以下の場合のみ、変数ペアの列番号と相関係数、p値を書き出す。
#母集団相関係数の無相関の検定（両側）を行い、「帰無仮説：相関係数は0」が有意水準5%で棄却できるかどうかをp値より判定する。
#相関係数の絶対値0.4以上

#結果を格納するmatrixを準備。
#nn 変数がn個の場合、n×n/2 が計算する最大のペア数。0.4以上しか出力しないので、そこから適当に減らした値を入れています。

nn<-3000
output<-array(0,dim=c(nn,4))
k<-0

#n 準備したデータdfの変数数
n<-ncol(df)

#n×nの組み合わせ全ては計算しない。自分自身の相関も計算しない。
#出力するのは相関係数とp値のみ
#
for(i in 1:(n-1)){
  for(j in (i+1):n){
      res<-cor.test(df[,i],df[,j])
　　　　# estimateオブジェクト　相関係数の値
     　if(abs(res$estimate) >=0.4 ){
      　　#p.valueオブジェクト　p値　0.05以下
      　　　　if(res$p.value <=0.05){
        　　　　k<-k+1
        　　　　output[k,1]<-i　
        　　　　output[k,2]<-j  
        　　　　output[k,3]<-res$estimate
        　　　　output[k,4]<-res$p.value
        　　　}　
      　}
  }
  colnames(output)<-c("i","j","cor","p_val")
}

#相関の高いペアは、結局kこ
#出力用に余計な行を消す
output.2<-output[1:k,]

#相関係数マトリクスと変数ラベルデータフレームを結合する。
#結合キーは、相関係数データでは i 、変数ラベルデータでは、var_no
#相関係数データのiに合わせ、ラベルを結合
soukan<merge(output.2,var_list,by.x="i",by.y="var_no",all.x=TRUE)
head(soukan)
#データの列順を変更　ラベルを変数連番の後に移動
soukan<-soukan[,c("i","label","j","cor","p_val")]

str(soukan)
#変数名(label)をlabel_iにリネームしておく。
colnames(soukan)<-c("i","label_i","j","cor","p_val")

#相関係数抽出データの変数jにラベルを結合する
soukan2<-merge(soukan,var_list,by.x="j",by.y="var_no",all.x=TRUE)

head(soukan2)
#列順を変更する
soukan2<-soukan2[,c("i","lalbel_i","j","label","cor","p_val")]

#列名を変更する
colnames(soukan2)<-c("i","label_i","j","label_j","cor","p_val")

#結果を出力する
#macのExcelでオープンすることを想定。
write.table(soukan2,"相関高いペア抽出.txt",sep="\t",fileEncoding = "CP932",row.names = F)

実行環境

R version 3.2.2 (2015-08-14)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.3 (El Capitan)