順位相関係数の信頼区間の算出

何度もやってるのにどうも覚えられないのでここにまとめてメモ書き。

自分の研究は、質問紙などの順序尺度データだったり、正規性のないデータだったりを使うことが多いので、相関分析をする場合には必然的にピアソンの積率相関係数ではなく、順位相関係数を使うことが多い。

基本的にはスピアマンを使っているが、ケンドールが報告されているのをあまり見かけないのでいつかドヤ顔で使ってみたいという厨二心もある。

順位相関係数を出すだけならすぐにできる。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

だがここにも落とし穴。

結局は中で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 replicates

95 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

このようなご指摘があったのでここで紹介させていただきます。

 

【追記2】2016.03.27

Nagoya.R #15にて、このブログの内容について発表しました。以下は発表スライドです。

【追記3】2016.03.27

スピアマンの順位相関係数の信頼区間について、以下のページをご紹介いただきました。

スピアマンの順位相関係数の信頼区間 – 裏 RjpWiki

cor.test関数とrank関数の組み合わせで算出した信頼区間と、ブートストラップで算出した信頼区間との間には大きな差はないようだ。

Nagoya.R #11,終わりました

こんばんは。

先日Nagoya.R #11が開催されました。

今回私は入門者講習を担当させていただき,なんとか無事(?)終えることができました。

拙い講習を優しく見守ってくださったオーディエンスの皆様ありがとうございました。

 

SlideShareに入門者講習にて使用したスライドを一部改変したものをアップしていますのでよろしければご覧ください。

このスライドの作成にあたって,天野修一先生,阪上辰也先生の以前の入門者講習の資料を参考とさせていただきました。この場をお借りして御礼申し上げます。

 

また,Ustreamにて入門者講習の中継(録画)が試聴できますので,こちらもよろしければご覧ください。
http://www.ustream.tv/recorded/41452159

 

これで今年(2013年)の発表が全て終了したことになります。

残すは来月上旬の修論の提出です。年末年始返上でがんばります。