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

アルコール摂取と食道がんの症例対照研究を実施したとする。当初，食道がん患者 200 人と対照 775 人のサン
プリングを行ったが，年齢が交絡している可能性を考え，年齢群別に集計したところ，55 歳未満では，患者群56 人のうちアルコール多量摂取者が 30 人，対照群 472 人のうちアルコール多量摂取者が 64 人で，55 歳以上では，患者群 144 人のうちアルコール多量摂取者が 66 人，対照群 303 人のうちアルコール多量摂取者が 45 人だったとする。
このデータから，アルコール摂取と食道がんには関連があるか，あるとしたらどの程度の関連か評価せよ。まず年齢層別に２つのクロス集計表を作り，別々に分析してから，その結果から判断して，適当な場合には，どの層でも共通した関連があるかどうか，あるとすればどの程度の関連か評価せよ。

µ´

まず array() と dimnames() を使って３次元のクロス集計表を定義する。55 歳未満群について「アルコール多量摂取と食道がんには関連がない」帰無仮説を検定し，55 歳以上群についても同様にする。その結果，どちらの群でも有意水準 1%としても有意な関連があるといえるので，mantelhaen.test() を使って，55 歳未満でも 55 歳以上でも共通した関連があるかどうかをみる。R のコードは下枠内のように
なる。

¶

³

x <- array(c(30,64,26,408,66,45,78,258),dim=c(2,2,2))dimnames(x) <- list(c("患者","対照"),c("アルコール多量","他"),c("55 歳未満","55 歳以上"))
fisher.test(x[,,1])
fisher.test(x[,,2])mantelhaen.test(x)µ´

結果として，55 歳未満群でフィッシャーの直接確率は p = 8.032 × 10−11 と小さく，オッズ比の点推定量は 7.3，95%信頼区間は [3.9,13.8] であった。この結果は，55 歳未満では，アルコール多量摂取者はそうでない者の 7.3 倍，食道がんになりやすいという解釈を与える（もちろん，これは症例対照研究なので，直接いえることは，食道がん患者のアルコール多量摂取曝露オッズは，そうでない者の曝露オッズの7.3 倍だということだけだが）。55 歳以上でも，p = 6.618 × 10−12 ，オッズ比の点推定量は 4.8，95%信
頼区間は [3.0,7.9] なので，55 歳未満群と同様に，アルコール多量摂取をしていると，そうでない場合の何倍も食道がんになりやすいとはいえそうである。
マンテルヘンツェルの検定の結果，χ2M H = 96.8(p < 2.2 × 10−16 ) となり，帰無仮説は棄却されるので，

55 歳未満でも 55 歳以上でも共通した関連があると判断される。共通オッズ比は 5.6，95%信頼区間は[3.9,8.0] なので，年齢調整するとアルコール多量摂取は 5.6 倍，食道がんのリスクを高めると示唆さ
れた。

µ

´

∗1 本資料は http://phi.med.gunma-u.ac.jp/medstat/it11-2005.pdf としてダウンロード可能である。

1

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

1. 標本データ X1 , X2 , ..., Xn が互いに独立に分布 F (x) に従い，別の標本データ Y1 , Y2 , ..., Yn が互いに独立かつ X とも独立で分布 G(y) に従う。

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

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

1. 変数 X のデータを x1 , x2 , ..., xm とし，変数 Y のデータを y1 , y2 , ..., yn とする。

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

2

2. まず，これらをまぜこぜにして小さい方から順に番号をつける∗5 。例えば，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

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

5. もし X と Y に差がなければ，X は N 個のサンプルから偶然によって m 個取り出したものであり，Yがその残りである，と考えることができる。順位についてみると，1, 2, 3, ..., N の順位から m 個の数値を取り出すことになる。同順位がなければ，ありうる組み合わせは，N Cm 通りある∗7 。

6. X > Y の場合には，N Cm 通りのうち，合計順位が RX と等しいかより大きい場合の数を k とする（X < Y の場合は，合計順位が RX と等しいかより小さい場合の数を k とする）。

7. k/N Cm が有意水準 α より小さいときに H0 を疑う。N が小さいときは有意になりにくいが，N が大きすぎると計算が大変面倒である∗8 。そこで，正規近似を行う（つまり，期待値と分散を求めて，統計量から期待値を引いて分散の平方根で割った値が標準正規分布に近似的に従うという関係を用いて検定する）。

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 )から，E(R2 ) = E((m
XRi )2 ) =

i=1

m
X
i=1

E(Ri2 ) + 2

X

E(Ri Rj )i<j∗5 同順位がある場合の扱いは後述する。
∗6 以下説明するように，順位和 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)) で両側検定の正確な有意確率が得られる。
∗7 R では choose(N,m) によって得られる。
∗8 もっとも，今ではコンピュータにやらせればよい。例えば R であれば，wilcox.test(X,Y,exact=T) とすれば，サンプル数の合計が 50 未満で同順位の値がなければ，総当りして正確な確率を計算してくれる。が，つい 15 年くらいまではコンピュータは誰もが使える道具ではなかったし，総当りをするには計算時間がかかりすぎた。今のコンピュータでもサンプルサイズが大きいと，総当りでは計算時間がかかりすぎて実用的でない。

3

となるので∗9 ，E(Ri2 ) = (12 + 22 + ... + N 2 )/N = (N + 1)(2N + 1)/6

と

N

N

k=1

k=1

