Rで相関分析する

2016-11-16 (水) 14:24:21 (165d) | Topic path: Top / バイオ・データ・マイニング / Rで相関分析する

はじめに

ここでは、Rを使って相関係数を求め、統計的有意性を検定する方法を説明します。

『Rによるバイオインフォマティクスデータ解析』の3.3節「基本統計関数」に少しだけ出てきます。

準備

Rのインストールについては、次のページを見てください。

ここでは、標準で使用できるirisデータセットを使います。

> data(iris)

このデータセットは、アヤメの種類(Species)を花びらの長さ(Sepal.Length)、幅(Sepal.Width)、がくの長さ(Petal.Length)、幅(Petal.Width)によって分類する問題です。 長さと幅は連続値、種類はsetosa, versicolor, virginicaのいずれかをとる離散値です。

このデータセットには、setosa, versicolor, virginicaという3種類のアヤメについて、それぞれ50個ずつ、合計150個のデータが含まれています。 ランダムに10個のデータを選択して、見てみましょう。

> iris[sort(sample(1:150,10)),]
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
4            4.6         3.1          1.5         0.2     setosa
22           5.1         3.7          1.5         0.4     setosa
65           5.6         2.9          3.6         1.3 versicolor
97           5.7         2.9          4.2         1.3 versicolor
100          5.7         2.8          4.1         1.3 versicolor
108          7.3         2.9          6.3         1.8  virginica
116          6.4         3.2          5.3         2.3  virginica
122          5.6         2.8          4.9         2.0  virginica
136          7.7         3.0          6.1         2.3  virginica
146          6.7         3.0          5.2         2.3  virginica

がく片(萼片)の長さと幅をプロットしてみます。

> plot(iris$Sepal.Length, iris$Sepal.Width)
iris.png

共分散

[math]X[/math] と [math]Y[/math] の共分散は、それぞれの平均値からのズレ(偏差)の積の和を [math]n - 1[/math] で割ったものです。 \[ \mathrm{Cov}(X,Y) = \frac{\sum_{i=1}^{n} (x_i - \bar{X})(y_i - \bar{Y})}{n-1} \] ここで、[math]n[/math] はデータ数、[math]\bar{X}[/math], [math]\bar{Y}[/math] は、それぞれ、[math]X[/math], [math]Y[/math] の平均です。

偏差の符号が同じとき、つまり、どちらも平均より大きいかどちらも平均より小さいとき、共分散は大きくなり、偏差の符号が異なるとき、つまり、いずれかが平均より大きくもう一方が平均より小さいとき、共分散は小さくなります。

cov.png

共分散を求めるには、cov関数を使います。

> cov(iris$Sepal.Length, iris$Sepal.Width)
[1] -0.042434

相関係数

[math]X[/math] と [math]Y[/math] の相関係数は、[math]X[/math] と [math]Y[/math] の共分散を [math]X[/math] の不偏標準偏差と [math]Y[/math] の不偏標準偏差の積で割ったものです。 \[ \mathrm{Cor}(X, Y) = \frac{\mathrm{Cov}(X, Y)}{\sigma(X) \sigma(Y)} \]

不偏標準偏差は正の値ですから、共分散が正なら相関係数も正、共分散が負なら相関係数も負になります。

相関係数の値は-1から1までの値をとります。 相関係数が1に近いほど強い正の相関があり、相関係数が-1に近いほど強い負の相関があります。 相関係数が0に近いときは相関がありません。

相関係数を求めるには、cor関数を使います。

> cor(iris$Sepal.Length, iris$Sepal.Width)
[1] -0.1175698

相関の検定

がく片(萼片)の長さと幅の相関係数は-0.12で、ほとんど相関がないことがわかりました。

本当に相関は無いのでしょうか。 相関係数の値がゼロに近かったのはこのサンプルだから偶然そうなったということは考えられないでしょうか。

そこで、「がく片(萼片)の長さと幅には相関がない」という仮説(これを帰無仮説といいます)を立て、これが成り立たないことを示します。 これを無相関検定といいます。 詳しくは統計の教科書で勉強してください。

