Bayesian Inference Using OpenBUGS

We will use the data set survey for our first demonstration of OpenBUGS. Although the example is elementary, it does contain all the essential steps. There are more advanced examples along with necessary background materials in the R Tutorial eBook.

The central concept of OpenBUGS is the BUGS model. As we will prove, it is not always necessary to create a BUGS model from scratch. Instead, we can build our models incrementally from simple ones. Hence our first task is to create our own library of basic BUGS models that we can reuse later.

The BUGS language bears a strong resemblance to R. We will introduce more BUGS syntax as we move along. If lack of patience, there is full detail in the WinBUGS online manual.

We begin with introducing the “~” operator, which describes the probability distribution of a random variable. For example, the following indicates that a random variable y fits a binomial distribution with probability of success p and size N.

Problem

The data set survey contains sample smoker statistics among university students. Denote the proportion of smokers in the general student population by p. With uniform prior, find the mean and standard deviation of the posterior of p using OpenBUGS.

Solution

The first step is to create a BUGS model. As we portrait the number of student smokers y as the outcome of a binomial experiment of size N and success probability p, we write the following in BUGS.

We assume the prior distribution of p to be uniform between 0 and 1, which is equivalent to Beta(1,1).

In summary, we have the following BUGS model. The sole purpose of the model function below is to serve as a packaging wrapper for OpenBUGS. The function itself is meaningless to R.

To transfer the model to OpenBUGS, we load the R2OpenBUGS extension and write the model to a temporary location using the method write.model. We denote the model file location by model.file.

Then we have to decide data parameters of the BUGS model. We find that there are 236 students in the survey, and 47 of them smoke, which we denote by N and y respectively.

We then identify data variables in a list called data.

And we identify the variable p to be monitored in a vector called params.

Lastly, we need to select some initial parameters for the simulation. A rule of thumb is to choose values as close to the expected result as possible. In this case, we initialize p to be 0.5. Notice how we wrap the initial values inside a list that is to be returned by a function.

Then we invoke OpenBUGS with the namesake method bugs and save the result in a variable out. We select 10,000 iterations per simulation chain.

Upon completion of the bugs method, we should check the Rhat component of the output. If we find any value higher than the 1.1 threshold, we should increase the n.iter argument in the bugs method and try again. However, even if all values in Rhat are less than 1.1, it still does not guarantee convergence. We shall demonstrate later how to perform further analysis using CODA.

Finally, we can retrieve the posterior mean and standard deviation of p from the output.

It is informative to print out the simulation result in full detail at this point.

MCMC Convergence

To use CODA for analyzing the MCMC convergence, we have to enable the codaPkg option, which allows us to convert the output for CODA using read.bugs.

Then we can print out the trace plot with xyplot in the coda extension and see if the simulation values stabilize.

Next, we should inspect if the density plot is well-defined.

We also need to check if the auto-correlation of the time series converge to zero using acfplot.

Lastly, we should perform the Gelman-Rubin convergence diagnostic with gelman.diag. The shrink factors should be below 1.05.

Actually, it is prudent to create the Gelman-Rubin-Brooks plot for visual confirmation of the shrink factor convergence as well.

With MCMC convergence assured, we can retrieve the point estimates and 95% credible intervals of p.