Convergence!

Check out these pretty tools:
mcmcplot()
ggmcmc()

 
# Example time!
library(R2jags)
?gelman.diag
?geweke.diag
library(car)
data(Prestige)
df <- Prestige[,1:4]

# model 1 weak priors
modelString = "
model {
  for(i in 1:N){
    prestige[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1 * income[i] + beta2 * education[i]
  }
  beta0 ~ dnorm(0, .01) 
  beta1 ~ dunif(-100, 100)
  beta2 ~ dunif(-100, 100)
  tau ~ dgamma(.01, .01) # note: JAGS uses shape and rate for gamma parameters
}
"
# plotting priors to play around with parameters
plot(density(rgamma(10000, shape=3, rate=2)), xlim=c(0,10))

# model 2 strong priors
modelString = "
model {
  for(i in 1:N){
    prestige[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1 * income[i] + beta2 * education[i]
  }
  beta0 ~ dnorm(0, .01) 
  beta1 ~ dgamma(1, .1)
  beta2 ~ dgamma(1, .1)
  tau ~ dgamma(.01, .01)
}
"
writeLines( modelString , con="model.jags" )

prestige <- df$prestige
income <- df$income
education <- df$education
N <- length(df$prestige)

data <- list(prestige=prestige, income=income, education=education, N=N)
params <- c("beta0", "beta1", "beta2")
inits1 <- list("beta0"=0, "beta1"=0, "beta2"=0)
inits2 <- list("beta0"=-7, "beta1"=1, "beta2"=4) # sorta from the OLS estimates
inits3 <- list("beta0"=10, "beta1"=10, "beta2"=10)
inits <- list(inits1, inits2, inits3)

set.seed(123)

fit <- jags(data = data, inits = inits,
                   parameters.to.save = params, n.chains = 3, n.iter = 9000,
                   n.burnin = 1000, model.file = "model.jags")

print(fit)
summary(lm(prestige ~ income + education, data=Prestige)) # to compare

fit.mcmc <- as.mcmc(fit)

summary(fit.mcmc)
plot(fit)
traceplot(fit)

xyplot(fit.mcmc)
xyplot(fit.mcmc, aspect="fill")
densityplot(fit.mcmc)
densityplot(fit.mcmc, aspect="fill")
autocorr.plot(fit.mcmc)

gelman.plot(fit.mcmc)
geweke.diag(fit.mcmc)
raftery.diag(fit.mcmc)
heidel.diag(fit.mcmc)

# one-stop shopping:
library(superdiag)
superdiag(fit.mcmc, burnin = 100)

library(mcmcplots)
denplot(fit.mcmc)
denplot(fit.mcmc, parms = c("beta0", "beta1", "beta2"))
traplot(fit.mcmc, parms = c("beta0", "beta1", "beta2"))
traplot(fit.mcmc[1], parms = c("beta0", "beta1", "beta2"))
traplot(fit.mcmc[2], parms = c("beta0", "beta1", "beta2"))
traplot(fit.mcmc[3], parms = c("beta0", "beta1", "beta2"))
mcmcplot(fit.mcmc) # whoa!!!

caterplot(fit.mcmc)
caterplot(fit.mcmc, parms = c("beta0", "beta1", "beta2"),
          labels = c("Intercept", "income", "education"))

library(ggmcmc) # pretty :)
fit.gg <- ggs(fit.mcmc)
ggs_density(fit.gg)
ggs_traceplot(fit.gg)
ggmcmc(fit.gg) # makes gorgeous compilation of plots saved as ggmcmc-output.pdf
# Produces density plots of the whole chain (black) with the density plot of only the latest part of the chain (green). By default, it uses the latest 10% of the chain as a comparison.