相関

相関係数

2008年9月24日に就任した中山成彬（なりあき）国土交通相は，「日教組（日本教職員組合）が強いところは学力が低い」などの発言で，5日後の28日には辞職してしまいます。

この「日教組が強いところは学力が低い」を反証するため，朝日新聞は27日朝刊で13道府県の全国学力調査の中3数学Aの点数をもとに「相関なし」と結論づけています。しかし，この13道府県はどういう基準で選んだのか，なぜ中3数学Aなのか，疑問の残る記事になってしまいました。紙面には「〈注〉科目は1位と最下位の県の得点差が最も大きいものを選びました」とありますが，恣意的に特定の科目を選んだととられないためには，全科目を調べるか，あるいは総合得点を使うべきでした。

同じ朝日新聞でも，版によっては13道府県の日教組組織率を数字で挙げています。文科省サイトで公開されている平成20年度全国学力・学習状況調査の都道府県別の全4科目正答率の合計と合わせて表にしておきます。組織率データの信頼性については，ここでは不問とします。

組織率（%）正答率合計
北海道50237.9
岩手県40238.8
秋田県50270.2
富山県50270.1
福井県90276.3
静岡県70259.2
愛知県60256.6
大阪府30231.4
香川県 1259.0
高知県10220.7
大分県60242.9
宮崎県10251.6
沖縄県40209.4

このデータの散布図（相関図）を描いてください。Rなら次のコマンドで描けます。

組織率 = c(50,40,50,50,90,70,60,30,1,10,60,10,40)
正答率合計 = c(237.9,238.8,270.2,270.1,276.3,259.2,256.6,231.4,259.0,220.7,242.9,251.6,209.4)
plot(組織率,正答率合計)

実際に描いてみると，なんとなく右肩上がりに見えるかもしれません。この傾向の度合，つまりこの場合は組織率と正答率合計の関連の度合を数字にしたものが，相関係数（correlation coefficient，長い名前はピアソンの積率相関係数 Pearson's product-moment correlation coefficient）です。相関係数は -1 から 1 までの値をとり，正負の関連が強いほど ±1 に近づき，関連が低ければ 0 に近づきます（詳しくは後で述べます）。

上のピアソンは Karl Pearson（1857〜1936年）で，ネイマン・ピアソンの理論で有名な E. S. Pearson の父です。

Rで相関係数を求める関数は cor() です。

> cor(組織率,正答率合計)
[1] 0.4251695

実は，相関係数と名の付くものはいくつかあり，Rでは後述のケンドールの順位相関係数，スピアマンの順位相関係数が同じ cor() で計算できます：

> cor(組織率,正答率合計,method="kendall")
[1] 0.3736324
> cor(組織率,正答率合計,method="spearman")
[1] 0.5076522

これらの有意性検定を行う関数は cor.test() です：

> cor.test(組織率,正答率合計)

	Pearson's product-moment correlation

data:  組織率 and 正答率合計 
t = 1.558, df = 11, p-value = 0.1475
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 -0.1643066  0.7908813 
sample estimates:
      cor 
0.4251695 

> cor.test(組織率,正答率合計,method="kendall")

	Kendall's rank correlation tau

data:  組織率 and 正答率合計 
z = 1.7298, p-value = 0.08366
alternative hypothesis: true tau is not equal to 0 
sample estimates:
      tau 
0.3736324 

Warning message:
In cor.test.default(組織率, 正答率合計, method = "kendall") :
   タイのため正確な p 値を計算することができません 
> cor.test(組織率,正答率合計,method="spearman")

	Spearman's rank correlation rho

data:  組織率 and 正答率合計 
S = 179.2146, p-value = 0.07656
alternative hypothesis: true rho is not equal to 0 
sample estimates:
      rho 
0.5076522 

Warning message:
In cor.test.default(組織率, 正答率合計, method = "spearman") :
   タイのため正確な p 値を計算することができません

このどれもが p ≦ 0.05 を満たしておらず，これだけのデータから何かを結論するのは早計ということになります。いずれにしても，この13道府県を選んだ理由が明らかでありませんので（たぶん朝日新聞が恣意的に選んだ？），仮に p ≦ 0.05 でも何の意味もないデータですが。

以下ではこれらの相関係数についてさらに詳しく説明します。

ピアソンの相関係数

互いに関連する（独立でない）二つの確率変数 X，Y を考えます。例えば X は数学の点数，Y は理科の点数だとすると，X が平均より大きいときは Y も平均より大きい傾向があり，X が平均より小さいときは Y も平均より小さい傾向がありそうです。このような二つの変数の間の関係を調べてみましょう。

平均よりどれくらい大きいか（小さいか）を調べるには，テストの点数なら偏差値に直すほうがわかりやすいでしょう。同じように，統計学でも，変数 X の平均値（期待値）を μX = E(X)，分散を σX2 = E((X - μX)2) とするとき，

\[ x = \frac{X - \mu_X}{\sigma_X} \]

で与えられる x に変換して考えると便利なことがあります。X を x に直すことを標準化する（standardize）ということにします。同様に Y を標準化したものを y とします。

