Rでとことん相関

研究会でせっかく資料作ったのでこちらにもアップしたいと思います。内容は「相関」について。

心理学をやっている人ならば恐らくは誰でも知っている相関係数ですが、あらためて調べてみると色々出来ることは多いですね。
なお、統計初心者向けの内容となっております。

1.  相関（correlation）とは

相関とは、2 つの変数間の関係性を記述するための統計指標を指します。例えば，変数 A の値が大きいほど，変数 B が大きければ正の相関。変数 A の値が大きいほど，変数 B が小さければ負の相関です。ここあたりは基本ですね。

特に Pearson の積率相関係数では共分散を用いて算出されます。

2.  共分散（covariance）

共分散は2 変数の散布度についての統計指標です。それぞれの観測変数が 2 つの平均値からどの程度離れているのか（平均からの偏差）を元に算出されて、直線状に並んだときが最大となります。分散と非常に近い指標で、2変数に基づく分散と考えるとわかりやすいでしょう。

・分散：平均からの偏差の 2 乗を平均した値
・共分散：2 変数それぞれの平均からの偏差の積を平均した値

式にするともっとわかりやすいです。

$${ S }_{ XY }=\frac { 1 }{ n } \sum _{ i=1 }^{ n }{ { (X_{ i }-\bar { X } ) }{ (Y_{ i }-\bar { Y } ) } } $$

　この数式に合わせて以下のプロットを見てみましょう。

これはあるサンプルデータにおける年齢と身長のプロット図です。真ん中に走っている点線が年齢、身長それぞれの平均値となります。

この図からもわかる通り、共分散の「2 変数それぞれの平均からの偏差の積」は第1象限・第3象限では正の値、第2象限・第4象限では負の値となります。共分散はこの平均ですので、観測点が第1象限・第3象限に伸びていくほど値は正となり、反対に観測点が第2象限・第4象限に伸びていくほど値は負となります。また、あくまで偏差の積の平均ですので外れ値があればその分、正になったり負になったりするでしょう。

「プロットの形状が右上になってるほど正の相関」とか「外れ値があるとそれにひっぱられる」とか漠然とした理解の人も多いかとは思いますが、実際にはこのような数式に基づいて算出されているわけですね。

3.  基本的な相関分析：Pearson の積率相関係数の算出

共分散は2変数の平均値と観測値の差の積ですので、観測値の分散が大きいほど共分散も大きくなります。これでは異なる指標間で比較することが難しくなってしまいます。そこで，比較できるように各変数の標準偏差で割ったものである Pearson の積率相関係数が登場します。

式にすると以下の通り。

$$r=\frac { { S }_{ XY } }{ { S }_{ X }{ S }_{ Y } } =\frac { \frac { 1 }{ n } \sum _{ i=1 }^{ n }{ { ({ X }_{ i }-\overline { X } ) }{ ({ Y }_{ i }-\overline { Y } ) } }  }{ \sqrt { \frac { 1 }{ n } \sum _{ i=1 }^{ n }{ { ({ X }_{ i }-\overline { X } ) }^{ 2 } }  } \sqrt { \frac { 1 }{ n } \sum _{ i=1 }^{ n }{ { (Y_{ i }-\overline { Y } ) }^{ 2 } }  }  } $$

基本としてはbaseパッケージにあるcor関数を使って算出されます。

#2 変数だけでもできるし
cor(data$age,data1$height)

#複数選択も可能
cor(data)
cor(data[1:3])

なお、cor関数では応答変数はリスト形式で返されます。マトリックスで表示されますので、そのまま行列指定で別の関数の引数に入れることもできますね。

ただし、cor関数は無相関検定を行ってくれません。つまり有意な相関かどうかはわかりません。これでは困ります。
この問題を全部解決してくれるのがpsychパッケージ内のcorr.test関数です。

#psych パッケージの corr.test 関数なら無相関検定もできる
corr.test(data[1:3])

#有効桁数を増やしたければ print 関数
print(corr.test(data[1:3]), digits=3) 

#欠損処理はペアワイズなのでリストワイズにしたければuseオプション
corr.test(data[1:3], use="complete")

#corr.test関数では検定の回数等に応じて修正したp値も同時に算出
#デフォルトはholmの方法。違うのを使いたければadjustオプション
corr.test(data[1:3], adjust="bonferroni")

#信頼区間を出したければ print 関数の short オプションを
#組み合わせる
print(corr.test(data[1:3]), digits=3, short = FALSE)

これを使えばばっちりp値も出してくれますし、信頼区間も出せるので便利ですね。

4.  順位相関係数：KendallとSpearman

我々が扱う尺度は間隔尺度、比例尺度とは限りません。時には順序尺度同士の相関を検討したいときもあるかと思います。そこで登場するのがKendallの順位相関係数、Spearmanの順位相関係数です。

この2つの順位相関係数は両方とも順位相関係数ではありますが、算出の方法が異なります。簡単に言えば、Kendallでは対応する2つの順位の一致率から算出しており、Spearmanでは順位の差から検討しています。

