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.

### Like this:

Like Loading...

*Related*