11. Regression with MCMC

Data set download

[1]:

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade iqplot colorcet bebi103 arviz cmdstanpy watermark"
    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    stdout, stderr = process.communicate()
    import cmdstanpy; cmdstanpy.install_cmdstan()
    data_path = "https://s3.amazonaws.com/bebi103.caltech.edu/data/"
else:
    data_path = "../data/"
# ------------------------------

import pandas as pd

import cmdstanpy
import arviz as az

import iqplot
import bebi103

import holoviews as hv
hv.extension('bokeh')
bebi103.hv.set_defaults()

import bokeh.io
bokeh.io.output_notebook()
Loading BokehJS ...

In a previous lesson, we plotted posterior distributions of models for microtubule size as a function of the size of droplets encapsulating them using data from Good, et al., Science, 342, 856-860, 2013. Here, we will perform a similar analysis using MCMC. We start by loading in the data set and making a quick plot.

[2]:

# Load in Data Frame
df = pd.read_csv(os.path.join(data_path, 'good_invitro_droplet_data.csv'), comment='#')

hv.Scatter(
    data=df,
    kdims='Droplet Diameter (um)',
    vdims='Spindle Length (um)',
)

Data type cannot be displayed:

[2]:

Updated generative model

We will estimate parameters for the conserved tubulin model. In this model, the spindle length \(l\) is a function of droplet diameter \(d\) according to

\begin{align} l(d;\gamma, \phi) = \frac{\gamma d}{\left(1+(\gamma d/\phi)^3\right)^{\frac{1}{3}}}. \end{align}

We will model the variation off of this curve as Normally distributed (assuming homoscedasticity), such that our likelihood is defined as

\begin{align} l_i \mid d_i, \gamma, \phi, \sigma \sim \text{Norm}(l(d_i;\gamma, \phi), \sigma) \;\forall i. \end{align}

We will choose weakly informative priors. This is a departure from our approach the previous time we studied this data set in which we used Jeffreys priors. The Jeffreys priors allowed for an analytical treatment, though they did not codify our prior knowledge. If we wish to properly encode our prior knowledge in the prior, we usually have to sacrifice analytical tractability. But that’s ok, since we have MCMC!

We will build the prior for the three parameters \(\phi\), \(\gamma\), and \(\sigma\), taking each to be independent of the others. Let’s start with \(\sigma\). It is possible that the spindle size is very carefully controlled. It is also possible that it could be highly variable. So, we can choose a Half Normal prior for \(\sigma\) with a large scale parameter of 10 µm; \(\sigma \sim \text{HalfNorm}(10)\). For \(\gamma\), we know the physically is must be between zero and one (as we have worked out in previous lessons), so we will assume a prior on is close to uniform that domain, but with probability density dropping right at zero and one; \(\gamma \sim \text{Beta}(1.1, 1.1)\). For \(\phi\), the typical largest spindle size, we assume a smallest spindle of about a micron (right where I would start to feel uncomfortable betting a year’s salary) and the largest spindle to be about a millimeter (the size of a very large cell, like a Xenopus egg). We therefore take \(\log_{10}\phi \sim \text{Norm}(1.5, 0.75)\), where \(\phi\) is in units of microns.

We thus have our complete generative model.

\begin{align} &\log_{10} \phi \sim \text{Norm}(1.5, 0.75),\\[1em] &\gamma \sim \text{Beta}(1.1, 1.1), \\[1em] &\sigma \sim \text{HalfNorm}(10),\\[1em] &\mu_i = \frac{\gamma d_i}{\left(1+(\gamma d_i/\phi)^3\right)^{\frac{1}{3}}}, \\[1em] &l_i \mid d_i, \gamma, \phi, \sigma \sim \text{Norm}(\mu_i, \sigma) \;\forall i. \end{align}

Using Stan to sample

The Stan code we will use is below.

functions {
  real ell_theor(real d, real phi, real gamma_) {
    real denom_ratio = (gamma_ * d / phi)^3;
    return gamma_ * d / (1 + denom_ratio)^(1.0 / 3.0);
  }
}


data {
  int N;
  real d[N];
  real ell[N];
}


parameters {
  real log10_phi;
  real gamma_;
  real<lower=0> sigma;
}


transformed parameters {
  real phi = 10^log10_phi;
}


model {
  log10_phi ~ normal(1.5, 0.75);
  gamma_ ~ beta(1.1, 1.1);
  sigma ~ normal(0.0, 10.0);

  for (i in 1:N) {
    ell[i] ~ normal(ell_theor(d[i], phi, gamma_), sigma);
  }
}

I pause to point out a few syntactical elements.

  1. Note how I used gamma_ as the variable name for gamma. This is because gamma is used as a distribution name in Stan.

  2. I have used the functions block in this Stan program. The definition of the function is real ell_theor(real d, real phi, real gamma_) { }, where the code to be executed in the function is between the braces. This is how you declare a function in Stan. The first word, real specifies what data type is expected to be returned. Then comes the name of the function, followed by its arguments in parentheses preceded by their data type.

  3. Within the function, I had to declare denom_ratio to be real. Remember, Stan is statically typed, so every variable needs a declaration of its type.

  4. Also within the function, note that the raise-to-power operator is ^, not ** as in Python.

  5. The return statement is like Python; return followed by a space and then an expression for the value to be returned.

  6. There is no Half-Normal distribution in Stan. It is implemented by putting a bound on the variable (in this case real<lower=0> sigma;) and then model sigma as being Normally distributed with location parameter zero.

  7. I use the transformed parameters block to conveniently handle the fact that I wanted to give my prior for \(\phi\) in terms of its base-ten logarithm.

All right! Let’s do some sampling!

[3]:

sm = cmdstanpy.CmdStanModel(stan_file='spindle.stan')

data = dict(
    N=len(df),
    d=df["Droplet Diameter (um)"].values,
    ell=df["Spindle Length (um)"].values,
)

samples = sm.sample(
    data=data,
    chains=4,
    iter_sampling=1000
)

samples = az.from_cmdstanpy(samples)

INFO:cmdstanpy:compiling stan program, exe file: /Users/bois/Dropbox/git/bebi103_course/2021/b/content/lessons/11/spindle
INFO:cmdstanpy:compiler options: stanc_options=None, cpp_options=None
INFO:cmdstanpy:compiled model file: /Users/bois/Dropbox/git/bebi103_course/2021/b/content/lessons/11/spindle
INFO:cmdstanpy:start chain 1
INFO:cmdstanpy:start chain 2
INFO:cmdstanpy:start chain 3
INFO:cmdstanpy:start chain 4
INFO:cmdstanpy:finish chain 3
INFO:cmdstanpy:finish chain 2
INFO:cmdstanpy:finish chain 1
INFO:cmdstanpy:finish chain 4

With samples in hand, let’s make a plot of the posterior samples. Remember that when we plotted the posterior by brute force, we had to analytically marginalize out \(\sigma\) (which we could do because we used a Jeffreys prior). Here, we can just plot our samples of \(\phi\) and \(\gamma\) and ignore the \(\sigma\) samples. So, we can look at a corner plot and check out the \(\phi\)-\(\gamma\) plane.

[4]:

bokeh.io.show(
    bebi103.viz.corner(
        samples,
        parameters=[("phi", "ϕ"), ("gamma_", "γ"), ("sigma", "σ")],
        plot_ecdf=True,
    )
)

The ECDFs of the samples depicted in the corner plot enable quick reference to the marginalized credible intervals.

We now have a complete picture of the posterior, computed directly from our samples. We will eschew plotting regression lines because this is not a preferred method of visualization and defer that to when we do posterior predictive checks.

[5]:

bebi103.stan.clean_cmdstan()

Computing environment

[6]:

%load_ext watermark
%watermark -v -p numpy,pandas,cmdstanpy,arviz,bokeh,iqplot,holoviews,bebi103,jupyterlab
print("cmdstan   :", bebi103.stan.cmdstan_version())

Python implementation: CPython
Python version       : 3.8.5
IPython version      : 7.19.0

numpy     : 1.19.2
pandas    : 1.2.0
cmdstanpy : 0.9.67
arviz     : 0.11.0
bokeh     : 2.2.3
iqplot    : 0.2.0
holoviews : 1.14.0
bebi103   : 0.1.2
jupyterlab: 2.2.6

cmdstan   : 2.25.0