The Rats growth data — hierarchical linear model (from scratch)¶
BiRats (single-normal hierarchy) and the mixture extension¶
References: Gelfand, Hills, Racine-Poon & Smith (1990, JASA); BUGS Vol I "Rats"/"BiRats".
| Section | Content |
|---|---|
| Data | 30 rats weighed at 5 ages — linear growth |
| Model | Random intercept + slope (hierarchical linear); BiRats = full-covariance single normal |
| Section 1 | Fit BiRats (K=1): population growth, shrinkage, intercept–slope correlation |
| Section 2 | Mixture extension (K≥2): do growth curves cluster into types? |
From-scratch engine: hlm_mixture_gibbs.py (rhier_linear_mixture) — K=1 is the single-normal hierarchy, K≥2 the mixture. (The outcome is continuous weight, so this is a hierarchical linear model — hier_mnl_gibbs.py is for multinomial choice and does not apply.)
The fitted model¶
Data. 30 rats, each weighed at five ages $x_j \in \{8, 15, 22, 29, 36\}$ days. $y_{ij}$ = weight (grams) of rat $i$ at age $x_j$.
Likelihood — one linear growth curve per rat (age centred at the midpoint $\bar x = 22$): $$y_{ij} = \alpha_i + \beta_i\,(x_j - 22) + \varepsilon_{ij}, \qquad \varepsilon_{ij}\sim N(0,\sigma^2).$$ Centring makes the parameters interpretable: $\alpha_i$ is rat $i$'s weight at day 22 (the middle of the study) and $\beta_i$ is its growth rate in g/day. The response is raw weight in grams (not log-weight) — the canonical BUGS Rats/BiRats specification, which reproduces the textbook population values (≈ 242 g at day 22, ≈ 6.2 g/day).
Hierarchy — the population of growth curves. The per-rat (intercept, slope) pairs share a population distribution: $$(\alpha_i,\beta_i)' \sim N(\mu,\ \Sigma)\quad\text{(BiRats: a single bivariate normal, Section 1)},$$ generalised in Section 2 to a $K$-component mixture of normals $(\alpha_i,\beta_i)'\sim\sum_{k=1}^{K}\pi_k\,N(\mu_k,\Sigma_k)$ to test for distinct growth types (segments).
Priors (weakly informative, conjugate): $\mu\sim N(\bar\mu, A_\mu^{-1})$ with $A_\mu=0.01$; $\Sigma\sim \mathrm{IW}(\nu=4,\ V=\mathrm{diag}(400,\,0.5))$; $\sigma^2_i\sim \mathrm{IG}(\nu_e/2,\cdot)$ with $\nu_e=3$ (per-rat error variance); mixture weights $\pi\sim\mathrm{Dirichlet}(a=5)$.
The data¶
30 young rats, each weighed at ages 8, 15, 22, 29, 36 days (a 30×5 table of weights in grams). Growth over this window is close to linear, so each rat is a small linear regression of weight on age. It is genuine panel data (5 repeated measures per unit) — what a hierarchical model needs. We centre age at $\bar x = 22$ so the intercept is weight at mid-study and decouples from the slope.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from hlm_mixture_gibbs import rhier_linear_mixture
Y = pd.read_csv('rats.csv', index_col=0).values # (30, 5) weights
ages = np.array([8, 15, 22, 29, 36.]); xc = ages - 22 # centred age
print(f'{Y.shape[0]} rats x {Y.shape[1]} weighings; mean weight by age: {Y.mean(0).round(0)}')
print(f'mean growth: {((Y[:,-1]-Y[:,0]).mean()/(ages[-1]-ages[0])):.2f} g/day')
plt.figure(figsize=(7,4.4))
for i in range(Y.shape[0]):
plt.plot(ages, Y[i], '-', color='steelblue', alpha=.35, lw=1)
plt.plot(ages, Y.mean(0), 'k-o', lw=2.5, label='population mean')
plt.xlabel('age (days)'); plt.ylabel('weight (g)'); plt.title('Rat growth curves (30 rats)')
plt.legend(); plt.tight_layout(); plt.savefig('rats_data.png', dpi=120, bbox_inches='tight'); plt.show()
30 rats x 5 weighings; mean weight by age: [152. 202. 245. 290. 325.] mean growth: 6.17 g/day
Model — hierarchical linear (random intercept + slope)¶
$$y_{ij} \sim N\!\big(\alpha_i + \beta_i\,(x_j-\bar x),\ \sigma^2\big),\qquad (\alpha_i,\beta_i)' \sim \text{MvN}(\mu,\ \Sigma),\quad \Sigma\sim\text{IW}(\nu,V).$$
Each rat has its own intercept $\alpha_i$ (weight at day 22) and slope $\beta_i$ (g/day); the population of $(\alpha_i,\beta_i)$ is a bivariate normal with mean $\mu$ and covariance $\Sigma$ — this single-normal version is BiRats (= rhier_linear_mixture with K=1). The mixture (K≥2) replaces that bivariate normal with a mixture, asking whether rats fall into discrete growth types.
regdata = [{'y': Y[i], 'X': np.column_stack([np.ones(5), xc])} for i in range(30)]
ols = np.array([np.linalg.lstsq(d['X'], d['y'], rcond=None)[0] for d in regdata]) # per-rat OLS
ssq = np.array([max(np.sum((d['y'] - d['X'] @ ols[i])**2)/3, 1e-2) for i, d in enumerate(regdata)])
Prior = dict(a=5., mubar=ols.mean(0), Amu=0.01, nu=4, V=np.diag([400., 0.5]), nu_e=3.0, ssq=ssq)
out = rhier_linear_mixture(regdata, K=1, Prior=Prior, Mcmc=dict(R=6000, burn=2000, keep=1, nprint=0), seed=1)
mu = out['mudraw'].mean(0)[0]; Sig = out['Sigmadraw'].mean(0)[0]; sig = np.sqrt(out['sigmasqdraw'].mean())
beta_b = out['betadraw'].mean(2) # (30,2) shrunken per-rat estimates
print(f'population intercept (weight at day 22) = {mu[0]:.1f} g')
print(f'population slope = {mu[1]:.2f} g/day')
print(f'SD across rats: alpha = {np.sqrt(Sig[0,0]):.1f} g, beta = {np.sqrt(Sig[1,1]):.2f} g/day')
print(f'intercept-slope correlation = {Sig[0,1]/np.sqrt(Sig[0,0]*Sig[1,1]):.2f}')
print(f'residual sigma = {sig:.2f} g')
print(f'slope shrinkage: OLS sd {ols[:,1].std():.2f} -> Bayes sd {beta_b[:,1].std():.2f}')
Done in 0.2 min | kept: 4000 population intercept (weight at day 22) = 242.3 g population slope = 6.21 g/day SD across rats: alpha = 14.0 g, beta = 0.51 g/day intercept-slope correlation = 0.55 residual sigma = 7.26 g slope shrinkage: OLS sd 0.57 -> Bayes sd 0.44
from matplotlib.patches import Ellipse
fig, ax = plt.subplots(1, 2, figsize=(13, 4.8))
# (a) shrinkage of per-rat slopes: OLS -> Bayesian (pulled toward population mean)
ax[0].plot([ols[:,1].min(), ols[:,1].max()], [ols[:,1].min(), ols[:,1].max()], 'k:', lw=1)
ax[0].scatter(ols[:,1], beta_b[:,1], color='firebrick', zorder=3)
ax[0].axhline(mu[1], color='gray', lw=.8, ls='--', label=f'pop. slope {mu[1]:.2f}')
ax[0].set_xlabel('per-rat OLS slope'); ax[0].set_ylabel('Bayesian (shrunken) slope')
ax[0].set_title('Shrinkage: noisy rats pulled toward the mean'); ax[0].legend()
# (b) (alpha_i, beta_i) with the fitted bivariate-normal 95% ellipse -> one elongated cluster
ax[1].scatter(beta_b[:,0], beta_b[:,1], color='steelblue', zorder=3)
vals, vecs = np.linalg.eigh(Sig)
o = vals.argsort()[::-1]; vals, vecs = vals[o], vecs[:, o] # descending: major axis first
width, height = 2 * np.sqrt(vals * 5.991) # 95% ellipse (major, minor)
ang = np.degrees(np.arctan2(vecs[1, 0], vecs[0, 0])) # angle of the major axis
ax[1].add_patch(Ellipse(mu, width, height, angle=ang, fill=False, edgecolor='firebrick', lw=2,
label='fitted 95% ellipse'))
ax[1].set_xlabel(r'intercept $\alpha_i$ (g at day 22)'); ax[1].set_ylabel(r'slope $\beta_i$ (g/day)')
ax[1].set_ylim(mu[1] - 4*np.sqrt(Sig[1,1]), mu[1] + 4*np.sqrt(Sig[1,1]))
ax[1].set_title('Per-rat growth params + fitted bivariate normal'); ax[1].legend()
plt.tight_layout(); plt.savefig('rats_birats.png', dpi=120, bbox_inches='tight'); plt.show()
Left — shrinkage. Each rat's noisy per-rat OLS slope (x-axis) against its Bayesian shrunken slope (y-axis), with the population slope (6.21 g/day) dashed. Points are pulled off the 45° line toward that mean, and the cross-rat SD of the slope tightens from 0.57 → 0.44 g/day. With only 5 weigh-ins per rat, an individual growth rate is barely identified on its own — so partial pooling does much of the work.
Right — the population of growth curves. The 30 rats in (intercept, slope) space with the fitted bivariate-normal 95% ellipse. Note the very different scales: the ellipse is wide in intercept (SD ≈ 14 g) but narrow in slope (SD ≈ 0.5 g/day), and it is tilted — an intercept–slope correlation of 0.55, so heavier rats also grow a little faster. Crucially it is one elongated cloud, not two — foreshadowing Section 2's finding that there are no discrete growth types.
# Fitted per-rat growth curves (shrunken), coloured by growth rate
from matplotlib import cm
from matplotlib.colors import Normalize
fig, ax = plt.subplots(figsize=(7.6, 4.7))
xg = np.linspace(ages.min(), ages.max(), 50)
slopes = beta_b[:, 1]; nrm = Normalize(slopes.min(), slopes.max())
for i in range(30):
ax.plot(ages, Y[i], '.', color='lightgray', ms=4, zorder=1)
ax.plot(xg, beta_b[i, 0] + beta_b[i, 1]*(xg - 22), '-', color=cm.viridis(nrm(slopes[i])), lw=1, alpha=.85, zorder=2)
ax.plot(xg, mu[0] + mu[1]*(xg - 22), 'k-', lw=2.6, label='population fit')
sm = cm.ScalarMappable(cmap='viridis', norm=nrm); sm.set_array([])
plt.colorbar(sm, ax=ax, label='rat growth rate (g/day)')
ax.set_xlabel('age (days)'); ax.set_ylabel('weight (g)')
ax.set_title('Fitted per-rat growth curves (shrunken estimates) over the data')
ax.legend(loc='upper left')
plt.tight_layout(); plt.savefig('rats_fitted.png', dpi=120, bbox_inches='tight'); plt.show()
1. Results — BiRats (single-normal hierarchy)¶
The from-scratch sampler reproduces the canonical numbers: population intercept ≈ 242 g (weight at day 22), slope ≈ 6.2 g/day, residual σ ≈ 7 g, with SD across rats ≈ 14 g (intercept) and ≈ 0.5 g/day (slope) and a positive intercept–slope correlation ≈ 0.5 (bigger rats grow a little faster). The shrinkage panel shows per-rat slopes pulled in from their noisy OLS values toward the population mean — the hallmark of partial pooling, and exactly why the hierarchy helps when each unit has only 5 observations. (Note the very different scales: the fitted ellipse is wide in intercept, ~±14 g, but narrow in slope, ~±0.5 g/day.)
2. Mixture extension: are there growth types?¶
Replace the single bivariate normal with a K-component mixture and compare by WAIC.
XtX = np.stack([d['X'].T @ d['X'] for d in regdata]); Xty = np.stack([d['X'].T @ d['y'] for d in regdata])
yty = np.array([d['y'] @ d['y'] for d in regdata]); Tn = np.full(30, 5.)
def waic(o, ss=2):
mud, Sigd, pid, s2d = o['mudraw'], o['Sigmadraw'], o['pidraw'], o['sigmasqdraw']
Rk, K = mud.shape[0], mud.shape[1]; idx = np.arange(0, Rk, ss); LL = np.empty((len(idx), 30))
for ri, r in enumerate(idx):
s2 = s2d[r]; llk = np.empty((30, K))
for k in range(K):
Sk = Sigd[r,k]; Si = np.linalg.inv(Sk); mu = np.tile(mud[r,k], (30,1))
Xte = Xty - np.einsum('ipq,iq->ip', XtX, mu)
ete = yty - 2*np.sum(mu*Xty,1) + np.einsum('ip,ipq,iq->i', mu, XtX, mu)
M = np.linalg.inv(Si[None] + XtX/s2[:,None,None]); quad = ete/s2 - np.einsum('ip,ipq,iq->i', Xte, M, Xte)/s2**2
B = np.eye(2)[None] + np.einsum('pq,iqr->ipr', Sk, XtX)/s2[:,None,None]
llk[:,k] = np.log(pid[r,k]+1e-300) - 0.5*(Tn*np.log(2*np.pi) + Tn*np.log(s2)+np.linalg.slogdet(B)[1] + quad)
m = llk.max(1); LL[ri] = m + np.log(np.exp(llk-m[:,None]).sum(1))
mm = LL.max(0); return -2*((np.log(np.exp(LL-mm).mean(0))+mm).sum() - LL.var(0).sum())
ks = [1, 2, 3]; waics = []
print(f"{'K':>2}{'WAIC':>10} component means [intercept, slope]")
for K in ks:
o = rhier_linear_mixture(regdata, K=K, Prior=Prior, Mcmc=dict(R=6000, burn=2000, keep=1, nprint=0), seed=1)
w = waic(o); waics.append(w)
print(f"{K:>2}{w:>10.1f} {np.round(o['mudraw'].mean(0),1).tolist()} pi={o['pidraw'].mean(0).round(2)}")
best = ks[int(np.argmin(waics))]
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(ks, waics, 'o-', color='steelblue', lw=2)
ax.scatter([best], [min(waics)], color='firebrick', s=90, zorder=5, label=f'best: K={best}')
ax.set_xticks(ks); ax.set_xlabel('number of mixture components K'); ax.set_ylabel('WAIC (lower = better)')
ax.set_title('Mixture model selection: K=1 wins → no distinct growth segments'); ax.legend()
plt.tight_layout(); plt.savefig('rats_waic.png', dpi=120, bbox_inches='tight'); plt.show()
K WAIC component means [intercept, slope]
Done in 0.2 min | kept: 4000 1 1089.1 [[242.3, 6.2]] pi=[1.]
Done in 0.2 min | kept: 4000 2 1090.6 [[238.4, 6.1], [247.1, 6.4]] pi=[0.53 0.47]
Done in 0.3 min | kept: 4000
3 1091.0 [[235.7, 6.0], [242.4, 6.2], [249.6, 6.5]] pi=[0.34 0.34 0.31]
2. Results — no growth segments¶
| K | WAIC | components |
|---|---|---|
| 1 | 1089.1 | [242, 6.2] |
| 2 | 1090.6 | [238, 6.1], [247, 6.4] |
| 3 | 1091.0 | three near-identical |
K=1 is best; the mixture does not improve fit and its components are nearly identical (all ≈ 242 g / 6.2 g/day). With 30 rats the population of growth curves is a single bivariate normal — there are no discrete growth types. This is the cleanest instance of the recurring lesson: the hierarchical normal earns its keep (shrinkage, partial pooling), while a mixture-as-segments adds nothing here. (To make the mixture machinery visibly bite you'd need genuinely clustered units; the bivariate normal is the right model for Rats.)
Takeaway. Rats is the canonical hierarchical-linear backbone — rhier_linear_mixture at K=1 recovers Gelfand et al.'s growth model (α≈242, β≈6.2, shrinkage), and the mixture extension confirms a single normal suffices.