Bayes-Example-HierarchicalModel.r

rm(list=ls(all=TRUE))

library(rjags)

# assuming the data is created from an ecological system that creates an
# exponential size distribution (e.g. you sample individuals from a population that can be 
# expected to follow this distribution), but this measurments are done with
# a considerable lognormal observation error
# for a realistic application see http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0058036

meanSize <- 10
trueLogSd <- 1
sampleSize <- 500
truevalues = rexp(rate = 1/meanSize, n = sampleSize)
observations  =  rlnorm(n = length(truevalues), mean = log(truevalues), sd = trueLogSd)


# plotting true and observed data
maxV <- ceiling(max(observations,truevalues))
counts <- rbind(
  obs = hist(observations, breaks = 0:maxV, plot = F)$counts,
  true = hist(truevalues, breaks = 0:maxV, plot = F)$counts
)
barplot(t(counts), beside=T)

################################################################
# Model specification of a non-hierarchical model in JAGS
# that does not account for the observation error
# textConnection avoids having to write the string to file (default option)
normalModel = textConnection('
    model {
    
    # Priors
    meanSize ~  dunif(1,100) 

    # Likelihood
    for(i in 1:nObs){
      true[i] ~ dexp(1/meanSize)
    }
    }
')

# Bundle data
positiveObservations <- observations[observations>0]
data = list(true = positiveObservations, nObs=length(positiveObservations))

# Parameters to be monitored (= to estimate)
params = c("meanSize")

jagsModel = jags.model( file= normalModel , data=data, n.chains = 3, n.adapt= 500)
results = coda.samples( jagsModel , variable.names=params,n.iter=5000)
plot(results)

# I explain how to interpret these plots in http://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/ 

# note that parameter estmiates are biased

#################################################################
# Model specification if hierarchical model that accounts for the 
# observation error in Jags

hierarchicalModel = textConnection('
    model {
    
    # Priors
    meanSize ~  dunif(1,100) 
    sigma ~ dunif(0,20) # Precision 1/variance JAGS and BUGS use prec instead of sd
    precision <- pow(sigma, -2)

    # Likelihood
    for(i in 1:nObs){
      true[i] ~ dexp(1/meanSize)
      observed[i] ~ dlnorm(log(true[i]), precision)
    }
    }
')
 
# Bundle data
data = list(observed = observations, nObs=length(observations))
 
# Parameters to be monitored (= to estimate)
params = c("meanSize", "sigma")

jagsModel = jags.model( file= hierarchicalModel , data=data, n.chains = 3, n.adapt= 500)
#update(jagsModel, 2500) # updating without sampling
results = coda.samples( jagsModel , variable.names=params,n.iter=5000)


plot(results)

# it's allways good to check the correlation structure in the posterior

# check correlations using
 pairs(as.matrix(results))
# or using https://gist.github.com/florianhartig/8423115 you can plot
# betterPairs(as.matrix(results))

# see http://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/ for comments on interpreting correlations