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

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

1. モデリング

 データサイエンスの分野では、観測データからノイズを取り除き、一定の法則を見つけ出して抽象化することをモデリングと呼ぶ。量的データの最も簡単なモデルは回帰分析である。本欄の第13回〜16回(2004年8月号〜11月号)で線形・非線形回帰モデルについて説明した。本稿では、カテゴリカルデータのモデリングについて説明する。
 モデリングには、応答変数が何らかの確率分布に従うという仮定の下で、モデルに必要となる係数・パラメータを推測する方法が最も多く用いられている。一般の線形回帰分析はデータが正規分布に従うという仮定の下で、モデルの推定を行う。
 カテゴリカルデータの場合は、観測データが2項分布、ポアソン分布、多項分布、などの確率分布に従うと見なし、モデルを推測する。
 しかし、何らかの仮定の下で構築したモデルが真のモデルであるかどうかは判断できない。同一のデータについて異なる仮定の下で推定した複数のモデルの中で、どのモデルが真のモデルに最も近似しているかを計測する統計量は幾つか提案されている。統計データ解析パッケージの中に最も多く実装されているのは赤池情報量規準AIC(Akaike Information Criterion)である。

2. 2項分布モデル

 観測データの応答変数(目的変数、あるいは被説明変数)が2とおりであると見なす場合は、2項分布に従うと仮定することができる。例えば、コインを投げる実験を行ったとき、表と裏が現れることのみを考えると実験結果は2項分布に従う。政権に対する支持率の調査結果を支持と不支持のみに分けて集計したデータ、薬物実験を行ったとき、効果ありと効果なしのみをカウントしたデータはいずれも2項分布に従う。
 2項分布に従う実験をn回行ったとき、事象Aがy回起こる確率はPy=nCypy(1-p)n-yとなる。pは事象Aが起こる確率、q=1-pは事象Aが起こらない確率である。
 2項分布に従うデータのモデリングには、次に示す事象 が起こる確率と、起こらない確率の対数オッズの線形回帰モデルが広く用いられている。

対数オッズの線形回帰モデル

 この式の左辺をロジット(logit)関数と呼び、次に示すその逆関数をロジスティック(logistic)関数と呼ぶ。

ロジスティック(logistic)関数ル

 この式を用いた回帰をロジスティック回帰と呼ぶ。本稿ではカテゴリカルデータのモデリング視点で解説することにする。

3. 関数glmによるモデルの推測

(1) 度数データのモデリング

 具体的なデータを用いて説明を行うことにする。データはパッケージfarawayの中のデータセットbabyfoodを用いることにする。パッケージfarawayはCRANミラーサイトからダウンロードできる。パッケージがインストールされている場合は次のコマンドでパッケージの中のデータを読み込むことができる。

> data(babyfood,package="faraway")
> babyfood

データbabyfood

 データbabyfoodは乳児の病気状況(あり、なし)を性別および3つの異なる形態の飲食物の与え方に分けて記録したデータである。これは表1に示す3元分割表に等しい。

表1 babyfoodの3元分割表
disease
Bottle Suppl Breast
Boy 77 19 47
Girl 48 16 31
nondisease
Bottle Suppl Breast
Boy 381 128 447
Girl 336 111 433

 ここでは病気の確率を応答変数とする。病気の確率は変数diseaseとnondiseaseの度数から求めることができる。しかし、Rの中の一般線形化モデルを推測する関数glmは、事象 が「起こる」と「起こらない」の度数データを用いてロジスティック回帰モデルを推測できるように設計されている。
 関数glmにおいて応答変数を指定する形式は、事象Aが起こる確率を指定する方法と、事象Aが起こる度数と起こらない度数を同時に指定する方法がある。後者の指定方法を用いたコマンドの例を下に示す。
 関数glmの書き式は、線形回帰関数lmに似ている。ただし、応答変数diseaseとnondiseaseの度数を用いる際には、cbind(disease, nondisease)のようにまとめる必要がある。また、関数glmを用いてロジスティック回帰分析を行う場合は、引数familyに2項分布(binomial)を指定しなければならない。そのコマンドを次に示す。

> mod1<-glm(cbind(disease,nondisease)~sex+food, family=binomial,data= babyfood)

 推測された係数および関連の統計量は関数summaryで返すことができる。

> summary(mod1)

推測された係数および関連の統計量

 返された結果の中の逸脱残差(Deviance Residuals)は用いたデータの各ケースの逸脱度である。値は以下に示す残差の逸脱度(Residual deviance)の平方根に、観測度数から予測度数を引いたときの符号(正、負)を付けている。Residual devianceは次の式で定義されている。式の中のyiは観測度数、yi^はモデルによる予測度数である。niはケースごとの総度数である。

Residual devianceの定義

