R, JAGS and ggplot2

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])

Density Plot - mcmcplots package

Trace Plot - mcmcplots package

Autoplot - mcmcplots package

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.

Density Plot - ggplot2 package
About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s