OpenBUGS in R
Rendered: Dec 31, 2019
OpenBUGS
This is an example of running R2OpenBUGS inside this project’s docker container. For more information on Bayesian inference Using Gibbs Sampling (BUGS), please visit this page.
Example
Library
Data
## $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 324Model Definition
The following model definition is taken from OpenBUGS tutorial.
# define model in strings
rats_model <- function() {
# likelihood
for( i in 1 : N ) {
for( j in 1 : T ) {
Y[i , j] ~ dnorm(mu[i , j],tau.c)
mu[i , j] <- alpha[i] + beta[i] * (x[j] - xbar)
# deleted unnecessary lines
}
alpha[i] ~ dnorm(alpha.c,alpha.tau)
beta[i] ~ dnorm(beta.c,beta.tau)
}
# prior
tau.c ~ dgamma(0.001,0.001)
sigma <- 1 / sqrt(tau.c)
alpha.c ~ dnorm(0.0,1.0E-6)
alpha.tau ~ dgamma(0.001,0.001)
beta.c ~ dnorm(0.0,1.0E-6)
beta.tau ~ dgamma(0.001,0.001)
alpha0 <- alpha.c - xbar * beta.c
}
# write model to disk
# model.file must be renamed with .txt rather than .bug
rats_model_file <- tempfile(pattern = "model_", fileext = ".txt")
write.model(rats_model, con = rats_model_file)
# print file
readLines(rats_model_file)## [1] "model"
## [2] "{"
## [3] " for (i in 1:N) {"
## [4] " for (j in 1:T) {"
## [5] " Y[i, j] ~ dnorm(mu[i, j], tau.c)"
## [6] " mu[i, j] <- alpha[i] + beta[i] * (x[j] - xbar)"
## [7] " }"
## [8] " alpha[i] ~ dnorm(alpha.c, alpha.tau)"
## [9] " beta[i] ~ dnorm(beta.c, beta.tau)"
## [10] " }"
## [11] " tau.c ~ dgamma(0.001, 0.001)"
## [12] " sigma <- 1/sqrt(tau.c)"
## [13] " alpha.c ~ dnorm(0.00000E+00, 1.00000E-06)"
## [14] " alpha.tau ~ dgamma(0.001, 0.001)"
## [15] " beta.c ~ dnorm(0.00000E+00, 1.00000E-06)"
## [16] " beta.tau ~ dgamma(0.001, 0.001)"
## [17] " alpha0 <- alpha.c - xbar * beta.c"
## [18] "}"Model Run
rats_fit <- bugs(
data = rats_data,
inits = list(
# source: http://www.openbugs.net/Examples/Ratsinits.html
list(
"alpha" = rep(250, times = rats_data[["N"]]),
"beta" = rep(6, times = rats_data[["N"]]),
"alpha.c" = 150,
"beta.c" = 10,
"tau.c" = 1,
"alpha.tau" = 1,
"beta.tau" = 1
)
),
model.file = rats_model_file,
parameters.to.save = c("alpha0", "beta.c", "sigma"),
# MCMC parameters to be consistent with
# http://www.openbugs.net/Examples/Rats.html
n.chains = 1,
n.burnin = 1000,
n.iter = 1000 + 10000,
# other params
digits = 9,
bugs.seed = sample(x = seq_len(14), size = 1) # getting some randomness
)## Inference for Bugs model at "/tmp/RtmpVlItOc/model_6173abbd8.txt",
## Current: 1 chains, each with 11000 iterations (first 1000 discarded)
## Cumulative: n.sims = 10000 iterations saved
## mean sd 2.5% 25% 50% 75% 97.5%
## alpha0 245.4 12.2 221.2 237.2 245.3 253.7 269.7
## beta.c -0.1 0.1 -0.3 -0.2 -0.1 0.0 0.1
## sigma 13.8 0.9 12.1 13.1 13.7 14.4 15.7
## deviance 1211.8 9.6 1195.0 1205.0 1211.0 1218.0 1233.0
##
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 34.0 and DIC = 1246.0
## DIC is an estimate of expected predictive error (lower deviance is better).Comparing to the results from OpenBUGS tutorial:
Note that the model did not specify constraints on alpha and beta’s hence creating identifiability issues.
Session Info
## 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] 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] coda_0.19-3 gridExtra_2.3 stringr_1.4.0 knitr_1.26
## [25] dplyr_0.8.3 stats4_3.6.2 grid_3.6.2 tidyselect_0.2.5
## [29] glue_1.3.1 inline_0.3.15 R6_2.4.1 processx_3.4.1
## [33] fansi_0.4.0 rmarkdown_2.0 rstan_2.19.2 ggplot2_3.2.1
## [37] callr_3.4.0 purrr_0.3.3 magrittr_1.5 htmltools_0.4.0
## [41] scales_1.1.0 ps_1.3.0 StanHeaders_2.19.0 matrixStats_0.55.0
## [45] assertthat_0.2.1 colorspace_1.4-1 stringi_1.4.3 lazyeval_0.2.2
## [49] munsell_0.5.0 crayon_1.3.4