jags: Run 'JAGS' from R

Description Usage Arguments Details Author(s) References Examples

View source: R/jags.R

Description

The jags function takes data and starting values as input. It automatically writes a jags script, calls the model, and saves the simulations for easy access in R.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
jags(data, inits, parameters.to.save, model.file="model.bug",
  n.chains=3, n.iter=2000, n.burnin=floor(n.iter/2),
  n.thin=max(1, floor((n.iter - n.burnin) / 1000)),
  DIC=TRUE, working.directory=NULL, jags.seed = 123,
  refresh = n.iter/50, progress.bar = "text", digits=5,
  RNGname = c("Wichmann-Hill", "Marsaglia-Multicarry",
              "Super-Duper", "Mersenne-Twister"),
  jags.module = c("glm","dic")
  )

jags.parallel(data, inits, parameters.to.save, model.file = "model.bug",
             n.chains = 2, n.iter = 2000, n.burnin = floor(n.iter/2),
             n.thin = max(1, floor((n.iter - n.burnin)/1000)),
             n.cluster= n.chains, DIC = TRUE,
             working.directory = NULL, jags.seed = 123, digits=5,
             RNGname = c("Wichmann-Hill", "Marsaglia-Multicarry",
              "Super-Duper", "Mersenne-Twister"),
             jags.module = c("glm","dic"),
             export_obj_names=NULL,
             envir = .GlobalEnv
             )

jags2(data, inits, parameters.to.save, model.file="model.bug",
  n.chains=3, n.iter=2000, n.burnin=floor(n.iter/2),
  n.thin=max(1, floor((n.iter - n.burnin) / 1000)),
  DIC=TRUE, jags.path="",
  working.directory=NULL, clearWD=TRUE,
  refresh = n.iter/50)

Arguments

data

(1) a vector or list of the names of the data objects used by the model, (2) a (named) list of the data objects themselves, or (3) the name of a "dump" format file containing the data objects, which must end in ".txt", see example below for details.

inits

a list with n.chains elements; each element of the list is itself a list of starting values for the BUGS model, or a function creating (possibly random) initial values. If inits is NULL, JAGS will generate initial values for parameters.

parameters.to.save

character vector of the names of the parameters to save which should be monitored.

model.file

file containing the model written in BUGS code. Alternatively, as in R2WinBUGS, model.file can be an R function that contains a BUGS model that is written to a temporary model file (see tempfile) using write.model

n.chains

number of Markov chains (default: 3)

n.iter

number of total iterations per chain (including burn in; default: 2000)

n.burnin

length of burn in, i.e. number of iterations to discard at the beginning. Default is n.iter/2, that is, discarding the first half of the simulations. If n.burnin is 0, jags() will run 100 iterations for adaption.

n.cluster

number of clusters to use to run parallel chains. Default equals n.chains.

n.thin

thinning rate. Must be a positive integer. Set n.thin > 1 to save memory and computation time if n.iter is large. Default is max(1, floor(n.chains * (n.iter-n.burnin) / 1000)) which will only thin if there are at least 2000 simulations.

DIC

logical; if TRUE (default), compute deviance, pD, and DIC. The rule pD=var(deviance) / 2 is used.

working.directory

sets working directory during execution of this function; This should be the directory where model file is.

jags.seed

random seed for JAGS, default is 123. This function is used for jags.parallell() and does not work for jags(). Use set.seed() instead if you want to produce identical result with jags()

.

jags.path

directory that contains the JAGS executable. The default is “”.

clearWD

indicating whether the files ‘data.txt’, ‘inits[1:n.chains].txt’, ‘codaIndex.txt’, ‘jagsscript.txt’, and ‘CODAchain[1:nchains].txt’ should be removed after jags has finished, default=TRUE.

refresh

refresh frequency for progress bar, default is n.iter/50

progress.bar

type of progress bar. Possible values are “text”, “gui”, and “none”. Type “text” is displayed on the R console. Type “gui” is a graphical progress bar in a new window. The progress bar is suppressed if progress.bar is “none”

digits

as in write.model in the R2WinBUGS package: number of significant digits used for BUGS input, see formatC. Only used if specifying a BUGS model as an R function.

RNGname

the name for random number generator used in JAGS. There are four RNGS supplied by the base moduale in JAGS: Wichmann-Hill, Marsaglia-Multicarry, Super-Duper, Mersenne-Twister

jags.module

the vector of jags modules to be loaded. Default are “glm” and “dic”. Input NULL if you don't want to load any jags module.

export_obj_names

character vector of objects to export to the clusters.

envir

default is .GlobalEnv

Details

To run:

  1. Write a BUGS model in an ASCII file.

  2. Go into R.

  3. Prepare the inputs for the jags function and run it (see Example section).

  4. The model will now run in JAGS. It might take awhile. You will see things happening in the R console.

BUGS version support:

Author(s)

Yu-Sung Su suyusung@tsinghua.edu.cn, Masanao Yajima yajima@stat.columbia.edu

References

Plummer, Martyn (2003) “JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.” http://citeseer.ist.psu.edu/plummer03jags.html.

Gelman, A., Carlin, J. B., Stern, H.S., Rubin, D.B. (2003) Bayesian Data Analysis, 2nd edition, CRC Press.

Sibylle Sturtz and Uwe Ligges and Andrew Gelman. (2005). “R2WinBUGS: A Package for Running WinBUGS from R.” Journal of Statistical Software 3 (12): 1–6.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
  # An example model file is given in:
  model.file <- system.file(package="R2jags", "model", "schools.txt")
  # Let's take a look:
  file.show(model.file)
  # you can also write BUGS model as a R function, see below:

#=================#
# initialization  #
#=================#

  # data
  J <- 8.0
  y <- c(28.4,7.9,-2.8,6.8,-0.6,0.6,18.0,12.2)
  sd <- c(14.9,10.2,16.3,11.0,9.4,11.4,10.4,17.6)


  jags.data <- list("y","sd","J")
  jags.params <- c("mu","sigma","theta")
  jags.inits <- function(){
    list("mu"=rnorm(1),"sigma"=runif(1),"theta"=rnorm(J))
  }

  ## You can input data in 4 ways
  ## 1) data as list of character
  jagsfit <- jags(data=list("y","sd","J"), inits=jags.inits, jags.params,
                n.iter=10, model.file=model.file)

  ## 2) data as character vector of names
  jagsfit <- jags(data=c("y","sd","J"), inits=jags.inits, jags.params,
                n.iter=10, model.file=model.file)

  ## 3) data as named list
  jagsfit <- jags(data=list(y=y,sd=sd,J=J), inits=jags.inits, jags.params,
                n.iter=10, model.file=model.file)

  ## 4) data as a file
  fn <- "tmpbugsdata.txt"
  dump(c("y","sd","J"), file=fn)
  jagsfit <- jags(data=fn, inits=jags.inits, jags.params,
                  n.iter=10, model.file=model.file)
  unlink("tmpbugsdata.txt")

  ## You can write bugs model in R as a function

  schoolsmodel <- function() {
    for (j in 1:J){                     # J=8, the number of schools
      y[j] ~ dnorm (theta[j], tau.y[j]) # data model:  the likelihood
      tau.y[j] <- pow(sd[j], -2)        # tau = 1/sigma^2
    }
    for (j in 1:J){
      theta[j] ~ dnorm (mu, tau)        # hierarchical model for theta
    }
    tau <- pow(sigma, -2)               # tau = 1/sigma^2
    mu ~ dnorm (0.0, 1.0E-6)            # noninformative prior on mu
    sigma ~ dunif (0, 1000)             # noninformative prior on sigma
  }

  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
                n.iter=10, model.file=schoolsmodel)


#===============================#
# RUN jags and postprocessing   #
#===============================#
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
    n.iter=5000, model.file=model.file)

  # Run jags parallely, no progress bar. R may be frozen for a while,
  # Be patient. Currenlty update afterward does not run parallelly
  #
   jagsfit.p <- jags.parallel(data=jags.data, inits=jags.inits, jags.params,
     n.iter=5000, model.file=model.file)

  # display the output
  print(jagsfit)
  plot(jagsfit)

  # traceplot
  traceplot(jagsfit.p)
  traceplot(jagsfit)

  # or to use some plots in coda
  # use as.mcmmc to convert rjags object into mcmc.list
  jagsfit.mcmc <- as.mcmc(jagsfit.p)
  jagsfit.mcmc <- as.mcmc(jagsfit)
  ## now we can use the plotting methods from coda
  #require(lattice)
  #xyplot(jagsfit.mcmc)
  #densityplot(jagsfit.mcmc)

  # if the model does not converge, update it!
  jagsfit.upd <- update(jagsfit, n.iter=100)
  print(jagsfit.upd)
  print(jagsfit.upd, intervals=c(0.025, 0.5, 0.975))
  plot(jagsfit.upd)

  # before update parallel jags object, do recompile it
  recompile(jagsfit.p)
  jagsfit.upd <- update(jagsfit.p, n.iter=100)



  # or auto update it until it converges! see ?autojags for details
  # recompile(jagsfit.p)
  jagsfit.upd <- autojags(jagsfit.p)
  jagsfit.upd <- autojags(jagsfit)

  # to get DIC or specify DIC=TRUE in jags() or do the following#
  dic.samples(jagsfit.upd$model, n.iter=1000, type="pD")

  # attach jags object into search path see "attach.bugs" for details
  attach.jags(jagsfit.upd)

  # this will show a 3-way array of the bugs.sim object, for example:
  mu

  # detach jags object into search path see "attach.bugs" for details
  detach.jags()

  # to pick up the last save session
  # for example, load("RWorkspace.Rdata")
  recompile(jagsfit)
  jagsfit.upd <- update(jagsfit, n.iter=100)

  recompile(jagsfit.p)
  jagsfit.upd <- update(jagsfit, n.iter=100)

#=============#
# using jags2 #
#=============#
  ## jags can be run and produces coda files, but cannot be updated once it's done
  ## You may need to edit "jags.path" to make this work,
  ## also you need a write access in the working directory:
  ## e.g. setwd("d:/")

  ## NOT RUN HERE
  ## Not run: 
    jagsfit <- jags2(data=jags.data, inits=jags.inits, jags.params,
      n.iter=5000, model.file=model.file)
    print(jagsfit)
    plot(jagsfit)
    # or to use some plots in coda
    # use as.mcmmc to convert rjags object into mcmc.list
    jagsfit.mcmc <- as.mcmc.list(jagsfit)
    traceplot(jagsfit.mcmc)
    #require(lattice)
    #xyplot(jagsfit.mcmc)
    #densityplot(jagsfit.mcmc)
  
## End(Not run)

R2jags documentation built on April 28, 2020, 1:08 a.m.