Summary
JAGSによるVarying Intercept and Slope Modelの書き方の備忘録。
Sponsored link
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
|
|
![[商品価格に関しましては、リンクが作成された時点と現時点で情報が変更されている場合がございます。]](https://hbb.afl.rakuten.co.jp/hgb/1797f0a4.6e20d489.1797f0a5.7bc9b355/?me_id=1346503&item_id=11420248&m=https%3A%2F%2Fthumbnail.image.rakuten.co.jp%2F%400_mall%2Fexcellentimports%2Fcabinet%2F06293380%2F200019952.jpg%3F_ex%3D80x80&pc=https%3A%2F%2Fthumbnail.image.rakuten.co.jp%2F%400_mall%2Fexcellentimports%2Fcabinet%2F06293380%2F200019952.jpg%3F_ex%3D240x240&s=240x240&t=picttext "[商品価格に関しましては、リンクが作成された時点と現時点で情報が変更されている場合がございます。]")
Sponsored link