Last week a question was asked on the ggplot2 list about using ggplot2 and jags in R (). Here’s what was my answer (a bit updated):

Using as an example the school dataset from R2WinBUGS package:

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, n.adapt = 500) school.out <- coda.samples(school.sim, variable.names = c("theta", "mu.theta", "sigma.theta"), n.iter = 1000) summary(school.out)

Than you can use the mcmcplots package which give a “feel” of ggplot2:

library(mcmcplots) denplot(school.out, parms = c('mu.theta', collapse = FALSE, auto.layout = TRUE)) traplot(school.out, parms = c('mu.theta')) autplot1(school.out[,"mu.theta", drop = FALSE])

If you really want to use ggplot2, you have to extract the information you need, like for example:

varnames(school.out) str(school.out) 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)) 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)

This would give you the density plot with its 95% credible interval.

Advertisements