5. Meta-analysis via Hierarchical Models*

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

rng = np.random.default_rng(6420)

5. Meta-analysis via Hierarchical Models*#

Adapted from Unit 7: rats_nocentering.odc.

This example is taken from Gelfand et al. (1990) and concerns 30 young rats whose weights were measured weekly for five weeks. Measurements were taken at 8, 15, 22, 29, and 36 days.

The observed weights (in grams) for each rat are arranged in the matrix \(\mathbf{Y}\), where \(y_{i,j}\) is the weight of the \(i\)-th rat at age \(x_j\). The first two and last rows of \(\mathbf{Y}\) are shown below:

Rat

\(y_{i,8} \)

\( y_{i,15} \)

\( y_{i,22} \)

\( y_{i,29} \)

\( y_{i,36} \)

1

151

199

246

283

320

2

145

199

249

293

354

…

…

…

…

…

…

30

153

200

244

286

324

Each row corresponds to a different rat (\( i = 1, \ldots, 30 \)), and each column corresponds to one of the five measurement times.

Data for \(x\) can be found here, and \(\mathbf{Y}\) here.

A hierarchical model lets us model each rat’s growth curve as a linear function of age, while also allowing the intercept (baseline weight) and slope (growth rate) to vary across rats. This structure is essentially “borrowing strength” across units in a meta-analysis. It allows us to estimate both group-level (population) and individual-level effects.

Note

The term meta-analysis was carried over from lectures here, but it’s kind of an odd choice of words. Typically meta-analysis refers to combining multiple studies.

The model is defined as:

\[\begin{split} \begin{align} y_{ij} &\sim \mathcal{N}(\mu_{ij}, \sigma^2)\\ \mu_{ij} &= \alpha_i + \beta_i x_j\\ \alpha_i &\sim \mathcal{N}(\mu_{\alpha}, \tau_{\alpha}^{-1}) \\ \beta_i &\sim \mathcal{N}(\mu_{\beta}, \tau_{\beta}^{-1})\\ \mu_{\alpha} &\sim \mathcal{N}(0, 10^6) \\ \tau_{\alpha} &\sim \text{Gamma}(0.001, 0.001) \\ \mu_{\beta} &\sim \mathcal{N}(0, 10^6) \\ \tau_{\beta} &\sim \text{Gamma}(0.001, 0.001) \\ \tau &\sim \text{Gamma}(0.001, 0.001) \\ \sigma &= \tau^{-1/2} \end{align} \end{split}\]

  • Line 1-2 are the likelihood

  • Lines 3 and 4 are the random effects

  • Lines 7-9 are hyperpriors

We will use a graphviz below to verify that we’ve specified our model correctly. This is particularly useful in models with complex hierarchical structures.

y = np.loadtxt("../data/rats_nocenteringY.txt")
x = np.loadtxt("../data/rats_nocenteringx.txt")
x.shape, y.shape

((5,), (30, 5))

with pm.Model() as m:
    mu_alpha = pm.Normal("mu_alpha", 0, tau=1e-6)
    tau_alpha = pm.Gamma("tau_alpha", 0.001, 0.001)
    mu_beta = pm.Normal("mu_beta", 0, tau=1e-6)
    tau_beta = pm.Gamma("tau_beta", 0.001, 0.001)

    tau = pm.Gamma("tau", 0.001, 0.001)
    sigma = pm.Deterministic("sigma", 1 / tau**0.5)

    alpha = pm.Normal("alpha", mu_alpha, tau=tau_alpha, shape=(30, 1))
    beta = pm.Normal("beta", mu_beta, tau=tau_beta, shape=(30, 1))

    mu = alpha + beta * x
    pm.Normal("y", mu, tau=tau, observed=y)

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

Hide code cell output

 
NUTS[nutpie]: [mu_alpha, tau_alpha, mu_beta, tau_beta, tau, alpha, beta]













az.summary(trace, var_names=["~alpha", "~beta"])










  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      mu_alpha
      106.58
      2.3
      100
      110
      12567
      9452
      1.00
      0.021
      0.015
    

    

      mu_beta
      6.186
      0.107
      6
      6.4
      13138
      9636
      1.00
      0.00093
      0.00069
    

    

      tau_alpha
      0.01017
      0.004
      0.0053
      0.017
      9026
      8851
      1.00
      4.5e-05
      6.8e-05
    

    

      tau_beta
      4.34
      1.5
      2.4
      6.9
      11327
      10206
      1.00
      0.014
      0.016
    

    

      tau
      0.02689
      0.00403
      0.021
      0.034
      7971
      8157
      1.00
      4.5e-05
      3.3e-05
    

    

      sigma
      6.15
      0.464
      5.5
      6.9
      7971
      8157
      1.00
      0.0052
      0.0038
    

  





If you have not had any courses in Regression or Design and Analysis of Experiments, take a moment and research how to interpret these coefficients: especially mu’s  and tau’s.






pm.model_to_graphviz(model=m)








../_images/4cc09b456cc2dbb32fdb6dee09e5d8a1c48d1484663902ee985aa3eb59695b18.svg









%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
pymc : 6.3.1