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
|
Related Posts
Sponsored link