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.
80% interval for each chain
,R-hat
, andmedians and 80% intervals
. That is what Abe's answer provides with your example above. All I added to your code above waslibrary(R2WinBUGS)
and I added a missing parenthesis to Abe'splot
statement (which I have now added to his post with a submitted edit). – Mattimcmc.list
object, perhaps for example one that was emailed to you, without running the model yourself then run Abe'splot
function? – Mattiplot(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