13. Multiple Regression: Brozek Index Prediction*

import pymc as pm
import pandas as pd
import numpy as np
import arviz as az

rng = np.random.default_rng(6420)

13. Multiple Regression: Brozek Index Prediction*#

Adapted from fat2d.odc and fatmulti.odc (they appear to be identical). This is the same dataset as the first regression example, available for download here.

This is not necessarily the best model for multiple regression; rather, it is just an attempt at implementing the lecture code, which is using Zellner’s G-prior Zellner [1986]. For multivariate regression in PyMC, we recommend pm.LKJCholeskyCov or pm.LKJCorr as your covariance matrix prior.

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

y = data["y"].to_numpy(copy=True)
X = data.iloc[:, 1:].to_numpy(copy=True)

# add intercept
X_aug = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1)
n, p = X_aug.shape

# Zellner's g
g = p**2

n, p, g

(252, 15, 225)

mu_beta = np.zeros(p)

with pm.Model() as m2d:
    tau = pm.Gamma("tau", 0.01, 0.01)
    variance = pm.Deterministic("variance", 1 / tau)

    tau_beta = (tau / g) * pm.math.dot(X_aug.T, X_aug)
    beta = pm.MvNormal("beta", mu_beta, tau=tau_beta)

    mu = pm.Deterministic("mu", pm.math.dot(X_aug, beta))
    pm.Normal("likelihood", mu=mu, tau=tau, observed=y)

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

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

Hide code cell output

 
NUTS[nutpie]: [tau, 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
    

  
  
    

      beta[0]
      -15.4
      16.6
      -48
      16
      3866
      4306
      1.00
      0.27
      0.19
    

    

      beta[1]
      0.0565
      0.0307
      -0.0033
      0.12
      10531
      9963
      1.00
      0.0003
      0.00021
    

    

      beta[2]
      -0.082
      0.051
      -0.18
      0.018
      4258
      5063
      1.00
      0.00079
      0.00057
    

    

      beta[3]
      -0.052
      0.106
      -0.26
      0.16
      7569
      7889
      1.00
      0.0012
      0.00088
    

    

      beta[4]
      0.067
      0.287
      -0.49
      0.63
      8294
      8198
      1.00
      0.0032
      0.0022
    

    

      beta[5]
      -0.443
      0.225
      -0.88
      0.0092
      8670
      8225
      1.00
      0.0024
      0.0017
    

    

      beta[6]
      -0.031
      0.099
      -0.22
      0.16
      6603
      7825
      1.00
      0.0012
      0.00086
    

    

      beta[7]
      0.874
      0.088
      0.7
      1
      9303
      8837
      1.00
      0.00092
      0.00065
    

    

      beta[8]
      -0.2
      0.142
      -0.48
      0.078
      5190
      6024
      1.00
      0.002
      0.0014
    

    

      beta[9]
      0.224
      0.141
      -0.049
      0.5
      9157
      8579
      1.00
      0.0015
      0.0011
    

    

      beta[10]
      0.003
      0.238
      -0.47
      0.47
      8097
      7650
      1.00
      0.0026
      0.0019
    

    

      beta[11]
      0.157
      0.211
      -0.27
      0.57
      12965
      9258
      1.00
      0.0019
      0.0013
    

    

      beta[12]
      0.148
      0.164
      -0.17
      0.46
      13114
      9677
      1.00
      0.0014
      0.001
    

    

      beta[13]
      0.43
      0.188
      0.062
      0.8
      12877
      9374
      1.00
      0.0017
      0.0012
    

    

      beta[14]
      -1.48
      0.509
      -2.5
      -0.48
      8408
      8224
      1.00
      0.0056
      0.004
    

    

      tau
      0.0595
      0.0053
      0.049
      0.07
      10988
      7533
      1.00
      5e-05
      3.6e-05
    

    

      variance
      16.94
      1.52
      14
      20
      10988
      7533
      1.00
      0.015
      0.011
    

    

      br2
      0.7337
      0.0239
      0.68
      0.78
      10988
      7533
      1.00
      0.00023
      0.00018
    

    

      br2_adj
      0.718
      0.0254
      0.66
      0.76
      10988
      7533
      1.00
      0.00024
      0.00019
    

  









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









bayesian_R2(mean=0.73, hdi_lb=0.68, hdi_ub=0.77)













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









Last updated: Sun, 18 Oct 2026

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

pytensor: 3.3.0

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