# Hierarchical Bayesian modeling of incidence rates

Kevin Murphy's book discusses a classical Hierarchical Bayesian problem (originally discussed in Johnson and Albert, 1999, p24):

Suppose that we are trying to estimate the cancer rate in $$N$$ cities. In each city, we sample a number of individuals $$N_i$$ and measure the number of people with cancer $$x_i \sim \text{Bin}(N_i, \theta_i)$$, where $$\theta_i$$ is the true cancer rate in the city.

We would like to estimate the $$\theta_i$$'s while allowing the data-poor cities to borrow statistical strength from data-rich cities.

To do so, he models $$\theta_i \sim \text{Beta}(a,b)$$ so that all cities share the same prior, so the final models look as follows:

$$p(\mathcal{D}, \theta, \eta|N)=p(\eta)\prod\limits^N_{i=1}\text{Bin}(x_i|N_i, \theta_i)\text{Beta}(\theta_i|\eta)$$

where $$\eta = (a,b)$$.

The crucial part about this model is of course (I quote), "that we infer $$\eta=(a,b)$$ from the data, since if we just clamp it to a constant, the $$\theta_i$$ will be conditionally independent, and there will be no information flow between them".

I am trying to model this in PyMC, but as far as I understand, I need a prior for $$a$$ and $$b$$ (I believe this is $$p(\eta)$$ above).

What would be one good prior for this model?

In case it helps, the code, as I have it now is:

bins = dict()
ps   = dict()
for i in range(N_cities):
ps[i]   = pm.Beta("p_{}".format(i), alpha=a, beta=b)
bins[i] = pm.Binomial('bin_{}'.format(i), p=ps[i],n=N_trials[i],  value=N_yes[i], observed=True)

mcmc = pm.MCMC([bins, ps])


A similar problem is discussed in Gelman, Bayesian Data Analysis, (2nd ed, p. 128; 3rd edition p. 110). Gelman suggests a prior $p(a,b)\propto (a+b)^{-5/2}$, which effectively constrains the "prior sample size" $a+b$, and therefore the beta hyperprior is not likely to be highly informative on its own. (As the quantity $a+b$ grows, the variance of the beta distribution shrinks; in this case, smaller prior variance constrains the "weight" of the observed data in the posterior.) Additionally, this prior does not set whether $a>b$, or the opposite, so appropriate distributions of pairs of $(a,b)$ are inferred from all the data together, as you would prefer in this problem.
Gelman also suggests reparameterizing the model in terms of the logit of the mean of $\theta$ and the "sample size" of the prior. So instead of doing inference on $(a,b)$ directly, the problem is about inference on the transformed quantities $\text{logit}\left(\frac{a}{a+b}\right)$ and $\log(a+b)$. This admits transformed prior values on the real plane, rather than untransformed prior values that must be strictly positive. Also, this accomplishes a posterior density that is more diffuse when plotted. This makes the accompanying graphs more legible, which I find helpful.