:::{post} Jan 10, 2023 :tags: generalized linear model, hierarchical model :category: intermediate :author: Demetri Pananos, Junpeng Lao, Raúl Maldonado, Farhan Reynaldo :::
This short tutorial demonstrates how to use PyMC to do inference for the rat tumour example found in chapter 5 of Bayesian Data Analysis 3rd Edition {cite:p}gelman2013bayesian. Readers should already be familiar with the PyMC API.
Suppose we are interested in the probability that a lab rat develops endometrial stromal polyps. We have data from 71 previously performed trials and would like to use this data to perform inference.
The authors of BDA3 choose to model this problem hierarchically. Let $y_i$ be the number of lab rats which develop endometrial stromal polyps out of a possible $n_i$. We model the number rodents which develop endometrial stromal polyps as binomial
allowing the probability of developing an endometrial stromal polyp (i.e. $\theta_i$) to be drawn from some population distribution. For analytical tractability, we assume that $\theta_i$ has Beta distribution
We are free to specify a prior distribution for $\alpha, \beta$. We choose a weakly informative prior distribution to reflect our ignorance about the true values of $\alpha, \beta$. The authors of BDA3 choose the joint hyperprior for $\alpha, \beta$ to be
For more information, please see Bayesian Data Analysis 3rd Edition pg. 110.
Our joint posterior distribution is
\propto p(\alpha, \beta) p(\theta \lvert \alpha,\beta) p(y \lvert \theta)$$which can be rewritten in such a way so as to obtain the marginal posterior distribution for $\alpha$ and $\beta$, namely
p(\alpha, \beta) \prod_{i = 1}^{N} \dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \dfrac{\Gamma(\alpha+y_i)\Gamma(\beta+n_i - y_i)}{\Gamma(\alpha+\beta+n_i)}$$See BDA3 pg. 110 for a more information on the deriving the marginal posterior distribution. With a little determination, we can plot the marginal posterior and estimate the means of $\alpha$ and $\beta$ without having to resort to MCMC. We will see, however, that this requires considerable effort.
The authors of BDA3 choose to plot the surface under the parameterization $(\log(\alpha/\beta), \log(\alpha+\beta))$. We do so as well. Through the remainder of the example let $x = \log(\alpha/\beta)$ and $z = \log(\alpha+\beta)$.
The plot shows that the posterior is roughly symmetric about the mode (-1.79, 2.74). This corresponds to $\alpha = 2.21$ and $\beta = 13.27$. We can compute the marginal means as the authors of BDA3 do, using
\sum_{x,z} \alpha p(x,z\lvert y) $$\sum_{x,z} \beta p(x,z\lvert y) $$Computing the marginal posterior directly is a lot of work, and is not always possible for sufficiently complex models.
On the other hand, creating hierarchical models in PyMC is simple. We can use the samples obtained from the posterior to estimate the means of $\alpha$ and $\beta$.
We can plot a kernel density estimate for $x$ and $y$. It looks rather similar to our contour plot made from the analytic marginal posterior density. That's a good sign, and required far less effort.
From here, we could use the trace to compute the mean of the distribution.
Analytically calculating statistics for posterior distributions is difficult if not impossible for some models. PyMC provides an easy way drawing samples from your model's posterior with only a few lines of code. Here, we used PyMC to obtain estimates of the posterior mean for the rat tumor example in chapter 5 of BDA3. The estimates obtained from PyMC are encouragingly close to the estimates obtained from the analytical posterior density.
:::{bibliography} :filter: docname in docnames :::
gelman2013bayesian by Demetri Pananos and Junpeng Lao (pymc#3054):::{include} ../page_footer.md :::