Brunner-Munzel 検定

Brunner-Munzel 検定は、対応のない 2 群の検定法であり、正規性を仮定しないノンパラメトリックな手法です。
同様のものにマン・ホイットニーの U 検定がありますが、マン・ホイットニーの U 検定は、等分散性を仮定しているため、不等分散の場合に検定がうまくいかないという問題があります。
それに対して、Brunner-Munzel 検定は等分散性を仮定しないため、マン・ホイットニーの U 検定よりも多くの場面で使用できる検定手法と言えます。

R で Brunner-Munzel 検定を行うには、lawstat パッケージの brunner.munzel.test() 関数を使います。

library(lawstat)
x = c(1,2,1,1,1,1,1,1,1,1,2,4,1,1)
y = c(3,3,4,3,1,2,3,1,1,5,4)
brunner.munzel.test(x,y)
## 
##  Brunner-Munzel Test
## 
## data:  x and y
## Brunner-Munzel Test Statistic = 3.138, df = 17.68, p-value =
## 0.005786
## 95 percent confidence interval:
##  0.5952 0.9827
## sample estimates:
## P(X<Y)+.5*P(X=Y) 
##            0.789

本記事では、不等分散の状況で Brunner-Munzel 検定がどれほどの検定精度を持つのかを調査します。
その際、比較対象として、t検定、Welch の t検定、マン・ホイットニーの U 検定についても検定精度を調査します。

それぞれの検定における仮定は次の通りです。

手法 正規性 等分散性
Student’s t ○ ○
Welch’s t ○ ×
Mann-Whitney × ○
Brunner-Munzel × ×

1. 正規分布

まずは、正規分布に対する検定精度を調査します。 標準偏差 1 の正規分布に対して、標準偏差を \(\frac{1}{4}\), \(\frac{1}{2}\), \(\frac{3}{4}\), 1, \(\frac{4}{3}\), 2, 4 と変化させた正規分布との検定をそれぞれ繰り返し、第一種の過誤(alpha error rate)を調査します。

sigma_ratio <- c("1/4", "1/2", "3/4", "1", "4/3", "2", "4")
sigma1 <- c(1/4, 1/2, 3/4, 1, 4/3, 2, 4)
sigma2 <- 1
iter_num <- 2000

library(pforeach)
library(lawstat)
pforeach(sigma=sigma1, .c=rbind)({
  npforeach(i=seq_len(iter_num), .c=rbind)({
    x1 <- rnorm(15, mean = 5, sd = sigma)
    x2 <- rnorm(45, mean = 5, sd = sigma2)
    student <- t.test(x1, x2, var.equal = TRUE)$p.value <= 0.05
    welch <- t.test(x1, x2, var.equal = FALSE)$p.value <= 0.05
    MH <- wilcox.test(x1,x2)$p.value <= 0.05
    BM <- brunner.munzel.test(x1, x2)$p.value <= 0.05
    data.frame(`Student's t`=student, `Welch's t`=welch, `Mann-Whitney`=MH, `Brunner-Munzel`=BM, check.names = FALSE)
  }) -> res
  colSums(res)/iter_num
}) -> result

library(tidyr)
library(dplyr)
data <- result %>% data.frame(check.names = FALSE) %>% 
  cbind(sigma=factor(sigma_ratio, levels=sigma_ratio)) %>% 
  gather(`Test Method`, `Alpha Error Rate`, 1:4)

library(ggplot2)
ggplot(data, aes(x=sigma, y=`Alpha Error Rate`, group=`Test Method`, color=`Test Method`)) + 
  geom_line() + geom_point() + ggtitle("Normal Dist.")

plot of chunk unnamed-chunk-3

有意水準 0.05 で検定を行っているため、alpha error rate はどれも 0.05 になるはずですが、異なる標準偏差に対する検定において、t検定とマン・ホイットニーは 0.05 から大きく外れており、検定精度が非常に悪いという結果になりました。
対して、Welch と Brunner-Munzel は異なる標準偏差に対しても alpha error rate は 0.05 を保っています。

