Brunner-Munzel検定

分散は等しくないし、正規分布でもない

　これまで、一番制約の強いStundentのt検定から始まり、検定のための標本の前提条件を徐々に緩めながら検定法の紹介を行ってきた。2標本の検定の最後は、おそらくこれまで紹介した中でも、最も制約が弱い状態を仮定した検定を紹介する。Brunner-Munzel検定は、2標本に関して正規性を仮定しないだけでなく、等分散性も仮定しない（分布が同じであることを仮定しない、と紹介されていることもある）。制約が緩いので様々な状況に対応できる可能性を秘めているが、一方であまり有名な検定ではないので検定の特徴の精査は一部を除いてまだまだである。この項ではBrunner-Munzel検定の特徴を調べ、適用可能な標本の特徴について考えてみる。

　今回、2つのパターンの標本について考える。1つ目は正規分布 vs 二峰分布のように2標本の母集団の分布が全く異なっている場合、2つ目は2標本がコーシー分布に由来する場合である。2つ目の場合はコーシー分布が分散を定義できない分布であるため、2標本は分散が異なる状況になりやすい。上記のような紹介からすると、いずれの場合でも利用可能に思える。結論から言うと、今回の場合、Brunner-Munzel検定はコーシー vs コーシーには利用可能だが、正規 vs 二峰に対しての利用は避けたほうが良い。

検出力

　Brunner-Munzel検定の検定統計量BMは下記のように算出する。

Brunner-Munzel検定は検定統計量BMに関連するpという値に関して、p = 1/2、すなわちR1 = R2という帰無仮説を検定する。R1、R2は2標本を合わせたときの、各標本の平均順位和であり、帰無仮説は平均順位和が等しくなる、ということを主張している。結局のところ、pはXとYの標本から1つずつ値を取り出したとき、Xi > Yjとなる確率と考えることができる（逆でもいいが）。検定統計量BMは、上記で計算した自由度dfのt分布に従うことを利用して検定を行う。サンプルサイズが大きいときは、平均0、標準偏差1の正規分布（標準正規分布）に従うと仮定してもよい。

　では、シミュレーションで検出力を確認してみる。まずは、正規分布と二峰分布の標本間での検定である。先にどんな標本の様子になるか確認してみよう。

------------------------------------------------------

library(plotn)

library(lawstat)

#install.packages("lawstat", repos="http://cran.ism.ac.jp/")

me1 <- 50

sd1 <- 10

m11 <- rnorm(10000, mean = me1, sd = sd1)

m22 <- NULL

for(i in 1:10000){

t <- rbinom(1,1,1/3)

if(t == 0){

m22 <- c(m22, rnorm(1, 45, 14.5))

} else {

m22 <- c(m22, rnorm(1, 90, 11))

}

}

mean(m11)

## [1] 49.93302

sd(m11)

## [1] 10.00159

mean(m22)

## [1] 60.41132

sd(m22)

## [1] 25.09921

m <- data.frame(m = c(m11, m22), g = c(rep("n", 10000), rep("b",10000)))

histn(m ~ g, data = m, xlab = "value", ylab = "頻度") #図1の描画

------------------------------------------------------

今回は検証を行う目的でこれらの2標本を比較対象にするのであるが、このような標本が取れたときに、まず考えてほしいのはこれらの標本は本当に比較するべきか、ということである。特に二峰分布のほうは、もしかしたら分けられるべき標本が混じっている可能性も考えるべきである。なにも考えずにただ検定するのは良い手とは言えない。それはさておき、以下のようなコードで検証を行う。Brunner-Munzel検定以外にも、比較対象としてWelchのt検定とMann-WhitneyのU検定も行う。

------------------------------------------------------

#正規分布 vs 二峰分布

n <- 60

t2 <- NULL

u3 <- NULL

BM <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

　m1 <- rnorm(n, 50, 10) #平均50の正規分布の標本を作る。

m2 <- NULL

for(i in 1:n){ #平均60の二峰分布の標本を作る。

b <- rbinom(1,1,1/3)

if(b == 0){

m2 <- c(m2, rnorm(1, 55, 14.5))

} else {

m2 <- c(m2, rnorm(1, 100, 11))

}

}

ue <- ((n - 1)*(sd(m1)^2) + (n - 1)*(sd(m2)^2))/(2*n - 2)

t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n + 1/n))) #統計量tを算出

