Varying Intercept Model

Summary

JAGSによるVarying Intercept Modelの書き方の備忘録。

Sponsored link

Varying Intercept Model

切片がグループごとに変化するモデルを指す。書き下すと、次のような数式になる。

\[
y_i \sim Normal(\mu_i, \sigma^2)\\
\mu_i = \beta_1 + \beta_2x_i + R_{j(i)}\\
R_j \sim Normal(0, \sigma_R^2)
\]

切片である\(\beta_1\)に対してグループ単位\(j\)で変化する変数\(R_j\)が加わった構造になっている。これにより、グループごとに切片の値が異なる(グループごとに応答変数\(y_i\)の平均水準が異なる)ことを表現する。意味としては同じだが、次のように書くこともできる。

\[
y_i \sim Normal(\mu_i, \sigma^2)\\
\mu_i = \beta_{1,j(i)} + \beta_2x_i\\
\beta_{1,j} \sim Normal(\mu_{\beta_1}, \sigma_{\beta_1}^2)
\]

ここで、切片のグループ間のばらつきを表す変数を、正規分布に従う確率変数とする。こうすることで、グループごとにそれぞれ切片を推定するのではなく、全体切片(\(\beta_1\)もしくは\(\mu_{\beta_1}\))と分散(\(\sigma_R^2\)もしくは\(\sigma_{\beta_1}^2\))という二つのパラメーターの推定で済ませることができる(グループ数が多くなるほどパラメーター数を節約できる)。ここでは次のテストデータに、このVarying Intercept Model(この場合GLMMと同じ)を当てはめてみる。

head(dat, 15)
##             Y         X groupID
## 1   0.7620270 5.1462547       1
## 2   0.6752850 4.3609549       1
## 3   0.6933545 4.4518697       1
## 4   1.4239270 7.7826422       1
## 5   1.5188387 7.8389356       1
## 6  -0.5569306 0.3245126       1
## 7   0.1546029 2.9847637       1
## 8   0.4130097 3.4364609       1
## 9   0.5041910 4.3507010       1
## 10  1.1575450 6.4613893       1
## 11  2.0210685 7.4991161       2
## 12  1.4163438 5.5429805       2
## 13  1.8921774 7.1133113       2
## 14  1.4546292 5.5301212       2
## 15  2.0675598 6.9715206       2
plot(Y ~ X, dat, pch = 21, col = NA, bg = COL[factor(dat$groupID)])

plot of chunk unnamed-chunk-1

JAGS script

JAGSの記法で上のモデルを書き下す。

model{
  ninfo <- 0.001

  ## Likelihood
  for(i in 1:Nsample){
    Y[i] ~ dnorm(mu[i], tau[1])
    mu[i] <- b[1] + b[2]*X[i] + R[groupID[i]]
  }

  for(j in 1:Ngroup){
    R[j] ~ dnorm(0, tau[2])
  }

  ## Prior
  ### prior distribution for residual and random-effect variance
  for(k in 1:2){
    tau[k] <- pow(sigma[k], -2)
    sigma[k] ~ dnorm(0, ninfo)T(0, 50)
  }

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

RからJAGSを走らせる。

# R script
# Read data
Y <- dat$Y
X <- c(scale(dat$X) ) # scaling X with a mean 0 and SD 1
groupID <- dat$groupID # must be numeric

# Dataframe for JAGS
Djags <- list( Y = Y, X = X, groupID = groupID,
               Nsample = nrow(dat), Ngroup = length(unique(groupID) ) )

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

# MCMC setting
n.ad <- 100 # number of adaptation
n.iter <- 3E+4 # 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/model_VIM.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.3.0 on Sun May 12 20:27:24 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: 14
##    Total graph size: 523
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 100
## -------------------------------------------------| 100
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 15000
## -------------------------------------------------| 15000
## ************************************************** 100%
## . . . Updating 30000
## -------------------------------------------------| 30000
## ************************************************** 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 4 variables....
## Finished running the simulation
summary(post)
##           Lower95    Median  Upper95      Mean          SD Mode
## sigma[1] 0.091691 0.1055765 0.122292 0.1061146 0.008165545   NA
## sigma[2] 0.353436 0.6218650 1.049360 0.6591029 0.194054183   NA
## b[1]     1.027640 1.4698300 1.931890 1.4691757 0.219905045   NA
## b[2]     0.812320 0.8366380 0.857253 0.8365962 0.011574117   NA
##                 MCerr MC%ofSD SSeff       AC.600      psrf
## sigma[1] 0.0002108335     2.6  1500  0.013629191 1.0002712
## sigma[2] 0.0069789269     3.6   773  0.015228704 1.0064985
## b[1]     0.0188719233     8.6   136  0.190227003 1.0261817
## b[2]     0.0003065785     2.6  1425 -0.008032277 0.9995096

sigma[1]は残差のSDを、sigma[2]はグループ間のばらつきのSDを表している。

Reference

Leave a Reply

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