# Deciphering output from R2WinBUGS

I have been struggling to figure out how to identify which output the parameter came from. For e.g. in the following code I run the model for 2 iterations. At the end of it I obtain 2 outputs named bugs.output[[1]], bugs.output[[2]]. Now when I type mean (p[,2]) it gives me the mean for p[2] from the second output. How do I specify that I need the mean for p[2] for the 1st or the ith output?

Here is how my output looks:

         mean       sd      2.5%       25%       50%       75%     97.5%
alpha      4.20338  1.14447   2.44975   3.39950   4.01550   4.82925   7.05405
p[1]       0.16691  0.04036   0.09974   0.13720   0.16445   0.19222   0.25308
p[2]       0.14286  0.03812   0.07665   0.11670   0.14020   0.16702   0.22630
p[3]       0.21238  0.04240   0.13740   0.18295   0.21050   0.23845   0.30360.......


And here is my code:

library("R2WinBUGS")
bugs.output <- list()
for(i in 1:2){
nausea <- as.integer(trial.data[i,])
bugs.output[[i]] <- bugs(
data=list(nausea=nausea, N=63),
inits=list(
list(theta=.300,mu=3),
list(theta=.350, mu=3)
),
model.file="conj_nausea_script_2.txt",
parameters.to.save = c("alpha","p"),
n.chains=2, n.iter=12000, n.burnin=5000,
bugs.directory="E://AChaudhuri/winbugs14/WinBUGS14",
working.directory=NULL)
attach.bugs(bugs.output[[i]],overwrite=NA)}


You complicated it a lot by calling bugs this way and using attach in a for loop! That's a big unclean thing. It's really hard to say what the result is. Attach should be used with care even in normal case, not even say in a for-loop. No wonder you got lost. Don't use attach this way!

• I completely agree; attach is dangerous, and best avoided. Dec 2, 2011 at 20:38
• without using attach if I type mean (p[,2]) it gives me the error Error in mean(p[, 2]) : object 'p' not found. How else can I read the output and do some analysis? Dec 5, 2011 at 22:23
• @anamika, you have to go directly to the structure: look at bugs.output[[i]]$mean$p. For more how to access the values see cran.r-project.org/web/packages/R2WinBUGS/vignettes/… Dec 8, 2011 at 9:39
• Thanks Tomas, This has been really helpful. I am new to this forum and I am pleased at how helpful everybody has been. I have seen this article earlier. I am actually trying to read in the quantiles. In normal case the 2.5% would be read as quantile(p[i],c(.025)) but in this case that doesnt seem to work. I dont see any reference on the article Dec 8, 2011 at 19:26
• @anamika, for quantiles use quantile function on bugs.output[[i]]\$sims.list (this is the posterior distribution) or something like that... or try to look into the bugs.output structure if there are the quantiles already computed. Dec 9, 2011 at 0:05

I suspect that the attach.bugs(bugs.output[[i]]...) when i=2 is masking the attach.bugs(bugs.output[[i]]...) when i=1. The overwrite=NA merely prevents deletion, not masking. Instead of "attach" you could refer directly to bugs.output[[1]] or bugs.output[[2]] in the call to mean. The results are stored in sims.array; reading the documentation of the bugs function return value in the R2WinBUGS manual should get you the rest of the way there.

I am not sure why you are making the two different bugs calls, each with two chains. Perhaps you could clarify your purpose? What I would have done is just run bugs once, specifying the number of chains, which in your case would seem to be four [(i in 1:2) x n.chains=2] and then print the result.

• My ultimate purpose is to be able to calculate coverage. I am making two different bugs calls because I am running the same model with 2 chains for 2 different datasets. I want to do this for say 300 datasets and then calculate coverage (out of 300 times how many times is the true value of p contained within the 95% prediction interval etc). This is the reason I was trying to identify mean or say any other descriptive specific to an output so that in the next loop I could call that. Dec 5, 2011 at 22:18