p2 <- c(p2, t.test(m1,m2,var.equal = F)$p.value) #Welchのt検定を行い、p値を取り出す

#以下Uの算出

d <- data.frame(

m = c(m1, m2),

sample = c(

rep(1, n),

rep(2, n)

)

)

d <- d[order(d$m),]

d$rank <- 1:nrow(d)

r1 <- sum(d[d$sample == 1,"rank"])

r2 <- sum(d[d$sample == 2,"rank"])

u1 <- n*n + n*(n + 1)/2 - r1

u2 <- n*n + n*(n + 1)/2 - r2

u3 <- c(u3, min(u1,u2))

#以上Uの算出

p3 <- c(p3, wilcox.test(m1,m2, exact = T)$p.value)

#以下BMの算出

rt1 <- rank(c(m1,m2))[1:n]

rt2 <- rank(c(m1,m2))[(n+1):(2*n)]

mrt1 <- mean(rt1)

mrt2 <- mean(rt2)

p <- (mrt1 - mrt2)/(2*n) + 1/2

r1 <- rank(m1)

r2 <- rank(m2)

mr1 <- mean(r1)

mr2 <- mean(r2)

S1 <- (1/(n-1)) * sum((rt1 - r1 - mrt1 + mr1)^2)

S2 <- (1/(n-1)) * sum((rt2 - r2 - mrt2 + mr2)^2)

BM <- c(BM, (n*n*(mrt2 - mrt1))/((n+n)*sqrt(n*S1+n*S2)))

#以上BMの算出

p4 <- c(p4, brunner.munzel.test(m1, m2)$p.value)

}

histn(t2, xlab = "統計量t", ylab = "頻度") #図2の描画

histn(p2, xlab = "P value", ylab = "頻度") #図3の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図4の描画

histn(p3, xlab = "P value", ylab = "頻度") #図5の描画

BM <- BM[!is.infinite(BM)]

histn(BM, xlab = "統計量BM", ylab = "頻度") #図6の描画

histn(p4, xlab = "P value", ylab = "頻度") #図7の描画

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 9967

p4 <- p4[!is.na(p4)] #BM検定はたまにうまく計算できずinfを返すことがあるのでそれを除く

sum(p4 < 0.05)

## [1] 9962

length(p4)

## [1] 10000

------------------------------------------------------

以上を確認すると、いずれの検定でも検出力にあまり差がないことがわかる。また、いずれの方法でも非常に検出力が高い。これを良いことと考えるかは危険率のところで議論しよう。

　次にコーシー分布の標本同士での比較を考える。Mann-WhitneyのU検定の検討の時と異なり、rcauchy()のscale引数を変えることで、外れ値の出やすさが変わっている。そうすると、分散も大きく変わるようになる。コードは以下の通り。

------------------------------------------------------

#コーシー分布

n <- 30

l1 <- 50

l2 <- 60

sc1 <- 1

sc2 <- 4

t2 <- NULL

u3 <- NULL

BM <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

m1 <- rcauchy(n, location = l1, scale = sc1)

m2 <- rcauchy(n, location = l2, scale = sc2)

ue <- ((n - 1)*(sd(m1)^2) + (n - 1)*(sd(m2)^2))/(2*n - 2)

t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n + 1/n))) #統計量tを算出

p2 <- c(p2, t.test(m1,m2,var.equal = F)$p.value) #Welchのt検定を行い、p値を取り出す

#以下Uの算出

d <- data.frame(

m = c(m1, m2),

sample = c(

rep(1, n),

rep(2, n)

)

)

d <- d[order(d$m),]

d$rank <- 1:nrow(d)

r1 <- sum(d[d$sample == 1,"rank"])

r2 <- sum(d[d$sample == 2,"rank"])

u1 <- n*n + n*(n + 1)/2 - r1

u2 <- n*n + n*(n + 1)/2 - r2

u3 <- c(u3, min(u1,u2))

#以上Uの算出

p3 <- c(p3, wilcox.test(m1,m2, exact = T)$p.value)

#以下BMの算出

rt1 <- rank(c(m1,m2))[1:n]

rt2 <- rank(c(m1,m2))[(n+1):(2*n)]

mrt1 <- mean(rt1)

mrt2 <- mean(rt2)

p <- (mrt1 - mrt2)/(2*n) + 1/2

r1 <- rank(m1)

r2 <- rank(m2)

mr1 <- mean(r1)

