How to use JAGS: Simple Linear Regression

Summary

JAGSの基本的な使い方の備忘録。ここでは、もっとも単純な単回帰モデルをJAGSで実装する。実装の手順は大きく二つに分けられる。(1)JAGSモデルの記述:JAGS言語を使ってJAGSのモデルファイルを記述する。(2)Rによる実装:RからJAGSのモデルファイルを走らせる。

Sponsored link

JAGSモデルの記述

サンプルデータ

ここではもっとも単純な例として、(GLMでも分析できる)単回帰モデルの分析をJAGSでする。ここで使用するサンプルデータは次のようになっている。

head(dat1)
##           Y        X
## 1 11.859762 6.635249
## 2  4.682198 2.290394
## 3 14.809464 7.288998
## 4  7.549894 3.507948
## 5  8.841815 3.150085
## 6 13.473219 8.793556
plot(Y ~ X, dat1)

plot of chunk code1

このデータに対し、応答変数\(y_i\)を説明変数\(x_i\)で説明する簡単な統計モデルを考える(添え字\(i\)はサンプルのID番号)。応答変数\(y_i\)の誤差は正規分布に従うと仮定し、以下の統計モデルを考える。
\[
y_i \sim Normal(\mu_i, \sigma^2)\\
\mu_i = \beta_1 + \beta_2 x_i
\]
ここで、\(\beta_1\)は切片(intercept)を、\(\beta_2\)は回帰係数(slope)を表しており、これらが推定対象となる主なパラメーターである。

JAGSの記法

上の統計モデルをJAGS言語で書き下すと以下のようになる。

# JAGS script
model{
  ninfo <- 0.0001

  ## Liklihood
  for(i in 1:Nsample){
    Y[i] ~ dnorm(mu[i], tau) # dnorm(mu[i], tau): normal distribution with a mean mu[i] and precision tau
    mu[i] <- b[1] + b[2]*X[i]
  }

  ## Prior
  ### prior distribution for residual variance
  tau <- pow(sigma, -2) # tau = 1/sigma^2
  sigma ~ dunif(0, 100) # dunif(0, 100): uniform distribution with a range of 0 to 100

  ### prior distribution for parameters
  for(j in 1:2){ b[j] ~ dnorm(0, ninfo) }
}

Liklihoodとされたブロックでは、上に述べたモデル式がそのまま書き下されている。JAGS言語ではベクター化した記法が実装されていないので、for構文を多用する。

Priorとされたブロックではパラメーターの事前分布を指定している。それぞれのパラメーターには、「とりうる値」の範囲を適切に表現できるよう事前分布が与えなければならない。上の例では、切片b[1]および傾きb[2]に対しては平均0、分散100の正規分布が仮定されている。これらのパラメーターはあらゆる実数をとりうるため、このような事前分布が指定されている。一方、標準偏差sigmaは負の値をとらないため、0から100の範囲の一様分布が仮定されている。

なお、頭文字が大文字の変数と小文字の変数は意図的に分けられている。大文字のものは「データ(Y, X)」である。一方、頭文字が小文字のものは「パラメーター(\(\sigma\), \(\beta_1\), \(\beta_2\))」であり、JAGSが推定する未知の定数である。この記法に従わないからといってエラーがでるわけではないが、解析者の頭を整理するという意味において有効な記法である。

ここまで書いたら、モデルファイルをRプロジェクト内の好きな場所に保存する。ここでは、ファイル名をmodel1.Rとし、作業フォルダ内のサブディレクトリBayesModelに保存する。

Rによる実装

データの読み込み

RからJAGSを起動し、model1.Rに記述したモデルのパラメーターを推定する。まず、以下の手順に沿ってデータを読み込む。

  • データの読み込み(ブロック"Read data")
  • データフレームのリスト化(ブロック"Dataframe for JAGS")
  • モニターするパラメーター名の宣言(ブロック"Claim parameters")
# Read data
Y <- dat1$Y
X <- c(scale(dat1$X) ) # scaling X with a mean 0 and SD 1

# Dataframe for JAGS
Djags <- list( Y = Y, X = X, Nsample = nrow(dat1))

# Claim parameters to be monitored
para <- c("sigma", "b")

この際、リストDjagsおよびベクターparaで与える変数名は、JAGSのモデルファイルの中で使われている変数名と一致するようにする。

パラメーター推定

ライブラリ"runjags"に入っている関数"run.jags"を使い、以下の手順でRからJAGSによるパラメーター推定を実行する。

  • MCMCの試行数、焼き捨て試行数、間引き数、サンプルするMCMC数を指定(ブロック"MCMC setting")
  • MCMCの初期値の指定(ブロック"Initial values")
  • JAGSによるパラメ―ター推定(ブロック"Run JAGS")
# MCMC setting
n.ad <- 100 # number of adaptation
n.iter <- 3E+3 # number of total iteration
n.thin <- max(3, ceiling(n.iter/500)) # number of thinning
burn <- ceiling(max(10, n.iter/2)) # number of burn-in
Sample <- ceiling(n.iter/n.thin) # number of sample per chain

# Intial values
inits <- replicate(3, list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = NA), simplify = FALSE)
for(i in 1:3){ inits[[i]]$.RNG.seed <- i }

# Run JAGS
library(runjags)
m <- read.jagsfile("BayesModel/model1.R") # read model file
post <- run.jags(model = m$model, monitor = para, data = Djags,
                 n.chains = 3, inits =inits, method = "parallel",
                 burnin = burn, sample = Sample, adapt = n.ad, thin = n.thin,
                 n.sims = 3)
## Calling 3 simulations using the parallel method...
## Following the progress of chain 1 (the program will wait for all
## chains to finish before continuing):
## Welcome to JAGS 4.1.0 on Tue Apr 09 23:13:36 2019
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 3
##    Total graph size: 414
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 100
## -------------------------------------------------| 100
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 1500
## -------------------------------------------------| 1500
## ************************************************** 100%
## . . . Updating 3000
## -------------------------------------------------| 3000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## . 
## All chains have finished
## Simulation complete.  Reading coda files...
## Coda files loaded successfully
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 3 variables....
## Finished running the simulation
summary(post)
##       Lower95   Median Upper95     Mean        SD Mode       MCerr MC%ofSD
## sigma 1.70331 1.956005 2.24335 1.964718 0.1403403   NA 0.003623571     2.6
## b[1]  8.33519 8.732565 9.11145 8.732600 0.2007176   NA 0.005182507     2.6
## b[2]  3.51754 3.886605 4.27558 3.890355 0.1940865   NA 0.005012817     2.6
##       SSeff         AC.60     psrf
## sigma  1500 -0.0145773216 1.000176
## b[1]   1500  0.0004248646 1.000246
## b[2]   1499 -0.0332069016 1.000495

function “run.jags”

run.jags関数の引数の説明。

  • model: モデルファイル(model1.R)をread.jagsfile関数を使って読み込んだもの
  • monitor: モニターするパラメーター名のベクター
  • data: モデルで使用するデータのリスト
  • n.chains: MCMCの鎖の数(デフォルトの3でよい)
  • inits: MCMCの初期値
  • method: 推定計算方法(parallelとすると並列計算を実行)
  • burnin: MCMCの焼き捨て試行数
  • sample: 事後分布の推定に用いるMCMC試行数(1鎖あたり)
  • adapt: 焼き捨て試行前の試運転MCMC試行数(1鎖あたり)
  • thin: 間引きするMCMC試行数(1鎖あたり)
  • n.sims: 並列試行の数(MCMCの鎖の数に合わせる)

Leave a Reply

Your email address will not be published. Required fields are marked *