- パッケージとデータの準備
- brunnermunzel.test( )を使ってBrunner-Munzel検定を行う
- brunnermunzel.permutation.test( )を使って並べ替えBrunner-Munzel検定を行う
- おわりに
- 参考資料
少し前に連続変数の代表値比較でどの検定方法を使えばよいのか考えた。
今回はそこで出てきたBrunner-Munzel検定をRで実行する方法を自分用に書いておく。
Brunner-Munzel検定を実行するための関数はlawstatパッケージにもあるみたいだけど、brunnermunzelパッケージだと小サンプルで使う並べ替え版（permutated BM test）も入っているみたい。
パッケージとデータの準備
brunnermunzelパッケージを読み込んで、模擬データを作る。このパッケージの説明で使われていた治療別（Y, N）の術後疼痛スケールのデータをそのまま使った。
library(brunnermunzel) Y <- c(1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1) N <- c(3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4)
データフレームを使う方法を確認するために、上のデータをデータフレーム型としておく。
dat <- data.frame(value = c(Y, N), group = factor(rep(c("Y", "N"), c(length(Y), length(N))), levels = c("Y", "N")) )
出来たデータフレームは次のとおり。
> head(dat) value group 1 1 Y 2 2 Y 3 1 Y 4 1 Y 5 1 Y 6 1 Y
brunnermunzel.test( )を使ってBrunner-Munzel検定を行う
比較したい変数の指定方法はt.test( )と同じ。他に指定できる主な引数は以下のとおり。
alternative："two.sided", "greater", "less"alpha：有意水準（デフォルトは0.05）perm：並べ替え法（permutated BM test）を使うかどうか。デフォルトは"FALSE"（並べ替え法を使わない）。est：何の信頼区間を推定するか。est = "original"のときは、を、
est = "difference"のときはを推定する。
全てデフォルトのまま実行してみる。
> brunnermunzel.test(Y,N) Brunner-Munzel Test data: Y and N Brunner-Munzel Test Statistic = 3.1375, df = 17.683, p-value = 0.005786 95 percent confidence interval: 0.5952169 0.9827052 sample estimates: P(X<Y)+.5*P(X=Y) 0.788961
P値は0.0058となった。
brunnermunzel.permutation.test( )を使って並べ替えBrunner-Munzel検定を行う
次に並べ替えBrunner-Munzel検定をやってみる。brunnermunzel.test( )の引数でperm=TRUEとしてもよいが、せっかくなのでbrunnermunzel.permutation.test( )を使ってみる。
データフレームを渡すときは、「比べたい値~群分け変数」の要領で指定すれば良い。このときは引数dataで使用するデータを指定するのを忘れないように。
estをデフォルトでない"difference"にしてみる。
brunnermunzel.permutation.test(value ~ group, data = dat, est = "difference")
結果は下のとおり。
permuted Brunner-Munzel Test data: value by group p-value = 0.008038 sample estimates: P(X<Y)-P(X>Y) 0.5779221
並べ替え版では信頼区間は計算されない。
サンプル数があまりにも多いと、すごく時間がかかってしまう。デフォルトでは40116600通り以上の組み合わせを計算しないといけないときは並べ替え法ではなく通常版を使う設定千種なっている（引数force=TRUEとすると組み合わせ数の上限がなくなり、強制的に並べ替え法を実行できる）。
おわりに
- そのうちMann-Whitney U検定に取って代わる日が来るんでしょうか。
参考資料
- 奥村先生の解説記事です。Brunner-Munzel検定の理論的背景を知りたい方は是非どうぞ。
- hoxo_mさんの代表値を検定方法を比較したシミュレーションの記事です。コードが分からなくても結果のグラフを見るだけでも面白かったです。