数式で表現しますと、Kendallの順位相関係数τでは、

$$\tau =\frac { 2P }{ \frac { 1 }{ 2 } n(n-1) } -1$$

となります。ここでのPは「大小関係が一致する観測値の数」です。先ほどのプロット図で言うのならば、第1象限・第3象限に入っている観測値の数となります。すごくシンプルですね。

これに対してSpearmanｎ順位相関係数ρは、

$$\rho =1-\frac { 6\sum { { D }^{ 2 } }  }{ { N }^{ 3 }-N } $$

となります。ここでのDは「対応する2つの変数の順位の差」です。差を2乗したものの合計を組み合わせの数で調整しています。やはりシンプルですね。どちらを使うべきかと言われれば、どちらでも、といった感じです。

なお、これをRで算出するのは非常に簡単です。

#method オプションで"kendall","spearman"を
#指定すれば順位相関も算出可能
cor(data[1:3],method = "kendall")
corr.test(data[1:3], method = "kendall")

たったこれだけ。簡単ですね。

5.  偏相関係数の算出

相関係数を算出するときに厄介なのが偽相関。偽相関とは変数 A と変数 B が高い相関関係にあるとき，本来的には変数 A と変数 Y のみが相関関係にあるのにも関わらず，変数 B と変数 Y の相関係数も高い値を示す現象を指します。例えば、身長とIQは恐らく正の相関を示すでしょうが、これは実際には年齢とIQが相関しているだけで、身長が直接IQと何か関わりがあるわけではありません。

数式は少しややこしくて

$$変数Bと変数Yの偏相関係数=\frac { { r }_{ BY }-({ r }_{ AY }\times { r }_{ AB }) }{ \sqrt { 1-{ r }_{ AY }^{ 2 } } \sqrt { 1-{ r }_{ AB }^{ 2 } }  }$$

となります。式を見れば、関連する2変数の相関係数（A・YとA・B）を調整しているのが一目瞭然ですね。

Rではこれを計算してくれる関数が存在します。

#base パッケージの partial.r 関数で偏相関係数自体は出せる
#3 列目の変数を統制した 1・2 列目の偏相関係数なら…
partial.r(data[1:5],c(1,2),3)
#4,5 列目の変数を統制した 1・2 列目の偏相関係数なら…
partial.r(data[1:5],c(1,2),c(4,5)) 

上記の通り、偏相関係数の算出は非常に面倒です。偏相関係数についての関数については、青木先生のHPで公開されておりますので、ご確認ください。

6.  相関係数の差の検討

例えば、男性と女性とで身長と体重の相関係数が異なるとしましょう。男性の方が鍛えていたりいなかったりの分散が大きいので女性よりも相関が小さいかもしれません。

しかし、性別間で異なる相関係数が算出されたとしても，それは統計的に異なるとは言い切れません。r=.70とr=.10ぐらいなら明らかに差はありますが、.50と.40とかならわかりませんよね。どうせなら層別で相関係数が異なるのか統計的に検討したいところ。もちろん，回帰分析に交互作用項として投入すればいいじゃないかという話もありますが、あくまでここでは相関のお話です。

算出する際には，相関係数を Z 変換することで正規分布に近似させたものを使用します。対応のある 2 つの相関係数の差の検定では，比較対象となる 2 つの相関係数だけでなく，対象外の 4 つの相関係数も計算に含めて調整していきます。

数式で表現するなら

$$z=\frac { { z }_{ 1 }^{ \prime  }-{ z }_{ 2 }^{ \prime  } }{ \sqrt { \frac { 1 }{ { n }_{ 1 }-3 } +\frac { 1 }{ { n }_{ 2 }-3 }  }  } $$

となります。

Rではデータセットから直接、相関係数の差の検定をしてくれる関数は存在しません。代わりに、直接相関係数とサンプルサイズを入力することでp値を計算してくれるr.test関数が存在します。

#分析には psych パッケージの r.test 関数を使用
#変数 1～4 まで存在するとき，変数 1 と変数 2 の
#相関係数を r12 と表現するとする
#独立な 2 群における相関係数の差(r12 vs r34)
r.test(r12=.3, r34=.5, n=100, n2=150)

#対応のある 2 つの相関係数の差(r12 vs r34)
r.test(r12=.3, r34=.5, r13=.1, r14=.1, r23=.1, r24=.1, n=150)

#１つの変数を共有する，対応のある 2 つの相関係数の差(r12 vs r13)
r.test(r12=.3, r13=.5, r23=.1, n=150)

#いちいち相関係数を手入力するのは非常に面倒
#なので cor 関数から引っ張って来よう
#人数は nrow 関数で引っ張る
data.r <- cor(data[data$sex==1,])
data.r2 <- cor(data[data$sex==2,])
r.test(r12=data.r1[4,5], r34=data.r2[4,5], n= nrow(data[data$sex==1,]),
n2=nrow(data[data$sex==2,])) 

