Embedded Files
目次
目次
・検出力
・危険率
・まとめ
・Brunner-Muzel検定の統計量の導出（工事中）
分散は等しくないし、正規分布でもない
分散は等しくないし、正規分布でもない
これまで、一番制約の強いStundentのt検定から始まり、検定のための標本の前提条件を徐々に緩めながら検定法の紹介を行ってきた。2標本の検定の最後は、おそらくこれまで紹介した中でも、最も制約が弱い状態を仮定した検定を紹介する。Brunner-Munzel（ブルンナー・ムンツェル）検定（Brunner and Muzel, 2000, Biometrical Journal）は、2標本に関して正規性を仮定しないだけでなく、等分散性も仮定しない（分布が同じであることを仮定しない、と紹介されていることもある）。制約が緩いので様々な状況に対応できる可能性を秘めているが、一方であまり有名な検定ではないので検定の特徴の精査は一部を除いてまだまだである。この項ではBrunner-Munzel検定の特徴を調べ、適用可能な標本の特徴について考えてみる。
今回、2つのパターンの標本について考える。1つ目は正規分布 vs 二峰分布のように2標本の母集団の分布が全く異なっている場合、2つ目は2標本がコーシー分布に由来する場合である。2つ目の場合はコーシー分布が分散を定義できない分布であるため、2標本は分散が異なる状況になりやすい。上記のような紹介からすると、いずれの場合でも利用可能に思える。結論から言うと、今回の場合、Brunner-Munzel検定はコーシー vs コーシーには利用可能だが、正規 vs 二峰に対しての利用は避けたほうが良い。
Brunner-Munzel検定
Brunner-Munzel検定
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の正規分布（標準正規分布）に従うと仮定してもよい。
検出力
検出力
・正規分布 vs 二峰分布
・正規分布 vs 二峰分布
では、シミュレーションで検出力を確認してみる。まずは、正規分布と二峰分布の標本間での検定である。先にどんな標本の様子になるか確認してみよう。
------------------------------------------------------
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の描画
------------------------------------------------------
図1 比較対象の2標本。橙は正規分布の標本m11、青は二峰分布の標本m22を示す。
今回は検証を行う目的でこれらの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
------------------------------------------------------
図2 正規分布 vs 二峰分布のWelchのt分布
図3 正規分布 vs 二峰分布のWelchのt検定のp値の分布
図4 正規分布 vs 二峰分布のMann-WhitneyのU分布
図5 正規分布 vs 二峰分布のMann-WhitneyのU検定のp値の分布
図6 正規分布 vs 二峰分布のBrunner-MunzelのBM分布
図7 正規分布 vs 二峰分布のBrunner-Munzel検定のp値の分布
以上を確認すると、いずれの検定でも検出力にあまり差がないことがわかる。また、いずれの方法でも非常に検出力が高い。これを良いことと考えるかは危険率のところで議論しよう。
・コーシー分布 vs コーシー分布
・コーシー分布 vs コーシー分布
次にコーシー分布の標本同士での比較を考える。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
------------------------------------------------------
図8 コーシー分布同士のWelchのt分布
図9 コーシー分布同士のWelchのt検定のp値の分布
図10 コーシー分布同士のMann-WhitneyのU分布
図11 コーシー分布同士のMann-WhitneyのU検定のp値の分布
図12 コーシー分布同士のBrunner-MunzelのBM分布
図13 コーシー分布同士のBrunner-Munzel検定のp値の分布
以上を確認すると、Welchのt検定では著しく検出力が低下し、45%程度となる。そもそも、統計量tの分布がt分布とならないので、信用のある結果は出ない。一方で、Mann-WhitneyのU検定やBrunner-Munzel検定の検出力はほぼ100%であり、高水準である。
危険率
危険率
・正規分布 vs 二峰分布
・正規分布 vs 二峰分布
続いて危険率について検討する。平均値が等しい状況を考えればよいので、以下のようなコードを書く。確認したいことがあるので、データを順位に変換したのち、Welchのt検定も行う。
------------------------------------------------------
#正規分布 vs 二峰分布
n <- 60
t2 <- NULL
tr <- NULL
u3 <- NULL
BM <- NULL
p2 <- NULL
p3 <- NULL
p4 <- NULL
p5 <- 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)
m1r <- d[d$sample == 1, "rank"]
m2r <- d[d$sample == 2, "rank"]
ue <- ((n - 1)*(sd(m1r)^2) + (n - 1)*(sd(m2r)^2))/(2*n - 2)
tr <- c(tr, (mean(m1r) - mean(m2r))/sqrt(ue*(1/n + 1/n))) #統計量tを算出
p3 <- c(p3, t.test(m1r,m2r,var.equal = F)$p.value) #Welchのt検定を行い、p値を取り出す
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の算出
p4 <- c(p4, 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の算出
p5 <- c(p5, 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(p5, xlab = "P value", ylab = "頻度") #図18の描画
sum(p2 < 0.05)
## [1] 8164
sum(p3 < 0.05)
## [1] 4163
sum(p4 < 0.05)
## [1] 4163
p5 <- p5[!is.na(p5)]
sum(p5 < 0.05)
## [1] 3752
length(p5)
## [1] 10000
------------------------------------------------------
図13 正規分布 vs 二峰分布のWelchのt分布
図14 正規分布 vs 二峰分布のWelchのt検定のp値の分布
図15 正規分布 vs 二峰分布のMann-WhitneyのU分布
図16 正規分布 vs 二峰分布のMann-WhitneyのU検定のp値の分布
図17 正規分布 vs 二峰分布のBrunner-MunzelのBM分布
図18 正規分布 vs 二峰分布のBrunner-Munzel検定のp値の分布
さて、今回の結果は非常に芳しくない。Welchのt検定、Mann-WhitneyのU検定、Brunner-Munzel検定のそれぞれの危険率は、81.6%、41.6%、37.5%であり、いずれの検定も大幅に5%を上回ってしまっている。ちなみに、生データではなく、順位に変更したのちにWelchのt検定を行っても、大した改善はなく、ほとんどMann-WhitneyのU検定と同様の傾向である。Welchのt検定は本当に壊滅的で、差がないにもかかわらず、80%以上の確率で差があると判断するということだ。このように全く分布の形が異なる標本同士での比較ではWelchのt検定は全く使えない。かといって、ノンパラメトリック検定のMann-WhitneyのU検定やBrunner-Munzel検定は、Welchのt検定よりはましであるが、決していい結果とはいえない。Brunner-Munzel検定はU検定よりもマシなようであるが、大幅に改善できているとは言えない。検出力の結果と併せて考えると、これらの検定は全く異なる分布で比較してしまうと、差があると判断してしまいやすい、ということである。ゆえに、全く異なる分布であるとわかった場合、これらのいずれの検定を使うことも推奨されない。
特に、危険率を確認する場合、p値の分布に関する累積分布関数(CDF)を確認するのが良い。標本に差がないとき、期待されるp値のCDFはy=xの直線となる。それに対するt検定とU検定のCDFも確認しよう。
------------------------------------------------------
cdf <- function(p, v){
sapply(v, function(x) sum(p < x))/length(p)
}
v <- seq(0, 1, length = 200)
p2c <- cdf(p2,v)
p3c <- cdf(p3,v)
p4c <- cdf(p4,v)
p5c <- cdf(p5,v)
cs <- data.frame(v,v,p2c,p3c,p4c,p5c)
plotn(cs[,1],cs[,2:6], mode = "m", xlab = "P-value", ylab = "CDF", type = "l",
col.line = c("black","red","blue","darkgreen","lightblue"),
pt.col.leg = c("black","red","blue","darkgreen","lightblue"),
legend = T, line = T, pt.cex.leg = 0,
leg.lab = c("expected", "Welch","ranked-Welch","Mann-Whitney","Brunner-Munzel"),
leg.sp = 8) #追加図1の描画
------------------------------------------------------
追加図1 期待されるp値のCDFとWelchのt検定、順位変換後Welchのt検定、U検定、BM検定のCDF
すると、すべての検定で。明らかに期待されるp値の分布より上側になっている。つまり、危険率が過剰に高まっていると言える。このなかでBrunner-Muzel検定が一番マシではある。
分布の形が同じであれば二峰分布であっても、Mann-WhitneyのU検定やBrunner-Munzel検定を使うことはできる。前述したとおり、標本の片方が二峰分布であるとしたら、標本が混ざってないかよく検討しよう。
・コーシー分布 vs コーシー分布
・コーシー分布 vs コーシー分布
幸い、2標本がともにコーシー分布に従う、つまり外れ値が出やすく、分散も違いやすい、という状況ならBrunner-Munzel検定は有効に使える。確認してみよう。検出力の時と同様にscale引数は変えてある。コードは以下の通り。
------------------------------------------------------
#コーシー分布
n <- 30
l1 <- 60
l2 <- 60
sc1 <- 1
sc2 <- 4
t2 <- NULL
tr <- NULL
u3 <- NULL
BM <- NULL
p2 <- NULL
p3 <- NULL
p4 <- NULL
p5 <- 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)
m1r <- d[d$sample == 1, "rank"]
m2r <- d[d$sample == 2, "rank"]
ue <- ((n - 1)*(sd(m1r)^2) + (n - 1)*(sd(m2r)^2))/(2*n - 2)
tr <- c(tr, (mean(m1r) - mean(m2r))/sqrt(ue*(1/n + 1/n))) #統計量tを算出
p3 <- c(p3, t.test(m1r,m2r,var.equal = F)$p.value) #Welchのt検定を行い、p値を取り出す
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の算出
p4 <- c(p4, 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の算出
p5 <- c(p5, 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
sum(p4 < 0.05)
## [1] 613
p5 <- p5[!is.na(p5)]
sum(p5 < 0.05)
## [1] 500
length(p5)
## [1] 10000
------------------------------------------------------
図19 コーシー分布同士のWelchのt分布
図20 コーシー分布同士のWelchのt検定のp値の分布
図21 コーシー分布同士のMann-WhitneyのU分布
図22 コーシー分布同士のMann-WhitneyのU検定のp値の分布
図23 コーシー分布同士のBrunner-MunzelのBM分布
図24 コーシー分布同士のBrunner-Munzel検定のp値の分布
Welchのt検定があてにならないのは今まで通りであるが、Mann-WhitneyのU検定とBrunner-Munzel検定はいい水準を保っている。この2つの検定で大きな差はないが、scaleが異なるコーシー分布同士であればBrunner-Munzel検定のほうがやや精度が良いことがわかる。
同様にCDFも確認しよう。
------------------------------------------------------
v <- seq(0, 1, length = 200)
p2c <- cdf(p2,v)
p3c <- cdf(p3,v)
p4c <- cdf(p4,v)
p5c <- cdf(p5,v)
cs <- data.frame(v,v,p2c,p3c,p4c,p5c)
plotn(cs[,1],cs[,2:6], mode = "m", xlab = "P-value", ylab = "CDF", type = "l",
col.line = c("black","red","blue","darkgreen","lightblue"),
pt.col.leg = c("black","red","blue","darkgreen","lightblue"),
legend = T, line = T, pt.cex.leg = 0,
leg.lab = c("expected", "Welch","ranked-Welch","Mann-Whitney","Brunner-Munzel"),
leg.sp = 8) #追加図2の描画
------------------------------------------------------
追加図2 期待されるp値のCDFとWelchのt検定、順位変換後Welchのt検定、U検定、BM検定のCDF
すると、Brunner-Munzel検定は、期待される直前上の乗るが、わずかにMann-WhitneyのU検定はその上側になる。つまり、危険率が高いことがわかる。
Rで行うBrunner-Muzel検定
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の描画
------------------------------------------------------
図25 2標本の描画。外れ値が出やすく、m2のほうが分散が大きい傾向がある。
分散が異なり、外れ値が存在しそうなので、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"))
------------------------------------------------------
図26 BMの近似分布。赤線は平均0、分散1の正規分布、青線は自由度dfのt分布
サンプルサイズが大きいときは正規近似でも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"))
------------------------------------------------------
図27 BMの近似分布。赤線は平均0、分散1の正規分布、青線は自由度dfのt分布
まとめ
まとめ
以上の検討から、上記の条件ではBrunner-Munzel検定はU検定よりも汎用性が高そうだ。とはいえ、どんな標本にでも活用できる、というほどではない。もし、そんな検定があるとしたら、とっくの昔に広まって、みんな使っているだろう。汎用性が高いBrunner-Munzel検定も、ちゃんと使いどころがあるのだ。それさえわきまえれば、非常に強力なツールとなる。
Brunner-Munzel検定の統計量の導出
Brunner-Munzel検定の統計量の導出
Page updated
Google Sites
Report abuse