文系のための「相関係数のt検定」

相関係数について「0.2以下だから相関は無い」とか、逆に、
「0.7以上だから相関があると言える」などという人が居る。
中には、詳しい解説をせずに、そのように記述する教科書もある。

本当にそれで大丈夫なのだろうか？一つの疑問が沸々と湧いてくる…。
対象の数が少ないと「偶然」に相関係数が高くなることがあるのではないか？
そうなると、対象の数が少ない場合を考慮した方法を考えなくてはならない。

ということで、今回は確率の力を借りてこの問題を考えてみる。

確率と言えば、正規分布というものがあったが、
正規分布以外にも、分布には様々な種類が存在する。

様々にありすぎて覚えきれないほど、多く存在するのであるが、
相関係数の有意性を評価する方法の一つに「t 分布」を使う方法がある。

t 分布は、ある値を基準にどのように他のデータが散らばっているか？
という状況を表すための分布関数と言える。

正規分布では、平均を基準にしたときの対象のバラツキを基準としていたが、
t 分布では、その基準そのものに自由度を与えることができる。これが重要。
この性質のおかげで、対象数に影響されにくい確率密度を得ることができる。

では、さっそく、頭が痛くなるような数式から。

定義上、t 分布によって推定される確率は次のように与えられる。
無理に、覚える必要は無いが、そのような定義が存在していることを留意しておく。



さて、正規分布の式と比較すると、やや難解に見える。いや、かなり、難解であろう。
「ν」という記号と「x」があって、「」 というのも...やはり複雑である。

とにかく、この関数は「f(x)」と書いてあるので、
どうやら「x」によって何かの値（確率）を得ることができる関数であることがわかる。

一方、「ν」と「  」は、最初から定められるべき値（定数）であることが解る。
この二つの定数部分というのは、様々な状況によって異なる状況が与えられる。
今回は、相関係数の場合の話なので、これに応じた値が定められる。

では、この「ν」と「」の部分は何であるかか？という話になるのであるが、
正直に言うと説明が難い。誤解を与えないような「良い説明」が見つからない…。
仕方無いので、とりあえず、暫定的に以下のように説明しておくことにする。

まずは、「ν」 の値について。

今回の場合、相関係数を出すために、二つのデータの真ん中、
つまり、一つ目の変数と二つ目の変数の平均を基準としていて、
この二つの平均は、全ての対象からの計算から導かれる。

この二つの値は、実は、標本全体から計算の結果として導かれるため、
実際に必要となるのは、この二つの平均分を抜いて「n-2」が検定の対象となる。

このように、ある何らかの関係式において対象の数から、
計算によって得られる値の数を引いた値のことを「自由度」と呼ぶ。
したがって、これから行う相関係数のt分布は、自由度「n-2」のt分布に従う。

次に、「」と呼ばれるものであるが、これはベータ関数と呼ばれるものである。
これも、説明が困難なのであるが、ここでは、二つのパラメータを与えてやると、
何か値が出てくるようなものと理解しておくことにする。

ちなみに、「ν」を 4〜100の間で変化させた状況を可視化してみると以下のようになる。
これを見ると解るが、要するに、自由度が増えると値が小さくなるような値が得られる。

# νの取る範囲を設定する。

nu <- seq(4, 100)

# νの範囲に相当するβ関数を格納するベクトルを初期化する。
b <- vector(length=length(nu))

# νを4〜100の間で変化させながらβ関数から値を得る
for(i in 1:length(nu)){
    b[i] <- sqrt(i+3) * beta(0.5, (i+3)/2)
}

# 結果を可視化する
plot(nu, b, type="l", col="red")


解ったような、解らないような...イマイチな説明。エレガントさに欠ける。美しくない。
しかし仕方が無い。ご愛嬌。ここは思い切って、そいういうものとして理解。
この辺りの話については詳細な説明が必要になるのだが、い

さて、気持ち悪い状態ではあるが、要するに、正規分布とは異なる、
別の関数「t 分布」によって確率密度を求めることができる。
とりあえず、上の式からどのような分布を描けるかを実験してみる。

# 作図領域を初期化する
dev.off()

# グラフの線色と凡例ラベルを準備する
led <- c("DF=4","DF=10","DF=50", "NORM")
col <- c("green","orange","red", "blue")

# 色を呼び出すためのパラメータを初期化する
i <- 1

# xの範囲を決める。これはテキトー。とりあえず、-4〜4の範囲に。
X.range <- seq(-4, 4, 0.01)

