相関の検定

目次

・相関とはなにか？

・相関係数：相関を示す量的な基準

　　・Pearsonの相関係数R

　　・決定係数R^2

　　・Spearmanの順位相関係数ρ

　　・Kendallの順位相関係数τ

・相関係数の検定

・どれくらいの相関係数なら「意味」があるか？

・Rで行う相関の検定

・相関と因果（工事中）

・追記：同順があるときのSpearmanのρ

・追記：同順があるときのKendallのτ

相関とはなにか？

　私たちは、2つの変数に「何かしらの関係」がないかを、しばしば探している。この「何かしらの関係」を足掛かりに、さらなる実験や解析を組んだり、結果の解釈に取り組んだりすることになるわけだ。統計をやっている者が考える「何かしらの関係」の代表的なものが、相関correlationである。相関とは、下の図に示すように「一方の変数が増加した時に、他方が減少するような関係」やその逆の関係がみられる状況を指す。

上記で示した、特に重要な相関が、線形関係の相関である。線形関係とは、2つの変数が直線関係にある場合を指す。他にも、下図のような非線形な相関も存在する。直線関係ではないが、明らかに2つの変数間では何かの関係がある。

相関とは、「何かの回帰式に当てはめたとき、回帰曲線と実測値の誤差が小さい状況」と言い換えることもできよう。本項では、線形関係の相関について解説をする。

　線形関係の相関を大まかに分ければ、以下の3パターンがある。

1. 正の相関：一方の変数が増加した時に、他方が増加するような関係

2. 無相関：一方の変数が増加した時に、他方が増加も減少もしないような関係

3. 負の相関：一方の変数が増加した時に、他方が減少するような関係

後述するような、ある量的な基準に基づき、正の相関や負の相関が観察された場合、その状況を「（正の or 負の）相関がある」と表現する。もし、無相関なら「相関がない」と表現する。

　相関と類似した指標として、線形回帰の時の傾きの係数が存在する。そのとき、係数が0に近ければ、ほとんど相関がないように見える。実際、傾きが小さいことを無相関と表現することもある。ただし、下の図のように「2変数が互いにランダムの場合」と「正負の関係は小さいが誤差が小さい場合」があることを覚えておくとよい。前者では、真の意味で2変数間に関係がない可能性が高いが、後者は関係がないとはいえない可能性が高い（ここの議論は間違っていたので2024.8.18に修正）。

データの状態としては明らかに異なっているが、必ずしも線形回帰の係数だけではその信頼度がわからないことがある。それに対して、後述する相関を示す量的基準は両者を区別できる可能性が高い。

相関係数：相関を示す量的な基準

・Pearsonの相関係数R

　前述してきた線形相関を示す量的な基準、は相関係数correlation coefficientと呼ばれる。まずは、Pearson（ピアソン）の積率相関係数Pearson's product-moment correlation coefficientを紹介しよう。この積率相関係数は、しばしば単に相関係数とも呼ばれる。Pearsonの相関係数はしばしばRと表記され、以下のように2変数の分散と共分散を使って定義される。

この時、分散と共分散は、基準がそろっていれば計算結果は変わらないので、標本でも不偏でもよい。つまり、不偏分散を使ったら不偏共分散を使えばよいし、標本分散を使ったら標本共分散を使えばよい。

　Pearsonの相関係数Rは、-1から1までの値をとる。定義から、相関係数Rの正負は、定義から分子の共分散に依存している。共分散の性質から、一方の変数が増加した時、他方が増加すれば正に、減少すれば負になる。ゆえに、Rが正になるなら正の相関があると言えるし、負なら負の相関があると言える。共分散は、変数の絶対値の大きさに依存して大きく変化するので、それを-1から1までに規格化するために分散で除している。

　線形回帰の傾きの推定値と比較すると、分母にy側の標準偏差が存在するか否かが異なっている。つまり、線形回帰の傾きはxのばらつきだけを考慮するが、相関係数はxとyの両方のばらつきを考慮する。これにより、上述した「2変数が互いにランダムの場合」と「正負の関係は小さいが誤差が小さい場合」について、前者の場合は傾きもRも0に近い値をとるが、後者の場合は傾きは小さい値をとるがRは大きい値をとる。なぜなら、後者は、yを考慮しない傾きには影響がないが、yのばらつきが小さいのでRの分母が小さくなり、それに合わせてRが大きくなるからである。逆に、Rが大きいからと言って必ずしも線形回帰の傾きが大きいとは限らず、ゆえに強い効果のある変数とは限らないこともあるので注意が必要である。

　また、この直後に解説する決定係数を見てもらう方がわかりやすいが、R = ±1となるときは、全くの誤差がなく回帰直線上にデータがのっている状況である。

・決定係数R^2

　Pearsonの相関係数と密接に関係した基準として決定係数coefficient of determinationが存在する。決定係数にはいくつかの定義（Wikipedia, 2024/8/15閲覧）があるが、ここで紹介するものは、その中でも特に一般的で、ソフトウェアRで算出される下記の定義を紹介する。

これは下記で示すように

(全変動) ＝ (回帰変動) + (残差変動)

とあらわせることを利用して、この式の全体を全変動で割り、

1 =  (回帰変動)/(全変動) + (残差変動)/(全変動)

と式変形し、

1 - (残差変動)/(全変動) = (回帰変動)/(全変動) = 決定係数

というように定義している。

ここで、線形回帰の結果を用いると

ゆえに、

すなわち、決定係数とは、全変動に対する回帰直線で説明できる変動の割合を表している。

　さらに、この定義の決定係数はR^2とあらわされるように、相関係数の2乗と同じ値となる。

決定係数R^2 = 1となるとき、誤差変動が0であるから、データは一直線状に並ぶ。ゆえに、相関係数が±1となるときは、誤差変動がないときに限られる。

・Spearmanの順位相関係数ρ

　Pearsonの相関係数および決定係数は、その定義から線形回帰と密接に関係する。線形回帰が、誤差が正規分布に従うことを仮定した解析であるため、同様にPearsonの相関係数や決定係数も、誤差が正規分布に従う前提のもとで有効な指標である。では、誤差がコーシー分布のように外れ値が出やすいデータの場合はどうするべきか。そんな時に使うのが、Spearman（スピアマン）の順位相関係数Spearman's rank correlation coefficientであり、ギリシャ文字ρ（ロー）で表されることが多い。その定義は単純で、2変数のデータを順位データに変換したのちに、Pearsonの相関係数と同様の計算をすることで求める（ただし、同順（タイ）データがない場合）（同順がある場合は追記）。2変数xとyを順位に変換したものをrx、ryすればρは以下のように求める。併せて同値な計算法を2つ示す。

ρの計算はサイトによって異なるように見える計算方法が紹介されているが、実はすべて同値である。以下に証明を示す。

・Kendallの順位相関係数τ

　Spearmanの順位相関係数と同様の傾向を示すものにKendall（ケンドール）の順位相関係数があり、こちらはギリシャ文字τ（タウ）で表されることが多い。同順がない場合、以下のように定義される（同順がある場合は追記）。

文字で書かれるとわかりにくいが、図にするとその意味は明白だ。

あるデータ（Xi, Yi）から見たときに、右上or左下にあるデータの数がPiで、左上or右下にあるデータの数がQiである。PとQはこれらの総和である。

　PとQはともに0からNまでの値をとり、N = P + Qという関係があるため、τは-1から1までの値をとる。τ = 1のとき、P = N、すなわち、2変数の順位の方向がすべて一致する場合である。逆にτ = -1のとき、Ｑ = N、すなわち、すべての順位の方向がすべて一致しない＝順位が逆方向の場合である。

相関係数の検定

　前述まで、相関係数のいくつかの定義について述べてきた。では、その相関の程度は、どれだけ信用に足る値かを知りたくなるだろう。これらの相関係数についても、平均値の差などと同様に検定を行い、偶然から考えられる以上に逸脱した値であるかを量的に示すことができる。

　サンプルサイズがnのとき、Pearsonの相関係数やサンプルサイズが十分に大きい場合のSpearmanの順位相関係数は、以下のような検定統計量tが自由度n - 2のt分布に従うことを利用して検定を行う。Rは相関係数を示す。

tは、相関係数の絶対値が1に近いほど、分子は大きく、分母は小さくなる。つまり、tは非常に大きくなり、有意になりやすい。また、サンプルサイズも大きいほどtは大きくなるので、同じ相関係数であってもサンプルサイズが大きいほど有意になりやすい。

　なぜ、相関係数の検定にtが登場するのだろうか。実は、このtは線形回帰のときの、傾きに関する検定から求めることができる。線形回帰では、傾きに関する有意性に関して、その係数をその係数の標準偏差で割った値がt分布に従うことを利用していたのだった。そのtは上記の相関係数のtと一致する。

つまりは、線形回帰を行ったとき、その傾きの係数が有意なら、相関係数も有意であると判断できる。また、p値も一致する。まあ、互いの定義と図形的状況を考えれば当然である。

　上記で、無相関の状況に2つのパターン「2変数が互いにランダムの場合」と「正負の関係はないが誤差が小さい場合」があることを紹介した。仮にRが傾きが0に近くても、前者の場合、線形回帰の傾きの係数の標準偏差が大きくなり、ほとんど有意になることはないだろう。一方で、後者なら標準偏差が小さくなって、有意になる可能性もある。

どれくらいの相関係数なら「意味」があるか？

　相関係数について、しばしば〇〇 < |R|なら強い相関とか、|R| < △△なら弱い相関とかの議論があるが、これにはどれくらい妥当性があるのだろうか。ここでは、有意になるか否か、という観点から、相関係数の大きさと、サンプルサイズの関係を示してみる。

　相関係数の検定統計量の計算式をRについて解きなおす。

　明らかにtの絶対値が大きくなれば、Rは±1に近づき、nが大きくなれば0に近づく。しかし、実際には、tもサンプルサイズの影響を受けているから、以下のようにこのRの計算について、サンプルサイズの影響を正しく評価しよう。今回は、正のほうだけを考慮して計算しているが、負の相関の場合はすべての符号が逆転しているものだと考えてもらってよい（あるいは、絶対値付き計算の結果だと思ってもらってもよい）。

------------------------------------------------------

library(plotn)

library(gtools)

n <- 3:100

R90 <- qt(p = 0.95, df = n - 2)/sqrt(qt(p = 0.95, df = n - 2)^2 + n - 2)#90%棄却域

R95 <- qt(p = 0.975, df = n - 2)/sqrt(qt(p = 0.975, df = n - 2)^2 + n - 2)#95%棄却域

R99 <- qt(p = 0.995, df = n - 2)/sqrt(qt(p = 0.995, df = n - 2)^2 + n - 2)#99%棄却域

R999 <- qt(p = 0.9995, df = n - 2)/sqrt(qt(p = 0.9995, df = n - 2)^2 + n - 2)#99.9%棄却域

Rs <- data.frame(n,R90,R95,R99,R999)

plotn(Rs[,1], Rs[,2:5], mode = "m", type = "l", #図1の描画

      xlab = "Sample size", ylab = "Pearson's R",

      xlim = c(0, 120),

      lty = c(2,1,2,2), col.dot = col_genelator()[c(1,4,1,1)], 

      lwd.dot = c(1, 1.5, 1, 1))

overdraw(text(x = 110, y = R90[98], labels = "90%"),

         text(x = 110, y = R95[98], labels = "95%"),

         text(x = 110, y = R99[98], labels = "99%"),

         text(x = 110, y = R999[98], labels = "99.9%"),

         text(x = 105, y = 0.4, labels = "Rejection region"))

