PyMC linear regression part 4: predicting actual height

At last, we have come to the end. This is the final post in a series of linear regression posts using PyMC3, from my reading of Statistical Rethinking. Part 1 set up the problem and explored the package’s objects. Part 2 interpreted the posterior distribution. Part 3 predicted average height which has its own uncertainty. Here, we’ll predict actual height which avoids the over-confidence that comes from predicting average height only.

Going deep into each step took some time, but it’s helped embed the concepts in my brain. As we’ve done before, we’ll start off duplicating code of the previous posts before we make the final predictions.

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import pymc3 as pm
import scipy.stats as stats
import seaborn as sns

%load_ext nb_black
%config InlineBackend.figure_format = 'retina'
%load_ext watermark
RANDOM_SEED = 8927
np.random.seed(RANDOM_SEED)
az.style.use("arviz-darkgrid")

Question

Ths question come’s from the winter 2020, week 2 homework.

The weights listed below were recorded in the !Kung census, but heights were not recorded for these individuals. Provide predicted heights and 89% compatibility intervals for each of these individuals. That is, fill in the table below, using model-based predictions.

Individual	weight	expected height	89% interval
1	45
2	40
3	65
4	31

Let’s quickly take a look at the data to get a handle on what we’re working with.

d = pd.read_csv("../data/a_input/Howell1.csv", sep=";", header=0)
d2 = d[d.age >= 18]  # filter to get only adults

Setting up the variables

# Get the average weight as part of the model definition
xbar = d2.weight.mean()

with pm.Model() as heights_model:

    # Priors are variables a, b, sigma
    # using pm.Normal is a way to represent the stochastic relationship the left has to right side of equation
    a = pm.Normal("a", mu=178, sd=20)
    b = pm.Lognormal("b", mu=0, sd=1)
    sigma = pm.Uniform("sigma", 0, 50)

    # This is a linear model (not really a prior or likelihood?)
    # Data included here (d2.weight, which is observed)
    # Mu is deterministic, but a and b are stochastic
    mu = a + b * (d2.weight - xbar)

    # Likelihood is height variable, which is also observed (data included here, d2.height))
    # Height is dependent on deterministic and stochastic variables
    height = pm.Normal("height", mu=mu, sd=sigma, observed=d2.height)

    # The next lines is doing the fitting and sampling all at once.
    trace_m2 = pm.sample(1000, tune=1000, return_inferencedata=False, progressbar=False)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, b, a]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 12 seconds.

trace_m2_df = pm.trace_to_dataframe(trace_m2)

The next set of code gives both the MAP as well as providing a sense of the variance of the data, which is more useful.

az.rcParams["stats.hdi_prob"] = 0.89  # sets default credible interval used by arviz
az.summary(trace_m2, round_to=3, kind="stats")

