Summary

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

Varying Intercept and Slope Model

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

\[
y_i \sim Normal(\mu_i, \sigma^2)\\
\mu_i = \beta_{j(i),1} + \beta_{j(i),2}x_i\\
\beta_{k,j} \sim MVN(\mu_{\beta_k}, \sum)
\]

切片および偏回帰係数である\(\beta_k\)がグループ単位\(j\)で変化すると仮定している。これにより、グループごとに切片および偏回帰係数の値が異なる（グループごとに応答変数\(y_i\)の平均水準と説明変数\(x_i\)の効果が異なる）ことを表現する。切片と偏回帰係数には多変量正規分布を事前分布とすることで、切片と偏回帰係数の相関も考慮する。ここでは次のテストデータに、このVarying Intercept and Slope Modelを当てはめてみる。

head(dat, 15)

##              Y        X groupID
## 1    1.6533548 3.715236       1
## 2    1.9452796 1.334316       1
## 3    1.0739305 7.079826       1
## 4   -0.9164404 9.527487       1
## 5    1.1134740 5.339768       1
## 6    1.3665161 5.336097       1
## 7    0.8860018 7.409342       1
## 8    0.7536042 6.252334       1
## 9    0.4763039 7.106482       1
## 10   1.5693835 4.488241       1
## 11  -8.5185040 4.881332       2
## 12 -15.3315302 9.889448       2
## 13  -8.4777794 5.365884       2
## 14  -6.9886498 4.510636       2
## 15  -8.9329260 5.492908       2

plot(Y ~ X, dat, pch = 21, col = NA, bg = COL[factor(dat$groupID)])

JAGS script

JAGSの記法で上のモデルを書き下す。多変量正規分布の分散共分散行列の事前分布には、逆ウィシャート分布を使う（Gelman and Hill 2007）。

model{
  ninfo <- 0.001

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

  for(j in 1:Ngroup){
      b[j,1:2] ~ dmnorm(mu_b[], TAU[,])
  }

  ## Prior
  ### prior distribution for residual
  tau <- pow(sigma, -2)
  sigma ~ dnorm(0, ninfo)T(0, 50)

  ### prior distribution for the variance-covariance matrix
  df <- 3 #df set to be "the number of parameters + 1"
  TAU[1:2,1:2] ~ dwish(W[,], df)
  SIGMA[1:2,1:2] <- inverse(TAU[,]) # SIGMA: variance-covariance matrix

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

RからJAGSを走らせる。

# R script
# Read data
# 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) ), W = diag(2) )

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

# 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_VISM.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:17:57 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: 550
## . 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 27 variables....
## Note: Unable to calculate the multivariate psrf
## Finished running the simulation

summary(post)

##               Lower95     Median   Upper95       Mean          SD Mode
## sigma       0.3751710  0.4404170  0.513176  0.4419768  0.03558219   NA
## mu_b[1]    -0.2577750  3.1641050  6.359250  3.1739754  1.70234159   NA
## mu_b[2]     0.0168568  1.7536100  3.697590  1.7594393  0.94669693   NA
## b[1,1]      0.9039170  1.1868550  1.452660  1.1845636  0.14414173   NA
## b[2,1]     -8.8103600 -8.5556750 -8.264360 -8.5552072  0.14109001   NA
## b[3,1]     -3.3198400 -2.9614350 -2.548480 -2.9587953  0.19426601   NA
## b[4,1]      5.4346200  5.6953950  5.971060  5.6946375  0.13866264   NA
## b[5,1]      4.4925200  4.7560500  5.058430  4.7549550  0.14239388   NA
## b[6,1]      5.7197000  6.0183900  6.288960  6.0173780  0.14507885   NA
## b[7,1]      2.2830700  2.5403650  2.832670  2.5394553  0.13927539   NA
## b[8,1]      6.1127500  6.3897700  6.699650  6.3934041  0.14455082   NA
## b[9,1]      9.2194300  9.5499750  9.840590  9.5499131  0.15528920   NA
## b[10,1]     6.5902000  6.8513250  7.121900  6.8459193  0.14074841   NA
## b[1,2]     -1.0614200 -0.7489370 -0.419042 -0.7453538  0.16837561   NA
## b[2,2]     -4.1007500 -3.7976550 -3.510640 -3.7956256  0.15268131   NA
## b[3,2]     -2.2949100 -1.8071200 -1.437390 -1.8149141  0.21748971   NA
## b[4,2]      2.8357300  3.1170800  3.406740  3.1114430  0.14635700   NA
## b[5,2]      3.0335100  3.2890450  3.550790  3.2880859  0.13166068   NA
## b[6,2]      3.3216000  3.5742500  3.821390  3.5759763  0.12517192   NA
## b[7,2]      0.5322530  0.8220305  1.115360  0.8253910  0.14718849   NA
## b[8,2]      3.9578600  4.3208950  4.620800  4.3215943  0.17057382   NA
## b[9,2]      5.1243500  5.3899700  5.675110  5.3897714  0.14032312   NA
## b[10,2]     3.0408600  3.2601250  3.497570  3.2599071  0.11711808   NA
## SIGMA[1,1] 10.5308000 25.2221000 60.410900 29.5145160 17.00228578   NA
## SIGMA[2,1]  4.6528500 13.5175000 30.932600 15.7947728  9.18789065   NA
## SIGMA[1,2]  4.6528500 13.5175000 30.932600 15.7947728  9.18789065   NA
## SIGMA[2,2]  2.7478300  7.8835750 17.601200  9.1295840  5.19271360   NA
##                   MCerr MC%ofSD SSeff       AC.600      psrf
## sigma      0.0009187281     2.6  1500  0.028654527 1.0014300
## mu_b[1]    0.0428430920     2.5  1579 -0.001540665 0.9997240
## mu_b[2]    0.0236876118     2.5  1597 -0.005380042 1.0004057
## b[1,1]     0.0040046288     2.8  1296  0.036204152 0.9995522
## b[2,1]     0.0037893253     2.7  1386  0.046546319 1.0013263
## b[3,1]     0.0050159268     2.6  1500 -0.008898043 1.0010946
## b[4,1]     0.0035802541     2.6  1500 -0.028980985 1.0001801
## b[5,1]     0.0036765943     2.6  1500 -0.013445106 1.0008722
## b[6,1]     0.0038762431     2.7  1401 -0.033237452 1.0048720
## b[7,1]     0.0035960751     2.6  1500  0.015975006 1.0003557
## b[8,1]     0.0039621617     2.7  1331  0.003186194 0.9998950
## b[9,1]     0.0041292621     2.7  1414 -0.014296008 1.0002408
## b[10,1]    0.0035827351     2.5  1543  0.010655330 1.0024681
## b[1,2]     0.0041429637     2.5  1652 -0.021483000 1.0037192
## b[2,2]     0.0036852979     2.4  1716 -0.010208375 1.0003075
## b[3,2]     0.0056155602     2.6  1500 -0.013438513 1.0027741
## b[4,2]     0.0036803256     2.5  1581 -0.010301315 1.0013358
## b[5,2]     0.0036376808     2.8  1310 -0.032447464 1.0009717
## b[6,2]     0.0030396787     2.4  1696  0.021798710 0.9999145
## b[7,2]     0.0038728773     2.6  1444  0.059380598 1.0003314
## b[8,2]     0.0047258325     2.8  1303 -0.026214046 1.0019232
## b[9,2]     0.0037061977     2.6  1434 -0.024624042 1.0008715
## b[10,2]    0.0030896924     2.6  1437 -0.027386486 1.0005015
## SIGMA[1,1] 0.4302102593     2.5  1562 -0.006253852 0.9995675
## SIGMA[2,1] 0.2446047599     2.7  1411 -0.011407742 1.0002483
## SIGMA[1,2] 0.2446047599     2.7  1411 -0.011407742 1.0002483
## SIGMA[2,2] 0.1388365512     2.7  1399 -0.012689870 1.0007052

パラメータb[j,1]はグループjの切片の推定値、b[j,2]は偏回帰係数の推定値を示している。図からわかるように、グループによって切片、偏回帰係数が大きく異なっている。

Reference

Data Analysis Using Regression and Multilevel / Hierarchical Models (Analytical Methods for Social Research)【楽天海外直送】【英語の本】【洋書】 価格：7501円（税込、送料無料) (2019/5/12時点) 楽天で購入