R　大量のカテゴリ変数から変数間の連関度を表すクラメールV指標を計算し、有意で連関があるペアのみ抽出する

##目的

大量のカテゴリ変数からなるデータセットから、関連の強いとみられる変数ペアを[クラメールV連関係数](http://software.ssri.co.jp/statweb2/gloss/glossary_k.html#k027)により探索します。カテゴリデータのクロス表の統計量には、χ2値、ファイ係数、クラメールV連関係数などがありますが、今回は、Rの[vcdパッケージのassocstat関数](https://www1.doshisha.ac.jp/~mjin/R/45/45.html)を使い、データフレーム内の変数間のクラメールV係数を一気に求め、連関度が強く有意な変数間のペアのみ抽出します。

##想定
連関度が高いと判断する閾値：クラメールV0.3以上([こちらを参考にしました](http://langstat.hatenablog.com/entry/20150103/1420256759))
帰無仮説：２つのカテゴリー変数間には連関がない(独立である）
有意性検定のα値：5％

##アウトプットイメージ
i列とj列は、入力データフレームの変数の連番です。
![image](https://qiita-image-store.s3.amazonaws.com/0/114127/a91e32c6-57f9-06af-3c5c-6c52d0f18805.png)


##インプットデータのイメージ
1)カテゴリーデータのデータセット
このようなダミー変数（０と１からなる）のデータか

 var_1,var_2,var_3
   1,      1,    0,
   0,      1,    1,
   0,      0,    1,
・・・

このようなカテゴリーデータを想定しています。
（数値にみえますが、カテゴリーと思って下さい）
 var_1,var_2,var_3
  1,  5,   3,
  2,  1,   2,
  1,  5,   2,
・・・


2)変数連番と変数名のデータセット

![image](https://qiita-image-store.s3.amazonaws.com/0/114127/22098f2a-8c6c-35ea-dbe4-1a513bcd064e.png)



##assocstats関数の結果
![image](https://qiita-image-store.s3.amazonaws.com/0/114127/f04acc37-c4a4-05fe-21dc-f151485afb79.png)

- Likelihood Ratio(尤度比統計量）
- Pearson (χ2値)
- Phi-Coefficient (ファイ係数)
- Contingency Coeff (ピアソンの一致性係数)
- Cramer'sV (クラメールのV係数)

が出力されます。

クラメールVが5％水準で有意（帰無仮説を棄却できる）かどうかは、ピアソンのχ2値のp値を使って判定できます[(根拠とした教科書 p.785)](http://uregina.ca/~gingrich/ch11a.pdf)。


##Rのサンプルコード
```R:cramer_v.R
#データを読み込む
df<read.table("testdata.dat",header=T)


#変数ラベルの読み込み
var_list<-read.table("var_list",header=T)
head(cat_list)

#vcdを読み込みます。
library(vcd)

#全ての変数ペアについてクロス表を作成し、クラメールVを求めます
#変数の数
n<-ncol(df)
#--計算する変数ペアの数
nn<-0
for(i in 1:(n-1)){
  nn<-nn+(n-i)
}

#アウトプット用のマトリクス
output<-matrix(0,nrow=nn,ncol=4)
k<-0

for(i in 1: (n-1)){
  for(j in (i+1):n){
    #クロス表
    cross<-xtabs(~df[,i]+df[,j],data=df)
    #assocstats関数より、クロス表の分析統計量をresに格納する
    res<-assocstats(cross)
    #chisq_testsオブジェクトの出力のうちピアソンのχ2値のp値（P(>X^2)のセルがp値
    p_val<- res$chisq_tests[2,3]
    #cramerオブジェクトが、クラメールV係数
    cramer_v<-res$cramer
    #クラメールVが0.3以上、かつp値が5％以下の場合のみ保持する
    if(cramer_v >=0.3 && p_val <=0.05){
      k<-k+1
      output[k,1]<-i
      output[k,2]<-j
      output[k,3]<-cramer_v
      output[k,4]<-p_val
    }
  }
}
#変数名をつける。
colnames(output)<-c("i","j","cramer_V","p_val")

head(output)

#余計な行を削除する
output<-output[1:k,]

#変数ペアの１こめについて、変数ラベルを結合する.
output.2<-merge(output,var_list,by.x="i",by.y="var_no",all.x=TRUE)

#変数番号の後に、変数ラベルにする。
output.2<-output.2[c("i","varname","j","cramer_V","p_val")]

#varnameの変数名を変える（あとに同じ名前のフィールドをさらに結合するため）
colnames(output.2)<-c("i","varname.i","j","cramer_V","p_val")

#変数ペアの２こめ(j)について、変数ラベルを結合する。
output.2<-merge(output.2,var_list,by.x="j",by.y="var_no",all.x=TRUE)

#変数番号ｊの後に変数ラベルという順番に入れ替える。
output.2<-output.2[,c("i","varname.i","j","varname","cramer_V","p_val")]

#ｊの変数ラベルをvarname.jに変更する。
colnames(output.2)<-c("i","varname.i","j","varname.j","cramer_V","p_val")

#macのExcelで開くことを想定した出力
write.table(output.2,"クラメールV有意なペア.txt",sep="\t",fileEncoding = "CP932",row.names=F)
```

##クラメールV係数の利点
- ０〜１の範囲をとり、１に近いほど「２つの変数に連関の度合いが強い」と判断できるため、直感的に分かりやすい。
- クラメールVは、χ2値と異なりサンプルサイズによる影響を受けない。

##ファイ係数とクラメールV係数との使い分け
2×２クロス表　：ファイ係数
r×ｃクロス表：　クラメールV

2×2クロス表のクラメールVは、ファイ係数と同じ。
ファイ係数もクラメールVもカイ二乗値を元に計算できる（[式はこちら](https://ja.wikipedia.org/wiki/%E5%88%86%E5%89%B2%E8%A1%A8)）

##ファイ係数と相関係数の使い分け
量的な変数の線形の関係の度合い：相関係数
ダミー変数の相関係数：ファイ係数

##参考）実行時間
およそ1,000変数で１時間10分程度

実行環境
- CPU  Intel Core i5 2.7GHz
- メモリ 16GB
- OS OS X El Capitan
- R version 3.2.2 (2015-08-14)
- Platform: x86_64-apple-darwin13.4.0 (64-bit)
- Running under: OS X 10.11.3 (El Capitan)