統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

相関係数の検定と信頼区間

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

統計ソフトRで相関係数の検定と推定はcor.test()でできるが、個々のデータが必要だ。

個々のデータを使わなくても、検定や推定はできないだろうか?

関数を作ってみたので、共有したい。

cor.test()の使い方はこちら。

toukeier.hatenablog.com

引用書籍

「新版医学への統計学」の6.1 相関係数の検定と信頼区間の例題を確認していく。

第3版も出てる。

相関係数  \rho の検定

相関係数  \rho がゼロかどうかの検定。スクリプトは以下の通り。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))
}

例題6.1のヘモグロビンー尿酸間の検定結果は、以下の通り計算される。

> pop.rho.test(r=-0.533,n=35)
[[1]]
            t            df       p-value 
-3.6187173384 33.0000000000  0.0009787127 

相関係数  \rho の信頼区間

相関係数  \rho の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 

Fisherの分散安定化変換は、hyperbolic arc-tangent (atanh())で計算できる

上述のスクリプト中の  Z の式は、Fisherの分散安定化変換と呼ばれるものである。これはhyperbolic tangenの逆関数のhyperbolic arc-tangentで、統計ソフトRではatanh()で計算できる。

> r <- -0.533
> 
> 1/2*log((1+r)/(1-r))
[1] -0.5943263
> 
> atanh(r)
[1] -0.5943263

二つの母相関係数の差の検定

二つの母相関係数  \rho_A \rho_B の差の検定のスクリプトは以下の通り。

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))
}

例題6.3に沿って、男性と女性の相関係数の差を検定した結果は以下の通り。

> 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 

まとめ

統計ソフトRで、相関係数の検定と信頼区間スクリプトを紹介した。

個々のデータが使える場合は、cor.test()で検定と信頼区間は簡単に得られるが、相関係数とサンプルサイズしかわからない場合は、cor.test()は使えないので、今回紹介したスクリプトが役に立つ。

また、2つの相関係数の差の検定は、デフォルトでは装備されていないので、今回のスクリプトが使える。よりスマートにまとまっているスクリプトはcor.diff.test()で検索すると見つかる。brainGraphパッケージに含まれている。