医学情報処理演習第 11 回「量的データのノンパラメトリックな分析」∗1
2006 年 12 月 18 日 中澤 港 (nminato@med.gunma-u.ac.jp)
演習サポート web ページ：http://phi.med.gunma-u.ac.jp/medstat/

前回の課題の回答例
まず 55 歳未満と以上で別々にクロス表を入力し，Fisher の直接確率を計算する。ともに有意であり（55 歳未満
p = 3.35 × 10−11 ，55 歳以上 p = 3.01 × 10−7 ），55 歳未満でも 55 歳以上でもアルコール多量摂取と食道がん発生には関連があった。
次に Woolf の検定で３次の交互作用がない（年齢層が違っても関連の向きは同じ）という帰無仮説を検討すると，有意水準 5%で帰無仮説は棄却された（χ2 = 5.47, p = 0.019）。したがって，年齢層によってアルコール多量摂取と食道がん発生の関連性は異なると考えられ，前提が満たされないのでマンテルヘンツェルの共通オッズ比は求められない∗2 。
年齢層別にオッズ比と 95%信頼区間をみると，55 歳未満群で 8.96 (4.43, 18.7)，55 歳以上群で 3.39 (2.07, 5.63) であり，55 歳未満ではアルコール多量摂取すると食道がん罹患リスクが約 9 倍に上昇するのに比べ，55 歳以上では 3 倍余りにとどまっていた。Fourfold プロットを見ても，55 歳未満群の方がアルコール多量摂取がより強く食道がんリス
クを高めると判断される。以上の解析をするコードは次の通り。

¶³

under55 <- matrix(c(30,16,64,308),nc=2)
fisher.test(under55)
over55 <- matrix(c(66,68,45,158),nc=2)
fisher.test(over55)
allbyage <- array(c(under55,over55),dim=c(2,2,2))dimnames(allbyage) <- list(c("大酒飲み","飲まない"),c("患者","対照"),c("55 歳未満","55 歳以上"))
allbyagelibrary(vcd)
woolf_test(allbyage)
fourfold(allbyage)
detach(package:vcd)µ´∗1 本資料は http://phi.med.gunma-u.ac.jp/medstat/it11-2006.pdf としてダウンロード可能である。
∗2 それでも強引にやってしまうと，mantelhaen.test(allbyage) の結果の共通オッズ比は 4.6 で 95%信頼区間は (3.13, 6.76) である。

1

２群の分布の位置の差に関するノンパラメトリックな検定
ノンパラメトリックな検定とは？
パラメータ（parameter）とは母数という意味である。これまで説明してきた検定法の多くは，母数，つまり母集団の分布に関する何らかの仮定をおいていた。その意味で，t 検定も F 検定もパラメトリックな分析法といえる。パラメ
トリックな分析法は，原則としては，その分析法が前提とする，母数に関する仮定を満たすようなデータに対して使うべきである。一方，フィッシャーの正確な確率は母数を仮定しないのでパラメトリックでない。ノンパラメトリックな分析とは，パラメトリックでない分析，つまり母数を仮定しない分析をさし，分布がひどく歪んでいたり，上限や下限があったりするようなデータに対しても使うことができるという利点がある∗3 反面，理想的な場合のパラメトリックな分析に比べると検出力は高くない∗4 。
２群の分布の位置の差に関して問題を定式化すると，次のようになる。

1. 標本データ X1 , X2 , ..., Xn が互いに独立に分布 F (x) に従い，別の標本データ Y1 , Y2 , ..., Yn が互いに独立かつX とも独立で分布 G(y) に従う。
2. F と G には連続分布であるという以外には制約をおかない。
3. このとき，「２つの分布に差はない」という帰無仮説（H0 : F (x) ≡ G(x)）を検定する。
つまり，２群の分布の位置の差についてのノンパラメトリックな検定では，「母数を仮定しない」とは言っても，連続分布であることだけは仮定する。もっとも理想的には分布の形が同じで位置だけがずれているという，「ズレのモデル」が仮定できると話は簡単である。具体的な方法としては，Wilcoxon の順位和検定，符号付順位和検定，符号検定などがある。得られたデータがある種の経験分布関数に一致するかどうかを調べるために良く使われる検定法としては
コルモゴロフ＝スミルノフ検定（KS 検定）がある∗5 。

Wilcoxon の順位和検定
Wilcoxon の順位和検定は，パラメトリックな検定でいえば，t 検定を使うような状況，つまり，独立２標本の分布の位置に差がないかどうかを調べるために用いられる。Mann-Whitney の U 検定と（これら２つほど有名ではないが，Kendall の S 検定とも）数学的に等価である。図示は層別箱ひげ図を用いるのが普通である。
データがもつ情報の中で，単調変換に対して頑健なのは順位なので，これを使って検定しようという発想である。以下，Wilcoxon の順位和検定の手順を箇条書きする。

1. 変数 X のデータを x1 , x2 , ..., xm とし，変数 Y のデータを y1 , y2 , ..., yn とする。
2. まず，これらをまぜこぜにして小さい方から順に番号をつける∗6 。例えば，x8 [1], y2 [2], y17 [3], ..., x4 [N ] のようになる（但し N = m + n）。
3. ここで問題にしたいのは，それぞれの変数の順位の合計がいくつになるかということである。ただし，順位の総合計は (N + 1)N/2 に決まっているので，片方の変数だけ考えれば残りは引き算でわかる。そこで，変数 X だ
け考えることにする。

4. X に属する xi （i = 1, 2, ..., m）の順位を Ri と書くと，X の順位の合計は
RX =

m
X

Ri

i=1
∗3 正規分布に従わないデータに出会ったときは，ノンパラメトリックな検定を行う他に，(1) 対数変換などによりデータの分布を正規分布に近づけるか，(2) 外れ値を吟味して，データそのものにエラーがありそうか，あるいは本質的に異なる母集団からの標本が混ざっていると判断されるなら，それを分析から除外することによって正規分布に近づける，といった方法もあるが，これらの操作が逆にデータを歪めてしまう
危険もあるので，データの性状を細かくみて，他の知見とも考え合わせ，慎重に行わねばならない。また，論文や報告書にまとめる際にも，これらの操作をしたならしたと明記せねばならない。
∗4 ただし，厳密に考えると区分はそれほど明確ではない。例えば，カイ二乗検定では母集団の分布には特定の仮定は置いていないので，定義からすると，実はノンパラメトリックな分析になる。ただし，カイ二乗統計量がカイ二乗分布に従うためにはデータ数が十分に多いことが必要である。もっとも，そう言ってしまえば正規近似する場合の順位和検定もデータ数が多いことが必要なので，問題は何を検定の本質と見なすかという話になってくる。一般には，量的な変数を分析するのに，量の情報を使わずに大小関係，即ち順位の情報だけを使う分析をノンパラ
メトリックな解析と呼ぶと考えておけば大過ない。なお，相関については，既にスピアマンとケンドールの順位相関を説明済みである。
∗5 R では ks.test(変数 1, 変数 2) で実行可能である。
∗6 同順位がある場合の扱いは後述する。

2

となる。RX があまり大きすぎたり小さすぎたりすると，X の分布と Y の分布に差がないという帰無仮説 H0が疑わしいと判断されるわけである。では，帰無仮説が成り立つ場合に，RX はどのくらいの値になるのだろうか？∗7

5. もし X と Y に差がなければ，X は N 個のサンプルから偶然によって m 個取り出したものであり，Y がその残りである，と考えることができる。順位についてみると，1, 2, 3, ..., N の順位から m 個の数値を取り出すことになる。同順位がなければ，ありうる組み合わせは，N Cm 通りある∗8 。
6. X > Y の場合には，N Cm 通りのうち，合計順位が RX と等しいかより大きい場合の数を k とする（X < Y の場合は，合計順位が RX と等しいかより小さい場合の数を k とする）。
7. k/N Cm が有意水準 α より小さいときに H0 を疑う。N が小さいときは有意になりにくいが，N が大きすぎると計算が大変面倒である∗9 。そこで，正規近似を行う（つまり，期待値と分散を求めて，統計量から期待値を引いて分散の平方根で割った値が標準正規分布に近似的に従うという関係を用いて検定する）。

8. 帰無仮説 H0 のもとでは，期待値はE(R) =m
X

E(Ri ) = m(1 + 2 + ... + N )/N = m(N + 1)/2i=1

（1 から N までの値を等確率 1/N でとるから）。分散はちょっと面倒で，var(R) = E(R2 ) − (E(R)2 )
から，

m
m
X
X
X
E(R2 ) = E((
Ri )2 ) =E(Ri2 ) + 2
E(Ri Rj )i=1となるので∗10 ，

i=1

i<j

E(Ri2 ) = (12 + 22 + ... + N 2 )/N = (N + 1)(2N + 1)/6とE(Ri Rj ) =

N

N

k=1

k=1

X
X
1
{(
k)2 −
k}
N (N − 1)N 2 (N + 1)2
N (N + 1)(2N + 1)
1
(
−
)
=
N (N − 1)
4
6(N + 1)(3N + 2)
12を代入して整理すると，結局，var(RX ) = m(N + 1)(N − m)/12 = mn(N + 1)/12 となる。
p
9. 標準化∗11 して連続修正∗12 し，z0 = {|RX − E(RX )| − 1/2}/ var(RX ) を求める。m と n が共に大きければこの値が標準正規分布に従うので，例えば z0 > 1.96 ならば，両側検定で有意水準 5%で有意である。R で有意確率を求めるには，z0 を z0 と書けば，2*(1-pnorm(z0,0,1)) とすればよい。
=∗7 以下説明するように，順位和 R をそのまま検定統計量として用いるのが Wilcoxon の順位和検定であり，R

X , RY の代わりに，UX =mn + n(n + 1)/2 − RY ，UY = mn + m(m + 1)/2 − RX として，UX と UY の小さいほうを U として検定統計量として用いるのが，Mann-Whitney の U 検定である。また，UX − UY を検定統計量とするのが Kendall の S 検定である。有意確率を求めるために参照する表が異なる（つまり帰無仮説の下で検定統計量が従う分布の平均と分散は，これら３つですべて異なる）が，数学的には等価な検定である。R では，Wilcoxon の順位和統計量の分布関数が提供されているので，例えばここで得られた順位和を RS と書くことにすると，
2*(1-pwilcox(RS,m,n)) で両側検定の正確な有意確率が得られる。
∗8 R では choose(N,m) によって得られる。
∗9 もっとも，今ではコンピュータにやらせればよい。例えば R であれば，wilcox.test(X,Y,exact=T) とすれば，サンプル数の合計が 50 未
満で同順位の値がなければ，総当りして正確な確率を計算してくれる。が，つい 15 年くらいまではコンピュータは誰もが使える道具ではな
かったし，総当りをするには計算時間がかかりすぎた。今のコンピュータでもサンプルサイズが大きいと，総当りでは計算時間がかかりすぎて実用的でない。
∗10 第１項が対角成分，第２項がそれ以外に相当する。m = 2 の場合を考えてやればわかるが，
2
X

E((

Ri )2 ) = E((R1 + R2 )2 ) = E(R12 + R22 + 2R1 R2 ) =

i=1

2
X
i=1

となる。
∗11 何度も出てくるが，平均（期待値）を引いて分散の平方根で割る操作である。
∗12 これも何度も出てくるが，連続分布に近づけるために 1/2 を引く操作である。

3

E(Ri2 ) + 2

X
E(Ri Rj )i<j10. ただし，同順位があった場合は，ステップ 2 の「小さい方から順に番号をつける」ところで困ってしまう。例えば，変数 X が {2, 6, 3, 5}，変数 Y が {4, 7, 3, 1} であるような場合には，X にも Y にも 3 という値が含まれる。こういう場合は，下表のように平均順位を両方に与えることで，とりあえず解決できる。
属する変数 YXX
Y Y X X Y
値

1

2

3

3

4

5

6

7

順位
1
2 3.5 3.5
5
6
7
8
11. ただし，このやり方では，正規近似をする場合に分散が変わる∗13 。帰無仮説の下で，E(RX ) = m(N + 1)/2 はステップ 8 と同じだが，分散がvar(RX ) = mn(N + 1)/12 − mn/{12N (N − 1)} ·TX(d3t − dt )t=1となる∗14 。ここで T は同順位が存在する値の総数であり，dt は t 番目の同順位のところにいくつのデータが重なっているかを示す。上の例では，T = 1，d1 = 2 となる。なお，あまりに同順位のものが多い場合は，この程度の補正では追いつかないので，値の大小があるクロス集計表として分析することも考慮すべきである（例えば
第 9 回に説明した Cochran-Armitage 検定などが考えられる）。

¶
例題 1

³

以下は，外来と入院患者に対して，それぞれ 3 ヶ月治療したときの胃潰瘍の面積変化割合の表であるa 。
人数 面積変化割合（肩数字はその割合を示した人数を意味する）−100(12) , −93, −92, −91(2) , −90, −85, −83, −81, −80, −78, −46, −40,
−34, 0, 29, 62, 75, 106, 147, 1321
Y: 外来患者
32−100(5) , −93, −89, −80, −78, −75, −74, −72, −71, −66, −59, −41,
−30, −29, −26, −20, −15, 20, 25, 37, 55, 68, 73, 75, 145, 146, 220, 1044
入院患者と外来患者の間で，面積変化割合に差があるといえるだろうか？ Wilcoxon の順位和検定を実行してみ
X：入院患者32

よう。
a 出典：Doll R, Pygott F: Factors influencing the rate of healing of gastric ulcers:

admission to hospital, phenobarbitone,
and ascorbic acid. Lancet, i: 171-175, 1952. から Armitage P, Berry G, Matthews JNS: Statistical Methods in Medical
Research 4th ed. Blackwell Publishing, 2002, pp.281. に引用されているものµ´

まず，2 群を併せて順位を計算する。

¶
it11-1-2006.R(1)³

X <- c(rep(-100,12), -93, -92, rep(-91,2), -90, -85, -83, -81, -80, -78, -46, -40,
-34, 0, 29, 62, 75, 106, 147, 1321)
Y <- c(rep(-100,5), -93, -89, -80, -78, -75, -74, -72, -71, -66, -59, -41, -30, -29,
-26, -20, -15, 20, 25, 37, 55, 68, 73, 75, 145, 146, 220, 1044)
d <- data.frame(gr=factor(c(rep(’X’,length(X)),rep(’Y’,length(Y)))), val=c(X,Y))
rnk <- rank(d$val)
dd <- data.frame(d,rk=rnk)
RX <- dd$rk[dd$gr==’X’]
RY <- dd$rk[dd$gr==’Y’]µ´

上枠内のように直接データ入力してもいいが，Excel で入院患者と外来患者を 1 行ずつ使って，横に 32 カラム使って面積変化割合を入力し，まず入院患者の全データを選択してコピーしてから，R に移って X <- scan("clipboard") として，次に再び Excel に戻って外来患者の全データを選択してコピーしてから R に移って Y <- scan("clipboard")
としてデータ入力の部分を済ませる方が便利だろう。RX に X の各値の順位が得られ，RY に Y の各値の順位が得られ∗13 正確な確率を求めることができれば問題ないけれども，同順位がある場合には，R では正確な確率は求められない。
∗14

var(RX ) = mn/{12N (N − 1)} × {N 3 − N −

T
X
t=1

とも書ける。

4(d3t − dt )}る。−100 はすべて 9 位となる（1 位から 17 位までが −100 なので平均順位の 9 位が与えられる）ので，RX や RY を表示させても最初の方の値は 9.0 ばかりである。

RX と RY それ ぞれの 合計は，sum(RX) より RX = 858，sum(RY) より RY = 1222 となる。RY の合計の方が大きいのでそちらを検定統計量として採用し，同順位がある場合の正規近似から，RY の期待値 E(RY ) は，E(RY )=32 × (32 + 32 + 1)/2 = 1040 となり，分散 V(RY ) は，V(RY )=32 × 32 × (32 + 32 + 1)/12 − 32 × 32/{12 ×(32 + 32 − 1) × (32 + 32)} × {173 − 17 + (23 − 2) × 5} = 5442.413 となるので，
RY − E(RY ) − 1/2p= 2.460259
V (RY )
より∗15 ，2*(1-pnorm(2.460259,0,1))=0.01388366 となるので∗16 ，有意水準 5%で帰無仮説は棄却され，2 群には有意な差があるといえる。この計算部分の R のプログラムは下枠内の通り。

¶
it11-1-2006.R(2)³

print(sum(RX))
print(sum(RY))
N <- 32+32+1
print(ERY <- 32*N/2)
print(VRY <- 32*32*N/12-32*32/(12*(N-2)*(N-1))*(17^3-17+(2^3-2)*5))
print(z0 <- (sum(RY)-ERY-0.5)/sqrt(VRY))
2*(1-pnorm(z0,0,1))

µ
¶
例題 2

´
³