mr2 <- mean(r2)

S1 <- (1/(n-1)) * sum((rt1 - r1 - mrt1 + mr1)^2)

S2 <- (1/(n-1)) * sum((rt2 - r2 - mrt2 + mr2)^2)

BM <- c(BM, (n*n*(mrt2 - mrt1))/((n+n)*sqrt(n*S1+n*S2)))

#以上BMの算出

p4 <- c(p4, brunner.munzel.test(m1, m2)$p.value)

}

histn(t2, xlab = "統計量t", ylab = "頻度") #図8の描画

histn(p2, xlab = "P value", ylab = "頻度") #図9の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図10の描画

histn(p3, xlab = "P value", ylab = "頻度") #図11の描画

BM <- BM[!is.infinite(BM)]

histn(BM, xlab = "統計量BM", ylab = "頻度") #図12の描画

histn(p4, xlab = "P value", ylab = "頻度") #図13の描画

sum(p2 < 0.05)

## [1] 4539

sum(p3 < 0.05)

## [1] 9992

p4 <- p4[!is.na(p4)]

sum(p4 < 0.05)

## [1] 9979

length(p4)

## [1] 10000

------------------------------------------------------

以上を確認すると、Welchのt検定では著しく検出力が低下し、45%程度となる。そもそも、統計量tの分布がt分布とならないので、信用のある結果は出ない。一方で、Mann-WhitneyのU検定やBrunner-Munzel検定の検出力はほぼ100%であり、高水準である。

危険率

　続いて危険率について検討する。平均値が等しい状況を考えればよいので、以下のようなコードを書く。

------------------------------------------------------

#正規分布 vs 二峰分布

n <- 60

t2 <- NULL

u3 <- NULL

BM <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

　m1 <- rnorm(n, 60, 10) #平均60の正規分布の標本を作る。

m2 <- NULL

for(i in 1:n){ #平均60の二峰分布の標本を作る。

b <- rbinom(1,1,1/3)

if(b == 0){

m2 <- c(m2, rnorm(1, 55, 14.5))

} else {

m2 <- c(m2, rnorm(1, 100, 11))

}

}

ue <- ((n - 1)*(sd(m1)^2) + (n - 1)*(sd(m2)^2))/(2*n - 2)

t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n + 1/n))) #統計量tを算出

p2 <- c(p2, t.test(m1,m2,var.equal = F)$p.value) #Welchのt検定を行い、p値を取り出す

#以下Uの算出

d <- data.frame(

m = c(m1, m2),

sample = c(

rep(1, n),

rep(2, n)

)

)

d <- d[order(d$m),]

d$rank <- 1:nrow(d)

r1 <- sum(d[d$sample == 1,"rank"])

r2 <- sum(d[d$sample == 2,"rank"])

u1 <- n*n + n*(n + 1)/2 - r1

u2 <- n*n + n*(n + 1)/2 - r2

u3 <- c(u3, min(u1,u2))

#以上Uの算出

p3 <- c(p3, wilcox.test(m1,m2, exact = T)$p.value)

#以下BMの算出

rt1 <- rank(c(m1,m2))[1:n]

rt2 <- rank(c(m1,m2))[(n+1):(2*n)]

mrt1 <- mean(rt1)

mrt2 <- mean(rt2)

p <- (mrt1 - mrt2)/(2*n) + 1/2

r1 <- rank(m1)

r2 <- rank(m2)

mr1 <- mean(r1)

mr2 <- mean(r2)

S1 <- (1/(n-1)) * sum((rt1 - r1 - mrt1 + mr1)^2)

S2 <- (1/(n-1)) * sum((rt2 - r2 - mrt2 + mr2)^2)

BM <- c(BM, (n*n*(mrt2 - mrt1))/((n+n)*sqrt(n*S1+n*S2)))

#以上BMの算出

p4 <- c(p4, brunner.munzel.test(m1, m2)$p.value)

}

histn(t2, xlab = "統計量t", ylab = "頻度") #図13の描画

histn(p2, xlab = "P value", ylab = "頻度") #図14の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図15の描画

histn(p3, xlab = "P value", ylab = "頻度") #図16の描画

BM <- BM[!is.infinite(BM)]

histn(BM, xlab = "統計量BM", ylab = "頻度") #図17の描画

histn(p4, xlab = "P value", ylab = "頻度") #図18の描画

sum(p2 < 0.05)