/Users/blacar/opt/anaconda3/envs/stats_rethinking/lib/python3.8/site-packages/arviz/data/io_pymc3.py:88: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(

	mean	sd	hdi_5.5%	hdi_94.5%
a	154.600	0.275	154.178	155.055
b	0.903	0.043	0.835	0.970
sigma	5.103	0.195	4.789	5.404

# Covariance matrix
trace_m2_df.cov().round(3)

	a	b	sigma
a	0.075	0.000	0.000
b	0.000	0.002	-0.000
sigma	0.000	-0.000	0.038

Predictions for actual height

We return now to the original problem where we are seeking predictions for actual height. We will produce compatability intervals where we’d expect to see any observed height (note how the uncertainty region for $\mu_i$ as shown above is much narrower than where the points are.)

We will do something similar as in the previous post but incorporating $\sigma$.

$\text{height}_i$ ~ Normal($\mu_i, \sigma$)

The repo initially confused me, but I saw that the solution was straightforward after more exploring of functions. I’ll spare you most of this exploration, but will put some of my key insights here.

Function exploration
One key realization for me was seeing how the stats function for drawing from a normal distribution can have varied mean (loc) and standard deviation (scale) values. Previously, I had typically used it to get multiple samples from the same mean and standard deviation like this:

stats.norm.rvs(loc=10, scale=2, size=10)

array([ 7.09403309, 12.255019  ,  7.32435791, 10.26594505, 11.79233188,
       12.09804331, 10.84097672, 14.51428249, 14.22362859,  9.25777341])

I hadn’t previously thought to use different mean and standard deviations for each draw. That is, loc and scale can be vectors. I’ll vary the vectors down below for each so that you can see what I mean in the resulting vector.

stats.norm.rvs(
    loc=[1, 1, 1, 10, 10, 10, 20, 20, 20], scale=[0.5, 1, 2, 0.5, 1, 2, 0.5, 1, 2]
)

array([ 0.73038967,  0.37624878,  1.87080349,  9.87713139, 10.7586794 ,
        9.36231624, 20.29745599, 19.36382859, 18.66087831])

This is when I had my a-ha moment about what to do next. Remember: the posterior distribution that pymc outputs produces a different mean and standard deviation in each line.

trace_m2_df.head()

	a	b	sigma
0	154.030585	0.853788	4.989442
1	153.969973	0.861619	5.128790
2	155.111070	0.957769	5.067568
3	154.106523	0.838728	4.961068
4	154.328210	0.863311	5.160797

We don’t need to use a multivariate normal distribution function but we can instead use the regular (non-multivariate?) version of the normal distribution function. This implementation is more consistent with R code 4.63 which I copied/posted below taken from McElreath’s repo.

## R code 4.63
post <- extract.samples(m4.3)
weight.seq <- 25:70
sim.height <- sapply( weight.seq , function(weight)
    rnorm(
        n=nrow(post) ,
        mean=post$a + post$b*( weight - xbar ) ,
        sd=post$sigma ) )
height.PI <- apply( sim.height , 2 , PI , prob=0.89 )

Making predictions for actual height ($h_i$) at a single weight value

When we input a weight, we can create a vector for height that results from the normal distribution that consists of the posterior’s parameters. This seems confusing at first, but let’s just look at one weight, 45.

Note how the parameters within the stats.norm.rvs function mimic McElreath’s R code above.

example_weight = 45
h_at_45 = stats.norm.rvs(
    loc=trace_m2_df["a"] + trace_m2_df["b"] * (example_weight - xbar),
    scale=trace_m2_df["sigma"],
)
print(h_at_45)
print("Length of vector: ", len(h_at_45))

[155.7948111  146.56221928 159.43749538 ... 155.33600385 149.83412841
 151.31287874]
Length of vector:  4000

We have a vector of values for which we can compute summary statistics, including getting a compatibility interval.

# Get 89% compatibility interval
az.hdi(np.array(h_at_45))

array([146.08867804, 162.48740096])

This seems repetitive to what we did for average height in the last post, but here we’ve incorporated sigma from our posterior distribution. Let’s repeat the code for average height and compare what we just did for actual height.

# Vector for average height - note that sigma is not used!
mu_at_45 = trace_m2_df["a"] + trace_m2_df["b"] * (45 - xbar)

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4), sharex=True)

az.plot_kde(mu_at_45, ax=ax1)
ax1.set_title("prediction of average height\nat weight of 45")
ax1.set_ylabel("density")
ax1.set_xlabel("average height ($\mu$)")
ax1.vlines(
    az.hdi(np.array(mu_at_45)),
    ymin=0,
    ymax=1.5,
    color="gray",
    linestyle="dashed",
    linewidth=1,
)

az.plot_kde(h_at_45, ax=ax2)
ax2.set_title("prediction of actual height\nat weight of 45")
ax2.set_ylabel("density")
ax2.set_xlabel("actual height ($h$)")
ax2.vlines(
    az.hdi(np.array(example_draw)),
    ymin=0,
    ymax=0.1,
    color="gray",
    linestyle="dashed",
    linewidth=1,
)

<matplotlib.collections.LineCollection at 0x7fe5ad6e2940>

I am intentionally plotting the x-axis with the same axis limits so that you can appreciate the difference in uncertainty. (Note the difference in the y-axis scale.) We would expect that the probability distribution for average height would be much tighter than the one for actual height. The plots drive that point home visually.

Making predictions for actual height ($h_i$) at a single weight value

