これまでは、独立２群に差がないか検定する方法として「２標本t検定」を学びました。
この検定法は母集団が正規分布に従っていることを前提としています。
母集団の正規性を前提とする検定法をパラメトリック検定といい、正規性を前提としない検定をノンパラメトリック検定といいます。

1 Wilcoxon順位和検定.

ここで扱うWilcoxonの順位和検定（rank sum test）は、データを「２群をあわせた順位」に変換し、順位をデータとして扱います。
2つのクラスの成績を比較する場合、素点を用いると外れ値の影響を受けやすいですが、データを順位にしてしまえば、近接する順位データの差は１であり、外れ値の影響を抑えることができます。

　帰無仮説は、「２群の順位分布に差がない」　です。

これは２群の分位点が等しいことを意味し、50% 分位点である中央値も等しくなるので、「２群の中央値の差の検定」と言われます。
順位分布に差がないという帰無仮説のもとでは、順位は２群の標本数に従って平等に配分され、標本数が同じならば順位和は同じになるはずです。（標本数が異なっていても順位平均は同じになります）

確率変数としての順位には注意すべき点があります。順位に変換されたデータは独立ではないからです。自分の順位が上がれば、他の順位が下がります。A群の順位和が決まると、B群の順位和が決まります。２群の順位和は独立ではないのです。（A群の順位平均が決まると、B群の順位平均も決まってしまうことになります）この順位データの独立性の問題は、分散を求める際にも注意する必要があります。順位平均差の分散も、独立性がないので簡単には求まりません。

ここでは、
帰無仮説：「２群の順位分布に差がない」のもとで、

各群の順位和がどのような分布になるのか.
対立仮説のもとで、各群の順位和がどのような分布になるのか.

シミュレーションで調べます。

まず、順位という確率変数の性質について見ていきましょう。
高校数学の数列の知識が役に立ちます。

2 順位分布

2.1 順位と一様分布.

大きさ10の乱数を発生させ、rank( )関数で順位をつけます。

  x <- rnorm(10,5,1)
  x

##  [1] 4.653982 5.448210 5.676536 4.397941 5.742300 4.651488 6.124564 4.530398
##  [9] 5.362304 6.697622

  r <- rank(x, ties.method = "min")
  r

##  [1]  4  6  7  1  8  3  9  2  5 10

1番目のデータの順位 r1は確率変数と考えられます。
rank( )関数が昇順の順位を返すことが分かります。
以下、順位は断らない限り昇順として話を進めます。

r1 を1000回出力したデータを S1 として分布を調べます。

S1=c()
n <- 10
for (i in 1:1000){
  x <- rnorm(n,5,1)
  r <- rank(x, ties.method = "min")
  r1 <- r[1]
  S1<- c(S1,r1)
}
hist(S1,prob=T,breaks=seq(0.5, 10.5, by = 1),main="",xlab="",col=0)

一様分布になっており、どの順位も同様に確からしいということになります。

2.2 平均と分散

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}\)

シミュレーションデータで確かめてみましょう。

m <- c(mean(S1),(n+1)/2)
m

## [1] 5.466 5.500

v <- c(var(S1),(n+1)*(n-1)/12)
v

## [1] 8.347191 8.250000

確かに一致しているようです。

2.3 順位和統計量

２群 \(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}\)

＜分散の導出＞.
　\(X_i\) が n 個の中での順位を表す確率変数として概略のみ記します。

\(\displaystyle V\left(\sum_i X_i\right)=\sum_i V(X_i)+2\sum_{i< j}Cov(X_i,X_j)\)
\(n=n_1+n_2\) とすると、\(Cov(X_i,X_j)=-\dfrac{n+1}{12}\)

これらを用いて順位和の分散を計算すると、

　　\(\begin{split} V(W_1)&=V\left(\sum_{i=1}^{n_1} R_{1i}\right)=\sum_{i=1}^{n_1} V(R_{1i})+2\sum_{i< j}^{n_1}Cov(R_{1i},R_{1j})\\ &=n_1\cdot \dfrac{(n+1)(n-1)}{12}+2\left(-\dfrac{n+1}{12}\right)\cdot \dfrac{n_1(n_1-1)}{2}\\ &=\dfrac{n_1n_2(n+1)}{12}=\dfrac{n_1n_2(n_1+n_2+1)}{12} \end{split}\).

＜注意＞.

導出の詳細についてはサイト「ひょっこり統計」の Wilcoxon を参照してください。
計算が煩雑なので、後で学習すると良いでしょう。
タイ（同順位）が存在する場合、分散を補正する必要がありますが、ここでは立ち入りません。
しかし、実際に検定をする場合は補正を行ってください。

3 小標本の場合

\(n_1=4,n_2=6\) の場合について、２群の順位分布を比較します。

3.1 同一分布の２群.

群G1（平均3, 標準偏差1）からの大きさ4の標本 \(x_1\)
群G2（平均3, 標準偏差1）からの大きさ6の標本 \(x_2\)

２つの標本を合わせて順位をつけ、順位データを r1,r2 とします。
同じ母集団からの標本なので、２群の順位分布に差はないはずです。

n1 <- 4
n2 <- 6
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)

ties.method = “average” は、同じ順位のとき平均順位にするというオプションです。
例えば、５位と６位が同じ点数なら 5.5位とします。

print(r1)

## [1] 7 3 2 6

print(r2)

## [1]  5  4  9  8  1 10

２群の順位分布を箱ひげ図で表示し比較します。

library(beeswarm)
list=list(A=r1, B=r2)
boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

同じ順位分布のはずですが、小標本だと変動が見られます。
順位和はどうでしょうか。

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

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

順位和は独立でなく、\(W_1+W_2=55\) という関係で結ばれているので、分散は等しくなります。

繰り返しプログラムで、順位和 \(W_1,~W_2\) を1000回出力します。
２群の順位和の分布を確かめましょう。

S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:1000){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
list=list(A=S1, B=S2)
par(family = "HiraKakuProN-W3") 
boxplot(list,main="順位和")

２群の標本数が異なるので、順位和の期待値は異なりますが、分散は同じになっています。

次は、順位平均の分布を比較します。

list=list(A=S1/n1, B=S2/n2)
par(family = "HiraKakuProN-W3") 
boxplot(list,main="順位平均")

順位平均の中央値は 5.5 で一致しています。
また、範囲や四分位範囲は相似形となっていることが観察されます。

　　　\(V(\bar{r_1})=\dfrac{V(r_1)}{4^2}=\dfrac{22}{16}\)　　　\(V(\bar{r_2})=\dfrac{V(r_1)}{6^2}=\dfrac{22}{36}\)

分散は \(\dfrac{36}{16}\)倍、標準偏差は \(\dfrac{6}{4}=1.5\)倍となり、範囲も 1.5 倍です。

max(S1/n1)-min(S1/n1)

## [1] 6

max(S2/n2)-min(S2/n2)

## [1] 4

3.2 ２群の順位和分布.

順位和 \(W_1,~W_2\) の分布の両方を求める必要はありません。

\(W_2=55-W_1\)より、
　　　\(E(W_2)=55-E(W_1)=55-22=33\)
　　　\(V(W_2)=V(55-W_1)=V(W_1)=22\)

順位平均の差も.

　　　\(\dfrac{W_1}{n_1}-\dfrac{W_2}{n_2}=\dfrac{W_1}{n_1}-\dfrac{55-W_1}{10-n_1}=\dfrac{10W_1−55 n_1}{n_1(10-n_1)}\)

\(W_1\)が分かれば、全てが決まってしまうからです。
このため、片方の順位和の分布を調べることが重要になります。

２群の順位和のヒストグラムを並置して比較します。

par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(S1,breaks = "Scott",main="W1",xlab="",col=0)
hist(S2,breaks = "Scott",main="W2",xlab="",col=0)

ヒストグラムはシンメトリー（対称）であることが分かります。
\(W_1\) のヒストグラムを裏返せば、\(W_2\) のヒストグラムに重ねることができます。

＜理由＞.

\(W_2=55-W_1\) なので、

\(W_1\)に \(W_1=x\)という値が出力されれば、
\(W_2\)に \(W_2=55-x\) という値が出力されます.

よって、\(W_1\)の \(x\)という値の度数と、\(W_2\)の \(55-x\) という値の度数は同じになります。
\(x,~55-x\)の中点は \(\dfrac{55}{2}\) なので、\(W_1,~W_2\) のヒストグラムは \(\dfrac{55}{2}=27.5\) について対称です。

hist(S1,main="",xlab="",col=0,xlim=c(8,50),breaks=seq(8.5, 47.5, by = 2))
hist(S2,main="",xlab="",col="lightblue",breaks=seq(8.5, 47.5, by = 2), add=T)
abline(v = 55/2,col="red")

赤の縦線は \(x=27.5\)であり、２つのヒストグラムは赤の縦線について左右対称になっています。
これで、統計量としては \(W_1,~W_2\) の片方を調べれば十分であることが分かります。
通常は標本数が小さいほうの順位和を統計量とします。

ここでは、\(W_1\) を統計量として使用します。

3.3 理論分布.

小標本なので、\(W_1\)の理論分布は組み合わせによる計算が可能です。

Rの関数 combn(x=c(1:10),m=4) は、1から10までの順位から4個の順位を取り出す組合せのすべてを出力します。
はじめの10個を表示します。

combn(x=c(1:10),m=4)[,1:10]

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    1    1    1    1    1    1    1    1    1     1
## [2,]    2    2    2    2    2    2    2    2    2     2
## [3,]    3    3    3    3    3    3    3    4    4     4
## [4,]    4    5    6    7    8    9   10    5    6     7

列の和が順位和ということになります。

列数は \({}_{10} C_4\) で、Rの関数choose(10,4) で計算されます。

N <- choose(10,4)
N

## [1] 210

次のコードで、210個の順位和が出力されます。

data <- combn(x=c(1:10),m=4)
N <- choose(10,4)
S=c()
for (i in 1:N){
  s <- sum(data[,i])
  S <- c(S,s)
}
S

##   [1] 10 11 12 13 14 15 16 12 13 14 15 16 17 14 15 16 17 18 16 17 18 19 18 19 20
##  [26] 20 21 22 13 14 15 16 17 18 15 16 17 18 19 17 18 19 20 19 20 21 21 22 23 16
##  [51] 17 18 19 20 18 19 20 21 20 21 22 22 23 24 19 20 21 22 21 22 23 23 24 25 22
##  [76] 23 24 24 25 26 25 26 27 28 14 15 16 17 18 19 16 17 18 19 20 18 19 20 21 20
## [101] 21 22 22 23 24 17 18 19 20 21 19 20 21 22 21 22 23 23 24 25 20 21 22 23 22
## [126] 23 24 24 25 26 23 24 25 25 26 27 26 27 28 29 18 19 20 21 22 20 21 22 23 22
## [151] 23 24 24 25 26 21 22 23 24 23 24 25 25 26 27 24 25 26 26 27 28 27 28 29 30
## [176] 22 23 24 25 24 25 26 26 27 28 25 26 27 27 28 29 28 29 30 31 26 27 28 28 29
## [201] 30 29 30 31 32 30 31 32 33 34

210個の順位和を table( ) 関数で整理し、相対度数を求めれば確率分布が分かります。

table(S)

## S
## 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 
##  1  1  2  3  5  6  9 10 13 14 16 16 18 16 16 14 13 10  9  6  5  3  2  1  1

sum(table(S)/N)

## [1] 1

table(S)/N

## S
##          10          11          12          13          14          15 
## 0.004761905 0.004761905 0.009523810 0.014285714 0.023809524 0.028571429 
##          16          17          18          19          20          21 
## 0.042857143 0.047619048 0.061904762 0.066666667 0.076190476 0.076190476 
##          22          23          24          25          26          27 
## 0.085714286 0.076190476 0.076190476 0.066666667 0.061904762 0.047619048 
##          28          29          30          31          32          33 
## 0.042857143 0.028571429 0.023809524 0.014285714 0.009523810 0.004761905 
##          34 
## 0.004761905

plot(table(S)/N,type="l",col="lightblue")

＜累積分布＞

分布関数は、cumsum( ) 関数で相対度数の累積和を求めればよく、

cumsum(table(S)/N)

##          10          11          12          13          14          15 
## 0.004761905 0.009523810 0.019047619 0.033333333 0.057142857 0.085714286 
##          16          17          18          19          20          21 
## 0.128571429 0.176190476 0.238095238 0.304761905 0.380952381 0.457142857 
##          22          23          24          25          26          27 
## 0.542857143 0.619047619 0.695238095 0.761904762 0.823809524 0.871428571 
##          28          29          30          31          32          33 
## 0.914285714 0.942857143 0.966666667 0.980952381 0.990476190 0.995238095 
##          34 
## 1.000000000

barplot( ) 関数を用いて、棒グラフで図示できます。
赤い横線は両側5% の棄却境界です。

par(family = "HiraKakuProN-W3") 
barplot(cumsum(table(S)/N),main="P(W1≦x)")
abline(h = c(0.025,0.975),col="red")

barplot(1-cumsum(table(S)/N),main="1-P(W1≦x)")
abline(h = c(0.025,0.975),col="red")

＜両側5% の棄却境界＞

　\(P(W_1\geq k)=1-P(W_1\leq k-1)\)より、
　　　\(P(W_1\geq 31)=1-P(W_1\leq 30)=1-0.9666667=0.33>0.25\).
　　　\(P(W_1\geq 32)=1-P(W_1\leq 31)=1-0.9809524=0.019<0.25\). 　　

　　　よって、右側境界 x=32　、　左側境界 x=12.

＜片側5% の棄却境界＞

　　　右側検定で x=31　、　左側検定で x=13.

となります。

3.4 帰無仮説.

２群に差がないの帰無仮説のもとで、\(W_1\)について

シミュレーションによる分布.
組合せ計算による理論分布.
正規近似による分布.

これらを比較します。

＜帰無仮説＞.

群G1（平均3, 標準偏差1）の正規分布からの大きさ4の標本 \(x_1\)
群G2（平均3, 標準偏差1）の正規分布からの大きさ6の標本 \(x_2\)

S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:1000){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,main="",ylab="",col=0,xlab="")
par(new=T)
plot(table(S)/N,type="l",col="blue", yaxt="n",ylab="", xaxt="n",xlab="")
curve(dnorm(x,22,sqrt(22)),add=T,col="green")#自由度

シミュレーションによる確率分布と理論分布（青）はほぼ一致しています。

また、同じ平均と標準偏差の正規曲線（緑）ともほぼ一致しますが、小標本の場合は理論分布によるEXACT検定が良いとされています。

3.5 対立仮説.

＜平均差がある場合＞.

２群の平均に差があるとき、２群の順位分布にも差が出てきます。

＜対立仮説＞.

群G1（平均3, 標準偏差1）の正規分布からの大きさ4の標本 \(x_1\)
群G2（平均4, 標準偏差1）の正規分布からの大きさ6の標本 \(x_2\)

n1 <- 4
n2 <- 6
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,4,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  list=list(A=r1, B=r2)
boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

「２群の平均に差がある」（平均差1）という対立仮説のもとで \(W_1\) を1000回出力し、有意差を検出できるかどうかシミュレーションします。

＜対立仮説＞.

群G1（平均3, 標準偏差1）の正規分布からの大きさ4の標本 \(x_1\)
群G2（平均2, 標準偏差1）の正規分布からの大きさ6の標本 \(x_2\)

S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:1000){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,2,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,main="",ylab="",col=0,xlab="")
par(new=T)
plot(table(S)/N,type="l",col="blue", yaxt="n",ylab="",cex.axis=0.5,xlab="")
abline(v = 31,col="blue")
curve(dnorm(x,22,sqrt(22)),add=T,col="green")
abline(v = qnorm(0.95,22,sqrt(22)),col="red")

赤の縦線は正規近似による右側棄却境界、青の縦線は理論分布による右側棄却境界です。赤と青の棄却境界が離れていますが、整数値補正を行わなかったためとも考えられます。
正確な分析のためには、理論分布によるEZACT検定が望ましいと考えられます。

平均差１での検出力は低いようです。

一般に、正規性を満たしていれば、t検定のほうが検出力が高いと言われています。
シミュレーションで検出力を比較することもできますが、ここでは立ち入りません。
探求課題として挑戦してください。

3.6 外れ値がある場合.

\(n_1=7,n_2=7\) の場合について、順位和の理論分布を求めます。

r1 の順位和\(W_1\) は確率変数で、期待値と分散は
　　\(E(W_1)=\dfrac{1}{2}\cdot 7\cdot15=52.5\) 　　\(V(W_1)=\dfrac{1}{12}\cdot7\cdot7\cdot15=61.25\)

data <- combn(x=c(1:14),m=7)
N <- choose(14,7)
S=c()
for (i in 1:N){
  s <- sum(data[,i])
  S <- c(S,s)
}
plot(table(S)/N,type="l",col="lightblue")
curve(dnorm(x,52.5,sqrt(61.25)),add=T,col="red")

標本数が増えると曲線は滑らかになり、正規分布との乖離も小さくなっています。

3.6.1 帰無仮説.

＜はずれ値を除くと平均差がない＞.

クラス成績を比較する際、ある優秀な生徒が常に1番であることがあります。
この1番の生徒の得点が外れ値に該当します。
いつも満点を取るとクラスの平均点が跳ね上がってしまいますが、統計量にはどのような形で現れるでしょうか。

＜帰無仮説＞.

平均5, 標準偏差2 の正規乱数6個に20を付加した \(x_1\)
平均5, 標準偏差2 の正規乱数7個 \(x_2\)

ともに大きさ7の標本ですが、\(x_1\) に外れ値があります。

n1 <- 7
n2 <- 7
x1 <- c(rnorm(n1-1,5,2),20);x2 <- rnorm(n2,5,2)
x1;x2

## [1]  2.875394  5.597399  4.673112  6.561924  4.368850  3.893005 20.000000

## [1] 3.287872 7.927668 3.180717 2.841752 3.724353 2.706374 5.014096

par(family = "HiraKakuProN-W3") 
list=list(x1=x1, x2=x2)
boxplot(list,main="生データの分布")
beeswarm(list,  pch = 16, add = TRUE)

x <- c(x1,x2)
r <- rank(x, ties.method = "average")
r1 <- head(r,n1);r1

## [1]  3 11  9 12  8  7 14

r2 <- tail(r,n2);r2

## [1]  5 13  4  2  6  1 10

list=list(r1=r1, r2=r2)
boxplot(list,main="順位データの分布")
beeswarm(list,  pch = 16, add = TRUE)

順位分布の差はあまりないようです。

順位和（昇順）の分布をみてみましょう。
\(x_1\) に外れ値がある状況では、昇順の順位和は \(W_1\) が少し大きくなるはずです。
外れ値の影響がどれくらいあるか見てみましょう。
同じ条件で1000回順位和を出力し、 \(W_1\) のヒストグラムと理論分布を比べます。

S1=c()
S2=c()
Y=c()
n1 <- 7
n2 <- 7
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,5,2),20);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,col=0, xaxt="n", yaxt="n",xlab="",main="順位和の分布")
par(new=T)
plot(table(S)/N,type="l",col="lightblue",ylab="",xlab="")

外れ値の影響はなく、順位和の分布は理論分布にほぼ当てはまっています。

3.6.2 対立仮説.

＜外れ値を除くと平均差がある＞.

平均3, 標準偏差2 の正規乱数6個に20を付加した \(x_1\)
平均5, 標準偏差2 の正規乱数7個 \(x_2\)

n1 <- 7
n2 <- 7
x1 <- c(rnorm(n1-1,3,2),20);x2 <- rnorm(n2,5,2)
x1;x2

## [1]  4.137425  1.578448  3.498493  1.243822  1.489587  3.452740 20.000000

## [1] 4.6114506 1.8000766 5.0139273 4.0645075 4.2183869 0.6588171 5.7338295

par(family = "HiraKakuProN-W3") 
list=list(x1=x1, x2=x2)
boxplot(list,main="生データの分布")
beeswarm(list,  pch = 16, add = TRUE)

x <- c(x1,x2)
r <- rank(x, ties.method = "average")
r1 <- head(r,n1);r1

## [1]  9  4  7  2  3  6 14

r2 <- tail(r,n2);r2

## [1] 11  5 12  8 10  1 13

list=list(r1=r1, r2=r2)
boxplot(list,main="順位データの分布")
beeswarm(list,  pch = 16, add = TRUE)

あまり差がないように見えるので、検出力は高くないかもしれません。
同じ条件で1000回順位和を出力し、ヒストグラムと理論曲線を重ねます。

S1=c()
S2=c()
Y=c()
n1 <- 7
n2 <- 7
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,3,2),20);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  y <- mean(x1)-mean(x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
  Y <- c(Y,y)
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,col=0, xaxt="n", yaxt="n",xlab="",main="順位和の分布")
par(new=T)
plot(table(S)/N,type="l",col="lightblue",ylab="",xlab="")
abline(v = qnorm(0.05,52.5,sqrt(61.25)),col="red")

赤の縦線は左側5%の棄却境界であり、その左側の値の検定会ようが棄却されます。
ある程度は平均差が検出されますが、検出力は高くありません。.

3.6.3 t検定との比較

順位和検定とヴェルチのt検定について、どちらが適切かという問題があります。
有意差が出やすいか、検出力はどうかなどです。
ここではP値の出力分布を調べます。

＜外れ値を含めて平均差がない場合＞.

G1に1番の生徒がいるとして、
外れ値を含めてクラス平均が同じであるときの２群の学力差を考えます。
外れ値を除くと平均差があることになるのですが、学力差は検出されるでしょうか。

平均4, 標準偏差2 の正規乱数6個に11を付加した \(x_1\)
平均5, 標準偏差2 の正規乱数7個 \(x_2\)

n1 <- 20
n2 <- 20
  x1 <- c(rnorm(n1-1,4,2),24);x2 <- rnorm(n2,5,2)
  mean <- c(mean(x1),mean(x2))
  mean

## [1] 4.654973 5.783670

標本数が小さいので、出力した平均には差がありますが、近い値になっています。

簡単のため Rの wilcox.exact()、t.test() を用いてP値を算出します。

library(exactRankTests)
wilcox.exact(x=x1,y=x2,paired=F)$p.value

## [1] 0.01547847

t.test(x1,x2)$ p.value

## [1] 0.369692

＜1000回繰り返し＞

S1=c()
S2=c()
PV=c()
PW=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,4,2),24);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
  PW=c(PW,wilcox.exact(x=x1,y=x2,paired=F)$p.value)
  PV=c(PV,t.test(x1,x2)$ p.value)
}
par(mfrow=c(1,2))
hist(PW,prob=T)
hist(PV,prob=T)

有意水準 0.05 付近では順位和検定の方が有利なようですが、さまざまなケースをシミュレーションしないと断言はできません。

S1=c()
S2=c()
PV=c()
PW=c()
n1 <- 10
n2 <- 10
for (i in 1:1000){
  x1 <- c(rnorm(n1,4,2));x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
  PW=c(PW,wilcox.exact(x=x1,y=x2,paired=F)$p.value)
  PV=c(PV,t.test(x1,x2)$ p.value)
}
par(mfrow=c(1,2))
hist(PW,prob=T)
hist(PV,prob=T)

外れ値のない場合は、t検定の方が有利なことが分かります。

4 大標本の場合.

これまで、小標本の順位和検定を考えてきました。
小標本であれば、組合せ計算で順位和の理論分布を計算することが可能です。
大標本になると、組合せ計算の負荷が大きく、筆者のノートパソコンではエラーを返します。
\(n_1=7,n_2=7\) の場合に、すでに正規分布に近かったので、大標本では順位和を正規近似できると考えられます。

２群 \(X_1=\left(x_{11},\cdots,x_{1n_1}\right)\:\:,\:\: X_2=\left(x_{21},\cdots,x_{2n_2}\right)\) について、
\(X_1\) の順位和\(W_1\)の期待値と分散は、
　　　\(\displaystyle E(W_1)=\dfrac{n_1(n_1+n_2+1)}{2}\)　　\(V(W_1)=\dfrac{n_1n_2(n_1+n_2+1)}{12}\)

\(z = \dfrac{W_1-E(W_1)}{\sqrt{Var(W_1)}}\)は標準正規分布に従うはずです。

4.1 \(n_1=20,n_2=20\)

S1=c()
S2=c()
Y=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- rnorm(n1,4,1);x2 <- rnorm(n2,4,1)
  #x1 <- c(rnorm(n1-1,5,2),20);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  #y <- mean(x1)-mean(x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
  #Y <- c(Y,y)
}
Z1  <- (S1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
par(family = "HiraKakuProN-W3") 
hist(Z1,prob=T,main="",ylab="",col=0,xlab="")
curve(dnorm(x,0,1),add=T,col="blue")
abline(v = qnorm(0.95),col="red")

正規分布で近似できると考えてよさそうです。

対立仮説のもとで、平均差を検出できるかシミュレーションします。

対立仮説.

S1=c()
S2=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- rnorm(n1,4,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
Z1  <- (S1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
par(family = "HiraKakuProN-W3") 
hist(Z1,prob=T,main="",ylab="",col=0,xlab="")
curve(dnorm(x,0,1),add=T,col="blue")
abline(v = qnorm(0.95),col="red")

平均差１の順位和検定では、順位和の正規化によって検出力が80%を超えている。

標本数によって検出力が変わる.
標本数が大きい方が有利.
分散が大きくなれば差の検出能力は上がる.
外れ値は不利にはならない.
どれくらいの標本数で、どれくらいの検出力があるか、十分シミュレーションすることが不可欠.

標本設計の重要性

4.2 \(n_1=10,n_2=10\)

S1=c()
S2=c()
Y=c()
n1 <- 10
n2 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,4,1);x2 <- rnorm(n2,4,1)
  #x1 <- c(rnorm(n1-1,5,2),20);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  #y <- mean(x1)-mean(x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
  #Y <- c(Y,y)
}
Z1  <- (S1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
par(family = "HiraKakuProN-W3") 
hist(Z1,prob=T,main="",ylab="",col=0,xlab="")
curve(dnorm(x,0,1),add=T,col="blue")
abline(v = qnorm(0.95),col="red")

S1=c()
S2=c()
n1 <- 10
n2 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,4,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
Z1  <- (S1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
par(family = "HiraKakuProN-W3") 
hist(Z1,prob=T,main="",ylab="",col=0,xlab="")
curve(dnorm(x,0,1),add=T,col="blue")
abline(v = qnorm(0.95),col="red")

検出力は、標本数が２０ずつの場合に比べて小さいようです。
一般に、次のようなことが言えます。

標本数は、検出力に影響する。
https://www.saecanet.com/2025/06/rでa-bテスト平均値の比較/
正規性と等分散があればt検定のほうが頑強である。
順位和検定では標本数を大きくしないと検出力は上がらない

要は、大標本であれば各群の標本平均は正規分布に近似できるので、パラメトリック検定が可能なはずです。
大標本の場合、２標本t検定と順位和検定のどちらが有利かを考えてみましょう。

4.3 P値の比較.

外れ値のある場合は統計検定を理解せずに使っている人のために II.
は多少検出力が落ちるらしいが，経験的に，たとえば， n = 6 程度で t 検定と順位和検定を平行して行ってみると，ほとんどの場合，検定結果は同じになり，分布の偏った場合にのみ異なる結果となる.したがって，研究者はパラメトリックとノンパラメトリック検定を上手く利用することを考えたほうが良い. また，外れ値があるとき，棄却検定で外すことが可能になる場合もあるが，安易に外すことには問題があるので，外さずにノンパラメトリック検定を行ってみる価値はある.ノンパラメトリック検定でも有意差が得られなければ，あきらめることができる.外れ値の取り扱いについては次回に述べるが，思い悩むデータがあるときは，まずは，パラメトリック検定とノンパラメトリック検定，さらには棄却検定の結果を見てから考えてみてはどうであろうか?

対応のないt検定における正規性の仮定とサンプルサイズ

4.4 P値の分布の比較.

＜1000回繰り返し＞

S1=c()
S2=c()
PV=c()
PW=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,4,2),11);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
  PW=c(PW,wilcox.test(x=x1,y=x2,paired=F)$p.value)
  PV=c(PV,t.test(x1,x2)$ p.value)
}
par(mfrow=c(1,2))
hist(PW,prob=T,ylim=c(0,3))
hist(PV,prob=T,ylim=c(0,3))

外れ値を含む場合、順位和検定の方が有効性が高いようである。

S1=c()
S2=c()
PV=c()
PW=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- rnorm(n1,4,2);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
  PW=c(PW,wilcox.test(x=x1,y=x2,paired=F)$p.value)
  PV=c(PV,t.test(x1,x2)$ p.value)
}
par(mfrow=c(1,2))
hist(PW,prob=T,ylim=c(0,5))
hist(PV,prob=T,ylim=c(0,5))

外れ値がなく正規性と等分散性が満たされる場合は、t検定の方がよい。

まとめ.

池田先生の言うように、

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

これまで順位和検定について説明してきました。

２群 \(X_1=\left(x_{11},\cdots,x_{1n_1}\right)\:\:,\:\: X_2=\left(x_{21},\cdots,x_{2n_2}\right)\) について、\(N=(n_1+n_2)\) 個の中での順位を考え、

\(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)\) とするとき、

\(X_1\:,\:X_2\) の順位和　\(\displaystyle W_1=\sum_{i=1}^{n_1} R_{1i}\:\:,\:\:W_2=\sum_{i=1}^{n_2} R_{2i}\) 　が Wilcoxonの順位和統計量です。

\(W_1+W_2=\dfrac{(n_1+n_2)(n_1+n_2+1)}{2}=\dfrac{N(N+1)}{2}\)　であるので、\(W_1,W_2\) は独立ではありません。

\(W_1,W_2\) のどちらを統計量として用いても検定の結果は同じですが、計算上、標本数の小さい方を統計量として用います。.

\(W_1,W_2\) は期待値は異なりますが分散は共通です。

\(\displaystyle E(W_1)=\dfrac{n_1(n_1+n_2+1)}{2}=\dfrac{n_1n_2}{2}+\dfrac{n_1(n_1+1)}{2}\).
\(E(W_2)=\dfrac{n_2(n_1+n_2+1)}{2}=\dfrac{n_1n_2}{2}+\dfrac{n_2(n_2+1)}{2}\).
\(V(W_1)=V(W_2)=\dfrac{n_1n_2(n_1+n_2+1)}{12}\).
２つのヒストグラムは、\(x=\dfrac{N(N+1)}{4}\) について対称です。

5.1 U統計量

U統計量は、どのようにして考え出されたのでしょうか。

奥村晴彦先生「Wilcoxon-Mann-Whitney検定」によれば、Wilcoxon が1945年に順位和統計量を提案し，1947年に Mannと Whitney が U統計量を導入しました。U統計量は順位和統計量の進化型ということになります。詳しく見ていきましょう。

「各群の中での順位と、２群を合わせた順位の関係」を考えます。
２群を合わせた順位のほうが大きくなりますが、順位和はどうでしょうか。

　\(X_1\) の中での順位を \(r_{X_1}=\left(r_{11},\cdots,r_{1n_1}\right)\)とすると、
　　　　その順位和は　\(\displaystyle \sum_{i=1}^{n_1} r_{1i}=\dfrac{n_1(n_1+1)}{2}\).

　\(X_2\) の中での順位を \(r_{X_2}=\left(r_{21},\cdots,r_{2n_1}\right)\)とすると.
　　　　その順位和は　\(\displaystyle \sum_{i=1}^{n_1} r_{2i}=\dfrac{n_2(n_2+1)}{2}\).

前節の \(E(W_1),E(W_2)\) の式の第２項は、各群の標本数に依存する値であり、本質的な部分は第１項の \(\dfrac{n_1n_2}{2}\) と考えられます。実際、

\(\displaystyle E\left[W_1-\dfrac{n_1(n_1+1)}{2}\right]=\dfrac{n_1n_2}{2}\)
\(E\left[W_2-\dfrac{n_2(n_2+1)}{2}\right]=\dfrac{n_1n_2}{2}\)

これら \(\dfrac{n_1n_2}{2}\) を期待値とする２つの統計量を U統計量といい、\(U_1,U_2\) と書きます。

＜U 統計量の定義＞

\(U_1=W_1-\dfrac{n_1(n_1+1)}{2}\).
\(U_2=W_2-\dfrac{n_2(n_2+1)}{2}\).

これらは、（２群を合わせた順位和）ー（各群の順位和） となっています。

＜注意＞.

　発見の順番としては、Wilcoxonの順位和統計量が最初です。（1945年）.
Wilcoxonの統計量から本質的な部分を抽出し改良したものが U統計量です。（1947年）
一見すると \(U_1,U_2\) は \(W_1,W_2\)で表されており、数理的には意味のないものに見えます。
しかし、\(U_1,U_2\)は同一の確率分布に従います。また、順位差を用いて統計量を計算することが可能になります。
これについては次節で説明します。

U統計量は、順位差の和で表すことができます。

\(U_1=\displaystyle\sum_{i=1}^{n_1}R_{1i}-\sum_{i=1}^{n_1} r_{1i}=\sum_{i=1}^{n_1}(R_{1i}-r_{1i})\).
\(U_2=\displaystyle\sum_{i=1}^{n_2}R_{2i}-\sum_{i=1}^{n_1} r_{2i}=\sum_{i=1}^{n_2}(R_{2i}-r_{2i})\).

「２群を合わせた順位」と「各群の中での順位」の差の和となっています。
次に示すように、この差は順位に頼ることなく計算できます。

＜順位差による計算法＞

そもそも順位とは何でしょうか。
テストの順位が気になるとき、自分より点数の高い人が何人いるかが気になります。
順位は自分より上位の人数に 1 を加えたものです。
成績が10位であることは、自分より上に９人いることを意味します。

\(X_1\)の要素 \(x_{1i}\) について、
\(X_1\) の中での順位を \(r_{1i}\) 、２群をあわせた中での順位を \(R_{1i}\) とします。

\(X_2\) の中に \(x_{1i}\) より大きいものが \(u_{i}\) 個あるとすると、

\(X_1\) の中での順位は \(r_{1i}\)
A,Bあわせた中での順位は、\(R_{1i}=r_{1i}+u_{i}\)

成績で例えると、他群の自分より上位の生徒の数だけ順位が下がる訳です。

\(R_{1i}-r_{1i}=u_{i}\) であるので、
　　　\(U_1=\displaystyle\sum_{i=1}^{n_1}(R_{1i}-r_{1i})=\sum_{i=1}^{n_1}u_i\).

よって、\(U_1= \displaystyle \sum_i(x_{1i}\) より大きい他方の群のデータ数）

\(X_1\) の各要素について、自分より大きい他群の要素をカウントし合計するだけで U統計量が計算できることになり、手計算でも容易に計算可能です。

5.2 計算例.

＜Wilcox.test( )関数＞

Rには wilcox.test( ) 関数があり、２群を比較するノンパラメトリック統計量が出力されます。
どのような計算をしているのか見ておきましょう。

n1 <- 4
n2 <- 6
x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
x <- c(x1,x2)
r <- rank(x, ties.method = "average")
r1 <- head(r,n1)
r2 <- tail(r,n2)
wilcox.test(x1, x2)$ statistic

##  W 
## 13

wilcox.exact(x=x1,y=x2,paired=F)$ statistic

##  W 
## 13

＜順位和統計量＞

W1=sum(r1)
W2=sum(r2)
c(W1,W2)

## [1] 23 32

Wilcoxonの提案した順位和統計量は、２群のどちらも wilcox.test( )の統計量と一致しません。
wilcox.test( ) で出力されるのは、進化型のU統計量です。

＜U統計量＞

\(U_1=W_1-\dfrac{n_1(n_1+1)}{2}\).
\(U_2=W_2-\dfrac{n_2(n_2+1)}{2}\).

U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
c(U1,U2)

## [1] 13 11

wilcox.test( ) では、U1が出力されています。

＜降順の場合＞.

デフォルトの順位は昇順です。降順にするとどうなるか見てみましょう。
rank(-x) とすれば、xの降順の順位が出力されます。

r <- rank(-x, ties.method = "average")
r1 <- head(r,n1)
r2 <- tail(r,n2)
U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
c(U1,U2)

## [1] 11 13

降順では U1,U2の値が入れ替わるだけです。

5.3 U1,U2 の関係

\(W_1,W_2\) の間には、

\(W_1+W_2=\dfrac{(n_1+n_2)(n_1+n_2+1)}{2}=\dfrac{N(N+1)}{2}\)　という関係があります。

\(U_1,U_2\) の間の関係は、

　　\(U_1=W_1-\dfrac{n_1(n_1+1)}{2}\)　　　\(U_2=W_2-\dfrac{n_2(n_2+1)}{2}\)　より、.

\(U_1+U_2=\dfrac{N(N+1)}{2}-\dfrac{n_1(n_1+1)}{2}-\dfrac{n_2(n_2+1)}{2}\)

\(\dfrac{N(N+1)}{2}=\dfrac{n_1(n_1+n_2+1)}{2}+\dfrac{n_2(n_1+n_2+1)}{2}\)　を用いると、

\(\dfrac{n_1(n_1+n_2+1)}{2}-\dfrac{n_1(n_1+1)}{2}=\dfrac{n_1n_2}{2}\)
\(\dfrac{n_2(n_1+n_2+1)}{2}-\dfrac{n_2(n_2+1)}{2}=\dfrac{n_1n_2}{2}\)

よって、　\(U_1+U_2=n_1n_2\) が成立します。

　　　 \(U_1,U_2\) の確率分布は \(x=\dfrac{n_1n_2}{2}\) について対称であり、データも\(\dfrac{n_1n_2}{2}\)について対称に出力します。

組合せ計算を用いた理論分布の図において、\(U_1\)を青の縦線、\(U_2\) を赤の縦線で表示します。

data <- combn(x=c(1:10),m=4)
N <- choose(10,4)
S=c()
for (i in 1:N){
  s <- sum(data[,i])
  S <- c(S,s)
}
plot(table(S-4*5/2)/N,type="l",col="lightblue")
abline(v = 4*6/2,col="lightblue")
abline(v = U1,col="blue")
abline(v = U2,col="red")

\(U_1\)が右側有意なら、\(U_2\) は左側有意で、P値は等く、どちらも統計量として使用可能です。
しかし、検定で両方の値を用いる必要はありません。

通常は\(U_1,U_2\) の小さほうを統計量とすることになっており、数表も左側境界値のみが記されています。

5.4 もう一つの公式.

統計書によって、U統計量の定義が異なっていることがあり、戸惑うことがあります。

\(W_1+W_2=\dfrac{N(N+1)}{2}\) を用いると、
　　　\(U_1=\dfrac{N(N+1)}{2}-W_2-\dfrac{n_1(n_1+1)}{2}\)

さらに \(\dfrac{N(N+1)}{2}=\dfrac{n_1(n_1+n_2+1)}{2}+\dfrac{n_2(n_1+n_2+1)}{2}\)　より、

　　　\(\dfrac{N(N+1)}{2}-\dfrac{n_1(n_1+1)}{2}=\dfrac{n_1n_2}{2}+\dfrac{n_2(n_1+n_2+1)}{2}\).

よって、\(U_1=\dfrac{n_1n_2}{2}+\dfrac{n_2(n_1+n_2+1)}{2}-W_2\).

さらに、\(U_1=n_1n_2+\dfrac{n_2(n_2+1)}{2}-W_2\).

同様に、\(U_2=n_1n_2+\dfrac{n_1(n_1+1)}{2}-W_1\).

これらを統計量とすることもできます。

一般には、上式を用いて両方を計算し、小さいほうを統計量として数表を調べることになっています。これは、便宜的なもので、どちらの統計量も同等です。

上式は、

\(U_1=n_1n_2-U_2\).
\(U_2=n_1n_2-U_1\).

と変形され、 \(U_1+U_2=n_1n_2\) であることを意味します。

順位和統計量の進化型であるU統計量の長所は \(U_1,U_2\) の確率分布が同じになることですが、次節で説明するように順位差を利用した計算が可能になることです。

この方法は、手計算でデータを処理する場合に有効と考えられます。

シミュレーションで確かめましょう。

5.5 順位差による計算.

\(X_1\) の中での順位を \(r_{1i}\) 、２群をあわせた中での順位を \(R_{1i}\) とすると、

　　　\(U_1=\displaystyle\sum_{i=1}^{n_1}(R_{1i}-r_{1i})\)　　が成立します。

\(\sum\) の中の \((R_{1i}-r_{1i})\) については、

順位が降順のとき、\(x_{1i}\)より大きい他方の群のデータ数
順位が昇順のとき、\(x_{1i}\)より小さい他方の群のデータ数

を意味します。

このことを用いると、順位和計算は簡単になります。

5.5.1 順位が降順のとき

n1 <- 4
n2 <- 6
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  list=list(A=x1, B=x2)
boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

箱ひげ図のドットを見て、順位差を根拠とした

\(U_1= \displaystyle \sum_i(x_{1i}\) より大きい他方の群のデータ数）.

による計算をしてください。
ほとんど黙視だけで、順位和が分かります。

同じ計算をRで行うと、

  s = 0
  for (j in 1:n1){
    s <- s+sum(x2 > x1[j])
  }
s

## [1] 18

集合Xと定数aに対し、sum(X > a) は、aより大きい X の要素の数をカウントします。

次に、U統計量の定義により、

\(U_1=W_1-\dfrac{n_1(n_1+1)}{2}\) を計算します。.
順位はr(-x)により、降順に直します。

  x <- c(x1,x2)
  r <- rank(-x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  sum(r1)-n1*(n1+1)/2

## [1] 18

どちらの方法でも同じ値になることが確認されました。　

5.5.2 順位が昇順のとき.

昇順の順位でもU統計量は計算できます。

箱ひげ図のドットを見て、順位差を根拠とした

\(U_1= \displaystyle \sum_i(x_{1i}\) より小さい他方の群のデータ数）.

を計算をしてください。

boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

同じ計算をRで行うと、

  s = 0
  for (j in 1:n1){
    s <- s+sum(x2 < x1[j])
  }
s

## [1] 6

集合Xと定数aに対し、sum(X < a) は、aより小さい X の要素の数をカウントします。

次に、U統計量の定義により、

\(U_1=W_1-\dfrac{n_1(n_1+1)}{2}\) を計算します。.
順位はr(x)で昇順のままです。

r <- rank(x, ties.method = "average")#昇順
r1 <- head(r,n1)
r2 <- tail(r,n2)
U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
U1

## [1] 6

やはり、どちらの方法でも同じ値になることが確認されました。　

注意

昇順と降順では、\(U_1,U_2\)の値が入れ替割っていることに注意してください。

U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
c(U1,U2)

## [1]  6 18

どちらを統計量としても問題ありませんが、このような仕組みを手計算せずブラックボックスのままにしておいたのでは心配になります。

5.6 参考.

\(U_1,U_2\) の対称性は箱ひげ図からも分かります。

S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:10){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1-rank(x1)))
  S2<- c(S2,sum(r2-rank(x2)))
}
list=list(A=S1, B=S2)
par(family = "HiraKakuProN-W3") 
boxplot(list,main="順位差の和")
beeswarm(list,  pch = 16, add = TRUE)

\(U_1,U_2\) の分布は、12を中心に上下反転した形状です。.

6 各統計量のヒストグラム.

最後に、これまで出てきた統計量について、ヒストグラムで確認しておきましょう。

S1=c()
S2=c()
WT=c()
U1=c()
U2=c()
n1 <- 4
n2 <- 6
for (i in 1:1000){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1 <- c(S1,sum(r1))
  S2 <- c(S2,sum(r2))
  U1 <- c(U1,sum(r1)-n1*(n1+1)/2)
  U2 <- c(U2,sum(r2)-n2*(n2+1)/2)
  WT <- c(WT,wilcox.test(x1, x2)$statistic)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(S1,breaks = "Scott",main="W1",xlab="",col=0)
hist(S2,breaks = "Scott",main="W2",xlab="",col=0)

hist(U1,breaks = "Scott",main="U1",xlab="",col=0)
hist(U2,breaks = "Scott",main="U2",xlab="",col=0)

hist(WT,breaks = "Scott",main="WT",xlab="",col=0)

１つ１つの公式を自分で確かめる。
順位和検定は、高校生の探究課題として格好の教材と考えます。

6.1 文献.

中外製薬久力洋「多重比較を正しく使うために -Dunnett 型多重比較を事例として-3.正規分布などの特定の分布が仮定できない場合.
リンクは第2期医薬安全研定例会のアーカイブです。

奥村晴彦先生「Wilcoxon-Mann-Whitney検定」

九州大学附属図書館.

https://sigma-eye.com/2019/05/16/mann-whitney/. 資料

https://corvus-window.com/excel_mann-whitney-test/ 化学者のF・ウィルコクソンと，経済学者のH・B・マンと大学院生のホイットニーが同年代にそれぞれで考え出された手法になります．

https://en.wikipedia.org/wiki/Frank_Wilcoxon アイルランド出身のアメリカ人化学者(2 September 1892 – 18 November 1965).

https://en.wikipedia.org/wiki/Henry_Mann.
Henry Mann(27 October 1905, Vienna – 1 February 2000) オーストリア出身、アメリカ人.

https://en.wikipedia.org/wiki/D._Ransom_Whitney Donald Ransom Whitney (November 27, 1915 – August 16, 2007[

ウィルコクソンの順位和検定の棄却域の境界.

統計検定のまとめブログ.

https://www.jstage.jst.go.jp/article/jamt/advpub/0/advpub_24J1-6/_pdf/-char/ja

浜田知久馬

https://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q12203549992

@kenkenvw様　「Wilcoxonの順位和検定（ノンパラメトリック検定）の使用について」

外れ値がある場合に平均差で比較するのは適切ではありません。

12章 ノンパラメトリック検定