はじめに
ある事象について観察された度数分布が、何らかの理由をもとに期待される度数分布に適合するかどうか調べたいときがある。例えば、サイコロが出す目が期待されている確率に従っているかどうかについて調べたい場合などが考えられる。あるいは、交配実験を行い、後代における形質の分離がメンデルの法則に従っているか否かを調べたい場合もある。期待される度数分布は、サイコロの例ように、階級（1〜6の目）による差がないと期待される場合は一様分布、サイコロを複数回ふってある目の出現回数を調べるようなベルヌーイ試行の場合には二項分布、多数のサイコロのゾロ目の出現など、まれにおこる現象の場合にはポアソン分布、多数のサイコロの目の和や平均のように、中心極限定理があてはまる状況であれば正規分布、といったように様々な分布が考えられる。また、交配実験の後代における分離のように、遺伝学的な推論から期待される分布がある場合もある。観察された分布と理論的に期待される分布を比較し、その適合度（goodness of fit）が明らかにすることができれば、観察されている現象の背後の仕組みを理解することにもつながる。本講義では、まず、適合度検定（goodness of fit test）とよばれる方法について説明する。
なお、観察分布と期待分布の適合度検定は、要因が2つ以上ある場合にも適応できる。例えば、ある薬剤による処理の有無と種子の発芽数に何らかの関係があるのか、あるいは、両者は独立なのかを調べたい場合などが考えられる。このように、ある要因と別の要因が独立なのかどうかを調べたい場合、これら要因について表にまとめた分割表（contingency table）を作り、この分割表における各セルの度数が期待されている度数に適合しているかどうかを検定する。こうした検定を行うことで要因間の独立性を明らかにし、何らかの関係がある（非独立的である）場合には、その関係についてさらなる研究を進めることも可能となる。本講義の後半では、分割表を用いた検定について説明する。
適合度検定
ある母集団から大きさ\(n\)の標本を抽出して、\(r\)種のカテゴリに分類したとき、観察度数（observed frequency）がそれぞれ\(n_1,n_2,…,n_r\) (\(n_1,+n_2+⋯+ n_r=n\))であったとする。いっぽう、理論的根拠からそれぞれのカテゴリは確率\(p_1,p_2,…,p_r\) (\(p_1,+p_2+⋯+ p_r=1\))をもっておこると期待されるとする。確率から期待される頻度\(np_1,np_2,…,np_r\)を理論度数（theoretical frequency）あるいは期待度数（expected frequency）とよぶ。ここでは、観察度数が期待度数に適合しているかどうかを検定する適合度検定について説明する。なお、カテゴリへの分類は、サイコロの目のように離散的な場合だけでなく、連続的な変数をある閾値で分けた場合にも適用できる。
適合度検定には、Karl Pearsonによって提案されたカイ二乗検定（chi-square test）がよく用いられる。カイ二乗検定では、観察度数と期待度数が適合している場合に、統計量 \[ \chi^2 = \sum_{i=1}^r \frac{(n_i - np_i)^2}{np_i} \] が、自由度\(r-1\)のカイ二乗分布に従うことを利用する。具体的には、計算された\(\chi^2\)値について、自由度\(r-1\)のカイ二乗分布における上側確率を求める。これが有意水準\(\alpha\)以下である場合には、「各カテゴリの事象が生じる確率は、理論などから期待される分布に適合している」という帰無仮説を棄却する。有意水準\(\alpha\)より大きい場合は、帰無仮説を採択する。
1つのサイコロを60回ふったときの以下の結果について、適合度検定を行ってみる。観察度数を見ると、1や5がたくさん生じており、サイコロの正しさに疑念が生じる。なお、サイコロが正しい場合には、1から6までの目がでる確率はそれぞれ等しく1/6になると期待されるため、サイコロを60回ふったときの期待度数は、いずれの目も\(60\times 1/6=10\)である。
表1. サイコロを60回ふる実験における各目の観察度数と期待度数
| 1 | 2 | 3 | 4 | 5 | 6 | 計 | |
|---|---|---|---|---|---|---|---|
| 観察度数(O) | 14 | 5 | 9 | 7 | 17 | 8 | 60 |
| 期待度数(E) | 10 | 10 | 10 | 10 | 10 | 10 | 60 |
期待度数を、Rを用いて計算し、度数分布図を描いてみよう。
# 60回サイコロを振った際の出目の集計
ni <- c(14, 5, 9, 7, 17, 8)
# 期待度数を計算
npi <- rep(10, 6)
# グラフを描く
data <- data.frame(ni, npi, row.names = 1:6)
barplot(t(data), beside = T)
図1. サイコロを60回ふったときの観察度数（濃）と期待度数（淡）
では、このデータをもとに、サイコロが正しいかどうか、いいかえると、帰無仮説「各カテゴリの事象が生じる確率が、想定した期待分布に適合する」が棄却されないかどうかを、カイ二乗検定で確認してみる。まずは、観察度数と期待度数の差 \(n_i-np_i\) を計算し、それをもとに \((n_i-np_i)^2/np_i\)を計算し、全カテゴリ分の和をとる。
表2. サイコロを60回ふったデータのカイ二乗検定の計算過程
| 1 | 2 | 3 | 4 | 5 | 6 | 計 | |
|---|---|---|---|---|---|---|---|
| 観察度数(O) | 14 | 5 | 9 | 7 | 17 | 8 | 60 |
| 期待度数(E) | 10 | 10 | 10 | 10 | 10 | 10 | 60 |
| O – E | 4 | -5 | -1 | -3 | 7 | -2 | |
| (O – E)2 / E | 1.6 | 2.5 | 0.1 | 0.9 | 4.9 | 0.4 | 10.4 |
これを、Rを用いて計算してみる。
ni - npi## [1] 4 -5 -1 -3 7 -2# （O - E)^2 / Eを計算
v <- (ni - npi)^2 / npi; v## [1] 1.6 2.5 0.1 0.9 4.9 0.4# 和をとって、自由度5のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 10.4得られた10.4という値について、自由度5のカイ二乗分布での上側確率を計算する。上側確率の計算には、カイ二乗分布の累積分布関数pchisqを用いる。
1 - pchisq(ch2, length(ni) - 1)## [1] 0.06466303帰無仮説は、有意水準5%で棄却されないことが分かる。いいかえると、観察度数の期待度数からのずれは有意でなく、サイコロはとくにどの目が出やすい（あるいは出にくい）という傾向があるとはいえない、と結論される。
なお、ここまでの一連の解析は、chisq.test関数を用いて行うこともできる。
# Rの備え付け関数を用いる
chisq.test(ni)##
## Chi-squared test for given probabilities
##
## data: ni
## X-squared = 10.4, df = 5, p-value = 0.06466次に、12個のサイコロを一度にふる実験を26,306回（！）繰り返し、5または6が出た回を数えたWeldonのサイコロデータについて解析を行ってみる。
表3. Weldonのサイコロデータ
| 12個のサイコロのうち、5または6が出たものの数 | 観察度数（O） |
|---|---|
| 0 | 185 |
| 1 | 1149 |
| 2 | 3265 |
| 3 | 5475 |
| 4 | 6114 |
| 5 | 5194 |
| 6 | 3067 |
| 7 | 1331 |
| 8 | 403 |
| 9 | 105 |
| 10 | 14 |
| 11 | 4 |
| 12 | 0 |
サイコロを1個ふって5または6がでる確率は1/3である。サイコロを同時に12個ふり、5または6となるサイコロの数は二項分布\(Bin(12,1/3)\)に従う。なお、\(Bin(n,p)\)は、パラメータ\(n, p\)をもつ二項分布 \[ P[X = k] = \left(\begin{array}{cc} n\\k\end{array}\right) p^k (1-p)^{n-k} \] を表す。したがって、期待度数は、実験の総繰り返し数26,306と\(Bin(12,1/3)\)の積として計算される。
Rを用いて期待度数を計算し、観察度数と期待度数の度数分布図を描いてみる。
# Weldonのサイコロデータ（12個のサイコロを振って、5または6が出た数)
ni <- c(185, 1149, 3265, 5475, 6114, 5194, 3067, 1331, 403, 105, 14, 4, 0)
# 期待値 (p = 1/3, n = 12の二項分布に従う)
pi <- dbinom(0:12, 12, 1/3)
(npi <- sum(ni) * pi)## [1] 2.027495e+02 1.216497e+03 3.345366e+03 5.575610e+03 6.272561e+03
## [6] 5.018049e+03 2.927195e+03 1.254512e+03 3.920351e+02 8.711891e+01
## [11] 1.306784e+01 1.187985e+00 4.949938e-02# グラフを描く
data <- data.frame(ni, npi, row.names = 0:12)
barplot(t(data), beside = T)
図2. Weldonのサイコロデータ (12個のサイコロを振って、5または6が出た頻度)
度数分布図をみると、観察度数は期待度数によく適合しているようにみえる。カイ二乗検定を用いて適合しているかどうか検定してみる。
# 観察度数(O)と期待度数(E)の差
ni - npi## [1] -17.74946043 -67.49676258 -80.36609708 -100.61016181 -158.56143203
## [6] 175.95085438 139.80466505 76.48771359 10.96491050 17.88109122
## [11] 0.93216368 2.81201488 -0.04949938# （O - E)^2 / Eを計算
(v <- (ni - npi)^2 / npi)## [1] 1.55385541 3.74502678 1.93064357 1.81547927 4.00820749 6.16946990
## [7] 6.67715753 4.66346196 0.30667985 3.67008067 0.06649373 6.65616728
## [13] 0.04949938# 和をとって、自由度r-1のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 41.312221 - pchisq(ch2, length(ni) - 1)## [1] 4.344864e-05観察度数の期待度数からのずれは高度に有意である。帰無仮説は0.1%水準でも棄却される。したがって、実験に用いられたサイコロの目の出やすさには偏りがある、と結論される。観察頻度と期待頻度の差を見直すと、小さな数のときに観察度数が期待度数を下回り、大きな数のときに上回っていることが分かる。これは、5や6が出やすいような偏りがあったことを示唆している。Pearsonは、Galtonの伝記の脚注の中で、目のひとつひとつが彫られているサイコロは、それによって、5や6が出やすいのではないかと推測している（Labby 2009, Chance 22: 6-13）。なお、特殊な機械を使って実験をした動画もある。
なお、観察度数と期待確率をもとに、chisq.testを使って、次のように計算することもできる。
chisq.test(ni, p = pi)## Warning in chisq.test(ni, p = pi): Chi-squared approximation may be incorrect##
## Chi-squared test for given probabilities
##
## data: ni
## X-squared = 41.312, df = 12, p-value = 4.345e-05ここまでの例では、理論分布のかたちがパラメータを含め既知である場合を扱ってきた。観察分布が理論的に従う分布が分かっている場合でも、パラメータが未知の場合がある。次の例は、このような場合にあてはまる。
表4. 20年10兵団において、兵団あたり年あたりの馬に蹴られて死亡した兵士の数
| 0 | 1 | 2 | 3 | \(\ge\) 4 | |
|---|---|---|---|---|---|
| 観察度数(O) | 109 | 65 | 22 | 3 | 1 |
馬に蹴られて死亡する事故は、あまり高い頻度で起こる事象ではなく、その頻度はポアソン分布にしたがうと期待される。ポアソン分布とは、自然数をとる確率変数\(X\)が、 \[ P[X = k] =\frac{\lambda^k \exp(-\lambda)}{k!} \] を満たすときに、確率変数\(X\)が従う分布である。ポアソン分布のパラメータは\(\lambda\)で、ポアソン分布の平均と分散はともに\(\lambda\)に等しいという性質をもつ。
馬に蹴られて死亡した兵士の数（年あたり、兵団あたり）が、ポアソン分布に従うかどうか検定してみる。まず、死亡した兵士数の平均を求めることで、未知であるポアソン分布のパラメータ\(\lambda\)を求める。つぎに、そのパラメータ\(\lambda\)のもとで、死亡した兵士数の期待度数を求める。
期待度数を計算し、度数分布図を描いてみる。なお、4人以上の兵士が死亡したという1例については、4人が死亡したとして平均をもとめる。また、4人以上になる確率については、3人以下になる確率を1から減ずることにより計算する
# 馬に蹴られて死亡した戦士数データ (10兵団×20年)
x <- 0:4
ni <- c(109, 65, 22, 3, 1)
# 兵団×年数
(n <- sum(ni))## [1] 200# 全死亡者数
(d <- sum(x * ni))## [1] 122# 1兵団1年あたりの死亡者数
(lmbd <- d / n)## [1] 0.61# ポアソン分布を仮定した期待頻度を考える
dpois(0:4, lmbd)## [1] 0.543350869 0.331444030 0.101090429 0.020555054 0.003134646# 4以上になる確率は3以下になる確率を1から引く
pi <- c(dpois(0:3, lmbd), 1 - sum(dpois(0:3, lmbd))); pi## [1] 0.543350869 0.331444030 0.101090429 0.020555054 0.003559618pi <- c(dpois(0:3, lmbd), 1 - ppois(3, lmbd)); pi## [1] 0.543350869 0.331444030 0.101090429 0.020555054 0.003559618npi <- sum(ni) * pi
data <- data.frame(ni, npi, row.names = x)
barplot(t(data), beside = T)
図3. 馬に蹴られて死亡した兵士数（年あたり、兵団あたり）
観察頻度と期待頻度は、よく一致しているようにみえる。では、カイ二乗検定を用いて、適合度検定を行ってみる。
# 観察度数(O)と期待度数(E)の差
ni - npi## [1] 0.3298262 -1.2888060 1.7819142 -1.1110108 0.2880765# （O - E)^2 / Eを計算
(v <- (ni - npi)^2 / npi)## [1] 0.00100106 0.02505734 0.15704840 0.30025340 0.11656877# 和をとって、自由度3のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 0.5999291 - pchisq(ch2, length(ni) - 1 - 1)## [1] 0.8964486観察度数の期待度数からのずれは有意でなく、馬に蹴られて死んだ兵士数の分布はポアソン分布に従っていないとはいえない、と結論される。
なお、観察度数と期待確率をもとに、chisq.testを使って計算すると、この例の場合は自由度が異なる（\(\lambda = 0.61\)をデータから計算しているために自由度が1減ってしまう）ために異なる\(p\)値が得られてしまうので注意する。
chisq.test(ni, p = pi)## Warning in chisq.test(ni, p = pi): Chi-squared approximation may be incorrect##
## Chi-squared test for given probabilities
##
## data: ni
## X-squared = 0.59993, df = 4, p-value = 0.9631ここまでは、統計学的に期待頻度が求められる例について扱った。先に述べたように、交配後代の形質分離のように、遺伝学的に期待される理論比や理論分布をもとに、期待度数を求め、観察度数と期待度数の適合度を検定することもできる。
では、有名なメンデルのエンドウの種子の色と形（丸としわ）のデータを用いて適合度検定を行ってみる。
表5. メンデルの実験におけるエンドウの種子の色と形の分離
| 丸・黄 | しわ・黄 | 丸・緑 | しわ・緑 | 計 | |
|---|---|---|---|---|---|
| 観察度数(O) | 315 | 108 | 101 | 32 | 556 |
メンデル遺伝をする形質では、雑種第2代（\(F_2\)）において、優性形質と劣性形質が3:1の分離を示すことが知られている。ここでは2つの形質、種子の色と形を同時に計測しており、それぞれが独立に遺伝しているとする。種子の形では、丸がしわに対して優性、黄色が緑色に対して優性のとき、丸・黄、しわ・黄、丸・緑、しわ・緑が、それぞれ、9:3:3:1の比で分離すると期待される。すなわち、確率で表すと、9/16, 3/16, 3/16, 1/16の確率で分離すると考えられる。
では、Rを用いて期待度数を計算し、度数分布図を描いてみましょう。
# メンデルのデータ
ni <- c(315, 108, 101, 32)
# 期待度数を計算
pi <- c(9, 3, 3, 1) / 16
npi <- sum(ni) * pi
# グラフを描く
data <- data.frame(ni, npi, row.names = c("RY", "WY", "RG", "WG"))
barplot(t(data), beside = T)
図4. メンデルのエンドウの種子のしわと色の分離データ
観察度数と期待度数はよく一致しているように見える。では、適合度検定を行ってみる。
観察度数と期待度数の差は有意ではなく、つまり、メンデルが観察したエンドウの形と色の分離比は、彼が提示した理論比に合わないとはいえない、と結論される。
# 観察度数(O)と期待度数(E)の差
ni - npi## [1] 2.25 3.75 -3.25 -2.75# （O - E)^2 / Eを計算
(v <- (ni - npi)^2 / npi)## [1] 0.01618705 0.13489209 0.10131894 0.21762590# 和をとって、自由度3のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 0.4700241 - pchisq(ch2, length(ni) - 1)## [1] 0.9254259なお、この例の場合はchisq.test関数を用いても同じ検定ができる。
# Rの備え付け関数を用いる
chisq.test(ni, p = pi)##
## Chi-squared test for given probabilities
##
## data: ni
## X-squared = 0.47002, df = 3, p-value = 0.9254分割表の検定
分割表は、通常、以下のようなかたちで記述される。
| 1 | 2 | … | c | 計 | |
|---|---|---|---|---|---|
| 1 | \(n_{11}\) | \(n_{12}\) | \(...\) | \(n_{1c}\) | \(n_{1\cdot}\) |
| 2 | \(n_{21}\) | \(n_{22}\) | \(...\) | \(n_{2c}\) | \(n_{2\cdot}\) |
| … | \(...\) | \(...\) | \(...\) | \(...\) | \(...\) |
| \(r\) | \(n_{r1}\) | \(n_{r2}\) | \(...\) | \(n_{rc}\) | \(n_{r\cdot}\) |
| 計 | \(n_{\cdot1}\) | \(n_{\cdot2}\) | \(...\) | \(n_{\cdot c}\) | \(n\) |
なお、\(r\)は1つめの要因Rのカテゴリの数を表し、 \(c\)は2つめの要因Cのカテゴリの数を表す。\((i,j)\)番目のセルの数値\(n_{ij}\)は、要因Rが\(i\)番目、要因Cが\(j\)番目のカテゴリに属する観察度数を表す。また、\(n_{i\cdot}\)は要因Rが\(i\)番目のカテゴリに属する観察度数、\(n_{\cdot j}\) を\(n\)と記す。このとき、セル\((i,j)\)に属する経験的確率は、\(n_{ij}/n\)と計算される。また、要因Rがカテゴリ\(i\)に属する経験的な周辺確率\(p_i\)は\(n_{i\cdot}/n\)、要因Cがカテゴリ\(j\)に属する経験的な周辺確率\(p_j\)は\(n_{\cdot j}/n\)と計算することができる。
ここで、要因RとCが「独立である」場合は、セル\((i,j)\)の期待度数\(e_{ij}\)は、\(np_i p_j\)と計算することができる。したがって、分割表内の数値を使って以下のように計算できる。 \[ e_{ij} =n \times \frac{n_{i\cdot}}{n} \times \frac{n_{\cdot j}}{n} = \frac{n_{i\cdot}\times n_{\cdot j}}{n} \]
こうして期待度数\(e_{ij}\)が計算されれば、それをもとに、カイ二乗検定によって観察度数\(n_{ij}\)が期待度数\(e_{ij}\)に適合しているかどうかを検定できる。期待度数は要因RとCの独立性を仮定した度数なので、この結果は、要因RとCの独立性を検定していることとなる。 要因RとCが独立のとき、観察度数\(n_{ij}\)と期待度数\(e_{ij}\)から計算される \[ \chi^2=\sum_{i=1}^r \sum_{j=1}^c \frac{(n_{ij}-e_{ij})^2}{e_{ij}} \] は、自由度\((r-1)\times (c-1)\)のカイ二乗分布に従う。仮説検定の仕方は、適合度検定の場合と同じとなる。
では、以下に示すRichardsonのアブラムシへの薬剤散布実験データについて、解析を行ってみる。この実験では、異なる濃度のオレイン酸ナトリウムをアブラムシ（Aphis rumicis L.）に散布し、それにともなうアブラムシの生死が調べられている。このデータをもとに、オレイン酸ナトリウムの散布濃度によって、死亡率に違いがあるかを検定してみる。
表5. Richardsonのアブラムシへのオレイン酸ナトリウムの散布実験データ
| 0.65 (%) | 1.10 (%) | 計 | |
|---|---|---|---|
| 死亡 | 55 | 62 | 117 |
| 生存 | 13 | 3 | 16 |
| 計 | 68 | 65 | 133 |
上述したように、周辺度数\(n_{i\cdot}\)および\(n_{\cdot j}\)を計算して、それをもとに期待度数を計算する。計算ができたら、観察度数との差をとり、カイ二乗値を計算して、検定をする。
# Richardsonのアブラムシ防除データ
n <- matrix(c(55, 13, 62, 3), 2)
# 周辺度数の計算
nx <- apply(n, 1, sum); nx## [1] 117 16ny <- apply(n, 2, sum); ny## [1] 68 65# 期待度数の計算
e <- nx %*% t(ny) / sum(n); e## [,1] [,2]
## [1,] 59.819549 57.180451
## [2,] 8.180451 7.819549# 観察度数(O)と期待度数(E)の差
n - e## [,1] [,2]
## [1,] -4.819549 4.819549
## [2,] 4.819549 -4.819549# （O - E)^2 / Eを計算
v <- (n - e)^2 / e; v## [,1] [,2]
## [1,] 0.388302 0.4062236
## [2,] 2.839458 2.9705104# 和をとって、自由度1のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 6.6044951 - pchisq(ch2, (length(nx) - 1)*(length(ny) - 1))## [1] 0.01017217検定の結果は5%水準で有意であり、オレイン酸ナトリウムの濃度とアブラムシの死亡率は独立ではない、と結論される。言い換えると、オレイン酸ナトリウムの濃度水準間でアブラムシの死亡率に5%水準で有意な差がある、といえる。ただし、この検定には、後述する点に注意する必要がある。
なお、分割表の検定も、chisq.test関数を用いて行うこともできる。行列のかたちの分割表をそのまま入力にする。
# chisq.test関数で検定
chisq.test(n)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: n
## X-squared = 5.3052, df = 1, p-value = 0.02126この結果を少し注意深く見てみよう。なお、いずれにしても5%水準で有意だが、先ほどの結果と比べると、\(p\)値が異なっていることに気づく。これは、結果の表示されている前に表示されているイェーツの補正（Yates’ correction）が行われているためである。分割表内に期待度数が5よりも小さい項目がある場合はカイ二乗分布による近似が悪くなることが知られており、その対処法として、以下のような補正が行われる。 \[ \chi^2=\sum_{i=1}^r \sum_{j=1}^c \frac{(|n_{ij}-e_{ij}|-0.5)^2}{e_{ij}} \]
では、さきほどの結果について手計算でイェーツの補正を行ってから、カイ二乗検定を行ってみる。
# Yatesの補正
(v <- (abs(n - e) - 0.5)^2 / e)## [,1] [,2]
## [1,] 0.3119131 0.3263091
## [2,] 2.2808647 2.3861354# 和をとって、自由度1のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 5.3052221 - pchisq(ch2, (length(nx) - 1)*(length(ny) - 1))## [1] 0.02126159chisq.test関数を用いた場合と同じ\(p\)値が得られる。
次に、カテゴリの数が3以上の分割表について、検定を行ってみる。次の例は、StrandとJessenの農地のデータである。
このデータでは、アイオア州オーデュボン群にける農地の肥沃土と農地の保有状況のデータである。ここで興味ある作業仮説は、所有、貸出、それらの混合、という土地の保有状況によって、肥沃度に違いがあるのではないか、というものである。つまり、ここで興味があるのは、保有状況間で、肥沃度に差があるかどうかということになる。言い換えると、保有状況は、肥沃度と独立であるか、ということになる。
表6. StrandとJessenの農地の肥沃度と保有状況のデータ
| 所有 | 貸出 | 混合 | 計 | |
|---|---|---|---|---|
| 肥沃度I | 36 | 67 | 49 | 152 |
| 肥沃度II | 31 | 60 | 49 | 140 |
| 肥沃度III | 58 | 87 | 80 | 225 |
| 計 | 125 | 214 | 178 | 517 |
では、各セルの期待度数を計算して、観察度数と比較してみる。
# Strand and Jessenの農地データ
n <- matrix(c(36, 31, 58, 67, 60, 87, 49, 49, 80), 3)
# 周辺度数の計算
(nx <- apply(n, 1, sum))## [1] 152 140 225(ny <- apply(n, 2, sum))## [1] 125 214 178# 期待度数の計算
(e <- nx %*% t(ny) / sum(n))## [,1] [,2] [,3]
## [1,] 36.75048 62.91683 52.33269
## [2,] 33.84913 57.94971 48.20116
## [3,] 54.40039 93.13346 77.46615# 観察度数(O)と期待度数(E)の差
n - e## [,1] [,2] [,3]
## [1,] -0.7504836 4.083172 -3.3326886
## [2,] -2.8491296 2.050290 0.7988395
## [3,] 3.5996132 -6.133462 2.5338491観察度数と期待度数の差をみると所有と貸出で少し系統的な違いがあるようにも見える。では、これが統計的に有意なものなのか、検定を行って確認してみる。
# （O - E)^2 / Eを計算
(v <- (n - e)^2 / e)## [,1] [,2] [,3]
## [1,] 0.01532566 0.26498944 0.21223471
## [2,] 0.23981531 0.07254031 0.01323919
## [3,] 0.23818240 0.40392957 0.08287996# 和をとって、自由度4のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 1.5431371 - pchisq(ch2, (length(nx) - 1)*(length(ny) - 1))## [1] 0.8189739観察度数と期待度数には有意な差は無く、したがって、農地の所有状況と農地の肥沃度は独立でないとはいえない、と結論される。言い換えると、農地の所有状況による農地の肥沃度に違いがあるとはいえない、と結論される。
chisq.test関数でも計算してみる。
# chisq.test関数で検定
chisq.test(n)##
## Pearson's Chi-squared test
##
## data: n
## X-squared = 1.5431, df = 4, p-value = 0.819同じ結果が得られる。分割表検定では、期待度数の計算が比較的単純なので、chisq.testを用いて計算することでほとんどの場合事足りると考えられる。
フィッシャーの正確確率検定
カイ二乗検定では、度数が5以下になるセルがある場合は、カイ二乗分布による良い近似が得られない。こうした場合、度数の少ないカテゴリを１つにまとめるなどして近似をよいものにする手続きがとられる場合がある。しかし、2×2の分割表では、カテゴリをこれ以上まとめることができない。先述したイェーツの補正は、こうした場合によく用いられてきた方法であるが、フィッシャーの直接（あるいは正確）確率検定（Fisher’s exact test）を用いるほうが精度が良いことが知られている。カイ二乗検定では、\(n\)の数が小さい場合にも近似が悪くなりますが、こうした場合にもフィッシャーの直接確率検定は有効である。
フィッシャーの直接確率検定の原理を説明するために、次の例について考えてみる。(1)ミルクを先に入れたミルクティーと、(2)紅茶を先に入れたミルクティーについて、それを飲めばどちらか識別できると主張する婦人がいた。この婦人が本当に識別できるかどうかを試すため、(1)、(2)それぞれ4杯ずつ、全8杯について、正しく識別できるか試験を行った。その結果、\(x\)杯について(1)を(1)として識別できた。分割表を以下に示す。
表7. 婦人は本当にミルクティーを識別できたのか？
| (1)と推測 | (2)と推測 | 計 | |
|---|---|---|---|
| (1) | \(x\) | \(4-x\) | 4 |
| (2) | \(4-x\) | \(x\) | 4 |
| 計 | 4 | 4 | 8 |
まず、8杯中4杯を無作為に選んだときに、(1)が\(x\)杯含まれている確率は、 \[ p[X=x]=\frac{\left(\begin{array}{cc} 4\\x\end{array}\right)\left(\begin{array}{cc} 4\\4-x\end{array}\right)}{\left(\begin{array}{cc} 8\\4\end{array}\right)} \] として与えられる。これは、超幾何分布（hypergeometric distribution）とよばれる分布であった。一般的には、\(m\)個の白玉と\(n\)個の赤玉が入った箱から、\(k\)個の玉を無作為に取り出す場合に、その中に含まれている白玉の個数\(X\)が従う分布で、その確率は、 \[ p[X=x]=\frac{\left(\begin{array}{cc} m\\x\end{array}\right)\left(\begin{array}{cc} n\\k-x\end{array}\right)}{\left(\begin{array}{cc} m+n\\k\end{array}\right)} \] で与えられる。
Rではdhyper関数を用いて超幾何分布の確率を計算することができる。xを0杯から4杯まで変化させて、それぞれに対応する確率を計算してみる。
# 超幾何分布
dhyper(0:4, 4, 4, 4) # dhyper(x, m, n, k) 今、m = 4, n = 4, k = 4, xは0〜4## [1] 0.01428571 0.22857143 0.51428571 0.22857143 0.01428571すなわち、でたらめ（無作為）に選んだ場合に、4杯ある(1)を4杯とも選ぶ確率は、0.014であることが分かる。なお、4杯ある(1)から3杯以上を選ぶ確率は\(0.229 + 0.014 = 0.243\)で、比較的高い確率であることが分かる。以上のことから、有意水準を5%とすると、婦人が(1)を4杯とも識別できた場合、すなわち、全8杯について正しく識別できた場合にのみ、婦人に識別能力があると結論することができる。
このように、周辺度数を固定した上で、全てのありうる組合せについて直接確率の計算を行い、それをもとに検定を行うのがフィッシャーの直接確率検定である。
では、先ほどイェーツの補正を行ったアブラムシのデータについてフィッシャーの直接確率検定を用いて検定をして、その結果をカイ二乗検定と比べてみる。
表8. Richardsonのアブラムシへのオレイン酸ナトリウムの散布実験データ
| 0.65 (%) | 1.10 (%) | 計 | |
|---|---|---|---|
| 死亡 | 55 | 62 | 117 |
| 生存 | 13 | 3 | 16 |
| 計 | 68 | 65 | 133 |
# Richardsonのアブラムシ防除データ
n <- matrix(c(55, 13, 62, 3), 2)
# Fisherの直接確率検定
fisher.test(n)##
## Fisher's Exact Test for Count Data
##
## data: n
## p-value = 0.01471
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.03595932 0.80781763
## sample estimates:
## odds ratio
## 0.2069404# カイ二乗検定
chisq.test(n)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: n
## X-squared = 5.3052, df = 1, p-value = 0.02126# カイ二乗検定（イェーツの補正なし）
chisq.test(n, correct = FALSE)##
## Pearson's Chi-squared test
##
## data: n
## X-squared = 6.6045, df = 1, p-value = 0.01017その結果、いずれの検定法においても5%水準で有意だが、フィッシャーの直接確率検定は、カイ二乗検定でイェーツの補正を行った場合とも、行わなかった場合とも異なる\(p\)値を示していることが分かる。先述したとおり、イェーツの補正を行うよりも、フィッシャーの直接確率検定を用いるほうが精度が高いと言われている。一昔前までは、計算に大変時間を要するのがフィッシャーの直接確率検定の弱点と言われてきたが、今日では、その弱点もコンピュータの演算能力の向上により克服されている。
なお、アブラムシの例におけるフィッシャーの直接確率検定は、両側検定になっている。そこでは、どちら向きに確率の低い事象が起こるのかは問うていない。いっぽう、先ほどのミルクティーの例は、片側検定になっている。すなわち、婦人が偶然8杯全部を間違える事象は\(p\)値の計算に含めていない。
このことについて、もう少しだけ詳しくみてみる。まずはミルクティーデータを、fisher.test関数を用いて検定してみる。
# ミルクティーデータ
n <- matrix(c(4, 0, 0, 4), 2)
# Fisherの直接確率検定
fisher.test(n) ##
## Fisher's Exact Test for Count Data
##
## data: n
## p-value = 0.02857
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.339059 Inf
## sample estimates:
## odds ratio
## Inf# セル(1,1)の要素が増える向きにのみ検定
fisher.test(n, alternative = "greater") ##
## Fisher's Exact Test for Count Data
##
## data: n
## p-value = 0.01429
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
## 2.003768 Inf
## sample estimates:
## odds ratio
## Inf\(x=4\)のとき、fisher.test関数をオプション無しで実行すると両側検定が行われ、先ほど計算されたより大きな\(p\)値が示される。これは、8杯全てを間違えるという事象も、8杯全て正解するのと同じ確率をもっており、両側検定を行うと、\(p\)値の計算に含められてしまうためである。
上述した説明をより明確にすべく、最後に、fisher.test関数を使わずに、アブラムシのデータについて、フィッシャーの直接確率検定を行ってみる。以降の説明が分かりやすくなるように、分割表の並び順を変え、もっとも度数の少ないセルを\((1,1)\)にもってきた。
表9. Richardsonのアブラムシへのオレイン酸ナトリウムの散布実験データ
| 1.10 (%) | 0.65 (%) | 計 | |
|---|---|---|---|
| 生存 | 3 | 13 | 16 |
| 死亡 | 62 | 55 | 117 |
| 計 | 65 | 68 | 133 |
アブラムシのデータを、超幾何分布として捉える場合、そのパラメータは、\(m=65,n=68,k=16\)となる。このとき、16頭のアブラムシのうち、\(x\)頭が高濃度で処理されたものだとする。すなわち、セル\((1,1)\)の度数が\(x\)頭だとする。dhyperを用いて、\(x=0,…,16\)の場合について、その確率を計算すると、
# アブラムシデータの超幾何分布
dhyper(0:16, 65, 68, 16)## [1] 8.207353e-06 1.610499e-04 1.431555e-03 7.652313e-03 2.753466e-02
## [6] 7.072060e-02 1.341253e-01 1.916075e-01 2.083732e-01 1.730750e-01
## [11] 1.094281e-01 5.210861e-02 1.831943e-02 4.596118e-03 7.759679e-04
## [16] 7.875495e-05 3.619253e-06ここで、セル(1,1)が3以下になる確率だけを考え、それを足し合わせる。
# このうち、生存個体が3個体以下になる確率は
dhyper(0:3, 65, 68, 16)## [1] 8.207353e-06 1.610499e-04 1.431555e-03 7.652313e-03# 足し合わせると
sum(dhyper(0:3, 65, 68, 16))## [1] 0.009253125次に、fisher.test関数で片側確率を計算する。なお、表9のように並び替えたので、改めてデータを準備して、それを用いて片側検定（この場合、セル(1,1)が減る方向なので、alternative = “less”というオプションを指定）をしてみる。
# 並べ替えた分割表を準備
n <- matrix(c(3, 62, 13, 55), 2)
# fisher.test関数で計算。ただし、3より少なくなる側のp値
fisher.test(n, alternative = "less")##
## Fisher's Exact Test for Count Data
##
## data: n
## p-value = 0.009253
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
## 0.0000 0.6801
## sample estimates:
## odds ratio
## 0.2069404計算される\(p\)値は、先ほどdhyperを用いて計算されたものに一致する。
最後に、両側検定をdhyperを用いて行ってみる。両側検定を行う場合は、\(x = 0,…,16\)のうち、現在観察されている状態（\(x = 3\)）と確率が同じか、それより小さくなる（より起こりにくい事象）ものについて、その確率を足し合わせる。すなわち、
# 両側検定のときは？
# x = 0, ..., 16のときの確率を計算
(pa <- dhyper(0:16, 65, 68, 16))## [1] 8.207353e-06 1.610499e-04 1.431555e-03 7.652313e-03 2.753466e-02
## [6] 7.072060e-02 1.341253e-01 1.916075e-01 2.083732e-01 1.730750e-01
## [11] 1.094281e-01 5.210861e-02 1.831943e-02 4.596118e-03 7.759679e-04
## [16] 7.875495e-05 3.619253e-06# x = 3のときの確率を計算
(pp <- dhyper(3, 65, 68, 16))## [1] 0.007652313# ｘ = 3のときの確率以下になるものを足し合わせる
(p.less <- pa[pa <= pp])## [1] 8.207353e-06 1.610499e-04 1.431555e-03 7.652313e-03 4.596118e-03
## [6] 7.759679e-04 7.875495e-05 3.619253e-06sum(p.less)## [1] 0.01470759この\(p\)値は、fisher.test関数で行う両側検定における\(p\)値と一致する。
# fisher.test関数で計算。両側検定
fisher.test(n)##
## Fisher's Exact Test for Count Data
##
## data: n
## p-value = 0.01471
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.03595932 0.80781763
## sample estimates:
## odds ratio
## 0.2069404参考文献
以下は、本テキストを作成するにあたり参考にした書籍です。より詳しく学ぶためにもお薦めの本です。
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.