R3. Choosing priors

Plot of Spindle Length vs. Droplet Diameter

[3]:

p = bokeh.plotting.figure(
    frame_width=300,
    frame_height=200,
    x_axis_label="droplet diamater (µm)",
    y_axis_label="spindle length (µm)",
)

p.circle(
    df_spindle["Droplet Diameter (um)"],
    df_spindle["Spindle Length (um)"],
    size=3,
    alpha=0.3,
)

bokeh.io.show(p)

ECDF of Spindle Length

[4]:

bokeh.io.show(iqplot.ecdf(df_spindle, q='Spindle Length (um)'))

Say the spindle length is independent of droplet size. We would then expect the spindle length to be given by

where $e_i$ is the noise component of the $i$the datum. So, we have a theoretical model for spindle length, $l=\varphi$, and to get a fully generative model, we need to model the errors $e_i$.

Let us use the following likelihood (for explanation on how this likelihood was chosen, please see the discussion of this data set from last term).

$$\begin{array}{r}{l}_i\sim \text{Norm}(\varphi ,\sigma )\mathrm{\forall }i.\end{array}$$

So, we will need to choose the priors $g(\varphi )$ and $g(\sigma )$.

Possibility #1:

Let’s think about what we know about mitotic spindles. I don’t know their exact range in size, but I would guess a reasonable choice for $\varphi$ might be around ~20 µm. Why? I know a eukaryotic cell is around 100 µm, so it’s going to be smaller than that. But I also know the mitotic spindle has to separate chromosomes during metaphase/anaphase, so they’re going to have some length probably greater than 1 µm. 20 um sounds reasonable, but I’m not super sure about it, so let’s be sure to keep it broad. I think the spindles could probably take on a large range of values. I also don’t think the mean length is going to be affected by any other parameters in my model. Based on this, I’ll make a rough sketch of my idea for $g(\varphi )$.

Great! I can use the Distribution Explorer to choose a reasonable prior based on my sketch. I will choose:

$$\begin{array}{r}g(\varphi )\sim \text{Norm}(20,20).\end{array}$$

We can plot the result.

[5]:

phi = np.linspace(-50, 100, 200)
p = bokeh.plotting.figure(
    frame_width=300, frame_height=200, x_axis_label="φ (µm)", y_axis_label="g(φ)"
)
p.line(phi, st.norm.pdf(phi, 20, 20), line_width=2)
bokeh.io.show(p)

For $g(\sigma )$, I might choose a normal prior again. I think spindle size may vary a fair amount. They are created spontaneously in a droplet, so I would guess there could be many sizes. Therefore, I might choose

[6]:

sigma = np.linspace(-10, 20, 100)
p = bokeh.plotting.figure(
    frame_width=300,
    frame_height=200,
    x_axis_label="sigma (µm)",
    y_axis_label="g(sigma)",
)
p.line(sigma, st.norm.pdf(sigma, 5, 5), line_width=2)
bokeh.io.show(p)

In this case, I chose the shape for the prior above because 1) I am speculating that the average length might vary by 25-70%, and 2) this standard deviation is likely to have a distribution that has the shape of a normal distribution.

Visualizing the Prior

The prior can be a bit confusing because it’s hard to think about what it means. When choosing priors, we think about our knowledge of the parameter before collecting data. From what we know, we can estimate the probability of the parameter taking on certain values, and therefore can choose a reasonable distribution for our prior probability. One question that came up in recitation a lot is: what does each part of the prior mean? For example, what does the µ in $g(\varphi )\sim \text{Normal}(\mu ,\sigma )$ represent? What about the $\sigma$?

Consider our prior $g(\varphi )\sim \text{Normal}(\mu ,\sigma )$. We want to choose a prior for $\varphi$ that represents the probability of $\varphi$ taking on possible values. We chose a normal prior, so we are concerned with prior parameters $\mu$ and $\sigma$. It is important to note that the $\mu$ in $g(\varphi )\sim \text{Normal}(\mu ,\sigma )$ is not the same parameter as the $\varphi$ in $l_i\sim \text{Norm}(\varphi ,\sigma )\mathrm{\forall }i$. In the likelihood, $l_i\sim \text{Norm}(\varphi ,\sigma )\mathrm{\forall }i$, $\varphi$ is the location parameter of the PDF of the likelihood, and therefore represents the center of the probability distribution of the values of $l_i$. In the prior for $\varphi$, $g(\varphi )\sim \text{Normal}(\mu ,\sigma )$, $\mu$ is the location parameter of the probability distribution of the possible values of $\varphi$, the center we were just talking about. Essentially, $g(\varphi )$ is giving the estimated probability distribution for where $\varphi$ lies.

What can be even more confusing is the $\sigma$ in $g(\varphi )\sim \text{Normal}(\mu ,\sigma )$. This is not related to the $\sigma$ in $l_i\sim \text{Norm}(\varphi ,\sigma )\mathrm{\forall }i$!! In $g(\varphi )\sim \text{Normal}(\mu ,\sigma )$, $\sigma$ is the scale parameter for our PDF for possible $\varphi$s. If we want to pick a broad normal prior probability distribution for $\varphi$, then $\sigma$ will be large.

Possibility #2:

You may be thinking: hold on, $\varphi$ and $\sigma$ could not physically have negative values, so why are we using normal priors that could potentially take on negative values? Totally fair point! We can resolve this in a few ways:

1. If the location parameter for our Normal prior is far away from 0, negative values have a very low probability anyway due to the light tails, so we won’t have to worry about them too much. We can just include a conditional statement in our model to avoid negative values for our prior.

2. We could choose a Log-Normal or Half-Normal prior to avoid this issue. Be sure to review those distributions on the distribution explorer to make sure they are reasonable choices.

Rethinking my prior $g(\varphi )$, I might consider a Log-Normal prior. 1) The Log-Normal prior will not take on negative values, which is physically true for $\varphi$, and 2) I think it is reasonable that $\text{ln}(\varphi )$ is normally distributed. I could expect the probability on the left side of $\mu$ to vanish more quickly than on the right. So, I will choose:

$$\begin{array}{r}g(\varphi )\sim \text{LogNorm}(\text{ln}20,0.75)\end{array}$$

Which looks like:

[7]:

phi = np.linspace(0, 100, 200)
p = bokeh.plotting.figure(
    frame_width=300, frame_height=200, x_axis_label="phi (µm)", y_axis_label="g(phi)"
)
p.line(phi, st.lognorm.pdf(phi, 0.75, loc=0, scale=20), line_width=2)
bokeh.io.show(p)

Looks reasonable!

Revisiting $g(\sigma )$, we know it doesn’t make sense for $\sigma$ to be negative. I think Half-Normal would be an appropriate choice here, especially considering our original plot of $g(\sigma )\sim \text{Norm}(5,5)$. It also doesn’t make sense that $\text{ln}(\sigma )$ would be normally distributed, so Half-Normal is a good choice. I might choose Half-Normal centered at 0 with a large scale parameter, such as

This would look like:

[8]:

sigma = np.linspace(0, 40, 200)
p = bokeh.plotting.figure(
    frame_width=300,
    frame_height=200,
    x_axis_label="sigma (µm)",
    y_axis_label="g(sigma)",
)
p.line(phi, st.halfnorm.pdf(sigma, 0, 10), line_width=2)
bokeh.io.show(p)

So, our final model would be:

Possibility #3: A New Model

Let’s use a new model for spindle length where droplet diameter is related to spindle length. To understand how we came up with this model, and its caveats, please refer to the discussion of the model from last term. Be sure to read through this first, otherwise the following will not make sense.

We have a theoretical model relating the droplet diameter to the spindle length. Let us now build a generative model. For spindle, droplet pair i, we assume

