How can I convert an mcmc.list to a bugs object?

问题:

I am using the rjags R library. The function coda.samples produces an mcmc.list, for example (from example(coda.samples)):

library(rjags)
data(LINE)
LINE$recompile()
LINE.out <- coda.samples(LINE, c("alpha","beta","sigma"), n.iter=1000)
class(LINE.out)
[1] "mcmc.list"

However, I would like to use the plot.bugs function, which requires a bugs object as input.

Is it possible to convert an object from an mcmc.list to a bugs object, so that plot.bugs(LINE.out)?

Note that there is a similar question on stats.SE that has been unanswered for over a month. That question had a bounty that ended on 08/29/2012.

More hints:

I have discovered that the R2WinBUGS package has a function "as.bugs.array" function - but it is not clear how the function can be applied to an mcmc.list.

回答1:

I do not know whether this will give you what you want. Note that the model code came from using your code and then typing LINE at the cursor. The rest is just standard bugs code, except I used tau = rgamma(1,1) for an initial value and do not know how standard that is. More than once I have seen tau = 1 used as an initial value. Perhaps that would be better.

In effect, I created an rjags object using the same model code you were using and added a jags statement to run it. I admit that is not the same thing as converting coda output to a bugs object, but it might result in you getting the desired plot.

If all you have is an mcmc.list and no model code and you simply want to plot the mcmc.list, then my answer will not help.

library(R2jags)

x <- c(1, 2, 2, 4, 4,  5,  5,  6,  6,  8) 
Y <- c(7, 8, 7, 8, 9, 11, 10, 13, 14, 13) 

N <- length(x)
xbar <- mean(x)

summary(lm(Y ~ x))

x2 <- x - xbar

summary(lm(Y ~ x2))

# Specify model in BUGS language

sink("model1.txt")

cat("

model  {
                for( i in 1 : N ) {
                        Y[i] ~ dnorm(mu[i],tau)
                        mu[i] <- alpha + beta * (x[i] - xbar)
                }
                tau ~ dgamma(0.001,0.001) 
                sigma <- 1 / sqrt(tau)
                alpha ~ dnorm(0.0,1.0E-6)
                beta ~ dnorm(0.0,1.0E-6)        
        }

",fill=TRUE)
sink()

win.data <- list(Y=Y, x=x, N=N, xbar=xbar)

# Initial values
inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), tau = rgamma(1,1))}

# Parameters monitored
params <- c("alpha", "beta", "sigma")

# MCMC settings
ni <- 25000
nt <-     5
nb <-  5000
nc <-     3

out1 <- jags(win.data, inits, params, "model1.txt", n.chains = nc, 
             n.thin = nt, n.iter = ni, n.burnin = nb)

print(out1, dig = 2)
plot(out1)

#library(R2WinBUGS)
#plot(out1)

EDIT:

Based on the comments perhaps something like this will help. The line str(new.data) suggests that a large amount of data are available. If you are simply trying to create variations of default plots then doing so may only be a matter of extracting and subsetting the data as desired. Here plot(new.data$sims.list$P1) is just one straight-forward example. Without knowing exactly what plot you want I will not attempt more specific data extractions. If you post a figure showing an example of the exact kind of plot you want perhaps someone can take it from here and post the code needed to create it.

By the way, I recommend reducing the size of the example data set to perhaps three chains and perhaps no more than 30 iterations until you have the exact code you want for the exact plot you want:

load("C:/Users/mmiller21/simple R programs/test.mcmc.list.Rdata")

class(test.mcmc.list)

library(R2WinBUGS)

plot(as.bugs.array(sims.array = as.array(test.mcmc.list)))

new.data <- as.bugs.array(sims.array = as.array(test.mcmc.list))

str(new.data)

plot(new.data$sims.list$P1)

EDIT:

Note also that:

class(new.data)
[1] "bugs"

whereas:

class(test.mcmc.list)
[1] "mcmc.list"

which is what the title of your post requests.

回答2:

This is not a solution to your question, but in response to your comment on @andybega's answer, here's a way to convert an mcmc.list object to typical coda text files.

mcmc.list.to.coda <- function(x, outdir=tempdir()) {
  # x is an mcmc.list object
  x <- as.array(x)
  lapply(seq_len(dim(x)[3]), function(i) {
    write.table(cbind(rep(seq_len(nrow(x[,,i])), ncol(x)), c(x[,,i])), 
                paste0(outdir, '/coda', i, '.txt'),
                row.names=FALSE, col.names=FALSE)
  })

  cat(paste(colnames(x), 
            1 + (seq_len(ncol(x)) - 1) * nrow(x),
            nrow(x) * seq_len(ncol(x))), 
      sep='\n', 
      file=file.path(outdir, 'codaIndex.txt'))
}

# For example, using the LINE.out from your question:
mcmc.list.to.coda(LINE.out, tempdir())

回答3:

Not an answer, but this blog post has the following wrapper function for converting coda output (.txt) to BUGS using R2WinBUGS:::bugs.sims:

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) 
}