This week I want to take a more detailed look at the use of JAGS with Stata and in particular I want to contrast JAGS with WinBUGS by analysing the biopsy data that I described last time. This posting really needs to be read in sequence with the previous two.<\/p>\n
Before I\u00a0start on the comparison,\u00a0let me say a little more about the JAGS syntax. I explained when I introduced JAGS that a WinBUGS program is likely to run in JAGS without alteration and this is more or less true, but there are a few potentially useful features of the JAGS language that derive from the original, sometimes called classic, version of BUGS that preceded WinBUGS.<\/p>\n
Some JAGS syntax<\/strong><\/p>\n In JAGS we have the option of declaring the dimensions of vectors and matrices in the model file. This is done by having a var<\/strong> block before the model starts. If the declaration is omitted, JAGS tries to work out the dimensions of the structures and it is usually able to do this successfully.<\/p>\n Second we can end an assignment or a distributional command with a semi-colon. In this way we can place multiple commands on a single line.<\/p>\n Finally we are allowed to implicitly define a calculated quantity (logical node) within a distributional statement so we could replace the WinBUGS commands<\/p>\n with a single JAGS command,<\/p>\n The only reason for using the former syntax in JAGS would be if we wanted to save the values of mu.<\/p>\n A JAGS program for analysing the biopsies<\/strong><\/p>\n When JAGS is installed a folder is created called classic-bugs<\/strong> that contains the files needed for fitting all of the WinBUGS examples in JAGS. The model file for the biopsy example is called biops.bug<\/strong> and contains<\/p>\n If this file is contrasted with the WinBUGS model file for the same problem you will see that here prior<\/strong> refers to a 4×4 matrix, which is defined in the data file to be,<\/p>\n 1\u00a0\u00a0 0\u00a0\u00a0 0\u00a0\u00a0 0 while previously prior was a vector of length 4 that was filled with 1\u2019s i.e. (1,1,1,1)<\/p>\n Another difference is that in JAGS, nbiops<\/strong>, the total number of biopsies for each heart, is read in as data while in the WinBUGS code it was calculated.<\/p>\n Other differences in the model file are that this JAGS model file defines the sizes of all of the vectors and matrices in a var<\/strong> block and the lines are ended with semi-colons.<\/p>\n The initial values are different in that JAGS does not give initial values for error[]<\/strong> but does give them for p<\/strong>[] and for true<\/strong>[]. In fact the way in which the structured zeros are handled is\u00a0quite different to the method used by WinBUGS. In JAGS it is not necessary to define error<\/strong>[] in both the data and the initial values files with NA in the appropriate places. In JAGS we can define error<\/strong>[]\u00a0implicitly and still estimate some of its terms as parameters. This seems to me to be convenient but a little illogical.<\/p>\n Running JAGS from Stata<\/strong><\/p>\n We can now prepare a Stata program that calls JAGS and fits the biopsy model.<\/p>\n\u00a0\u00a0\u00a0\u00a0 mu[i] <- a + b*x[i]<\/span><\/strong><\/pre>\n\u00a0\u00a0\u00a0\u00a0 y[i] ~ dnorm(mu[i],tau)<\/span><\/strong><\/pre>\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 y[i] ~ dnorm(a+b*x[i],tau)<\/span><\/strong><\/pre>\nvar\r\n\u00a0\u00a0 biopsies[ns,4], # grades observed in ith session (multinomial)\r\n\u00a0\u00a0 nbiops[ns],\u00a0\u00a0\u00a0\u00a0 # total number of biopsies in ith session\r\n\u00a0\u00a0 true[ns],\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 # true state in ith \r\n\u00a0\u00a0 error[4,4],\u00a0\u00a0\u00a0\u00a0 # error matrix in taking biopsies\r\n\u00a0\u00a0 prior[4,4],\u00a0\u00a0\u00a0\u00a0 # prior parameters for rows of error[,]\r\n\u00a0\u00a0 p[4];\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 # underlying\u00a0\u00a0 incidence of true states\r\nmodel {\r\n\u00a0\u00a0 for (i in 1:ns)\r\n\u00a0\u00a0\u00a0\u00a0 true[i]\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 ~ dcat(p[]);\r\n\u00a0\u00a0\u00a0\u00a0 biopsies[i,] ~ dmulti(error[true[i],],nbiops[i]); \r\n\u00a0\u00a0 }\r\n\u00a0\u00a0 for (j in 1:4)\r\n error[j,]\u00a0\u00a0\u00a0\u00a0 ~ ddirch(prior[j,]);\r\n }\r\n\u00a0\u00a0 p[]\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 ~ ddirch(prior[4,]);\u00a0\u00a0\u00a0\u00a0 # prior for p\r\n}<\/strong><\/span><\/pre>\n
\n1\u00a0\u00a0 1\u00a0\u00a0 0\u00a0\u00a0 0
\n1\u00a0\u00a0 1\u00a0\u00a0 1\u00a0\u00a0 0
\n1\u00a0\u00a0 1\u00a0\u00a0 1\u00a0\u00a0 1<\/p>\n*---------------------------------\r\n* write data file\r\n*---------------------------------\r\nuse biopsy.dta, clear\r\nmatrix prior = J(4,4,1)\r\nforvalues r=1\/4 {\r\n\u00a0\u00a0 forvalues c=2\/4 {\r\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if `c' > `r' matrix prior[`r',`c'] = 0\r\n}\r\n }\r\n egen nbiops = rsum(grade*)\r\n gen true = 4\r\n scalar ns = 157\r\n wbslist_v2 (scalar ns , format(%5.0f)) (matrix prior , format(%4.0f)) \/\/\/\r\n (structure grade1 grade2 grade3 grade4 , name(biopsies) format(%3.0f)) \/\/\/\r\n (vector nbiops , format(%3.0f) ) \/\/\/\r\n using dataj.txt , replace jags\r\n *---------------------------------\r\n * write initial values file\r\n *---------------------------------\r\n matrix p = J(1,4,0.25)\r\n wbslist_v2 (matrix p) (vec true, format(%3.0f) ) \/\/\/\r\n using initj.txt , replace jags\r\n *---------------------------------\r\n * write model file\r\n *---------------------------------\r\n wbsmodel biopsy_jags.do modelj.txt\r\n \/*\r\n var\r\n biopsies[ns,4], # grades observed in ith session (multinomial)\r\n nbiops[ns], # total number of biopsies in ith session\r\n true[ns], # true state in ith session\r\n error[4,4], # error matrix in taking biopsies\r\n prior[4,4], # prior parameters for rows of error[,]\r\n p[4]; # underlying incidence of true states\r\n model {\r\n for (i in 1:ns){\r\n true[i] ~ dcat(p[]);\r\n biopsies[i,] ~ dmulti(error[true[i],],nbiops[i]);\r\n }\r\nfor (j in 1:4) {\r\n error[j,] ~ ddirch(prior[j,]);\r\n }\r\np[] ~ ddirch(prior[4,]); # prior for p\r\n }\r\n *\/\r\n *---------------------------------\r\n * write script file\r\n *---------------------------------\r\n wbsscript_v2 using scriptj.txt , replace \/\/\/\r\n model(modelj.txt) data(dataj.txt) init(initj.txt) \/\/\/\r\n burnin(1000) update(10000) \/\/\/\r\n set(p error) coda(bsj) jags\r\n *---------------------------------\r\n * run and read the results\r\n *---------------------------------\r\n wbsrun_v2 using scriptj.txt , jags\r\n type jagslog.txt\r\n wbscoda_v2 using bsj , clear jags\r\n *---------------------------------\r\n * inspect the results\r\n *---------------------------------\r\n mcmcstats *\r\n mcmctrace p*<\/span><\/strong><\/pre>\n