# 自由度の異なるt分布を重ねて描く
for(nu in c(4,10,50)){
    fx <- (1/(sqrt(nu)*beta(0.5, nu/2)))*(1+(X.range^2/nu))^(-1*((nu+1)/2))
    par(new=TRUE)
    plot(X.range, fx, type="l", xlab="", ylab="", ylim=c(0, 0.4), xaxt="n", yaxt="n", col=col[i])
    i <- i+1
}

# ついでに正規分布も描く
par(new=TRUE)
plot(X.range, dnorm(X.range), type="l", xlab="", ylab="", ylim=c(0, 0.4), col=col[4])

# 図の体裁を整える
abline(h=0, v=0, col="grey", lty=2)
title(main="Student's t distribution", xlab="t-value", ylab="density")
legend("topleft", legend=led, col=col, lty=1)

この図において、DFが自由度を表し、異なる自由度の分布が描かれている。
なお、参考のために青い線で平均「0」標準偏差「1」の正規分布を重ねている。
なお、自由度が無限に大きくなると、t分布は正規分布に近似するのである。

ここでは、最初に示した定義に基づいてt分布を描いているが（「fx <- ・・・」の部分）、
実は、R には、正規分布と同様にt分布に関する関数を持っていて、
以下の関数のうち、dt()関数を使うことで簡単に描くことができる。

• dt()関数：ある値xの確率密度fxの値を返す。
• pt()関数：確率密度fxに対応する累積分布（p値）を返す。
• qt()関数：ある確率（p値）における値（t値）を返す。
• rt()関数：自由度νにおけるt分布に従う乱数をn個発生させる。

さてさて、いよいよ本題。相関係数の相関の有無を確認したい。

まずは、相関係数の式を思い出してみる。ちょっと、難しいか。
二つの変数の共分散を、それぞれの標準偏差で掛けたようなものになるのだった。



そして、この相関係数から、検定に用いる値を次のように与えられる。
これを「検定統計量 t」と呼ぶ。実は、今回の話で一番重要な式がコレ。



この値を計算すると、先程説明した自由度「n-2」の t分布に従うことが知られている。
したがって、逆に、この検定統計量 t を使って、
ある確率「p」における相関係数を得ることができる。



この式において、「 」は、ある確率pにおける相関係数という意味。
では、「ある確率」には、どのような基準を与えるかというと、
一般的には、1%、5%、10%のいずれか状況を考える。

この値を「有意水準」と呼ぶ。この水準の上下をめぐって、一喜一憂するのである。
もちろん、絶対的な基準ではないのであるが。

さて、相関係数の場合は正負の可能性があるので上の式に当てはめて、

負の相関の状況は、　となり、

正の相関の状況は、　となる。

ここで注意が必要。t検定を行う場合には、釣り鐘状の分布の上側の「袖」の部分、
下側の「袖」の部分、両袖の部分を取る場合の3つの場合があるが、
両袖を取る場合には、有意水準分を二つに分ける。

解りにくいかもしれないが、本ブログにおいては、
とりあえずは、両袖を取る場合を想定して話を進める。

ふむ。なるほど。イマイチ解りにくい。

ということで、作図して確認。とりあえず、自由度は「50」。

# xの範囲を決める。これはテキトー。とりあえず、-4〜4の範囲に。
X.range <- seq(-4, 4, 0.01)

# 自由度50のt分布を描く
plot(X.range, dt(X.range, 50), type="l", xlab="", ylab="", ylim=c(0, 0.4))

# 面積の部分を塗りつぶす
## 上限の片袖 2.5%
poly.x <- seq(qt(0.975,50),4,0.01)
poly.y <- c(0,dt(poly.x,50),0)
polygon(c(qt(0.975,50),poly.x,4), poly.y, col="red")

## 下限の片袖 2.5%
poly.x <- seq(-4, qt(0.025,50),0.001)
poly.y <- c(0,dt(poly.x,50),0)
polygon(c(-4,poly.x,qt(0.025,50)), poly.y, col="red")

## 信頼区間
poly.x <- seq(qt(0.025,50),qt(0.975,50),0.001)
poly.y <- c(0,dt(poly.x,50),0)
polygon(c(qt(0.025,50),poly.x,qt(0.975,50)), poly.y, col="lightblue")

# 図の体裁を整える
abline(h=0, v=0, col="grey", lty=2)
title(main="Student's t distribution with DF=50", xlab="range of value", ylab="density")

# 解りやすいようにラベルを配置
text(-3, 0.04, "2.5%", cex=1.5)
text(3, 0.04, "2.5%", cex=1.5)
text(0, 0.04, "95%", cex=1.5)

既に述べたように、相関係数のt分布は、自由度がn-2に従うので、
この場合は、対象数が52存在していることになる。