------------------------------------------------------

結果、サンプルサイズが大きくなるにつれて、有意になる相関係数の値は小さくなってゆく。特に、3～20にかけてのRの減少はかなり大きい。上の図では、赤線よりも上側の相関係数となれば、検定では有意になる。

　20、30、40のサンプルサイズについて確認してみると、

------------------------------------------------------

Rs[Rs$s == 20,]

##     s       R90       R95       R99      R999

## 18 20 0.3783409 0.4437634 0.5614354 0.6787811

Rs[Rs$s == 30,]

##     s       R90       R95       R99      R999

## 28 30 0.3060566 0.3610069 0.4628923 0.5703174

Rs[Rs$s == 40,]

##     s       R90       R95      R99      R999

## 38 40 0.2638092 0.3120064 0.402641 0.5007004

------------------------------------------------------

　それぞれ、R = 0.444、0.361、0.312よりも大きければ、有意な相関となることがわかる。この値をみてどう思っただろうか。私は、思ったよりも小さな相関係数で有意になるものだな、と感じている。

　同じ相関係数R = 0.4だとしても、その意味は、サンプルサイズによって劇的に解釈が変わる。ゆえに、私の見解としては、Rの絶対値の大きさによって強い相関とか、弱い相関とかを区別する意味はあまりないように感じている。相関係数にしろ、線形回帰にしろ、t値なり、サンプルサイズなりを記載することに変わりはない。それに、そこに科学的観点からどのような解釈をするべきかは、また別問題である。

Rで行う相関の検定

　では、Rを使って相関の検定を行ってみよう。以下のように正規分布に従うデータを生成する。

------------------------------------------------------

library(plotn)

library(gtools)

n <- 30

b0 <- 10

b1 <- -1

x <- round(runif(n, min = 10, max = 30), digits = 2)

y <- round(rnorm(n, mean = b1 * x + b0, sd = 2), digits = 2)

plotn(y ~ x) #図2の描画

------------------------------------------------------

負の相関があることは明らかだろう。相関係数の検定はcor.test関数を使うことで実行できる。今回は、外れ値などはなさそうなのでPearsonの相関係数を求めることにする。

------------------------------------------------------

cor.test(x, y, method = "pearson")

## 

##  Pearson's product-moment correlation

## 

## data:  x and y

## t = -16.077, df = 28, p-value = 1.138e-15

## alternative hypothesis: true correlation is not equal to 0

## 95 percent confidence interval:

##  -0.9761085 -0.8963415

## sample estimates:

##        cor 

## -0.9498735

------------------------------------------------------

この結果から、Pearsonの相関係数は-0.95であり、負の相関を示している。相関係数だけなら、以下のようにcor関数を使ったり、手計算で求めることもできる。

------------------------------------------------------

cor(x, y, method = "pearson")

## [1] -0.9498735

cov(x, y)/sqrt(var(x)*var(y))

## [1] -0.9498735

------------------------------------------------------

t値とp値は以下のように求めることができる。自由度はサンプルサイズが30だから、28だ。

------------------------------------------------------

cor(x, y)*sqrt(n-2)/sqrt(1 - cor(x, y)^2)

## [1] -16.07711

pt(q = cor(x, y, method = "pearson")*sqrt(n-2)/sqrt(1 - cor(x, y, method = "pearson")^2), df = n - 2) * 2

## [1] 1.138395e-15

------------------------------------------------------

この結果は線形回帰と一致することを示そう。

------------------------------------------------------

fit1_1 <- lm(y ~ x)

fit1_2 <- summary(fit1_1)

fit1_2

## 

## Call:

## lm(formula = y ~ x)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -4.9231 -1.4703  0.3305  1.1604  4.0463 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept) 12.91618    1.50774   8.567 2.61e-09 ***

## x           -1.14754    0.07138 -16.077 1.14e-15 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 2.02 on 28 degrees of freedom

