せいびし

専門性を持たない自分の正体を求めて

【各種検定の理論】2標本の検定・位置母数の推定(ノンパラメトリック検定)

導入

本記事では、ノンパラメトリック検定における2標本の検定について取り扱いたいと思います。パラメトリック検定で言う所の、対応なしt検定に相当します。今回はウィルコクソン順位和検定とマン・ホイットニーのU検定を紹介します。後で述べていきますが、ウィルコクソン順位和検定とマン・ホイットニーのU検定は実質的に同じことをしており、ウィルコクソン順位和検定の検定統計量からマン・ホイットニーの検定統計量に数理的に繋げることができます。ここで次のような疑問が生まれます。本当に何もかも完全に一致するのか。僕は初めて同じである事を知った時に本当か?と疑問に思いました。だって同じことをするなら、どっちかの存在意義が無くなりませんか?なぜ2つ残っているんでしょうかね?気になったら即手を動かして遊んでみましょう。


(※追記:本当はどちらもぴったり一致する事を確かめようと思ったのですが、マン・ホイットニーのU検定が意味分からんなので、主にマン・ホイットニーのU検定の話をしていきます…勉強不足ですみません…とりあえずやろうとしていたことはメモとして残しておきます。できなかったところは線分を引っ張っておきます。)

(※追記2:後で示しますが、期待値の謎問題については現在進行形で検証中です。深い理論の話になってしまうので医療従事者のための統計学には書きませんが、また改めて数理統計学の枠を設けて遊んでいきたいと思います。理論が完全に明らかになり次第、しれっと修正しておきます(2019年11月9日 追記))

(※追記3:期待値の謎問題は解決しました。帰無仮説 \displaystyle F_1(x)=F_2(x)の下での期待値である事をすっかり忘れていました。従って、2つ山ができるのは、単にwilcox.test(data1, data2)とwilcox.test(data2, data1)で順番を変えただけであってそりゃそうという結論になりました。ちなみに、このブログで述べているMann-WhitneyやMoodの検定などは、帰無仮説の下での検定統計量の分布は左右対称になっています。従って、例えばwilcox.test(data1, data2)とwilcox.test(data2, data1)で算出される検定統計量を足してみると期待値の2倍になる事が確認できます。不勉強のために一人で騒いでしまいました…とりあえず勉強の後を残すために、一人で騒いでいる様子はそのまま残しておきたいと思います(2019年12月10日 追記))

1.ウィルコクソン順位和検定の概要

しくみ

1標本の時と同じように、初めに仕組みの話をします。まずは簡単に前提条件と文字の定義をしておきましょう。
 \displaystyle X_1 = (X_{11}, ..., X_{1m})を分布関数 \displaystyle F_1(x)からの確率標本、 \displaystyle X_2 = (X_{21}, ..., X_{2n})を分布関数 \displaystyle F_2(x)からの確率標本とします。ここで、分布関数 \displaystyle F_1(x), F_2(x)は連続である事以外には何も分からないものとします。次の図を見てみましょう。

f:id:seibibibi:20191107191020p:plain

母集団の分布関数は点線で、データは丸で表しています。データを見て、左の状態と右の状態で、どちらがデータに差が表れているように見えるでしょうか。例えば、この2群のデータに対して順位をつけてみましょう。左側の状態では各群で順位がばらばらであり、右の状態では全体的に青の分布の方が順位が高い事が分かるのではないでしょうか。これが本記事で述べる検定の仕組みです。もし、2群間でデータに差が表れているなら、どちらかの群のデータはもう一方の群のデータよりも順位が高くなるはずです。では、この仕組みをどうやって実現しているのでしょうか。

仮説と検定統計量

仮説

まず初めに仮説を示します。

 \displaystyle H_0:F_1(x)=F_2(x)

 \displaystyle H_1:F_2(x)=F1(x-\Delta)

ここで、謎の対立仮説が現れました。この仮説を図で表すなら、以下のようになります。
f:id:seibibibi:20191107191231p:plain

 \displaystyle \Deltaは移動量のようなものであり、位置母数と呼ばれています。つまり、一方の分布を \displaystyle \Deltaだけ移動させると完全に重なります。この事から、ウィルコクソン順位和検定は、データに差が表れているなら分布関数が大きくシフトしているはず(差が大きい)という事を検証します。

検定統計量

