Support for jags and (open)bugs?

[Bill Harris]

I'm beginning to do more Bayesian analysis using jags and occasionally OpenBUGS.  There are various packages that help one set up, run, and then analyze those results inside of R; I've tried rjags and am currently using R2jags.  As you can see in the example in the R2jags manual (http://cran.r-project.org/web/packages/R2jags/index.html) for the jags command, they use plot, traceplot, xyplot, and densityplot to analyze results that are of class rjags or mcmc.list.  They seem to use base graphics and lattice.

Has anyone produced similar functions that use ggplot2 for working with these data types or other types produced by jags or perhaps openbugs?  I haven't found anything searching in this mailing list.

Thanks,

Bill

[Denis Haine]

Hello Bill,

Using as an example the school dataset from R2WinBUGS package (see its vignette):

library(rjags)
library(coda)
library(ggplot2)
library(R2WinBUGS)

data(schools)
schoolmodel <- function(){
  for (j in 1:J)
    {
      y[j] ~ dnorm (theta[j], tau.y[j])
      theta[j] ~ dnorm (mu.theta, tau.theta)
      tau.y[j] <- pow(sigma.y[j], -2)
    }
  mu.theta ~ dnorm (0.0, 1.0E-6)
  tau.theta <- pow(sigma.theta, -2)
  sigma.theta ~ dunif (0, 1000)
}
write.model(schoolmodel, "~/schoolmodel.bug")
J <- nrow(schools)
y <- schools$estimate
sigma.y <- schools$sd
forJags <- list(J = nrow(schools), y = schools$estimate, sigma.y = schools$sd)
inits <- function(){
  list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),
       sigma.theta = runif(1, 0, 100))
}

school.sim <- jags.model(file = "schoolmodel.bug", data = forJags,
                         inits = inits, n.chains = 3)
school.out <- coda.samples(school.sim, variable.names =
                           c("theta", "mu.theta", "sigma.theta"),
                           n.iter = 1000)
summary(school.out)

Than you can either use the mcmcplots packages which gives a "feel" of ggplot2:

library(mcmcplots)
denplot(school.out, parms = c('mu.theta', collapse = FALSE, auto.layout = TRUE))

Or go with ggplot2 after extracting the info you want from your output, like for example:

mu.theta.out <- cbind(school.out[[1]][,1], school.out[[2]][,1], school.out[[3]][,1])
attributes(mu.theta.out) <- NULL
dens.mu.theta <- density(mu.theta.out)
q25.mu.theta <- quantile(mu.theta.out, .025)
q975.mu.theta <- quantile(mu.theta.out, .975)
dd.mu.theta <- with(dens.mu.theta, data.frame(x,y))

To get a density plot for mu.theta with its 95% credible interval:

qplot(x, y, data = dd.mu.theta, geom = "line", ylab = "", xlab = "") +
  geom_ribbon(data = subset(dd.mu.theta, x>q25.mu.theta & x<q975.mu.theta),
              aes(ymax = y), ymin = 0, fill = "red", colour = NA, alpha = 0.5)

Denis

2011/5/14 Bill Harris <bill_...@facilitatedsystems.com>

--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: http://gist.github.com/270442
 
To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

[Bill Harris]

Thanks, Denis.  I didn't know about mcmcplots.  I'll check all this out.

Bill