Code for Basic Hierarchical Bayesian Analysis

Question

I have this code written in winBUGS:

n <- 100

x1 <- rbinom(n,1,.7)

x2 <- rbinom(n,1,.5)

sum(x1)

sum(x2)

model{

x1 ~ dbin(p1, n) x2 ~ dbin(p2, n) p1 ~ dbeta(a1, b1) p2 ~ dbeta(a2,b2)

diff <- p1 - p2 p.value <- step(diff)

} list(n = 100, x1 = 70, x2 = 54, a1 = 1, b1 = 1, a2 = 1, b2 = 1)

I am having trouble on how to execute this in R/JAGS. As a matter of fact, I am not even entirely sure what this code is trying to do (I think to calculate the posteriors?). I have never used winBUGS before, and I am new to R. This is also my first Bayes class, and I am quite lost once the code was introduced.

In addition, how would I calculate the posterior mean and standard deviation for the difference between the proportions? or how to find the posterior probability that p1 is larger than p2, and whether or not that is significant?

Answer 1

If you need help getting rjags set up, I'd expect your TA or classmates could help you with that. This answer focuses on explaining what the code does.

For pedagogical purposes, let's impose a simple narrative on this data. Say we have two coins and we want to know if one is more likely to flip a head than the other. The code is in two parts.

Generative/Forward Simulation

The first is a forward simulation, where we already know what the probability of flipping a head is for the two coins.

## how many times we will flip each coin
n <- 100

## coin 1 has 70% chance to land on heads
## simulate n flips
x1 <- rbinom(n, 1, .7)

## coin 2 has 50% chance to land on heads
## simulate n flips
x2 <- rbinom(n, 1, .5)

## count number of coin 1 heads
sum(x1) # 70

## count number of coin 2 heads
sum(x2) # 54

Posterior Probabilities

Now we take this simulated data and try to reverse the experiment. That is, if we observed 70 heads from one coin and 54 heads from another, can we say something about the probabilities that each coin has to flip a head? Specifically, the code is going to ask, "Under a uniform assumption about the true probability that each coin has, what is the (posterior) probability that coin 1 is more likely to flip a head than coin 2?"

Computing Likelihoods

The way that JAGS will do this is with MCMC, where samples are taken from the space of all possible parameter configurations (in this case, values of p1 and p2), weighted by how likely each configuration is to generate the observed data. Hence, the main thing JAGS code needs to accomplish is define how to generate parameter values, and how to compute the likelihood of the data, given those values. This is the first part of the model code.

Computing Values of Interest

The second part of the code in the model section is testing the question we posed. This is done by computing some additional variables that don't contribute to the likelihood of the configurations. Instead they tell us information about the configurations being sampled. Specifically, diff will keep track of what the distribution of differences between the two coin probabilities are, and p.value will track how often p1 is greater than p2.

model{

## COMPUTE LIKELIHOODS
## compute likelihood that heads resulted from coins
## with given probabilities after `n` coin flips
x1 ~ dbin(p1, n) 
x2 ~ dbin(p2, n) 

## SAMPLE PARAMETERS
## randomly select coin probabilities from Beta distributions
## Note that since these are all 1, it's really just a Uniform[0,1]
p1 ~ dbeta(a1, b1)
p2 ~ dbeta(a2, b2)

## HYPOTHESIS TEST
## compute how often coin1's probability is greater than coin2's
diff <- p1 - p2 
p.value <- step(diff)
} 

## INPUT VALUES
list(n = 100, x1 = 70, x2 = 54, a1 = 1, b1 = 1, a2 = 1, b2 = 1)