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

正規近似による方法
マン・ホイットニーのU検定

正規近似による方法

データの準備と要約

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))
df <- data.frame(data, group)
print(df)

   data group
1    12     X
2     4     X
3     9     X
4     6     X
5     2     X
6    10     X
7    12     X
8    11     X
9     7     X
10   10     X
11   10     Y
12   10     Y
13   25     Y
14   13     Y
15   18     Y
16    7     Y
17   15     Y
18    9     Y
19   18     Y
20   19     Y

順位を追加

# Add a ranking column (rank) to the data
df2 <- df %>%
  mutate(rank = rank(data, ties.method = "average"))
# Sort by rank
df3 <- df2 %>%
  arrange(rank)
print(df3)

   data group rank
1     2     X  1.0
2     4     X  2.0
3     6     X  3.0
4     7     X  4.5
5     7     Y  4.5
6     9     X  6.5
7     9     Y  6.5
8    10     X  9.5
9    10     X  9.5
10   10     Y  9.5
11   10     Y  9.5
12   11     X 12.0
13   12     X 13.5
14   12     X 13.5
15   13     Y 15.0
16   15     Y 16.0
17   18     Y 17.5
18   18     Y 17.5
19   19     Y 19.0
20   25     Y 20.0

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

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

groupをファクター変数に変更してYを基準に変更

df$group <- factor(df$group)

yを基準にして片側検定

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

        Asymptotic Wilcoxon-Mann-Whitney Test

data:  data by group (Y, X)
Z = 2.2798, p-value = 0.01131
alternative hypothesis: true mu is greater than 0

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

wilcox_testの Z値の求め方

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

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

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

# A tibble: 2 × 2
  group sum_rank
  <chr>    <dbl>
1 X           75
2 Y          135

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]}}$

# Calculate the sum of ranks for each group
Wx = 75
Wy = 135
# Sample sizes
nx <- sum(group == "X")
ny <- sum(group == "Y")
n <- nx + ny
# Calculation of U statistic
# Using group Y as the reference, so here we use Wy
# Expected value
EW <- ny*(n+1)/2
# Variance
VW <-  nx*ny*(n+1)/12
# Calculation of Z value
Zw <- (Wy - EW) / sqrt(VW)
# Output the result
Zw

> Zw
[1] 2.267787

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

> pnorm(Zw, lower.tail = FALSE)
[1] 0.0116711

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]}}$

# Calculate the sum of ranks for each group
Wx = 75
Wy = 135

# Sample sizes
nx <- sum(group == "X")
ny <- sum(group == "Y")
n <- nx + ny

# Calculation of U statistic
# Using group Y as the reference, so here we use Wy
U <- nx*ny + ny*(ny+1)/2 - Wy

# Expected value
EU <- nx*ny / 2

# Count of tied data
t <- table(df3$rank) 
# Correction for ties
T <- sum(t^3 - t) 
# Variance
VU <- (nx*ny / 12) * ((n+1) - T / (n*(n-1)))

# Calculation of Z value
Zu <- abs(U - EU) / sqrt(VU)

# Output the result
Zu

> Zu
[1] 2.279818

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

> pnorm(Zu, lower.tail = FALSE)
[1] 0.01130925

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

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

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

        Wilcoxon rank sum test

data:  data by group
W = 80, p-value = 0.01131
alternative hypothesis: true location shift is greater than 0

警告メッセージ:
wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...) で:
  タイがあるため、正確な p 値を計算することができ

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