BLOG

BUGSファイルの記述

まずはBUGS言語による統計モデルの記述から説明を進めます。BUGSファイルは「統計モデルの構造」のみを記述するファイルです。ここでは、森林パッチiの鳥の個体数Yiを、森林面積Xiで説明するポアソン線形モデルを例に挙げて考えます。

Eq. 1)
Yi ~ Pois(λi)
log(λi) = α + β*Xi

上記のモデルは、RのGLMでは次のように記述されます。

###R code###
data <- read.csv(“data.csv”)
glm(y ~ X, data, family=poisson)

ここで気を付けなければならないのは、

  • GLMのようにRでパッケージ化されているものだと、eq.1で書かれているほとんど部分が自動化されており、リンク関数(log)などを自分で考える必要がない
  • 特に指定しなくとも、サンプルの繰り返しが“i”として認識されている

ことです。これら2点がBUGSとR言語では大きく異なります。このことを念頭におき、上記のモデルをBUGS言語で書いてみます。

###BUGS code###
model{
#Likelihood
 for(i in 1:N.sample){ # loop for i
   y[i] ~ dpois(lambda[i]) # poisson model
   log(lambda[i]) <- alpha + beta * X[i] # linear model structure
 }

#Prior
 alpha ~ dnorm(0, 0.001)
 beta ~ dnorm(0, 0.001)
}

model{}で囲われている範囲がモデル構造を指定するエリアになります({}の外にテキストを書くとエラーがでるので注意)。「~」は確率論的ノード(Stochastic node)と呼ばれ、左辺が確率変数である場合に使われます。確率論的ノードの右辺には、必ず確率分布の記述が必要です(上記ではdpoisが確率分布、Poisson distributionを指すもの)。「<-」は決定論的ノード(Deterministic node)と呼ばれ、=(イコール)と対応します。こうしてみると、BUGS言語ではEq.1をすべて自分で書き下す必要があると分かります。 Prior以降の記述は「事前分布」と呼ばれるものです。ベイズ統計の枠組みでは、すべてパラメーターに対し、事前分布を設定しなければなりません(事前分布を設定し忘れるとエラーが出る)。ここでは切片αと傾きβに対し、事前の情報は何も得られないので、「無情報事前分布」というものを設定しています(ほとんどケースはこれにあたる)。ここでは、平均0、分散1000の正規分布を事前分布として設定することで、切片αと傾きβは「かなり広い範囲の値をとる可能性がある」としています。「dnorm」が正規分布であることを示しており、一つ目の引数が「平均」、二つ目の引数が「精度(分散の逆数=1/1000)」を表します。ここまで書き終えたら、ファイルを.Rもしくは.txtファイルとして保存します(ここでは“model_PoissonRegression.R”として保存することにします)。続きはこちら

SPONSOR LINK

BUGSの使い方

最近ではRパッケージも実装され、容易にRからBUGSソフトウェアを使うことができます。その一連の流れは以下のようになっています。

  • 扱いたい統計モデルをBUGS言語で記述し、.Rもしくは.txtファイルとして保存
  • 解析に使うデータをRで準備(データ整形のしやすさはRのほうが優れているため)
  • RからBUGSソフトウェアにデータを渡し、指定した統計モデルのもとで計算を実行
  • BUGSからRに渡される計算結果を保存

この際、「どんな統計モデルを計算させるのか」という部分だけ、RではなくBUGS言語で記述する必要があります。したがって、統計モデルの実装の際には、RおよびBUGSスクリプトをそれぞれ別個に用意することになります。

SPONSOR LINK

ベイズ統計とBUGS

ベイズ統計、と聞くとなにやらすべてが新しいものに聞こえますが、扱うモデルの「構造」自体は大きく変わりません(線形回帰、ANOVAなど)。ただし、回帰係数のような“パラメーター”の扱いと推定手法が大きく異なっており、ベイズ統計の実装にはマルコフ連鎖モンテカルロのような逐次的な推定法が必要です。自前でアルゴリズムを構築することも可能ですが、そんな煩雑なプログラムを自分で書くというのは現実的な選択肢ではないでしょう。

そこで、BUGSソフトウェアの登場です。BUGSとはBayesian inference Using Gibbs Samplerの略で、ギブスサンプリングによるベイズ統計の実装を指します。このギブスサンプリングを実行してくれるのがWinBUGSをはじめとするBUGSソフトウェアで、Rから動かすことができます。そのため、昨今のベイズ統計を用いた研究では、RとBUGSソフトウェアの連携が一般的な手法になっています。なお、WinBUGSのほかにも、OpenBUGSとJAGS(Stanもベイズ実装のためのソフトですが、ギブスサンプリングではない)があります。WinBUGSは開発が終了していますが、OpenBUGSとJAGSは最新バージョンが随時ネット上でダウンロードできるようになっており、計算効率もそれらのほうがより良いものとなりつつあります。ここでは、BUGの中ではもっともクラッシュしにくいJAGS(Just Another Gibbs SamplerなのでBUGSというには語弊がありますが。。。)について解説します。

JAGSのダウンロードは次のページから可能です(JAGS)。ダウンロードはProgram files(デフォルト)もしくはCドライブ直下がお勧めです。