![[PukiWiki]](https://xn--p8ja5bwe1i.jp:443/wiki/image/banner.png "[PukiWiki]")
はじめに †
ここでは、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)),]
> 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
論文に使用するグラフはカッコイイ方がいいので、グラフ作成にはggplot2を使います。
install.packages("ggplot2") library(ggplot2)
まず、がく片（萼片）の長さと幅をプロットしてみます。
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width))+ geom_point()+ theme(aspect.ratio=1)
共分散 †
[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関数を使います。
cov(iris$Sepal.Length, iris$Sepal.Width)
> 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)}{s(X) s(Y)} \]
標準偏差は正の値ですから、共分散が正なら相関係数も正、共分散が負なら相関係数も負になります。
相関係数の値は-1から1までの値をとります。 相関係数が1に近いほど強い正の相関があり、相関係数が-1に近いほど強い負の相関があります。 相関係数が0に近いときは相関がありません。
相関係数を求めるには、cor関数を使います。
cor(iris$Sepal.Length, iris$Sepal.Width)
> cor(iris$Sepal.Length, iris$Sepal.Width) [1] -0.1175698
相関係数の値によって、相関を次のように表現します。（ただし、この値の範囲や表現は厳格に決まっているわけではありません。）
- 0.9以上のとき「極めて強い正の相関がある」
- 0.7以上0.9未満のとき「強い正の相関がある」
- 0.4以上0.7未満のときは「（中程度の）正の相関がある」
- 0.2以上0.4未満のとき「弱い正の相関がある」
- -0.2以上0.2未満のとき「ほとんど相関がない」
- -0.4以上-0.2未満のとき「弱い負の相関がある」
- -0.7以上-0.4未満のとき「（中程度の）負の相関がある」
- -0.9以上-0.7未満のとき「強い負の相関がある」
- -0.9未満のとき「極めて強い負の相関がある」
がく片（萼片）の長さと幅の相関係数は-0.12なので、ほとんど相関はありません。
相関の検定 †
本当にがく片（萼片）の長さと幅に相関は無いのでしょうか。 相関係数の値がゼロに近かったのは偶然そうなったとは考えられないでしょうか。
そこで、「がく片（萼片）の長さと幅には相関がない」という仮説（これを帰無仮説といいます）を立て、これが成り立たないことを示します。 これを無相関検定といいます。 詳しくは統計の教科書で勉強してください。
具体的には、p値を求め、有意水準と比べます。 p値が有意水準より小さいとき、その水準において統計的に有意な相関があると考えます。 そうでないとき、その水準において有意な相関はないと考えます。 有意水準には、5%や1%がよく用いられます。
p値を求めるには、cor.test関数を使います。
cor.test(iris$Sepal.Length, iris$Sepal.Width)
> 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) を見ます。 p値が 0.15 なので、1%有意水準でも、5%有意水準でも、「ほとんど相関がない」というのは統計的に有意ではありません。
次に、相関係数の95%信頼区間（95 percent confidence interval）を見ます。 相関係数は、95%の確率で、約-0.27から約0.04の区間に入ります。 この区間に 0 が含まれているので、相関係数が 0 になる可能性もあります。 つまり、「相関がない（相関係数が 0 である）」という帰無仮説は棄却できません。
そこで、setosaのデータだけを取り出して、がく片（萼片）の長さと幅の相関を調べてみましょう。
setosa <-iris[1:50,] plot(setosa$Sepal.Length, setosa$Sepal.Width)
cor.test(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
相関係数が0.74、p値は6.7e-10、つまり0.00000000067だったので、setosaのがく片（萼片）の長さと幅の間には1%有意水準で（0.0000001%水準でも）統計的有意に強い正の相関があると言えます。
サンプルが少ないとp値が大きくなります。（サンプルが多いとp値が小さくなります。）
たとえば、setosaの1番目から5番目のデータだけを取り出してsetosa1とし、がく片（萼片）の長さと幅の相関を調べてみます。
setosa1 <- setosa[1:5,] ggplot(data=setosa1, aes(x=Sepal.Length, y=Sepal.Width))+ geom_point()+ theme(aspect.ratio=1)
cor.test(setosa1$Sepal.Length, setosa1$Sepal.Width)
> cor.test(setosa1$Sepal.Length, setosa1$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
すると、相関係数は0.68ですが、p値が0.21なので、有意な正の相関があるとは言えません。
同じように、setosaの21番目から25番目までのデータだけを取り出してsetosa2とし、相関を調べてみます。
setosa2 <- setosa[21:25,] ggplot(data=setosa2, aes(x=Sepal.Length, y=Sepal.Width))+ geom_point()+ theme(aspect.ratio=1)
cor.test(setosa2$Sepal.Length, setosa2$Sepal.Width)
> cor.test(setosa2$Sepal.Length, setosa2$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
setosaに含まれる50個のデータからは有意な正の相関が見られたのに、setosa2では相関係数が負になってしまいますが、p値が0.63ですから、これも統計的に有意な負の相関があるとは言えません。
このように、サンプルの取り方によっては相関係数が大きく変わってしまいますので、サンプルの相関係数を求めるだけでなく、無相関検定のp値を調べることによってその相関係数が統計的に有意かどうかを調べる必要があります。
演習 †
irisデータのPetal.LengthとPetal.Widthの間の相関係数を調べてみよう。
また、その統計的有意性を検定してみよう。
参考文献 †
- 樋口千洋: Rによるバイオインフォマティクスデータ解析 第2版－Bioconductorを用いたゲノムスケールのデータマイニング－, 共立出版 (2011)
cov.png
setosa2.png
setosa1.png
setosa.png
iris.png
