[連載]フリーソフトによるデータ解析・マイニング第49回

Rとカテゴリカルデータのモデリング(3)

1. 分割表のモデリング

 カテゴリカルデータを分割表にまとめ、そのセルの度数をモデリングする方法として、分割表の対数線形モデル(log linear model)という方法がある。説明を簡単にするために、まず表1に示す2元分割表を考えよう。

表1 2元分割表の一般形式と周辺度数
b1 b2 ・・・ bj ・・・ bc
a1 n11 n12 ・・・ n1j ・・・ n1c n1+
a2 n21 n22 ・・・ n2j ・・・ n2c n2+
・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
ai ni1 ni2 ・・・ nij ・・・ nic ni+
・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
ar nr1 nr2 ・・・ nrj ・・・ nrc nr+
n+1 n+2 ・・・ n+j ・・・ n+c n++

 2元分割表のセルの度数nijの期待度数2元分割表のセルの期待度数を対数変換すると、

2元分割表のセルの期待度数対数変換後
となる。式の中の項目-log(n++)は総度数の効果、log(ni+)は第i行の周辺度数の効果、log(n+j)は第j列の周辺度数の効果である。分割表が独立である場合は、この式で観測データを表現することができる。しかし、観測データが独立ではない場合は、行と列の相互効果を加えて表現することが必要である。グッドマン(Goodman)は1970年代に、次に示すような式で2元分割表の観測度数をモデリングする方法を提案した。
2元分割表の観測度数をモデリングする方法
 ただし、2元分割表の観測度数をモデリングする方法の条件を条件とする。式の中のλは総度数の効果、λiAは第i行の周辺度数の効果、λjBは第j列の周辺度数の効果、λijABは第i行と第j列の交互作用効果を表すパラメータである。このモデルを2元分割表の飽和モデル(saturated model)と呼ぶ。飽和モデルのパラメータの数は分割表のセルの数に等しい。交互作用効果のパラメータを持たないモデルを独立モデルと呼ぶ。
 モデルの中のパラメータは観測データの対数値log(nij)=vijを用いて、次のように求めることができる。
モデルの中のパラメータ
 パラメータを求めるプロセスを表2に示すデータを用いて説明する。表2の対数変換の結果を表3に示す。

表2 喫煙と肺癌の分割表
肺癌あり 肺癌なし
喫煙者 182 156
非喫煙者 72 98

表3 表2の対数変換の分割表
肺癌あり 肺癌なし 平均
喫煙者 5.204 5.050 5.127
非喫煙者 4.277 4.585 4.431
平均 4.740 4.817 4.779

 2×2の分割表におけるモデルの制約条件は、

2×2の分割表におけるモデルの制約条件
である。表1のセルn11=182を表現するモデルのパラメータは、
2×2の分割表におけるモデルの制約条件(表1を例に)
となり、モデルによるセルn11の推測値n11^
モデルによる推測値
となる。同様のモデルに基づいて求めた結果を表4に示す。

表4 表2の対数線形モデルの推測値
肺癌あり 肺癌なし
喫煙者 182 156
非喫煙者 72 98

 このモデルによる推測値は、観測値と完全に一致する。

2. 関数loglinによるモデリング

 Rには対数線形モデルのパラメータの推測や検定統計量を求める関数loglinがある。関数loglinの書き式を次に示す。

loglin(table, margin, fit = FALSE, param = FALSE, print =TRUE…)

 引数tableにはデータの分割表を、引数marginではリスト形式で分割表の範囲を、引数fitでは分割表の推測値を返すか否かを、引数paramでは推測したパラメータを返すか否かを論理値(FALSE, TRUE)で指定する。これ以外にも不完備表(incomplete tables、表の中に論理的にゼロとなるセルを含む分割表)に対応する引数startがある。
 表2のデータを用いた例を次に示す。まず次のように分割表を作成する。

> tab1<-matrix(c(182,72,156,98),nc=2)
> rownames(tab1)<-c("喫煙者","非喫煙者")
> colnames(tab1)<-c("肺癌あり","肺癌なし")
> (tab1<-as.table(tab1))

肺癌あり 肺癌なし
喫煙者 182 156
非喫煙者 72 98

 データtabを用いた関数loginの使用例を次に示す。コマンドの中のlist(c(1,2))は2元分割表の表側と表頭を表す。

> log.m1<-loglin(tab1,margin=list(c(1,2)), param=T,fit=T)

 返す項目は関数summaryを用いて確認できる。

> summary(log.m1)
データtabを用いた関数loginの使用例

 $lrtには尤度比の検定統計量、$pearsonにはピアソンのカイ2乗検定統計量、$dfにはモデルの当てはめの自由度、$marginにはモデルの当てはめのために指定した値が記録されている。$fitはモデルによって当てはめた結果、$paramには推測されたモデルのパラメータの値が記録されている。

