Title :¶

MCMC and PyMC3 Hierarchical Models

Description :¶

Download this notebook on your local environment. Some parts of it are time-consuming, so we are providing the output cells as well. Once you have gone through them, make changes and run the cells again.

CS109B Data Science 2: Advanced Topics in Data Science

Lecture 13 - MCMC and Hierarchical Models in PyMC3¶

Harvard University
Spring 2021
Instructors: Pavlos Protopapas, Mark Glickman, and Chris Tanner
Additional Instructor: Eleni Angelaki Kaxiras
Content: Eleni Angelaki Kaxiras

---

Learning Objectives¶

By the end of this lab, you should be able to understand how to:

• run a PyMC3 bayesian model using MCMC.
• implement hierarchical models in PyMC3.

Table of Contents¶

1. Markov Chain Monte Carlo (MCMC) Simulations in PyMC3.
2. Hierarchical models.

Top

1. Markov Chain Monte Carlo (MCMC) Simulations¶

PyMC3 is a Python library for programming Bayesian analysis, and more specifically, data creation, model definition, model fitting, and posterior analysis. It uses the concept of a model which contains assigned parametric statistical distributions to unknown quantities in the model. Within models we define random variables and their distributions. A distribution requires at least a name argument, and other parameters that define it. You may also use the logp() method in the model to build the model log-likelihood function. We define and fit the model.

PyMC3 includes a comprehensive set of pre-defined statistical distributions that can be used as model building blocks. Although they are not meant to be used outside of a model, you can invoke them by using the prefix pm, as in pm.Normal.

Monte Carlo methods are rooted in the 1940s when Enrico Fermi, John von Neumann, Stan Ulam, Nicolas Metropolis, and others, employed the use of random numbers to study physics from a stochastic point of view. The name is attributed to Metropolis 1.

For a history of the Monte Carlo Method first used in Physics see: The beginning of the Monte Carlo Method

PyMC3 uses the No-U-Turn Sampler (NUTS) and the Random Walk Metropolis, two Markov chain Monte Carlo (MCMC) algorithms for sampling in posterior space. Monte Carlo gets its name because when we sample in posterior space, we choose our next move via a pseudo-random process. NUTS is a sophisticated algorithm that can handle a large number of unknown (albeit continuous) variables.

Bayesian Linear Regression¶

We will artificially create the data to predict on. We will then see if our model predicts them correctly.

Now let's see if our model will correctly predict the values for our unknown parameters, namely $b_0$, $b_1$, $b_2$ and $\sigma$.

Defining the Problem¶

Our problem is the following: we want to perform multiple linear regression to predict an outcome variable $Y$ which depends on variables $\bf{x}_1$ and $\bf{x}_2$.

We will model $Y$ as normally distributed observations with an expected value $mu$ that is a linear function of the two predictor variables, $\bf{x}_1$ and $\bf{x}_2$.

\begin{equation} Y \sim \mathcal{N}(\mu,\,\sigma^{2}) \end{equation}

\begin{equation} \mu = \beta_0 + \beta_1 \bf{x}_1 + \beta_2 x_2 \end{equation}

where $\sigma^2$ represents the measurement error (in this example, we will use $\sigma^2 = 10$)

We also choose the parameters to have normal distributions with those parameters set by us.

\begin{eqnarray} \beta_i \sim \mathcal{N}(0,\,10) \\ \sigma^2 \sim |\mathcal{N}(0,\,10)| \end{eqnarray}

Defining a Model in PyMC3¶

Note: If our problem was a classification for which we would use Logistic regression see below

Python Note: pm.Model is designed as a simple API that abstracts away the details of the inference. For the use of with see Compounds statements in Python..

Fitting the Model with Sampling - Doing Inference¶

See below for PyMC3's sampling method. As you can see it has quite a few parameters. Most of them are set to default values by the package. For some, it's useful to set your own values.

pymc3.sampling.sample(draws=500, step=None, n_init=200000, chains=None,
                      cores=None, tune=500, random_seed=None)

Parameters to set:

• draws: (int): Number of samples to keep when drawing, defaults to 500. Number starts after the tuning has ended.
• tune: (int): Number of iterations to use for tuning the model, also called the burn-in period, defaults to 500. Samples from the tuning period will be discarded.
• target_accept (float in $[0, 1]$). The step size is tuned such that we approximate this acceptance rate. Higher values like 0.9 or 0.95 often work better for problematic posteriors.
• (optional) chains (int) number of chains to run in parallel, defaults to the number of CPUs in the system, but at most 4.

pm.sample returns a pymc3.backends.base.MultiTrace object that contains the samples. We usually name it a variation of the word trace. All the information about the posterior is in trace, which also provides statistics about the sampler.

How does a good trace plot look like? Is this a good one?¶

Model Plotting¶

PyMC3 provides a variety of visualizations via plots: https://docs.pymc.io/api/plots.html. arviz is another library that you can use.

$\hat{R}$ is a metric for comparing how well a chain has converged to the equilibrium distribution by comparing its behavior to other randomly initialized Markov chains. Multiple chains initialized from different initial conditions should give similar results. If all chains converge to the same equilibrium, $\hat{R}$ will be 1. If the chains have not converged to a common distribution, $\hat{R}$ will be > 1.01. $\hat{R}$ is a necessary but not sufficient condition.

For details on the $\hat{R}$ see Gelman and Rubin (1992).

This linear regression example is from the original paper on PyMC3: Salvatier J, Wiecki TV, Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2:e55 https://doi.org/10.7717/peerj-cs.55

Top

2. Hierarchical Models¶

Gelman et al.'s famous radon dataset is a classic for hierarchical modeling. In this dataset the amount of the radioactive gas radon has been measured among different households in counties in several states. Radon gas is known to be the highest cause of lung cancer in non-smokers. It is believed to be more strongly present in households containing a basement and to differ in amount present among types of soil. Here we'll investigate this differences and try to make predictions of radonlevels in different county's based on the county itself and the presence of a basement. We’ll look at Minnesota, a state that contains 85 counties in which different measurements are taken, ranging from 2 to 116 measurements.

Each row in the dataframe represents the radon measurements for one house in a specific county including whether the house has a basement (floor = 0) or not (floor = 1). We are interested in whether having a basement increases the radon measured in the house.

To keep things simple let's keep only the following three variables: county, log_radon, and floor

Let's check how many different counties we have. We also notice that they have a different number of houses. Some have a large number of houses measured, some only 1.

Pooling: Same Linear Regression for all¶

We can just pool all the data and estimate one big regression to asses the influence of having a basement on radon levels across all counties. Our model would be:

\begin{equation} y_{i} = \alpha + \beta*floor_{i} \end{equation}

Where $i$ represents the measurement (house), and floor contains a 0 or 1 if the house has a basement or not. By ignoring the county feature, we do not differenciate on counties.

Remember, with the pooled model we have only one intercept, $\alpha$, and only one slope, $\beta$ for all the counties. Let's plot the regression lines.

Unpooling: Separate Linear Regression for each county¶

We believe that different counties have different relationships of radon and basements. Our model would be:

\begin{equation} \textbf{radon}_{i,c} = \alpha_c + \beta_c*\textbf{floor}_{i,c} \end{equation}

Where $\textbf{i}$ represents the measurement, $\textbf{c}$ the county, and floor contains a 0 or 1 if the house has a basement or not.

Notice we have separate coefficients for each county in $a_c$ and $b_c$. They are totally different, they could even come from different distributions.

We will do this for only one county, since this is very time consuming, as an example.

Print the regression line for our chosen county alone.

Partial pooling: Hierarchical Regression (Varying-Coefficients Model)¶

Counties, of course, have similarities, so there is a middle ground to both of these extremes. Specifically, we may assume that while $\alpha_c$ and $\beta_c$are different for each county as in the unpooled case, the coefficients are all drawn from the same distribution:

\begin{equation} radon_{i,c} = \alpha_c + \beta_c*floor_{i,c} \end{equation}

\begin{equation} a_c \sim \mathcal{N}(\mu_a,\,\sigma_a^{2}) \end{equation}

\begin{equation} b_c \sim \mathcal{N}(\mu_b,\,\sigma_b^{2}) \end{equation}

where the common parameters are: \begin{eqnarray} \mu_a \sim \mathcal{N}(0,\,10) \\ \sigma_a^2 \sim |\mathcal{N}(0,\,10)| \\ \mu_b \sim \mathcal{N}(0,\,10) \\ \sigma_b^2 \sim |\mathcal{N}(0,\,10)| \end{eqnarray}

The different counties are effectively sharing information through the common priors. We are thus observing what is known as shrinkage; modeling the groups not as independent from each other, neither as a single group but rather as related.

TALKING POINTS¶