R の組み込みデータ sleep は，20 人の患者を 10 人ずつ 2 群に分けて，それぞれ異なる催眠剤を与えたときに睡眠時間が何時間長くなったかという値をもつ。変数 extra が睡眠時間の増分を示し，group が催眠剤の異なる 2群を示す要因型の変数である。2 群間で睡眠時間の増加量に差はあるか Wilcoxon の順位和検定をせよ。
µ´
手順は以下の通り。実際に入力して結果をみてみよう。

¶³

attach(sleep)
boxplot(extra~group)wilcox.test(extra~group,exact=F)
detach(sleep)µ´

ここでも exact=F オプションをつけないと，同順位があるので正確な p 値が計算できないという警告メッセージが表示される。サンプルサイズが 50 未満なのでデフォルトで正確な確率を計算しようとするのだけれども同順位の値があるので正確な確率が計算できず，正規近似で計算されているためである。W = 25.5, p = 0.069 より，2 群間に睡眠時間の増加量の差がないという帰無仮説は有意水準 5%で棄却できないので，差があるとはいえない。

メディアン検定
検定統計量を計算するために，順位そのものでなくても，大小関係を反映するスコアならば使うことができる。順位が i 番目のオブザーベーションのスコア s(Ri ) として，順位そのものを使う代わりに，i ≥ [(N + 1)/2] のとき 1，

i < [(N + 1)/2] のとき 0 を用いるのがメディアン検定である。次のように言い換えることもできる。
m 個のデータからなる X と n 個のデータからなる Y を合わせた N = m + n 個のデータを，全体のメディアン以上
かメディアンより小さいかによって分類すると，以下の２×２クロス集計表が得られる。

メディアン以上
メディアンより小さい
合計

X
H
m−H
m

Y
(m + n/2) − H
H + (n − m)/2n

合計

(m + n)/2(m + n)/2m+n∗15 この場合は 2 群のサンプルサイズが等しいので，期待値や分散は RX と RY のどちらを検定統計量にしても変わらないが，連続修正の効く方向が異なるので合計順位の大きい方を検定統計量とする。
∗16 wilcox.test(X,Y,exact=F) の結果の p-value と一致する。exact=F は，有意確率の計算を正確な確率でなく正規近似で行うオプション指定である。このオプションがないと，サンプルサイズが 50 未満の場合，順位の組み合わせを使った正確な確率が計算される。

5

帰無仮説の下では H は m/2 の周りに分布する（超幾何分布）ので，

P r(H = h0 ) =n Ch0 ·n C(m+n)/2−h0 /m+n C(m+n)/2 より，P r(H ≥ h0 ) をすべて合計して２倍すれば，両側検定での有意確率が得られる。下枠内のように関数定義しておけば，例えば，median.test(rnorm(100),rnorm(100)) のよ
うに２つのベクトルを与えることによって，簡単にメディアン検定をすることができる。ただし，通常はデータが歪んでいても Wilcoxon の順位和検定ができれば十分なので，メディアン検定はあまり使われない。

¶³

median.test <- function(X,Y) {M <- median(c(X,Y))
fisher.test(cbind(table(X>=M),table(Y>=M))) }

µ
¶
例題 3

´
³

例題 2 で取り上げたデータ sleep の２群間でメディアン検定をせよ。

µ´

上で定義した関数では X~C 型の引数は与えられないので，上のように定義した後でmedian.test(sleep$extra[sleep$group==1],sleep$extra[sleep$group==2]) とするか，または，下枠内のようにする。得られる p-value が約 0.18 なので，有意水準 5%では有意差はないと判断できる。

¶
it11-2-2006.R³

attach(sleep)
X <- subset(extra,group==1,drop=T)
Y <- subset(extra,group==2,drop=T)
M <- median(extra)
fisher.test(cbind(table(X>=M),table(Y>=M)))
detach(sleep)µ´

符号付き順位和検定
２群間の各サンプルに対応がある場合には，単純な順位和検定よりも切れ味がよい方法がある。符号化順位検定とも呼ばれるこの方法は，対応のある t 検定の場合と同じような考え方に基づく。
変数 X の任意の i 番目（i は 1 から n までの整数値）のデータが xi = ei + θi のように，誤差変動 ei と真の効果 θiの和であると捉えれば，もし X と Y が同じ母集団からのサンプルであるならば X − Y により X と Y に共通する真の効果を打ち消すことができて，Ui = xi − yi = ei − e0i が得られる。このとき帰無仮説は，ei と e01 の分布が同じということなので，Ui は原点に対して対称になるはずである。そこで，Ui の絶対値が小さい方から順に順位 Ri をつける。さらに，εi = 1(Ui > 0), εi = −1(Ui < 0) とすれば，帰無仮説の下で P r(εi = 1) = P r(εi = −1) = 1/2 となる。
いま，

