ノンパラメトリック検定

必要なパッケージ

require(BSDA)
require(PMCMR)
require(coin)

ノンパラメトリック検定

ここまでに紹介してきた検定方法のほとんどは、パラメトリックな仮説に基づいていた。データは、既知のよく知られた分布（正規分布、ポアソン分布など）に従っており、それら分布は、1つまたはそれ以上のパラメータで規定され（例えば、正規分布では平均と分散）、それらうち少なくとも1つは未知で推定する必要があった。しかし、より複雑な実験や標本抽出を行う場合には、得られるデータが既知の分布に従わない場合もある。

データが従う分布が明確でない場合には、データのもつ分布がどのようなものであろうと適用できる検定手法が必要となる。こうした手法は、分布に依存しない（distribution free）検定、あるいは、ノンパラメトリック（nonparametric）検定とよばれる。

　\(t\)分布、ピアソンの相関係数、分散分析などの基本的な統計手法の多くは、標本が正規分布から得られていると仮定されているか、標本サイズが十分に大きく、中心極限定理により正規性を仮定できるとされている（パラメトリック検定）。

　次は、よく用いられるパラメトリック検定とノンパラメトリック検定の対応表である。

パラメトリック検定	ノンパラメトリック検定
1標本の\(t\)検定（平均の検定）	符号検定、ウィルコクソンの符号順位検定
対応のある\(t\)検定	符号検定、ウィルコクソンの符号順位検定
2標本の\(t\)検定	ウィルコクソンの順位和検定、ウィルコクソン-マン-ホイットニー検定
1要因の分散分析	クラスカル-ウォリス検定
ブロックのある分散分析	フリードマン検定

符号検定

　ノンパラメトリック検定のうち、もっとも単純なのが符号検定である。まず、連続的な累積分布関数をもつ母集団が中央点として\(m_0\)をもつという仮説\(H_0\)を、対立仮説\(H_1: med > \ne < m_0\)に対して検定することを考える。

　ここで、\(X_i > m_0\)のとき、すなわち、\(X_i - m_0\)が正のとき、その標本に+の符号を割りふり、そうでなければ-の符号を割り振る。\(X_i = m_0\)（タイ）は連続分布では理論的にはありえないが、限られた精度で観察値が計測される場合には、タイが生じる可能性がある。なお、次の説明では、タイは起こらないと仮定している。

　もし、\(m_0\)が中央値で、\(H_0\)が真であれば、中央値の定義より、\(P(X_i > m_0) = P(X_i < m_0) = \frac{1}{2}\)が成り立つ。したがって、+の数を表す検定統計量\(T\) \[ T = \sum_{i=1}^n \boldsymbol{1}(X_i > m_0) \] は、パラメータが\(n, 1/2\)の二項分布に従う。

　ここで、検定の有意水準を\(\alpha\)とする。対立仮説が\(H_1: med > m_0\)のとき、\(T\)の棄却限界値は、\(k_\alpha\)以上の整数となる。ここで、\(k_\alpha\)は、以下の関係を満たす最小の整数である。 \[ \sum_{t=k_\alpha}^n \left( \begin{matrix} n \\ t \end{matrix}\right) 2^{-n} < \alpha \]

　同様に、対立仮説が\(H_1: med < m_0\)のとき、\(T\)の棄却限界値は、\(k'_\alpha\)以下の整数となる。ここで、\(k'_\alpha\)は、以下の関係を満たす最大の整数である。 \[ \sum_{t=0}^{k'_\alpha} \biggl( \begin{matrix} n \\ t \end{matrix}\biggr) 2^{-n} < \alpha \]

　もし、対立仮説が両側、すなわち、\(H_1: med \ne m_0\)であった場合、\(T\)の棄却限界値は、、\(k'_\alpha\)以下で、\(k_\alpha\)以上の整数となる。ここで、\(k'_\alpha\)、\(k_\alpha\)は、それぞれ、以下の不等式を満たす最大、最小の整数である。 \[ \sum_{t=0}^{k'_\alpha} \biggl( \begin{matrix} n \\ t \end{matrix}\biggr) 2^{-n} < \alpha/2 \\ \sum_{t=k_\alpha}^n \biggl( \begin{matrix} n \\ t \end{matrix}\biggr) 2^{-n} < \alpha/2 \]

　\(T\)の値が観察されているとき、対立仮説\(H_1: med > m_0\)に対する検定は、\(T\)が大きくなる向きに棄却限界値をとることとなり、\(p\)値は、 \[ p = \sum_{t=T}^n \biggl( \begin{matrix} n \\ t \end{matrix}\biggr) 2^{-n} = \sum_{t=0}^{n-T} \biggl( \begin{matrix} n \\ t \end{matrix}\biggr) 2^{-n} \] となる。

また、対立仮説\(H_1: med < m_0\)に対する検定は、\(T\)が小さくなる向きに棄却限界値をとることとなり、\(p\)値は、 \[ p = \sum_{t=0}^T \biggl( \begin{matrix} n \\ t \end{matrix}\biggr) 2^{-n} \] となる。

さらに、対立仮説が両側である場合は、\(p\)値は、 \[ p = 2 \sum_{t=0}^{T'} \biggl( \begin{matrix} n \\ t \end{matrix}\biggr) 2^{-n} \] となる。ここで、\(T'=\min\{T, n-T\}\)である。

　次に、対応のある2標本について考える。標本\(X_1,...,X_n\)と\(Y_1,...,Y_n\)が観察されており、対応間の差\(X_1-Y_1,X_2-Y_2,...,X_n-Y_n\)の母集団の中央値が\(0\)に等しいかどうかを検定することを考える。この場合、 \[ T = \sum_{i=1}^n \boldsymbol{1}(X_i > Y_i) \] は、差(\(X_i - Y_i\))が正である標本のペア数となる。

　上述した説明はタイがない場合であるが、タイが存在する場合には、以下のような方法が提案されている。

1. 無視する。もし\(s\)個のタイがあった場合は、タイでない観察値だけを用いる。標本数は\(n-s\)に減少する。
2. タイに多い方の符号をつける。例えば、-が2、タイが2、+が6の場合は、2つのタイを+にする。これにより、+と-に最も極端な偏りを生じさせる。
3. 無作為化する。もし2つのタイがあれば、2回コインを投げて、コインが表を向けば+を割り振り、裏を向けば-を割り振る。

例 枯葉剤による影響（Vidakovic 2011より引用）

ベトナム戦争の退役軍人の多くは、枯葉剤Agent Orangeへの暴露の結果、ダイオキシン2,3,7,8-TCDDの血中・脂肪組織中濃度が危険な水準にあるといわれている。Chemosphere (vol. 20, 1990)で発表された研究では、枯葉剤への暴露があった20人のマサチューセッツ州のベトナム戦争退役軍人のTCDDレベルについて報告されている。血漿と脂肪組織内のTCDD（一兆分率、PPT）の濃度は以下のとおりである。

tcddpla <- c(2.5, 3.1, 2.1, 3.5, 3.1, 1.8, 6.8, 3.0, 36.0, 4.7, 6.9, 3.3, 4.6, 1.6, 7.2, 1.8, 20.0, 2.0, 2.5, 4.1)
tcddfat <- c(4.9, 5.9, 4.4, 3.5, 7.0, 4.2, 10.0, 5.5, 41.0, 4.4, 7.0, 2.9, 4.6, 1.4, 7.7, 1.8, 11.0, 2.5, 2.3, 2.5)

血漿中と脂肪組織中のTCDDの濃度の分布に違いはみられるか、符号検定を用いて両側検定を行う。

まずは、タイを確認してみる。

tcddpla - tcddfat

##  [1] -2.4 -2.8 -2.3  0.0 -3.9 -2.4 -3.2 -2.5 -5.0  0.3 -0.1  0.4  0.0  0.2 -0.5
## [16]  0.0  9.0 -0.5  0.2  1.6

3つのタイがある。

BSDAパッケージのSIGN.test関数を用いて検定を行う。検定は両側で行う。

SIGN.test(x = tcddpla, y = tcddfat, alternative = "two.sided")

## 
##  Dependent-samples Sign-Test
## 
## data:  tcddpla and tcddfat
## S = 6, p-value = 0.3323
## alternative hypothesis: true median difference is not equal to 0
## 95 percent confidence interval:
##  -2.4000000  0.1767059
## sample estimates:
## median of x-y 
##          -0.3 
## 
## Achieved and Interpolated Confidence Intervals: 
## 
##                   Conf.Level L.E.pt U.E.pt
## Lower Achieved CI     0.8847   -2.4 0.0000
## Interpolated CI       0.9500   -2.4 0.1767
## Upper Achieved CI     0.9586   -2.4 0.2000

違いは5%水準でも有意ではない。

なお、タイの扱い方には3通りあった。これらを順に試してみる。

まずは、無視する。

tie <- tcddpla - tcddfat == 0
SIGN.test(x = tcddpla[!tie], y = tcddfat[!tie], alternative = "two.sided")

## 
##  Dependent-samples Sign-Test
## 
## data:  tcddpla[!tie] and tcddfat[!tie]
## S = 6, p-value = 0.3323
## alternative hypothesis: true median difference is not equal to 0
## 95 percent confidence interval:
##  -2.498985  0.200000
## sample estimates:
## median of x-y 
##          -0.5 
## 
## Achieved and Interpolated Confidence Intervals: 
## 
##                   Conf.Level L.E.pt U.E.pt
## Lower Achieved CI     0.8565 -2.400    0.2
## Interpolated CI       0.9500 -2.499    0.2
## Upper Achieved CI     0.9510 -2.500    0.2

結果は変わらない。SIGN.test関数は、タイを無視して計算をしているようである。

次に、無作為化して計算してみます。正規乱数をタイの数だけ発生させて、タイのデータと入れ替える（ただし、順番は入れ替わる）。

tcddpla2 <- c(tcddpla[!tie], rnorm(3)) # 正規乱数を代わりにして入れ替える
tcddfat2 <- c(tcddfat[!tie], rnorm(3)) # 正規乱数を代わりにして入れ替える
SIGN.test(x = tcddpla2, y = tcddfat2, alternative = "two.sided")

## 
##  Dependent-samples Sign-Test
## 
## data:  tcddpla2 and tcddfat2
## S = 7, p-value = 0.2632
## alternative hypothesis: true median difference is not equal to 0
## 95 percent confidence interval:
##  -2.405394  0.200000
## sample estimates:
## median of x-y 
##    -0.7212507 
## 
## Achieved and Interpolated Confidence Intervals: 
## 
##                   Conf.Level  L.E.pt U.E.pt
## Lower Achieved CI     0.8847 -2.4000    0.2
## Interpolated CI       0.9500 -2.4054    0.2
## Upper Achieved CI     0.9586 -2.4061    0.2

無作為化の結果に依存するが、多くの場合、有意でない結果が得られると思われる。

最後に、タイに多い方の符号をつけてみる。

table(tcddpla - tcddfat > 0)

## 
## FALSE  TRUE 
##    14     6

table(tcddpla - tcddfat < 0)

## 
## FALSE  TRUE 
##     9    11

tcddpla3 <- c(tcddpla[!tie], 0, 0, 0)
tcddfat3 <- c(tcddfat[!tie], 1, 1, 1)
SIGN.test(x = tcddpla3, y = tcddfat3, alternative = "two.sided", conf.level = 0.90)

## 
##  Dependent-samples Sign-Test
## 
## data:  tcddpla3 and tcddfat3
## S = 6, p-value = 0.1153
## alternative hypothesis: true median difference is not equal to 0
## 90 percent confidence interval:
##  -2.40000000 -0.03783901
## sample estimates:
## median of x-y 
##            -1 
## 
## Achieved and Interpolated Confidence Intervals: 
## 
##                   Conf.Level L.E.pt  U.E.pt
## Lower Achieved CI     0.8847   -2.4 -0.1000
## Interpolated CI       0.9000   -2.4 -0.0378
## Upper Achieved CI     0.9586   -2.4  0.2000

最も極端な場合でも、有意にならない。従って、両組織間にTCDD濃度に有意な差はみられないと考えられる。

順位の検定

　\(X_1, X_2, ..., X_n\)を連続する分布\(F\)をもつ母集団からの標本とする。分布に依らない検定法の多くは、パラメータ\(\theta\)に対して、あるいは、他のサンプルに対して、標本内の観察値の順位がどのようになるか、ということに基づいている。標本\(X_1,X_2, ..., X_n\)の順位を、\(r(X_1),r(X_2),...,r(X_n)\)として定義する。

　今、連続分布からの標本\(X_1,...,X_n\)の順位が、それぞれ、\(R_i = r(X_i)\)であるとする。このとき、順位\(R_i\)は離散一様分布から得られた確率変数と考えられる。すなわち、 \[ P(R_i = j) = \frac 1n,\; 1 \le j \le n \]

このとき、整数の和の特性から、以下の期待値が得られる。 \[ E(R_i) = \sum_{j=1}^n \frac jn = \frac {n+1}2 \] また、 \[ E(R_i^2) = \sum_{j=1}^n \frac {j^2}n = \frac {n(n+1)(2n+1)}{6n}=\frac {(n+1)(2n+1)}6 \]

したがって、順位の分散\(Var(R_i)\)は、 \[ Var(R_i) = E(R_i^2) - E(R_i)^2 = \frac {n^2-1}{12} \] となる。

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

　ウィルコクソンの符号順位検定は、符号検定よりも検出力の高い検定法である。ウィルコクソンの符号順位検定では、その名の通り、符号と順位を考慮に入れた検定が行われる。

　対になった標本\((X_i, Y_i)\)が観察され、その差を\(D_i = X_i - Y_i\)とする。これを2標本の対比較の問題ととらえ、差\(D\)の真の値が0かどうかを検定することを考える。あるいは、1標本で、中央値の検定に興味があり、\(H_0: med = m_0\)を片側あるいは両側の対立仮説に対して検定することを考える。すなわち、観察値\(X_i\)を\(m_0\)に比較して、差\(D_i = X_i - m_0\)の真の値が0かどうかを検定する。

　ウィルコクソンの符号順位検定では、差の絶対値（\(|D_1|,|D_2|,...,|D_3|\))の順位、\(r(|D_1|),r(|D_2|),...,r(|D_n|)\)をもとに検定を行う。帰無仮説のもとでは、差\(D_i\)は、0に対して対称に分布すると期待され、正の値をとる\(D_i\)の順位和の期待値\(W^+\)と、負の値をとる\(D_i\)の順位和の期待値\(W^-\)は一致する。 \[ W^+ = \sum_{i=1}^n S_i r(|D_i|) \\ W^- = \sum_{i=1}^n (1-S_i) r(|D_i|) \] ここで、\(D_i > 0\)のとき\(S_i = 1\)、\(D_i < 0\)のとき\(S_i = 0\)である。\(D_i = 0\)のときはタイであることを意味しており、無視する。\(W^+ + W^-\)は全順位の和\(\sum_{i=1}^n r(|D_i|)\)であり、タイがない場合は\(\sum_{i=1}^n i = n(n+1) / 2\)に等しくなる。ウィルコクソンの符号順位検定のための統計量は、正の差をもつ順位和と、負の差をもつ順位和間の差として、以下のように表される。 \[ W = W^+ - W^- = 2 \sum_{i=1}^n r(|D_i|)S_i - \sum_{i=1}^n r(|D_i|) = 2 \sum_{i=1}^n r(|D_i|)S_i - n(n + 1) /2 \]

　なお、ウィルコクソンの符号順位検定では、標本数が多い場合には、大数の法則に基づき正規近似が利用される場合がある。標本が同じ集団から得られている場合、正の差をもつ順位の和は負の差をもつ順位の和に近くなる。したがって、帰無仮説\(H_0\)のもとで、タイが無ければ\(E(W) = 0\)となり、\(Var(W) = \sum_i (r(|D_i|)^2) = \sum_i i^2 = n(n+1)(2n+1)/6\)となると期待される。統計量　 \[ Z = \frac W {\sqrt{Var(W)}} \] は、標準正規分布で近似され、正規累積分布関数を用いて、ある特定の対立仮説に対する観察された統計量\(W\)の\(p\)値を計算できる。Rのwilcox.text関数では、サンプル数が50以上になると正規近似が用いられる。

例 一卵性双生児 （Vidakovic 2011より引用）

　本データセットは、12ペアの一卵性双生児が心理学テストをうけ、個人の攻撃性が評価されたものである。双子を互いに比較して、最初に生まれた双子が、あとで生まれた双子よりも、より攻撃的になる傾向があるかどうかを検定したい。心理学テストの結果は以下のとおり（高い値は高い攻撃性を表す）である。

#first-born
fb <- c(86, 71, 77, 68, 91, 72, 77, 91, 70, 71, 88, 87)
# second-born
sb <- c(88, 77, 76, 64, 96, 72, 65, 90, 65, 80, 81, 72)

　仮説は、\(H_0\): 2人の双子の攻撃性得点の平均は同じ、すなわち\(E(X_i) = E(Y_i)\)である。一方、対立仮説\(H_1\)は、最初に生まれた双子は他よりもより攻撃的、すなわち\(E(X_i) > E(Y_i)\)である。\(D_i = X_i - Y_i\)はそれぞれ独立で、対称であると仮定すると、ウィルコクソンの符号順位検定を用いて検定ができる。

wilcox.test(x = fb, y = sb, paired = T, alternative = "greater")

## Warning in wilcox.test.default(x = fb, y = sb, paired = T, alternative =
## "greater"): cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = fb, y = sb, paired = T, alternative =
## "greater"): cannot compute exact p-value with zeroes

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  fb and sb
## V = 41.5, p-value = 0.2382
## alternative hypothesis: true location shift is greater than 0

タイで警告が出る場合は、パッケージcoinの関数を使うと良い。

wilcoxsign_test(fb ~ sb, alternative = "greater")

## 
##  Asymptotic Wilcoxon-Pratt Signed-Rank Test
## 
## data:  y by x (pos, neg) 
##   stratified by block
## Z = 0.78567, p-value = 0.216
## alternative hypothesis: true mu is greater than 0

検定により\(H_0\)は棄却されない。

　ウィルコクソンの符号順位検定は、\(t\)検定のように、\(H_0: \mu = \mu_0\)の検定にも用いられる。ウィルコクソンの符号順位検定で扱われる差\(X_1-Y_1, X_2-Y_2, ..., X_n-Y_n\)の代わりに、\(X_1 - \mu_0, X_2 - \mu_0, ..., X_n - \mu_0\)が用いられる。

例 月の錯覚（Vidakovic 2011より引用）

　Kaufman and Rock (1962)は、地平線付近の月が天頂（頭上で最も高い位置）の月よりも大きく見えるという一般的に観察される事象について検証するために、地平線と天頂に2つの人工的な月を表示する装置を考案した。被験者は、天頂の月の大きさに合わせて地平線の月の大きさを調節するように、またはその逆に調節するように求められた。その結果に基づく、各被験者について、地平線の月の大きさと天頂の月の大きさの比が記録された。1.00であれば錯覚が無いことを意味し、1.00以外であれば錯覚があることを意味する。例えば、1.50という値は、地平線の月は、天頂の月の直径の1.50倍に見えたことを意味する。この錯覚を検定するために、帰無仮説\(H0 : \mu = 1.00\)を棄却し、対立仮説\(H1 : \mu > 1.00\)をとることができるか検定する。

moon <- c(1.73, 1.06, 2.03, 1.40, 0.95, 1.13, 1.41, 1.73, 1.63, 1.56)
mu0 <- 1
wilcox.test(x = moon - mu0, alternative = "greater")

## Warning in wilcox.test.default(x = moon - mu0, alternative = "greater"): cannot
## compute exact p-value with ties

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  moon - mu0
## V = 54, p-value = 0.004002
## alternative hypothesis: true location is greater than 0

警告がでるときは、パッケージcoinの関数を使って計算するとよい。

wilcoxsign_test(moon ~ rep(mu0, length(moon)), alternative = "greater")

## 
##  Asymptotic Wilcoxon-Pratt Signed-Rank Test
## 
## data:  y by x (pos, neg) 
##   stratified by block
## Z = 2.7029, p-value = 0.003437
## alternative hypothesis: true mu is greater than 0

なお、以下のように\(t\)検定を行うと、\(p\)値は、0.0009988となり、帰無仮説\(H_0\)が棄却される。

t.test(x = moon - mu0, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  moon - mu0
## t = 4.2976, df = 9, p-value = 0.0009988
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.26551     Inf
## sample estimates:
## mean of x 
##     0.463

ウィルコクソンの符号順位検定では、\(t\)検定よりも\(p\)値が高いが、いずれの検定においても\(H_0\)は棄却される。

ウィルコクソンの順位和検定、ウィルコクソン・マン・ホイットニー検定

　ウィルコクソンの順位和検定、ウィルコクソン・マン・ホイットニー検定は同等の検定のため、ここでは、ウィルコクソンの順位和検定として紹介する。ウィルコクソンの順位和検定は、比較しようとする母集団が独立で、それぞれが正規分布に従っていない可能性が高い場合に、2標本の\(t\)検定に代わりに用いられる。

　この検定を用いることができるタイプのデータの例は、例えば、アンケートなどで使われる心理検査的回答尺度であるリッカート尺度（例えば、1 = より悪い、2 = 悪い、3 = 変わらない、4 = 良い、5 = より良い、など）のデータである。このようなデータでは、正規分布の仮定が成り立っていないと考えられるため、\(t\)検定を適用するのは適切ではない。ウィルコクソンの順位和検定は、グループが一様かどうか、あるいは、一方のグループが他方に比べて大きいかどうか、という点について検定できる。これを、より一般的な言い方をすると、ウィルコクソンの順位和検定の帰無仮説は、2つの母集団が等しい、というものである。すなわち、\(H_0: F_X(x) = F_Y(x)\)。この検定では、2つの母集団の分布の形が似ていることを仮定しているが、その形状に対する厳密な仮定はなされない。

　ここで、\(\boldsymbol {X} = X_1, ...,X_{n1}\)、\(\boldsymbol {Y} = Y_1, ...,Y_{n2}}\)を、比較しようとしている母集団からのサイズ\(n_1, n_2\)の標本とする。これら標本をまとめて、\(n = n_1 + n_2\)の順位が割り振られるとする。検定統計量\(W_n\)は、最初の標本\(\boldsymbol {X}\)に対応する順位和である。もし\(X_1 = 1, X_2 = 13, X_3 = 7, X_4 = 9, Y_1 = 2, Y_2 = 0, Y_3 = 18\)であるとすると、\(W_n\)の値は、2 + 4 + 5 + 6 = 17となる。

X <- c(1, 13, 7, 9)
Y <- c(2, 0, 18)
is.X <- c(T, T, T, T, F, F, F)
rank(c(X, Y))

## [1] 2 6 4 5 3 1 7

sum(rank(c(X, Y))[is.X])

## [1] 17

　2つの集団が同じ分布をもつとすると、1番目の集団の標本と2番目の集団の標本の順位和は互いに似たような値になり、それぞれの標本サイズに依存する値となりる。ウィルコクソンの順位和検定の統計量は、 \[ W_n = \sum_{i = 1}^n iS_i(\boldsymbol {X}, \boldsymbol {Y}) \] ここで、\(S_i(\boldsymbol {X}, \boldsymbol {Y})\)は、\(i\)番目に順位付けされる観察値が1つめの集団の標本である場合は1、2つめの集団の標本であれば0になる指示関数である。

　例えば、\(X_1 = 1, X_2 = 13, X_3 = 7, X_4 = 9, Y_1 = 2, Y_2 = 0, Y_3 = 18\)では、\(S_1 = 0, S_2 = 1, S_3 = 0, S_4 = 1, S_5 = 1, S_6 = 1, S_7 = 0\)となる。したがって、 \[ W_n = 1\times0 + 2\times1 + 3\times 0 + 4 \times 1 = 5 \times 1 + 6 \times 1 + 7 \times 0 = 2 + 4 + 5 + 6 =17 \] となる。

　統計量\(W_n\)は、タイがなければ、帰無仮説\(H_0\)のもとで \[ E(W_n) = \frac {n_1(n+1)}{2} \] \[ Var(W_n) = \frac {n_1 n_2 (n + 1)}{12} \] となる。

　なお、1番目の集団の標本がいずれも2番目の集団の標本よりも小さな値をもつ場合は、統計量\(W_n\)は最小値をとり、また、その逆の場合には最大値をとる。 \[ \min W_n = \sum_{i=1}^{n_1} i = \frac {n_1(n_1 + 1)}{2} \\ \max W_n = \sum_{i=n - n_1 + 1}^{n} i = \frac {n_1(2n - n_1 + 1)}{2} \]

　統計量\(W_n\)の正規近似は、 \[ W_n \sim N\biggl(\frac{n_1(n+1)}2, \frac {n_1 n_2 (n+1)}{12}\biggr) \] となる。

例 アワノメイガの卵数（Snedecor and Cockran 1989より引用）

　以下は、トウモロコシの植物体20個体に産み付けられていたアワノメイガの卵の数のデータである。草丈（植物体の高さ）が高い植物体ほど、より多くの卵が産み付けられることが知られている。

草丈（植物体の高さ）	卵の数
23インチ未満	0	14	18	0	31	0	0	0	11	0
23インチ以上	36	42	12	32	105	84	15	47	51	65

　ウィルコクソンの順位和検定を用いて、集団の同質性の検定を行ってみる。

low <- c(0, 14, 18, 0, 31, 0, 0, 0, 11, 0)
high <- c(36, 42, 12, 32, 105, 84, 15, 47, 51, 65)
hist(c(low, high))

産み付けられた卵の数が0の植物体も多く、また、多くの卵を産み付けられた植物体もあり、分布は正規分布に従っているとは思えない。そこで、ウィルコクソンの順位和検定を行う。

wilcox.test(x = low, y = high, paired = F)

## Warning in wilcox.test.default(x = low, y = high, paired = F): cannot compute
## exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  low and high
## W = 5, p-value = 0.0006519
## alternative hypothesis: true location shift is not equal to 0

タイがあり警告が出る場合は、coinの関数を使う。

class <- factor(c(rep("low", length(low)), rep("high", length(high))))
wilcox_test(c(low, high) ~ class)

## 
##  Asymptotic Wilcoxon-Mann-Whitney Test
## 
## data:  c(low, high) by class (high, low)
## Z = 3.4473, p-value = 0.0005661
## alternative hypothesis: true mu is not equal to 0

草丈が高い植物体と低い植物体でのアワノメイガの卵数の違いは有意である。

クラスカル・ウォリス検定

クラスカル・ウォリス検定は、ウィルコクソンの順位和検の3つ以上の集団への一般化である。クラスカル・ウォリス検定では、全ての集団が同じ分布関数をもっているという帰無仮説を、少なくとも全集団のうち2つは異なる中央値をもつという対立仮説に対して検定する。

クラスカル・ウォリス検定は、一元配置の分散分析における\(F\)検定に似ている。分散分析では、比較を行っている全ての母集団が、独立の正規分布に従っているという仮定に依存しているが、クラスカル・ウォリス検定では分布の形状に関する仮定がない。

　データは、\(k\)組の独立した無作為標本で、標本サイズは、それぞれ\(n_1, ..., n_k\)で、\(n = n_1 + \cdots + n_k\)とする。帰無仮説のものでは、全\(k\)集団の標本が同じ母集団から得られていると考える。\(X_{ij}\)が集団\(i\)の標本\(j\)の観察値とする。

　集団\(i\)の標本の順位和の期待値\(E(R_i)\)は、1つの観察値の期待順位\((n+1)/2\)の\(n_i\)倍になる。すなわち、 \[n_i(n+1)/2\] である。また、分散は、 \[Var(R_i) = n_i(n+1)(n-n_i)/12\] となる。ここで\(H_0\)を検定するひとつの方法は、\(R_i = \sum_{j = 1}^{n_i} r(X_{ij})\)と 集団\(i\)の標本の順位和の期待値の違いを計算することである。すなわち、統計量 \[ \sum_{i=1}^k \biggl[ R_i - \frac{n_i(n+1)}{2} \biggr]^2 \] は、標本が異なる母集団から得られていれば大きくなると考えられる。したがって、帰無仮説\(H_0\)から考えてこの値が大きすぎる場合には、\(H_0\)を棄却するという方法が考えられる。

　こうした考えをもとに、KruskalとWallis（1952）は、次の統計量を提案した。 \[ H' = \frac{1}{S^2}\biggl[ \sum_{i=1}^k \frac{R_i^2}{n_i} - \frac{n(n+1)^2}{4} \biggr] \] ここで、 \[ S^2 = \frac{1}{n-1}\biggl[ \sum_{i=1}^k \sum_{j=1}^{n_i} r(X_{ij})^2-\frac{n(n+1)^2}{4}\biggr] \] タイがない場合は、上式は、 \[ H = \frac{12}{n(n+1)} \sum_{i=1}^k \frac{1}{n_i} \biggl[ R_i - \frac{n_i(n+1)}{2} \biggr]^2 \] と簡略化される。 \(H'\)は帰無仮説のもとで自由度\(k-1\)の\(\chi^2\)分布で近似される。

　Rでは、kruskal.test関数を用いて検定を行うことができる。今、3つの処理グループがあり、以下の反応が得られていると考える。 \[ (1,3,4),(3,4,5),(4,4,4,6,5) \] このとき、以下のようにしてクラスカル・ウォリス検定を行うことができる。

data   <- c(1,3,4, 3,4,5, 4,4,4,6,5)
belong <- c(1,1,1, 2,2,2, 3,3,3,3,3)
kruskal.test(data, belong)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  data and belong
## Kruskal-Wallis chi-squared = 3.8923, df = 2, p-value = 0.1428

例 プロット間の収量の差（Vidakovic 2011）

　以下のデータは、ある作物の栽培試験における、異なる処理のもとでの収量を計測したものである。ここでは、簡単のため、処理を整数｛1,2,3,4｝で表す。

data <- c(83, 91, 94, 89, 89, 96, 91, 92, 90, 84, 
          91, 90, 81, 83, 84, 83, 88, 91, 89, 
          101, 100, 91, 93, 96, 95, 94, 81, 
          78, 82, 81, 77, 79, 81, 80)
belong <- c(rep(1, 10), rep(2, 9), rep(3, 8), rep(4, 7))
boxplot(data ~ belong)

3番目の処理は残りの処理に比べて大きな値を示し、一方、4番目の処理は残りの処理に比べて小さな値を示す。

では、この収量性のデータについてクラスカル・ウォリス検定を行ってみる。

kruskal.test(data, belong)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  data and belong
## Kruskal-Wallis chi-squared = 20.337, df = 3, p-value = 0.0001445

帰無仮説は0.1%水準で棄却される。

フリードマン検定

　フリードマン（Friedman）検定は、分散分析における乱塊法（randomized block design: RBD）のノンパラメトリック版である。フリードマン検定は、正規性の仮定が疑わしい場合や、母集団間で分散が異なる可能性がある場合に、乱塊法の代わりとして用いられる。この検定でも、検定統計量の計算にはデータの順位を用いられる。なお、フリードマン検定では分布に関する仮定をおかないために、集団が正規分布にしたがっている場合には分散分析に比べて検出力が低くなる。

　乱塊法では、処理の各水準において、各ブロックで反復のある計測が行われる。\(X_{ij}\)が対象\(i\)、処理\(j\)の実験において得られた観察値とする。ここで、\(i = 1, ...,b\)および\(j = 1,...,k\)である。

　検定統計量を構成するために、順位\(\{1,2,...,k\}\)を、表の各行に割り振る。なお、それぞれの観察値にとって、\(H_0\)のもとでの順位の期待値は\((k+1)/2\)である。次に、列ごとに全順位の和をとり、\(R_j = \sum_{i=1}^b r(X_{ij})\), \(1\le j\le k\)を得る。\(H_0\)が真の場合、\(R_j\)の期待値は、\(E(R_j) = b(k+1)/2\)となる。そこで、 \[ \sum_{j=1}^k \biggl[ R_j - \frac{b(k+1)}{2}\biggr]^2 \] は、処理間の違いを評価するための統計量として利用できる。このとき、統計量 \[ S = \frac{12}{bk(k+1)}\sum_{j=1}^k \biggl[ R_j - \frac{b(k+1)}{2}\biggr]^2 \] は、帰無仮説\(H_0\)のもとで自由度\(k- 1\)の\(\chi^2\)分布で近似される。

例 車の能力評価（Vidakovic 2011より引用）

3種類の車（A, B, C）の能力を評価するために、6名のプロドライバー（1-6）が、車A, B, Cを無作為な順番で評価した。評価は車の能力に基づいて行われ、車のメーカー名や、他の外的な因子には影響されないとした。ここに、1から10のスコアをつけられた結果がある。

評価者	1	2	3	4	5	6
車A	7	6	6	7	7	8
車B	8	10	8	9	10	8
車C	9	7	8	8	9	9

Rではfriedman.test関数を用いて検定できる。入力データは、行をブロック、列を処理とした行列で与えられる。

data <- t(matrix(c(7,8,9, 6,10,7, 6,8,8, 7,9,8, 7,10,9, 8,8,9), nrow = 3))
colnames(data) <- c("A", "B", "C")
data

##      A  B C
## [1,] 7  8 9
## [2,] 6 10 7
## [3,] 6  8 8
## [4,] 7  9 8
## [5,] 7 10 9
## [6,] 8  8 9

boxplot(data)

フリードマン検定を行ってみる。

friedman.test(data)

## 
##  Friedman rank sum test
## 
## data:  data
## Friedman chi-squared = 8.2727, df = 2, p-value = 0.01598

車の能力評価の結果は、３種の車の間で有意に異なる。

参考文献

以下は、本テキストを作成するにあたり参考にした書籍です。より詳しく学ぶためにもお薦めの本です。

1. Vidakovic B (2011) Statistics for Bioengineering Sciences With MATLAB and WinBUGS Support. Springer.
2. Snedecor GW and Cochran WG (1989) Statistical Methods. Eighth Edition. Iowa State University Press.