独立性の検定をやってみる

性別と進学予定ありなしの架空クロス表．

ex<-matrix(c(60,30,50,60),2,2)
rownames(ex)<-c("male","female")
colnames(ex)<-c("yes","no")

ex_margin<-addmargins(ex)
ex_margin
##        yes  no Sum
## male    60  50 110
## female  30  60  90
## Sum     90 110 200
prop.table(ex,margin = 1)
##              yes        no
## male   0.5454545 0.4545455
## female 0.3333333 0.6666667
barplot(prop.table(t(ex),margin = 2),
        col=c("lightgreen","darkgray"),horiz = TRUE, legend.text = TRUE,
        args.legend = list(x ='bottomleft', bty='n'))

期待度数

ex_expected<-matrix(c(ex_margin[3,1]*ex_margin[1,3]/ex_margin[3,3],
  ex_margin[3,1]*ex_margin[2,3]/ex_margin[3,3],
  ex_margin[3,2]*ex_margin[1,3]/ex_margin[3,3],
  ex_margin[3,2]*ex_margin[2,3]/ex_margin[3,3]),2,2)
ex_expected
##      [,1] [,2]
## [1,] 49.5 60.5
## [2,] 40.5 49.5

カイ二乗値を計算する．

ex_chisq<-(ex_expected-ex)^2/ex_expected
ex_chisq
##             yes       no
## male   2.227273 1.822314
## female 2.722222 2.227273
sum(ex_chisq)
## [1] 8.999082

5%, 1%水準棄却域を算出し，これによって帰無仮説「クロス表は独立」を棄却するかどうかを判断する．

qchisq(0.95,df=1)
## [1] 3.841459
qchisq(0.99,df=1)
## [1] 6.634897

p値は以下の通り．

1-pchisq(sum(ex_chisq),df=1)
## [1] 0.002701153

ちなみに，カイ二乗値はRのbase関数chisq.testで算出できる（デフォルトでは2x2クロス表の時，「イェーツの補正」を勝手にするので，correct = FALSEでオフにしておく）．

chisq.test(ex,correct = FALSE)
## 
##  Pearson's Chi-squared test
## 
## data:  ex
## X-squared = 8.9991, df = 1, p-value = 0.002701

独立性の検定の理論的基礎（一般）

多項分布による説明が一般的であるが，ここでは直感的理解を優先してポワソン分布による説明を紹介する．

\(I \times J\)クロス集計表を仮定する．

\[ \begin{array}{|cccc|c|} \hline X_{11} & X_{12} & \cdots & X_{1J} & n_{1*}\\ X_{21} & X_{22} & \cdots & X_{2J} & n_{2*}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ X_{I1} & X_{I2} & \cdots & X_{IJ} & n_{I*}\\ \hline n_{*1} & n_{*2} & \cdots & n_{*J} & N \\ \hline \end{array} \]

それぞれの\(X_{ij}\)が，それぞれ平均・分散\(\mu_{ij}\)のポワソン分布にしたがうとする． \[X_{ij}\sim Pois(\mu_{ij})\] サンプルサイズが十分に大きいとき，中心極限定理より\(Z\)値は標準正規分布に近似する．つまり， \[Z_{ij}=\frac{X_{ij}-\mu_{ij}}{\sqrt{\mu_{ij}}}\approx N(0,1)\] すると，\(Z\)の2乗は，自由度1のカイ2乗分布に近似的にしたがう．つまり， \[Z^2_{ij}=\frac{(X_{ij}-\mu_{ij})^2}{\mu_{ij}}\approx \chi^2(1)\] 結局，\(Z^2\)を各セルで計算して総和した\(\chi^2\)は自由度\(IJ\)のカイ2乗分布にしたがう．つまり， \[\chi^2=\sum_{i=1}^I \sum_{j=1}^J Z^2_{ij}=\sum_{i=1}^I \sum_{j=1}^J \frac{(X_{ij}-\mu_{ij})^2}{\mu_{ij}}\approx \chi^2(IJ)\]

さて，じっさいはポワソン分布のパラメータは未知である．そこで，データから推定する．

