This is an example of running rjags inside this project’s docker container. For more information on Just Another Gibbs Sampler (JAGS), please visit this page.
## Loading required package: coda## Linked to JAGS 4.3.0## Loaded modules: basemod,bugs## $x
## [1] 8 15 22 29 36
##
## $xbar
## [1] 22
##
## $N
## [1] 30
##
## $T
## [1] 5
##
## $Y
## [,1] [,2] [,3] [,4] [,5]
## [1,] 151 141 154 157 132
## [2,] 199 189 200 212 185
## [3,] 246 231 244 259 237
## [4,] 283 275 289 307 286
## [5,] 320 305 325 336 331
## [6,] 145 159 171 152 160
## [7,] 199 201 221 203 207
## [8,] 249 248 270 246 257
## [9,] 293 297 326 286 303
## [10,] 354 338 358 321 345
## [11,] 147 177 163 154 169
## [12,] 214 236 216 205 216
## [13,] 263 285 242 253 261
## [14,] 312 350 281 298 295
## [15,] 328 376 312 334 333
## [16,] 155 134 160 139 157
## [17,] 200 182 207 190 205
## [18,] 237 220 248 225 248
## [19,] 272 260 288 267 289
## [20,] 297 296 324 302 316
## [21,] 135 160 142 146 137
## [22,] 188 208 187 191 180
## [23,] 230 261 234 229 219
## [24,] 280 313 280 272 258
## [25,] 323 352 316 302 291
## [26,] 159 143 156 157 153
## [27,] 210 188 203 211 200
## [28,] 252 220 243 250 244
## [29,] 298 273 283 285 286
## [30,] 331 314 317 323 324The following model definition is taken from this document.
# define model in strings
rats_model_str <- "model {
# likelihood
for (i in 1:n) {
for (j in 1:k) {
Y[i, j] ~ dnorm(alpha + beta[i] * (x[j] - mean(x)), tauy)
}
beta[i] ~ dnorm(mub, taub)
}
## flat prior for (alpha, tauy, mub, taub)
alpha ~ dnorm(0, 0.000001)
tauy ~ dgamma(0.001, 0.001)
mub ~ dnorm(0, 0.000001)
taub ~ dgamma(0.001, 0.001)
}
"
# write model to disk
rats_model_file <- tempfile(pattern = "model_", fileext = ".bug")
cat(rats_model_str, file = rats_model_file)
# print file
readLines(rats_model_file)## [1] "model {"
## [2] " # likelihood"
## [3] " for (i in 1:n) {"
## [4] " for (j in 1:k) {"
## [5] " Y[i, j] ~ dnorm(alpha + beta[i] * (x[j] - mean(x)), tauy)"
## [6] " }"
## [7] " beta[i] ~ dnorm(mub, taub)"
## [8] " }"
## [9] " "
## [10] " ## flat prior for (alpha, tauy, mub, taub)"
## [11] " alpha ~ dnorm(0, 0.000001)"
## [12] " tauy ~ dgamma(0.001, 0.001)"
## [13] " mub ~ dnorm(0, 0.000001)"
## [14] " taub ~ dgamma(0.001, 0.001)"
## [15] "}"rats_model <- jags.model(
file = rats_model_file,
data = list(
"n" = as.integer(rats_data[["N"]]),
"k" = as.integer(rats_data[["T"]]),
"x" = as.numeric(rats_data[["x"]]),
"Y" = as.matrix(rats_data[["Y"]])
)
)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 150
## Unobserved stochastic nodes: 34
## Total graph size: 501
##
## Initializing modelN_ITER <- 10000
rats_sample <- coda.samples(
model = rats_model,
variable.names = c("beta"),
n.iter = N_ITER,
progress.bar = "none"
)##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta[1] -0.09426 0.5108 0.005108 0.04628
## beta[2] -0.09443 0.5067 0.005067 0.04631
## beta[3] -0.09263 0.5121 0.005121 0.04750
## beta[4] -0.08245 0.5107 0.005107 0.04606
## beta[5] -0.07303 0.5266 0.005266 0.04357
## beta[6] -0.08485 0.5161 0.005161 0.04555
## beta[7] -0.08623 0.5069 0.005069 0.04511
## beta[8] -0.08369 0.5185 0.005185 0.04671
## beta[9] -0.08962 0.5072 0.005072 0.04404
## beta[10] -0.09651 0.5138 0.005138 0.04559
## beta[11] -0.08771 0.5120 0.005120 0.04745
## beta[12] -0.09470 0.5117 0.005117 0.04370
## beta[13] -0.10293 0.5119 0.005119 0.04457
## beta[14] -0.10937 0.5045 0.005045 0.04459
## beta[15] -0.09802 0.5118 0.005118 0.04297
## beta[16] -0.08776 0.5109 0.005109 0.04313
## beta[17] -0.09006 0.5130 0.005130 0.04555
## beta[18] -0.08458 0.5103 0.005103 0.04475
## beta[19] -0.08014 0.5110 0.005110 0.04537
## beta[20] -0.08063 0.5127 0.005127 0.04431
## beta[21] -0.08972 0.5177 0.005177 0.04306
## beta[22] -0.09808 0.5232 0.005232 0.04520
## beta[23] -0.09577 0.5173 0.005173 0.04406
## beta[24] -0.10693 0.5081 0.005081 0.04415
## beta[25] -0.10825 0.5117 0.005117 0.04538
## beta[26] -0.08926 0.5132 0.005132 0.04474
## beta[27] -0.08600 0.5120 0.005120 0.04386
## beta[28] -0.09002 0.5168 0.005168 0.04494
## beta[29] -0.09194 0.5046 0.005046 0.04446
## beta[30] -0.09303 0.5123 0.005123 0.04439
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta[1] -1.107 -0.3965 -0.08910 0.2017 0.9364
## beta[2] -1.096 -0.3906 -0.08995 0.2049 0.9210
## beta[3] -1.099 -0.3955 -0.08581 0.1973 0.9329
## beta[4] -1.058 -0.3848 -0.08360 0.2102 0.9592
## beta[5] -1.077 -0.3725 -0.06825 0.2169 0.9953
## beta[6] -1.080 -0.3910 -0.08579 0.2105 0.9760
## beta[7] -1.087 -0.3911 -0.08609 0.2110 0.9295
## beta[8] -1.081 -0.3914 -0.08250 0.2112 0.9661
## beta[9] -1.083 -0.3905 -0.08507 0.1994 0.9242
## beta[10] -1.113 -0.3979 -0.08773 0.2000 0.9347
## beta[11] -1.073 -0.3927 -0.08582 0.2066 0.9331
## beta[12] -1.090 -0.3984 -0.08379 0.2064 0.9098
## beta[13] -1.097 -0.4050 -0.09775 0.1961 0.9181
## beta[14] -1.122 -0.4095 -0.10564 0.1907 0.8867
## beta[15] -1.124 -0.3953 -0.09264 0.1980 0.9130
## beta[16] -1.071 -0.3905 -0.08715 0.2054 0.9302
## beta[17] -1.104 -0.3943 -0.09157 0.2013 0.9401
## beta[18] -1.089 -0.3869 -0.08499 0.2039 0.9718
## beta[19] -1.076 -0.3857 -0.07719 0.2115 0.9334
## beta[20] -1.072 -0.3761 -0.07637 0.2110 0.9582
## beta[21] -1.101 -0.3934 -0.08499 0.2050 0.9601
## beta[22] -1.148 -0.4045 -0.09322 0.2007 0.9519
## beta[23] -1.097 -0.4022 -0.09424 0.1993 0.9458
## beta[24] -1.127 -0.4030 -0.09817 0.1919 0.8980
## beta[25] -1.137 -0.4019 -0.09662 0.1900 0.9072
## beta[26] -1.084 -0.3973 -0.08863 0.2044 0.9360
## beta[27] -1.093 -0.3883 -0.08130 0.2099 0.9256
## beta[28] -1.088 -0.3945 -0.08422 0.2050 0.9495
## beta[29] -1.066 -0.3965 -0.08555 0.1953 0.9101
## beta[30] -1.085 -0.3952 -0.08738 0.2007 0.9396## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.3.5.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rjags_4-10 coda_0.19-3 R2OpenBUGS_3.2-3.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 pillar_1.4.3 compiler_3.6.2 prettyunits_1.0.2
## [5] tools_3.6.2 boot_1.3-23 digest_0.6.23 pkgbuild_1.0.6
## [9] lattice_0.20-38 evaluate_0.14 lifecycle_0.1.0 tibble_2.1.3
## [13] gtable_0.3.0 pkgconfig_2.0.3 rlang_0.4.2 cli_2.0.0
## [17] yaml_2.2.0 parallel_3.6.2 xfun_0.11 loo_2.2.0
## [21] gridExtra_2.3 stringr_1.4.0 knitr_1.26 dplyr_0.8.3
## [25] stats4_3.6.2 grid_3.6.2 tidyselect_0.2.5 glue_1.3.1
## [29] inline_0.3.15 R6_2.4.1 processx_3.4.1 fansi_0.4.0
## [33] rmarkdown_2.0 rstan_2.19.2 ggplot2_3.2.1 callr_3.4.0
## [37] purrr_0.3.3 magrittr_1.5 htmltools_0.4.0 scales_1.1.0
## [41] ps_1.3.0 StanHeaders_2.19.0 matrixStats_0.55.0 assertthat_0.2.1
## [45] colorspace_1.4-1 stringi_1.4.3 lazyeval_0.2.2 munsell_0.5.0
## [49] crayon_1.3.4