X
X
1
k}
E(Ri Rj ) ={(k)2 −N (N − 1)
2

=

2

1
N (N + 1)
N (N + 1)(2N + 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. 標準化∗10 して連続修正∗11 し，z0 = {|RX − E(RX )| − 1/2}/ var(RX ) を求める。m と n が共に大=

きければこの値が標準正規分布に従うので，例えば z0 > 1.96 ならば，両側検定で有意水準 5%で有意である。R で有意確率を求めるには，z0 を z0 と書けば，2*(1-pnorm(z0,0,1)) とすればよい。

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

∗9 第１項が対角成分，第２項がそれ以外に相当する。m = 2 の場合を考えてやればわかるが，
2
X

E((

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

i=1

2
X

E(Ri2 ) + 2

X
E(Ri Rj )
i<j

i=1

となる。
∗10 何度も出てくるが，平均（期待値）を引いて分散の平方根で割る操作である。
∗11 これも何度も出てくるが，連続分布に近づけるために 1/2 を引く操作である。
∗12 正確な確率を求めることができれば問題ないけれども，同順位がある場合には，R では正確な確率は求められない。
∗13

var(RX ) = mn/{12N (N − 1)} × {N 3 − N −TX(d3t − dt )}t=1とも書ける。

4

¶
例題 1

³

以下は，外来と入院患者に対して，それぞれ 3 ヶ月治療したときの胃潰瘍の面積変化割合の表であるa 。
人数 面積変化割合（肩数字はその割合を示した人数を意味する）

X：入院患者32−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 の順位和検定を実行してみよう。
a 出典：Doll R, Pygott F: Factors influencing the rate of healing of gastric ulcers: admission to hospital, pheno-barbitone, 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.R は，次の枠内も含んでいるファイルである）のように直接データ入力してもいいが，Excel で入院患者と外来患者を 1 行ずつ使って，横に 32
カラム使って面積変化割合を入力し，まず入院患者の全データを選択してコピーしてから，R に移ってX <- scan("clipboard") として，次に再び Excel に戻って外来患者の全データを選択してコピーしてからR に移って Y <- scan("clipboard") としてデータ入力の部分を済ませる方が便利だろう。
¶
it11-1.R³

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’]µ´RX に X の各値の順位が得られ，RY に Y の各値の順位が得られる。−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 )
より∗14 ，2*(1-pnorm(2.460259,0,1))=0.01388366 となるので∗15 ，有意水準 5%で帰無仮説は棄却され，2∗14 この場合は 2 群のサンプルサイズが等しいので，期待値や分散は RX と RY のどちらを検定統計量にしても変わらないが，連続修正の効く方向が異なるので合計順位の大きい方を検定統計量とする。
∗15 wilcox.test(X,Y,exact=F) の結果の p-value と一致する。

5

群には有意な差があるといえる。この計算部分の R のプログラムは下枠内の通り。

¶³

print(sum(RX))
print(sum(RY))# 見比べると sum(RY) が大きいとわかるので E(RY) と V(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 の順
位和検定をせよ。

µ´

手順は以下の通り。実際に入力して結果をみてみよう。

¶³

data(sleep)
wilcox.test(sleep$extra~sleep$group,exact=F)µ´

exact=F オプションをつけないと，同順位があるので正確な p 値が計算できないという警告メッセージが表示される。サンプルサイズが 50 未満なのでデフォルトで正確な確率を計算しようとするのだけれども同順位の値があるので正確な確率が計算できず，正規近似で計算されているためである。

メディアン検定
検定統計量を計算するために，順位そのものでなくても，大小関係を反映するスコアならば使うことができる。順位が 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

帰無仮説の下では 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 の順位和検定ができれば十分なので，メディアン検定はあまり使われない。

6

¶

³

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.R³

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

符号付き順位和検定
２群間の各サンプルに対応がある場合には，単純な順位和検定よりも切れ味がよい方法がある。符号化順位
検定とも呼ばれるこの方法は，対応のある 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 Ri

i=1

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

Ri E(εi ) =nX

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

i=1

分散は

var(R∗ ) =nX

Ri 2 var(εi )i=1∗16 この正確な確率の計算法は，R. A. Fisher が考案した「並べ換え検定」(permutation test) と呼ばれている。R で実行するにはexactRankTests ライブラリの perm.test() 関数を用いる。

7

=

n
X

Ri 2 (12 × 1/2 + (−1)2 × 1/2)

i=1

=

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.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 の検定の仕組みを箇条書きで説明する。

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

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

• 次に，各群ごとに順位を足し合わせて，順位和 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 を求める。同順位を含むときは，すべての同順位の値について，その個数に個数の２乗から１を引いた値を掛けたものを計算し，その総和を 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.R³

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

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

り１少ないカイ二乗分布に従う統計量 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 オプショ
ンをつけないと警告メッセージがでる）。

¶
it11-5.R³

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

µ

´

課題
10 人の健康な日本人成人男性ボランティアを募り，同じ日の 9:00 と 21:00 に採血をして血清鉄濃度 (mg/L)を測定した結果が下表のように得られたとする（注：架空のデータである）。9:00 と 21:00 の血清鉄濃度に有意差があるといえるか？

ノンパラメトリックな方法で有意水準 5%で検定せよ。結果は配布する紙に学籍番

号，氏名と共に自筆して提出すること。結果の提出をもって出席確認とする。
時刻＼対象者
1
2
3
4
5
6
7
8

9

10

9:00

0.98

0.87

1.12

1.34

0.88

0.91

1.04

1.21

1.17

1.09

21:00

1.03

0.78

1.04

1.52

0.97

0.84

1.32

1.12

1.09

1.32

（注）データは第 6 回の例題と同じなので，http://phi.med.gunma-u.ac.jp/medstat/it06-2.R の先頭 2行をコピーして実行すれば，改めて入力しなくてもよい。

10