次に、検定統計量を述べていきます。しくみの所でも述べたように、基本的には全てのデータを1つにまとめた後で順位で並べ替えて、各群の順位を見るという処理を行います。式を書く前にイメージを図で示します。

f:id:seibibibi:20191107200523p:plain

ごくシンプルに、順位と割り当てた数の積和になっています。この事を踏まえて検定統計量を眺めてみましょう。検定統計量は次の式で表されます。
 \displaystyle W = \sum_{i=1}^{m}R_i

ここで、 \displaystyle R_iは群1に所属するデータの順位です。0を割り当てた物は総和しても仕方が無いので、1を割り当てた物だけを総和しています。とてもシンプルですよね。

2.マン・ホイットニーのU検定

ウィルコクソン順位和検定との違い

基本的に、ウィルコクソン順位和検定と同じ仕組み、同じ仮説を設定しています。従って、特にこれと言って説明する事は無いのですが、強いて言うなら検定統計量の導出が少し異なっています。

検定統計量

初めに図を示します。

f:id:seibibibi:20191107200412p:plain

マン・ホイットニーのU検定の場合は、網羅的にデータの大きさを比較します。群1のデータの方が大きければ1を割り当て、そうでなければ0を割り当てます。最終的に割り当てた数字を全て総和してあげれば、マン・ホイットニーのU検定の検定統計量が導出されます。式は以下の通りです。
 \displaystyle MW = \sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)

ここで、 \displaystyle U(i,j)は、 \displaystyle X_{1i} \geq X_{2j}の時に1、それ以外が0となる量です。要するに上で示した処理を行っています。

何が同じ

マン・ホイットニーのU検定の検定統計量はウィルコクソン順位和検定の検定統計量から導出できます。具体的には、

 \displaystyle W=\sum_{i=1}^{m}R_i = \sum_{i=1}^{m} (\sum_{j=1}^nU(i,j)+i) = \sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)+\frac{ m(m+1) }{2}

イメージ的には、ぜんぶ1つにまとめて並べた時に、母集団1の任意のサンプルから次のサンプルまでに入っている母集団2のサンプルの個数を足し合わせている感じでしょうか。図で言えば以下のようになります。


f:id:seibibibi:20191108113648p:plain


やってることが同じな気がしてきましたね、なんだかとてもつまらない予感がしてきました…

3.シミュレーション

やること

さて、ここからシミュレーションに入ります。導入でも述べた通り、確かめたい事はウィルコクソン順位和検定とマン・ホイットニーのU検定は完全に一致するのかのはずでした。以下にやろうとしていた実験の手順を示します。

f:id:seibibibi:20191107202142p:plain

母集団分布を2つ用意し、各群でデータをサンプルします。サンプルした後で各検定の検定統計量を求めます。この手順を複数回繰り返し、検定統計量に関するヒストグラムを作成します。ここで、今回完全に一致するのかどうかは検定統計量を正規近似した下で行う事とします*1。どちらの検定統計量も、サンプル数が大きい時には正規近似できることが知られています。従って、この2つの検定統計量を標準化し、標準正規分布上で比較したいと思います。もし本当に何もかも一致するなら、ヒストグラムは完全に一致するはずです。

ここで、用いる母集団分布を示しておきます。今までは分布関数について話をしてきましたが、以下で示す分布は確率密度関数である事に注意です。

f:id:seibibibi:20191107203845p:plain

2つの分布*2から、サンプルするデータの個数は赤の分布から100個、青の分布から50個、標本抽出は10000回行います。

正規近似のシミュレーション

理論的背景

1標本のノンパラ検定でも述べましたが、この検定においても正規近似を適用する事が出来ます。ある検定の検定統計量を \displaystyle Uとした時、標準化後の検定統計量 \displaystyle Zは次の式で正規近似できます。

 \displaystyle Z = \frac{U-E[U]}{\sqrt{V[U]}}~N(0, 1)

ここで、それぞれの検定量の期待値と分散は次の式で示されます。
 \displaystyle E[W]=\frac{m(m+n+1)}{2}, V[W]=\frac{mn(m+n+1)}{12}

 \displaystyle E[MW]=\frac{mn}{2}, V[MW]=\frac{mn(m+n+1)}{12}

という事で、この期待値と分散により正規近似していきます。

実験結果1 自作した検定統計量とパッケージとの対応

