Calculating Gelman and Rubins convergence statistic for only a subset of iterations (coda package)

Question

I am trying to calculate Gelman and Rubin's convergence diagnostic for a JAGS analysis I am currently running in R using the R package rjags.

For example, I would like to assess the convergence diagnostic for my parameter beta. To do this, I am using the library coda and the command:

library(coda)
gelman.diag(out_2MCMC$beta)

with out_2MCMC being an MCMC list object with more than one chain, resulting in correct output with no error messages or whatsoever. However, as I used a high amount of iterations as burn-in, I would like to calculate the convergence diagnostic for only a subset of iterations (the part after the burn in only!).

To do this, I tried:

gelman.diag(out_2MCMC$beta[,10000:15000,])

This resulted in the following error:

Error in mcmc.list(x) : Arguments must be mcmc objects

Therefore, I tried:

as.mcmc(out_2MCMC$beta[,10000:15000,]))

But, surprisingly this resulted in the following error:

Error in gelman.diag(as.mcmc(out_2MCMC$beta[, 10000:15000,])) 
You need at least two chains

As this is the same MCMC list object as I get from the JAGS analysis and the same as I am using when I am assessing the convergence diagnostic for all iterations (which works just perfect), I don't see the problem here.

The function itself only provides the option to use the second half of the series (iterations) in the computation of the convergence diagnostic. As my burn in phase is longer than that, this is unfortunately not enough for me.

I guess that it is something very obvious that I am just missing. Does anyone have any suggestions or tipps?

As it is a lot of code, I did not provide R code to run a full 2MCMC-JAGS analysis. I hope that the code above illustrates the problem well enough, maybe someone encountered the same problem before or recognizes any mistakes in my syntax. If you however feel that the complete code is necessary to understand my problem I can still provide example code that runs a 2MCM JAGS analysis.

Answer 1

I was looking for a solution to the same problem and found that the window() function in the stats package will accomplish the task. For your case:

window(out_2MCMC, 100001,15000)

will return an mcmc.list object containing the last 5,000 samples for all monitored parameters and updated attributes for the new list.