符号検定とWilcoxonの符号付き順位検定

符号検定

10人の患者にある睡眠薬を飲ませたところ，睡眠時間がそれぞれ次の時間だけ増えました (Arthur R. Cushny and A. Roy Peebles, The Journal of Physiology 32, 501-510 (1905)):

1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4

つまり，10個のうち負の値は1個だけで，残り9個は正です。正・負の符号の付き方は全部で $2^{10} = 1024$ 通りあり，そのうちで

すべて正になる場合は $_{10}C_0 = 1$ 通り
一つだけ負になる場合は $_{10}C_1 = 10$ 通り
二つだけ負になる場合は $_{10}C_2 = (10 \times 9) / (2 \times 1) = 45$ 通り

等々のように場合分けできます。すべての場合を合計すると，当然ですが $2^{10} = 1024$ になります。

もし正になる確率と負になる確率が同じなら，

すべて正になる確率は 1/1024
一つだけ負になる確率は 10/1024
二つだけ負になる確率は 45/1024

となるはずです。実際のデータは10個のうち一つだけ負ですので，このようになる確率と，もっと極端な（すべて正になる）確率を合計すれば，10/1024 + 1/1024 = 11/1024 です。逆に，すべて負になる確率と，一つだけ正になる確率を合計すれば，やはり 10/1024 + 1/1024 = 11/1024 です。そこで，10個のうち1個以下の符号が他と異なる確率は，22/1024 で，$p$ 値は約0.02です。つまり，偶然では50回に1回しか起きない事象です。

このような検定法を符号検定（sign test）といいます。

符号検定では，差が 0 のデータは外して考えます。

一般に，0でない数が n 個あって，そのうち m 個が正（または負）であるなら，binom.test(m, n) で符号検定できます。例えば上の例の場合は

> binom.test(1, 10)

	Exact binomial test

data:  1 and 10
number of successes = 1, number of trials = 10, p-value = 0.02148
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.002528579 0.445016117
sample estimates:
probability of success 
                   0.1

で，$p = 0.02148$ となります。

Wilcoxonの符号付き順位検定

さきほどのデータで，順序関係

\[ -0.1 < 0.1 < 0.8 < 1.1 < 1.6 < 1.9 < 3.4 < 4.4 < 4.6 < 5.5 \]

しか意味がないのであれば，符号検定以上のことはできません。

ここで，もし絶対値の順序関係

\[ 0.1 = 0.1 < 0.8 < 1.1 < 1.6 < 1.9 < 3.4 < 4.4 < 4.6 < 5.5 \]

に意味があるのであれば，順位（同順位があれば順位の平均をとります）

\[ 1.5, 1.5, 3, 4, 5, 6, 7, 8, 9, 10 \]

に変換してから，符号を戻して，合計した値

\[ -1.5 + 1.5 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 = 52 \]

を考えます。この符号の付け方は $2^{10} = 1024$ 通りありますが，順位和が 52 以上になるのは全部 $+$ にした場合と二つの $1.5$ のどちらかに $-$ を付けた場合の3通りです。同様に，順位和が $-52$ 以下になるのも3通りですので，偶然で順位和が $\pm 52$ またはそれより極端な値になる確率は $p = 6/1024 = 0.005859$ です。これがWilcoxon（ウィルコクソン）の符号付き順位検定（符号付き順位和検定，Wilcoxon signed-rank test）です。

一般にさきほどの単なる符号検定より小さい $p$ 値が得られますが，こちらのほうは，絶対値が比較可能でなければなりません。別の言い方をすると，分布の密度関数が正と負で対称であるという帰無仮説を検定することになります。

なお，すべて符号が同じ（またはゼロ）なら，符号検定と同じことになります。

Rによる計算

Rの標準の関数でWilcoxonの符号つき順位検定をするには次のようにします。

x = c(1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4)
> wilcox.test(x)

	Wilcoxon signed rank test with continuity correction

data:  x 
V = 53.5, p-value = 0.009298
alternative hypothesis: true location is not equal to 0 

Warning message:
In wilcox.test.default(x) :
   タイがあるため、正確な p 値を計算することができません

このように，残念ながら同じ値 1.5 が複数ありますので，近似計算になります。データの個数が十分多ければこれで問題ないのですが，この場合には，正確な値 $p = 0.005859$ とやや外れてしまいます。

厳密な答えを求める一つの方法は，exactRankTests パッケージの wilcox.exact() を使うものです：

> library(exactRankTests)
 Package ‘exactRankTests’ is no longer under development.
 Please consider using package ‘coin’ instead.

> wilcox.exact(x)

	Exact Wilcoxon signed rank test

data:  x 
V = 53.5, p-value = 0.005859
alternative hypothesis: true mu is not equal to 0

しかし，上のメッセージにあるように，exactRankTests パッケージは開発が止まっています（このメッセージは2006年に出るようになりましたが，その後も小さな修正はなされているようです）。代替としては coin（conditional inference）パッケージの wilcoxsign_test() があります。使い方は，次のように，この場合なら0を10個並べただけのデータ rep(0,10) と比較します。比較は ~ という演算子で行います：

> library(coin)
> wilcoxsign_test(x ~ rep(0,10), distribution="exact")

        Exact Wilcoxon-Signed-Rank Test

data:  y by
         x (neg, pos)
         stratified by block
Z = 2.6519, p-value = 0.005859
alternative hypothesis: true mu is not equal to 0

要素ごとに対応するベクトル x と y を比較するなら，wilcoxsign_test(x ~ y, distribution="exact") とします。このとき，0の処理法のデフォルトは "Pratt" です。wilcox.exact() と同じにしたいなら，wilcoxsign_test(..., zero.method="Wilcoxon") とします。

奥村晴彦

Last modified: 2012-08-14 21:21:36