「R Mann Whitney」で調べてもwilcox.testを使っている物ばかりだし、「R Mann Whitney documentation」で調べても意味不明な物が出てきて萎えたので、それなら一から作ってしまおうという事で作りました。以下に作成したコードの1部分を載せます。

# 検定統計量の計算
calc <- function(data1, data2){
	# 両方向の比較を書いてみる(data1→data2とdata2→data1)
	# 1. data1→data2
	# Mann-Whitney
	MW1 <- 0
	d1 <- sort(data1)
	d2 <- sort(data2)
	for (i in 1:length(d1)){
		MW1 = MW1 + length(d2[d2 < d1[i]])
	}
	# たしかめ
	wilcox.test(data1, data2)
	# 2. data2→data1
	MW2 <- 0
	for (i in 1:length(d2)){
		MW2 = MW2 + length(d1[d1 < d2[i]])
	}
	MW2
	# たしかめ
	wilcox.test(data2, data1)
	return (c(MW1, MW2))
}

自分で作成したは良いものの、合ってるかどうかが肝心なので、たしかめの部分を実行するとこんな感じになりました。

> # Mann-Whitney
> MW1 <- 0
> d1 <- sort(data1)
> d2 <- sort(data2)
> for (i in 1:length(d1)){
+ MW1 = MW1 + length(d2[d2 < d1[i]])
+ }
> # たしかめ
> wilcox.test(data1, data2)

        Wilcoxon rank sum test with continuity correction

data:  data1 and data2
W = 1199, p-value = 2.163e-07
alternative hypothesis: true location shift is not equal to 0

> MW1
[1] 1199
> 

良い感じに一致していました。で、作った後で気付いたのですが、wilcox.testのdocumentationに次のような記述がありました。

R's value can also be computed as the number of all pairs (x[i], y[j]) for which y[j] is not greater than x[i], the most common definition of the Mann-Whitney test.

wilcox.test function | R Documentation
おう、まじか…ちゃんと調べれば良かった…と思いました(萎えないでちゃんと調べよう)。でも、逆に言えばウィルコクソンの検定統計量が計算されている訳では無いという事です。マン・ホイットニーの検定統計量とウィルコクソンの検定統計量の関係式を用いて計算したかったのですが、次に示す謎により計算できませんでした。

実験結果2 謎

とりあえずマン・ホイットニーのU検定の検定統計量のヒストグラムを書いてみます。以下の様になりました。

f:id:seibibibi:20191108160552p:plain

赤が群1→群2で比較した検定統計量、青が群2→群1で比較した検定統計量です。上で示したマン・ホイットニーの検定統計量の期待値の式に \displaystyle m=100, n=50を入れると、2500になります。しかし、どちらの山も2500にはならず、なんなら2つの山の真ん中が2500なのかなという感じがします。期待値が山のピークになる事は必ずしもそうではない(指数分布とか連続一様分布と考えると)ので、2つの山の真ん中が期待値と言われても別に違和感がない(むしろ期待値の式に突っ込めば真ん中になるのでしょう)と言えばないですが、うーん…わからない。両方考えなきゃなのかなあ、定義は片方だけなのに。こんな結果なので、ウィルコクソンの検定統計量を計算しようにも、 \displaystyle nを入れるのか \displaystyle mを入れるのか分からないので、残念ながら比較はできませんでした。

考察

つらい

まとめ

今回は、ウィルコクソン順位和検定とマン・ホイットニーのU検定で遊びました。本当にどっちも同じなのか?という問いに対しては、これも後から気付いたのですが、以下の式で証明できそうです。
 \displaystyle E[W] = E[\sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)+\frac{ m(m+1) }{2}] = E[MW+\frac{ m(m+1) }{2}]
 \displaystyle = \frac{mn}{2}+\frac{ m(m+1) }{2} = \frac{m(m+n+1)}{2} = E[MW]
 \displaystyle V[W] = V[\sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)+\frac{ m(m+1) }{2}] = V[MW]
同じじゃんか…

後から気付いて証明っぽいのもできてしまった上に、検定統計量の振る舞いも良く分からないし……なんもわからんです…がんばりたい…

*1:サンプル数が少ない時の帰無分布を求めるのがとても面倒くさいので、サンプル数が少ない時の検証は魔が差した時にします。

*2:赤、青ともに混合比が0.3と0.7の混合正規分布です。赤は \displaystyle N(10, 5), N(20, 5)の2つ、青は \displaystyle N(15, 5), N(50, 5)の2つを混ぜています。