［R］［データ分析］階層ベイズモデルのサンプルコード bayesmパッケージを利用
Rの階層ベイズモデルのサンプルコードが全然見当たらなかったので、自分で書くことにします。詳細を説明しだすとかなり面倒な領域なので、取り合えず使えるというレベルを目指します。
利用するパッケージは「bayesm」です。
階層ベイズに限らずベイズ推定用MCMCの実行はWinBUGSが一般的だと思いますが、Rのみで利用可能かつ事前分布に関する知識なしで利用可能なのが魅力的なので。
階層ベイズモデルについて
階層ベイズモデルは簡単に説明すると個体差を取り入れた統計モデルです。イメージとしては回帰モデルを作成した際の回帰係数が個体ごとに異なっているようなモデルで、最尤法に基づく重回帰モデルやロジスティック回帰モデルより高い表現力を持ちます。
もちろん単純に人ごとに回帰係数を変えるとデータ数より係数の方が多くなり推定できないのですが、係数は個体ごとに大きく異ならないという仮定を入れて問題を解きます。このとき回帰係数自体もあるパラメータに依存しているとし、パラメータ⇒回帰係数⇒予測値という階層関係を持つモデルを構築するのが大雑把な概要です。
階層ベイズモデルの詳細は、例えば以下の文献が参考になります。
http://hosho.ees.hokudai.ac.jp/~kubo/stat/2006/b/kubostat2006b.pdf
http://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/39717/2/kubo2009IEICE.pdf
階層ベイズモデルを使う上で取り合えず以下の2つの式を知っていれば良いと思います。
①個体内モデル
これは分析対象とする個体の反応構造を直接説明するモデルで、hを個体の添え字として個体ごとに通常の回帰モデルが定義されています。Xは説明変数でεは誤差です。
②個体間モデル
こちらは個体ごとの回帰パラメータβを説明変数Zで回帰するモデルになっており、Δが個体によらない共通パラメータ、ηが誤差項になってます。
階層ベイズモデルでは、MCMC法を用いて回帰係数β、共通パラメータΔなどをパラメータの事後確率が最大になるように推定します。ここでインプットとなるのは個体ごとの説明変数XとZですので、モデルを構築する際は基本的にこの2変数を設定することが必要になります。
階層ベイズモデルのサンプルコード
では、階層ベイズ回帰モデルを作ってみたいと思います。利用するパッケージは先ほども述べたように「bayesm」、利用するデータは「LearnBayes」に含まれているbirdextinctを用います。
birdextinctは野鳥の絶滅期間に関するデータで、平均絶滅期間time、つがいの平均数nesting、大きさsize、渡り鳥か否かstatusからなります。
> library(LearnBayes)
> library(bayesm)
>
> # dataのロード
> data(birdextinct)
> attach(birdextinct)
>
> dim(birdextinct)
[1] 62 5
> head(birdextinct)
species time nesting size status
1 Sparrowhawk 3.030 1.000 0 1
2 Buzzard 5.464 2.000 0 1
3 Kestrel 4.098 1.210 0 1
4 Peregrine 1.681 1.125 0 1
5 Grey_partridge 8.850 5.167 0 1
6 Quail 1.493 1.000 0 0
階層ベイズモデルを構築する前に、後ほど比較を行うための通常の重回帰モデルを作成します。ここではtimeの対数をnesting、size、statusを用いて説明する回帰モデルを作成します。
> # 目的変数のlog化
> logtime <- log(time)
> # 重回帰モデルの構築
> fit.lm <- lm(logtime~ nesting + size + status, data=birdextinct)
> summary(fit.lm)
Call:
lm(formula = logtime ~ nesting + size + status, data = birdextinct)
Residuals:
Min 1Q Median 3Q Max
-1.8410 -0.2932 -0.0709 0.2165 2.5167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.43087 0.20706 2.081 0.041870 *
nesting 0.26501 0.03679 7.203 1.33e-09 ***
size -0.65220 0.16667 -3.913 0.000242 ***
status 0.50417 0.18263 2.761 0.007712 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.6524 on 58 degrees of freedom
Multiple R-squared: 0.5982, Adjusted R-squared: 0.5775
F-statistic: 28.79 on 3 and 58 DF, p-value: 1.577e-11 nestingとsizeが絶滅までの時間logtimeに有意に効いていることがわかります。このモデルを元データに適用した場合の予実プロットが以下になります。
> # 実測値と適用値のプロット > pred.lm <- predict(fit.lm,birdextinct[,3:5]) > plot(logtime,pred.lm,xlim=c(0,4),ylim=c(0,4),xlab="実測値",ylab="適用値") > abline(0,1)
実線が実測値と適用値が一致した場合の直線なので、実測値が大きな所で結構ずれていることがわかります。
では、次に同様のデータで階層ベイズモデルを作成します。今回はrhierLinearModelを用いて上記と同様のモデルを作成します。個体差を表すベクトルはZ=1とします。（Z=(1,1,1,1,...) n=hの縦ベクトル)
まず、rhierLinearModelに渡すデータを作成します。
> # 階層ベイズモデルのインプットデータの作成
> regdata <- NULL
> for(i in 1:nrow(birdextinct)){
+ X <- as.matrix(cbind(1,birdextinct[i,3:5]))
+ y <- logtime[i]
+ regdata[[i]] <- list(X=X,y=y)
+ }
> Data <- list(regdata=regdata)Dataに説明変数データが入っており、Zを指定する場合はこの引数にZも渡します。
（Zを渡さなかった場合は、デフォルトでZ=(1,1,1,...)で計算します。）
次にMCMCのパラメータを設定します。
> # MCMCのパラメータ設定 > R <- 10000 # MCMCのイタレーション数 > Mcmc <- list(R=R)
ここではMCMCのイタレーションパラメータのRのみ指定しています。MCMCのパラメータ推定は以下の関数で実行します。
> # MCMCによるパラメータ推定 > out <- rhierLinearModel(Data=Data,Mcmc=Mcmc)
このコマンドによりGibbsサンプリングによるパラメータ推定が行われ、結果がoutに入ります。個々の個体の回帰係数βは以下のコマンドで確認できます。
> # betaのプロット > plot(out$betadraw)
個体によってβの分布が異なるのが確認できると思います。
次に個々の個体のβと全体の平均値を計算します。ここでは2000ステップをbarninとして、2000～10000をサンプリングデータとして利用します。
> # 個々のbetaの平均値を抽出
> beta <- data.frame()
> for (i in 1:nrow(birdextinct)){
+ tmp <- rowMeans(out$betadraw[i,,seq(2000,10000)])
+ beta <- rbind(beta,tmp)
+ }
> colnames(beta) <- c("I",colnames(birdextinct)[3:5])
> head(beta) # 個々のbeta
I nesting size status
1 0.4169747 0.2639782 -0.6531243 0.4630047
2 0.4873114 0.2954235 -0.6862281 0.5263316
3 0.4824536 0.2759155 -0.6713109 0.5169199
4 0.2843932 0.2430212 -0.5800434 0.3086360
5 0.4507167 0.2428417 -0.6514790 0.4898842
6 0.3188920 0.2521075 -0.6259157 0.5079930
> # betaの平均
> apply(beta,2,mean)
I nesting size status
0.4447607 0.2711495 -0.6663089 0.4819408 個体によって回帰係数の平均値が異なっていることが見て取れます。また、全体のβの平均値が重回帰モデルで求めた値と近い値になっていることが確認できます。
次に階層ベイズモデルを元データに適用させて、重回帰モデルと比較します。
> # 元データをモデルに適用
> pred.blm <- c()
> for(i in 1:nrow(birdextinct)){
+ pred.blm <- c(pred.blm,sum(beta[i,] * cbind(1,birdextinct[i,3:5])))
+ }
>
> # 実測値と適用値のプロット
> plot(logtime,pred.lm,xlim=c(0,4),ylim=c(0,4),xlab="実測値",ylab="適用値")
> par(new=T)
> plot(logtime,pred.blm,,xlim=c(0,4),ylim=c(0,4),xlab="",ylab="",col=2)
> abline(0,1)
上図において、重回帰モデルを適用したデータが黒の点、階層ベイズモデルを適用したデータが赤の点で示されています。階層ベイズモデルを用いることによって、全体的に予実差が小さくなっていることがわかります。
次にR^2値を比較します。
> R2.lm <- 1 - sum((logtime - pred.lm)^2) / sum((logtime - mean(logtime))^2) > R2.blm <- 1 - sum((logtime - pred.blm)^2) / sum((logtime - mean(logtime))^2) > R2.lm [1] 0.5982319 > R2.blm [1] 0.9552781
上記のように重回帰モデルで0.6だったR^2値が、階層ベイズモデルを用いると0.96と非常に良くなっています。
ということで、階層ベイズモデルにより個体差を考慮することでモデルの精度が改善されるようです。重回帰モデルを構築した際、結果が予測と大幅にずれるような場合に階層ベイズモデルを用いるのはありなのではないでしょうか。ただし計算コストと相談することにはなりますが。。。
本当はZの指定の仕方（今回は指定なし）や階層ベイズ2項ロジットモデルについても説明したかったのですが、長くなったのでここまでにします。
メモという名の注意点
・階層ベイズモデルの特徴は個体内情報（例えばトランザクション）を取り入れた柔軟なモデルを構築することが可能な点であるので、今回のように個体集約情報のみからモデルを作るのは微妙。
・元データに対する当てはめは良くなるが、どの程度の個体差があるか不明な未知データに対する予測力は重回帰モデルと同程度。（係数として平均値を利用するしかないため）
最後に、今回利用したコードを以下に載せます。
library(LearnBayes)
library(bayesm)
# dataのロード
data(birdextinct)
attach(birdextinct)
dim(birdextinct)
head(birdextinct)
# 目的変数のlog化
logtime <- log(time)
# 重回帰モデルの構築
fit.lm <- lm(logtime~ nesting + size + status, data=birdextinct)
summary(fit.lm)
# 実測値と適用値のプロット
pred.lm <- predict(fit.lm,birdextinct[,3:5])
plot(logtime,pred.lm,xlim=c(0,4),ylim=c(0,4),xlab="実測値",ylab="適用値")
abline(0,1)
# 階層ベイズモデルのインプットデータの作成
regdata <- NULL
for(i in 1:nrow(birdextinct)){
X <- as.matrix(cbind(1,birdextinct[i,3:5]))
y <- logtime[i]
regdata[[i]] <- list(X=X,y=y)
}
Data <- list(regdata=regdata)
# MCMCのパラメータ設定
R <- 10000 # MCMCのイタレーション数
Mcmc <- list(R=R)
# MCMCによるパラメータ推定
out <- rhierLinearModel(Data=Data,Mcmc=Mcmc)
# betaのプロット
plot(out$betadraw)
# 個々のbetaの平均値を抽出
beta <- data.frame()
for (i in 1:nrow(birdextinct)){
tmp <- rowMeans(out$betadraw[i,,seq(2000,10000)])
beta <- rbind(beta,tmp)
}
colnames(beta) <- c("I",colnames(birdextinct)[3:5])
head(beta) # 個々のbeta
# betaの平均
apply(beta,2,mean)
# 元データをモデルに適用
pred.blm <- c()
for(i in 1:nrow(birdextinct)){
pred.blm <- c(pred.blm,sum(beta[i,] * cbind(1,birdextinct[i,3:5])))
}
# 実測値と適用値のプロット
plot(logtime,pred.lm,xlim=c(0,4),ylim=c(0,4),xlab="実測値",ylab="適用値")
par(new=T)
plot(logtime,pred.blm,,xlim=c(0,4),ylim=c(0,4),xlab="",ylab="",col=2)
abline(0,1)
R2.lm <- 1 - sum((logtime - pred.lm)^2) / sum((logtime - mean(logtime))^2)
R2.blm <- 1 - sum((logtime - pred.blm)^2) / sum((logtime - mean(logtime))^2)
R2.lm
R2.blm