この図は、その時の有意水準5%（両袖）の場合を図式化したものである。
この図を見て解るように、「-0.025から下」と「0.975から上」を取っている。
つまり、両方合わせると、丁度全体の「5%」になるようになっている。

何を表しているかというと、手元の52の対象というのは、
母集団から抽出された標本であり、取り出し方を変える度に、
相関係数を算出すると、その相関係数の出現確率はこのように分布する。

ここからは、説明するよりも実際に値を見て検討した方が解りやすい。
自由度「50」の場合に、相関が起こり得る基準はどのような状況であるのか？
以下の三行の処理は、検定統計量 t から相関係数を計算したもの。

# 相関があると言って良い基準
abs(qt((0.01)/2, 50)/sqrt(50+qt((0.01)/2, 50)^2))    # 1%
abs(qt((0.05)/2, 50)/sqrt(50+qt((0.05)/2, 50)^2))    # 5%
abs(qt((0.10)/2, 50)/sqrt(50+qt((0.10)/2, 50)^2))    # 10%

今回は両袖で状況を観察するので、「 」を二つに分けている。
この実行結果は以下の通り。abs()関数は、絶対値を取る関数。

> abs(qt((0.01)/2, 50)/sqrt(50+qt((0.01)/2, 50)^2))    # 1%
[1] 0.3541529
> abs(qt((0.05)/2, 50)/sqrt(50+qt((0.05)/2, 50)^2))    # 5%
[1] 0.2732435
> abs(qt((0.10)/2, 50)/sqrt(50+qt((0.10)/2, 50)^2))    # 10%
[1] 0.2306199

したがって、対象数が50の場合には、

1%の有意水準の場合には0.35以上、

5%の有意水準の場合には0.27以上、

10%の有意水準の場合には0.23以上、

おっ？意外に低くても良い？

がおおよその基準となる。これ以下の場合には相関があるとは言い難い。

重要なポイントは、相関係数が「対象数に依存」するということであって、

一部の入門書で述べられているように、絶対的な基準は存在しないということ。

さて、もう少しこの状況を詳しく見るために、
簡単な実験を行なってみることにする。

以下の例では、異なる二つのサイコロを別々に振って、
二つのサイコロの結果に相関があるかを実験してみる。
当然、相関があるハズの無い実験であるが…結果を予測できるか？

パソコンによっては若干時間がかかるかもしれない。
さすがに、サイコロを「10,000回」振るのは大変。

# サイコロを振る回数
n <- 6

# 自由度はn-2
nu <- n-2

# 相関の有無の基準
r_99 <- abs(qt((0.01)/2, nu)/sqrt(nu+qt((0.01)/2, nu)^2))
r_95 <- abs(qt((0.05)/2, nu)/sqrt(nu+qt((0.05)/2, nu)^2))
r_90 <- abs(qt((0.10)/2, nu)/sqrt(nu+qt((0.10)/2, nu)^2))

# 結果を格納する行列を作成する
res <- matrix(nrow=1000, ncol=4)

for(i in seq(1, 10000)){
    # サイコロを二つ用意する
    d1 <- rep(0, n)
    d2 <- rep(0, n)

    # サイコロを振る（ゾロ目は計算上問題あるので回避）
    while(length(table(d1))==1) {d1 <- round(runif(n,1,6))}
    while(length(table(d2))==1) {d2 <- round(runif(n,1,6))}

    # 相関係数を算出する
    r <- abs(round(cor(d1,d2),3))

    # 検定統計量tを算出する
    t <- round((r*sqrt(n-2))/sqrt(1-(r^2)),3)

    # p値を計算する（t値が正の場合は、全体からp値を引く）
    p <- 2 * min(pt(t, nu), 1 - pt(t, nu))

    # 相関の有無の判定
    if(r>r_99){res_99<-"correlated"}else{res_99<-"uncorrelated"}
    if(r>r_95){res_95<-"correlated"}else{res_95<-"uncorrelated"}
    if(r>r_90){res_90<-"correlated"}else{res_90<-"uncorrelated"}
    if(r>0.7){res_gen<-"correlated"}else{res_gen<-"uncorrelated"}

    # 判定結果を格納する
    res[i, 1] <- res_99
    res[i, 2] <- res_95
    res[i, 3] <- res_90
    res[i, 4] <- res_gen
}

さて、この実験では、10,000回ほど同じ処理を行ったわけであるが、
「相関有り」とみなされた場合は、何回含まれているのであろうか？