We’ll start off by generating a matrix that contains the posterior predictions for actual heights over a range of weights, between 25 and 70 kg. This is similar to what we did for average heights in the last post.

# Input a range of weight values
weight_seq = np.arange(25, 71)

# This is generating a matrix where the predicted h values will be kept
# Each weight value will be its own row
h_pred = np.zeros((len(weight_seq), len(trace_m2_df)))

# Fill out the matrix in this loop
for i, w in enumerate(weight_seq):
    h_pred[i] = stats.norm.rvs(
        loc=trace_m2_df["a"] + trace_m2_df["b"] * (w - xbar), scale=trace_m2_df["sigma"]
    )

# Output the matrix in a dataframe for viewing
df_h_pred = pd.DataFrame(h_pred, index=weight_seq)
df_h_pred.index.name = "weight"
df_h_pred.head()

	0	1	2	3	4	5	6	7	8	9	...	3990	3991	3992	3993	3994	3995	3996	3997	3998	3999
weight
25	146.819087	131.950216	129.793433	141.425616	134.363193	134.193045	138.625503	127.385217	135.971128	132.278403	...	128.438620	128.187695	140.142773	136.868872	132.564979	140.662170	135.461373	129.179569	133.174905	139.440426
26	136.456664	139.880879	138.150449	144.660997	130.087085	133.466438	134.668864	130.054509	137.989162	144.520142	...	140.724558	122.746829	141.763728	140.855240	138.828098	125.438146	139.518013	145.103823	144.708165	137.155993
27	133.427578	146.873131	148.040623	134.562287	143.913001	141.570235	138.758449	126.249154	139.976170	136.232402	...	132.069314	142.616321	139.039138	144.439965	140.008615	132.783098	136.572733	131.534774	135.076601	145.095471
28	143.316194	130.229996	135.166804	140.417059	143.595620	146.727791	142.071108	135.161632	138.845736	140.713138	...	131.457780	131.502991	140.776164	144.107353	150.544606	144.113096	133.674960	137.908407	131.267668	141.673441
29	142.093570	142.195453	145.842916	144.098927	146.044805	138.547015	148.580323	134.374695	139.463730	141.527401	...	142.074329	145.405100	144.384839	131.004655	150.856133	139.352867	133.089377	139.672434	139.466480	133.085828

5 rows × 4000 columns

To re-iterate, each row are a weight’s sampling of heights as dictated by the posterior distribution of parameters. We can illustrate that graphically.

# Plotting the data as a scatter plot
plt.scatter(d2["weight"], d2["height"], alpha=0.5)
plt.plot(weight_seq, h_pred, "C0.", color="gray", alpha=0.005)
plt.xlabel("weight")
plt.ylabel("height")
plt.title("Uncertainty of $h_i$")

Text(0.5, 1.0, 'Uncertainty of $h_i$')

PyMC method for predictions of actual height

The above plot helps show the uncertainty for a range of weights after some step-by-step calculations. Not surprisingly, pymc and arviz can help us get to those steps more quickly. Here’s one way to do that.

# Predict heights
# set parameter to 200 to match the notebook example
# could set to 400 to match the samples and nchains
height_pred = pm.sample_posterior_predictive(trace_m2, 200, heights_model)
height_pred

