Plot posterior densities of MCMC chains, by scenario, parameter and replicas, together with priors and true values

ggplot-posteriors_by_scenario.R

### ggplot MCMC posteriors and priors by scenario, parameter and replicas ###
### Author : Facundo Muñoz (facundo.munoz@orleans.inra.fr)
### Licence: GPL-3
### Date   : 2014-09-10


### Required packages
library(rlist)    # data-manipulation
library(plyr)
library(reshape2)
library(actuar)   # inverse gamma distribution
library(ggplot2)  # plotting

## Input: text files with MCMC chains
## each file holds one chain, corresponding to one single replica
## for one scenario and one parameter

## Assume files are in a directory 'mcmc-output', and are named
## like scni-parj-repk.txt, and they contain one single number per line
## without any header
datdir <- 'mcmc-output'
filepat <- '^scn[[:digit:]]+-par[[:digit:]]+-rep[[:digit:]]+.txt$'

# Simulate some fake chains
dir.create(datdir)
for( scn in 1:10 ) {
  for( par in 1:3 ) {
    for( rep in 1:6 ) {
      chain <- rgamma(1e3, shape = scn, scale = rep)
      fn <- paste('scn', scn, '-par', par, '-rep', rep, '.txt', sep = '')
      write.table(chain,
                  file = file.path(datdir, fn),
                  row.names = FALSE,
                  col.names = FALSE)
    }
  }
}

flist <- list.files(datdir, pattern = filepat)

# List of chains
# use the pattern as the names, removing the extension .txt
chainl <- sapply(file.path(datdir, flist), read.table)
names(chainl) <- gsub('\\.txt', '', flist)

# List of densities
denl <- list.select(lapply(chainl, density, cut = 1), x, y)

# Data.frame of densities
dendf <- ldply(denl, as.data.frame)

# Separate scenario from replica
dendf <- cbind(dendf, colsplit(dendf$.id, '-', c('scenario', 'parameter', 'replica')))

# remove scn, par and rep tags, and transform categorical variables to factors
rmtag <- function(x, tag) as.factor(as.numeric(gsub(tag, '', x)))
dendf <- transform(dendf,
                   scenario  = rmtag(scenario, 'scn'),
                   parameter = rmtag(parameter, 'par'),
                   replica   = rmtag(replica, 'rep'))



### True parametric values
### Each parameter has a true value at each scenario
### These are nonsense. Use true values here!
truedf <- data.frame(scenario = factor(1:10),
                     par1 = 2*(1:10),
                     par2 = 3*(1:10),
                     par3 = 4*(1:10))
truedf <- melt(truedf, variable.name = 'parameter', value.name = 'x')
truedf <- transform(truedf,
                    parameter = rmtag(parameter, 'par'))



### Prior distribution
### Scaled-inverted chi-squared distribution with parameters
### degrees of freedom: nu, scaling: tau2
### We use the fact that :
### Scale-Inv-Chisq(nu, tau2) = Inv-Gamma(nu/2, nu*tau2/2)
### We always use nu = 5, while tau2 equals the true value
### of the parameter for the corresponding scenario
build_prior <- function(true) {
  nu   <- 5
  tau2 <- true$x
  limits <- qinvgamma(c(.01, .99), shape = nu/2, scale = nu*tau2/2)
  x <- seq(limits[1], limits[2], length = 101)
  y <- dinvgamma(x, shape = nu/2, scale = nu*tau2/2)
  data.frame(scenario  = true$scenario,
             parameter = true$parameter,
             replica   = 0,
             x, y)
}
priordf <- list.stack(lapply(1:nrow(truedf), function(x) build_prior(truedf[x,])))


### Plot curves and true value
p <- ggplot(dendf, aes(x, y, group = replica)) + 
  geom_line() +
  geom_line(data = priordf, colour = 'darkred') +
  geom_vline(aes(xintercept = x), colour = 'darkgray', size = 1, data = truedf) +
  facet_grid(scenario ~ parameter, scales = 'free') + 
  theme_bw()

### Save plot (a4 size)
pdf(file = 'result.pdf', width = 8.27, height = 11.7)
print(p)
dev.off()


### Improvements
# Use proper labels for parameters and scenarios
# Automatically highlight those posteriors whose mass don't contain the true value