gelman-rubin

[d23yuipp0]

I want to obtain gelman-rubin statistics after running models with multiple chains. I believe this can be done with the gelman.diag() function in R.

The syntax is:

Code: Select all

gelman.diag(x, confidence = 0.95, transform=FALSE, autoburnin=TRUE, multivariate=TRUE)

where x is a mcmc.list object with more than one chain.

However, for the life of me, I cannot figure out how chains are stored and how to get them to be mcmc.list objects in order to analyze with gelman.diag().

Any insight would be greatly appreciated!

[ChrisCharlton]

You can find the MCMC chains in the chains slot of the mlwinfitMCMC object returned from the runMLwiN function. An example of this would be as follows:

Code: Select all

# for runMLwiN function
library(R2MLwiN)

# for gelman.diag function
library(coda)

# load the data
data(tutorial, package = "R2MLwiN")

# run the model
(mymodel <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, mcmcMeth = list(nchains = 3)), data = tutorial))

# calculate gelman-rubin diagnostics
gelman.diag(mymodel@chains)

[Norudomo]

Thank you for reminding about runMLwiN function, I had a similar issue and was thinking how to fix it yesterday at my Munich home. I will try to do as you wrote today and hope it's going to work, thank you!

[d23yuipp0]

Thanks, Chris. This was super helpful.