この2つの逸脱度(Null deviance、Residual deviance)は、説明変数を用いていないモデルと、説明変数を用いたモデルとの比較分析を行う際によく用いられる。逸脱度については、関数anovaを用いて次のように分散分析を行うことができる。

> anova(mod1,test="Chi")

関数anovaを用いた分散分析の結果

 コマンドsummary(mod1)が返す結果の逸脱度と自由度は、それぞれコマンドdeviance(mod1)、df.residual(mod1)で呼び出すことができる。また、Coefficientsの部分は、モデルに用いた係数の推測値、標準誤差などである。次のように係数だけを呼び出すこともできる。

> round(coefficients(mod1),4)

(Intercept) sexGirl foodBreast foodSuppl
-1.6127 -0.3126 -0.6693 -0.1725

 返された結果から分かるように、モデルに用いる説明変数の係数はsexGirl、foodBreast、 foodSupplである。推測された係数を用いたモデルは、

推測された係数を用いたモデル1
推測された係数を用いたモデル2

となる。このようなモデルでは、説明変数の性別がGirlの場合は、式の中の は1,そうではない場合は0とする。説明変数foodがBreastであると式の中のBreastは1、Supplは0とする。説明変数foodがBottleである場合はBreast、Suppl両変数を0とする。このように、カテゴリカルデータの場合はダミー変数0,1を用いる。
 データbabyfoodの第1行のデータはsex=boy, food=Bottleである。この値を用いた、第1行のdiseaseの予測確率は、

第1行のdiseaseの予測確率1
第1行のdiseaseの予測確率2

であり、その予測度数は0.1662×(77+381)≒76である。
 このように計算されたロジット値と予測確率は、関数predictとfittedで返すことができる。

> round(predict(mod1),3)

1 2 3 4 5 6
-1.61 -1.79 -2.28 -1.93 -2.10 -2.59

 返された値はロジット関数値であり、ロジスティック関数に代入すると応答変数の予測の確率値が求まる。予測の確率値を求めるには、関数predictに引数type=”response”を用いる方法と、関数fittedを用いる2とおりの方法がある。

> round(fitted.values(mod1),3)

1 2 3 4 5 6
0.166 0.144 0.093 0.127 0.109 0.069

 diseaseの予測度数は次のように、求めることができる。

> round(fitted(mod1)*(babyfood[,1]+babyfood[,2]),0)

1 2 3 4 5 6
76 21 46 49 14 32

(2) 生データのモデリング

 データbabyfoodのように度数表に集計していない2値データとして記録したデータのロジスティック回帰の例を示す。
 パッケージISwRの中に健康診断で収集したデータjuulがある。パッケージISwRはCRANミラーサイトからダウンロードできる。データjuulは1339の個体、6変数のデータフレームである。本稿で用いる変数のみを次に示す。

age 年を単位とした年齢。
menarche 初経に関する量的データ。なしが1、ありが2
sex 性別に関する量的データ。boyが1, girlが2
igf1 量的データ。インシュリン様成長因子の測定値

 データを読み込み、ランダムに4行を取り出したデータを示す。

> data(juul,package="ISwR")
> set.seed(10)
> juul[sample(nrow(juul),4),]

ランダムに4行を取り出したjuulデータ

 カテゴリカルデータとしてモデルの推測を行うため、変数menarcheの1,2を次のコマンドで"no","yes"に変換する。

> juul$menarche<-factor(juul$menarche, label=c("no","yes"))

 初経に関するモデルを構築するので、次のコマンドで、年齢が8歳から18歳までのデータを抽出し、その中の1列から4列までに欠損値を含む個体を取り除いたデータセットjuul.girを作成する。

> temp<-subset(juul,age>8 & age<18)
> juul.gir<-na.omit(temp[,1:4])
> head(juul.gir,n=3)

age menarche sex igf1
343 13.01 no 2 682
745 8.13 no 2 210
746 8.17 no 2 564

 変数menarcheを応答変数(2値データ)、変数ageを説明変数としたロジスティック回帰モデルを推測する。

> m1<-glm(menarche~age,data=juul.gir, binomial)
> summary(m1)

ロジスティック回帰モデルの推測

 応答変数yiが2値データの場合の逸脱度は

逸脱度

となる。ただし、yiは0,1を取る。式の中のpi^はモデルによる予測値(確率)である。
 モデルに用いた説明変数は1つであり、かつ量的データであるので、予測値と観測値は図1のように図示することができる。横軸は変数ageであり、縦軸は応答変数の確率である。図の中の実線は、ロジスティック回帰の予測値、白丸点は観測値である。まず実線を作成するコマンドを次に示す。