$$\begin{array}{r}{l}_i=\frac{\gamma d_i}{{(1+(\gamma d/\varphi {)}^3)}^{\frac{1}{3}}}+e_i.\end{array}$$

We will assume that $e_i$ is Normally distributed with variance ${\sigma }^2$. This leads us to our statistical model.

$$\begin{array}{rl}& {\mu }_i=\frac{\gamma d_i}{{(1+(\gamma d_i/\varphi {)}^3)}^{\frac{1}{3}}},\\ & l_i\sim \text{Norm}({\mu }_i,\sigma )\mathrm{\forall }i,\end{array}$$

where

$$\begin{array}{rl}\gamma & ={\left(\frac{T_0-T_{\min }}{k^2{T}_{\mathrm{s}}}\right)}^{\frac{1}{3}}\\ \varphi & =\gamma {\left(\frac{6\alpha \beta }{\pi }\right)}^{\frac{1}{3}},\end{array}$$

which is equivalently stated as

Yikes, we have a lot more parameters in this model! Not to worry, we will choose priors for one parameter at a time. Starting with $\varphi$, this still represents the location parameter for the distribution of spindle size. So, we will use the same prior:

$$\begin{array}{r}\varphi \sim \text{LogNorm}(\ln 20,0.75)\end{array}$$

Now let’s consider $\gamma$. $\gamma$ determines the relationship between $d$, the droplet diameter, and $\varphi$, the spindle length. We can already think of some bounds on $\gamma$. $\gamma$ must be greater than 0. Also, $\gamma$ must be less than or equal to 1, or else the spindle length exceeds the droplet diameter. Further, we can see that for large droplets, with $d\gg \varphi /\gamma$, the spindle size becomes independent of $d$, with

$$\begin{array}{r}l\approx \varphi .\end{array}$$

Conversely, for $ d \ll `:nbsphinx-math:phi`/:nbsphinx-math:`gamma`$, the spindle length varies approximately linearly with diameter.

$$\begin{array}{r}l(d)\approx \gamma d.\end{array}$$

Of course, these cases are physically possible. Please see the notes I referenced from BE/Bi 103a to review the limitations of this model. So, we know our model will only work in cases of droplets with intermediate diameter. Thus, we will choose:

$$\begin{array}{r}\gamma \sim \text{Beta}(2,2).\end{array}$$

Let’s take a look.

[9]:

gamma = np.linspace(0, 1, 100)
p = bokeh.plotting.figure(
    frame_width=300,
    frame_height=200,
    x_axis_label="γ (µm)",
    y_axis_label="γ(gamma)",
    x_range=[0, 1],
)
p.line(gamma, st.beta.pdf(gamma, 2, 2), line_width=2)
bokeh.io.show(p)

Notice how useful the beta distribution can be!

Next, we will consider $\sigma$. We could use our previous choice, $\sigma \sim \text{HalfNorm}(0,10)$. We could also consider the case where we assume $\sigma$ scales with $\varphi$. This is not an unreasonable assumption. As $\varphi$ increases, we may expect the standard deviation $\sigma$ to increase as well. We can write the expression:

$$\begin{array}{r}\sigma ={\sigma }_0\varphi \end{array}$$

We have now introduced a new prior parameter, ${\sigma }_0$. ${\sigma }_0$ is our constant of proportionality for the relationship between $\sigma$ and $\varphi$. We choose a gamma prior for ${\sigma }_0$:

$$\begin{array}{r}{\sigma }_0=\text{Gamma}(2,10)\end{array}$$

Gamma makes sense in this case because, as stated in the Distribution Explorer: “it imparts a heavier tail than the Half-Normal distribution (but not too heavy; it keeps parameters from growing too large), and allows the parameter value to come close to zero.” We want to leave room for the case where ${\sigma }_0$ is small, so a Gamma distribution is appropriate.

