# PyMC linear regression part 4: predicting actual height

11 minute read

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 ofh_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


Categories:

Updated: