3

I am new to both Bayes and JAGS, so please forgive my ignorance.

I have been sent an R script (by a colleague) which is written using JAGS code.

The author of this code has defined the set of coda samples as below:

codaSamples = coda.samples( jagsModel , variable.names=parameters , 
                            n.iter=nPerChain , thin=thinSteps )

I wish to obtain the following, and have had limited success:

Gelman diagnostics: I have used "show(gelman.diag(codaSamples))" which is fine for a single simulation. However, how do I output to a file each of the gelman diagnostics, per parameter, for every simulation of interest? Of more interest, is it possible just to record the proportion of simulations where the Rhat value is >1.1?

Density plot: I have used "show(densplot(codaSamples))". However, this produces each plot on a separate plot (I have 96 parameters in the model). Is there an equivalence to "autocorr.plot", which places several plots per page?

Quantiles: I used "show (summary(codaSamples))", but although this gave the mean, SD, and specific centiles for each parameter (which is what I wanted), it also gave the MCMC matrix. Is there anyway in which to just specify the basic statistical properties for each parameter?

Posterior distribution: Is there a way to calculate the centile at which a given value (say zero), of each parameter lies below/above? Then to summarise across all simulations?

Thank you in advance for any help that you can offer.

agstudy
  • 119,832
  • 17
  • 199
  • 261
user3628889
  • 271
  • 1
  • 5
  • 10

1 Answers1

2
  • record gelman output with something like:

write.csv(gelman.diag(codaSamples)$psrf, file = "gelman.csv")

  • to plot density of mcmc: mcmc object has a specific S3 class mcmc.list so class(codaSamples) should return mcmc.list. There is a plot S3 method for this class of object with arguments: trace and density.

plot(codaSamples, trace = FALSE, density = TRUE)

That should correspond to your expectations.

  • summary(codaSamples) gives statistic of posteriors for the parameters. Not sure I understand the question, but codaSamples is a list of length: the number of chains, in columns: your estimeted parameters, and in rows: each iterations of the chain. So you can calculate every statistics you want with this object or sub-object (quantiles, ...). For example the quantile for the first parameter :

quantile(codaSamples[[1]][,1])

[[1]] because I take the first chain and [,1] for the first parameters.

To look at the structure of an object you need to use str(). With this function you can select all sub-object you want.

Philippe
  • 194
  • 1
  • 12