n
X

R∗ =

εi R ii=1とおけば，R∗ の大きさによって検定ができる。
すべての場合（εi の値が各 i について 2 通りあるので，2n 通り）を計算してやれば正確な確率が計算できるが∗17 ，n が大きくなると計算が大変なので，n ≥ 15 ならば近似を行ってよいことになっている。R∗ の期待値はE(R∗ ) =nX

Ri E(εi ) =

i=1

n
X

Ri (1 × 1/2 + (−1) × 1/2) = 0

i=1

分散は

var(R∗ ) =nX

Ri 2 var(εi )

i=1

=

n
X

Ri 2 (12 × 1/2 + (−1)2 × 1/2)i=1
∗17 こ の 正 確 な 確 率 の 計 算 法 は ，R.

A. Fisher が 考 案 し た「 並 べ 換 え 検 定 」(permutation test) と 呼 ば れ て い る 。R で 実 行 す る に はexactRankTests ライブラリの perm.test() 関数を用いる。

6

n
X

=

Ri 2 = n(n + 1)(2n + 1)/6i=1となるので，標準化と連続性の補正をして，|R∗ | − 1/2pvar(R∗ )が標準正規分布に従うことを利用して検定する。なお，R では，対応のある２群の生のデータを X と Y に付値しておき，wilcox.test(X,Y,paired=T) とすればこの検定ができる。

¶
例題 4

³

以下は，躁鬱の患者９人についてのハミルトンの抑うつ尺度の測定値である。治療前と精神安定剤で治療を開始した後の値を示すa 。治療によって抑うつが改善されたといえるか，Wilcoxon の符号付き順位和検定で検定せよ。
治療前
1.83
0.50
1.62 2.48 1.68 1.88 1.55 3.06 1.30
治療後 0.878 0.647 0.598 2.05 1.06 1.29 1.06 3.14 1.29
a 出典: Hollander M, Wolfe DA (1973) Nonparametric statistical inference, New York: John Wiley & Sons, pp.27-33.

µ´

下枠内のように入力すればよい。注意すべきことは，治療後の方が抑うつ尺度は改善していると期待されるので，対立仮説を「X が Y より大きい」とする片側検定にすべきだということである。なお，対応のある t 検定と同じく，差を計算して，その期待値がゼロより大きいという対立仮説を検定する形にしても同等である。

¶
it11-3-2006.R³

X <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30)
Y <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)
wilcox.test(X, Y, paired=T, alt="greater")# wilcox.test(X-Y, alt="greater") と同じµ´

多群間の分布の位置の差の検定
クラスカル＝ウォリス (Kruskal-Wallis) の検定
ノンパラメトリックな分析でも，多群間の分布の位置の差を調べるためのアプローチとしては，パラメトリックな場
合と同じく大別して２つある。１つは，一元配置分散分析と同じく，群分け変数が量的な変数に与える効果をみるという形にする方法で，もう１つは，２群ずつすべての組み合わせについて検定の多重性を調整して Wilcoxon の順位和検
定を繰り返す方法である。
前者を行う代表的な方法がクラスカル＝ウォリス (Kruskal-Wallis) の検定である。量的な変数 X，群分け変数 C にデータが付値されているとして，kruskal.test(X~C) とすれば実行できる。なお，一元配置分散分析を行うときに前もってバートレットの検定によって分散の同等性をしておく必要があったのと同様，厳密に言えば，Fligner-Killeen の検定 (fligner.test(X~C)) によって，前もってばらつきの同等性を検定しておかねばならない。以下，Kruskal-Wallisの検定の仕組みを箇条書きで説明する。

•「少なくともどれか１組の群間で大小の差がある」という対立仮説に対する「すべての群の間で大小の差がない」という帰無仮説を検定する。

• まず２群の比較の場合の順位和検定と同じく，すべてのデータを込みにして小さい方から順に順位をつける（同順位がある場合は平均順位を与える）。

• 次に，各群ごとに順位を足し合わせて，順位和 Ri (i = 1, 2, ..., k; k は群の数) を求める。
• 各群のオブザーベーションの数をそれぞれ ni とし，全オブザーベーション数を N としたとき，各群について統計量 Bi を Bi = ni {Ri /ni − (N + 1)/2}2 として計算し，
B=

