スピアマンの順位相関係数とその検定

　スピアマンの順位相関係数の検定と、並びかえ検定について書いていきます。

1.スピアマンの順位相関係数

　スピアマンの順位相関係数は二つの変数の関連の強さをあらわす統計量です。次のような変数(x, y)のペアからなるn個の標本が得られたとき、

(x_1, y_1), (x_2, y_2), ...(x_i, y_i), ...(x_n, y_n)

　スピアマンの順位相関係数は以下の式から計算できます。

ρ = \frac{12}{n(n^2-1)}\sum_{i=1}^{n}(R_{xi}-R_xmean)(R_{yi}-R_ymean)

　「Rxi」はすべての「x」を大きさの順に並べた際の「xi」の順位で、「Rxmean」はその順位の平均値です。順位の平均はサンプル数「n」を使って以下のように計算できます。

R_xmean=\frac{n+1}{2}

　「Ry」についても同様です。Wikipediaの計算とは形が違いますが、結果は一緒です。

2.順位相関の検定

　スピアマンの順位相関係数は、二つの変数「xi」と「yi」が、それぞれに独立した分布に従うという仮説の下で、漸近的に正規分布に従うという統計量です。
　検定はこの正規分布と標本から計算される順位相関係数をもとに、二つの変数「xi」と「yi」がそれぞれに独立した分布に従う、という帰無仮説を一定の有意水準のもと棄却できるか、という問題になります。

　正規分布以外にも、順位相関係数を利用した以下のような検定統計量を利用する手もあります。

t=\frac{ρ}{\sqrt{\frac{1-ρ^2}{n-2}}}

　この検定統計量「t」は、サンプル数「n」が十分に大きいとき、上記の仮説「二つの変数が独立～」の下で、自由度「n-2」のスチューデントのt分布に近似的に従います。二つのパラメータを考えなければならない正規分布より楽ですね。

3.サンプル数に応じた帰無分布の導出

　スピアマンの順位相関係数のように、標本間の観測値の順位だけが問題になるような統計量は、同順がないとき、その統計量がある帰無仮説の下で従う分布（以下帰無分布と呼ぶことにします）を、サンプルの数に基づいて特定することができます。

　小さなサンプルを例に見ていくとわかりやすいと思います。サンプル数「n=3」の場合を考えていきます。

(x_1, y_1), (x_2, y_2), (x_3, y_3)

　このような標本が得られたとき、「xi」と「yi」の順位数「Rxi」と「Ryi」は(1, 2, 3)のいずれかになるはずです。
　スピアマンの順位相関係数は計算を見ればわかるように、サンプルのペアのちがいは問題にされていないので、同順を考えない時、順位数「Rxi」と「Ryi」の付き方は次の６パターンを考えれば十分です。

(R_{x1}=1, R_{y1}=1), (R_{x2}=1, R_{y2}=2), (R_{x3}=3, R_{y3}=3)

(R_{x1}=1, R_{y1}=1), (R_{x2}=1, R_{y2}=3), (R_{x3}=3, R_{y3}=2)

(R_{x1}=1, R_{y1}=2), (R_{x2}=1, R_{y2}=1), (R_{x3}=3, R_{y3}=3)

(R_{x1}=1, R_{y1}=2), (R_{x2}=1, R_{y2}=3), (R_{x3}=3, R_{y3}=1)

(R_{x1}=1, R_{y1}=3), (R_{x2}=1, R_{y2}=1), (R_{x3}=3, R_{y3}=2)

(R_{x1}=1, R_{y1}=3), (R_{x2}=1, R_{y2}=2), (R_{x3}=3, R_{y3}=1)

　ここでスピアマンの順位相関の検定の帰無仮説「二つの変数が独立～」が真であるとき、これらの順位の分布のパターンは、どれも「同様に確からしい」と言えます。特定の順位パターンの出現確率が高くなると言える合理的な根拠はどこにもありません。
　それぞれの順位パターンから計算されるスピアマンの順位相関係数は以下のようになります。

1.0, 0.5, 0.5, -0.5, -0.5, -1.0

　同順を考えなければ、これ以外の値はとりようがないので、サンプル数「n=3」の時、帰無仮説「二つの変数が独立～」の下でスピアマンの順位相関係数は以下の分布に従います。

ρ	-1.0	-0.5	0.5	1.0
出現確率	1/6	2/6	2/6	1/6
累積確率	1/6	6/3	5/6	6/6

　これはサンプル数が増えても同じことで、(xi, yi)に可能なすべての順位の組み合わせから相関係数を計算し、その値をカウントしてやれば分布が得られます。次のようなスクリプトを使えば、帰無分布を計算できます。

import itertools

def spearmanr_null_dis(sample_size):
    mean = (sample_size + 1) / 2
    cnst =  12 / (sample_size * ((sample_size * sample_size) - 1))
    null_dis = []
    for i in itertools.permutations(list(range(1, sample_size+1))):
        sum = 0
        for j, k in enumerate(i):
            sum += (j + 1 - mean) * (k - mean)
        null_dis.append(sum * cnst)
    return null_dis

　組み合わせの総数は、サンプル数「n」の階乗個になります。「n=3」から「n=8」までの帰無分布は大体以下のようになります。

　「n」が増えていくのに従って、きれいな左右対称の山型に近づいていくのがわかります。

4.scipy.stats.permutation_test()

　上のスクリプトの難点は、サンプル数の増大にともなって帰無分布の要素数が跳ね上がっていくことです。「n=10」の時点で「10!=3,628,000」、ここからサンプルが一つ増えるごとに最低でも桁が一つずつ上がっていきます。
　現実的な範囲での近似分布を与えてくれる関数がscipy.statsのpermutation_test()関数です。

scipy.stats.permutation_test(data, statistic, *, permutation_type='independent', vectorized=None, n_resamples=9999, batch=None, alternative='two-sided', axis=0, random_state=None)

　permutation_test()はノンパラメトリック検定の一種である並び替え検定を行います。この関数から帰無分布を得るには、具体的なデータを与えてやる必要があるので次のような適当なデータを用意します。

import scipy.stats as stats

x = stats.norm.rvs(loc=5, size=10)
y = stats.norm.rvs(size=10)

　さらに検定統計量を計算する関数を変数：statisticに渡す必要があるので、データを受け取ってスピアマンの順位相関係数を計算する関数を作成します。

def statistic(x, y):
    res, _ = stats.spearmanr(x, y)
    return res

　stats.pearsonr()は相関係数と正規近似した場合のp値を返しますが、今回はp値は必要ないので無視しています。
　材料をpermutation_test()に渡すと次のような帰無分布が得られます。比較として同じサンプルサイズで前述のspearmanr_null_dis()で作った帰無分布ものせておきます。

ref = stats.permutation_test((x, y), statistic, permutation_type='pairings')
plt.subplot(121)
plt.hist(ref.null_distribution)
plt.xlabel('permutation_test()')

null_dis = spearmanr_null_dis(10)
plt.subplot(122)
plt.hist(null_dis)
plt.xlabel('spearmanr_null_dis()')

plt.show()

　近似としては微妙な感じですが、度数に注目すれば、まあ扱いやすいのではないだろうかと思えます。

　permutation_test()はサンプル数が多く帰無分布が巨大になってしまう場合、自動的にランダムサンプリングによる近似を返すようになります。そのさい帰無分布はキーワード引数：n_ressamplesで指定した要素数をもつリストとして得られます。
　この近似はサンプル数「n」がn! ≧ n_ressamplesとなるとき行われます。n_ressamplesの初期値は9999なので、「n≧8」であればランダムサンプリングに自動的に切り替わります。

　n! < n_ressamplesであれば、scipyのリファレンスによれば、正確検定が行われるのですが、こちらもやや癖があります。
　permutation_test()の行う並べかえ検定は、3.で見たようなサンプルの並べかえから帰無分布を導出するのですが、並べかえ検定ではこの並び替えがさらに徹底したものになります。

(x_1, y_1), (x_2, y_2), (x_3, y_3)

　このような標本が得られたとき、3.で行った並び替えは以下のように「yi」の並びだけを問題にしたもので事足りました。

(x_1, y_1), (x_2, y_2), (x_3, y_3)

(x_1, y_1), (x_2, y_3), (x_3, y_2)

...

　しかし並び替え検定では以下のように、「xi」と「yi」を入れ替えるパターンも問題にします。

(x_1, y_1), (x_2, y_2), (x_3, y_3)

...

(y_3, x_3), (y_2, x_2), (y_1, x_1)

　並び替え検定はこのすべてのパターンについて検定統計量を計算して帰無分布を導出するので、その分布はn!×2^n個のリストとして得られます。

　permutation_test()は正確検定を行う場合でもn_ressamples個のリストで帰無分布を返します。初期値の場合は要素数9999のリストです。
　これはサンプル数「n=5」のように、正確検定の帰無分布の表現に9999個も必要ない（5!×2^5=3840）時であっても、帰無分布として9999のリストを返してきます。
　9999に足りない分は、おそらく同じ計算を繰り返すことで埋めてると思うのですが、正確な分布を得るには5!×2^5=3840の定数倍の数が必要で、9999個のリストでは分布が若干ゆがみます。

　まぁサンプル数的に現実的になにか問題になるとは思えませんが。

5.おわりに

　ずいぶん長くなってしまったので、この辺で〆ようと思います。本当は正規近似との比較や、もう少しくわしく並べかえ検定について書ければと思っていたのですが、前準備だけがだらだらと続く結果になってしまいました。
　並び替え検定については、とても分かりやすい資料を獨協大学の先生が作ってくださっているので、詳しく知りたい方は以下の参考のリンクから当たってみてください。

6.参考

朝倉書店　統計解析スタンダード　ノンパラメトリック法　村上 秀俊(著)
https://www.asakura.co.jp/detail.php?book_code=12852
獨協大学学術リポジトリ　パーミュテーション検定について　 松井 敬（著）
https://dokkyo.repo.nii.ac.jp/records/312
Wikipedia(英)　Permutation test
https://en.wikipedia.org/wiki/Permutation_test
（日本のWikipediaには並び替え検定（パーミュテーション検定）の記事はないようです）