Germán Rodríguez
Multilevel Models Princeton University

Running WinBUGS in Batch Mode

Pointing and clicking following the instructions in a compound document such as hospBUGs is fine for becoming familiar with the software, but for serious work you will want to create a script that documents all the steps.

Model and Data Files

This turns out to be relatively easy. You need to save the model into a text file. For the hospital data I created hospModel.txt with the following contents

model {
    # N observations
    for (i in 1:N) {
        hospital[i] ~ dbern(p[i])
        logit(p[i]) <- bcons + blonginc*loginc[i] + bdistance*distance[i] + 
            bdropout*dropout[i] + bcollege*college[i] + u[group[i]] 
    }
    # M groups
    for (j in 1:M) {
        u[j] ~ dnorm(0,tau)
    }
    # Priors
    bcons     ~ dnorm(0.0,1.0E-6)
    bloginc   ~ dnorm(0.0,1.0E-6)
    bdistance ~ dnorm(0.0,1.0E-6)
    bdropout  ~ dnorm(0.0,1.0E-6)
    bcollege  ~ dnorm(0.0,1.0E-6)
    # Hyperprior
    tau ~ dgamma(0.001,0.001)
}

Next I saved the list with the number of units per level in a file called hospList.txt which looks like this:

list(N=1060,M=501)

The data themselves, of course, are already in a file called hospital.txt. Click here to see it.

At this point we add a file with the initial values, which I called hospInits.txt and which looks like this:

list(bcons=-2.69517150,bloginc=0.45852520,bdistance=-0.07627492,
	bdropout=-1.57007674,bcollege=0.82072827,tau=1)

The Script File

Wih these four files, all we need is the following script:

display('log')
check('d:/pop510/hospModel.txt')
data('d:/pop510/hospList.txt')
data('d:/pop510/hospital.txt')
compile(1)
inits(1,'d:/pop510/hospInits.txt')
gen.inits()
update(500)
set(bcons)
set(bloginc)
set(bdistance)
set(bdropout)
set(bcollege)
set(tau)
update(5000)
coda(*,'d:/pop510/hospOut')

Each of the commands listed above corresponds to one of the actions we performed before: checking the model, reading the data, compiling, reading initial values for the parameters, generating initial values for the random effects, running a burn-in, setting the nodes to be traced, updating, and finally writing the output to a file for further analysis.

Note that all file names are specified to WinBUGS using forwards slashes, although the path separator in Windows is a backslash. That's why we write d:/pop510 rather than d:\pop510.

You can open the script file in WinBUGS (make sure you say it is of type "text" before you open it) and then select Model | Script in the menu or press Alt-M-C, or even better, press Start and type cmd in the search box to open a command window, and then type

D:\WinBUGS\WInBUGS14.exe /PAR d:/pop510/hospScript.txt

Make sure you type the path and name of the WinBUGS executable in your system (I installed my copy in a folder called D:\WinBUGS. The executable itself is called WinBUGS14.exe) and the name of the script file.

The CODA Output

CODA is a popular program for convergence diagnostics written for R. WinBUGS produces output that can be read by CODA in a simple format. The command coda(*,'d:/pop510/hospOut') in our script generates two files: an index and a data file.

The index file has a name consisting of the stem we specified and the postfix "Index.txt", so in our case the file is called "D:\pop510\hospOutIndex.txt" and looks like this

bcons     1        5000
bloginc   5001    10000
bdistance 10001   15000
bdropout  15001   20000
bcollege  20001   25000
tau      25001   30000

The data file has a name consisting of the stem we specified and the postfix "1.txt" where 1 stands for chain 1, so in our case the file is called :D:\pop510\hospOut1.txt" and looks like this (I show only the first two lines and the last one)

501     -3.038
502     -2.795
...
5500    0.7913

This file simply has the sample number and the value of the parameter for each draw, in the order specified in the index file. In our case lines 1 to 5000 are the values of bcons, lines 5001 to 10000 are the values of bloginc, and so on.

Obviously reading this file into Stata is not hard, but a user-contributed command makes it very easy.

Running WinBUGS from Stata

Thompson, Palmer and Moreno (2006) "Bayesian Analysis in Stata using WinBUGS", The Stata Journal,6:530-549, have written a series of commands to facilitate using WinBUGS from Stata.

The interface commands are rather low level. You need to prepare external text files with a script as outlined in the previous section, the model in WinBUGS language, the data using lists and vectors, and the initial values.

Two helpers are the commands wbarray to write data in array format, and wbvector to write a list. The command wbrun is a shortcut for running WinBUGS with a script.

Importing CODA Output

Once WinBUGS has ran the command wbcoda is extremely useful for importing the samples into Stata. All you need to do is specify the root name given to the coda() command in the script.