We saw that some counties had only one sample, so if that house is a really old with old lead pipes, our prediction will be that all houses in this county have radon. On the other extreme, if we have a newer house with no radon then again we will have missleading results. In one case, you will overestimate the bad quality and in the other underestimate it. Under a hierarchical model, the miss-estimation of one group will be offset by the information provided by the other groups. As always, gathering more data helps, if this is an option.

Inference¶

Adapted from PyMC3 docs: https://docs.pymc.io/notebooks/GLM-hierarchical.html

Appendix A: Is this a fair coin?¶

Can we find out if this a fair coin without using sampling?¶

Let's say you visit the casino. You want to test your theory that casinos are dubious places where coins have been manipulated to have a larger probability for tails. So you will try to estimate how fair a coin is based on a certain amount of flips. You have no prior opinion on the coin's fairness (i.e. what $\theta$ might be), and begin flipping the coin. You get either Heads ($H$) or Tails ($T$) as our observed data and want to see if your posterior probabilities change as you obtain more data, that is, more coin flips. A nice way to visualize this is to plot the posterior probabilities as we observe more data. Your data is the pair ($n,k$) where $n$ is the number of flips and $k$ is the number of successes.

We will be using Bayes rule. $\textbf{D}$ is our data.

\begin{equation} P(\theta|\textbf{D}) \propto P(\theta) \times P(\textbf{D} |\theta) \end{equation}

a) Let's say we believe that our experiment is governed by a $\text{Binomial}$ distribution (this will be our Likelihood).

\begin{equation} P(\theta|\textbf{D}) \propto \text{Beta(a,b)} \times \text{Binomial}(\textbf{k},\textbf{n}|\theta) \end{equation}

b) We start with a non-informative prior, a $\text{Beta}$ distribution with (a=b=1.)

\begin{equation} P(\theta|1.,1.) = \text{Beta(1., 1.)} \end{equation}

c) We enter the data as pairs of ($\textbf{k},\textbf{n}$) (we observe $\textbf{k}$ heads in $\textbf{n}$ tosses), and we update our Beta with new a,b as follows $^*$:

\begin{equation} P(\theta|\textbf{k}) = Beta(\alpha + \textbf{k}, \beta + (n - \textbf{k})) \end{equation}

d) We are done! That is why using a congugate prior helps you find your parameters analytically.

$^*$ (the proof for this formula is beyond our scope, if interested, see this Wikipedia article)

Appendix B: MAP¶

Fitting the Model with MAP (FYI, we will not directly use this method)¶

In Bayesian analysis we have our prior(s), we define our likelihood, and, having specified our model, we try to calculate posterior estimates for the unknown variables in the model. We could try to calculate the posterior estimates analytically, but for most the models, this is not feasible. What we do then is compute summaries based on samples drawn from the posterior distribution using Markov Chain Monte Carlo (MCMC) sampling methods.

\begin{equation} P(\theta|\textbf{D}) \rightarrow \{\theta_1,....\theta_n\} \end{equation}

Then we can find any estimate we want by using these samples, for example:

\begin{equation} \mathbb{E}[f(\theta] = \int d\theta{p(\theta) f(\theta)} \end{equation}

So we calculate the maximum a posteriori (MAP) point using optimization methods.

\begin{equation} f(\hat{\theta}), \hat{\theta} = argmax ({p(\theta))} \end{equation}

The maximum a posteriori (MAP) estimate for a model, is the mode of the posterior distribution and is generally found using numerical optimization methods. This is often fast and easy to do, but only gives a point estimate for the parameters and can be biased if the mode isn’t representative of the distribution. PyMC3 provides this functionality with the find_MAP function.

MAP estimate is not always reasonable, especially if the mode is at an extreme or we have a multimodal distribution, or we have high dimensional posteriors. This will often occur in hierarchical models with the variance parameter for the random effect. If the individual group means are all the same, the posterior will have near infinite density if the variance parameter for the group means is almost zero. Most techniques for finding the MAP estimate only find a local optimium (which is often good enough), and can therefore fail badly for multimodal posteriors, as mentioned above.

To solve these issues we turn to sampling as our method for finding the posterior.

You do not have to worry about MAP in our problems. Our pyMC3 models use the MAP method to initialize the variables under the hood and we do not have to explicitly set this.

References:

• Salvatier J, Wiecki TV, Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2:e55 (https://doi.org/10.7717/peerj-cs.55)
• Distributions in PyMC3
• More Details on Distributions

Information about PyMC3 functions including descriptions of distributions, sampling methods, and other functions, is available via the help command.