Rで相関係数を求める方法とそのまま使える自作関数例

R言語入門

ライター：高橋光太郎

今回は相関係数の求め方についてです。R言語で相関係数は”cor()”という関数を使って行います。また、相関係数に関する検定はcor.test()で簡単に行えます。相関係数

この記事では、R言語を用いた実際の相関係数の算出方法について記述していますが、相関係数の意味と数式上の定義は相関係数の意味と定義（公式）でまとめてあります。

※R言語入門のトップページはこちら

1 単に相関係数を求めるだけのプログラム
- 1.1 ピアソンの積率相関係数
- 1.2 スピアマンの順位相関係数
2 相関係数に関する検定を行う関数、cor.test()
- 2.1 使用例：ピアソンの積率相関係数の場合
- 2.2 使用例：スピアマンの順位相関係数の場合
3 ピアソン、スピアマンの相関係数に関する検定を行い、フィアルへの出力を行う自作関数
- 3.1 使用例

単に相関係数を求めるだけのプログラム

説明のための、ベクトルデータを３種類用意して、それをそれぞれ変数a,b,cに代入しておきます。

プログラム

a <- c(1:10)
b <- c(1,5,23,6,4,2,9,8,10,6)
c <- c(2,4,6,8,10,13,8,2,19,29)

さらにこれらをひとつのdatという変数にまとめて格納ひとつのデータとして格納します。

プログラム

dat <- cbind(a,b,c)

ピアソンの積率相関係数

a,b,c、互いの、ピアソンの相関係数は関数”cor()”を用いて以下のように求められます。

プログラム

> cor(dat)
a b c
a 1.0000000 0.0355953 0.7208661
b 0.0355953 1.0000000 -0.0456211
c 0.7208661 -0.0456211 1.0000000

また、aとbの相関係数だけならば

プログラム

cor(a,b)
[1] 0.0355953

とすることで求めることが出来ます。

スピアマンの順位相関係数

cor()で何も指定をしなければ、ピアソンの積率相関係数が結果として出力されます。スピアマンの相関係数を求める場合は、次のようにオプションで

method=”s”　もしくは　method=”spearman”

を指定します。

プログラム

#どちらでも可

cor(a,b,method=”s”)
cor(a,b,method=”spearman”)

相関係数に関する検定を行う関数、cor.test()

相関係数に関する検定を行う場合は、cor.test()を使います。使い方は、

cor.test(ベクトル1,ベクトル2)

として、ベクトル1とベクトル2の相関係数に関する検定を行います。

使用例：ピアソンの積率相関係数の場合

プログラム

cor.test(dat[,1],dat[,3])

Pearson’s product-moment correlation
data: dat[, 1] and dat[, 3]
t = 1.2949, df = 45, p-value = 0.202
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1032543 0.4520929
sample estimates:
cor
0.1895324

結果は、

p-value = 0.202　の部分がp値
t = 1.2949　の部分が統計量
最後の　cor 0.1895324　の部分がピアソンの相関係数

を示しています。

使用例：スピアマンの順位相関係数の場合

プログラム

cor.test(dat[,1],dat[,3],method = “s”)

Spearman’s rank correlation rho
data: dat[, 1] and dat[, 3]
S = 18428, p-value = 0.6621
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.06543193

結果の見方は、ピアソンの場合と同様です。

ピアソン、スピアマンの相関係数に関する検定を行い、フィアルへの出力を行う自作関数

相関係数を大量に算出したい場合に、そのまま使える自作関数を作成しました。解析用データは1列目に目的変数、その他の列を説明変数とします。

以下の関数は解析用データの1列目と、そのほかの列のピアソンの相関係数とスピアマンの相関係数を出力します。また、相関係数に関する仮説検定を行い、p値も出力するようにしてあります。

プログラム

cor_fun <-function(dat,filename) #引数は（解析用データ,出力先ファイル名)
{
for(i in 2:ncol(dat))
{
ans_c = cor.test(dat[,1],dat[,i])
c= ans_c$estimate #ピアソンの相関係数
c_p= ans_c$p.value #そのp値
ans_s = cor.test(dat[,1],dat[,i],method=“s”)
s= ans_s$estimate #スピアマンの相関係数
s_p= ans_s$p.value #そのp値
result = c(paste(colnames(dat)[1],“VS”,colnames(dat)[i]),c,c_p,s,s_p)
names(result)[1] = “Name”
names(result)[3] = “ピアソンの相関係数のp値”
names(result)[5] = “スピアマンの相関係数のp値”
if(i == 2)
{
write.table(matrix(c(colnames(t(result))),nrow=1),filename,append=T,quote=F,sep=“,”,row.names=F,col.names=F,fileEncoding=“CP932”)
}
write.table(t(result),filename,append=T,quote=F,sep=“,”,row.names=F,col.names=F,fileEncoding=“CP932”)
}
}