Brunner-Munzel 検定（モンテカルロ法）

Brunner-Munzel 検定 ッテ，もはや，参照する価値もない記事です。パチもん検定です。

 

いくら，並べ替え検定が正確（？）といっても，サンプルサイズによる縛りがあるので，いつでも使えるわけではない。

そんなとき，サンプルサイズが大きい場合でも，任意の並べ替えに基づく検定を何回も行い，そのときの検定統計量と観察データにおける検定統計量の関係から完全な並べ替え検定の近似値を求めるのがモンテカルロ法（まあ，円周率を求めるときに，[0,1) の一様乱数を2組発生させて，原点からの距離が 1 未満になる確率を求めるようなことと同じ）

対象とする検定統計量は lawstat パッケージの brunner.munzel.test 関数の戻り値 estimate つまり P(X<Y)+.5*P(X=Y)

しつこいけど，t 値を対象としてはいけない。

bm.perm = function(x, y, trials=5000) {
  s = c(brunner.munzel.test(x, y)$estimate, brunner.munzel.test(y, x)$estimate)
  s.lo = min(s) + .Machine$double.eps
  s.hi = max(s) - .Machine$double.eps
  xy = c(x, y)
  nx = length(x)
  ny = length(y)
  nxy = nx + ny
  results = replicate(trials, {
    xy = sample(xy, nxy)
    brunner.munzel.test(xy[1:nx], xy[nx+1:ny])$estimate
  })
  mean(results <= s.lo | s.hi <= results)
}

これを使って，モンテカルロ法による p 値の推定値と，brunner.munzel.test の戻り値の p.value を比較する。

状況としては，第 2 群のサンプルサイズが 第 1 群のサンプルサイズの 2 倍 で 標準偏差も 2 倍

m = 200
p1 = p2 = numeric(m)
nx = 100; meanx = 50; sdx = 10
ny = 200; meany = 50; sdy = 20
for (i in 1:m) {
  cat(i, " ")
  x = rnorm(nx, meanx, sdx)
  y = rnorm(ny, meany, sdy)
  p1[i] = brunner.munzel.test(x, y)$p.value
  #t.test(x, y)$p.value
  p2[i] = bm.perm(x, y)
}
plot(p1, p2, xlab="brunner.munzel.test", ylab="permutation test", pch=19, cex=0.7, col="#66000020")
abline(0, 1, col=2)
abline(v=0.05, h=0.05, col=3, lty=2)

結果は以下の図のようになる。brunner.munzel.test()$p.value は，モンテカルロ法による p 値より小さめである。brunner.munzel.test()$p.value ＜ 0.05 となるのは 200 回中 9 回であるが，モンテカルロ法の場合には 200 回中 4  回だった。

trials=10000, m=1000 でやってみたら，同様な傾向であった。brunner.munzel.test()$p.value ＜ 0.05 となるのは 1000 回中 59 回であるが，モンテカルロ法の場合には 1000 回中 40  回だった。