Monday, February 11, 2008

A wraper function to convert coda files into a BUGS object

I used to fit Bayesian model using WinBUGS via R a lot. But now I am more and more prone to run models on OpenBUGS directly. I have document the reason why I like OpenBUGS and wrote a auto OpenBUGS function here. In short, I like to be able to know what's going on why my MCMC is running. However, exporting MCMC results from OpenBUGS to R is the necessary step to do the post analysis. 'coda' package has functions that can read coda files generated from WinBUGS or OpenBUGS. One annoying thing we need to take care of is the filenames for the coda files:

For WinBUGS, file names are like `coda1.txt', `coda2.txt', `coda3.txt' and `codaIndex.txt'. For OpenBUGS, file names are like `CODAchain1.txt', `CODAchain2.txt', `CODAchain3.txt' and `CODAindex.txt'. So to read coda files into R:

coda.WinBUGS <- read.coda(stem)
coda.OpenBUGS <- read.openbugs(stem)

where stem is the path you store the coda files. This will convert coda files into `mcmc.list' object. You can do the following tests for convergence:

gelman.diag(coda.WinBUGS)
geweke.diag(coda.WinBUGS)
heidel.diag.diag(coda.WinBUGS)


However, if you want to do more analysis and extract the values from the coda files, converting coda files into `bugs' object is one of the ways to go. Here is the wraper function to do this:

coda2bugs <- function(path, para, n.chains=3, n.iter=5000,
 n.burnin=1000, n.thin=2){   
 setwd(path)   
 library(R2WinBUGS)   
 fit <- R2WinBUGS:::bugs.sims(para, n.chains=n.chains, 
 n.iter=n.iter, n.burnin=n.burnin, n.thin=n.thin, 
 DIC = FALSE)   class(fit) <- "bugs"   
 return(fit) 
} 
To see the summary statistics:

fit <- coda2bugs("c:/")
print(fit)
plot(fit)

To extract values of parameters:

attach.bugs(fit)


Note: To use this function, you need to save coda files as the WinBUGS file formate as stated above.

1 comment:

David said...

Do you know how to convert an mcmc.list to a bugs object?