第5章 クロス集計表
2変数の関係を確認する方法として，クロス集計表と独立性の検定であるχ2検定について学習する。
データの読み込みと確認
chap5 <- read.csv("chap5.csv")table()を用いてシンプルなクロス表を作成する。
attach(chap5) # データフレーム名を無視
table(sex,regu) # クロス表を作成
regu
sex 0 1
0 10 4
1 3 8addmargins()を用いて，クロス表に合計欄を付け加える。
t1 <- table(sex,regu)
addmargins(t1)
regu
sex 0 1 Sum
0 10 4 14
1 3 8 11
Sum 13 12 25prop.table()関数で相対度数にできる。 デフォルトでは，全体相対度数となっているが，オプションでmargin = 1を指定すると行相対度数，margin = 2を指定すると列相対度数となる。
prop.table(t1)
regu
sex 0 1
0 0.40 0.16
1 0.12 0.32テキストで紹介されているgmodels関数を用いると，キレイな表ができるがコードが長くなるので，各自で学習する。
パッケージを読み込む。
# install.packages("gmodels")
library(gmodels)CrossTable()関数で作表する。
CrossTable(sex,regu,
expected = F, prop.r = T, prop.c = F,
prop.t = F, prop.chisq = F,
chisq = T, asresid = T, format = "SPSS")
Cell Contents
|-------------------------|
| Count |
| Row Percent |
| Adj Std Resid |
|-------------------------|
Total Observations in Table: 25
| regu
sex | 0 | 1 | Row Total |
-------------|-----------|-----------|-----------|
0 | 10 | 4 | 14 |
| 71.429% | 28.571% | 56.000% |
| 2.194 | -2.194 | |
-------------|-----------|-----------|-----------|
1 | 3 | 8 | 11 |
| 27.273% | 72.727% | 44.000% |
| -2.194 | 2.194 | |
-------------|-----------|-----------|-----------|
Column Total | 13 | 12 | 25 |
-------------|-----------|-----------|-----------|
Statistics for All Table Factors
Pearson's Chi-squared test
------------------------------------------------------------
Chi^2 = 4.811855 d.f. = 1 p = 0.02826461
Pearson's Chi-squared test with Yates' continuity correction
------------------------------------------------------------
Chi^2 = 3.205388 d.f. = 1 p = 0.07339608
Minimum expected frequency: 5.28 クラメールのVとχ2検定
\(r\)行かける\(c\)列のクロス表の行要素と列要素の関係の強さを表す尺度の1つにクラメールのV(Cremer’s V)がある。
\[ V = \sqrt{\frac{\chi^2}{N \times \min (r-1,c-1)}} \]
\(r\)行\(c\)列のクロス集計表の\(i\)行\(j\)列の観測度数を\(O_{ij}\)，期待度数を\(E_{ij}\)， \(i\)行の合計を\(n_{i \cdot}\)，\(j\)列の合計を\(n_{\cdot j}\)，全データの合計を\(n\)とおくと，分子の\(\chi ^2\)値は次式で計算される。
\[ \chi ^2 = \displaystyle \sum_{i = 1}^r {\displaystyle \sum_{j = 1}^c} \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]
ここで，期待度数\(E_{ij}\)は次式で求める。
\[ E_{ij} = \frac{n_{i\cdot} \times n_{\cdot j}}{n} \]
クロス表からクラメールのVと\(\chi ^2\)統計量を出力するため，vcdパッケージのassocstats()を用いる。
# install.packages("vcd") # first time only
library(vcd)
assocstats(t1)
X^2 df P(> X^2)
Likelihood Ratio 4.9748 1 0.025719
Pearson 4.8119 1 0.028265
Phi-Coefficient : 0.439
Contingency Coeff.: 0.402
Cramer's V : 0.439 クラメールのVは0.439と比較的強い関連性が確認できた。
せっかくなので，上記の関数を使わずに，定義通りに計算して\(\chi ^2\)値とクラメールの\(V\)を計算してみる。 クロス表は，
| あり | なし | |
|---|---|---|
| 男 | 10 | 4 |
| 女 | 3 | 8 |
ここから，期待度数を計算する。
\[ E_{\text{男},\text{あり}} = \frac{14 \times 13}{25} = 7.28\\ E_{\text{男},\text{なし}} = \frac{14 \times 12}{25} = 6.72\\ E_{\text{女},\text{あり}} = \frac{11 \times 13}{25} = 5.72\\ E_{\text{女},\text{なし}} = \frac{11 \times 12}{25} = 5.28 \]
つまり期待度数のクロス表は次のようになる。
| あり | なし | |
|---|---|---|
| 男 | 7.28 | 6.72 |
| 女 | 5.72 | 5.28 |
これをRで計算する。 期待度数を計算して，行列emを作成する。
m <- matrix(t1, nrow = 2, ncol = 2) # 実際の度数
em <- matrix(nrow=2,ncol=2) # 期待度数
for (i in 1:ncol(m)){
for (j in 1:nrow(m)) {
em[i,j] <- (sum(m[i,]) * sum(m[,j])) / sum(m)
}
}
em
[,1] [,2]
[1,] 7.28 6.72
[2,] 5.72 5.28ここから\(\chi^2\)値を計算し，chisq.test()関数の結果と比較する。
(chi2 <- sum( (m - em)^2 / em ))
[1] 4.811855
(res <- chisq.test(t1))
Pearson's Chi-squared test with Yates' continuity correction
data: t1
X-squared = 3.2054, df = 1, p-value = 0.0734あれ，手計算の4.812と，関数3.205の結果が違う。。
モザイク図を作成する。
mosaic(t1, shade=TRUE)