ここで，帰無仮説として，周辺分布を固定したとき，行と列が独立であるとする．つまり， \[\hat{\mu}_{ij}=\frac{n_{i*}\times n_{*j}}{N}\] としたときに， \[H_0: \forall i, j,\ \ \mu_{ij}=\hat{\mu}_{ij}\] と仮定する．この帰無仮説のもとで，\(\chi^2\)を計算する．

ただし，データから帰無仮説ものとでのパラメータを推定するために，制約式の数だけ自由度が失われる．

\[ \begin{array}{lcc} \hline \mathrm{constraint} & \mathrm{n.\ of\ constraints} & \mathrm{ex.\ of\ unmoved\ cell} \\ \hline N=\sum_{i}\sum_{j}\hat{\mu}_{ij} & 1 & \hat{\mu}_{IJ} \\ n_{i*}=\sum_{j}\hat{\mu}_{ij}\ \mathrm{for}\ i=1,\ldots,I-1 & I-1 & \hat{\mu}_{1J},\ldots , \hat{\mu}_{(I-1)J} \\ n_{*j}=\sum_{i}\hat{\mu}_{ij}\ \mathrm{for}\ j=1,\ldots,J-1 & J-1 & \hat{\mu}_{I1},\ldots , \hat{\mu}_{I(J-1)} \\ \hline \end{array} \]

\[ \begin{array}{ccc:c} \hat{\mu}_{11} & \cdots & \hat{\mu}_{1(J-1)} & \hat{\mu}_{1J} \\ \vdots & \ddots & \vdots & \vdots \\ \hat{\mu}_{(I-1)1} & \cdots & \hat{\mu}_{(I-1)(J-1)} & \hat{\mu}_{(I-1)J} \\ \hdashline \hat{\mu}_{I1} & \cdots & \hat{\mu}_{I(J-1)} & \hat{\mu}_{IJ} \end{array} \]

結局，\(IJ\)から制約式の合計数を引くと，自由度は \[d.f.=IJ-1-(I-1)-(J-1)=(I-1)(J-1)\] となる．

ゆえに，\(\chi^2\)は，帰無仮説のもとで，自由度\((I-1)(J-1)\)のカイ2乗分布にしたがう．つまり， \[\chi^2=\sum_{i=1}^I \sum_{j=1}^J \frac{(X_{ij}-\hat{\mu}_{ij})^2}{\mu_{ij}}\approx \chi^2((I-1)(J-1))\]

\(2\times 2\)クロス集計表に特化した説明

Ref. Agresti & Finlay, SMSS, Ch. 7 and 8.

各セルの値．

\[ \begin{array}{cc|c} a & b & n_{1*}\\ c & d & n_{2*} \\ \hline n_{*1} & n_{*2} & N \end{array} \]

背後の確率．

\[ \begin{array}{cc|c} \ \ \ \pi_1\ \ \ & 1-\pi_1 & 1\\ \ \ \ \pi_2\ \ \ & 1-\pi_2 & 1\\ \hline \ \ \ \pi\ \ \ & 1-\pi & 1 \end{array} \]

帰無仮説 \[H_0: \pi_1 = \pi_2\]

推定量．

\[\begin{align} \hat{\pi}_1 &= \frac{a}{a+b}=\frac{a}{n_{1*}} \\ \hat{\pi}_2 &= \frac{c}{c+d}=\frac{c}{n_{2*}} \\ \hat{\pi} &= \frac{a+c}{a+b+c+d} \end{align}\]

\(Z\)値． \[Z=\frac{\hat{\pi}_1-\hat{\pi}_2}{\sqrt{\hat{\pi}(1-\hat{\pi})\left(\frac{1}{n_{1*}}+\frac{1}{n_{2*}}\right)}}\] \(Z^2\)値． \[Z^2=\frac{(\hat{\pi}_1-\hat{\pi}_2)^2}{\hat{\pi}(1-\hat{\pi})\left(\frac{1}{n_{1*}}+\frac{1}{n_{2*}}\right)}\approx \chi^2(1)\] ちなみに，\(Z^2\)を変形すると， \[Z^2=\frac{N(ad-bc)^2}{(a+b)(a+c)(b+d)(c+d)}=\chi^2\] ゆえに，\(Z^2\)は\(2\times 2\)のカイ2乗値に等しい．\(Z^2\)は自由度1のカイ2乗分布にしたがうので，\(2\times 2\)のカイ2乗値も自由度1のカイ2乗分布にしたがう．