2. 連続一様分布

次に、連続一様分布に対して、正規分布と同様に標準偏差比を変えて検定を繰り返し、第一種の過誤について調査しました。

sigma_ratio <- c("1/4", "1/2", "3/4", "1", "4/3", "2", "4")
sigma1 <- c(1/4, 1/2, 3/4, 1, 4/3, 2, 4)
a1 <- 5 - sqrt(3) * sigma1
b1 <- 5 + sqrt(3) * sigma1
a2 <- 5 - sqrt(3)
b2 <- 5 + sqrt(3)
iter_num <- 2000

library(pforeach)
library(lawstat)
pforeach(a=a1, b=b1, .c=rbind)({
  npforeach(i=seq_len(iter_num), .c=rbind)({
    x1 <- runif(15, min = a, max = b)
    x2 <- runif(45, min = a2, max = b2)
    student <- t.test(x1, x2, var.equal = TRUE)$p.value <= 0.05
    welch <- t.test(x1, x2, var.equal = FALSE)$p.value <= 0.05
    MH <- wilcox.test(x1,x2)$p.value <= 0.05
    BM <- brunner.munzel.test(x1, x2)$p.value <= 0.05
    data.frame(`Student's t`=student, `Welch's t`=welch, `Mann-Whitney`=MH, `Brunner-Munzel`=BM, check.names = FALSE)
  }) -> res
  colSums(res)/iter_num
}) -> result

library(tidyr)
library(dplyr)
data <- result %>% data.frame(check.names = FALSE) %>% 
  cbind(sigma=factor(sigma_ratio, levels=sigma_ratio)) %>% 
  gather(`Test Method`, `Alpha Error Rate`, 1:4)

library(ggplot2)
ggplot(data, aes(x=sigma, y=`Alpha Error Rate`, group=`Test Method`, color=`Test Method`)) + 
  geom_line() + geom_point() + ggtitle("Uniform Dist.")

plot of chunk unnamed-chunk-4

一様分布に対する結果は、正規分布のときと同様になりました。
異なる標準偏差に対しては、t検定とマン・ホイットニーは検定精度が非常に悪く、Welch と Brunner-Munzel は良い精度を保っています。

3. 対数正規分布(平均値)

上記二つの分布は対称分布でしたが、歪んだ分布として対数正規分布について調査しました。
対数正規分布は平均値と中央値が異なるため、まずは同じ平均値を持つ対数正規分布について検定を繰り返し、第一種の過誤を調査しました。

sigma_ratio <- c("1/4", "1/2", "3/4", "1", "4/3", "2", "4")
sigma1 <- c(1, 2, 3, 4, 4, 4, 4)
sigma2 <- c(4, 4, 4, 4, 3, 2, 1)
sigmas1 <- sqrt(log((sigma1^2)/(5^2) + 1))
sigmas2 <- sqrt(log((sigma2^2)/(5^2) + 1))
mus1 <- log(5) - (sigmas1^2)/2
mus2 <- log(5) - (sigmas2^2)/2
iter_num <- 2000

library(pforeach)
library(lawstat)
pforeach(mu1=mus1, sigma1=sigmas1, mu2=mus2, sigma2=sigmas2, .c=rbind)({
  npforeach(i=seq_len(iter_num), .c=rbind)({
    x1 <- rlnorm(15, meanlog = mu1, sdlog = sigma1)
    x2 <- rlnorm(45, meanlog = mu2, sdlog = sigma2)
    student <- t.test(x1, x2, var.equal = TRUE)$p.value <= 0.05
    welch <- t.test(x1, x2, var.equal = FALSE)$p.value <= 0.05
    MH <- wilcox.test(x1,x2)$p.value <= 0.05
    BM <- brunner.munzel.test(x1, x2)$p.value <= 0.05
    data.frame(`Student's t`=student, `Welch's t`=welch, `Mann-Whitney`=MH, `Brunner-Munzel`=BM, check.names = FALSE)
  }) -> res
  colSums(res)/iter_num
}) -> result

