R で相関係数の検定と推定は cor.test() でできるが、個々のデータが必要だ。
個々のデータを使わなくても、検定や推定はできないだろうか?
相関係数の検定
母相関係数 がゼロかどうかの検定。
スクリプトは以下の通り。
r がサンプルの相関係数で、n がサンプルサイズ。
pop.rho.test <- function(r,n){ T <- r/(sqrt(1-r^2)/sqrt(n-2)) p <- pt(abs(T), df=n-2,lower.tail=F)*2 list(c(t=T, df=n-2, "p-value"=p)) }
r = -0.533, n = 53 のときの検定結果は、以下の通り計算される。
> pop.rho.test(r=-0.533,n=35) [[1]] t df p-value -3.6187173384 33.0000000000 0.0009787127
p = 0.0009787127 で統計学的有意に母相関係数は 0 と異なると言える。
相関係数の信頼区間
母相関係数 の95%信頼区間を計算するスクリプトは以下の通り。
pop.rho.ci <- function(r,n){ Z <- 1/2*log((1+r)/(1-r)) Zl <- Z-1/sqrt(n-3)*qnorm(1-0.05/2) Zu <- Z+1/sqrt(n-3)*qnorm(1-0.05/2) rl <- (exp(2*Zl)-1)/(exp(2*Zl)+1) ru <- (exp(2*Zu)-1)/(exp(2*Zu)+1) list(c(Estimate=r, LL=rl, UL=ru)) }
上述の相関係数の信頼区間は以下のように計算される。
> pop.rho.ci(r=-0.533,n=35) [[1]] Estimate LL UL -0.5330000 -0.7355906 -0.2428969
下限が -0.7355906, 上限が -0.2428969 であった。
相関係数の信頼区間を計算する際のFisherの分散安定化変換
上述のスクリプト中の の式は、Fisherの分散安定化変換と呼ばれるものである。
これは hyperbolic tangent の逆関数の hyperbolic arc-tangent で、R では atanh() で計算できる。
> r <- -0.533 > > 1/2*log((1+r)/(1-r)) [1] -0.5943263 > > atanh(r) [1] -0.5943263
相関係数の差の検定
二つの母相関係数 と の差の検定のスクリプトは以下の通り。
pop.rho.diff <- function(rA,nA,rB,nB){ T <- abs((atanh(rA)-atanh(rB)))/sqrt(1/(nA-3)+1/(nB-3)) p <- 2*(1-pnorm(T)) list(c(rA=rA,nA=nA,rB=rB,nB=nB,T=T,"p-value"=p)) }
男性 -0.459 と女性 -0.097 の相関係数の差を検定した結果は以下の通り。
> pop.rho.diff(rA=-0.459,nA=30,rB=-0.097,nB=24) [[1]] rA nA rB nB T p-value -0.4590000 30.0000000 -0.0970000 24.0000000 1.3704342 0.1705514
p = 0.1705514 で統計学的有意ではなかったが、-0.459 と -0.097 はかなり大きな違いである。
サンプルサイズが小さすぎたことによる検出力不足と考えられる。
まとめ
Rで、相関係数の検定と信頼区間のスクリプトを紹介した。
個々のデータが使える場合は、cor.test() で検定と信頼区間は簡単に得られるが、相関係数とサンプルサイズしかわからない場合は、cor.test()は使えないので、今回紹介したスクリプトが役に立つ。
また、2つの相関係数の差の検定は、デフォルトでは装備されていないので、今回のスクリプトが使える。
よりスマートにまとまっているスクリプトは cor.diff.test() で検索すると見つかる。
brainGraphパッケージに含まれている。
cor.diff.test function – RDocumentation
関連記事
cor.test()の使い方はこちら。
引用書籍
出典は、「新版医学への統計学」の6.1 相関係数の検定と信頼区間の例題
第3版も出てる。
コメント