2
$\begingroup$

This question is based on question 1 of the week 2 Statistical Rethinking problems, i.e. q1 here: https://github.com/rmcelreath/stat_rethinking_2022/blob/main/homework/week02.pdf

I have a pandas data frame df with observed weights and heights (df.weight and df.height resp.)

I am trying to use linear regression (predicting weight from height) with Pymc3 to find an 89% confidence interval for the weight of an 140cm tall individual.

This is my setup so far

with pm.Model() as q1_model:
    # specify the model 
    alpha = pm.Normal('alpha', mu=45, sd=100)
    beta = pm.Normal('beta', mu=0, sd=10)
    sigma = pm.Uniform('sigma', lower=0, upper=50)
    weight = pm.Normal('weight', mu=alpha + beta * df.height_c, sd=sigma, observed=df.weight)
    
    # find the posterior distribution of the weight
    trace = pm.sample(1000)

where df.height_c is df.height - df.height.mean(). But I'm not sure how to get the interval. I tried:

pm.hdi(trace.alpha + (140 - df.height.mean()) * trace.beta, 0.89)

which gives [35.0, 36.6] but I think this is just an interval for the mean height mu of an 140cm person, not for height itself. The actual interval is supposedly meant to be [29.1, 42.8] confirming my suspicion here.

$\endgroup$
7
  • $\begingroup$ What is the function of 'sigma' in your model and is it involved in computing the posterior distribution for the weight of a person of 140cm? $\endgroup$ Jan 10, 2023 at 21:46
  • $\begingroup$ Sigma is a parameter of the normal distribution that the weight follows. I.e. if weight and height are $w$ and $h$, then $w$ follows distribution $N(\mu, \sigma^2)$, $\mu = \alpha + \beta h$. So sigma is just the standard deviation. $\endgroup$ Jan 10, 2023 at 22:18
  • $\begingroup$ Is this sigma used in the computation of the interval? Or do you only use the posterior of $\mu$? $\endgroup$ Jan 10, 2023 at 22:26
  • $\begingroup$ Regarding terminology, it is a credible interval you wish to calculate. $\endgroup$ Jan 10, 2023 at 23:43
  • $\begingroup$ @SextusEmpiricus Yes sigma is used in the computation of the interval $\endgroup$ Jan 11, 2023 at 9:35

1 Answer 1

1
$\begingroup$

This exercise asks for a model of weight given height and a 89% compatibility interval for the weight of an individual with height = 140cm.

Aside: @GrahamBornholt comments that "it is a credible interval you wish to calculate". So let's clarify the terminology. McElreath prefers the term "compatibility interval" to highlight that the "the interval indicates a range of parameter values compatible with the model and data." The crucial thing to keep in mind is not the name of the interval but its interpretation, which is fundamentally different from the meaning of a frequentist confidence interval. See What's the difference between a confidence interval and a credible interval?

I'll use a bit of math to go with the python code. Let $x$ denote height (in cm) and $y$ denote weight (in kg). Then the regression model for $y$ given $x$ is: $$ \begin{aligned} y = \alpha + \beta(x - \bar{x}) + e, \quad e \sim \operatorname{N}(0,\sigma^2) \end{aligned} $$

The weight of a 140cm-tall individual ($x = 140$) is modeled as the sum of two quantities: the average weight of 140cm-tall persons $\operatorname{E}(y) = \alpha + \beta(140 - \bar{x})$ plus the difference $e$ between the individual's weight and the average weight at that height; this difference is called the error.

So to construct a compatibility interval for the weight of a person (rather than a compatibility interval for the average weight) at a specific height, you have to simulate both $\operatorname{E}(y)$ and the error from the posterior. As McElreath explains in Section 4.4.3.5: sample from the posterior to simulate $y$s, not a distribution of plausible average $y$.

First set up and fit the model.

import numpy as np
import pandas as pd
import pymc3 as pm
import scipy.stats as stats

np.random.seed(1234)

howell = pd.read_csv("Howell1.csv", sep=";")
adults = howell.query("age >= 18")
x = adults["height"].values
y = adults["weight"].values
xbar = np.mean(x)

with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=45, sd=100)
    beta = pm.Normal("beta", mu=0, sigma=10)
    sigma = pm.Uniform("sigma", lower=0, upper=50)
    weight = pm.Normal("weight", mu=alpha + beta * (x - xbar), sd=sigma, observed=y)
    trace = pm.sample(1000, tune=1000)

The exercise asks about the weight of a 140cm-tall person.

x0 = 140

First we construct a 89% compatibility interval to represent the uncertainty about the average weight of 140cm-tall individuals.

pp_Ey_at_x0 = trace.alpha + trace.beta * (x0 - xbar)
pm.hdi(pp_Ey_at_x0, hdi_prod=0.89)
#> [34.8, 36.7]

Then we construct a 89% compatibility interval to represent the uncertainty about the weight of a person who is 140cm tall. [This is also called a prediction interval.]

pp_y_at_x0 = pp_Ey_at_x0 + stats.norm.rvs(0, trace.sigma)
pm.hdi(pp_y_at_x0, hdi_prod=0.89)
#> [27.8, 43.6]

PS. It would be worth it to revisit the prior specification. Does the prior weight ~ $\operatorname{N}(45,100)$ at the average adult height seem plausible? Careful consideration about the prior is a big part of Statistical Rethinking. And it's a big part of doing Bayesian analysis well. See Section 4.4.1.3 on choosing priors and especially Figure 4.5 for an illustration of prior predictive simulation.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.