> newage<-data.frame(age=seq(8,18,0.1))
> m1.pr<-predict(m1,newage,type= "response")
> plot(newage$age,m1.pr,type="l")

 図1の中の白丸点は、観察データを年ごとにまとめた比率データである。そのデータを作成し、実線プロットに追加するコマンドを次に示す。

> age.group<-cut(juul.gir$age,c(8:18))
> m1.pr<-predict(m1,newage,type="response")
> tb<-table(age.group, juul.gir$menarche)
> rel.freq<-prop.table(tb,1)[,2]
> points(c(8:17)+0.5,rel.freq,pch=1)

モデルm1の予測値と観測値

図1 モデルm1の予測値と観測値

 変数menarcheを応答変数、変数ageとigf1を説明変数としたロジスティック回帰モデルの推測の例を次に示す。

> m2<-glm(menarche~age+igf1,data=juul.gir, binomial)
> summary(m2)

ロジスティック回帰モデルの推測の例

 両モデル(m1,m2)のAICを比べると、明らかに変数ageとigf1を用いたm2の方が小さい。
 モデルの観測データへの当てはめの良さは混同マトリクスからも考察することができる。
 応答変数menarcheは、カテゴリカルの2値データであるが、ロジスティック回帰の予測結果は0から1までの確率値である。予測された確率値と2値データとの対応関係は、確率値0.5を境界として二分して用いることが1つの方法である。混同マトリクスと正解率を求めるコマンドを次に示す。

> m1.p<-round(predict(m1,type="resp"))
> (m1.t<-table(juul.gir$menarche,m1.p))

m1.p
0 1
menarche no 168 15
yes 17 166

> sum(diag(m1.t))/sum(m1.t)
[1] 0.9125683

> m2.p<-round(predict(m2,type="resp"))
> (m2.t<-table(juul.gir$menarche,m2.p))

m2.p
0 1
menarche no 168 15
yes 8 175

> sum(diag(m2.t))/sum(m2.t)
[1] 0.9371585

 モデルm2の正解率がm1より約2.5ポイント高い。このようなロジスティック回帰モデルは非線形2群判別分析の方法として用いることが可能である。

4.2変数ロジスティック回帰モデル

 パッケージVGAMにはデータcoalminersがある。パッケージVGAMはCRANミラーサイトから入手できる。データcoalminersは、炭鉱夫の呼吸器に関して、次に示す項目に分けて整理した9行5列のデータフレームである。

BW 息切れあり、喘鳴ありの度数
BnW 息切れあり、喘鳴なしの度数
nBW 息切れなし、喘鳴ありの度数
nBnW 息切れなし、喘鳴なしの度数
age 炭鉱夫の年齢について5歳きざみで年を1つのカテゴリーとした代表値

> library(VGAM);data(coalminers)
> coalminers

coalminersのデータ

 ロジスティック回帰モデルのアイディアに基づいて、「息切れあり・なし」と「喘鳴あり・なし」のような2つの2項分布を応答変数とした回帰を2変数ロジスティック回帰(Bivariate Logistic Regression)と呼ぶ。2変数ロジスティック回帰分析の応答変数は、2つ(Y1,Y2)の2値データであるため、データcoalminersのように4とおりの可能性がある。この2変数は説明変数を指すことではない。
 パッケージVGAMの中の関数vlgmとリンク関数binom2.orを用いて2変数ロジスティック回帰分析を行うことができる。
 関数vlgmの書き式は、関数glmの場合と似ている。データcoalminers の中の変数BW、BnW、nBW、nBnWを応答変数とし、変数ageを説明変数とした2変数ロジスティック回帰の例を次に示す。

> coa.vglm<-vglm(cbind(BW,BnW,nBW,nBnW)~age,family=binom2.or, data=coalminers)
> (coa.fit<-fitted(coa.vglm))

2変数ロジスティック回帰の例
<後略>

 現段階のバージョンでは、返される結果の列ラベルが入力した応答変数と対応していない。ラベルを対応させる1つの方法は、次のように、返された結果の列のラベルを、用いた応答変数のラベルで入れ替えることである。

> lab<- colnames(coalminers[,1:4])
> colnames(coa.fit)<-lab
> coa,fit

2変数ロジスティック回帰の例ラベル入れ替え
<後略>

 モデルにより予測される度数は次のコマンドで求めることができる。

> round(coa.fit*apply(coalminers,1,sum))
<結果は省略>

 ここで用いた説明変数は1つであり、かつ量的な変数として見なすことが可能であるのでは、予測値を次のように図示することができる。

> x<-coalminers[,5]
> matplot(x,coa.fit,type="b",xlab="年齢",ylab="予測値(確率)")
> legend(24,0.6,lab,lty=1:4,pch=c("1","2","3","4"),col=1:4)

回帰モデルの予測値

図2 回帰モデルの予測値