【統計】Wilcoxonの順位和検定とマンホイットニーのU検定

Wilcoxonの順位和検定（Wilcoxon rank sum test）とマンホイットニーのU検定（Mann–Whitney U test）について。

【目次】

帰無仮説、対立仮設
計算式等
計算例
プログラムコード
W統計量とU統計量の関係
参考

両方とも２標本の分布の位置のズレを評価するノンパラメトリックな手法。
これらは等価な手法であり、同じ結果が得られる。

一方、パラメトリックな手法としては t検定が代表的。

cochineal19.hatenablog.com

帰無仮説、対立仮設

帰無仮説 $H_{0}:$ ２標本の確率分布に差がない（ズレていない）。
対立仮設 $H_{1}:$ ２標本の確率分布に差がある（ズレている）。

※２標本は同じ形状の確率分布であるという前提で、分布の位置のズレを見ている。正規分布に依存しないのが特徴。

計算式等

Wilcoxonの順位和検定とマンホイットニーのU検定は等価であり、同じ結果が得られるが、統計量と期待値に違いがある。

順位和

統計量の算出前に、先ず、比較したい２標本の観測値に全体順位（Rank）をつけて、標本ごとの合計（順位和）を算出する。観測値にタイ（同一順位）がある場合は、順位の平均値をつける。
※以降、本記事では標本をA、B、順位和を $R_{A}$ 、 $R_{B}$ とする。

Wilcoxonの順位和検定

$\quad W = min \left( R_{A},\ R_{B} \right)$

$\quad E[ W] =\dfrac{n_{A}\left( n_{A}+n_{B}+1\right) }{2}\qquad※R_{B}\verb|の方が小さければ |\ \dfrac{n_{B}\left( n_{B}+n_{A}+1\right) }{2}$

$\quad Var[ W] =\dfrac{n_{A}n_{B}\left( n_{A}+n_{B}+1\right) }{12} - \dfrac{n_{A}n_{B}}{12\left( n_{A}+n_{B}\right) \left( n_{A}+n_{B}-1\right) } \sum \left( t^{3}_{i}-t_{i}\right)$

$\quad z = \dfrac{W-E[ W] }{\sqrt{Var [ W] }}$

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

$\quad U = min \left( R_{A}-\dfrac{n_{A}\left( n_{A}+1\right) }{2} ,\quad R_{B}-\dfrac{n_{B}\left( n_{B}+1\right) }{2} \right)$

$\quad E[ U] =\dfrac{n_{A} n_{B}}{2}$

$\quad Var[ U] =\dfrac{n_{A}n_{B}\left( n_{A}+n_{B}+1\right) }{12} - \dfrac{n_{A}n_{B}}{12\left( n_{A}+n_{B}\right) \left( n_{A}+n_{B}-1\right) } \sum \left( t^{3}_{i}-t_{i}\right)$

$\quad z = \dfrac{U-E[ U] }{\sqrt{Var [ U] }}$

Wilcoxonの順位和検定もマンホイットニーのU検定も分散は同じ。
分散の第２項の $\dfrac{n_{A}n_{B}}{12\left( n_{A}+n_{B}\right) \left( n_{A}+n_{B}-1\right) } \sum \left( t^{3}_{i}-t_{i}\right)$ はタイ（同一順位）の補正項であり、式中の t が同一値の個数を示す。
例えば、{10,10,11,12,12,12} のデータがあるとき、10が２つ、11が１つ、12が３つになるため、それぞれ $t=2、t=1、t=3$ となり、 $t^{3}-t$ の計算で $2^{3}-2=6、1^{1}-1=0、3^{3}-3=24$ になる。
タイがない数値は解が０になり、タイにだけ働くことが分かる（タイがまったくない場合は第２項の解全体が０）。

計算例

今回は次の架空データAとBを比較してみる。

　A {4, 9, 10, 11, 12, 13, 13, 15, 18, 18, 20}
　B {16, 16, 17, 19, 19, 22, 22, 23, 25}

Excelでの計算。

f:id:cochineal19:20201215223427p:plain — Wilcoxon-Mann–Whitney test

Excelの簡単な説明。
・A列：標本のグループ名
・B列：観測値。小さい順にならべている。
・C列：同値の数を上から数ている（セルC4：=IF(B3<>B4,1,D3+1)）。
・D列：同値の最終行のみ表示している（セルD4：=IF(B4=B5,"",C4)）。
・E列：分散の補正項用（セルE4：=IF(D4<>"",D4^3-D4,"")）
・G列：標本Aの順位（セルG4：=IF($A4=G$3, RANK.AVG($B4,$B:$B,1), "")）
・H列：標本Bの順位（セルH4：=IF($A4=H$3, RANK.AVG($B4,$B:$B,1), "")）

Wilcoxonの順位検定の場合

W：順位和 $R_{A}=77,\ R_{B}=133$ だったため、 $R_{A}=77$ （小さい方）をWとする。
期待値： $R_{A}$ の方が小さいため、次の式で計算する。
$E[ W] =\dfrac{n_{A}\left( n_{A}+n_{B}+1\right) }{2} =\dfrac{11\left( 11+9+1\right) }{2}=115.5$
分散：タイがあるため補正（E列およびK15参照）した値を導出する。
$Var[ W] = 172.59868$
z値：z= (77-115.5)/sqrt(172.6) = -38.5/sqrt(172.6) = -2.93050。
p値：標準正規分布を使って p=0.00338。
（セルN23：=2*(1-NORM.S.DIST(ABS(N22),TRUE)) ）

マンホイットニーのU検定の場合

U： $U_{A}=11,\ U_{B}=88$ だったため、 $U_{A}=11$ （小さい方）をUとする。
期待値：計算式のとおり (11×9)/2 = 49.5。
分散：Wilcoxonと同じ。タイを補正した値を導出する。 $Var[ W] = 172.59868$
z値：z= (11-49.5)/sqrt(172.6) = -38.5/sqrt(172.6) = -2.93050。
p値：標準正規分布を使って p=0.00338。
（セルK23：=2*(1-NORM.S.DIST(ABS(K22),TRUE))）

z値の導出過程で、Wilcoxonもマンホイットニーも「-38.5/sqrt(172.6)」となり、等価になることが分る。どちらの検定も今回の架空データでは２標本に差があり（p<0.01）。

Rでの実行：

> ads <- data.frame(
+   TRT01P <- c("A","A","A","A","A","A","A","A","B","B","B","A","A","B","B","A","B","B","B","B")
+   ,AVAL  <- c(4,9,10,11,12,13,13,15,16,16,17,18,18,19,19,20,22,22,23,25)
+ )
> ads$TRT01PF <- factor(ads$TRT01P)
> wilcox_test(AVAL ~ TRT01PF, data=ads, method="exact")

	Asymptotic Wilcoxon-Mann-Whitney Test

data:  AVAL by TRT01PF (A, B)
Z = -2.9305, p-value = 0.003384
alternative hypothesis: true mu is not equal to 0

Rのデフォルトで使えるwilcoxon.testは値が全然一致しない。詳細不明。
coinパッケージのwilcoxon_testを使った方が良い。

SASでの実行：

data ads;
  do AVAL = 4,9,10,11,12,13,13,15,18,18,20; TRT01PN=1; output; end;
  do AVAL = 16,16,17,19,19,22,22,23,25; TRT01PN=2; output; end;
run;
proc sort data=ads: by AVAL TRT01PN; run;
proc npar1way data=ads wilcoxon correct=no /*今回は連続修正なし*/;
  class TRT01PN;
  var   AVAL;
run;

f:id:cochineal19:20201219133544p:plain

z値の符号が逆転しているが、Pr > |z| なのでp値は同じ。

SAS社のサイトに詳細がないので不明だが、
T近似のp値は、おそらくz統計量と自由度N-1で算出している。

= T.DIST.2T(2.9305, 20-1) ≒ 0.008583004

プログラムコード

■ Rのコード

library(coin)
wilcox_test(AVAL ~ TRT01PF, data=ads, method="exact")

■ SASのコード

proc npar1way data=[InputDS] wilcoxon correct=no /*連続修正無の場合correct=noを追加*/;
 class TRT01PN;
 var AVAL;
 output out = [OutDS] wilcoxon; /* Z_WIL, P2_WILを取得。*/
run;

■ Pythonのコード

整備中

W統計量とU統計量の関係

W統計量とU統計量には次の関係があり等価。

$\qquad W-\dfrac{n_{A}\left( n_{A}+1\right) }{2}=U$

　<再掲>
　Wilcoxonの順位和検定
$\qquad W = min \left( R_{A},\ R_{B} \right)$
　マンホイットニーのU検定
$\qquad U = min \left( R_{A}-\dfrac{n_{A}\left( n_{A}+1\right) }{2} ,\quad R_{B}-\dfrac{n_{B}\left( n_{B}+1\right) }{2} \right)$