Brunner-Munzel test with confidence interval estimation
@fronori
July 17, 2016
{lawstat}パッケージのbrunner.munzel.test()関数でBrunner Munzel検定を行う場合
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.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この例はLumley (1996) のデータだそうです。
{nparcomp}パッケージのnpar.t.test()関数を使ってBrunner Munzel検定を行う場合
library(nparcomp)npar.t.test()関数を用いるには、引数にデータフレームを指定する必要があります。１列目にスコア（計測値など）、 ２列目に群（“male”/“female”などでもいいし、0/1などのダミー変数でもよい）のデータフレームを作ります。
dat.Lumley <- data.frame(score=c(x,y), group=c(rep("group1", length(x)), rep("group2", length(y))))データフレームの最初の３行を表示
head(dat.Lumley, n=3)## score group
## 1 1 group1
## 2 2 group1
## 3 1 group1データフレームの最後の３行を表示
tail(dat.Lumley, n=3)## score group
## 23 1 group2
## 24 5 group2
## 25 4 group2箱ヒゲ図と蜂群図で可視化
par(mfrow=c(1,2))
boxplot(score~group, data=dat.Lumley)
library(beeswarm)
beeswarm(score~group, data=dat.Lumley)
npar.t.test()関数のオプションでmethod=“t.app”
を指定すると、B-M検定を行ってくれます。
a <- npar.t.test(score~group, data = dat.Lumley, alternative = "two.sided",
method="t.app", info=FALSE, plot.simci=FALSE)後でサマリーを表示するので、info=FALSE。後でプロットさせるので、plot.simci=FALSE
ちなみに、method = “permu”を指定すると、permutation testによる結果を計算してくれます。
より詳しい説明は公式ドキュメント参照。 https://cran.r-project.org/web/packages/nparcomp/nparcomp.pdf
検定結果のサマリーを表示
summary(a)##
## #-----Nonparametric Test Procedures and Confidence Intervals for relative effects-----#
##
## - Alternative Hypothesis: True relative effect p is less or equal than 1/2
## - Confidence level: 95 %
## - Method = Brunner - Munzel - T - Approx with 17.683 DF
## #---------------------------Interpretation---------------------------------------------#
## p(a,b) > 1/2 : b tends to be larger than a
## #--------------------------------------------------------------------------------------#
##
## #----Data Info-------------------------------------------------------------------------#
## Sample Size
## group1 group1 14
## group2 group2 11
##
## #----Analysis--------------------------------------------------------------------------#
## Effect Estimator Lower Upper T p.Value
## 1 p(group1,group2) 0.789 0.595 0.983 3.137 0.006“p(a,b) > 1/2 : b tends to be larger than a”" という結果から、二番目の群、この場合はy（group2)の方が大きいということが言えます。
信頼区間をプロット
plot(a)
brunner.munzel.test()とnpar.t.test()のどちらの関数を使っても、Brunner-Munzel検定で推定している確率はP(X<Y)+.5*P(X=Y) = 0.789、95%信頼区間 = (0.595, 0.983)となり、結果は一致してますね！
ちなみに、対応のあるデータの場合
nparcomp::npar.t.test.paired()で計算できます。