## [1] 8164

sum(p3 < 0.05)

## [1] 4163

p4 <- p4[!is.na(p4)]

sum(p4 < 0.05)

## [1] 3752

length(p4)

## [1] 10000

------------------------------------------------------

さて、今回の結果は非常に芳しくない。Welchのt検定、Mann-WhitneyのU検定、Brunner-Munzel検定のそれぞれの危険率は、81.6%、41.6%、37.5%であり、いずれの検定も大幅に5%を上回ってしまっている。Welchのt検定は本当に壊滅的で、差がないにもかかわらず、80%以上の確率で差があると判断するということだ。このように全く分布の形が異なる標本同士での比較ではWelchのt検定は全く使えない。かといって、ノンパラメトリック検定のMann-WhitneyのU検定やBrunner-Munzel検定は、Welchのt検定よりはましであるが、決していい結果とはいえない。Brunner-Munzel検定はU検定よりもマシなようであるが、大幅に改善できているとは言えない。検出力の結果と併せて考えると、これらの検定は全く異なる分布で比較してしまうと、差があると判断してしまいやすい、ということである。ゆえに、全く異なる分布であるとわかった場合、これらのいずれの検定を使うことも推奨されない。分布の形が同じであれば二峰分布であっても、Mann-WhitneyのU検定やBrunner-Munzel検定を使うことはできる。前述したとおり、標本の片方が二峰分布であるとしたら、標本が混ざってないかよく検討しよう。

　幸い、2標本がともにコーシー分布に従う、つまり外れ値が出やすく、分散も違いやすい、という状況ならBrunner-Munzel検定は有効に使える。確認してみよう。検出力の時と同様にscale引数は変えてある。コードは以下の通り。

------------------------------------------------------

#コーシー分布

n <- 30

l1 <- 60

l2 <- 60

sc1 <- 1

sc2 <- 4

t2 <- NULL

u3 <- NULL

BM <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

m1 <- rcauchy(n, location = l1, scale = sc1)

m2 <- rcauchy(n, location = l2, scale = sc2)

ue <- ((n - 1)*(sd(m1)^2) + (n - 1)*(sd(m2)^2))/(2*n - 2)

t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n + 1/n))) #統計量tを算出

p2 <- c(p2, t.test(m1,m2,var.equal = F)$p.value) #Welchのt検定を行い、p値を取り出す

#以下Uの算出

d <- data.frame(

m = c(m1, m2),

sample = c(

rep(1, n),

rep(2, n)

)

)

d <- d[order(d$m),]

d$rank <- 1:nrow(d)

r1 <- sum(d[d$sample == 1,"rank"])

r2 <- sum(d[d$sample == 2,"rank"])

u1 <- n*n + n*(n + 1)/2 - r1

u2 <- n*n + n*(n + 1)/2 - r2

u3 <- c(u3, min(u1,u2))

#以上Uの算出

p3 <- c(p3, wilcox.test(m1,m2, exact = T)$p.value)

#以下BMの算出

rt1 <- rank(c(m1,m2))[1:n]

rt2 <- rank(c(m1,m2))[(n+1):(2*n)]

mrt1 <- mean(rt1)

mrt2 <- mean(rt2)

p <- (mrt1 - mrt2)/(2*n) + 1/2

r1 <- rank(m1)

r2 <- rank(m2)

mr1 <- mean(r1)

mr2 <- mean(r2)

S1 <- (1/(n-1)) * sum((rt1 - r1 - mrt1 + mr1)^2)

S2 <- (1/(n-1)) * sum((rt2 - r2 - mrt2 + mr2)^2)

BM <- c(BM, (n*n*(mrt2 - mrt1))/((n+n)*sqrt(n*S1+n*S2)))

#以上BMの算出

p4 <- c(p4, brunner.munzel.test(m1, m2)$p.value)

}

histn(t2, xlab = "統計量t", ylab = "頻度") #図19の描画

histn(p2, xlab = "P value", ylab = "頻度") #図20の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図21の描画

histn(p3, xlab = "P value", ylab = "頻度") #図22の描画

BM <- BM[!is.infinite(BM)]

histn(BM, xlab = "統計量BM", ylab = "頻度") #図23の描画

histn(p4, xlab = "P value", ylab = "頻度") #図24の描画

sum(p2 < 0.05)

## [1] 216

sum(p3 < 0.05)

## [1] 613

