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

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

1. ポアソン分布

 実験・観察で事象Aが、ある一定の時間内で生起する回数(あるいは確率)を考える問題のモデルとして、次に示すポアソン分布(Poisson distribution)がある。

ポアソン分布の式

ポアソン分布は、2項分布のnを限りなく大きくすることによって得られることが知られている。ポアソン分布のパラメータμと2項分布のn,pは、μ=npの関係を持っている。
 ポアソン分布は、与えられた単位時間内で事象Aがy回起こる確率を示す。事象がポアソン分布に適応するためには、次に示す条件を満たすことが必要である。

 (1) 事象Aが同時に2回起こらない。
 (2) 事象の生起は独立である。
 (3) 単位時間内の事象の平均生起の数は一定である。

 大きな集団の中で起こる偶発的事故や病気の頻度、コールセンターに掛かってくる電話の回数などはポアソン分布に従うと仮定して解析することが多い。
 事象Aの生起するパラメータがμiであるとき、事象Aがyi回生起する確率は、

ポアソン分布の式2
 このポアソン分布を用いた回帰をポアソン回帰(Poisson regression)と呼ぶ。ポアソン回帰において、推定すべきパラメータはμである。その推定方法としては、次に示す対数線形モデルが多く用いられている。
対数線形モデル
 式の中のパラメータは対数尤度を最大化することで求めることができる。
 パッケージfarawayの中のデータgalaを用いてポアソン回帰を説明することにする。

> data(gala,package=”faraway”)
> dim(gala);head(gala)

ポアソン回帰

 データgalaは、ガラパゴス(Galapagos)諸島の30の島と亀の種類との関連の6つの変数を記録した30行7列のデータフレームである。その7つの変数を次に示す。