少し面倒は面倒ですが、こうすることで相関係数の差を統計的に検討できます。

7.  ポリコリック相関・ポリシリアル相関

我々はよく 5 件法あたりを使用するものの，そもそも連続変量として扱えるのは 7 件法以降であるという話もあります。つまり 3 件法，さらには5 件法ですら順序尺度となるんですね。ならば Kendall の順位相関を使用すればよいという話もありますが，Kendall だけでは汎用性に欠けます。そこで登場するのがポリコリック相関（多分相関係数; polycholic correlation coefficient），ポリシリアル相関係数（多分系列相関係数; polyserial correlation coefficient）です。

基本的な考え方として，「現在は順序尺度として表現されているものの，実際には連続変量である」という仮定が存在しています。その背後の連続変量が正規分布であるとするときに，どの程度の相関係数が得られるのが尤もらしいかという最尤推定法の形で算出されます。数式的には小杉考司先生の資料が詳しいです。
（※ただし，R における psych パッケージ内の polychoric 関数では最尤推定法以外を使用している模様。色々調べてはみましたが最終的に理解できず。どうやらタウ値を使っているようではあるのですが…）

なお，順序尺度同士の相関係数をポリコリック相関係数，順序尺度と連続変量との相関係数をポリシリアル相関係数と呼びます。これらを用いれば，例え，3 水準の尺度であろうとも因子分析を行うことができます。

では実際にどの程度しっかり推定できるのかシミュレーションしてみましょう。やり方としては清水裕士先生がHADで行ったシミュレーションをRでやってみようかと思います。

まずはサンプルデータを用意します。今回はMASSパッケージ内のmvrnorm関数を使ってデータ生成を行います。詳しい使い方は平川真さんのスライドがわかりやすいです。

#まず相関行列を定義
#今回は変数 a と変数 b の相関係数が.70 ちょうどになるように設定
x <- matrix(c(1.0, 0.7,
              0.7, 1.0),ncol=2)

#そしてMASSパッケージ内のmvrnorm関数を使ってデータ生成
#なおmvrnorm関数では標準偏差が1になるように乱数が発生します
#また、muオプションで平均値を定義
ex1 <- mvrnorm(n=1000, mu=c(10,10), Sigma=x, empirical=T)

#標準偏差1だと分けにくいので標準偏差が3になるように変換
a <- ex1[,1]
b <- ex1[,2]
ex1 <- data.frame(a,b)
ex1$a <- (ex1$a-mean(ex1$a))/sd(ex1$a)*3+10
ex1$b <- (ex1$b-mean(ex1$b))/sd(ex1$b)*3+10

#さらに変数 a,b に対して，～7，7～10，10～13，13～に
#それぞれ 1,2,3,4 を変数 ar,brに割り当てる
ex11 <- ex1[ex1$a<7,]
ex11$ar <- 1
ex12 <- ex1[ex1$a<10 & ex1$a>=7,]
ex12$ar <- 2
ex13 <- ex1[ex1$a<13 & ex1$a>=10,]
ex13$ar <- 3
ex14 <- ex1[ex1$a>=13,]
ex14$ar <- 4

ex1<-rbind(ex11,ex12,ex13,ex14)

ex11 <- ex1[ex1$b<7,]
ex11$br <- 1
ex12 <- ex1[ex1$b<10 & ex1$b>=7,]
ex12$br <- 2
ex13 <- ex1[ex1$b<13 & ex1$b>=10,]
ex13$br <- 3
ex14 <- ex1[ex1$b>=13,]
ex14$br <- 4

data.p <- rbind(ex11,ex12,ex13,ex14)

#まずはPearsonの積率相関係数、Kendallの順位相関係数を算出
cor(data.p)
cor(data.p[3:4], method="kendall")

なんとも美しくないスクリプトですね。適度に変更してください。

そしてここからポリコリック相関・ポリシリアル相関を算出します。今回はpsychパッケージ内にあるpolycholic関数・polyserial関数を使用します。

#ポリコリック相関では psych パッケージの polycholic 関数を使用
polychoric(data.p[3:4])

#ポリシリアル相関では psych パッケージの polyserial 関数を使用
#初めに連続変量，後に順序尺度
#このサンプルデータであれば1列目が連続変量、4列目が順序尺度なので
polyserial(data.p[1], data.p[4])

恐らく、ポリコリック相関やポリシリアル相関のほうが0.7、つまり真値に近い値となったのではないでしょうか？

なお、このpolycholic関数ではrhoにデータフレーム形式でポリコリック相関係数のマトリックスを返します。ですので、

poly <- polychoric(data.p[3:4])
poly$rho

とすることで因子分析にこの相関係数を用いることができます。いわゆるカテゴリカル因子分析と呼ばれるものです。

8. 最後に

いかがでしたでしょうか。あくまでここでは相関係数をRで扱っていくための初歩的な話を掲載しています。もっと詳しく知りたい方は随所随所に載せている資料等を参考にしてくださいませ。