統計と情報の専⾨誌「エストレーラ」 ESTRELA 参考

2014年4月号の連載第129回「bootstrap ― car、fpc ―」において紹介したスクリプトrandmz、randmzpを用意しました。ダウンロードして、記事と併せてご活用ください。

ブートストラップ2標本検定

基本的には X, Y を用意し、randmzp(X,Y)とすれば、2000回のブートストラップ統計量を元にしてのp値を求めることを100回行い、そのmedianとして最終的なp値推定を行う。グラフも描画するブートストラップそのものはrandmz( )が行う。randomzp( )は それを呼び出す。
compile すれば2割程度早くなる（実行時間が 0.8倍になる）。
（例えば library(compiler); randmzpc <- cmpfun(randmzp) とする）。
compile が単純なように randmzp の内部に randmz を入れ込んだが、 本来はrandmz は独立して定義していた。
対応ありの場合はパラメータ paired=TRUE を加える。
p値算出の例のグラフを左に、p値の分布とmedian確認のグラフを右に描く。
median確認のグラフは、p値分布の箱ひげ図と、平均値の位置に赤丸を付す。

[実例]

http://aoki2.si.gunma-u.ac.jp/R/u-test.html の実例を引用する。

X <- rep(1:4, c(9, 12, 6, 3))
Y <- rep(1:4, c(4, 9, 11, 5))
randmzp(X,Y) ## median(p-values)=0.054
t.test(X,Y) ## p-value = 0.05486

分布が特異であるとはいえ、標準偏差は同程度であり、t検定と本件とはほぼ同じp値となっている。p > 0.05 により、帰無仮説は棄却できない。
一方、ノンパラメトリックでは p < 0.05で棄却となる。

wilcox.text(X,Y, correct=FALSE) ## p=0.04883 NOT EXACT
wilcox.text(X,Y, correct=FALSE, exact=TRUE) ## p=0.04883 NOT EXACT
exactRankTests::wilcox.exact(X,Y) ## p=0.04965 EXACT

このX,Yは整数から構成されており、分付を気にすると通常ではノンパラメトリック検定を行いたくなる。p値は 0.05 の前後で微妙に動く。