裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

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

2016年03月14日 | ブログラミング

R で,cor.test(x, y, method="spearman") とやれば,スピアマンの順位相関係数が表示されるが,信頼区間は表示されない。さて,どうしたものかというのがスタート地点。

上のように呼んだら,R は,x, y の代わりに rank(x), rank(y) を引数としてピアソンの積率相関係数を計算する。そして(理論的に)それは,スピアマンの順位相関係数に他ならない(実際,計算結果は一致する)。

では,cor.test(rank(x), rank(y)) としたらどうなるか。cor.test は x, y の測定値そのままではなくて,順位をデータだと思ってピアソンの積率相関係数を計算する。R は,計算されるのはスピアマンの順位相関係数だとわかっているからそのように計算する。しかし,ピアソンの積率相関係数の信頼区間は計算しない。

なぜ,計算しないか?順位をデータだと思ってピアソンの積率相関係数を計算するが,そのまま信頼区間の計算をしてよいかどうかについては「否」としているんだろう。

なぜか?元のデータは少なくとも正規分布するけど,順位に変換したデータは正規分布しないじゃないか。だったら,信頼区間は計算できないだろうということかな?

 

以下にシミュレーション・プログラムを示す。アイリス・データファイルからデータを取って,ブートストラップにより信頼区間を求めるプログラムである。(プログラムが間違っていたら元も子もないが)

当然ながら,両者が完全に一致するわけはないが,近似計算としては十分な意味を持つと思うがいかが?

CL.app というのは,データの代わりにデータの順位を与えてピアソンの積率相関係数の信頼区間を計算する方式で頸sんした信頼区間

CL.boot というのブートストラップによる信頼区間

まあ,結論は,どっちでもいいけどね(^_^;)

# ブートストラップ回数
N = 10000

# iris データセットの iris[1:50, 1:2] から,サンプルサイズ 30 で無作為抽出
j = sample(50, 30)
x = iris[j, 1]
y = iris[j, 2]

# rank(x), rank(y) で cor.test を呼べば,cor.test は順位相関係数を計算する
# おまけに,信頼区間も
a = cor.test(rank(x), rank(y))
a$estimate # スピアマンの順位相関係数(ただし,表示はピアソンの積率相関係数)
a$conf.int # 信頼区間
# ブートストラップ
z = sort(atanh(sapply(seq_len(N), function(i) {
    i = sample(30, 30, replace=TRUE)
    cor.test(x[i], y[i])$estimate
    })))
tanh.z = tanh(z)
# ブートストラップによる信頼区間
boot = tanh.z[N*c(0.025, 0.975)]
layout(1:2)
par(mar=c(3, 3, 1.2, 0.5), mgp=c(1.8, 0.8, 0))
hist(z, breaks=30, probability=TRUE, main="Fisher's z transformation")
x = seq(0.2, 1.6, 0.01)
lines(x, dnorm(x, mean(z), sd(z)), col="red")
hist(tanh.z, breaks=30, probability=TRUE, main="Correlation coefficients")
legend("topleft", legend=c(
sprintf("Spearman's corr. =%.5f", a$estimate),
sprintf("CL.app =[%.5f, %.5f]", a$conf.int[1], a$conf.int[2]),
sprintf("CL.boot=[%.5f, %.5f]", boot[1], boot[2])),
bty="n")
layout(1)

コメント   この記事についてブログを書く
« 平成27年3月24日の金星は(そ... | トップ | 効率のよい荷物配送法 »
最新の画像もっと見る

コメントを投稿

ブログ作成者から承認されるまでコメントは反映されません。

ブログラミング」カテゴリの最新記事