/Users/blacar/opt/anaconda3/envs/stats_rethinking/lib/python3.8/site-packages/pymc3/sampling.py:1687: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(




<div>
    <style>
        /* Turns off some styling */
        progress {
            /* gets rid of default border in Firefox and Opera. */
            border: none;
            /* Needs to be in here for Safari polyfill so background images work as expected. */
            background-size: auto;
        }
        .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {
            background: #F44336;
        }
    </style>
  <progress value='200' class='' max='200' style='width:300px; height:20px; vertical-align: middle;'></progress>
  100.00% [200/200 00:00<00:00]
</div>






{'height': array([[152.44015719, 143.8660879 , 139.40217284, ..., 156.3123799 ,
         162.99350218, 160.15982149],
        [161.48691291, 152.72929135, 138.60784035, ..., 162.5252433 ,
         167.12620516, 160.8162708 ],
        [154.18712134, 148.66960604, 140.74062707, ..., 165.56777632,
         160.19434841, 157.95943753],
        ...,
        [160.39855129, 147.75720026, 140.01191018, ..., 159.5460152 ,
         153.48386574, 163.72311138],
        [149.698651  , 142.51656198, 145.54450618, ..., 167.74115779,
         165.92839171, 164.0738555 ],
        [159.73732783, 153.56864439, 141.05003146, ..., 166.80582294,
         161.04500287, 160.92227836]])}

f, ax1 = plt.subplots()
# ax = az.plot_hdi(weight_seq, mu_pred.T)
az.plot_hdi(d2["weight"], height_pred["height"], ax=ax1)
ax1.scatter(d2["weight"], d2["height"])
ax1.set_xlabel("weight")
ax1.set_ylabel("height")
ax1.set_title("Uncertainty of $h_i$\n(plotting with the help of arviz)")

/Users/blacar/opt/anaconda3/envs/stats_rethinking/lib/python3.8/site-packages/arviz/stats/stats.py:493: FutureWarning: hdi currently interprets 2d data as (draw, shape) but this will change in a future release to (chain, draw) for coherence with other functions
  warnings.warn(





Text(0.5, 1.0, 'Uncertainty of $h_i$\n(plotting with the help of arviz)')

The answer to the problem

After all that, we come to the answer, which can be written in a few lines of code.

# Initialize the table
df_answer = pd.DataFrame()
df_answer["individual"] = range(1, 5)
df_answer["weight"] = [45, 40, 65, 31]
df_answer["expected_height"] = None
df_answer["89perc_interval"] = None

# Fill out the table
for i, weight in enumerate(df_answer["weight"]):

    draw = stats.norm.rvs(
        loc=trace_m2_df["a"] + trace_m2_df["b"] * (weight - xbar),
        scale=trace_m2_df["sigma"],
    )

    df_answer.loc[i, "expected_height"] = "{0:0.2f}".format(draw.mean())
    # .at is for insertion of alist
    # https://stackoverflow.com/questions/26483254/python-pandas-insert-list-into-a-cell
    df_answer.at[i, "89perc_interval"] = "{0:0.2f}, {1:0.2f}".format(
        az.hdi(draw)[0], az.hdi(draw)[1]
    )

# Output the table
df_answer

	individual	weight	expected_height	89perc_interval
0	1	45	154.55	146.06, 162.41
1	2	40	150.22	141.98, 157.86
2	3	65	172.79	164.58, 180.99
3	4	31	141.96	133.63, 149.96

Summary: A lesson about two kinds of uncertainty

In a series of posts, I took a deep dive into this problem. I initially had not appreciated the two kinds of uncertainty that this problem illustrates. It took me some time to absorb the concept. This was in the McElreath book.

Rethinking: Two kinds of uncertainty. In the procedure above, we encountered both uncertainty in parameter values and uncertainty in a sampling process. These are distinct concepts, even though they are processed much the same way and end up blended together in the posterior predictive simulation. The posterior distribution is a ranking of the relative plausibilities of every possible combination of parameter values. The distribution of simulated outcomes, like height, is instead a distribution that includes sampling variation from some process that generates Gaussian random variables. This sampling variation is still a model assumption. It’s no more or less objective than the posterior distribution. Both kinds of uncertainty matter, at least sometimes. But it’s important to keep them straight, because they depend upon different model assumptions. Furthermore, it’s possible to view the Gaussian likelihood as a purely epistemological assumption (a device for estimating the mean and variance of a variable), rather than an ontological assumption about what future data will look like. In that case, it may not make complete sense to simulate outcomes.

Appendix: Environment and system parameters

%watermark -n -u -v -iv -w

Last updated: Wed May 19 2021

Python implementation: CPython
Python version       : 3.8.6
IPython version      : 7.20.0

numpy     : 1.20.1
scipy     : 1.6.0
arviz     : 0.11.1
seaborn   : 0.11.1
json      : 2.0.9
pandas    : 1.2.1
matplotlib: 3.3.4
pymc3     : 3.11.0

Watermark: 2.1.0