p4 <- p4[!is.na(p4)]

sum(p4 < 0.05)

## [1] 500

length(p4)

## [1] 10000

------------------------------------------------------

Welchのt検定があてにならないのは今まで通りであるが、Mann-WhitneyのU検定とBrunner-Munzel検定はいい水準を保っている。この2つの検定で大きな差はないが、scaleが異なるコーシー分布同士であればBrunner-Munzel検定のほうがやや精度が良いことがわかる。

Rで行うBrunner-Muzel検定

　上記のような背景のもと、Brunner-Munzel検定を行ってみよう。以下のようなサンプルサイズ20の独立ではないデータがあるとする。有意水準は5%としよう。まずは、データを描画してみて様子を確認する。

------------------------------------------------------

m1 <- c(49, 50, 55, 49, 50, 45, 49, 49, 20, 49,

51, 56, 44, 49, 49, 57, 51, 47, 52, 49,

50, 54, 50, 49, 47, 49, 49, 59, 47, 50)

m2 <- c(60, 58, 51, 56, 60, 60, 59, 75, 62, 58,

57, 61, 62, 52, 65, 111, 56, 53, 71, 54,

75, 58, 68, 53, 55, 6, 64, 29, 73, 57)

n1 <- length(m1)

n2 <- length(m2)

sd1 <- sd(m1)

sd2 <- sd(m2)

sd1

## [1] 6.404381

sd2

## [1] 16.27808

d <- data.frame(d = c(m1, m2), g = c(rep("m1", n1), rep("m2", n2)))

boxplotn(d$d ~ d$g, ylab = "data", jitter.method = "swarm") #図25の描画

------------------------------------------------------

分散が異なり、外れ値が存在しそうなので、Brunner-Munzel検定を行う。lawstatパッケージのbrunner.munzel.test()関数を使おう。

------------------------------------------------------

fit <- brunner.munzel.test(m1, m2)

fit

##

## Brunner-Munzel Test

##

## data: m1 and m2

## Brunner-Munzel Test Statistic = 7.8857, df = 39.954, p-value =

## 1.139e-09

## 95 percent confidence interval:

## 0.7883879 0.9871677

## sample estimates:

## P(X<Y)+.5*P(X=Y)

## 0.8877778

------------------------------------------------------

Brunner-Munzel Test Statisticは検定統計量のことで以下のように計算する。

------------------------------------------------------

rt1 <- rank(c(m1,m2))[1:n1]

rt2 <- rank(c(m1,m2))[(n1+1):(n1+n2)]

mrt1 <- mean(rt1)

mrt2 <- mean(rt2)

r1 <- rank(m1)

r2 <- rank(m2)

mr1 <- mean(r1)

mr2 <- mean(r2)

S1 <- (1/(n1-1)) * sum((rt1 - r1 - mrt1 + mr1)^2)

S2 <- (1/(n2-1)) * sum((rt2 - r2 - mrt2 + mr2)^2)

BM0 <- (n1*n2*(mrt2 - mrt1))/((n1+n2)*sqrt(n1*S1+n2*S2))

BM0

## [1] 7.885673

------------------------------------------------------

dfは自由度のことで、検定統計量をt分布に近似するときの計算に用いる。以下のように計算する。

------------------------------------------------------

df <- ((n1*S1+n2*S2)^2)/(((n1*S1)^2)/(n1-1)+((n2*S2)^2)/(n2-1))

df

## [1] 39.95415

------------------------------------------------------

p valueは検定統計量が自由度dfのt分布に従うと考えて算出されており、以下のように計算する。

------------------------------------------------------

pt(df = df, q = BM0, lower.tail = F)*2

## [1] 1.13922e-09

------------------------------------------------------

95 percent confidence intervalは検定統計量の95%信頼区間ではなく、p（p値ではない）の信頼区間である。検定統計量がt分布に従うことと仮定して、以下のように計算する。下記の計算は、上記の検定法に記した検定統計量BMとpの関係BM = (p - 1/2) × nm/sqrt(nΣ1 + mΣ2)から、p = BM × sqrt(nΣ1 + mΣ2)/nm + 1/2となるので、BMの値をt分布の95%信頼区間で置き換えればよい。

------------------------------------------------------

(BM0 - qt(df = df, p = 0.975))/(n1*n2/sqrt(n1*S1+n2*S2)) + 1/2

## [1] 0.7883879

BM0/(n1*n2/sqrt(n1*S1+n2*S2)) + 1/2

## [1] 0.8877778

(BM0 + qt(df = df, p = 0.975))/(n1*n2/sqrt(n1*S1+n2*S2)) + 1/2

## [1] 0.9871677

------------------------------------------------------

最後のP(X<Y)+.5*P(X=Y) がpのことであり、下記のように計算することができる。直上の2番目の計算と一致する。

------------------------------------------------------

p <- (mrt2 - mrt1)/(n1+n2) + 1/2

## [1] 0.8877778

------------------------------------------------------

今回はBrunner-Munzel検定の関数が明確にt分布を利用していることで、一致したp値を計算できた。念のため、これまで通り、検定統計量BMが正規分布近似やt分布近似が可能であることを確認しておこう。

------------------------------------------------------

n1 <- 30

n2 <- 30

BM <- NULL

for(i in 1:300000){

m1 <- rnorm(n1, mean = 0, sd = 20)

m2 <- rnorm(n2, mean = 0, sd = 20)

rt1 <- rank(c(m1,m2))[1:n1]

rt2 <- rank(c(m1,m2))[(n1+1):(n1+n2)]

mrt1 <- mean(rt1)

mrt2 <- mean(rt2)

r1 <- rank(m1)

r2 <- rank(m2)

mr1 <- mean(r1)

mr2 <- mean(r2)

S1 <- (1/(n1-1)) * sum((rt1 - r1 - mrt1 + mr1)^2)

S2 <- (1/(n2-1)) * sum((rt2 - r2 - mrt2 + mr2)^2)

BM <- c(BM, (n1*n2*(mrt2 - mrt1))/((n1+n2)*sqrt(n1*S1+n2*S2)))

df <- ((n1*S1+n2*S2)^2)/(((n1*S1)^2)/(n1-1)+((n2*S2)^2)/(n2-1))

}

histn(BM, freq = F, xlab = "統計量BM") #図26の描画

xx <- seq(-10,10, length = 200)

yy1 <- dnorm(xx, mean = 0, sd = 1)

yy2 <- dt(xx, df = df)

overdraw(points(xx, yy1, type = "l", col = "red"),

points(xx, yy2, type = "l", col = "blue"))

------------------------------------------------------

サンプルサイズが大きいときは正規近似でもt分布近似でも問題がない。一方でサンプルサイズが小さいときは正規近似よりもt分布近似のほうが精度が良い。

------------------------------------------------------

n1 <- 5

n2 <- 5

BM <- NULL

for(i in 1:300000){

m1 <- rnorm(n1, mean = 0, sd = 20)

m2 <- rnorm(n2, mean = 0, sd = 20)

rt1 <- rank(c(m1,m2))[1:n1]

rt2 <- rank(c(m1,m2))[(n1+1):(n1+n2)]

mrt1 <- mean(rt1)

mrt2 <- mean(rt2)

r1 <- rank(m1)

r2 <- rank(m2)

mr1 <- mean(r1)

mr2 <- mean(r2)

S1 <- (1/(n1-1)) * sum((rt1 - r1 - mrt1 + mr1)^2)

S2 <- (1/(n2-1)) * sum((rt2 - r2 - mrt2 + mr2)^2)

BM <- c(BM, (n1*n2*(mrt2 - mrt1))/((n1+n2)*sqrt(n1*S1+n2*S2)))

df <- ((n1*S1+n2*S2)^2)/(((n1*S1)^2)/(n1-1)+((n2*S2)^2)/(n2-1))

}

BM <- BM[!is.infinite(BM)]

histn(BM, freq = F, xlab = "統計量BM") #図27の描画

xx <- seq(-5,5, length = 200)

yy1 <- dnorm(xx, mean = 0, sd = 1)

yy2 <- dt(xx, df = df)

overdraw(points(xx, yy1, type = "l", col = "red"),

points(xx, yy2, type = "l", col = "blue"))

------------------------------------------------------

まとめ

　以上の検討から、上記の条件ではBrunner-Munzel検定はU検定よりも汎用性が高そうだ。とはいえ、どんな標本にでも活用できる、というほどではない。もし、そんな検定があるとしたら、とっくの昔に広まって、みんな使っているだろう。汎用性が高いBrunner-Munzel検定も、ちゃんと使いどころがあるのだ。それさえわきまえれば、非常に強力なツールとなる。