具体的には、p値を求め、有意水準と比べます。 p値が有意水準より小さいとき、その水準において有意な相関があると考えます。 そうでないとき、その水準において有意な相関はないと考えます。 有意水準には、5%や1%がよく用いられます。

p値を求めるには、cor.test関数を使います。

> cor.test(iris$Sepal.Length, iris$Sepal.Width)

	Pearson's product-moment correlation

data:  iris$Sepal.Length and iris$Sepal.Width
t = -1.4403, df = 148, p-value = 0.1519
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.27269325  0.04351158
sample estimates:
       cor 
-0.1175698 

まず、p値 (p-value) を見ます。 0.15なので、1%有意水準でも、5%有意水準でも、相関があるとは言えません。

次に、相関係数の95%信頼区間(95 percent confidence interval)を見ます。 相関係数は、95%の確率で、約-0.27から約0.04の区間に入ります。 この区間に 0 が含まれているので、相関係数が 0 になる可能性もあります。 つまり、「相関がない(相関係数が 0 である)」という帰無仮説は棄却できません。

そこで、setosaのデータだけを取り出して、がく片(萼片)の長さと幅の相関を調べてみましょう。

> setosa <- iris[iris['Species']=='setosa',]
> plot(setosa$Sepal.Length, setosa$Sepal.Width)
> cor.test(setosa$Sepal.Length, setosa$Sepal.Width)

	Pearson's product-moment correlation

data:  setosa$Sepal.Length and setosa$Sepal.Width
t = 7.6807, df = 48, p-value = 6.71e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5851391 0.8460314
sample estimates:
      cor 
0.7425467 
setosa.png

相関係数が0.74、p値は6.7e-10、つまり0.00000000067だったので、setosaのがく片(萼片)の長さと幅の間には1%水準で(0.0000001%水準でも)有意にかなり強い正の相関があると言えます。

サンプル数が少ないと、p値が大きくなる傾向にあります。

たとえば、setosaの1番目から5番目のデータだけを取り出して、がく片(萼片)の長さと幅の相関を調べてみます。

> s1 <- setosa[1:5,]
> s1
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
> plot(s1$Sepal.Length, s1$Sepal.Width)
> cor.test(s1$Sepal.Length, s1$Sepal.Width)

	Pearson's product-moment correlation

data:  s1$Sepal.Length and s1$Sepal.Width
t = 1.6064, df = 3, p-value = 0.2065
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.5055650  0.9764541
sample estimates:
      cor 
0.6800193 
s1.png

すると、相関係数は0.68ですが、p値が0.21なので、有意な正の相関があるとは言えません。

また、21番目から25番目までのデータだけを取り出して、同じように相関を調べてみると、(setosaに含まれる50個のデータからは有意な正の相関が見られたのに)相関係数が-0.30と負になってしまいますが、p値が0.63ですから、これも有意な負の相関があるとは言えません。

> s2 <- setosa[21:25,]
> s2
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
21          5.4         3.4          1.7         0.2  setosa
22          5.1         3.7          1.5         0.4  setosa
23          4.6         3.6          1.0         0.2  setosa
24          5.1         3.3          1.7         0.5  setosa
25          4.8         3.4          1.9         0.2  setosa
> plot(s2$Sepal.Length, s2$Sepal.Width)
> cor.test(s2$Sepal.Length, s2$Sepal.Width)

	Pearson's product-moment correlation

data:  s2$Sepal.Length and s2$Sepal.Width
t = -0.5371, df = 3, p-value = 0.6285
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9343030  0.7934151
sample estimates:
       cor 
-0.2961744 
s2.png

このように、サンプルの取り方によっては相関係数が大きく変わってしまいますので、p値を調べることによって、その相関が統計的に有意な相関であるかどうかを調べる必要があります。

演習

irisデータのPetal.LengthとPetal.Widthの間の相関係数を調べてみよう。

また、その統計的有意性を検定してみよう。

参考文献

この本の3.3節「基本統計関数」に少しだけ出てきます。

添付ファイル: files2.png 70件 [詳細] files1.png 65件 [詳細] filesetosa.png 63件 [詳細] fileiris.png 64件 [詳細] filecov.png 58件 [詳細]
トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS