ウィルコクソンの順位和検定、マン・ホイットニーのU検定

ウィルコクソンの順位和検定（Wilcoxon rank-sum test）、またはマン・ホイットニーU検定（Mann-Whitney U test）は、2つの独立したサンプルが同じ母集団から得られたものかどうかを検定するためのノンパラメトリック検定です。データが正規分布に従わない場合や、サンプルサイズが小さい場合に有用です。両検定は実質的に同じ方法なので、まとめてマン・ホイットニー・ウィルコクソン検定(Wilcoxon-Mann-Whitney test)とも呼ばれています。

適用
パッケージ
ウィルコクソンの順位和検定
P値を直接求める方法
正規近似による方法
マン・ホイットニーのU検定

適用

ChatGPT様のお陰で複雑な図も簡単に描けるようなりました・・・

正規分布から大きく逸脱している場合にはｔ検定の結果が信頼できなくなる可能性があるので、ウィルコクソンの順位和検定のようなノンパラメトリック検定が適用されます

データの分布の形状に依存しないため、データの分布が偏っている場合やサンプルサイズが小さい場合に適用できます

location shift model（2つのグループのデータ分布の形は同じで位置が異なる場合）のもとでは、ｔ検定とおなじ平均値の差の検定が可能になります

パッケージ

この記事で使用するパッケージです（パッケージのインストール方法）

library(ggplot2)
library(dplyr)
library(coin)

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

サンプルの作成

ｐ値の算出方法を解説するためにサンプルサイズを小さくしています。解析対象４例という非現実的な例題ですが、ご理解ください。

data <- c(2, 5, 1, 3)
group <- c("A", "A", "B", "B")
dat <- data.frame(data, group)
#View(dat)

Aの平均値＝3.5、Bの平均値＝2

順位を追加

順位を表示 (群に関係なく小さい順に順位をつけます)

#dataの順位列 (rank) を追加
dat2 <- dat %>%
  mutate(rank = rank(data, ties.method = "average"))
#順位でソート
dat3 <- dat2 %>%
  arrange(rank)
# 順位列に「位」を追加
dat4 <- dat3 %>%
  mutate(rank = paste0(rank, "位"))
#View(dat4)

図でイメージしてみましょう

# 一直線上にデータポイントをプロットするグラフの描画
ggplot(dat4, aes(x = data, y = 1, color = group)) + 
  geom_point(size = 4, position = position_dodge(width = 0.2)) +  # データポイント
  geom_text(aes(label = rank), vjust = -1, position = position_dodge(width = 0.2)) +  # 順位のテキスト
  theme_minimal() + 
  labs(x = "data", y = "", title = "dataの順位") +
  scale_color_manual(values = c("A" = "blue", "B" = "red")) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

仮説（重要！）

目的：A群2名の得点が(2点, 5点)、B群２名の得点が(1点, 3点)　となった場合に、２群間に差があるかを検定したい（有意水準5%）

A群の母集団の得点分布とB群の分布が統計的に有意に異なるかを検定することになります

帰無仮説: A群とB群の得点分布に差がない

帰無仮説として以下の3つの仮説が考えられます

対立仮説１: A群の得点分布は、B群よりも高い（A群 > B群）
対立仮説２: A群の得点分布は、B群よりも低い（A群 < B群）
対立仮説３: A群とB群の得点分布には差がある（A群 $\neq$ B群）

この検定は、中央値ではなく得点分布の差を比較する手法なので、敢えて中央値という言葉を使用していません。また対立仮説３は、「A群がB群より高いか、もしくは低いか・・・」という仮説になっています。研究目的にマッチしているかどうか十分吟味して、どの仮説にするかを決定してください。

Rでウィルコクソンの順位和検定

基準（reference）の設定

Rは、どちらの群が高いか低いか・・・という答えは出しません

したがって、基準の設定が重要になります

基準を設定するためには、名義変数を factor変数 に変換します

必ず下記の操作を実行してください

dat$group <- factor(dat$group)

Rはアルファベット順に基準が設定されます

A群とB群がある場合、自動的にA群が基準になりま

基準をB群に変更したい場合は以下の記事をご参照ください

Rにデフォルトで含まれているwilcox.testを使用します。サンプルサイズが小さく、タイ（順位が同じデータ）がない場合に適用可能です。ただし、サンプルサイズが小さいので検出力は低下します。
後半の例題で使用するcoinパッケージのwilcox_testは、タイがある場合やサンプルサイズが大きい場合にも適用可能です。

対立仮説1の場合（alternative = “greater”）

基準がA群なので、対立仮説（alternative = “greater”）は、A群 > B群を意味します

したがってp<0.05となった場合、統計的に稀な現象と見なされ、帰無仮説（A群＝B群）が棄却され、A群＞B群が支持されます

wilcox_greater <- wilcox.test(
  data ~ group,
  data = dat,
  alternative = "greater"
)
print(wilcox_greater)

p=0.333ですので、統計的に有意な差（A群 > B群）があるとは言えず、帰無仮説（A群とB群には差がない）を棄却することができません

対立仮説2の場合（alternative = “less”）

基準がA群なので、対立仮説（alternative = “less”）は、A群 < B群を意味します

p<0.05になった場合、統計的に稀な現象と見なされ、帰無仮説（A群とB群には差がない）が棄却され、A群 < B群が支持されます

wilcox_less <- wilcox.test(
  data ~ group,
  data = dat,
  alternative = "less"
)
print(wilcox_less )

p=0.8333ですので、統計的に有意な差（A群 < B群）があるとは言えず、帰無仮説（A群＝B群）を棄却することができません

対立仮説3の場合（alternative = “two.sided”）

p<0.05になった場合、統計的に稀な現象と見なされ、帰無仮説（A群とB群には差がない）が棄却されます

wilcox_two.sided <- wilcox.test(
  data ~ group,
  data = dat,
  alternative = "two.sided"
)
print(wilcox_two.sided )

注意）wilcox.testではデフォルトになっているので、alternative = ” “ の記述がない場合は、この両側検定になるので注意してください

p=0.6667ですので、統計的に有意な差（A群 $\neq$ B群）があるとは言えず、帰無仮説（A群$=$B群）を棄却することができません

P値を直接求める方法

サンプルサイズが小さい場合、組合せからｐ値を算出できます

Rの wilcox.test もサンプルサイズが小さい場合には、この方法で算出しているようです（すみません、ChatGPTの情報ですので不確かです。。。）

wilcox.testで求められる p 値は、特にサンプルサイズが小さい場合には、全ての可能なデータの順位の組み合わせから直接計算されることがあります。これは「正確検定」（exact test）と呼ばれ、統計的検定の精度を高めるために行われます。
ChatGPT

（再掲）

A群(2位, 4位)、A群の順位和 $ W_A=6$
B群(1位, 3位)、B群の順位和 $ W_B=4$

対立仮説1の場合（A群＞B群）

片側検定

全ての可能な組み合わせの中で $W_A≧6$ となる確率

$\dfrac{2}{6}=0.3333$

この考え方は比率の検定をご確認ください

対立仮説2の場合（A群＜B群）

反対側の片側検定

全ての可能な組み合わせの中で $W_A≦6$ となる確率

$\dfrac{5}{6}=0.8333$

対立仮説3の場合（A群$\ne$B群）

両側検定のp値

A群の順位和≧6または反対側のA群の順位和≦4となる確率

$\dfrac{4}{6}=0.6666$　

両側検定には、以下のような検定統計量を利用した考え方もあるようですが、ご助言いただける方は、下のコメント欄にご記入いただければ助かります。

ウィルコクソンの順位和検定に使う検定統計量(どちらか一方の順位和)

$W_A=\sum_{i=1}^{n_A}(R_i)$

基準（A群）の順位和の期待値

$E(W_A)=\frac{n_A(n_A+n_B+1)}{2}$

$E(W_A)=\frac{2(2+2+1)}{2}=5$

A群の順位和（実測値）と期待値との差：$6-5=1$

差が $1$ になるような、反対側の順位和は、$4$ (5-4=1)

両側検定のｐ値＝ A群の順位和≧6または反対側のA群の順位和≦4となる確率

正規近似による方法

データの準備と要約

data <- c(
  12, 4, 9, 6, 2, 10, 12, 11, 7, 10,
  10, 10, 25, 13, 18, 7, 15, 9, 18, 19)
group <- c(rep("X",10), rep("Y",10))
dat_xy <- data.frame(data, group)
#View(dat_xy)

順位を追加

#dataの順位列 (rank) を追加
dat_xy2 <- dat_xy %>%
  mutate(rank = rank(data, ties.method = "average"))
#順位でソート
dat_xy3 <- dat_xy2 %>%
  arrange(rank)
#View(dat_xy3)

片側検定（帰無仮説Y群>X群）

ここからは、coinパッケージのwilcox_testを使用します

