0
$\begingroup$

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")
trial.data <- read.table("simuldatBB_6_30.csv", header=T)
p_true<- read.table("p_BB_6_30.csv",header=F)
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)}
Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

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!

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ I completely agree; attach is dangerous, and best avoided. $\endgroup$
    – jbowman
    Commented Dec 2, 2011 at 20:38
  • $\begingroup$ 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? $\endgroup$
    – Anamika
    Commented Dec 5, 2011 at 22:23
  • $\begingroup$ @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/… $\endgroup$
    – Tomas
    Commented Dec 8, 2011 at 9:39
  • $\begingroup$ 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 $\endgroup$
    – Anamika
    Commented Dec 8, 2011 at 19:26
  • $\begingroup$ @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. $\endgroup$
    – Tomas
    Commented Dec 9, 2011 at 0:05
1
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Anamika
    Commented Dec 5, 2011 at 22:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.