## Multiple R-squared:  0.9023, Adjusted R-squared:  0.8988 

## F-statistic: 258.5 on 1 and 28 DF,  p-value: 1.138e-15

------------------------------------------------------

確かに、xの係数に関するt値とp値は一致している。また、Multiple R-squaredは決定係数のことであり、以下の値と一致する。

------------------------------------------------------

cor(x, y, method = "pearson")^2

## [1] 0.9022597

------------------------------------------------------

　続いて、外れ値が出やすいコーシー分布に従うようなデータを生成しよう。

------------------------------------------------------

x <- round(runif(n, min = 10, max = 30), digits = 2)

y <- round(rcauchy(n, location = b1 * x + b0, scale = 3), digits = 2)

plotn(y ~ x) #図3の描画

------------------------------------------------------

さて、図1同様、負の相関がありそうだが、外れ値が気になるところだ。この時にPearsonの相関係数を求めると以下のようになる。

------------------------------------------------------

cor.test(x, y, method = "pearson")

## 

##  Pearson's product-moment correlation

## 

## data:  x and y

## t = -3.7315, df = 28, p-value = 0.0008591

## alternative hypothesis: true correlation is not equal to 0

## 95 percent confidence interval:

##  -0.7755488 -0.2726393

## sample estimates:

##        cor 

## -0.5763021

------------------------------------------------------

有意ではあるが、R = -0.58程度と、やや低めである。これは外れ値の影響を受けている結果だろう。対して、Spearmanの順位相関係数では外れ値の影響を受けにくい。以下ではまず、正確検定を行う。

------------------------------------------------------

cor.test(x, y, method = "spearman", exact = T)

## 

##  Spearman's rank correlation rho

## 

## data:  x and y

## S = 7666, p-value = 2.261e-05

## alternative hypothesis: true rho is not equal to 0

## sample estimates:

##        rho 

## -0.7054505

------------------------------------------------------

ρ = -0.71となり、より負の相関が際立った形となった。シミュレーションしてみるとわかるが、外れ値が多いときには、SpearmanのρやKendallのτなどを使う方が、検出力が高い。

　相関係数だけを求めるなら、以下のようにcor関数を使ったり、データを順位に変換したのちにいくつかの方法を使ったりして、求めることができる。

------------------------------------------------------

cor(x, y, method = "spearman")

## [1] -0.7054505

Rx <- rank(x)

Ry <- rank(y)

plotn(Ry ~ Rx)#図4の描画

------------------------------------------------------

すると、外れ値の影響はかなり軽減される。

------------------------------------------------------

cov(Rx, Ry)/sqrt(var(Rx)*var(Ry))

## [1] -0.7054505

1 - 6*sum((Rx - Ry)^2)/(n * (n^2 - 1))

## [1] -0.7054505

12*sum((Rx - mean(Rx))*(Ry - mean(Ry)))/(n * (n^2 - 1))

## [1] -0.7054505

------------------------------------------------------

このとき、Sは以下の値である。

------------------------------------------------------

sum((Rx - Ry)^2)

## [1] 7666

------------------------------------------------------

このときのp値は、簡単に計算できるわけではないようだ。もし、漸近検定なら、Pearsonの相関係数と同様に、t分布を使ってp値を求められる。その場合は、exact = Fとする。

------------------------------------------------------

cor.test(x, y, method = "spearman", exact = F)

## 

##  Spearman's rank correlation rho

## 

## data:  x and y

## S = 7666, p-value = 1.338e-05

## alternative hypothesis: true rho is not equal to 0

## sample estimates:

##        rho 

## -0.7054505

pt(q = cor(x, y, method = "spearman")*sqrt(n-2)/sqrt(1 - cor(x, y, method = "spearman")^2), df = n - 2) * 2

## [1] 1.337728e-05

------------------------------------------------------

