ランダム効果入りのベイズ推定モデル

ここでは一般的なGLMMをJAGSで実装するためのスクリプトを紹介します。GLMM、すなわちグループ(ランダム効果)によって切片の値が異なるとするモデルですが、その中身は以下のようなものになります。

SPONSOR LINK

Yi ~ Normal(μi, σ2)
μi = α + β*Xi + Rj(i)

上の式ではRj(i)がランダム効果ですが、これは平均0、分散σR2の正規分布に従う確率変数になります。j(グループID)ごとに正規分布に従う乱数が生成され、その値に応じて切片(グループごとのYの平均値)の値が変化します。平均は0に固定されており、問題となるσR2はデータから推定することになります(グループ間のばらつき)。このようなモデルをJAGSで実装する場合、ランダム効果も数字表記(Numeric)でなければならないので、その点に注意が必要です(lmerなどRの関数では因子型として指定;BUGSでは因子型変数は利用不可)。

データセットを用意します。

> head(D,10)
             x group         y
1   1.89336207     1 3.1029784
2   1.96143183     1 3.0784519
3   1.46765994     1 2.8607315
4  -1.61391038     1 2.1712377
5  -0.57316145     1 2.3753531
6   0.70379683     2 0.6382498
7   2.72372835     2 1.0363066
8   0.97064235     2 0.7026342
9   0.02245384     2 0.6256619
10  1.63627604     2 0.8633311

モデルスクリプトを書きます(model_random.Rとして保存)。ランダム効果は入れ子として表記(R[group[i]])されます。データセットはiをユニットとした繰り返しになるのに対し、ランダム効果はjをユニットとした繰り返しになります。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}

上記のモデルは正規分布になってますが、該当箇所をdpoisやdnegbinとすることで、ポアソンや負の二項分布のモデルでもランダム効果を加えることができます。とてもフレキシブルです。下にいくつかほかの確率分布の場合の例を挙げてみます。モデルを走らせるRコードも下につけますが、詳細はこちらを参照してください。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dpois(mu[i])
    log(mu[i]) <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}
model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  size ~ dnorm(0,ninfo)T(0,1000)

  # Likelihood
  for(i in 1:N){
    y[i] ~ dnegbin(p[i],size)
    p[i] <- size/(mu[i]+size)
    log(mu[i]) <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}
###MCMC setting
n.iter <- 3E+4
n.thin <- max(3, ceiling(n.iter/500))
burn <- ceiling(max(10, n.iter/2))
n.ad <- 100
Sample <- ceiling(n.iter/n.thin)

# data for JAGS
Djags <- list(N = nrow(D), Ng = length(unique(D$group)), y = D$y, x = c(scale(D$x, scale=F)), group=D$group)

# parameters to be monitored
para <- c("b", "sigma", "sigma.R")
  
# set initial values
inits <- replicate(3, list(.RNG.name="base::Wichmann-Hill", .RNG.seed=NA), simplify = FALSE)
for(i in 1:3){ inits[[i]]$.RNG.seed <- 2*i+2 }

# Run JAGS
library(runjags)
m <- read.jagsfile("model_random.R")
post <- run.jags(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)
Spread the love

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です