「統計学関連なんでもあり」の過去ログ--- 048

No.22724　Brunner-Munzel検定　　【青木繁伸】　2019/05/07(Tue) 18:23

https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
でも紹介されているが，R の
lawstat パッケージの brunner.munzel.test() も
brunnermunzel パッケージの brunnermunzel.test() も
alternative="greater", "less" の場合の信頼区間がおかしい（"two.sided" の場合とまったく同じになっている）
nparcomp パッケージの npar.t.test で method="t.app" を指定して得られるのが正しい信頼区間のようだ。

No.22725　Re: Brunner-Munzel検定　　【青木繁伸】　2019/05/09(Thu) 11:38

> nparcomp パッケージの npar.t.test で method="t.app" を指定して得られるのが正しい信頼区間のようだ。

と書いたが，そもそも npar.t.test は片側検定で "greater" と "less" を取り違えているようだ。p 値が t.test と比べて逆になっている。これは困ったことになる。
df = data.frame(
  d = c(1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1, 3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4),
  g = rep(1:2, c(14, 11)))
library(nparcomp)
cl = 0.95
alt = "greater"
a = npar.t.test(d ~ g, data=df, round=5, info=FALSE, method="t.app", alternative=alt, conf.level=cl)
a$Analysis
b = t.test(d ~ g, data=df, alternative=alt, conf.level=cl)
b
を実行すると，npar.test の結果は
  Effect Estimator  Lower Upper       T p.Value
1 p(1,2)   0.78896 0.6291     1 3.13747 0.00289
だが，t.test の結果は
t = -2.9486, df = 15.916, p-value = 0.9953
p 値が 0.00289 と 0.9953 と 180 度異なっている。
alt="less" では 0.99711 と 0.004739 である。
ちなみに，wilcox.test の p 値は t.test と同一方向だ。

No.22727　Re: Brunner-Munzel検定　　【青木繁伸】　2019/05/10(Fri) 18:00

引き続き

並べ替え検定については，奥村先生が
https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
でも独自のプログラム例を示しているが，
brunnermunzel パッケージには brunnermunzel.permutation.test がある。
しかし，これは場合によって正しい答えを出さない。というか，奥村先生のプログラムと結果が異なることがある。
そのようなデータの例は，x = c(36, 45, 51, 49, 57)，y = c(48, 42, 49, 39, 45) の場合。

奥村先生のプログラムでは 0.365079365079365 になるのだが，brunnermunzel.permutation.test では 0.341269841269841 になる。

前者は 92/252，後者は 86/252 だ。92 の内訳は，データを並べ替えたデータにより計算される検定統計量の絶対値が，観察値における検定統計量の絶対値と同じかより極端な場合の数 45*2 と，検定統計量が +Inf, -Inf になる 2 通りの計 92 である。

観察値における検定統計量の絶対値と同じになるのは 4*2 = 8 通りあるが brunnermunzel.permutation.test では，計算誤差のためか 8 個全てが同じとなっていないのであろう。つまり，数え落としている。

また，Inf, -Inf となる場合の扱いもチェックする必要があるだろう。
前に書いて消したが，検定統計量が Inf になるのは R の場合だけで，Python の場合には 0 による割り算エラーで計算が中止される。

いろいろ問題ありで，それぞれのプログラムどれも信頼して使うわけにはいかな状態である。

自分で書くしかないか。

No.22728　Re: Brunner-Munzel検定　　【青木繁伸】　2019/05/10(Fri) 21:15

Brunner-Munzel検定最強！！

なんて，書かれているページが多くて，ですね。ちょっと，立ち止まって検証してみませんか？という趣旨で書きました。

ほかの人の意見も聞きたいので，「拡散希望」（^_^;）

No.22729　Re: Brunner-Munzel検定　　【it may be that】　2019/05/19(Sun) 19:00

npar.t.test の greater/less の扱いですが，取り違えているわけではないと思いますよ。
summary(a)
に
#---------------------------Interpretation---------------------------------------------#
 p(a,b) > 1/2 : b tends to be larger than a 