library(tidyr)
library(dplyr)
data <- result %>% data.frame(check.names = FALSE) %>% 
  cbind(sigma=factor(sigma_ratio, levels=sigma_ratio)) %>% 
  gather(`Test Method`, `Alpha Error Rate`, 1:4)

library(ggplot2)
ggplot(data, aes(x=sigma, y=`Alpha Error Rate`, group=`Test Method`, color=`Test Method`)) + 
  geom_line() + geom_point() + ggtitle("Log Normal Dist.(Mean)")

plot of chunk unnamed-chunk-5

異なる標準偏差に対して、t検定とマン・ホイットニーの精度が悪いのは今までと同様ですが、Brunner-Munzel も精度が悪くなっています。
意外なのは、正規性を仮定しているはずの Welch が善戦しているということです。

4. 対数正規分布(中央値)

次は、同じ中央値を持つ対数正規分布について、標準偏差比を変更して検定を繰り返し、第一種の過誤を調査します。

sigma_ratio <- c("1/4", "1/2", "3/4", "1", "4/3", "2", "4")
sigma1 <- c(1, 2, 3, 4, 4, 4, 4)
sigma2 <- c(4, 4, 4, 4, 3, 2, 1)
sigmas1 <- sqrt(log((1 + sqrt((4 * sigma1^2)/25 + 1))/2))
sigmas2 <- sqrt(log((1 + sqrt((4 * sigma2^2)/25 + 1))/2))
mu1 <- log(5)
mu2 <- log(5)
iter_num <- 2000

library(pforeach)
library(lawstat)
pforeach(sigma1=sigmas1, sigma2=sigmas2, .c=rbind)({
  npforeach(i=seq_len(iter_num), .c=rbind)({
    x1 <- rlnorm(15, meanlog = mu1, sdlog = sigma1)
    x2 <- rlnorm(45, meanlog = mu2, sdlog = sigma2)
    student <- t.test(x1, x2, var.equal = TRUE)$p.value <= 0.05
    welch <- t.test(x1, x2, var.equal = FALSE)$p.value <= 0.05
    MH <- wilcox.test(x1,x2)$p.value <= 0.05
    BM <- brunner.munzel.test(x1, x2)$p.value <= 0.05
    data.frame(`Student's t`=student, `Welch's t`=welch, `Mann-Whitney`=MH, `Brunner-Munzel`=BM, check.names = FALSE)
  }) -> res
  colSums(res)/iter_num
}) -> result

library(tidyr)
library(dplyr)
data <- result %>% data.frame(check.names = FALSE) %>% 
  cbind(sigma=factor(sigma_ratio, levels=sigma_ratio)) %>% 
  gather(`Test Method`, `Alpha Error Rate`, 1:4)

library(ggplot2)
ggplot(data, aes(x=sigma, y=`Alpha Error Rate`, group=`Test Method`, color=`Test Method`)) + 
  geom_line() + geom_point() + ggtitle("Log Normal Dist.(Median)")

plot of chunk unnamed-chunk-6

今度は、Welch の精度が悪く、Brunner-Munzel の精度が良くなりました。

5. まとめ

4つの検定手法に対して、不等分散の状況下で検定精度を調査しました。

t検定は不等分散の状況では全く使い物にならないことが分かりました。
2群のノンパラメトリック検定で有名なマン・ホイットニーの U 検定も、不等分散の状況では全く使えないことが分かりました。

Welch の t検定については、不等分散の状況において、対称分布に対しては良い精度を保ちました。
また、歪んだ分布に対しても、平均値の検定としてはある程度は精度を保ちました。

Brunner-Munzel 検定については、不等分散の状況において、対称分布に対して良い精度を保ちました。
また、歪んだ分布に対しては、中央値の検定として良い精度を保ちました。

6. 参考