Databases
StackOverflow2013
Postsdbo
PO
 

Note that there are some explanatory texts on larger screens.

plurals
  1. PO
    primarykey
    data
    text
    <p>When you get error running your BUGS model from R, one option is to try a mock run of the model in OpenBUGS or WinBUGS itself. It can help you (via the cursor placement after you hit check model button) to locate problematic lines.</p> <p>I did this with your BUGS model. I found problems in the definition of <code>mn</code>, <code>prec</code> and <code>R</code> in the BUGS model. You can drop these as they are already defined in the data (which, looks like the appropriate place to define them). Once I dropped these from your BUGS model everything ran fine. </p> <p>Note, to run a model in OpenBUGS you have to edit the format of your data, for example the script I ran was:</p> <pre><code>model{ #likelihood for(j in 1 : Nf){ p1[j, 1:2 ] ~ dmnorm(gamma[1:2], T[1:2 ,1:2]) for (i in 1:2){ logit(p[j,i]) &lt;- p1[j,i] Y[j,i] ~ dbin(p[j,i],n) } } #priors gamma[1:2] ~ dmnorm(mn[1:2],prec[1:2 ,1:2]) expit[1] &lt;- exp(gamma[1])/(1+exp(gamma[1])) expit[2] &lt;- exp(gamma[2])/(1+exp(gamma[2])) T[1:2 ,1:2] ~ dwish(R[1:2 ,1:2], 2) sigma2[1:2, 1:2] &lt;- inverse(T[,]) rho &lt;- sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2]) } #data list(Y=structure(.Data=c(1,11,6,1,8,5,1,25,13,1,1,13,1,8,22),.Dim=c(5,3)), Nf=5, n=60, mn=c(-1.59,-2.44), prec=structure(.Data=c(0.0001,0,0,0.0001),.Dim=c(2,2)), R=structure(.Data=c(0.001,0,0,0.001),.Dim=c(2,2))) #inits list(gamma=c(0,0), T=structure(.Data=c(0.9,0,0,0.9),.Dim=c(2,2))) </code></pre> <p>where the data and inits need a bit work to convert from your R script. </p> <p>A couple of other points: 1) I am not sure you have the right format for Y as it has 3 columns, your distribution only considers the first two (X and Y1). 2) you had an unnecessary set of curly brackets in the likelihood.</p> <p>To run the code in BUGS via R you can use the following R syntax...</p> <pre><code>#BUGS code as a character string bugs1&lt;- "model{ #likelihood for(j in 1 : Nf){ p1[j, 1:2 ] ~ dmnorm(gamma[1:2], T[1:2 ,1:2]) for (i in 1:2){ logit(p[j,i]) &lt;- p1[j,i] Y[j,i] ~ dbin(p[j,i],n) } } #priors gamma[1:2] ~ dmnorm(mn[1:2],prec[1:2 ,1:2]) expit[1] &lt;- exp(gamma[1])/(1+exp(gamma[1])) expit[2] &lt;- exp(gamma[2])/(1+exp(gamma[2])) T[1:2 ,1:2] ~ dwish(R[1:2 ,1:2], 2) sigma2[1:2, 1:2] &lt;- inverse(T[,]) rho &lt;- sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2]) }" #write the BUGS code to a txt file in current working directory writeLines(bugs1, "bugs1.txt") #create data Y&lt;-data.frame(X=1,Y1=c(11,8,25,1,8),Y2=c(6,5,13,13,22)) #run BUGS from R library("R2OpenBUGS") mcmc1 &lt;- bugs(data = list(Y=as.matrix(Y), Nf=5, n=60, mn=c(-1.59, -2.44), prec=matrix(c(.0001,0,0,.0001),nrow=2,ncol=2), R=matrix(c(.001,0,0,.001),nrow=2,ncol=2)), inits = list(list(gamma=c(0,0), T=matrix(c(0.9,0,0,0.9),nrow=2,ncol=2))), param = c("gamma", "sigma2"), model = "bugs1.txt", n.iter = 11000, n.burnin = 1000, n.chains = 1) </code></pre> <p>A couple of points to note here. 1) This uses OpenBUGS not WinBUGS. 2) If you use R2WinBUGS you might hit a trap if you are not running R (or Rstudio, or whatever you are using) as an administrator.</p> <p>To run the above code a 1000 times you could put it within a loop, something like....</p> <pre><code>#create and write the BUGS code to a txt file in current working directory (outside the loop) bugs1&lt;-... #loop for(i in 1:1000){ Y &lt;- read.csv(file=paste0("MVN",i,".csv")) #run BUGS from R library("R2OpenBUGS") mcmc1 &lt;- bugs(data = list(Y=as.matrix(Y), Nf=5, n=60, mn=c(-1.59, -2.44), prec=matrix(c(.0001,0,0,.0001),nrow=2,ncol=2), R=matrix(c(.001,0,0,.001),nrow=2,ncol=2)), inits = list(list(gamma=c(0,0), T=matrix(c(0.9,0,0,0.9),nrow=2,ncol=2))), param = c("gamma", "sigma2"), model = "bugs1.txt", n.iter = 11000, n.burnin = 1000, n.chains = 1) #save mcmc write.csv(mcmc1$sims.matrix,paste0("mcmc",i,".csv")) } </code></pre>
    singulars
    1. This table or related slice is empty.
    1. POAutomation of WinBUGS from R
      singulars
      1. PTQuestion
    1. PTAnswer
    1. USgjabel
    1. USgjabel
    plurals
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. POAutomation of WinBUGS from R
      singulars
      1. PTQuestion
    1. This table or related slice is empty.
    1. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    3. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTAcceptedByOriginator
    1. CO@gjable,Thankyou for your response. The model ran fine with BUGS even earlier as I had mentioned. There was an error in this window when I pasted the code. I didnot mention mn, prec and R in code as they were listed under data. The problem is in how R is reading in the datafile. I have edited my code above for R and if I use your model code I still get the trap error. I am looking for some guidance in the R code as I am not so familiar with R calling WinBUGS.
      singulars
      1. PO
      1. USuser1560215
    2. CO@user1560215 updated answer with R code. did this help?
      singulars
      1. PO
      1. USgjabel
    3. COIt runs fine with 1 dataset, I just tested it. But how can I automate this for 1000 datastes. The first part is that how do I read in the 1000 simulated datasets? Thanks for all the help.
      singulars
      1. PO
      1. USuser1560215
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload