How to specify a Bayesian binomial model with shrinkage to the population?

Question

Problem

I’m currently working on a problem where I have count data for $n$ items in the following form:

    Item   Count   Total
    --------------------
    1          1      30
    2         10     340
    3          0       0
    4          1      10
             ...

I want to calculate the posterior of the success probability $p_i$ for each item. Additionally, I suspect that all the items have a similar or even the same success probability but I want the model to determine that.

My expectations for the resulting model are:

• Estimate the underlying global parameter of all the items and
• Shrink the success probabilities $p_i$ towards this global parameter in order to deal with items which have few observations.
• Items with observations indicating a different success probability should be able to “override” the global success probability.

Stan code

This is the Stan code I implemented so far:

    data {
      int<lower=0> J; // amount of items
      int count[J];
      int total[J];
    }
    parameters {
      real<lower=0,upper=1> mu;
      real<lower=0.001> ss;
      real<lower=0,upper=1> prob[J];
    }
    model {
      mu ~ uniform(0, 1);
      ss ~ gamma(0.01, 0.01);
      prob ~ beta(mu*ss, (1-mu)*ss);
      count ~ binomial(total, prob);
    }

Here I achieve the shrinkage by assuming that the success probabilities of the items are sampled from a beta distribution (parametrized using mean and sample size) where the mean is the global success probability.

The problem occurs when specifying the prior of the sample size. If I choose a very vague prior, Stan requires a lot of iterations until the effective sample size is large enough. Also I’m not sure if this is the right approach since it is not clear if the amount of shrinkage this model provides is always appropriate.

I’m interested if there is an elegant way to implement my expected behavior or if there is some flaw in my overall model.

Answer 1

Try this: $p_i\stackrel{iid}{\sim} Be(\alpha,\beta)$ and $p(\alpha,\beta)\propto (\alpha+\beta)^{-5/2}$.

I believe the issue you are running into is that $ss\sim Ga(\gamma,\gamma)$ with $\gamma\to 0$ results in the improper $p(ss)\propto 1/ss$ prior. This prior, together with your uniform prior on mu, results in an improper posterior. Despite the fact that you are using a proper prior, it is close enough to this improper prior to cause issues.

You might want to take a look at page 110 of Bayesian Data Analysis (3rd ed) for a discussion of priors for this model as well as the prior suggested above.

Answer 2

This Blog entry here http://lingpipe-blog.com/2009/09/23/bayesian-estimators-for-the-beta-binomial-model-of-batting-ability/ shows another possibility for modelling the prior on the parameters of the Beta population distribution, namely

$\frac{\alpha}{\alpha + \beta}$ ~ Uniform(0,1) = Beta(1,1)

$(a+b) \sim $ Pareto(1.5, 1)

Here is another resource regarding your question (slide 5): http://www.stat.cmu.edu/~brian/724/week06/lec15-mcmc.pdf.