Species:その島の亀の種類の数
Endemics:亀固有種の数。
Area:島の面積(km2)
Elevation:島の標高(m)。
Nearest:最近隣の島との距離(km)。
Scruz:Santa Cruz島との距離(km2)
Adjacent:近隣の島のエリア(km2

 ポアソン回帰モデルは関数glmとリンク関数poissonを用いて次のように推測することができる。

> gala.pm1<-glm(Species~ ., data = gala, family = poisson)
> summary(gala.pm1)

関数glmとリンク関数poissonを用いた推測

 ポアソン回帰の残差の逸脱度(Residual deviance)は次の式で定義される。

ポアソン回帰の残差の逸脱度

 データの行う数をn、モデルに用いた変数の数をkとすると、Dは自由度(n-1-k)のχ2分布に従う。式の中のμiは回帰モデルによる予測値であり、fitted(gala.pm1)で呼び出すことができる。
 モデルgala.pm1のsummaryで分かるように用いた変数Elevation、Scruz、AdjacentのPr(>|z|)は十分小さくない。そこで、関数stepを用いて変数を自動的に選択して、モデルを作成してみる。

> gala.pm2<-step(gala.pm1)
<返される結果は省略>
> summary(gala.pm2)

pm2のsummary

 関数stepを用いて、変数の選別を行ったモデルgala.pm2では、変数Scruz、Adjacentが取り除かれ、両モデルのAICの値を比較するとgala.pm2の方が若干減少している。小数点以下4桁まで丸めたモデルgala.pm2の係数を呼び出すコマンドとその結果を次に示す。

> coe<- coefficients(gala.pm2)
> round(coe,4)

モデルgala.pm2の係数を呼び出した結果

 この結果を用いたモデルを次に示す。

最終モデル

 データセットgalaの第1行の値は、

> gala[1,2:5]

Endemics Area Elevation Nearest
Baltra 23 25.09 346 0.6

 である。この値をモデルに代入すると予測値3.69246が得られる。

> coe[1]+coe[-1]%*%t(gala[1,2:5])

Baltra
[1,] 3.69246

 すべての個体に対する予測値は関数predictを用いて呼び出すことができる。

> predict(gala.pm2)

Baltra Bartolome
3.692460 3.538945

 μの予測値は指数変換μ=exp(3.69246)で簡単に得られる。予測されたμの値は関数predictに引数type=”respon”を用いて返すことができる。また関数fiitedを用いて呼び出すことも可能である。

> predict(gala.pm2,type="respon")

Baltra Bartolome
40.14349 34.43059

 予測値の考察の便利ため、モデルgala.pm2による予測結果と観測度数を1つの図にプロットするコマンドとその結果を示す(図1)。

> plot(fitted(gala.pm2),pch=3,xlab="個体番号",ylab="予測値と観測値")
> points(gala$Species,col=2,cex=1.4)
> legend(5,400,c("予測値","観測値"), pch=c(3,1),col=1:2,cex=1.5)

モデルgala.pm2による予測結果と観測度数

図1 予測値と観測値のプロット

 残差の診断図は関数plotで次のように作成することができる(図2)。

> par(mfrow=c(2,2),mar=c(4.5,4.5,2,2))
> plot(gala.pm2)

回帰の残差プロット

図2 回帰の残差プロット

2. 負の2項分布

 2項分布と同じく、事象が起こるか、起こらないかの応答変数について、事象が第k回目に起こる(あるいは成功する)確率を示す次の分布を負の2項分布(Negative binominal distribution)と呼ぶ。

負の2項分布
 式の中のpは事象が起こる(成功する)確率である。負の2項分布を用いた回帰を負の二項回帰(Negative binominal regression)と呼ぶ。
 変換y=z-k,p=(1+p)-1を施すと負の2項分布は次の式で表される。
負の2項分布変換後
 この場合の期待値はμ=θρ、分散はθρ+θρ2=μ+μ2/θである。負の二項回帰には、次に示す対数線形モデルを用いることができる。
負の2項分布対数線形モデル
 負の2項回帰モデルの推測は、関数glmにリンク関数negative.binomialを引数として用いて行うことができる。リンク関数 negative.binomialはパッケージMASSの中に含まれている。パッケージfarawayの中のデータsolderを用いて、負の2項回帰の手順を示す。まずデータを読み込み、その構造を概観する。

> data(solder,package=”faraway”)
> dim(solder) [1] 900 6
> head(solder)

solderのデータ

 データsolderの変数skipsは度数データであるので応答変数として用いることにする。

> library(MASS)
> mod.nb<-glm(skips~.,family= negative.binomial(1),data=solder)
> summary(mod.nb)

mod.nbのsummary

 アソン回帰の場合と同じく、関数summary、coefficients, predict, fittedなどで関連情報を呼び出すことができる。
 回帰係数の項目が多いので、具体的なモデルの書き式は省略する。
 リンク関数negative.binomial(1)に用いた1は、自由に指定することができる。このパラメータをパラメータθと呼ぶことにする。モデルの推測結果はパラメータ に依存する。従って、負の2項回帰分析を行うとき、パラメータθを幾つにし、どのように決めるかが1つの課題である。パッケージMASSの中には、パラメータθを最尤推定して関数glmを用いる関数glm.nbがある。関数glm.nbを用いた例を次に示す。

> mod.gnb<-glm.nb(skips~.,data=solder)
> summary(mod.gnb)
<結果は省略>

3.多項ロジットモデル

 ある実験・観察で、起こり得る結果が m種類であり、それぞれの事象が発生する確率を確率としたとき,n回の試行でそれぞれの事象がyk回起きる確率は次式で定義される多項分布(Multinomial distribution)で表される。m=2の場合は、2項分布となる。

多項分布の式

 応答変数が2値である2項分布の回帰には、ロジットモデルを用いた。同様なアプローチで、多項分布のカテゴリーの中から任意に1つのカテゴリーを選択し、それを基準としてロジットモデルを表現することが可能である。選択されたカテゴリーをベースラインカテゴリー(baseline -category)と呼ぶ。説明の便利のため、最後のカテゴリーをベースラインカテゴリーとした場合のロジットモデルを次に示す。

多項分布の式2

 このモデルを多項ロジットモデル(Multinomial logit model) と呼ぶ。多項ロジットモデルに対応するロジスティック関数を次に示す。

多項ロジットモデル

 多項ロジットモデルのパラメータは、一般的に使用されている尤度を最大化する方法で推測する。多項ロジットモデルの推測には、パッケージneetの中の関数multinom、パッケージVGAMの中の関数vlgmを用いることが可能である。パッケージVGAMはCARNミラーサイトからダウンロードできる。
 花菖蒲のデータirisを用いて、多項ロジスティック回帰の例を示すことにする。データirisの第5列Speciesは3つの異なる種類を示すカテゴリカルデータである。変数Speciesを応答変数、その他を説明変数とし、パッケージVGAMを用いることにする。関数vglmを用いて多項ロジットモデルを推測するためには、リンク関数multinomialを引数として指定することが必要である。

> library(VGAM)
> iris.fit = vglm(Species ~ ., multinomial, data=iris)
> summary(iris.fit)

花菖蒲のデータirisを用いた多項ロジスティック回帰の例

 ここではirisの3つの異なる種類を応答変数としている。返された結果から分かるように、このプログラムでは最後のカテゴリー(3番目)をベースラインとしている。推測された係数を小数点以下1桁まで丸めたモデルを次に示す。式の中の変数SLはSepal.Length、SWはSepal.Width、PLはPetal.Length、PWはPetal.Widthである。また式の中のpkはsummary(iris.fit)の中のmu[,k]に対応する。

花菖蒲のデータirisモデル

 データセット第1行の値を

> iris[1,1:4]

データセット第1行の値

とし、logp1/p3logp2/p3に代入して得られる値は、90.02517、61.73516となる。その計算のコマンドと結果を次に示す。

> (coe<-coefficients(iris.fit))
計算のコマンド
> coe13<- coe[c(1,3,5,7,9)]
> coe23<-coe[c(2,4,6,8,10)]
> coe13[1]+coe13[-1]%*%t(iris[1,1:4])
1
[1,] 90.02517
> coe23[1]+coe23[-1]%*%t(iris[1,1:4])
1
[1,] 61.73516

 このように計算された値は関数predictで呼び出すことが可能である。すべて返すと150行になるので、その一部のみを呼び出して次に示す。データirisの変数Speciesは3つのカテゴリーからなる。1〜50行がsetosa、51〜100行がversicolor、101~150行がvirginicaである。次に3種類の代表として1, 50, 51, 100, 101行の予測値を示す。

> predict(iris.fit)[c(1,50,51,100,101),]

1, 50, 51, 100, 101行の予測値

logp1/p2logp1p3-logp2p3で求めることができる。応答確率は次の式で求まる。
応答確率

 例えば、データirisの第1行の応答確率の予測値はそれぞれ

応答確率の予測値

となる。計算された応答確率は関数fiitedで呼び出すことができる。整数に丸めた2値データ(0,1)を次に示す。

> (temp<-round(fitted(iris.fit)))

setosa versicolor virginica
1 1 0 0
<中略>
50 1 0 0
51 0 1 0
<中略>
100 0 1 0
101 0 0 1
<中略>
150 0 0 1

 結果が正しく予測されているか否かは次のように混同行列を用いると簡単に確認できる。

> fit.tb<-rbind(apply(temp[1:50,],2,sum), apply(temp[51:100,],2,sum),apply(temp[101:150,],2,sum))
> fit.tb

setosa versicolor virginica
[1,] 50 0 0
[2,] 0 49 1
[3,] 0 1 49

 多項ロジットモデルは、非線形重判別法として用いることができる。