と出ます。
したがって，群1と2 (g1とg2) の観測量をX1, X2として，対立仮説は
'greater'はP(X1 < X2) > 1/2, 
'less'はP(X1 < X2) < 1/2
です。
Estimate が 0.78896 > 1/2 であることから，P(1,2)はP(X2 - X1 > 0)であると思われます。
よって，書き直すと
'greater'はP(X2 - X1 > 0) > 1/2,
'less'はP(X2 - X1 > 0) < 1/2 あるいは P(X2 - X1 < 0) > 1/2
となります。

一方，g1とg2の位置母数をそれぞれm1とm2とすると，t.testの対立仮説は
'greater' :m1 - m2 > 0
'less' : m1 - m2 < 0
です。
つまり，npar.t.testの対立仮説は
m2 - m1 > 0
m2 - m1 < 0
に相当し，基準が逆になっています。
確かに紛らわしいですが，lm, glmなどでcontr.treatmentに馴染んでしまった身としては，むしろt.testやwilcox.testのfactorの扱いに違和感を覚えます。
普段contr.SASを使っていればt.test,wilcox.testが自然でしょう。

No.22730　Re: Brunner-Munzel検定　　【青木繁伸】　2019/05/21(Tue) 11:01

「npar.t.test の greater/less の扱い」についてのコメント，ありがとうございます。
確かに，説明の全体を見ればその通りですね。

No.22731　Re: Brunner-Munzel検定　　【青木繁伸】　2019/05/21(Tue) 11:04

追加の考察

http://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html

では，p = Pr(x＜y)+0.5*Pr(x=y) ではなく，t 分布による漸近近似を使っている。

n1, n2 が小さいときに行う並べ替え検定で，n1, n2 が大きいことを仮定する漸近近似を行うという，相矛盾する状況である。

p そのままを使って並べ替え検定を行うと，
foo2 = function(X) {
  brunner.munzel.test(xandy[X], xandy[-X])$estimate
}
z2 = combn(1:N, n1, foo2)
bm2 = brunner.munzel.test(x, y)$estimate
mean(abs(z2-0.5) >= abs(bm2-0.5))
0.006690223 が得られる。

なお，この値の正当性としては，マン・ホイットニーの U 検定統計量と Brunner-Munzel の p との間に p=1-U/n1/n2 という単純な関係があるので，マン・ホイットニーの U 検定統計量を計算して使う並べ替え検定と同じ値になることを示せば十分。

または，coin パッケージの wilcox_test で正確な p 値を求めると，同じ値 0.006690223 が得られる。結局は wilcox_test が並べ替え検定をやっているということだ。
> df = data.frame(d=c(x, y), g=factor(rep(1:2, c(length(x), length(y)))))
> wilcox_test(d ~ g, data=df, distribution="exact")

    Exact Wilcoxon-Mann-Whitney Test

data:  d by g (1, 2)
Z = -2.6934, p-value = 0.00669
alternative hypothesis: true mu is not equal to 0
ということは，マン・ホイットニーの U 検定と Brunner-Munzel 検定って，前者が正規分布，後者が t 分布（など）で漸近近似を行うだけで，本質は同じということだ...
> npar.t.test(d ~ g, data=df, info=FALSE, method="perm", nperm=400000, round=5)$Analysis
       Estimator Statistic   Lower   Upper p.value
id       0.78896   3.13747 0.59708 0.99057 0.00708
logit    0.78896   2.38394 0.57650 0.91561 0.00694
probit   0.78896   2.51949 0.58109 0.92446 0.00694
となれば，奥村先生のページの冒頭の説明，

> 二つの確率変数 X1，X2 が同じ分布に従うという帰無仮説を検定するには，有名なWMW検定（Wilcoxon-Mann-Whitney test）が使えます。しかし，分布が異なる（例えば分散が異なる）場合には，これでは正確に検定できません。

> そこで，分布が同じことは仮定せず，両群から一つずつ値を取り出したとき，どちらが大きい確率も等しいという帰無仮説を検定することにします。

という，Brunner-Munzel 検定の優位性というのもなくなるのか。（前提が異なるといえども，検定統計量が同じなんだからねぇ）

No.22732　Re: Brunner-Munzel検定　　【青木繁伸】　2019/05/24(Fri) 18:57

検定統計量ではなく，その近似（t分布に近似する）のときに「ベーレンス・フィッシャー問題を考慮した」ということですね。