#groupをファクター変数に変更してYを基準に変更
dat_xy$group <- factor(dat_xy$group)

yを基準にして片側検定

dat_xy$group <- relevel(dat_xy$group, ref="Y") 
coin::wilcox_test(
  data ~ group,
  data = dat_xy,
  alternative = "greater"
)

有意水準5％のもとで、YはXより優位に高いことが分かりました

wilcox_testの Z値の求め方

Z=2.2798 の算出方法について解説します

X群とY群の順位和 $(W)$

rank_sum <- dat_xy3 %>%
  group_by(group) %>%
  summarise(sum_rank = sum(rank))
print(rank_sum)

X群の順位和: $W_x=75$
Y群の順位和: $W_y=135$
X群のサイズ: $n_x$
Y群のサイズ: $n_y$
サンプルサイズ: $N=n_x+n_y$

Z値の求め方

Y群を基準として記載

順位和の期待値 $E[W]=\dfrac{n_y(n+1)}{2}$

順位和の分散 $V[W]=\dfrac{n_xn_y(ｎ+1)}{12}$

Z値

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

# グループ別のランク和の計算
Wx = 75
Wy = 135

# サンプルサイズ
nx <- sum(group == "X")
ny <- sum(group == "Y")
n <- nx + ny

# U統計量の計算
# Y群を基準にしているので、ここではWyを使用します


#期待値
EW <- ny*(n+1)/2
#分散
VW <-  nx*ny*(n+1)/12
# Z値の計算
Zw <- (Wy - EW) / sqrt(VW)

# 結果の出力
Zw

#p値
pnorm(Zw, lower.tail = FALSE)

coinパッケージのwilcox_testの答えと微妙に違うようです・・・？
タイデータを考慮したマン・ホイットニーU検定を実行してみましょう

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

タイデータのがあるので補正したZ値を求めます

マン・ホイットニーのU統計量　$U=n_xn_y+\dfrac{n_y(n_y+1)}{2}-W_y$

順位和の期待値 $E[U]=\dfrac{n_xn_y}{2}$

補正項 $T=\sum{t_i^3-t_i}$ （t=タイデータの個数）

分散　$V[U]=\dfrac{n_xn_y(N+1)}{12}-\dfrac{n_xn_y}{12N(N-1)}T$

検定統計量

$z=\dfrac{|U-E[U]|}{\sqrt{V[U]}}$

# グループ別のランク和の計算
Wx = 75
Wy = 135

# サンプルサイズ
nx <- sum(group == "X")
ny <- sum(group == "Y")
n <- nx + ny

# U統計量の計算
# Y群を基準にしているので、ここではWyを使用します
U <- nx*ny+ny*(ny+1)/2-Wy

#期待値
EU <- nx*ny/2

# タイデータの個数
t　<- table(dat_xy3$rank) 
# タイの修正項
T <- sum(t^3 - t) 
#分散
VU <- (nx*ny/12)*((n+1)-T/(n*(n-1)))

# Z値の計算
Zu <- abs(U - EU) / sqrt(VU)

# 結果の出力
Zu

#p値
pnorm(Zu, lower.tail = FALSE)

これでwilcox_testのZ値、p値と同値になりました

前半で使用した　wilcox.test　を使用する場合には、correct =FALSE とすることで連続性の補正が省かれてマン・ホイットニーU検定を実行することができます

wilcox.test(
  data ~ group,
  data = dat_xy,
  alternative = "greater",
  correct =FALSE
)

ｐ値は同値となりました

タイがあるので正確検定ができず”警告メッセージ”が記載されますが問題ありません（correct =FALSEにしているので、このメッセージは不要かと思うのですが、、、）

適用

パッケージ

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

サンプルの作成

順位を追加

仮説（重要！）

Rでウィルコクソンの順位和検定

基準（reference）の設定

対立仮説1の場合（alternative = “greater”）

対立仮説2の場合（alternative = “less”）

対立仮説3の場合（alternative = “two.sided”）

P値を直接求める方法

対立仮説1の場合（A群＞B群）

対立仮説2の場合（A群＜B群）

対立仮説3の場合（A群$\ne$B群）

正規近似による方法

データの準備と要約

片側検定（帰無仮説Y群>X群）

wilcox_testの Z値 の求め方

X群とY群の順位和 $(W)$

Z値の求め方

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

コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください

wilcox_testの Z値の求め方

コメント欄『間違い』や『分かりにくい部分』などのご意見もお寄せください