備忘ログ

チラシの裏的備忘録&メモ

Rでウィルコクソンの順位和検定・符号順位和検定をwilcox.testとwilcox.exactでやったときのp値の差

Rでウィルコクソンの順位和検定・符号順位和検定を行うときには、Rにはじめからインストールされているstatsのwilcox.test()を使うことが多いと思いますが、

例えば

wilcox.test(sample(1:9,10,replace = T),sample(1:9,10,replace = T))

    Wilcoxon rank sum test with continuity correction

data:  sample(1:9, 10, replace = T) and sample(1:9, 10, replace = T)
W = 47.5, p-value = 0.878
alternative hypothesis: true location shift is not equal to 0

 警告メッセージ: 
 wilcox.test.default(sample(1:9, 10, replace = T), sample(1:9,: 
   タイがあるため、正確な p 値を計算することができません 

となり、警告メッセージがでます。

Rの解説書やサイトに寄っては無視してもいいよ、とのたまっていることがありますが、この関数はデータにタイが存在したときに正確な値を返せないため、exactRankTestsのwilcox.exact()を使った方が良いと言われることがあります。

結局のところケース・バイ・ケースでデータの内容によっては両者で大きな差が出ることがありますが、乱数で作ったデータの場合にはどの程度p値に差がでるのか検証してみました。

(ちなみにexactRankTestsを読み込むと古いからcoin使ってねと出るけど、wilcox.testと同じ使い方ができるのと、枯れた手法なので問題ないと色んな所で見かけたので今回はこちらで検証した。)

あくまで乱数なので、やるたび微妙に結果は異なると思いますことはあしからず。

library("exactRankTests")
wilcox.exact.normal <- function(n = 10, times = 100, paired = C(F,T)){
  p.wilcox.test <- NULL
  p.exact.wilcox.test <- NULL
  paired.wilcox <- paired
  for(i in 1:times){
    x <- round(rnorm(n))
    y <- round(rnorm(n))
    p.wilcox.test[i] <- wilcox.test(x, y, paired = paired.wilcox)$p.value
    p.exact.wilcox.test[i] <- wilcox.exact(x, y, paired = paired.wilcox)$p.value
  }
  data.frame(p.wilcox.test, p.exact.wilcox.test)
}



ウィルコクソンの順位和検定の場合

ウィルコクソンの順位和検定とウィルコクソンの符号順位和検定では同じwilcox.test()、wilcox.exact()を関数を使いますが、ウィルコクソンの順位和検定は対応のない検定であり、ウィルコクソンの符号順位和検定は対応のある検定となります。

まずはウィルコクソンの順位和検定で検証します。

dat <- wilcox.stats.exact(paired = F)
 50 件以上の警告がありました (最初の 50 個の警告を見るには warnings() を使って下さい) 

で警告が出ていますが、警告の内容はタイがあるので正確な値が計算ができませんでしたということなのでOK。

plot(dat)
abline(0,1)

f:id:indenkun:20200315225657p:plain
fig.1

0.1以下の範囲で見ると

plot(dat,ylim = c(0,0.1),xlim = c(0,0.1))
abline(0,1)

f:id:indenkun:20200315225724p:plain
fig.2

となります。

まぁまぁばらついている様に見えますが、どちらが小さい傾向にあるのかが微妙に分かりづらいので、ウィルコクソンの符号順位和検定で検証しついでに95%信頼区間も算出。

wilcox.test(dat$p.wilcox.test,dat$p.exact.wilcox.test,paired = T,conf.int = 0.95)

    Wilcoxon signed rank test with continuity
    correction

data:  dat$p.wilcox.test and dat$p.exact.wilcox.test
V = 953, p-value = 3.051e-07
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -0.02585789 -0.01184374
sample estimates:
(pseudo)median 
   -0.01858252 

と有意な差があり(当たり前)、wilcox.test()で算出されるp値の方が95%信頼区間で0.025~0.011くらい小さいとのこと。

個人的には設定した有意水準に微妙に届かないときには絶妙にほしい差くらいに感じた。

だからといって、一方は正確な値ではないとわかっている上に正しい値で計算する方法を知っているので(道義的に)使うわけにはいかない(涙……汗?)。


ウィルコクソンの符号順位和検定の場合

ウィルコクソンの符号順位和検定の場合は

dat <- wilcox.stats.exact(paired = T)
 50 件以上の警告がありました (最初の 50 個の警告を見るには warnings() を使って下さい) 

でプロットすると

plot(dat)
abline(0,1)

f:id:indenkun:20200315231804p:plain
fig.3

となりfig.1よりもばらついてみえます。

同様に0.1以下の範囲のみでプロットすると

plot(dat,ylim = c(0,0.1),xlim = c(0,0.1))
abline(0,1)

f:id:indenkun:20200315231940p:plain
fig.4

となります。

同様にウィルコクソンの符号順位和検定を行うと

wilcox.exact(dat$p.wilcox.test,dat$p.exact.wilcox.test,paired = T,conf.int = 0.95)

    Asymptotic Wilcoxon signed rank test

data:  dat$p.wilcox.test and dat$p.exact.wilcox.test
V = 204, p-value = 2.915e-13
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
 -0.0727962 -0.0457684
sample estimates:
(pseudo)median 
    -0.0585168 

となり、wilcox.test()で出力するp値の方が0.072~0.045小さいとなり、ウィルコクソンの順位和検定のときよりも差が大きいという結果になりました。


データを増やすとwilcox.test()とwilcox.exact()のp値は近づいていく

ウィルコクソンの順位和検定、ウィルコクソンの符号順位和検定の両者ともにデータの数が増えるとwilcox.test()とwilcox.exact()の出力するp値は近づいていくような感じになっていました。

ウィルコクソンの順位和検定の場合にデータを50個ずつにすると

dat <- wilcox.stats.exact(n = 50,paired = F)
plot(dat)
abline(0,1)

f:id:indenkun:20200315234707p:plain
fig.5

wilcox.exact(dat$p.wilcox.test,dat$p.exact.wilcox.test,paired = T,conf.int = 0.95)

    Asymptotic Wilcoxon signed rank test

data:  dat$p.wilcox.test and dat$p.exact.wilcox.test
V = 5050, p-value < 2.2e-16
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
 0.001908317 0.002330849
sample estimates:
(pseudo)median 
   0.002202397 

ウィルコクソンの符号順位和検定の場合にデータを50個にすると

dat <- wilcox.stats.exact(n = 50,paired = T)
 50 件以上の警告がありました (最初の 50 個の警告を見るには warnings() を使って下さい) 
plot(dat)
abline(0,1)

f:id:indenkun:20200315235018p:plain
fig.6

wilcox.exact(dat$p.wilcox.test,dat$p.exact.wilcox.test,paired = T,conf.int = 0.95)

    Asymptotic Wilcoxon signed rank test

data:  dat$p.wilcox.test and dat$p.exact.wilcox.test
V = 760, p-value = 1.29e-09
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
 -0.005156391 -0.002866855
sample estimates:
(pseudo)median 
  -0.003988123 

となりfig.5はfig.1、fig.6はfig.3と比較してp値が近づいている感じになり、それぞれのウィルコクソンの符号順位和検定の結果の95%信頼区間もデータ10個(fig.1とfig.3のとき)のときよりも狭くなっている。

しかし、wilcox.test()がせっかくタイでエラーを出していて正確ではありませんと言っているので、わかっていれば正確に計算してくれるwilcox.exact()で計算し直すのが妥当かと思います。