カテゴリカルデータの解析 (その2)
その1からのつづきです．
変数間の関係性を調べる
集計したカテゴリカルデータを解析する方法を紹介します．独立性の検定
「変数のあいだに関連性があるか」を調べる検定です．先に紹介したArthritisで，処置と改善度合いのあいだに関連性があるかを調べてみましょう．この場合には，カイ二乗検定 (データ数が少ない場合にはフィッシャーの直接確率検定のほうがよい) を使います．
帰無仮説「処置と改善度合いに関連性はない」を，カイ二乗検定で検討してみます．
> chisq.test(arthritis.imp.tre) Pearson's Chi-squared test data: arthritis.imp.tre X-squared = 13.055, df = 2, p-value = 0.001463
結果，p値 < 0.0015で，対立仮説「処置と改善度合いには関連性がある」をとることになります．
カイ二乗検定の用途には，他に，適合度の検定があり，「観察値と理論値が同じか」を調べることができます．
カイ二乗検定の原理については，以下などがわかりやすいかもしれません．
『カテゴリカルデータ解析 (Rで学ぶデータサイエンス1)』には，以下に示したように，ほかにもいろいろな場合に使う検定法が記されていますので，参考にしてみてください．
- マンテル検定: カテゴリー間の順序関係をスコアに変換/li>
- ウィルコクソン順位和検定: スコアではなく順位の状態で検定
- コクラン・アーミテージ検定: ロジスティック回帰分析のように順序カテゴリーを対象とする
- クラスカル・ワリス検定: 3群以上の比較に拡張したウィルコクソン検定
- マクネマー検定: 変数間に対応関係がある場合に使う
回帰分析
ひとつの変数を目的変数として，目的変数を，その他の変数で説明するかたちのモデルを用いて分析を行なう手法です．これについては以前まとめました．
ほかにも，目的変数が個数データのときに用いるポアソン回帰モデルや，目的変数を特に定めないで変数の関係性を調べるときに用いる対数線形モデルがありますが，ここでは大胆に省略します:(
対応分析
カテゴリカル変数間の関係をうまく表そうとする方法．複数変数のカテゴリーに，相関が高くなるような数値を割り当てて，関係を分析します．たとえば，DanishWelfareのデータで，婚姻状態と居住地域の関係を表してみましょう．
> danish.sta.urb <- xtabs(Freq ~ Status + Urban, data = DanishWelfare)
> round(prop.table(danish.sta.urb, margin = 2), 3)
Urban
Status Copenhagen SubCopenhagen LargeCity City Country
Widow 0.210 0.148 0.093 0.121 0.064
Married 0.522 0.700 0.705 0.705 0.773
Unmarried 0.268 0.151 0.202 0.173 0.163
> chisq.test(danish.sta.urb)
Pearson's Chi-squared test
data: danish.sta.urb
X-squared = 158.1145, df = 8, p-value < 2.2e-16独立性に関するカイ二乗検定を行なうと，p値 < 0.001で，これらふたつの変数間に何らかの関係性がありそうです．
対応分析では，それぞれのカテゴリーに数値を割り当ててみて，変数どうしの相関が高くなる最適な数値の組み合わせを探索します．
最適化問題はcorresp()関数で解くことができます．それぞれのカテゴリーを2次元の値で表現するように指定 (nf = 2) して，結果のベクトルを2次元平面上に表してみましょう．
> danish.corresp <- corresp(danish.sta.urb, nf = 2) > plot(danish.corresp)
割り当てられた数値が近いほど関連が近いと考えられ，平面上でもカテゴリー名が近くに表示されます．Copenhagen以外の地域は比較的まとまっており，これらは同じような傾向を示していることがみてとれます．また，この4地点の近くにMarriedが表示されており，これらの地域では婚姻者の割合が大きこともわかります．
このようにして，変数のカテゴリー間の対応関係を調べる方法が対応分析です．
決定木
説明変数を利用して各個体を段階的に分類していき，どのカテゴリーに属するかを決定する方法です．フィッシャーのアヤメのデータを，決定木を使って解析してみます．mvpartまたはrpartパッケージを使います．
> install.packages("mvpart")
> library(mvpart)
> (iris.part <- rpart(Species ~ ., data = iris))
n= 150
node), split, n, loss, yval, (yprob)
* denotes terminal node
1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)
2) Petal.Length< 2.45 50 0 setosa (1.00000000 0.00000000 0.00000000) *
3) Petal.Length>=2.45 100 50 versicolor (0.00000000 0.50000000 0.50000000)
6) Petal.Width< 1.75 54 5 versicolor (0.00000000 0.90740741 0.09259259)
12) Petal.Length< 4.95 48 1 versicolor (0.00000000 0.97916667 0.02083333) *
13) Petal.Length>=4.95 6 2 virginica (0.00000000 0.33333333 0.66666667) *
7) Petal.Width>=1.75 46 1 virginica (0.00000000 0.02173913 0.97826087) *結果の意味は以下のとおりになります．
最初のnは，対象としたデータの行数を表しています．
次の行は，出力値の意味を示しています．nodeは分岐する点の順番，splitはデータの分け方，nはそのグループに用いた行数，lossは予想が外れた数，yvalはグループの予想値，(yprob)はそれぞれのカテゴリーの割合 (アルファベット順に配置) を表すそうです．*は，分岐の終末点を表します．
このままではわかりづらいので図にしましょう．
> plot(iris.part) > text(iris.part, use.n = TRUE, all = TRUE)
ノード1) のPetal.Lengthによる分類で，setosaと{versicolor, virginica}が100%分類できていることがわかります．
後者は，ノード3) のPetal.Widthによる分類でさらに分けられていますが，versicolorに分類されてしまったvirginicaが5/50個体，virginicaに分類されてしまったversicolorが1/50個体いることがわかります．
ノード 6)では，Petal.Lengthによる分類で，versicolor 49個体とvirginica 5個体を分けています．
irisデータの場合，説明変数はカテゴリーではなく量的変数です．そうした場合，さまざまな可能性が探索されながら，もっとも良い値が各変数の分割点に採用されるようになっています．
説明変数がカテゴリーとして表されている場合でも決定木の方法は使えます．
また，predict()関数を使った予測や，rpartOrdinalパッケージを使った順序カテゴリカルデータに対する方法もありますが，ここでは省略します:( 詳細は『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』をご覧ください．
みせかけの相関
最後に，みせかけの相関の問題について．母集団全体と，いずれかの変数に着目して分割した集団とで，相関関係が異なっている場合があるので気をつけましょう，というものです．
例によって，ここでは詳しく述べませんが，以下などを参考にしてみてください．
シンプソンのパラドックス (Wikipedia)
シンプソンのパラドックスの図解 (社会学者の研究メモ)
ふたつの変数間に相関が見られたとしても，両方に関連する別な因子がないかを常に意識することが大切ですね．
参考
Kashiwa.R#3で発表する内容を文章にしたものです．内容のほぼすべては，『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』 藤井良宜 2010年 共立出版 を大いに参考にしています．