table(res[,1])    # 1%有意水準で相関係数を判定した場合
table(res[,2])    # 5%有意水準で相関係数を判定した場合
table(res[,3])    # 10%有意水準で相関係数を判定した場合
table(res[,4])    # 相関係数を0.7で判定した場合

ここでは、table()関数を用いている。この関数を使うと同じ値ごとに集計してくれる。
非常に便利なのでExcel のデータをコピーしてこの関数で集計することもある。
話が逸れたが、とにかく、結果を表示してみる。

> table(res[,1])    # 1%有意水準で相関係数を判定した場合
  correlated uncorrelated 
          92         9908 

> table(res[,2])    # 5%有意水準で相関係数を判定した場合
  correlated uncorrelated 
         527         9473 

> table(res[,3])    # 10%有意水準で相関係数を判定した場合
  correlated uncorrelated 
        1072         8928 

> table(res[,4])    # 相関係数を0.7で判定した場合
  correlated uncorrelated 
        1285         8715 

この結果から解るように、有意水準「1%」というのは、
繰り返し標本抽出した場合に、誤判別が1%ほど含まれることを意味する。
同様に5%と10%の場合にも、同様のことが言えるのが確認できる。

では、一般的に言われているように相関係数の基準を「0.7」とした場合、
12.85%も「誤判別」が生じていることが解る。この結果を可視化してみる。

# 結果をパーセンテージで表示
cls.F <- round(table(res[,4])[1]/10,3)
cls.T <- round(table(res[,4])[2]/10,3)

# 円グラフを描画する
pie(table(res[,4]),  col=c("red", "green"), clockwise=TRUE)

# 図を整える
title(main="Misclassification of correlation (r=0.7)")
text(0.2,0.5,cls.F, cex=1.8)
text(-0.2,-0.3,cls.T, cex=1.8)
mtext("Total size: 10 000", side=1, adj=1)

この結果は乱数によるものなので、実行の度に値が変わるが、
おおよそ、似た結果になるであろう。

この実験は「n=6」 の場合であるが、「n=100」の場合はどうなるであろうか？
おそらく、この実験の結果では、誤判別は「0」となるであろう。

先程の実験結果の教訓とは、全く逆のことが言えるのである。
t分布は、n-2の自由度に依存するため、対象数が増えると、
相関係数は「0.7」よりもずっと低くても良いのである。「n=52」の例も低かった。

どういう意味か？ちょっとした、「早見表」を作ってみる。

r_table <- matrix(nrow=50, ncol=3)
rownames(r_table) <- c(seq(4,50),100,1000,10000)
colnames(r_table) <- c("0.01","0.05","0.10")
counter <- 1
for(n in c(seq(4,50),100,1000,10000)){
    # 相関の有無の基準
    r_table[counter, 1] <- abs(qt((0.01)/2, n-2)/sqrt(n-2+qt((0.01)/2, n-2)^2))
    r_table[counter, 2] <- abs(qt((0.05)/2, n-2)/sqrt(n-2+qt((0.05)/2, n-2)^2))
    r_table[counter, 3] <- abs(qt((0.10)/2, n-2)/sqrt(n-2+qt((0.10)/2, n-2)^2))
    counter <- counter + 1
}

この処理の実行結果は以下の通り。

> r_table
            0.01       0.05       0.10
4     0.99000000 0.95000000 0.90000000
5     0.95873500 0.87833945 0.80538364
6     0.91719970 0.81140135 0.72929928
7     0.87452638 0.75449223 0.66943947
8     0.83434163 0.70673440 0.62148925
9     0.79768120 0.66638361 0.58220560
10    0.76459250 0.63189686 0.54935683
11    0.73478634 0.60206878 0.52140436
12    0.70788755 0.57598299 0.49726475
13    0.68352763 0.55294266 0.47615599
14    0.66137560 0.53241280 0.45750017
15    0.64114481 0.51397748 0.44086084
16    0.62259073 0.49730904 0.42590199
17    0.60550592 0.48214602 0.41236048
18    0.58971445 0.46827731 0.40002705
19    0.57506679 0.45553051 0.38873305
20    0.56143540 0.44376340 0.37834086
21    0.54871103 0.43285756 0.36873700
22    0.53679962 0.42271350 0.35982694
23    0.52561988 0.41324703 0.35153123
24    0.51510117 0.40438632 0.34378257
25    0.50518184 0.39606973 0.33652351
26    0.49580785 0.38824400 0.32970471
27    0.48693164 0.38086286 0.32328346
28    0.47851116 0.37388591 0.31722267
29    0.47050913 0.36727768 0.31148989
30    0.46289233 0.36100691 0.30605660
31    0.45563108 0.35504589 0.30089765
32    0.44869879 0.34937001 0.29599073
33    0.44207154 0.34395729 0.29131600
34    0.43572774 0.33878805 0.28685573
35    0.42964790 0.33384462 0.28259402
36    0.42381429 0.32911104 0.27851658
37    0.41821085 0.32457292 0.27461052
38    0.41282290 0.32021717 0.27086418
39    0.40763704 0.31603193 0.26726695
40    0.40264101 0.31200637 0.26380923
41    0.39782355 0.30813060 0.26048222
42    0.39317433 0.30439558 0.25727790
43    0.38868380 0.30079300 0.25418891
44    0.38434318 0.29731521 0.25120851
45    0.38014434 0.29395520 0.24833047
46    0.37607975 0.29070645 0.24554908
47    0.37214243 0.28756298 0.24285905
48    0.36832589 0.28451923 0.24025548
49    0.36462410 0.28157003 0.23773384
50    0.36103143 0.27871059 0.23528994
100   0.25648345 0.19655120 0.16542980
1000  0.08142144 0.06199742 0.05204468
10000 0.02575724 0.01960021 0.01644948

