背景
ノンパラメトリックな多重比較であるクラスカル・ウォリス検定について調べていて、気になる記事を見つけました。
Kruskal-Wallis 検定を使えば U 検定は不要：漸近と正確検定
この記事の中で
つまり，通常使われる Kruskal-Wallis 検定は，2 群の場合，漸近 U 検定であり，漸近 Steel-Dwass 検定なのである。
とあります。
この近似は中心極限定理によるものですが（中心極限定理の概要についてはヨビノリさんの動画が分かりやすいです）、この近似ができるのはサンプルサイズがある程度大きいときです（一例としてn > 20：統計ソフトで，小標本でのノンパラメトリックな検定をしたときには注意が必要より）。生物実験でサンプルサイズが20以上なのは稀で3や5の場合も多いです。
そのため、近似ではなく正確なp値を計算する必要がありますが、Rのstatsパッケージにおいてクラスカル・ウォリス検定を行うkruskal.test関数は「method = exact」と指定しない限りは近似値になるようです（最初の記事参照）。そこで、他のノンパラメトリック検定の関数で求められるp値がデフォルトで近似か正確かを確認してみました。
注意
統計に関する知識が十分でないため、誤った解釈をしている可能性があります。
方法
下記のノンパラメトリック検定関数をRで実行しp値を比較する。
※ドキュメントに記載がある場合もあるのですが、「method = NULL」などデフォルトでどちらが計算されているか分からないものも一部あったため。
- kruskal.test (stats)：クラスカル・ウォリス検定
- kruskal_test (coin)：クラスカル・ウォリス検定
- pSDCFlig (NSM3)：スティール・ドゥワス検定
- wilcox.test (stats)：ウィルコクソンの順位和検定（マンホイットニーのU検定）
- wilcox.exact (exactRankTests)：ウィルコクソンの順位和検定
- wilcox_test (coin)：ウィルコクソンの順位和検定
# 環境
sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS 10.14.1
結果
まず最初の記事で紹介されているkSamplesパッケージのqn.test関数でクラスカル・ウォリス検定の近似と正確なp値を求めます。
group1 <- c(1, 2, 3)
group2 <- c(100, 200, 300)
library(kSamples)
qn.test(group1, group2, test="KW", method = "exact")
Kruskal-Wallis k-sample test.
Number of samples: 2
Sample sizes: 3, 3
Number of ties: 0
Null Hypothesis: All samples come from a common population.
test statistic asympt. P-value exact P-Value
3.85700 0.04953 0.10000
Warning: At least one sample size is less than 5,
asymptotic p-values may not be very accurate.
近似のp値は0.04953、正確なp値は0.10000になりました。けっこう差あるな。。。
kruskal.test (stats)
kruskal.test(list(group1, group2))
Kruskal-Wallis rank sum test
data: list(group1, group2)
Kruskal-Wallis chi-squared = 3.8571, df = 1, p-value = 0.04953
statsパッケージのクラスカル・ウォリス検定はデフォルトで近似のp値を出力するようです。
kruskal_test (coin)
library(coin)
valueVec <- c(group1, group2)
groupVec <- factor(c(rep("A", 3), rep("B", 3))) #A群・B群という名前をつけました
kruskal_test(valueVec ~ groupVec)
Asymptotic Kruskal-Wallis Test
data: valueVec by groupVec (A, B)
chi-squared = 3.8571, df = 1, p-value = 0.04953
coinパッケージのクラスカル・ウォリス検定はデフォルトで近似のp値を出力するようです。
pSDCFlig (NSM3)
library(NSM3)
pSDCFlig(valueVec, groupVec)
Group sizes: 3 3
Using the Exact method:
For treatments A - B, the Dwass, Steel, Critchlow-Fligner W Statistic is 2.7775.
The smallest experimentwise error rate leading to rejection is 0.1 .
NSM3パッケージのスティール・ドゥワス検定はデフォルトで正確なp値を出力するようです。なお、本関数は「method = 」で近似か正確かを選択できます。
# 正確
pSDCFlig(valueVec, groupVec, method = "Exact")
Group sizes: 3 3
Using the Exact method:
For treatments A - B, the Dwass, Steel, Critchlow-Fligner W Statistic is 2.7775.
The smallest experimentwise error rate leading to rejection is 0.1 .
# 近似
pSDCFlig(valueVec, groupVec, method = "Asymptotic")
Group sizes: 3 3
Using the Asymptotic method:
For treatments A - B, the Dwass, Steel, Critchlow-Fligner W Statistic is 2.7775.
The smallest experimentwise error rate leading to rejection is 0.0496 .
# 0.049531は四捨五入しても0.0496にはなりませんが、切り上げ処理が行われているのでしょうか？
wilcox.test (stats)
wilcox.test(group1, group2)
Wilcoxon rank sum test
data: group1 and group2
W = 0, p-value = 0.1
alternative hypothesis: true location shift is not equal to 0
statsパッケージのウィルコクソンの順位和検定はデフォルトで正確なp値を出力するようです。
なお、本関数で近似のp値を計算したい場合は下記の通りです（最初の記事参照）。
wilcox.test(group1, group2, exact = F, correct = F)
Wilcoxon rank sum test
data: group1 and group2
W = 0, p-value = 0.04953
alternative hypothesis: true location shift is not equal to 0
wilcox.exact (exactRankTests)
library(exactRankTests)
wilcox.exact(group1, group2)
Exact Wilcoxon rank sum test
data: group1 and group2
W = 0, p-value = 0.1
alternative hypothesis: true mu is not equal to 0
exactRankTestsパッケージのウィルコクソンの順位和検定はデフォルトで正確なp値を出力するようです（関数名もwilcox.exactですしね）。
wilcox_test (coin)
df <- data.frame(value = valueVec, group = groupVec)
wilcox_test(value ~ group, data = df)
Asymptotic Wilcoxon-Mann-Whitney Test
data: value by group (A, B)
Z = -1.964, p-value = 0.04953
alternative hypothesis: true mu is not equal to 0
coinパッケージのウィルコクソンの順位和検定はデフォルトで近似のp値を出力するようです。
なお、本関数で正確なp値を算出するには「distribution = "exact"」を指定すれば良いようです（coinパッケージの公式ドキュメントより）。
wilcox_test(value ~ group, data = df, distribution = "exact")
Exact Wilcoxon-Mann-Whitney Test
data: value by group (A, B)
Z = -1.964, p-value = 0.1
alternative hypothesis: true mu is not equal to 0
まとめると
- kruskal.test (stats)：クラスカル・ウォリス検定 → 近似
- kruskal_test (coin)：クラスカル・ウォリス検定 → 近似
- pSDCFlig (NSM3)：スティール・ドゥワス検定 → 正確
- wilcox.test (stats)：ウィルコクソンの順位和検定 → 正確
- wilcox.exact (exactRankTests)：ウィルコクソンの順位和検定 → 正確
- wilcox_test (coin)：ウィルコクソンの順位和検定） → 近似
結論
クラスカル・ウォリス検定とcoinパッケージのウィルコクソンの順位和検定（マンホイットニーのU検定）はデフォルトで近似のp値を計算するので、サンプルサイズの小さい実験で使用するときは注意した方が良さそうです。
補足
近似と正確でどの程度計算速度が変わるのが検証してみました。
# 正確
t1 <- Sys.time()
for (i in 1:1000) {
wilcox.test(group1, group2)
}
t2 <- Sys.time()
t2 - t1
Time difference of 0.09528399 secs
# 近似
t1 <- Sys.time()
for (i in 1:1000) {
wilcox.test(group1, group2, exact = F, correct = F)
}
t2 <- Sys.time()
t2 - t1
Time difference of 0.240438 secs
なぜか近似の方が時間かかったのですが、現在のPCならどちらも時間的な問題はなさそうです。そのため、サンプルサイズの小さい生物実験では基本的に正確なp値を求めれば良さそうです。