#目的
Rの標準で入っているライブラリlibrary(stats)の関数cor.testを使い、変数間の相関行列を求めます。
Rでは、変数のデータフレームを与えて一気に変数間の相関行列を出力させることが可能です。ですが変数が膨大にある場合、相関行列が出力されてもその読みとりは困難です。そこで、相関が強いと判定されたペアのみデータとして出力するスクリプトを作りました。
##想定
「相関が強い」の閾値は、絶対値で0.4としています。
母相関係数の無相関の検定の有意水準は5%です。
##アウトプットイメージ
こんな感じになります(結果はダミーです)。
i列とj列:変数の列番号
cor: 相関係数
p_val : p値
##データ準備
-
全ての列が量的な変数のデータ。
(factorの変数があれば、あらかじめ除去します) -
変数ラベルリスト
データとは別に変数ラベルリストを用意。
下記のvar_no は連番
ファイル名:labelfile.txt
##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)