3日目：クロス集計表

2変数の関係を確認する方法として，クロス集計表と独立性の検定である \chi^2 検定について学習します。以下では、前回学習したtidyverseのパッケージを用いて、データの読み込みやデータの抽出・加工、グラフの作成を行います。

pacman::p_load(tidyverse)

データの読み込みと確認

chap5 <- read_csv("data/chap5.csv")

Rows: 25 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (4): ID, sex, satis, regu

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

table()関数を用いてシンプルなクロス表を作成します。 table()関数は引数として2つのベクトルを受け取り、クロス集計表を作成します。このとき、ベクトルの要素がカテゴリカル変数となっている必要があります。

度数のみのシンプルなクロス表はtable()関数を用いて作成しますが、連続変数の場合は、cut()関数を用いてカテゴリカル変数に変換する必要があります。

table(chap5$sex, chap5$regu) # クロス表を作成

addmargins()を用いて，クロス表に合計欄を付け加えます。

周辺合計を追加するならaddmargins()

t1 <- table(chap5$sex, chap5$regu)
addmargins(t1)

     
       0  1 Sum
  0   10  4  14
  1    3  8  11
  Sum 13 12  25

列合計と行合計が追加されていることが確認できました。

prop.table()関数で相対度数にしてみます。デフォルトでは，全体相対度数となっているが，オプションでmargin = 1を指定すると行相対度数，margin = 2を指定すると列相対度数となります。

相対度数のクロス表はprop.table()関数で，margin = 1で行相対度数，margin = 2で列相対度数のクロス表になる

prop.table(t1)

   
       0    1
  0 0.40 0.16
  1 0.12 0.32

テキストで紹介されているgmodels関数を用いると，キレイな表ができるがコードが長くなるので，各自で学習する。

パッケージを読み込む。

pacman::p_load(gmodels)

CrossTable()関数で作表する。

gmodelsパッケージのCrossTable()関数で情報量の多いクロス表が作成できる

CrossTable(chap5$sex, chap5$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 

             | chap5$regu 
   chap5$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と \chi^2 検定

r 行 \times c 列のクロス表の行要素と列要素の関係の強さを表す尺度の1つにクラメールのV(Cremer’s V)があります。クラメールの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()を用いる。

pacman::p_load(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

ここから，期待度数を計算する。

\begin{align*} 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 \end{align*}

つまり期待度数のクロス表は次のようになる。

	あり	なし
男	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)  
  }
}
print(em) # 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の結果が違う。。