[R] Coda, HPDinterval and multiple chains

[NORRIS Paul]

Dear all,

I'm not 100% sure if this question is best directed at the r-list, or a mailing list concerned with Bayesian analysis, so please accept my apologies if another audience may be more appropriate.

I have been using the rjags package to run Jags models with multiple chains and store the results in a Coda based mcmc list.  For instance, having created a jags model and done initial adapting and updating, I run the following command:

coda_odp_gini_only <- coda.samples(odp_gini_only,
                                   variable.names=c("beta","sigma2.u2", "deviance"),
                                   n.iter=itercount, thin=thincount)

This create an object with 4 separate chains as requested in my initial call to Jags which created the object "odp_gini_only".  I then use the  Coda package to look at the results stored in "coda_odp_gini_only".  Prior to running any analysis with Coda, I used the command, 

coda.options(combine.plots=TRUE, combine.stats=TRUE)

to ask for results that combine the four separate chains.  Sure enough, if I enter "summary(coda_odp_gini_only)", I am given a single set of output combining the four chains.  However, if I enter "HPDinterval(coda_odp_gini_only)" I receive 4 sets of HPD figures, one for each chain.  Is it possible to combine the four chains together to receive a single set of HPD estimates? 

In a similar vein, is it possible to use the Coda object to estimate the proportion of a given parameters distribution which is above (or below) a given value, for instance, the proportion of the distribution of beta[1] greater than zero?  Again, in doing so, is it possible to combine the results of the four chains into a single estimate?

Kind regards, and many thanks in advance for any advice anyone can offer me.

Paul
Dr Paul Norris
Lecturer in Social Policy
University of Edinburgh
-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.