k
X

Bi

i=1

として B を求め，H = 12 · B/{N (N + 1)} として H を求める。同順位を含むときは，すべての同順位の値に7

ついて，その個数に個数の２乗から１を引いた値を掛けたものを計算し，その総和を A として，

H0 =H
1 − N (NA2 −1)

により H を補正した値 H 0 を求める。

• H または H 0 から表を使って（データ数が少なければ並べかえ検定によって）有意確率を求めるのが普通だが，k ≥ 4 で各群のオブザーベーション数が最低でも 4 以上か，または k = 3 で各群のオブザーベーション数が最低でも 5 以上なら，H や H 0 が自由度 k − 1 のカイ二乗分布に従うものとして検定できる。
¶
例題 5

³

R の組み込みデータ chickwts（詳細は第７回の資料参照）を使って，餌の種類 (feed) によって，６週間飼育後の鶏の体重 (weight) に差がでるかを，ノンパラメトリックな方法で分析せよ。
µ´
上で説明したとおり，次の枠内のように入力すればよい。ノンパラメトリックな位置母数の比較のためのグラフ表示
としては，層別箱ヒゲ図を使うのが普通である（２群の場合と同様）。

¶
it11-4-2006.R³

attach(chickwts)
boxplot(weight~feed)
fligner.test(weight~feed)kruskal.test(weight~feed)
detach(chickwts)µ´

フリードマン (Friedman) の検定
対応のある多群間の差をノンパラメトリックな方法で調べるには，フリードマン (Friedman) の検定と呼ばれる手法
を用いた方が切れ味がよい。簡単に説明すると，まず同じ個体について群間で順位をつける（群といっても，対応がある場合だから，例えば 2005 年の予測値と 2010 年の予測値と 2025 年の予測値というように，個々の個体について順位をつけることが可能である）。次に，群ごとに順位の合計（順位和）を計算する。順位和の二乗和から順位和の平均の二乗を引いた値を統計量 S として，サンプル数が少ない場合は表によって（コンピュータシミュレーションによってもよい）有意確率を計算し，サンプル数が多い場合は自由度が群数より１少ないカイ二乗分布に従う統計量 Q を S の

12 倍を個体数と群数と「群数＋１」の積で割った値として計算して有意確率を計算する。ただし同順位がある場合は調整が必要であり，煩雑なので，通常はコンピュータソフトウェアに計算させる。R では friedman.test() 関数を用いる。使い方は example(friedman.test) を参照されたい。

多重比較多群があるときに２群ずつのすべての組み合わせについて分布の位置の差を検定するには，パラメトリックな分析でやったのと同様，多重比較を行って第一種の過誤を調整する必要がある。第一種の過誤の調整法としては，Holm の方
法と Bonferroni の方法は問題なく使えるが，Tukey の HSD のように母集団の分布を仮定した方法は使えない。R では量的変数を X，群分け変数を C として，pairwise.wilcox.test(X,C) 関数を用いればよい（X~C ではなく X,C であることに注意）。なお，サンプルサイズが小さいときは正確な確率を求めるのがデフォルト動作だが，同順位の値があると正確な確率を求めることができず，警告メッセージが表示される。警告を見たくないときは exact=F オプションをつければよい。

¶
例題 6

³

例題 5 の chickwts データを使って，どの餌とどの餌の間で体重が異なるかを，検定の多重性を調整して検定
せよ。

µ

´

上で説明したとおり，下枠内のように入力すればよい（同順位の値が含まれているので exact=F オプションをつけないと警告メッセージがでる）。

8¶it11-5-2006.R³

attach(chickwts)
pairwise.wilcox.test(weight,feed,exact=F)
detach(chickwts)

µ

´

課題
http://phi.med.gunma-u.ac.jp/medstat/p11.txt は，ある途上国の３つの地域（変数 GRP，離島にある P 村，
首都から車で約 1 時間離れた M 村，首都 H 市）の，再生産をおえたカップル 30 組ずつが生涯に産んだ子供の数（変数 PARITY）のデータである（架空のものである）。これら３地域間で，生涯子供数には互いに有意差があるか，ノンパ
ラメトリックな方法で調べよ。作図をしてパワーポイントに貼り付け，検定結果とその解釈をテキストボックスで打ち込み，印刷して氏名を自筆して提出せよ。

9