Runjags is reporting very low psrf = 1.0047 for a chain that apparently has convergence problems:
> print(o, vars = "q_date2")
JAGS model summary statistics from 3000 samples (chains = 3; adapt+burnin = 1000):
Lower95 Median Upper95 Mean SD MCerr MC%ofSD SSeff AC.10 psrf
q_date2 -0.17611 -0.10467 0.0053844 -0.087376 0.063296 0.023726 37.5 7 0.022495 1.0047
When I try to compute the psrf using coda, I get a result which looks much more reasonable:
> gelman.diag(as.mcmc.list(o)[,'q_date2'], transform=FALSE, autoburnin=FALSE)
Potential scale reduction factors:
Point est. Upper C.I.
[1,] 3.54 7.94
So why is the psrf reported by runjags so low? Is it a problem of runjags, or am I doing something wrong?
I use current version of runjags (1.2.1-0) in R 3.1.0.
EDIT: during creating summaries I got warnings - sorry for not mentioning them before:
Warning messages:
1: In autocorrs[x$stochastic] <- x$autocorr[4, ] :
number of items to replace is not a multiple of replacement length
2: In psrfs[x$stochastic] <- x$psrf$psrf[, 1] :
number of items to replace is not a multiple of replacement length
It appears (from information sent to me off-line) that the psrf is being calculated correctly, but is being reported out of sequence because some information is not available for some semi-stochastic monitored variables. The fact that this isn't being picked up by the software is a bug that I will fix!
In the meantime, you can either (a) ignore the psrf listed in the output for summary() and use RJout$psrf (or your own code) instead, or (b) remove the monitored variables (M in this case) that are causing the problem. The development version of runjags has even better solutions (summary statistics and plots are (re-)calculated after the model has been returned), and should be on CRAN in the next couple of months.
This is also a good reminder that manual checking of trace plots is an essential part of MCMC analysis :)