こうしておくと，

\[ E(x) = 0, \quad E(x^2) = 1, \quad E(y) = 0, \quad E(y^2) = 1 \]

となります。

ここでもし X，Y が独立なら積の期待値は E(xy) = 0 ですが，一般には E(xy) は必ずしも 0 になりません：

この

\[ \rho = E(xy) = E\biggl(\frac{X - \mu_X}{\sigma_X} \cdot \frac{Y - \mu_Y}{\sigma_Y}\biggr) \]

が，X と Y の相関係数（ピアソンの積率相関係数）です。証明は省きますが，相関係数 ρ は必ず

\[ -1 \leqq \rho \leqq 1 \]

の範囲に入ります。

上の ρ は母集団 X，Y の相関係数でしたが，標本については

\begin{align*} r &= \frac{1}{n-1} \sum_{i=1}^n \frac{X_i - \bar{X}}{s_X} \frac{Y_i - \bar{Y}}{s_Y} \\ &= \frac{\displaystyle \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})(Y_i - \bar{Y})}{\displaystyle \sqrt{\frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2} \sqrt{\frac{1}{n-1} \sum_{i=1}^n (Y_i - \bar{Y})^2}} \end{align*}

で計算します。この最後の式では 1/(n - 1) はすべて約分して消すことができますので，「n で割るか n - 1 で割るか」の話はここでは影響しません。この r も必ず -1 ≦ r ≦ 1 の範囲に入ります。

数学的には，r は二つの n 次元の単位ベクトルの内積にほかならず，このことがわかれば -1 ≦ r ≦ 1 は自明です（（コーシー・）シュワルツの不等式）。

この最後の式の分子 $\frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})(Y_i - \bar{Y})$ を X と Y の共分散（covariance）と呼びます。

Rで例えば X1 = 1, X2 = 2, X3 = 3 と Y1 = 1, Y2 = 3, Y3 = 2 の相関係数を求めるには，次のようにします。

> x = c(1,2,3)    # x = 1:3 でも同じ
> y = c(1,3,2)
> cor(x,y)

これで r = 0.5 が求まります。

X，Y が独立に正規分布に従うなら，

\[ t = \frac{r \sqrt{n-2}}{\sqrt{1 - r^2}} \]

は自由度 n - 2 の t 分布に従います。

さきほどの例で計算すると，

> x = c(1,2,3)    # x = 1:3 でも同じ
> y = c(1,3,2)
> r = cor(x,y)    # r = 0.5 になる
> n = 3
> t = r * sqrt(n-2) / sqrt(1 - r^2)  # t = 0.5773503
> 2 * pt(-t,n-2)  # 0.6666667 と表示される

となります。より簡単には cor.test(x,y) と打ち込めば相関係数と p 値などが出力されます。

> cor.test(x,y)

	Pearson's product-moment correlation

data:  x and y 
t = 0.5774, df = 1, p-value = 0.6667
alternative hypothesis: true correlation is not equal to 0 
sample estimates:
cor 
0.5

相関係数 r の不偏推定量 G(r) は次の式で近似できます（Ingram Olkin and John W. Pratt, ``Unbiased Estimation of Certain Correlation Coefficients'', The Annals of Mathematical Statistics, Vol.29, No.1 (1958), pp.201-211.）。

\[ G(r) = \left( 1 + \frac{1 - r^2}{2(n - 3)} \right) r\]

この式の誤差は n ≧ 8 で 0.01 以下，n ≧ 18 で 0.001 以下とのことです。

順位相関係数

t 分布を使って求めたピアソンの積率相関係数の p 値は，二つの変数が正規分布をするときには正確ですが，一般の場合には必ずしも正しくありません。そのため，データの数値そのものではなく，その順位だけによる方法がいくつか考えられました。

一つは，単にデータの順位からピアソンの相関係数を求める方法です。この方法による相関係数を，スピアマンの順位相関係数またはスピアマンの ρ（ロー）（Spearman's rank correlation coefficient，Spearman's rho）といいます。

等しい値（タイ，tie）が現れるときは，その順位は，等しくなかったときの順位の平均値にします。たとえば，実際の値が 5, 7, 7, 9, 10 のとき，順位は 1, 2.5, 2.5, 4, 5 とします（あるいは大きい順に 5, 3.5, 3.5, 2, 1 としても同じことです）。

スピアマンの ρ の有意性検定は，n が十分大きければピアソンの相関係数と同じ方法（t 検定）で可能です。

現バージョンのRでは，帰無仮説を ρ = 0 としたときのスピアマンの ρ の有意性検定は，n ≦ 1290 のときは D. J. Best and D. E. Roberts, ``Algorithm AS 89: The Upper Tail Probabilities of Spearman's Rho,'' Applied Statistics, Vol.24, No.3 (1975), pp.377-379 にほぼ従って計算します。これは，n ≦ 9（元論文では6）のときは片方のデータの n! 通りの並べ替えを行って，そのうち何通りが，観測された |ρ| 以上の値をとるかを調べ，それを n! で割ったものを求めます。n > 9（元論文では n > 6）のときは 1/n についての級数展開で，少なくとも2桁の精度があるとのことです。n > 1290 のときは t 検定を使います。

もう一つのよく使われるものは，ケンドールの順位相関係数またはケンドールの τ（タウ）（Kendall's rank correlation coefficient，Kendall's tau）と呼ばれるものです。ケンドールが1938年に流行らせたのでこう呼ばれますが，ケンドールによれば，1900年ごろからいろいろな人が使っていたとのことです。

これの求めようとするものは，ランダムに二つを選んだとき，その X の順序と Y の順序が同じになる確率（例えば，A君がB君より数学ができたとき，A君がB君より英語もできる確率）です。実際には，この確率を2倍して1を引くことにより，-1 から 1 の範囲に収めたものが，ケンドールの τ です。スピアマンの ρ より具体的に求めようとしているものがはっきりしているのと，正規分布で近似しやすいことが特長です。

より具体的には，すべてのペアについて，両変数が同順のペアの数から逆順のペアの数を引き，ペアの総数 n(n - 1)/2 で割ったものがケンドールの τ です。タイがある場合も含めてもっと厳密に言うと，S = mx = my = 0 を初期値として，1 ≦ i < j ≦ n を満たすすべての整数のペア (i,j) について，

を行い，最後に

\[ \tau = \frac{S}{\sqrt{\vphantom{m} \smash{m_x m_y}}} \]

とすると τ が求まります。分子 S は同順のペア数から逆順のペア数を引いたもので，分母は基本的にはペアの総数 m = n(n - 1)/2 ですが，タイがある場合は各変数のタイでないペアの総数の相乗平均としています。

ケンドールの τ とスピアマンの ρ は，いずれも -1 から 1 までの値をとり，二つの変数の順序関係がまったく同順であれば 1，まったく逆順であれば -1 になるという点では同じです。両者に1対1の対応はありませんが，近似的に非線形な関係があり，中程度の相関では τ のほうが絶対値が小さくなります。どちらが統計的に有意になりやすいということはありません。τ のほうが正規分布に近いので扱いやすい反面，n が小さいと τ は ρ に比べて飛び飛びの値しか取らないことが目立ちます。

下の図は10ペアの乱数を何回も作ってケンドールの τ（横軸）とスピアマンの ρ（縦軸）を計算した結果の散布図です。

2変数が同じ順に並んでいても，タイの位置が異なれば，どちらの順位相関係数も 1 になりません。

ケンドールの τ のタイの処理の仕方はいくつか考えられますが，上で述べたものは τb と呼ばれる方法です。

現バージョンのRでは，帰無仮説を τ = 0 としたときのケンドールの τ の有意性検定は，n < 50 またはオプション exact=TRUE 指定時には，タイがなければ厳密な方法を使います。これ以外の場合は正規分布で近似する方法で行います。τ の分子 S が単なる ±1 の和であることを考えれば，n が大きければ中心極限定理により正規分布に近づくことが理解できます。具体的には母集団が独立のとき τ はタイがなければ正規分布 N(0, (4n+10)/9n(n-1)) に近づきます。

タイがある場合はRでは τ の分子 S が分散

\[ V(S) = \frac{n (n - 1) (2 n + 5) - V_x - V_y}{18} + \frac{T_x T_y}{2n(n - 1)} + \frac{U_x U_y}{9n(n - 1)(n - 2)} \]

の正規分布に近づくことを使っています。ここで

\[ T_k = \sum t_k (t_k - 1), \quad U_k = \sum t_k (t_k - 1) (t_k - 2), \quad V_k = \sum t_k (t_k - 1) (2 t_k + 5) \]

です（k = x, y で，tx, ty はそれぞれ X, Y の個々のタイの長さのベクトルです。例えば X = (1, 2, 3, 3, 4, 5, 5, 5) なら，tx = (2, 3) で，$\sum t_x (t_x - 1) = 2(2-1) + 3(3-1) = 8$ となります）。

タイのある場合も含めて正確な p 値を求めるには，並べ替え検定の考え方を使います。例えば1万回のシミュレーションでは

t = cor(X,Y,method="kendall")
a = replicate(10000,cor(X,sample(Y),method="kendall"))
mean(abs(a) >= abs(t)) # 両側確率

このページの最初に挙げた例では p = 0.08366 でしたが，100万回のシミュレーションでは p = 0.095 ほどです。ただし，気をつけなければならないのは，ケンドールの τ が飛び飛びの値をとることです。このため，値をコピー＆ペーストして使っても四捨五入のために正しくない結果を生じることがあります：

> t
[1] 0.3736324
> mean(abs(a) >= t)
[1] 0.095084
> mean(abs(a) >= 0.3736324)
[1] 0.07212
> mean(abs(a) >= 0.3736323)
[1] 0.095084

この場合は p = 0.095 のほうが正しいのですが，少しずらすと p = 0.072 になります。正規分布による近似 p = 0.08366 はこのほぼ真ん中の値になっています。ちなみに，100万個の中にユニークな値は67個しかありません：

> length(unique(a))
[1] 67

奥村晴彦

Last modified: 2008-12-14 08:46:59