JAGS with Stata II

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.

Before I start on the comparison, let 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.

Some JAGS syntax

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 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.

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.

Finally we are allowed to implicitly define a calculated quantity (logical node) within a distributional statement so we could replace the WinBUGS commands

     mu[i] <- a + b*x[i]

     y[i] ~ dnorm(mu[i],tau)

with a single JAGS command,

           y[i] ~ dnorm(a+b*x[i],tau)

The only reason for using the former syntax in JAGS would be if we wanted to save the values of mu.

A JAGS program for analysing the biopsies

When JAGS is installed a folder is created called classic-bugs 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 and contains

var
   biopsies[ns,4], # grades observed in ith session (multinomial)
   nbiops[ns],     # total number of biopsies in ith session
   true[ns],       # true state in ith 
   error[4,4],     # error matrix in taking biopsies
   prior[4,4],     # prior parameters for rows of error[,]
   p[4];           # underlying   incidence of true states
model {
   for (i in 1:ns)
     true[i]       ~ dcat(p[]);
     biopsies[i,] ~ dmulti(error[true[i],],nbiops[i]); 
   }
   for (j in 1:4)
       error[j,]     ~ ddirch(prior[j,]);
   }
   p[]       ~ ddirch(prior[4,]);     # prior for p
}

If this file is contrasted with the WinBUGS model file for the same problem you will see that here prior refers to a 4×4 matrix, which is defined in the data file to be,

1   0   0   0
1   1   0   0
1   1   1   0
1   1   1   1

while previously prior was a vector of length 4 that was filled with 1’s i.e. (1,1,1,1)

Another difference is that in JAGS, nbiops, the total number of biopsies for each heart, is read in as data while in the WinBUGS code it was calculated.

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 block and the lines are ended with semi-colons.

The initial values are different in that JAGS does not give initial values for error[] but does give them for p[] and for true[]. In fact the way in which the structured zeros are handled is quite different to the method used by WinBUGS. In JAGS it is not necessary to define error[] in both the data and the initial values files with NA in the appropriate places. In JAGS we can define error[] implicitly and still estimate some of its terms as parameters. This seems to me to be convenient but a little illogical.

Running JAGS from Stata

We can now prepare a Stata program that calls JAGS and fits the biopsy model.

*---------------------------------
* write data file
*---------------------------------
use biopsy.dta, clear
matrix prior = J(4,4,1)
forvalues r=1/4 {
   forvalues c=2/4 {
       if `c' > `r' matrix prior[`r',`c'] = 0
}
 }
 egen nbiops = rsum(grade*)
 gen true = 4
 scalar ns = 157
 wbslist_v2 (scalar ns , format(%5.0f)) (matrix prior , format(%4.0f)) ///
 (structure grade1 grade2 grade3 grade4 , name(biopsies) format(%3.0f)) ///
 (vector nbiops , format(%3.0f) ) ///
 using dataj.txt , replace jags
 *---------------------------------
 * write initial values file
 *---------------------------------
 matrix p = J(1,4,0.25)
 wbslist_v2 (matrix p) (vec true, format(%3.0f) ) ///
 using initj.txt , replace jags
 *---------------------------------
 * write model file
 *---------------------------------
 wbsmodel biopsy_jags.do modelj.txt
 /*
 var
 biopsies[ns,4], # grades observed in ith session (multinomial)
 nbiops[ns], # total number of biopsies in ith session
 true[ns], # true state in ith session
 error[4,4], # error matrix in taking biopsies
 prior[4,4], # prior parameters for rows of error[,]
 p[4]; # underlying incidence of true states
 model {
 for (i in 1:ns){
 true[i] ~ dcat(p[]);
 biopsies[i,] ~ dmulti(error[true[i],],nbiops[i]);
 }
for (j in 1:4) {
 error[j,] ~ ddirch(prior[j,]);
 }
p[] ~ ddirch(prior[4,]); # prior for p
 }
 */
 *---------------------------------
 * write script file
 *---------------------------------
 wbsscript_v2 using scriptj.txt , replace ///
 model(modelj.txt) data(dataj.txt) init(initj.txt) ///
 burnin(1000) update(10000) ///
 set(p error) coda(bsj) jags
 *---------------------------------
 * run and read the results
 *---------------------------------
 wbsrun_v2 using scriptj.txt , jags
 type jagslog.txt
 wbscoda_v2 using bsj , clear jags
 *---------------------------------
 * inspect the results
 *---------------------------------
 mcmcstats *
 mcmctrace p*

Notice that I am using the modified versions of some of the wbs programs, for example, wbsrun_v2. These were introduced in an earlier posting   https://staffblogs.le.ac.uk/bayeswithstata/2014/11/21/jags-with-stata/ and can be downloaded from there.

You might like to look back and compare this program with the version that called WinBUGS that was given in last week’s posting. The only real difference in structure is the use of the jags option in wbslist, sbsscript, wbsrun and wbscoda. All other differences relate to the model file itself and the data and initial values that are supplied and stem from the small differences between WinBUGS and JAGS.

Here are the results calculated by JAGS

----------------------------------------------------------------------------
 Parameter     n  mean    sd    sem median    95% CrI
 ----------------------------------------------------------------------------
 p_1       10000 0.152 0.051 0.0012 0.154 ( 0.046, 0.247 )
 p_2       10000 0.313 0.057 0.0011 0.308 ( 0.213, 0.435 )
 p_3       10000 0.388 0.043 0.0005 0.387 ( 0.306, 0.475 )
 p_4       10000 0.147 0.030 0.0003 0.145 ( 0.093, 0.210 )
 error_1_1 10000 1.000 0.000      . 1.000 ( 1.000, 1.000 )
 error_2_1 10000 0.588 0.067 0.0013 0.589 ( 0.456, 0.715 )
 error_3_1 10000 0.341 0.045 0.0006 0.340 ( 0.258, 0.432 )
 error_4_1 10000 0.099 0.042 0.0005 0.094 ( 0.033, 0.196 )
 error_1_2 10000 0.000 0.000      . 0.000 ( 0.000, 0.000 )
 error_2_2 10000 0.412 0.067 0.0013 0.411 ( 0.285, 0.544 )
 error_3_2 10000 0.036 0.018 0.0002 0.034 ( 0.009, 0.077 )
 error_4_2 10000 0.023 0.023 0.0003 0.015 ( 0.001, 0.087 )
 error_1_3 10000 0.000 0.000      . 0.000 ( 0.000, 0.000 )
 error_2_3 10000 0.000 0.000      . 0.000 ( 0.000, 0.000 )
 error_3_3 10000 0.622 0.047 0.0007 0.624 ( 0.527, 0.710 )
 error_4_3 10000 0.206 0.060 0.0008 0.201 ( 0.103, 0.337 )
 error_1_4 10000 0.000 0.000      . 0.000 ( 0.000, 0.000 )
 error_2_4 10000 0.000 0.000      . 0.000 ( 0.000, 0.000 )
 error_3_4 10000 0.000 0.000      . 0.000 ( 0.000, 0.000 )
 error_4_4 10000 0.672 0.073 0.0010 0.676 ( 0.522, 0.805 )
 ----------------------------------------------------------------------------

Notice that, unlike WinBUGS, JAGS allows us to have items in the error matrix that are not updated. So the results file will contain items like error_1_4 that do not change. In fact, if you look back you will see that the error matrix is not defined in either the data file or the initial values file and so takes it initial values of zero by default when it is listed in the var block.

Both this week’s JAGS code and last week’s WinBUGS code were run for a burn-in of 1,000 followed by 10,000 updates so ought to give comparable results. To make comparison easier I have collected together some of the summary statistics. There is no systematic benefit in using either program.

 --------------------------------------------------
                 WinBUGS              JAGS
 Parameter  mean    sd    sem  mean    sd    sem
 --------------------------------------------------
 error_2_1 0.587 0.066 0.0013 0.588 0.067 0.0013
 error_2_2 0.413 0.066 0.0013 0.412 0.067 0.0013
 error_3_1 0.342 0.046 0.0007 0.341 0.045 0.0006
 error_3_2 0.037 0.018 0.0002 0.036 0.018 0.0002
 error_3_3 0.621 0.048 0.0007 0.622 0.047 0.0007
 error_4_1 0.099 0.042 0.0005 0.099 0.042 0.0005
 error_4_2 0.022 0.023 0.0003 0.023 0.023 0.0003
 error_4_3 0.204 0.061 0.0008 0.206 0.060 0.0008
 error_4_4 0.675 0.073 0.0010 0.672 0.073 0.0010
 p_1       0.153 0.050 0.0011 0.152 0.051 0.0012
 p_2       0.311 0.055 0.0010 0.313 0.057 0.0011
 p_3       0.389 0.044 0.0006 0.388 0.043 0.0005
 p_4       0.147 0.030 0.0003 0.147 0.030 0.0003
 ---------------------------------------------------

Next time I will return to this example and experiment to see how many of the differences between the WinBUGS and JAGS models files are actually essential if we are to get the JAGS program to run.