# 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)
# Can also compute the DIC using pD (=Dbar-Dhat), via dic.samples(), which
# is a closer approximation to the original formulation of Spiegelhalter et
# al (2002), instead of pV (=var(deviance)/2), which is the default in JAGS
jagsfit.pD <- jags(data=jags.data, inits=jags.inits, jags.params,
n.iter=5000, model.file=model.file, pD=TRUE)
# 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
if (FALSE) {
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)
}
Run the code above in your browser using DataLab