以下のように、サンプルサイズnのとき、Sは0からn(n^2 - 1)/3の値をとる。正確検定では、Sがn(n^2 - 1)/6より大きければ、Sより大きい値をとる確率の総和の2倍、Sが小さければSより小さい値をとる確率の総和の2倍をp値とする。

小標本であれば、以下のように正確なp値を手計算できる。

------------------------------------------------------

n <- 7

x <- round(runif(n, min = 10, max = 30), digits = 2)

y <- round(rcauchy(n, location = b1 * x + b0, scale = 3), digits = 2)

plotn(y ~ x)#図5の描画

Rx <- rank(x)

Ry <- rank(y)

plotn(Ry ~ Rx)#図6の描画

------------------------------------------------------

------------------------------------------------------

cor.test(x, y, method = "spearman", exact = T)

## 

##  Spearman's rank correlation rho

## 

## data:  x and y

## S = 106, p-value = 0.0123

## alternative hypothesis: true rho is not equal to 0

## sample estimates:

##        rho 

## -0.8928571

------------------------------------------------------

以下のように、gtoolsライブラリのpermutations関数を使い、すべての場合の順位を求める。そこから、すべての場合のSを求める。

------------------------------------------------------

head(permutations(n, n), 10)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7]

##  [1,]    1    2    3    4    5    6    7

##  [2,]    1    2    3    4    5    7    6

##  [3,]    1    2    3    4    6    5    7

##  [4,]    1    2    3    4    6    7    5

##  [5,]    1    2    3    4    7    5    6

##  [6,]    1    2    3    4    7    6    5

##  [7,]    1    2    3    5    4    6    7

##  [8,]    1    2    3    5    4    7    6

##  [9,]    1    2    3    5    6    4    7

## [10,]    1    2    3    5    6    7    4

vS <- apply(permutations(n, n), MARGIN = 1, FUN = function(x){

  sum((1:n - x)^2)

  })

histn(vS, xlab = "S", freq = F)#図7の描画

------------------------------------------------------

------------------------------------------------------

tS <- table(vS)

tS

## vS

##   0   2   4   6   8  10  12  14  16  18  20  22  24  26  28  30  32  34  36  38 

##   1   6  10  14  29  26  35  46  55  54  74  70  84  90  78  90 129 106 123 134 

##  40  42  44  46  48  50  52  54  56  58  60  62  64  66  68  70  72  74  76  78 

## 147  98 168 130 175 144 168 144 184 144 168 144 175 130 168  98 147 134 123 106 

##  80  82  84  86  88  90  92  94  96  98 100 102 104 106 108 110 112 

## 129  90  78  90  84  70  74  54  55  46  35  26  29  14  10   6   1

pS <- tS/length(vS)

S <- sum((Rx - Ry)^2)

S

## [1] 106

cols <- rep("grey", length = length(pS))

if(S > median(as.numeric(names(pS)))){

  cols[as.numeric(names(pS)) >= S] <- "red"

} else {

  cols[as.numeric(names(pS)) <= S] <- "red"

}

barplotn(as.vector(tS/length(vS)), ylab = "Proportion", col.fill = cols)#図8の描画

------------------------------------------------------

------------------------------------------------------

if(S > median(as.numeric(names(pS)))){

  sum(pS[as.numeric(names(pS)) >= S])*2

} else {

  sum(pS[as.numeric(names(pS)) <= S])*2

}

## [1] 0.01230159

------------------------------------------------------

確かにp値は一致した。サンプルサイズが大きいときに、どのようにして計算を省略しているかはわからないのだが、上記のような考えに基づいていることは確かなようだ。

相関と因果

　（工事中）

追記：同順があるときのSpearmanのρ

　同順があるときのρは以下のように求めることができる。

これは以下のように求まる。同順があるときの順位の分散は、Mann-WhitneyのU検定の時の導出を参考にしてほしい。

追記：同順があるときのKendallのτ

　同順があるときのKendallのτは以下のように計算する。

図的に解釈すると以下のようになる。