Let’s plot it up to take a look.

[10]:

sigma_0 = np.linspace(0, 2, 100)
p = bokeh.plotting.figure(
    frame_width=300, frame_height=200, x_axis_label="σ₀", y_axis_label="g(σ₀)"
)
p.line(sigma_0, st.gamma.pdf(sigma_0, 2, loc=0, scale=0.1), line_width=2)
bokeh.io.show(p)

We have now finished building the model:

$$\begin{array}{rl}& \varphi \sim \text{LogNorm}(\ln 20,0.75),\\ & \gamma \sim \text{Beta}(2,2),\\ & {\sigma }_0\sim \text{Gamma}(2,10),\\ & \sigma ={\sigma }_0\varphi ,\\ & {\mu }_i=\frac{\gamma d_i}{{(1+(\gamma d_i/\varphi {)}^3)}^{\frac{1}{3}}},\\ & l_i\sim \text{Norm}({\mu }_i,\sigma )\mathrm{\forall }i.\end{array}$$

ECDF of Catastrophe Times for Labeled Tubulin

[12]:

p = iqplot.ecdf(labeled, q='t (s)', conf_int=True)

bokeh.io.show(p)

Last term, we explored a Gamma model for this dataset, which gives the amount of time we have to wait for $\alpha$ arrivals of a Poisson process. This models a multistep process where each step happens at the same rate. So $\alpha$ is the number of arrivals, and $\beta$ is the rate of arrivals.

In the distribution explorer, it is easier to see the impacts of changing the parameters. When $\beta$ is small, the gamma distribution gets broad. This makes sense because the rate is small and thus the waiting time will likely be larger or more variable. When $\beta$ rate is large, the gamma distribution gets sharp; the waiting time will be short with a large rate.

We have two parameters - $\alpha$ and $\beta$. Let’s try to pick reasonable priors for $\alpha$ and $\beta$.

For $\beta$, we know that it is not physically possible for it to be negative or zero. We also know it will not be faster than the speed of light, ~ 3e8 m/s. In scipy, the lognormal pdf takes in a scale parameter, which is 1/$\beta$. It might be easier to think of a prior for 1/$\beta$, instead of $\beta$ itself. We can define $\theta$ = 1/$\beta$. If $\beta$ is very large, say approaching the speed of light, then $\theta$ will be very small, approaching zero.

I am not a biologist, but I would bet to say it is unlikely that microtubule catastrophe will take longer than one hours = 3600 seconds. 1/3600 ~ 2e-4 $\beta$ rate. So we could use 1/2e-4 = 5000 for a $\theta$ upper bound.

The range 0 – 5000 is a pretty wide, but we could use a lognormal, since I expect more probability towards faster rates than really slow ones. Let’s use a mean of 500.

[13]:

theta = np.linspace(0, 5000, 200)
p = bokeh.plotting.figure(
    frame_width=400, frame_height=300, x_axis_label="θ (s)", y_axis_label="g(θ)"
)
p.line(theta, st.lognorm.pdf(theta, 0.7, loc=0, scale=500), line_width=2)
bokeh.io.show(p)

For $\alpha$, we do have a rough bound. We know $\alpha$ must not be negative, or zero. We also know that if we use a mean $\theta$ of 500, then $\beta$ = 1/500 = .002. And, looking at the data, a mean catastrophe time is about 500. If we then use the mean moment of the gamma distribution, and set 500 = $\alpha$/ $\beta$ = $\alpha$/ .002, we get an order of magnitude estimate for $\alpha$ of 1. Based on this, it is not likely to be larger than 10.

[14]:

alpha = np.linspace(0, 20, 200)
p = bokeh.plotting.figure(
    frame_width=300, frame_height=200, x_axis_label="α", y_axis_label="g(α)"
)
p.line(alpha, st.halfnorm.pdf(alpha, 0, 5), line_width=2)
bokeh.io.show(p)