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

3

13

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.

Capitalistic answered 22/8, 2012 at 17:10 Comment(8)
What is wrong with the answer Abe provided to your question on Cross Validated? Could you post a figure showing the plot you want for the above example? You posted a figure on Cross Validated, but it does not appear to be for the example above.Matti
@MarkMiller the answer at Cross Validated is incomplete.Capitalistic
What specific addition results do you want? Abe's answer runs on my computer. Knowing what addition output you want would help people provide the necessary code. That is why I suggest you provide a figure for the example above showing exactly what you want.Matti
In your post on Cross Validated you provided a figure showing graphs of 80% interval for each chain, R-hat, and medians and 80% intervals. That is what Abe's answer provides with your example above. All I added to your code above was library(R2WinBUGS) and I added a missing parenthesis to Abe's plot statement (which I have now added to his post with a submitted edit).Matti
Your figure on Cross Validated includes plots of additional parameters perhaps because that figure is from a different example or perhaps because it is from a different model using the same data set and more parameters were monitored than in the example above. That is why I ask what additional results you want.Matti
Are you saying that you want to import an mcmc.list object, perhaps for example one that was emailed to you, without running the model yourself then run Abe's plot function?Matti
@MarkMiller Here is a link to an example mcmc.list, test.mcmc.list.Rdata. When I plot it (plot(as.bugs.array(sims.array = as.array(test.mcmc.list))), I get four plots on the right, labeled P1,P2,P3,P4. I expect one plot for each of the parameters; in the case of beta.site, with three lines (one for each level of the effect). Does that help?Capitalistic
I edited my answer to try to show how data extraction (and perhaps data subsetting) could be done to create custom plots, if that is an objective.Matti
M
4

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.

Matti answered 6/5, 2014 at 6:3 Comment(1)
Thanks. Its not so much that all I have is an mcmc.list so much as that this is embedded inside code that requires an mcmc.list for downstream processing (starting from line 101 here.Capitalistic
D
1

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) 
} 
Dolhenty answered 22/8, 2012 at 17:35 Comment(2)
This could be useful if I could figure out how to save an mcmc.list as a coda fileCapitalistic
See my answer below for one method of converting mcmc.list objects to coda. However, I believe @andybega's (Yu-Sung's) solution involves re-fitting and sampling.Guardroom
G
1

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())
Guardroom answered 8/5, 2014 at 6:14 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.