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.

## Recent Comments