11. Brozek index prediction*

import pymc as pm
import pandas as pd
import numpy as np
import arviz as az
from pymc.math import dot
import pytensor.tensor as pt

rng = np.random.default_rng(6420)

11. Brozek index prediction*#

This example goes over linear regression and Bayesian \(R^2\).

Adapted from Unit 7: fat1.odc, fat2d.odc, and fatmulti.odc.

Data can be found here.

Percentage of body fat, age, weight, height, and ten body circumference measurements (e.g., abdomen) were recorded for 252 men. Percentage of body fat is estimated through an underwater weighing technique.

The data set has 252 observations and 15 variables. Brozek index (Brozek et al., 1963) is obtained by the underwater weighing while other 14 anthropometric variables are obtained using scales and a measuring tape.

  • y = Brozek index

  • X1 = 1 (intercept)

  • X2 = age

  • X3 = weight

  • X4 = height

  • X5 = adipose

  • X6 = neck

  • X7 = chest

  • X8 = abdomen

  • X9 = hip

  • X10 = thigh

  • X11 = knee

  • X12 = ankle

  • X13 = bicep

  • X14 = forearm

  • X15 = wrist

These anthropometric variables are less intrusive but also less reliable in assessing the body fat index.

Set a linear regression to predict the Brozek index from these body measurements.

Single predictor (X8)#

This is a recreation of fat1.odc.

data = pd.read_csv("../data/fat.tsv", sep="\t")

y = data["y"].to_numpy(copy=True)
X = data["X8"].to_numpy(copy=True)

# p will be the number of predictors + intercept (1 + 1 in this case)
n, p = X.shape[0], 2

with pm.Model() as m:
    tau = pm.Gamma("tau", 0.001, 0.001)
    beta0 = pm.Normal("beta0_intercept", 0, tau=0.001)
    beta1 = pm.Normal("beta1_abdomen", 0, tau=0.001)
    variance = pm.Deterministic("variance", 1 / tau)

    mu = pm.Deterministic("mu", beta0 + beta1 * X)
    pm.Normal("likelihood", mu=mu, tau=tau, observed=y)

    # Bayesian R2 from fat1.odc
    # no need to use this anymore, use az.bayesian_r2 instead
    sse = (n - p) * variance
    cy = y - y.mean()
    sst = dot(cy, cy)
    br2 = pm.Deterministic("br2", 1 - sse / sst)

    trace = pm.sample(2000, random_seed=rng)

Hide code cell output

 
NUTS[nutpie]: [tau, beta0_intercept, beta1_abdomen]













az.summary(trace, var_names="~mu", ci_prob=0.95)










  
    

      
      mean
      sd
      eti95_lb
      eti95_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      beta0_intercept
      -35
      2.49
      -40
      -30
      1266
      1258
      1.00
      0.07
      0.049
    

    

      beta1_abdomen
      0.583
      0.027
      0.53
      0.64
      1268
      1303
      1.00
      0.00075
      0.00052
    

    

      tau
      0.049
      0.0044
      0.041
      0.058
      3995
      3516
      1.00
      7e-05
      5.1e-05
    

    

      variance
      20.57
      1.88
      17
      25
      3995
      3516
      1.00
      0.03
      0.024
    

    

      br2
      0.659
      0.0312
      0.59
      0.71
      3995
      3516
      1.00
      0.0005
      0.00039
    

  









az.bayesian_r2(
    trace,
    pred_mean="mu",
    scale="variance",
    scale_kind="var",
    ci_prob=0.95,
    ci_kind="hdi",
)









bayesian_R2(mean=0.66, hdi_lb=0.6, hdi_ub=0.71)













All predictors#


This is not in the lecture, just adding an extra PyMC example with all predictors.






y = data["y"].to_numpy()
X = data.drop(columns="y")













with pm.Model(coords={"predictor": X.columns}) as m_multiple:
    tau = pm.Gamma("tau", 0.001, 0.001)
    alpha = pm.Normal("intercept", 0, tau=0.001)
    beta = pm.Normal("beta", 0, tau=0.001, dims="predictor")
    variance = pm.Deterministic("variance", 1 / tau)

    mu = pm.Deterministic("mu", alpha + pm.math.dot(X, beta))

    pm.Normal("likelihood", mu=mu, tau=tau, observed=y)

    trace = pm.sample(3000, random_seed=6420)









NUTS[nutpie]: [tau, intercept, beta]















az.summary(trace, var_names="~mu", ci_prob=0.95)










  
    

      
      mean
      sd
      eti95_lb
      eti95_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      intercept
      -12.2
      14.6
      -41
      16
      3700
      4990
      1.00
      0.24
      0.18
    

    

      beta[X2]
      0.0568
      0.0304
      -0.0021
      0.12
      10942
      10301
      1.00
      0.00029
      0.00021
    

    

      beta[X3]
      -0.073
      0.046
      -0.17
      0.016
      4108
      5302
      1.00
      0.00072
      0.00053
    

    

      beta[X4]
      -0.064
      0.102
      -0.26
      0.14
      7290
      8004
      1.00
      0.0012
      0.00086
    

    

      beta[X5]
      0.064
      0.28
      -0.48
      0.62
      7969
      8146
      1.00
      0.0031
      0.0022
    

    

      beta[X6]
      -0.456
      0.221
      -0.89
      -0.02
      7706
      8078
      1.00
      0.0025
      0.0019
    

    

      beta[X7]
      -0.039
      0.097
      -0.23
      0.16
      6869
      7952
      1.00
      0.0012
      0.00084
    

    

      beta[X8]
      0.879
      0.086
      0.71
      1
      8102
      8096
      1.00
      0.00095
      0.00068
    

    

      beta[X9]
      -0.217
      0.132
      -0.47
      0.048
      4763
      6044
      1.00
      0.0019
      0.0014
    

    

      beta[X10]
      0.224
      0.138
      -0.047
      0.5
      7669
      7633
      1.00
      0.0016
      0.0011
    

    

      beta[X11]
      -0.01
      0.228
      -0.45
      0.44
      8501
      8341
      1.00
      0.0025
      0.0018
    

    

      beta[X12]
      0.15
      0.209
      -0.26
      0.56
      13769
      9745
      1.00
      0.0018
      0.0012
    

    

      beta[X13]
      0.146
      0.163
      -0.17
      0.47
      13077
      9531
      1.00
      0.0014
      0.001
    

    

      beta[X14]
      0.428
      0.189
      0.05
      0.8
      13739
      9203
      1.00
      0.0016
      0.0012
    

    

      beta[X15  ]
      -1.49
      0.494
      -2.5
      -0.53
      7642
      7930
      1.00
      0.0056
      0.004
    

    

      tau
      0.0627
      0.00577
      0.052
      0.074
      12074
      7100
      1.00
      5.2e-05
      3.7e-05
    

    

      variance
      16.09
      1.5
      13
      19
      12074
      7100
      1.00
      0.014
      0.01
    

  









az.bayesian_r2(
    trace,
    pred_mean="mu",
    scale="variance",
    scale_kind="var",
    ci_prob=0.95,
    ci_kind="hdi",
)









bayesian_R2(mean=0.74, hdi_lb=0.69, hdi_ub=0.78)













%load_ext watermark
%watermark -n -u -v -iv









Last updated: Sun, 18 Oct 2026

Python implementation: CPython
Python version       : 3.14.7
IPython version      : 9.16.1

arviz   : 1.3.0
numpy   : 2.4.6
pandas  : 3.0.5
pymc    : 6.3.1
pytensor: 3.3.0