これも，非公開にされたのね。リンク切れになちゃった。
確かに，不都合な事実は消去したいというのはもっともだとは思うのだけど。不十分でしたと潔く振る舞うのもブシノスガタではないかな?
http://my-notes.hatenablog.com/entry/2017/09/15/173257
| My Notes
| 統計学とかR(R言語)とかPython3とかプログラミングとかの覚え書きとか走り書きとか。 座右の銘にしたい: All work and no play makes Jack a dull boy.
| 2017-09-15
| R(R言語) ノンパラメトリック検定(独立サンプルの比較、独立した2群の中心位置の比較、Mann-Whitney (マン・ウイットニー) 検定(U検定))、wilcox.test()
で，Mann-Whitney (マン・ウイットニー) 検定(U検定) -- wilcox.test において，同順位を含むデータで
> 患者_fac <- factor(c(rep("男性患者", 20), rep("女性患者", 16)))
> 評価_vec <- c(1, 4, 2, 3, 5, 3, 2, 3, 2, 5, 6, 7, 4, 2, 5, 3, 2, 2, 2, 3,
+ 5, 4, 6, 4, 3, 7, 5, 6, 1, 2, 5, 6, 7, 5, 3, 5)
> df <- data.frame(患者 = 患者_fac, 評価 = 評価_vec)
> res <- wilcox.test(df$評価 ~ df$患者, data = df)
警告メッセージ:
wilcox.test.default(x = c(5, 4, 6, 4, 3, 7, 5, 6, 1, 2, 5, 6, で:
タイがあるため、正確な p 値を計算することができません
> res
Wilcoxon rank sum test with continuity correction
data: df$評価 by df$患者
W = 231, p-value = 0.02253
alternative hypothesis: true location shift is not equal to 0
> res$p.value
[1] 0.0225325
「タイがあるため、正確な p 値を計算することができません」に対して，「これは無視していい」とある。
wilcox.test を紹介するページの多くで「無視してよい」とあるのが無責任。せっかくノンパラメトリック検定を採用するなら，最後までちゃんとしてほしいものだ。
正確な p 値を計算する手段はある。
ひとつは exactRankTests パッケージの wilcox.exact 関数を使うやり方であった。
> library(exactRankTests)
Package ‘exactRankTests’ is no longer under development.
Please consider using package ‘coin’ instead.
> wilcox.exact(評価 ~ 患者, data = df)
Exact Wilcoxon rank sum test
data: 評価 by 患者
W = 231, p-value = 0.02071
alternative hypothesis: true mu is not equal to 0
もうひとつは，library(exactRankTests) したときのメッセージ
Package ‘exactRankTests’ is no longer under development.
Please consider using package ‘coin’ instead.
にしたがうことである。
> library(coin)
要求されたパッケージ survival をロード中です
> wilcox_test(評価 ~ 患者, data=df, distribution="exact")
Exact Wilcoxon-Mann-Whitney Test
data: 評価 by 患者 (女性患者, 男性患者)
Z = 2.2974, p-value = 0.02071
alternative hypothesis: true mu is not equal to 0
wilcox.exact も wilcox_test も同じ p 値を示す。
なお，wilcox_test で表示される Z 値は p 値の計算に使われるものではないことに注意。
> pnorm(2.2974, lower.tail=FALSE)*2
[1] 0.02159596
W = 231 という統計量 W は紛らわしいが，マン・ホイットニー検定の統計量 U の大きい方である。
> n1 = 20
> n2 = 16
> (r = rank(df$評価))
[1] 1.5 19.5 6.5 14.0 25.5 14.0 6.5 14.0 6.5 25.5 31.5 35.0 19.5 6.5 25.5 14.0 6.5 6.5
[19] 6.5 14.0 25.5 19.5 31.5 19.5 14.0 35.0 25.5 31.5 1.5 6.5 25.5 31.5 35.0 25.5 14.0 25.5
> (W1 = sum(r[1:n1])-n1*(n1+1)/2)
[1] 89
> (W2 = sum(r[n1+1:n2])-n2*(n2+1)/2)
[1] 231
> (W = max(W1, W2))
[1] 231
普通は，マン・ホイットニーの U 検定の検定統計量は W1, W2 の小さい方とされる（大きい方でもかまわないのであるが），W1+W2 = n1*n2 なので，U = n1*n2 - W で求めることができる。
> (U = n1*n2 - W)
[1] 89
示される Z を計算するためには W の期待値 E(W) と分散 V(W) が必要である。
期待値は E(U) = n1*n2/2
> (E.W = n1*n2/2) # or (W1+W2)/2
[1] 160
同順位がある場合の分散は n = n1+n2，同順位の大きさを ti として，V(W) = n1*n2/12/(n^2-n) * sum(n^3-n-Σ(ti^3-ti)
> table(r) # 同順位の個数
r
1.5 6.5 14 19.5 25.5 31.5 35
2 8 7 4 8 4 3
> sum(t^3-t)
[1] 1494
> n = n1+n2
> (V.W = n1*n2/12/(n^2-n)*(n^3-n-sum(t^3-t)))
[1] 955.0476
Z 値は Z= (W-E(W)) / sqrt(V(W))
> (Z = (W-E.W) / sqrt(V.W))
[1] 2.297449
しかし，wilcox_test の両側 p 値はこのZ 値から計算されているのではない。
なぜなら，
> pnorm(Z, lower.tail=FALSE)*2
[1] 0.02159318
であるから。
pnorm(Z, lower.tail=FALSE)*2 で計算される値は，連続性の補正をしないマン・ホイットニー検定の結果で示される p 値と同じものである。wilcox.test 関数を correct=FALSE で使う。
> wilcox.test(評価 ~ 患者, data = df, correct=FALSE)
Wilcoxon rank sum test
data: 評価 by 患者
W = 231, p-value = 0.02159
alternative hypothesis: true location shift is not equal to 0
警告メッセージ:
wilcox.test.default(x = c(5, 4, 6, 4, 3, 7, 5, 6, 1, 2, 5, 6, で:
タイがあるため、正確な p 値を計算することができません
これらとは対極的に，wilcox_test では，正確な p 値が表示される。
> (a = wilcox_test(評価 ~ 患者, data=df, distribution="exact"))
Exact Wilcoxon-Mann-Whitney Test
data: 評価 by 患者 (女性患者, 男性患者)
Z = 2.2974, p-value = 0.02071
alternative hypothesis: true mu is not equal to 0
p.value だけを取り出すためには，以下のようにする。
> pvalue(a)
[1] 0.02071037
その計算方法は，フィッシャーの正確検定 fisher.test と同じやり方で，途中の検定統計量をカイ二乗統計量ではなく，マン・ホイットニー統計量（U または W）を使う。
関数の実装例は，以下にもある。
マン・ホイットニーの U 検定(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-mw.html
> exact.mw(df$評価[1:20], df$評価[21:36])
U = 89, Z = 2.29745, P 値 = 0.0215932
正確な P 値 = 0.02071036886
査察した分割表の個数は 14094 個
結論
詳細は以上のようであるが，どんな場合でも，coin ライブラリの wilcox_test を使おうということだ。
wilcox.test は使わない。