Here's all we need to (1) run WinBUGS from Stata using the files prepared in the previous section and (2) import the results. (I have modified wbrun to default to the path to the WinBUGS14 executable in my system.)

. wbrun, script(D:/pop510/hospScript.txt)

. wbcoda, root(D:/pop510/hospOut)

. sum bcons-tau

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       bcons |      5000   -3.516448    .4539867     -5.264     -2.152
     bloginc |      5000    .5862913    .0694852      .3525       .813
   bdistance |      5000   -.0776984    .0332205     -.1857     .03305
    bdropout |      5000   -2.054222    .2577039     -3.099      -1.15
    bcollege |      5000    1.058578    .3996331     -.2771      2.754
-------------+--------------------------------------------------------
         tau |      5000    .6226772     .207884      .2314      2.256

We'll print a different summary below.

Convergence Diagnostics

Once the samples are in Stata you can analyze them using standard tools, plus a few specialized commands. We can obtain trace plots for all the parameters using the wbtrace command

. wbtrace bcons-tau, cgoptions(title(Trace Plots for Hospital Data))
. graph export hospTrace.png, width(500) replace

The command wbsection Divides a chain into subsections and plots a smooth density estimate for each section. If the chain has converged the superimposed densities will be similar.

. wbsection bcons-tau, m(5) cgoptions(title(Density Plots for Hospital Data))

. graph export hospDen.png, width(500)

These graphs would like much tighter if we droppped the first couple of thousand samples. Convergence was particularly slow for bloginc.

The command wbgeweke performs a test comparing means in early and later parts of a chain, and wbbgr produces the Brooks-Gelman-Rubin plot that compares parallel chains. Two other commands are wbac for autocorrelations and wbintervals to plot intervals by section. Here are the Geweke tests for our run, raising no concerns

. wbgeweke bcons-tau

Parameter: bcons first 10.0% (n=500) vs last 50.0% (n=2500)
Means (se)    -3.5872 (    0.1041)   -3.5297 (    0.0649)
Autocorrelations  0.9518   0.9633
Mean Difference (se)    0.0575 (    0.1227) z =  0.468 p =  0.6395

Parameter: bloginc first 10.0% (n=500) vs last 50.0% (n=2500)
Means (se)     0.5967 (    0.0170)    0.5869 (    0.0092)
Autocorrelations  0.9519   0.9605
Mean Difference (se)   -0.0098 (    0.0193) z =  0.510 p =  0.6099

Parameter: bdistance first 10.0% (n=500) vs last 50.0% (n=2500)
Means (se)    -0.0758 (    0.0035)   -0.0775 (    0.0016)
Autocorrelations  0.7105   0.7108
Mean Difference (se)   -0.0017 (    0.0039) z =  0.441 p =  0.6590

Parameter: bdropout first 10.0% (n=500) vs last 50.0% (n=2500)
Means (se)    -2.0415 (    0.0240)   -2.0471 (    0.0133)
Autocorrelations  0.6766   0.7410
Mean Difference (se)   -0.0057 (    0.0274) z =  0.206 p =  0.8366

Parameter: bcollege first 10.0% (n=500) vs last 50.0% (n=2500)
Means (se)     1.0329 (    0.0268)    1.0589 (    0.0141)
Autocorrelations  0.4710   0.5187
Mean Difference (se)    0.0260 (    0.0303) z =  0.859 p =  0.3906

Parameter: tau first 10.0% (n=500) vs last 50.0% (n=2500)
Means (se)     0.6308 (    0.0400)    0.6128 (    0.0232)
Autocorrelations  0.9372   0.9503
Mean Difference (se)   -0.0180 (    0.0462) z =  0.390 p =  0.6962

The command wbstats computes a standard BUGS summary for each chain. Two other commands for summarizing data are wbdensity and wbhull.

. wbstats bcons-tau
Parameter        n    mean      sd      se  median         95% CrI
bcons         5000  -3.516   0.454  0.0476  -3.504 (  -4.416,  -2.655 )
bloginc       5000   0.586   0.069  0.0072   0.584 (   0.458,   0.721 )
bdistance     5000  -0.078   0.033  0.0012  -0.078 (  -0.143,  -0.011 )
bdropout      5000  -2.054   0.258  0.0095  -2.044 (  -2.593,  -1.580 )
bcollege      5000   1.059   0.400  0.0102   1.047 (   0.269,   1.876 )
tau           5000   0.623   0.208  0.0202   0.591 (   0.326,   1.106 )

The mean is not too different from median except for the precision, which has a more skewed posterior distribution. If you want to report the standard deviation of the random effect I suggest you calculate sigma=1/sqrt(tau) and summarize that. YOu'll get a mean of 1.315, not quite the same as 1/sqrt(0.623)=1.267. The median and credible interval are not affected by transformations.