> log.m1$fit

肺癌あり 肺癌なし
喫煙者 182 156
非喫煙者 72 98

 返された結果から分かるように、推測値は用いた観測値と全く同じである。

> log.m1$param

推測されたモデルのパラメータの値

 パラメータλ、λiA、λjB、λijABの推測値は上から順番にそれぞれ項目$`(Intercept)`、$`1`、$`2`、$`1.2`に返されている。  関数loglinは多元分割表のデータを扱うことも可能である。パッケージdatasetの中に髪の色(4種類)と目の色(4種類)を性別に分けた3元分割表データHairEyeColorがある。

> class(HairEyeColor)
[1] "table"
> dim(HairEyeColor)
[1] 4 4 2
> HairEyeColor

3元分割表データHairEyeColor

 この3元分割表のカテゴリカル変数である髪色、目の色、性別をそれぞれA=Hair、B=Eye、C=Sexにした飽和モデルの一般式を次に示す。

HairEyeColor飽和モデルの一般式
 ただし、パラメータの制約条件は、各カテゴリカル変数のパラメータの合計はゼロである。
HairEyeColor飽和モデルのパラメータの制約条件
 データHairEyeColorの3つのカテゴリカル変数の飽和モデルの推測例を次に示す。コマンド中の引数list(c(1,2),c(1,3), c(2,3))は、3つのカテゴリカル変数を組み合わせるモデルの書き式である。

> fm <- loglin(HairEyeColor,list(c(1,2), c(1,3),c(2, 3)),fit=T,param=T)

 比較のために、モデルにより推測された分割表の値を整数に丸めて次に示す。

> round(fm$fit)

モデルにより推測された分割表の値

 推測値の残差を次に示す。

> HairEyeColor-round(fm$fit)

推測値の残差

3. 関数loglmによるモデリング

 パッケージMASSには分割表を対数線形モデルで当てはめる関数loglmがある。関数loglmの書式は、線形回帰関数lmや一般線形回帰関数glmと似ている。
 関数loglmでは、分割表(表2)のデータを表5のような度数分布表の形式にして用いる。

表5 表2の度数分布表
喫煙歴 肺癌歴 度数
あり あり 182
あり なし 156
なし あり 72
なし なし 98

> tab2<-data.frame(喫煙=c(rep("あり",2),rep("なし",2)),
             肺癌=c(rep(c("あり","なし"),2)),
             度数=c(182,156,72,98))
> tab2

喫煙 肺癌 度数
1 あり あり 182
2 あり なし 156
3 なし あり 72
4 なし なし 98

 作成したデータセットtab2を飽和モデルに当てはめる例を次に示す。

> library(MASS)
> tab.m2<-loglm(度数~喫煙+肺癌+喫煙:肺癌,data=tab2)
> tab.m2
Call:
loglm(formula = 度数 ~ 喫煙 + 肺癌 + 喫煙:肺癌, data = tab2)

Statistics:

X^2 df P(> X^2)
Likelihood Ratio 0 0 1
Pearson 0 0 1

 このように関数loglmは、関数loglinと同じく尤度比とピアソンカイ2乗統計量を返す。次のように関数fittedとresidualsで当てはめ値と残差を返すことができる。

> fitted(tab.m2)
Re-fitting to get fitted values

肺癌
あり なし
喫煙 あり 182 156
なし 72 98


> residuals(tab.m2)
Re-fitting to get frequencies and fitted values

肺癌
あり なし
喫煙 あり 0 0
なし 0 0

 この結果は、関数loglinを用いた場合と同じであることが確認できる。
 Rでは関数updateを用いて作成したモデルの更新を行うことができる。モデルtab.m2の交互作用のパラメータを除いたモデルの書き式を例として次に示す。コマンドの中の.~.の右に追加、あるいは除くパラメータを指定する。追加する場合は+,除く場合は-を付ける。

> (tab.m3<-update(tab.m2,.~.-喫煙:肺癌))
Call:
loglm(formula = 度数 ~ 喫煙 + 肺癌, data = tab2)

Statistics:

X^2 df P(> X^2)
Likelihood Ratio 5.994097 1 0.01435383
Pearson 5.976471 1 0.01449799


> fitted(tab.m3)
Re-fitting to get fitted values

肺癌
あり なし
喫煙 あり 169 169
なし 85 85


> residuals(tab.m3)
Re-fitting to get frequencies and fitted values

肺癌
あり なし
喫煙 あり 0.9875737 -1.013250
なし -1.4484958 1.376219

4. モデルの選択

 分割表(あるいは度数分布表)のモデリングを行う際、飽和モデルを用いると当てはめはよいが、用いるパラメータの数が多いのが問題である。データをモデリングする際に重要なのは、少ないパラメータで簡潔なモデルを構築することである。
同一の分割表の対数線形モデルは唯一ではない。例えば、2元分割表のモデルとしては、次に示すようなモデルが考えられる。

2元分割表のモデル
 これらのモデルの中で、どのモデルを用いるべきであるかが1つの問題となる。
 通常、モデルの評価にはAICやBICなどの情報量規準がよく用いられている。これらは、返された尤度比とカイ2乗統計量を用いて簡単に計算することができる。関数loglmが返す尤度比(G2: Likelihood Ratio)を次の式の「-2(対数尤度)」に代入するとAICが計算できる。
AICの一般式
 2×2の分割表の飽和モデルのパラメータの数は4(松田 p.86)であるため、表2(あるいは表5)の飽和モデルのAICは0+2×4=8になり、交互作用効果を除いたモデル(独立モデル)のAICは5.994+2×3=11.994になる。
 AICを求める式から分かるように、AICは分割表の総度数の影響を無視している。しかし、分割表のカイ2乗値および尤度比は総度数の影響を受ける。各セルの比率は全く同じであっても総度数の増加に伴い分割表のカイ2乗値、尤度比が増大する。例えば、tab2の総度数は508であり、独立モデルの尤度比は5.994097である。各セルの比率が同じであり、総度数が2000の場合の独立モデルの尤度比は23.8992になる。

> tab3<-tab2
>tab3[,3]<-round(2000*tab3[,3]/sum(tab3[,3]),0)
> loglm(度数~喫煙+肺癌,data=tab3)

loglm(度数~喫煙+肺癌,data=tab3)

 分割表の総度数の影響を考慮して、モデルを評価する情報量規準としてBIC(Bayesian Information Criterion)がある。分割表のBICは次の式で定義されている。

BICの一般式
 データtab2の飽和モデルのBICは0+log(508)×4=24.9219、独立モデルのBICは5.994097+log(508)*3=24.6855となる。
 情報量規準を用いてモデルの選択を行う場合には、値が小さいモデルを選択する。上記のデータに対してどのモデルを用いるかに関しては、用いる情報量規準によって結果が異なる。このデータにおいては、AICを用いると飽和モデルが選択され、BICを用いると独立モデルが選択される。
 対数線型モデルのAIC, BICは関数extractAICで求めることができる。引数k=2のときはAIC、k=log(n++)のときはBICを返す。デフォルトにはk=2になっている。モデルtab.m2を用いた例を次に示す。

> extractAIC(tab.m2)
[1] 4 8        #AIC
> extractAIC(tab.m2,k=log(sum(tab2[,3])))
[1] 4.00000 24.92193 #BIC

5. 多元度数表と対数線形モデル

 パッケージMASSの中に、住宅環境に対する満足度に関する調査結果を4元度数表にまとめたデータセットhousingがある。

> data(housing)
> head(housing,n=3)

Sat Infl Type Cont Freq
1 Low Low Tower Low 21
2 Medium Low Tower Low 21
3 High Low Tower Low 28

 変数Satは現在の住宅状況に対する満足度(高、中、低)、変数Inflは不動産のマネジメントにおける影響力の認知度(高、中、低)、変数Typeは賃貸住宅のタイプ(高層ビール、中庭付き、アパート、テラス)、変数Contは他の住民とのふれあいの度合い(高、低)、変数Freqは度数である。変数のすべての組み合わせは3×3×4×2=72通りであるので、データセットは72行5列である。
 データhousingの分割表形式は次のコマンドで確認できる。

> xtabs(Freq~Sat+Infl+Type+Cont, data=housing)

データhousingの分割表形式

<後略>

 関数loglmを用いてモデルを作成する例を次に示す。

> hous.log<-loglm(Freq~Sat+Infl+Type+ Cont,data=housing)
> hous.log
loglm(formula = Freq ~ Sat + Infl + Type + Cont, data = housing, fitted = T, param = T)

Statistics:

X^2 df P(> X^2)
Likelihood Ratio 295.3518 63 0
Pearson 305.9267 63 0

> hous.log$param
<省略>

 関数glmでは、リンク関数をポアソン(poisson)に指定することにより対数線形モデルに当てはめることができる。その例を次に示す。

> hous.glm<-glm(Freq~Sat+Infl+Type+Cont, data=housing,family=poisson)

 計算された尤離度は関数loglmの尤度比に等しい。

> deviance(hous.glm)
[1] 295.3518

 関数summaryは回帰係数と標準誤差などを返す。

> summary(hous.glm)
<結果は省略>

 関数dummy.coefは、すべてのパラメータを返す。

> dummy.coef(hous.glm)
<結果は省略>

[1] 松田 紀之 著(1988). 質的情報の多変量解析、朝倉出版 (http://www.sci. kagoshima-u.ac.jp/~ebsa/matsuda01/index.html)
[2] http://www.stat.ufl.edu/~aa/cda/cda.html (Categorical Data Analysis by Alan Agresti (Wiley, 2002)のRコード)