何度もやってるのにどうも覚えられないのでここにまとめてメモ書き。
自分の研究は、質問紙などの順序尺度データだったり、正規性のないデータだったりを使うことが多いので、相関分析をする場合には必然的にピアソンの積率相関係数ではなく、順位相関係数を使うことが多い。
基本的にはスピアマンを使っているが、ケンドールが報告されているのをあまり見かけないのでいつかドヤ顔で使ってみたいという厨二心もある。
順位相関係数を出すだけならすぐにできる。Rだとこう(ダジャレではない)。
> cor(dat, method=”spearman”) #スピアマン
> cor(dat, method=”kendall”) #ケンドール
無相関検定とか信頼区間の算出をする場合はcor.test関数を使う(引数が2つ必要なことに注意)。
> cor.test(dat[,1],dat[,2],method=”spearman”)
Spearman’s rank correlation rho
data: dat[, 1] and dat[, 2]
S = 459.55, p-value = 0.001742
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.6544704警告メッセージ:
cor.test.default(dat[, 1], dat[, 2], method = “spearman”) で:
タイのため正確な p 値を計算することができ
はいでてきt、あ、あれ?
調べてみるとcor.test関数では、ピアソンの積率相関の場合のみにしか信頼区間を出してくれないらしい。
http://www.inside-r.org/r-doc/stats/cor.test
次はこれを試してみた。corrgramパッケージのcorrgram関数。信頼区間が一気出しできるpanel.confオプションが便利。
https://cran.r-project.org/web/packages/corrgram/corrgram.pdf
だがここにも落とし穴。
母相関係数の信頼区間を一気出しするためにcorrgramを使ったら、出てきた相関係数(スピアマンの順位相関)が、psychパッケージのpairs.panels(method=”spearman”)で出した相関係数と一致せず、おかしいなあと思って調べてみた。
— Kawaguchi Yusaku (@kwsk3939) 2015年7月31日
マニュアルを読みなおしてみると、どうもcorrgramで信頼区間を出すpanel.confオプションは、オプションで順位相関を選んだ場合(cor.method=”spearmman/kendall”)機能しないそうだ。 https://t.co/wjfBsh6rSG
— Kawaguchi Yusaku (@kwsk3939) 2015年7月31日
結局は中でcor.testが動いているので、順位相関では信頼区間が出せないいうことだそうだ。
はてさて、困った。
【解決策1】
langtestを使う。スピアマンでもケンドールでも使える。一番簡単かつ手っ取り早い。
http://langtest.jp/shiny/cor/
【解決策2】
スピアマンの場合、DescToolsパッケージのSpearmanRho関数を使う。
https://cran.r-project.org/web/packages/DescTools/DescTools.pdf
1,2番目の引数には変数、3番めの引数では
> library(DescTools)
> SpearmanRho(dat[,1],dat[,2],conf.level=0.95)
rho lwr.ci ups.ci
0.6544704 0.2983596 0.8506335
この結果はlangtestの結果と一致する。
【解決策3】
これもスピアマンの場合。RVAideMemoireパッケージのspearman.ci関数を使う。
https://cran.r-project.org/web/packages/RVAideMemoire/RVAideMemoire.pdf
これは上の2つとは違い、スピアマン順位相関のブートストラップ信頼区間を出してくれる。
> library(RVAideMemoire)
*** Package RVAideMemoire v 0.9-52 ***> spearman.ci(dat[,1], dat[,2], nrep=1000, conf.level=0.95)
Spearman’s rank correlation
data: dat[, 1] and dat[, 2]
1000 replicates95 percent confidence interval:
0.2263722 0.9238737
sample estimates:
rho
0.6544704
【解決策4】
これはケンドールの場合のみ。NSM3パッケージのkendall.ci関数を使う。
https://cran.r-project.org/web/packages/NSM3/NSM3.pdf
どうもこの本の計算方法に基づいて信頼区間を算出するようだ。
> kendall.ci(dat[,1], dat[,2], alpha=0.05, bootstrap=F) #漸近解析を使う場合
1 – alpha = 0.95 two-sided CI for tau:
0.187, 0.867> kendall.ci(dat[,1], dat[,2], alpha=0.05, bootstrap=T, B=1000) #ブートストラップを使う場合
1 – alpha = 0.95 two-sided CI for tau:
0.109, 0.838
ブートストラップを使用する場合は、Bの部分で抽出回数を指定する。
いろいろ種類があって、どれを使うのがいいかはよくわからないし、もう全部langtestでOKじゃなかろうかと思ったりもするが、手数は多いに越したことはないだろう。
【追記1】2016.03.13
このようなご指摘があったのでここで紹介させていただきます。
@kwsk3939 cor.test(rank(x), rank(y)) により,スピアマンの順位相関係数と信頼区間が求まる(表示はピアソンの積率相関係数だが)。
— foobar (@anonymous124816) 2016年3月13日
【追記2】2016.03.27
Nagoya.R #15にて、このブログの内容について発表しました。以下は発表スライドです。
【追記3】2016.03.27
スピアマンの順位相関係数の信頼区間について、以下のページをご紹介いただきました。
cor.test関数とrank関数の組み合わせで算出した信頼区間と、ブートストラップで算出した信頼区間との間には大きな差はないようだ。