この結果を見てみると、対象数が「100」以上だと「0.25」でも相関があると言えてしまう。
したがって、十分に対象数が存在する場合には「0.7」にこだわる必要が無いと言える。

相関が「0.3〜0.7」とする表記を見かけることがあるが、
この基準というのは、対象数が「10件〜100件」が目安ということであろうか。

さて、今回は、少々難しい話が多かったかもしれない。
実際、説明する側としても、かなりの苦労を要した。
最後に、少しだけ、重要な補足を加えることにする。本当に、少しだけ。

今回の方法というのは、実際には少々、粗い方法であって、
実際には、もう少し気の利いた補正というのが必要になるし、
また、t分布以外の分布によって、相関係数の妥当性を評価する方法もある。

Rには、相関係数の妥当性を評価する関数にcor.test()関数というのがある。
この関数は、二つのベクトルを与えると、上記のt分布を用いた方法に加えて、
巧妙な補正をかけてより精度の高いに相関係数の検定を行なってくれる。

試しに、サイコロの実験の例を用いてみる。

# サイコロを振る回数
n <- 100


# サイコロを二つ用意して、サイコロを振る
d1 <- round(runif(n,1,6))
d2 <- round(runif(n,1,6))

# 相関係数の検定を行う
cor.test(d1, d2, alternative="two.sided", method="pearson", conf.level = 0.95)

実行結果は以下の通り。

> cor.test(d1, d2, alternative="two.sided", method="pearson", conf.level = 0.95)

Pearson's product-moment correlation

data:  d1 and d2 
t = -1.3222, df = 98, p-value = 0.1892
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 -0.32046720  0.06574572 
sample estimates:
       cor 
-0.1323819 

cor.test()関数では、相関を見るための二つの変数に格納されたベクトルを入力とし、
それぞれのベクトル自身の分布を使って相関の基準となるの範囲を調整している。
実際のデータで相関係数の検定を行うのであれば、この関数を使うと良い。

それにしても、この結果、どのようにして見るのか解りづらい。


この関数では、「相関が0である」という仮説（帰無仮説）と、
この仮説に対する「相関が0では無い」という仮説（対立仮説）を用意している。
重要な問題は、帰無仮説は棄却され、対立仮説が採択されるかどうか？である。

帰無仮説というのは「無に帰する仮説」という意味。実は、言いたくないこと。
では、本当に言いたいことは対立仮説の方になる。どのように判断するのか？

まず、この結果において注意するべき値は「95 percent confidence interval:」部分と、
p-value（p値）と呼ばれる値が書いてある部分。この二つを観察する。

まず、「95 percent confidence interval:」と書いてある部分は、
95%信頼区間と言われる範囲。
対象となっている変数の、その信頼度での属する区間を表している。


元の相関係数は、95%の信頼度で「-0.32〜0.065」の間に属していることになる。
この範囲に「0」が含まれているので、相関係数が「０」の可能性もあり得る。
つまり、相関無しの可能性を否定はできないということになる。

帰無仮説は「相関が0である」であり、この仮説は棄却したいのだけれど、
今回の検定の結果は、相関が0である可能性が95%信頼区間に入ってしまっているので、
残難ながら、帰無仮説を棄却することはできない。

cor.test()関数の場合は、既定の設定では95%信頼区間で見るので、
p値を見てみると0.05よりも大きいから、棄却できないことは解る。
しかし、相関係数がどのくらいに見積もるかは信頼区間を観察する必要がある。