カイ二乗検定
削除 （当サイト分割表の解析 参照。）
http://hyoccori-toukei.la.coocan.jp/_/id/id3.html#4_分割表 高校情報
http://hyoccori-toukei.la.coocan.jp/id-2/id-3/id/id1.html 高校数学
ピアソンは、様々な統計量がどのような分布に従うかを研究したと言われています。
（人間環境大学
三中信宏様「統計の落とし穴と蜘蛛の糸」）参照
通常は統計量を天下り的に覚え適用するのが普通で、筆者も何の疑問ももたずにいました。しかし、理解を深めるには「何故、如何にして」（Why & How）を考えてみることも有意義であると思います。
高等学校では、情報Ⅰでクロス集計が扱われ、独立性の検定としてカイ二乗検定が紹介されています。
ここでは、四分表について、シミュレーションを用いて、カイ二乗統計量を考えます。
1 四分表
新薬の効果を検証するには、新薬を投与するグループとブラセポ（偽薬）を投与するグループをランダムに２群に分け、次のような四分表をつくります。
| 効果あり | 効果なし | 計 | ||||
|---|---|---|---|---|---|---|
| 新薬 | a | b | a+b | |||
| ブラセポ | c | d | c+d | |||
| 合計 | a+c | b+d | a+b+c+d | |||
新薬に効果がないなら、
効果の比率は同じで、 \(a:b=c:d~~~~~~~\therefore~~ad-bc=0\)
\(T=ad-bc\) は 誤差変動を考慮しても 0
に近い値をとります。
このとき「薬の種類」と「効果」は独立であると言われます。
新薬に効果があるならば、
新薬の群の方が「効果あり」の比率が高いはずです。
\(T=ad-bc\) は 0
から離れた値をとると予想されます。
薬の種類により、効果に差があるかを調べることを独立性の検定といいます。
独立であれば、薬の種類という情報は、効果の判断に影響を与えません。
四分表において「２群が独立であるかどうか」の統計量としては、
\(T=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)} \sim \chi^2(1)\)
が用いられています。ピアソンは、どのようにしてこの統計量にたどり着いたのでしょうか？
（ ad-bc は、比例式 a:b=c:d と関係がありそうです）
探究的学習として、統計量の発見に至るまでのピアソンの試行錯誤を追体験してみましょう。このような経験は遠回りに見えますが、これから研究活動を進める生徒にとって有意義なことだと考えます。
1.1 独立.
「薬の種類」と「効果」が独立という帰無仮説のもとで、統計量
\(X=ad-bc\) の分布を調べます。
薬の種類によって効果に差がないモデルは、効果の有無の比率が同じモデルです。
調査総数を N とします。実験モデルに指定するパラメータは、
- 新薬を投与される人の割合 \(\pi_1\)
- 共通の効果の割合 \(\pi_2\)
独立という仮説のもとでの四分表の確率分布は次のようになります。.
| 確率分布 | 効果あり | 効果なし | 計 |
|---|---|---|---|
| 新薬 | \(\pi_1\pi_2\) | \(\pi_1(1-\pi_2)\) | \(\pi_1\) |
| ブラセポ | \((1-\pi_1)\pi_2\) | \((1-\pi_1)(1-\pi_2)\) | \(1-\pi_1\) |
| 合計 | \(\pi_2\) | \(1-\pi_2\) | 1 |
データの出力には、多項乱数
rmultinorm(乱数の数, 調査数, 比率) を用います。
＜注意＞.
簡単のため、行の周辺度数を固定して考えます。
各行の中での効果の人数は二項分布に従い、考えやすくなります。
データ数が十分大きければ二項分布は正規分布で近似できるため、探究的学習の可能性が広がります。
N=500、\(\pi_1=3/5\)
とすると、行の比率は \(3:2\)
なので、
１行目の合計は300、２行目の合計は200 となります。
１行目の合計 300 を、\(\pi_2=3/10\) に従って \(3:7\) に分けると、
ab <- rmultinom(n=1, size = 300, prob = c(3/10,7/10))
t(ab)#転置## [,1] [,2]
## [1,] 71 229新薬を300人に投与したとき、(効果のあった人数, なかった人数)
です。
１番目の数字が a 、２番目の数字が b になります。
２行目の合計 200 を、\(\pi_2=3/10\) に従って \(3:7\) に分けると、
cd <- rmultinom(n=1, size = 200, prob = c(3/10,7/10))
t(cd)#転置## [,1] [,2]
## [1,] 61 139ブラセポ（偽薬）を200人に投与したとき、(効果のあった人数,
なかった人数) です。
１番目の数字が c 、２番目の数字が d になります。
四分表の形式で表示するには、matrix()関数を使います。
a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
t(matrix(c(a,b,c,d),2,2))#転置## [,1] [,2]
## [1,] 71 229
## [2,] 61 139これで ad-bc が1つ計算されます。
n=1000 とすれば、1000個の ad-bc
が計算され、ヒストグラムを描くことができます。
ab <- rmultinom(n=1000, size = 300, prob = c(3/10,7/10))
cd <- rmultinom(n=1000, size = 200, prob = c(3/10,7/10))
a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
x <- a*d-b*c
hist(x,prob=T)
正規分布に従っているかどうかをみるため、データをscal()関数で標準化したヒストグラムと標準正規分布の理論曲線を比べます。
hist(scale(x),prob=T)
curve(dnorm(x, 0,1),add=T)
統計量 \(X=ad-bc\) は正規分布に従っているようです。
1.2 ad-bc の分散
独立という仮説のもとで \(X=ad-bc\)
は平均０の正規分布に従うとして、
\(X\)
の分散が分かれば、標準化したものを統計量とすることができます。
ピアソンの統計量は、
\(T=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)} \sim \chi^2(1)\)
\(\chi^2(1)\) は標準正規分布の２乗なので、
\(t=\dfrac{(ad-bc)\sqrt{N}}{\sqrt{(a+b)(c+d)(a+c)(b+d)}} \sim N(0,1)\)
\(X=ad-bc\) の分散は、
\(v=\dfrac{(a+b)(c+d)(a+c)(b+d)}{N}\)
となりそうですが、
このシミュレーションモデルでは \(v\)
は変動する量です。
行の周辺度数 \((a+b)=300,~ (c+d)=200\)
は定数ですが、列の周辺度数は変動するからです。
代替として、列の周辺度数を、 \(500\times
\dfrac{3}{10}=150,~~500\times \dfrac{7}{10}=350\) と考え、
\(V=\dfrac{300\times 200 \times
150\times 350}{500}\) としてみましょう。
m <- c(mean(a+c),mean(b+d)) #列の周辺度数
m## [1] 149.634 350.366v <- c(300*200*150*350/500,mean(300*200*(a+c)*(b+d)/500),var(x))
v## [1] 6300000 6279236 66312471.3 ピアソンの予想
筆者の想像ですが、ピアソンは上記の例で
\(X=ad-bc\) の分散を \(V=\dfrac{300\times 200 \times 150\times
350}{500}\) と考え、
\(z=\dfrac{ad-bc}{\sqrt{V}}\sim
N(0,~1)\) と予想したのではないでしょうか。
実際、
V <- 300*200*150*350/500
z <- x/sqrt(V)
hist(z,prob=T)
curve(dnorm(x, 0,1),add=T)
数理的証明は難しそうですが、標準正規分布に近いようです。
さらにピアソンは、
- 実際の調査では、パラメーター\(\pi_1,\pi_2\) が未知であること、
- \(a+c\fallingdotseq \pi_2N,~ b+d\fallingdotseq (1-\pi_2)N\) が成立することを勘案し、
\(t=\dfrac{(ad-bc)\sqrt{N}}{\sqrt{(a+b)(c+d)(a+c)(b+d)}}\)
を近似的に標準正規分布に従う統計量として利用することを考えたのではないかというのが筆者の予想です。
1.3.1 数理的考察
数学Bでは、二項分布が正規分布で近似できることを扱っているので、この観点から考えてみましょう。
２群の効果が等しいという帰無仮説のもと、全体で効果がある確率を \(p\) とします。
新薬投与 \(n_1\)
人の中で効果がある人が \(a\)
人とすると、\(a\sim B(n_1,~p)\).
\(\hat{p_1}=\dfrac{a}{n_1}\)
とすると、 \(E(a)=n_1p~,~~V(a)=n_1p(1-p)\) だから、
\(E(\hat{p_1})=p\) \(V(\hat{p_1})=\dfrac{p(1-p)}{n_1}\).
\(n_1\)
が大きければ正規近似により、\(\hat{p_1}\sim
N\left(p,~\dfrac{p(1-p)}{n_1}\right)\)
ブラセポ投与 \(n_2\)
人の中で効果がある人が \(c\)
人とすると、\(c\sim B(n_2,~p)\).
\(\hat{p_2}=\dfrac{c}{n_2}\)
とすると同様に、
\(n_2\)
が大きければ正規近似により、\(\hat{p_2}\sim
N\left(p,~\dfrac{p(1-p)}{n_2}\right)\).
標本比率の差 \(w=\hat{p_1}-\hat{p_2}\) は、
\(w \sim
N\left(0,~~\dfrac{p(1-p)}{a+b}+\dfrac{p(1-p)}{c+d}\right)=N\left(0,~~\dfrac{p(1-p)N}{(a+b)(c+d)}\right)\)
\(p=\dfrac{a+c}{N}\) で近似すると、
\(\dfrac{p(1-p)N}{(a+b)(c+d)}=\dfrac{(a+c)}{N}\dfrac{(b+d)}{N}\dfrac{N}{(a+b)(c+d)} =\dfrac{1}{N}\dfrac{(a+c)(b+d)}{(a+b)(c+d)}\).
\(\hat{p_1}=\dfrac{a}{a+b}\:,\hat{p_2}=\dfrac{c}{c+d}\) であるので、その差は、
\(w=\hat{p_1}-\hat{p_2}=\dfrac{a}{a+b}-\dfrac{c}{c+d}=\dfrac{ad-bc}{(a+b)(c+d)}\).
標準化すると、
\(z=\dfrac{(ad-bc)}{(a+b)(c+d)}\sqrt{N\dfrac{(a+b)(c+d)}{(a+c)(b+d)}}=\dfrac{(ad-bc)\sqrt{N}}{\sqrt{(a+b)(c+d)(a+c)(b+d)}}\sim N(0,~1)\).
二乗すると自由度１のカイ二乗統計量が得られ、
\(\chi^2=z^2=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\sim
\chi^2(1)\)
但し、近似が使われており、近似が正当性をもつには、周辺度数が大きいことが必要になります。
まとめ
周辺度数が十分大きいとき、
\(X=ad-bc\) の標準化は \(t\fallingdotseq \dfrac{(ad-bc)\sqrt{N}}{\sqrt{(a+b)(c+d)(a+c)(b+d)}} \sim N(0,1)\)
\(T=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)} \sim\chi^2(1)\)
第二の式は、偏差平方和という乖離度を表す統計量につながっていきます。
1.4 独立検定
対立仮説「新薬に効果がある」のもとでシミュレーションを行い、独立性を検定できるか確かめます。
1.4.1 z検定
帰無仮説のもとでの理論曲線を青で表示しておきます。
N=500
ab <- rmultinom(n=1000, size = 300, prob = c(3/10,7/10))
cd <- rmultinom(n=1000, size = 200, prob = c(2/10,8/10))
a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
x <- a*d-b*c
t <- x*sqrt(N)/sqrt((a+b)*(c+d))/sqrt((a+c)*(b+d))
hist(t,prob=T)
curve(dnorm(x, 0,1),add=T,col="blue")
abline(v = qnorm(0.95,0,1),col="red")
赤の縦線は帰無仮説のもとでの棄却境界であり、ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。つまり２群の比率の差が検出されます。
この例では調査数が大きいので、わずかの比率の差でも検出できるようです。
1.4.2 カイ二乗検定
Z値の二乗であるカイ二乗値を用いて、比率の差を検定します。
N=500
ab <- rmultinom(n=1000, size = 300, prob = c(3/10,7/10))
cd <- rmultinom(n=1000, size = 200, prob = c(2/10,8/10))
a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
x <- a*d-b*c
k <- x^2*N/(a+b)/(c+d)/(a+c)/(b+d)
hist(k,prob=T)
curve(dchisq(x,1),add=T,col="blue")
abline(v = qchisq(0.95,1),col="red")
ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。
z検定でもカイ二乗検定でも、比率の差を検出できない可能性があることも知っておく必要があります。
2 適合度検定
適合度検定とは「分割表の観測度数が理論比率にあっているか」を検証する検定方法です。
検定統計量は \(\displaystyle
\chi_0^2=\sum_{i=1}^{k}\dfrac{(o_i-e_i)^2}{e_i}\)（理論比率からの乖離度を表す統計量）
ピアソンは、どのようにしてこの検定量にたどり着いたのでしょうか？
四分表の独立性を評価する検定量と、理論度数からの乖離度を表す統計量を比較することで調べてみましょう。
2.1 乖離度(適合度)
- 適合度とは、観測度数が理論度数にどの程度一致しているかをいいます。
- 乖離度とは、観測度数と理論度数の差の大きさの程度をいいます。
乖離度が小さければ、理論通りの度数が観測されていることになります。
つまり、適合度を測る統計量が乖離度なのです。
帰無仮説（独立）のもとでの四分表の理論度数は、
| 度数分布 | 効果あり | 効果なし | 計 |
|---|---|---|---|
| 新薬 | \(\pi_1\pi_2N\) | \(\pi_1(1-\pi_2)N\) | \(\pi_1N\) |
| ブラセポ | \((1-\pi_1)\pi_2N\) | \((1-\pi_1)(1-\pi_2)N\) | \((1-\pi_1)N\) |
| 合計 | \(\pi_2N\) | \((1-\pi_2)N\) | N |
観測度数の理論度数からの乖離度を計算しましょう。
四分表は、1つのセルの値が増えれば、隣のセルは同じだけ減じます。
構造上、各セルの理論値からの乖離度は全て同じはずです。
\(\pi_1=\dfrac{a+b}{N},~~\pi_2=\dfrac{a+c}{N}\) で近似すると、
\(a-N\pi_1\pi_2=a-\dfrac{(a+b)(a+c)}{N}=\dfrac{a(a+b+c+d)-(a+b)(a+c))}{N}=\dfrac{ad-bc}{N}\)
同様に、
\(b-N\pi_1(1-\pi_2)=c-N(1-\pi_1)\pi_2=-\dfrac{ad-bc}{N}\)
\(d-N(1-\pi_1)(1-\pi_2)=\dfrac{ad-bc}{N}\) となることが確かめられます。
＜考察＞
- 独立ならば、観測度数はこの理論度数にほぼ一致しており、乖離度は０に近いはずです。
- 独立ならば、カイ二乗統計量も０に近い値をとります。
各セルの乖離度とカイ二乗統計量の関係を調べましょう。
統計学では平方和を統計量とすることがよくあります。
乖離度の平方和から統計量を作ることを考えます。
- 各セルの乖離度の平方 \(\left(\dfrac{ad-bc}{N}\right)^2\)
- \(\chi^2=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\)
両者の関係に注目しましょう。
\(\chi^2=\left(\dfrac{ad-bc}{N}\right)^2~\dfrac{N^3}{(a+b)(c+d)(a+c)(b+d)}=\left(\dfrac{ad-bc}{N}\right)^2~\dfrac{N^4}{(a+b)(c+d)(a+c)(b+d)}\dfrac{1}{N}\)
\(\pi_1=\dfrac{a+b}{N},~~\pi_2=\dfrac{a+c}{N},~~1-\pi_1=\dfrac{c+d}{N},~~1-\pi_2=\dfrac{b+d}{N}\) であるので、
\(\chi^2=\left(\dfrac{ad-bc}{N}\right)^2 \dfrac{1}{\pi_1\pi_2 (1-\pi_1)(1-\pi_2)}\dfrac{1}{N}\)
ところで、
\(\dfrac{1}{\pi_1(1-\pi_1)\pi_2
(1-\pi_2)}=\left(\dfrac{1}{\pi_1}+\dfrac{1}{1-\pi_1}\right)\left(\dfrac{1}{\pi_2}+\dfrac{1}{1-\pi_2}\right)\) を展開すると、
\(\dfrac{1}{\pi_1\pi_2 (1-\pi_1)(1-\pi_2)}=\dfrac{1}{\pi_1\pi_2}+\dfrac{1}{\pi_1(1-\pi_2)}+\dfrac{1}{(1-\pi_1)\pi_2 }+\dfrac{1}{(1-\pi_1)(1-\pi_2)}\)
\(\dfrac{1}{\pi_1\pi_2 (1-\pi_1)(1-\pi_2)N}=\dfrac{1}{\pi_1\pi_2N}+\dfrac{1}{\pi_1(1-\pi_2)N}+\dfrac{1}{(1-\pi_1)\pi_2N }+\dfrac{1}{(1-\pi_1)(1-\pi_2)N}\)
従って、
\(\chi^2=\dfrac{\left(\frac{ad-bc}{N}\right)^2}{N\pi_1\pi_2}+\dfrac{\left(\frac{ad-bc}{N}\right)^2}{N\pi_1(1-\pi_2)}+\dfrac{\left(\frac{ad-bc}{N}\right)^2}{N(1-\pi_1)\pi_2}+\dfrac{\left(\frac{ad-bc}{N}\right)^2}{N(1-\pi_1)(1-\pi_2)}\)
この４つの項の分母は四分表の理論度数そのものであり、
カイ二乗値は、各セルの乖離度の平方を理論度数で除した値の和
に等しいことがわかります。
独立性のカイ二乗統計量は、乖離度の平方和で求めることが可能になったわけです。
ピアソンは、この乖離度の平方和を一般の分割表のカイ二乗統計量として定義し、独立性の検定や分布の適合度検定に使おうとしたのではないかと筆者は考えます。
四分表の計算は、高校生にも確かめることは可能であり、探究的学習の良い題材となるでしょう。
実際に調査をしてクロス集計をする前に、
- どのくらいの調査数があれば理論曲線に近くなるか
- どのくらいの比の差があれば検出できるか
シミュレーションで事前に確かめておくことが必要です。
2.2 ピアソンの定理
データをある区分に従って集計し次のような区分表を得たとしましょう。
| 区分 | \(A_1\) | \(A_2\) | \(\cdots\) | \(A_k\) | 合計 |
|---|---|---|---|---|---|
| 観測度数 | \(o_1\) | \(o_2\) | \(\cdots\) | \(o_k\) | N |
| 理論度数 | \(e_1\) | \(e_2\) | \(\cdots\) | \(e_k\) | N |
四分表のカイ二乗値の一般化と考えれば、区分表のカイ二乗値は、
\(\displaystyle
\chi_0^2=\sum_{i=1}^{k}\dfrac{(o_i-e_i)^2}{e_i}\) と考えることができます。
（カイ二乗値は、各セルの乖離度の平方を理論度数で除した値の和）
＜ピアソンの定理＞
「帰無仮説：観測度数は理論度数の分布から得られた」が真であるとき、
即ち、観測度数が理論比率に従っているとき、
\(\chi_0^2\) は漸近的に自由度 𝑘 −
𝑟 − 1のカイ二乗分布に従います。
（𝑟は期待度数を求めるために観測値から推定した未知母数の数）
四分表
四分表は適合度検定の k=4 のときの特殊な場合です。
「2.1 乖離度」の表では、
\(\pi_1=\dfrac{a+b}{N},~~\pi_2=\dfrac{a+c}{N}\)
によって2つの母数を推定しているので、
自由度は \(4-2-1=1\) となることに注意しましょう。
高校でピアソンの定理の証明を扱うことはありません。
しかし、実験も証明もなしで使うことには抵抗があります。
そこで探究学習の一環として、
ピアソンの定理を、乱数シミュレーションで体験的に理解する.
という試みを行ってみましょう。
2.3 二項分布.
二項分布 \(X \sim B(10,0.5)\)に従う二項乱数を100個出力し、データと理論分布を比較します。
a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(-0.5, 10.5, by = 1),plot = FALSE)
par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(mfrow=c(1,2))
barplot(h$counts,names.arg=h$mids,ylim=c(0,30),main="データ分布")
x=0:10
m=dbinom(x,10,0.5)*100
barplot(m,names.arg=c(0:10),ylim=c(0,30),main="理論分布")
適合度検定量 \(\chi_0^2\) を計算します。
x <- h$counts
K <- sum((x-m)^2/m)
K## [1] 10.18971ピアソンの定理によれば、
\(N=100\)、各区分の確率を \(p_0,~p_1,~\cdots p_{10}\) とすると、
n区分の適合度検定量は自由度(n-1)のカイ二乗分布に従います。即ち、
\(T=\displaystyle \sum_{i=0}^{10}\dfrac{(X_i-Np_i)^2}{Np_i}\sim \chi^2(10)\)
これをシミュレーションで確かめましょう。
削除.
2.3.1 乖離度の計算.
乱数による出現分布と理論確率による分布との乖離度をカイ二乗値によって計算する.
削除.
a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(-0.5, 10.5, by = 1),plot=FALSE)
x <- h$counts
a=0:10
m <- dbinom(a,10,0.5)*100
par(mfrow=c(1,2))
plot(x,type="l")
plot(m,type="l")
K <- sum((x-m)^2/m)
K## [1] 9.7224132.3.2 シミュレーション.
1000回繰り返して、自由度１０のカイ二乗分布の理論曲線と比較します。
K <- c()
for(i in 1:1000){
a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(-0.5, 10.5, by = 1), plot = FALSE)
x <- h$counts
m <- dbinom(0:10,10,0.5)*100
k <- sum((x-m)^2/m)
K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,breaks="scott")
curve(dchisq(x,10),xlim=c(0,20),add=T)
mean(K)## [1] 10.40606mean(K) の値から見て、自由度１０のカイ二乗分布に従いますが、「理論分布と少しずれがあるように感じられます。
ピアソンの定理は、\(Np_i\geq 5\) （𝑖
= 1,2, ⋯ , 𝑘）
であれば、カイ二乗分布による近似がよいと言われています。
度数５未満のセルが多くあり、近似の精度に影響していると考えられます。
そこで、コクランの規則に従って補正を行います。
2.3.3 コクランの規則.
コクランの規則（Cochran’s rule）とは、カイ二乗検定を適用するには、期待値が 5 未満のセルが全体の20%以上にならないようにしなければならないというものです。
［コクランの規則］Cochran’s rule.
- 期待値が 5 未満のセルが全体の20％以上にならないようにする.
具体的には、
期待度数５以下の区分が
20％を超えないよう左端の区間と右端の区間を合併します.
コクランの規則に従ってセルを合併すると近似の精度が上がります。
2.3.4 補正シミュレーション.
右端の３区分と、左端の３区分を合併し、全部で７区分とすれば、適合度検定の自由度は 7-1=6 となるはずです。
K <- c()
for(i in 1:1000){
a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(-0.5, 10.5, by = 1), plot = FALSE)
x <- c(sum(h$ counts[1:3]),h$ counts[4:8],sum(h$ counts[9:11]))
m <- dbinom(0:10,10,0.5)
n <- c(sum(m[1:3]),m[4:8],sum(m[9:11]))*100
k <- sum((x-n)^2/n)
K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,breaks="scott")
curve(dchisq(x,6),xlim=c(0,20),add=T)
mean(K)## [1] 5.930898かなり近似が良くなっており、自由度は６です。
2.3.5 推定母数があるとき.
期待度数を求めるために観測値からr個の母数を推定していれば、自由度は r 減ります。
\(T=\displaystyle \sum_{i=1}^{k}\dfrac{(X_i-Np_i)^2}{Np_i}\rightarrow \chi^2(k-r-1)\).
p=0.5 の代わりに、推定値 p = mean(a)/10 を用います
理論累積確率 dbinom(h$breaks, 10, 0.5) のところを.
推定値を用いて dbinom(h$breaks, 10, mean(a)/10).
としてシミュレーションすると自由度は１減り、自由度５となるはずです。
K <- c()
for(i in 1:1000){
a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(-0.5, 10.5, by = 1), plot = FALSE)
x <- c(sum(h$ counts[1:3]),h$ counts[4:8],sum(h$ counts[9:11]))
m <- dbinom(0:10,10,mean(a)/10)
n <- c(sum(m[1:3]),m[4:8],sum(m[9:11]))*100
k <- sum((x-n)^2/n)
K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T)
curve(dchisq(x,5),xlim=c(0,20),add=T)
mean(K)## [1] 4.876545統計量の平均は５となっており、自由度が５となっていることが分かります。
自由度５のカイ二乗分布の理論曲線はヒストグラムによく当てはまっています。
以下は、削除.
2.4 四分表.
四分表のカイ二乗統計量については、そのモデルによって自由度が異なるので注意が必要です。
- \(\chi^2=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\) について、
\(\pi_1=\dfrac{a+b}{N},~~\pi_2=\dfrac{a+c}{N}\) という推定値を用いており、自由度は 4-1-2=1.
推定値を用いずに４つのセルの乖離度を合計した場合、
- \(T=\displaystyle \sum_{i=1}^{4}\dfrac{(X_i-Np_i)^2}{Np_i}\sim \chi^2(3)\) について、
４つの乖離度の平方和として求めるので、自由度は 4-1=3.
カイ二乗分布 \(\chi^2(n)\) の期待値は n であるので、 乱数シミュレーションでカイ二乗統計量の平均を求めることで自由度を確かめてみましょう。
2.4.1 多項乱数で生成.
当サイト 計数値データの検定 参照.
市原清志先生の「バイオサイエンスの統計学」p125の例題５−４を題材とします.
白血病の患者140名をランダムに２群に分け、A,B２種の方法で治療するとき、治療後の結果を比較します。 結果をクロス集計し次のような表が得られたとします。
| 生存 | 死亡 | 計 | |
|---|---|---|---|
| 治療 A | 48 | 32 | 80 |
| 治療 B | 24 | 36 | 60 |
| 合計 | 72 | 68 | 140 |
N=140 は固定されていますが、周辺度数は調査ごとに変動します。。
目的は「治療法A,Bに差があるか」を検定することです。
治療Aの割合を p1=72/140 、生存確率を p2=80/140 とし、多項乱数 rmultinorm( ) を用いて 1000回分のクロス集計データを生成します。
p1=72/140;p2=80/140
#c(p1*p2,p2*(1-p1),p1*(1-p2),(1-p1)*(1-p2))
d <- rmultinom(n=1000, size = 140, prob = c(p1*p2,p2*(1-p1),p1*(1-p2),(1-p1)*(1-p2)))
# 乱数を発生
data <- as.data.frame(t(d))
colnames(data) <- c("a", "b", "c", "d")
#data[1,]
a=data[,1]
b=data[,2]
c=data[,3]
d=data[,4]
kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
mean(kai)## [1] 0.958564X = c()
for (i in 1:1000) {
a=data[i,1]
b=data[i,2]
c=data[i,3]
d=data[i,4]
kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
X <-c(X,kai)
}
mean(X)## [1] 0.958564同じデータでも、次のように母比率を用いると自由度は３になる。
これは、
検定量 T を用いた検定は、
- 母比率からの乖離度ではなく.
- 独立性からの乖離度を表している.
つまり、独立かどうかだけを検定していることになる。
p1=72/140;p2=80/140
d <- rmultinom(n=1000, size = 140, prob = c(p1*p2,p2*(1-p1),p1*(1-p2),(1-p1)*(1-p2)))
data <- as.data.frame(t(d))
colnames(data) <- c("a", "b", "c", "d")
X = c()
for (i in 1:1000) {
a=data[i,1]
b=data[i,2]
c=data[i,3]
d=data[i,4]
x <- c(a,b,c,d)
#print(x)
y <- c(80*c(72/140,68/140),60*c(72/140,68/140))
s <- sum((abs(x-y))^2/y)
X <-c(X,s)
}
mean(X)## [1] 2.958621hist(X,prob=T)
curve(dchisq(x,3),add=T)
#abline(v = s1)
#1-pchisq(s1, 3)自由度が３になっている.
しかも、乖離度は有意ではない。
四分表のカイ二乗値についてシミュレーションする。
詳しくは、当サイト 番外編
分割表の解析 を参照のこと。.
| 肺がんあり | 肺がんなし | 計 | |
|---|---|---|---|
| A:喫煙あり | a | b | a+b |
| B:喫煙なし | c | d | c+d |
| 合計 | a+c | b+d | a+b+c+d |
四分表の検定は、
\(T=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\)
が自由度１のカイ二乗分布 \(\chi^2(1)\)に従うことに基づいて行われる。.
この検定量は、観測値から周辺分布の確率を、
\(\pi_1=\dfrac{a+b}{N}\)
\(\pi_2=\dfrac{a+c}{N}\). によって推定し、
独立という仮説のもとで
理論度数と実現値との乖離度を計算したものである。
＜独立という仮説のもとでの確率分布＞.
| 確率分布 | 肺がんあり | 肺がんなし | 計 |
|---|---|---|---|
| A：喫煙あり | \(\pi_1\pi_2\) | \(\pi_1(1-\pi_2)\) | \(\pi_1\) |
| B：喫煙なし | \((1-\pi_1)\pi_2\) | \((1-\pi_1)(1-\pi_2)\) | \(1-\pi_1\) |
| 合計 | \(\pi_2\) | \(1-\pi_2\) | 1 |
この統計量 T は、「四分表の要因が独立である」 という状態からの乖離度を表している。 このことについて、乱数シミュレーションで実験してみよう。
2.5 例５−４.
市原清志先生の「バイオサイエンスの統計学」p125の例題５−４を参照.
白血病の患者140名をランダムに２群に分け、A,B２種の方法で治療するとき、治療後の結果を比較する。 これまでの説明に合わせ、次のような表にする。
| 生存 | 死亡 | 計 | |
|---|---|---|---|
| 治療 A | 48 | 32 | 80 |
| 治療 B | 24 | 36 | 60 |
| 合計 | 72 | 68 | 140 |
目的は「治療法A,Bに差があるか」を検定することである。
検定は、
- 帰無仮説「差がない」のもとで、調査を無限回行い統計量の分布を描いたとき、
- 上記の表から得られる統計量が、統計的に有意な位置にあるかどうかをみる。
無限回の乱数をどのように出力すれば良いだろうか？
\(N=140\) は固定しておこう。
2.5.1 二項乱数を２種類使用（１）.
周辺分布（合計の分布）から、生存確率を \(p_1=\dfrac{72}{140}\) とし、.
治療A 80人、治療B 60人 を固定しておくことにすると、
a1 <- rbinom(n=1, size = 80, prob = c(72/140))
a2 <- rbinom(n=1, size = 60, prob = c(72/140))
x <- c(a1,80-a1,a2,60-a2)
x## [1] 35 45 30 30期待度数と乖離度.
#期待度数
y <- c(80*c(72/140,68/140),60*c(72/140,68/140))
round(y,0)## [1] 41 39 31 29abs(x-y)^2/y## [1] 0.91716270 0.97111345 0.02380952 0.02521008#乖離度
sum(abs(x-y)^2/y)## [1] 1.937296これらは手計算で確かめることができる。
問題に与えられた調査データ (48,32,24,68) の、期待度数からの乖離度は.
x <- c(48,32,24,36)
s1 <- sum(abs(x-y)^2/y)
s1## [1] 5.490196シミュレーション.
1000回出力してみよう。
X = c()
for (i in 1:1000) {
a1 <- rbinom(n=1, size = 80, prob = c(72/140))
a2 <- rbinom(n=1, size = 60, prob = c(72/140))
x <- c(a1,80-a1,a2,60-a2)
#print(x)
y <- c(80*c(72/140,68/140),60*c(72/140,68/140))
s <- sum((abs(x-y))^2/y)
X <-c(X,s)
}
mean(X)## [1] 1.917144hist(X,prob=T)
curve(dchisq(x,2),add=T)
abline(v = s1)
#1-pchisq(s1, 2)自由度が２になっている理由は、
- ２種類の二項乱数からデータを出力している.
- 各二項乱数は、2つの区分を生成するので自由度は１.
- 検定量は各二項乱数の乖離度の和であるので、
- カイ二乗分布の再生性から、自由度は 1+1=2.
2.5.2 二項乱数を２種類使用（２）.
治療A 80人、治療B 60人 を固定しておくが、
生存確率を \(A:p_1=\dfrac{24}{60}~~~B:p_2=\dfrac{48}{80}\)
のように異なる値に設定しても、得られる統計量の分布は \(\chi^2(2)\)である。.
a1 <- rbinom(n=1, size = 80, prob = c(24/60))
a2 <- rbinom(n=1, size = 60, prob = c(48/80))
x <- c(a1,80-a1,a2,60-a2)
x## [1] 29 51 34 26期待度数と乖離度.
#期待度数
y <- c(80*c(24/60,36/60),60*c(48/80,32/80))
round(y,0)## [1] 32 48 36 24abs(x-y)^2/y## [1] 0.2812500 0.1875000 0.1111111 0.1666667#乖離度
sum(abs(x-y)^2/y)## [1] 0.7465278これらは手計算で確かめることができる。
調査データ (48,32,24,68) の乖離度は.
x <- c(48,32,24,36)
s1 <- sum(abs(x-y)^2/y)
s1## [1] 23.33333理論確率を変えたので、与えられたデータの乖離度は大きくなっている。
シミュレーション.
1000回出力してみよう。
X = c()
for (i in 1:1000) {
a1 <- rbinom(n=1, size = 80, prob = c(24/60))
a2 <- rbinom(n=1, size = 60, prob = c(48/80))
x <- c(a1,80-a1,a2,60-a2)
#print(x)
y <- c(80*c(24/60,36/60),60*c(48/80,32/80))
s <- sum((abs(x-y))^2/y)
X <-c(X,s)
}
mean(X)## [1] 2.036059hist(X,prob=T)
curve(dchisq(x,2),add=T)
abline(v = s1)
1-pchisq(s1, 2)## [1] 8.574939e-06２組の２分表のカイ二乗値の和なので、カイ二乗分布の再生性により自由度２のカイ二乗分布になる。
2.5.3 生存確率に推定値を用いる.
生存確率に推定値を用いれば自由度は１になるはずである。
x <- c(48,32,24,36)
a1=48;a2=24
y <- c(80*c((a1+a2)/140,1-(a1+a2)/140),60*c((a1+a2)/140,1-(a1+a2)/140))
s1 <- sum(abs(x-y)^2/y)
s1## [1] 5.490196X = c()
for (i in 1:1000) {
a1 <- rbinom(n=1, size = 80, prob = c(72/140))
a2 <- rbinom(n=1, size = 60, prob = c(72/140))
x <- c(a1,80-a1,a2,60-a2)
#print(x)
y <- c(80*c((a1+a2)/140,1-(a1+a2)/140),60*c((a1+a2)/140,1-(a1+a2)/140))
s <- sum((abs(x-y))^2/y)
X <-c(X,s)
}
mean(X)## [1] 1.037hist(X,prob=T)
curve(dchisq(x,1),add=T)
abline(v = s1)
#1-pchisq(s1, 2)確かに、自由度は１になっている。
#削除
X = c()
for (i in 1:10) {
a1 <- rbinom(n=1, size = 80, prob = c(72/140))
#a2 <- rbinom(n=1, size = 60, prob = c(72/140))
#x <- c(a1,80-a1,a2,60-a2)
x <- c(a1,80-a1)
#print(x)
#y <- c(80*c((a1+a2)/140,1-(a1+a2)/140),60*c((a1+a2)/140,1-(a1+a2)/140))
y <- c(80*c((a1)/80,1-(a1)/80))
s <- sum((abs(x-y))^2/y)
X <-c(X,s)
}
X## [1] 0 0 0 0 0 0 0 0 0 0mean(X)## [1] 0hist(X,prob=T)
curve(dchisq(x,1),add=T)
#abline(v = s1)
#1-pchisq(s1, 2)2.6 検定量 T を用いる.
四分表の検定量 \(T=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\) を用いると、自由度は１になる。.
注目すべきことは、
- 二項乱数を２種類用いて生成したデータでも.
- 多項乱数を用いて生成したデータでも.
自由度が１になっていることを確かめることができる。
x <- c(48,32,24,36)
a=x[1];b=x[2];c=x[3];d=x[4]
s1 <- (a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
s1## [1] 5.4901962.6.1 二項乱数２種類で生成.
X = c()
for (i in 1:1000) {
a1 <- rbinom(n=1, size = 80, prob = c(72/140))
a2 <- rbinom(n=1, size = 60, prob = c(72/140))
x <- c(a1,80-a1,a2,60-a2)
a=x[1];b=x[2];c=x[3];d=x[4]
kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
X <-c(X,kai)
}
mean(X)## [1] 1.052966hist(X,prob=T)
curve(dchisq(x,1),add=T)
abline(v = s1)
1-pchisq(s1, 1)## [1] 0.01912343 感想.
推定した母数を用いて計算すると自由度が１減じるのは、不偏分散の計算で母平均を標本平均で推定すると自由度が１減じるのと似ている. このため不偏分散の計算では、nで割るところを n-1 で割るのであった。
帰無仮説のもとでは、二項乱数でも多項乱数でも簡易計算値は自由度１のカイ二乗分布に従っていることは驚きであり、びっくり ピアソン である。
K = c()
X = c()
for (i in 1:1000) {
a1 <- rbinom(n=1, size = 80, prob = c(72/140))
a2 <- rbinom(n=1, size = 60, prob = c(72/140))
x <- c(a1,80-a1,a2,60-a2)
a=x[1];b=x[2];c=x[3];d=x[4]
kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
K <-c(K,kai)
y <- c(80*c(72/140,68/140),60*c(72/140,68/140))
s <- sum((abs(x-y))^2/y)
X <-c(X,s)
}
mean(K)## [1] 1.010992mean(X)## [1] 1.936924X = c()
for (i in 1:1000) {
a1 <- rbinom(n=1, size = 80, prob = c(72/140))
a2 <- rbinom(n=1, size = 60, prob = c(72/140))
x <- c(a1,80-a1)
#print(x)
y <- c(80*c(72/140,68/140))
s <- sum((abs(x-y))^2/y)
X <-c(X,s)
}
mean(X)## [1] 0.9826726