Wilcoxon

はじめに

Wilcoxonの符合順位和検定は、独立２群の差を順位の和から判断する方法である。
順位は離散変量であり数列を学んだ高校生であれば、期待値や分散などの統計量を数理的に導出することもできる。これから高校でデータサイエンスが学ばれるようになれば、高校生が統計的検定を学ぶ格好の題材になる。
このような考えのもと、Wilcoxonの符合順位和検定を高生校向け教材として考えてみたいと思う。
rのプログラムによるシミュレーションも加えて、理解の補助としたい。

ノンパラメトリック検定

独立２群に差がないか検定する方法として「２標本t検定」があった。
この検定法は母集団が正規分布に従っていることを前提としている。（母集団の正規性の仮定）
母集団の正規性を前提とする検定法をパラメトリック検定といい、正規性を前提としない検定をノンパラメトリック検定という。
ここで扱うWilcoxonの符合順位和検定はデータを順位に変換し、順位をデータとして扱う。
例えば、2つのクラスの成績を比較する場合、平均点を用いると外れ値の影響を受けやすい。
データを順位にしてしまえば、近接するデータの差は１であり、外れ値の影響を抑えることができる。
具体的にみてみよう。
10人の順位が２クラスA,Bにランダムに散らばっているとする。
rで離散型一様分布の乱数を生成する runif() を用いて１０個の乱数を作り、順位関数rankを適用する。
10個の順位をグループ x（４個）,グループy（６個） に分ける。

r <- rank(runif(10) )
x <- r[1:4];y <-r[5:10]

print(x)

## [1]  3  1 10  4

print(y)

## [1] 2 8 9 7 5 6

xの順位和\(W_x\) は確率変数で、期待値と分散は
　　\(E(W_x)=\dfrac{1}{2}\cdot 4\cdot11=22\) 　　\(V(W_x)=\dfrac{1}{12}\cdot4\cdot6\cdot11=22\)

導出は以下の通り。（数列の和の公式を用いる）

順位の期待値と分散

n 個から１個選んだものの順位をRとすると、Rは一様分布に従う確率変数で、
　　　\(P(R=i)=\dfrac{1}{n}\)　から、　\(\displaystyle E(R)=\sum_{i=1}^n i\cdot \dfrac{1}{n}=\dfrac{1}{n}\cdot\dfrac{n(n+1)}{2}=\dfrac{n+1}{2}\)
　　　\(\displaystyle E(R^2)=\sum_{i=1}^n i^2\cdot \dfrac{1}{n}=\dfrac{1}{n}\cdot\dfrac{n(n+1)(2n+1)}{6}=\dfrac{(n+1)(2n+1)}{6}\)
　　　\(\therefore V(R)=\dfrac{(n+1)(2n+1)}{6}-\left(\dfrac{n+1}{2}\right)^2=\dfrac{(n+1)(n-1)}{12}\)
このように、数列の和の公式を用いて計算ができる。

Wilcoxon順位和検定の統計量

２群 \(X_1=\left(x_{11},\cdots,x_{1n_1}\right)\:\:,\:\: X_2=\left(x_{21},\cdots,x_{2n_2}\right)\) について
\((n_1+n_2)\) 個の中での順位Rを考える。
\(X_1\:,\:X_2\) の要素の順位を、\(R_{X_1}=\left(R_{11},\cdots,R_{1n_1}\right)\:\:,\:\: R_{X_2}=\left(R_{21},\cdots,R_{2n_2}\right)\) とすると、
これらは順位 \(1,\cdots ,(n_1+n_2)\) からの \(n_1,n_2\) 個の標本である。（独立ではない）
\(X_1\:,\:X_2\) の順位和を、\(\displaystyle W_1=\sum_{i=1}^{n_1} R_{1i}\:\:,\:\:W_2=\sum_{i=1}^{n_2} R_{2i}\) とすると期待値は、
　　　\(\displaystyle E(W_1)=\sum_{i=1}^{n_1} E(R_{1i})=\dfrac{n_1(n_1+n_2+1)}{2}\:\:,\:\:E(W_2)=\sum_{i=1}^{n_2} E(R_{2i})=\dfrac{n_2(n_1+n_2+1)}{2}\)

次に分散であるが、順位は独立ではないので簡単には行かない。結果のみ示すと、
　　　\(V(W_1)=V(W_2)=\dfrac{n_1n_2(n_1+n_2+1)}{12}\)
詳細については当サイト Wilcoxon を参照のこと。

さて、繰り返しプログラムで\(W_x\) を1000回出力し、期待値と分散を確かめてみよう。

s = c()
for (i in 1:1000) {
  r <- rank(runif(10) )#離散型一様分布の乱数
  x <- r[1:4];y <-r[5:10]
  s <-c(s,sum(x))
}
hist(s)

mean(s)

## [1] 21.807

var(s)

## [1] 20.88263

理論値とヒストグラムは、ほぼ合致しているようである。

このように正規分布に近い曲線になることは、\(W_x\) の値を用いて２群に差があるかどうかを検定できることを示唆している。
P値は組み合わせ計算によって求める。
10個の順位を、グループ x（４個）,グループy（６個）に分ける方法は、
　　　\({}_{10} C_4=210\) 通り。

データの順位が
　　x=(1,2,3,4)　y=(5,6,7,8,9,10)　とすると、 \(W_x=10\)
　　\(P(W_x\leq 10)=P(W_x=10)=1/210=0.0047\)

このP値は有意であり、順位が偏っていることがわかる。
Rでは次のように計算される。

library(exactRankTests)

##  Package 'exactRankTests' is no longer under development.
##  Please consider using package 'coin' instead.

x <- c(1,2,3,4)
y <- c(5,6,7,8,9,10)
#片側検定
wilcox.exact(x=x,y=y,paired=F,alternative = "less")

## 
##  Exact Wilcoxon rank sum test
## 
## data:  x and y
## W = 0, p-value = 0.004762
## alternative hypothesis: true mu is less than 0

組合せ計算から求めた値と一致している。

マン・ホイットニー

Wilcoxonの順位和検定とマンホイットニーのU検定は同等であると言われている。
検定量Uは両群で
　　　U= \(\displaystyle \sum_i\)( \(x_i\) より大きい他方の群のデータ数 )
rを用いてデータ数を数えてみよう。

r <- runif(10)
x <- r[1:4];y <-r[5:10]

x

## [1] 0.61418627 0.50932384 0.79393178 0.05221353

y

## [1] 0.2066053 0.9437817 0.1714336 0.7796532 0.5673343 0.5125095

sum(y >= x[1])　#x[1]より大きいyの要素の個数を数える

## [1] 2

M = c()
for (i in 1:1000) {
  r <- rank(runif(10) )#離散型一様分布の乱数
  x <- r[1:4];y <-r[5:10]
  s = 0
  for (j in 1:4){
    s <- s+sum(y > x[j])
  }
  M <-c(M,s)
}
hist(M)

W <- M+10
hist(W)

\(U=W_1-\dfrac{n_1( n_{1}+1) }{2}\) から、 　　　\(W=M+\dfrac{4( 4+1) }{2}=U+10\) 　であることと一致している。