Hierarchical Linear Model — Python Gibbs sampler (hlm_gibbs.py)

Python companion to hlm_bayesm.ipynb (rhierLinearModel), matching it section-by-section.

Model. Unit $i$ has its own regression, and the unit coefficients are tied together by a population prior: $$y_i = X_i\,\beta_i + \varepsilon_i,\ \ \varepsilon_i\sim N(0,\sigma_i^2 I),\qquad \beta_i = \Delta'z_i + u_i,\ \ u_i\sim N(0,V_\beta).$$ The hierarchy lets each $\beta_i$ borrow strength (shrink toward the population surface $\Delta'z_i$) — the cure for noisy per-unit OLS.

Section	Content
Sections 2–3	Synthetic data — parameter recovery + trace plots
Section 4	Cheese scanner data — Exploratory data analysis + OLS baseline (why pool)
Section 5	Model A — hierarchy, no unit covariates
Section 6	Model B — large-market covariate $Z$
Section 7	OLS vs Model A vs Model B

1. Imports

import os, subprocess
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from hlm_gibbs import rhier_linear_model, simulate_hlm, posterior_summary
np.set_printoptions(suppress=True)

2. Synthetic data

Match the bayesm setup: $n_{\text{reg}}=60$ units, $n_{\text{obs}}=80$ each, $k=3$ regressors (intercept, $x_1$, $x_2$), $n_z=2$ unit covariates (constant + one continuous characteristic).

rng = np.random.default_rng(1)
nreg, nobs, nvar, nz = 60, 80, 3, 2
u  = rng.uniform(0, 5, nreg); zc = (u - u.mean())/u.std()
Z  = np.column_stack([np.ones(nreg), zc])                      # const + unit_char
Delta_true = np.array([[2.0, -1.5, 0.8],                       # const -> (intercept, x1, x2)
                       [0.5,  0.3, -0.2]])                     # unit_char effect
Vbeta_true = np.array([[0.5, 0.2, 0.0],
                       [0.2, 0.4, 0.1],
                       [0.0, 0.1, 0.3]])
tau_true = 0.5
regdata_s, Beta_true = simulate_hlm(nreg, nvar, Z, Delta_true, Vbeta_true,
                                    true_sigma2=tau_true, ni=nobs, seed=2)
print('Simulated %d units x %d obs; k=%d, nz=%d' % (nreg, nobs, nvar, nz))

Simulated 60 units x 80 obs; k=3, nz=2

3. rhier_linear_model — synthetic data

Prior_s = dict(Deltabar=np.zeros((nz, nvar)), A=0.01*np.eye(nz),
               nu=nvar+3, V=(nvar+2)*np.eye(nvar), nu_e=3.0, ssq=tau_true)
out_s = rhier_linear_model(regdata_s, Z, Prior=Prior_s,
                           Mcmc=dict(R=10000, burn=2000, nprint=0), seed=11)
Delta_pm = out_s['Deltadraw'].mean(0).reshape(nz, nvar, order='F')
Vbeta_pm = out_s['Vbetadraw'].mean(0).reshape(nvar, nvar, order='F')
beta_pm_s = out_s['betadraw'].mean(2)                          # (nreg, k)
print('Delta posterior mean:\n', np.round(Delta_pm, 3))
print('Delta truth:\n', Delta_true)
print('\nunit-beta recovery: corr(posterior mean, truth) = %.3f' %
      np.corrcoef(beta_pm_s.ravel(), Beta_true.ravel())[0,1])

Delta posterior mean:
 [[ 1.953 -1.513  0.813]
 [ 0.587  0.443 -0.222]]
Delta truth:
 [[ 2.  -1.5  0.8]
 [ 0.5  0.3 -0.2]]

unit-beta recovery: corr(posterior mean, truth) = 0.999

Recovery — trace plots and unit-beta scatter

fig, ax = plt.subplots(2, 2, figsize=(12, 6))
names_z = ['const','unit_char']; names_v = ['intercept','x1','x2']
for iz in range(2):
    for iv in range(2):
        a = ax[iz, iv]; idx = iv*nz + iz
        a.plot(out_s['Deltadraw'][:, idx], color='steelblue', lw=.4)
        a.axhline(Delta_true[iz, iv], color='firebrick', lw=2)
        a.set_title('Delta[%s, %s]' % (names_z[iz], names_v[iv]), fontsize=9)
plt.suptitle('Synthetic Delta trace (red = truth)'); plt.tight_layout()
plt.savefig('hlm_syn_trace.png', dpi=120, bbox_inches='tight'); plt.show()

fig, ax = plt.subplots(1, 3, figsize=(13, 4))
for j in range(3):
    ax[j].scatter(Beta_true[:, j], beta_pm_s[:, j], s=18, color='steelblue', alpha=.7)
    lim=[min(Beta_true[:,j].min(),beta_pm_s[:,j].min()), max(Beta_true[:,j].max(),beta_pm_s[:,j].max())]
    ax[j].plot(lim, lim, 'k--', lw=1); ax[j].set_title('beta recovery: %s'%names_v[j], fontsize=10)
    ax[j].set_xlabel('true'); ax[j].set_ylabel('posterior mean')
plt.tight_layout(); plt.savefig('hlm_syn_recovery.png', dpi=120, bbox_inches='tight'); plt.show()

Top — traces for four entries of the population matrix $\Delta$, with the truth in red: the chains sit right on the true values. Bottom — unit-level recovery, each unit's posterior-mean $\beta_i$ against its true value, scattered tightly along the 45° line (correlation 0.999). Both levels of the hierarchy — the population surface $\Delta'z_i$ and the 60 unit-level regressions — are recovered correctly before the model meets real data.

4. Cheese scanner data

bayesm::cheese (exported once to cheese_raw.csv): weekly VOLUME, shelf PRICE, and DISP (feature/display intensity) for 88 retailers. We model $\log(\text{VOLUME})\sim$ intercept $+\log(\text{PRICE})+\text{DISP}$ per store.

R_EXE = r'C:\Program Files\R\R-4.6.0\bin\x64\Rscript.exe'
if not os.path.exists('cheese_raw.csv'):
    open('_x.R','w').write("userLib<-file.path(Sys.getenv('USERPROFILE'),'R','win-library','4.6'); .libPaths(c(userLib,.libPaths())); library(bayesm); data(cheese); write.csv(cheese,'cheese_raw.csv',row.names=FALSE)")
    subprocess.run([R_EXE,'_x.R'], capture_output=True, text=True); os.remove('_x.R')
cheese = pd.read_csv('cheese_raw.csv')
retailers = sorted(cheese['RETAILER'].unique()); nreg_c = len(retailers)
cnt = cheese['RETAILER'].value_counts()
print('Retailers: %d   total obs: %d   obs/store min=%d max=%d median=%.0f' %
      (nreg_c, len(cheese), cnt.min(), cnt.max(), cnt.median()))

Retailers: 88   total obs: 5555   obs/store min=52 max=68 median=61

Exploratory data analysis

lv = np.log(cheese['VOLUME']); lp = np.log(cheese['PRICE']); dp = cheese['DISP']
fig, ax = plt.subplots(2, 3, figsize=(14, 7))
ax[0,0].hist(lv, bins=50, color='#4477aa'); ax[0,0].set_title('log(VOLUME)')
ax[0,1].hist(cheese['PRICE'], bins=50, color='#aa4444'); ax[0,1].set_title('PRICE (shelf $)')
ax[0,2].hist(dp, bins=50, color='#44aa77'); ax[0,2].set_title('DISP (feature intensity)')
ax[1,0].plot(lp, lv, '.', ms=2, alpha=.2, color='#333333'); ax[1,0].set_title('log(VOL) vs log(PRICE)')
b=np.polyfit(lp, lv, 1); xs=np.array([lp.min(),lp.max()]); ax[1,0].plot(xs, b[0]*xs+b[1],'firebrick',lw=2)
ax[1,0].set_xlabel('log(PRICE)'); ax[1,0].set_ylabel('log(VOLUME)')
ax[1,1].plot(dp, lv, '.', ms=2, alpha=.2, color='#333333'); ax[1,1].set_title('log(VOL) vs DISP')
b2=np.polyfit(dp, lv, 1); xs2=np.array([dp.min(),dp.max()]); ax[1,1].plot(xs2, b2[0]*xs2+b2[1],'steelblue',lw=2)
psd = cheese.groupby('RETAILER')['PRICE'].std()
ax[1,2].hist(psd.dropna(), bins=30, color='#aa8844'); ax[1,2].set_title('within-store PRICE sd')
plt.suptitle('Cheese Exploratory data analysis'); plt.tight_layout(); plt.savefig('hlm_eda.png', dpi=120, bbox_inches='tight'); plt.show()

The bayesm::cheese panel: 88 retailers, 5,555 store-weeks (52–68 observations each). Top row — the marginals of log(VOLUME), shelf PRICE and DISP (feature/display intensity). Bottom row — log-volume against log-price (the downward-sloping demand relation we will estimate, red fit) and against DISP (promotions lift volume), and then the panel that matters most: the within-store price SD. Each store sees only a narrow spread of its own prices, so there is very little variation from which to identify a store-specific elasticity — which is exactly why pooling will be needed.

# build per-store regression data + OLS baseline
regdata_c = []; ols_c = np.zeros((nreg_c, 3)); ssq_c = np.zeros(nreg_c)
for i, ret in enumerate(retailers):
    sub = cheese[cheese['RETAILER']==ret]
    y = np.log(sub['VOLUME'].values); X = np.column_stack([np.ones(len(sub)), np.log(sub['PRICE'].values), sub['DISP'].values])
    regdata_c.append({'y': y, 'X': X})
    b = np.linalg.lstsq(X, y, rcond=None)[0]; ols_c[i] = b
    r = y - X@b; ssq_c[i] = max(float(r@r)/max(len(y)-3,1), 1e-3)
print('OLS log-price elasticity across stores: mean=%.3f sd=%.3f range=[%.2f, %.2f]' %
      (ols_c[:,1].mean(), ols_c[:,1].std(), ols_c[:,1].min(), ols_c[:,1].max()))

OLS log-price elasticity across stores: mean=-2.266 sd=1.765 range=[-12.65, 2.86]

OLS baseline — why pool?

fig, ax = plt.subplots(1, 3, figsize=(13, 4))
labs = ['Intercept','log-Price elasticity','Display effect']
for j in range(3):
    d = ols_c[:, j]; d = d[np.isfinite(d)]
    ax[j].hist(d, bins=30, color='#aaaaaa'); ax[j].set_title(labs[j])
    ax[j].axvline(d.mean(), color='firebrick', lw=2); ax[j].axvline(np.median(d), color='steelblue', lw=2, ls='--')
    ax[j].set_xlabel('OLS estimate')
ax[2].legend(['mean','median'], fontsize=8)
plt.suptitle('Per-store OLS estimates — wide spread motivates pooling'); plt.tight_layout()
plt.savefig('hlm_ols.png', dpi=120, bbox_inches='tight'); plt.show()

Run a separate OLS regression per store and the estimates fly apart. The log-price elasticity has mean −2.27 but an SD of 1.77 and a range of [−12.65, +2.86] — some stores look absurdly elastic, and others come out positively sloped, which demand theory forbids. With only ~60 collinear observations apiece (and little within-store price variation, as the previous figure showed), unpooled OLS is largely fitting noise. This is the motivation for the hierarchy.

5. Model A — hierarchy without unit covariates

$Z=\mathbf 1$ (one population mean $\Delta$ shared by all stores). Each store's $\beta_i$ shrinks toward that pooled mean by an amount set by $V_\beta$ and the store's own precision.

Z_A = np.ones((nreg_c, 1))
Prior_A = dict(Deltabar=np.zeros((1,3)), A=0.01*np.eye(1), nu=3+3, V=(3+2)*np.eye(3), nu_e=3.0, ssq=ssq_c)
out_A = rhier_linear_model(regdata_c, Z_A, Prior=Prior_A, Mcmc=dict(R=10000, burn=2000, nprint=0), seed=21)
beta_pm_A = out_A['betadraw'].mean(2)                          # (88, 3)
Delta_A = out_A['Deltadraw'].mean(0).reshape(1, 3, order='F')  # (1,3) pooled mean
print('pooled mean (Delta_A): intercept=%.3f  log-price=%.3f  disp=%.3f' % tuple(Delta_A[0]))

pooled mean (Delta_A): intercept=10.299  log-price=-2.152  disp=0.995

Shrinkage — OLS vs Bayes posterior mean

fig, ax = plt.subplots(1, 3, figsize=(13, 4.2)); labs=['intercept','log_price','disp']
for j in range(3):
    ax[j].scatter(ols_c[:, j], beta_pm_A[:, j], s=22, color='steelblue', alpha=.7)
    lim=[np.nanpercentile(ols_c[:,j],1), np.nanpercentile(ols_c[:,j],99)]
    ax[j].plot(lim, lim, 'grey', ls='--'); ax[j].set_xlim(lim)
    ax[j].axvline(Delta_A[0,j], color='firebrick', ls=':'); ax[j].axhline(Delta_A[0,j], color='firebrick', ls=':')
    ax[j].set_xlabel('OLS'); ax[j].set_ylabel('posterior mean'); ax[j].set_title(labs[j])
plt.suptitle('Shrinkage: OLS vs Bayes posterior mean (red dotted = pooled prior mean)'); plt.tight_layout()
plt.savefig('hlm_shrinkage.png', dpi=120, bbox_inches='tight'); plt.show()
print('SD of log-price elasticity:  OLS %.3f  ->  Bayes A %.3f  (shrunk toward pool)' % (ols_c[:,1].std(), beta_pm_A[:,1].std()))

SD of log-price elasticity:  OLS 1.765  ->  Bayes A 0.764  (shrunk toward pool)

The shrinkage picture: each store's Bayesian posterior mean (y-axis) against its own OLS estimate (x-axis), with the pooled population mean marked by the red dotted lines. Points are pulled off the 45° line toward the pooled mean — and the further an OLS estimate strays into the noisy extremes, the harder it is pulled back. The cross-store SD of the log-price elasticity collapses accordingly (printed below the figure), with the whole distribution drawn toward the pooled value of −2.15. This is partial pooling doing exactly its job: borrow strength where a store's own data is thin, defer to it where the data is strong.

6. Model B — hierarchy with a large-market covariate

Add $Z=[\,1,\ \text{large-market dummy}\,]$: the population mean is allowed to differ between large metros and the rest. $\Delta$ row 1 = small-market mean, row 2 = large-market increment.

cities = [r.split(' - ')[0] for r in retailers]
large_pat = ['LOS ANGELES','CHICAGO','NEW YORK','NEW ENGLAND','ATLANTA','PHILADELPHIA','SAN FRANCISCO','DETROIT','DALLAS','HOUSTON','MIAMI']
large_flag = np.array([int(any(p in c for p in large_pat)) for c in cities])
Z_B = np.column_stack([np.ones(nreg_c), large_flag])
print('Large-market stores: %d / %d' % (large_flag.sum(), nreg_c))
Prior_B = dict(Deltabar=np.zeros((2,3)), A=0.01*np.eye(2), nu=3+3, V=(3+2)*np.eye(3), nu_e=3.0, ssq=ssq_c)
out_B = rhier_linear_model(regdata_c, Z_B, Prior=Prior_B, Mcmc=dict(R=10000, burn=2000, nprint=0), seed=22)
beta_pm_B = out_B['betadraw'].mean(2)
Delta_B = out_B['Deltadraw'].mean(0).reshape(2, 3, order='F')

Large-market stores: 24 / 88

Model B — Δ credible intervals

names=['intercept','log_price','disp']; rows=[]
Dd = out_B['Deltadraw']                                        # (R_kept, nz*k) col-major, nz=2,k=3
for iz,zl in enumerate(['small-market mean','large-market increment']):
    for iv in range(3):
        d = Dd[:, iv*2 + iz]; rows.append((zl, names[iv], d.mean(), np.percentile(d,2.5), np.percentile(d,97.5)))
dfB = pd.DataFrame(rows, columns=['Z row','param','mean','q2.5','q97.5'])
print(dfB.round(3).to_string(index=False))
print('\nlarge-market log-price elasticity = %.3f (small %.3f + increment %.3f)' %
      (Delta_B[0,1]+Delta_B[1,1], Delta_B[0,1], Delta_B[1,1]))

                 Z row     param   mean   q2.5  q97.5
     small-market mean intercept 10.083  9.798 10.380
     small-market mean log_price -2.151 -2.393 -1.911
     small-market mean      disp  0.906  0.686  1.138
large-market increment intercept  0.775  0.222  1.336
large-market increment log_price -0.001 -0.461  0.448
large-market increment      disp  0.336 -0.102  0.783

large-market log-price elasticity = -2.152 (small -2.151 + increment -0.001)

7. OLS vs Model A vs Model B

fig, ax = plt.subplots(1, 3, figsize=(14, 4))
ax[0].hist(ols_c[:,1], bins=25, color='#aaaaaa'); ax[0].axvline(ols_c[:,1].mean(), color='black', lw=2); ax[0].set_title('OLS (per store)')
ax[1].hist(beta_pm_A[:,1], bins=25, color='#4477aa'); ax[1].axvline(Delta_A[0,1], color='firebrick', lw=2, ls='--'); ax[1].set_title('Bayes A (no Z)')
ax[2].hist(beta_pm_B[:,1], bins=25, color='#aa4444'); ax[2].axvline(Delta_B[0,1], color='steelblue', lw=2, ls='--'); ax[2].axvline(Delta_B[0,1]+Delta_B[1,1], color='firebrick', lw=2, ls='--'); ax[2].set_title('Bayes B (with Z)')
for a in ax: a.set_xlabel('log-price elasticity')
plt.suptitle('Log-price elasticity: OLS vs Bayes A vs Bayes B'); plt.tight_layout()
plt.savefig('hlm_compare_hist.png', dpi=120, bbox_inches='tight'); plt.show()

ordr = np.argsort(beta_pm_A[:,1]); xs = np.arange(nreg_c)
fig, ax = plt.subplots(figsize=(11, 4.6))
ax.scatter(xs, ols_c[ordr,1], s=10, color='#aaaaaa', label='OLS')
ax.scatter(xs, beta_pm_A[ordr,1], s=14, color='steelblue', label='Bayes A')
ax.scatter(xs, beta_pm_B[ordr,1], s=14, color='firebrick', label='Bayes B')
ax.axhline(0, color='grey', ls='--'); ax.axhline(Delta_A[0,1], color='steelblue', ls=':')
ax.set_xlabel('store (sorted by Bayes A)'); ax.set_ylabel('log-price elasticity')
ax.set_title('Per-store log-price elasticity: OLS vs Bayes A vs Bayes B (shrinkage of OLS extremes)'); ax.legend(fontsize=8)
plt.tight_layout(); plt.savefig('hlm_compare_stores.png', dpi=120, bbox_inches='tight'); plt.show()

Top — the log-price elasticity across the 88 stores under the three specifications: OLS (grey) with its heavy tails and impossible positive values; Model A (blue) tightened around the pooled mean −2.15; Model B (red) with the large-market covariate added. Bottom — the same store by store, sorted, showing the OLS extremes being reeled in.

The honest finding from Model B: the large-market increment on the elasticity is −0.001, with a 95% CrI of [−0.46, 0.45] — straddling zero, so big metros are not more price-sensitive than small ones. What does differ is the baseline: the intercept increment is +0.78 [0.22, 1.34], comfortably above zero. Large markets simply sell more cheese at any given price — they do not respond differently to price.

Takeaways

• Synthetic recovery passes: posterior-mean $\Delta$ matches the truth and the Δ traces hover around their true values; unit-level $\beta_i$ are recovered with high correlation.
• Why pool: per-store OLS log-price elasticities are wildly spread (some positive, some hugely negative) because each store has few, collinear observations. The hierarchy shrinks these noisy estimates toward the pooled mean — the shrinkage plot shows posterior means pulled off the 45° line toward the red prior-mean lines, and the cross-store SD of the elasticity collapses.
• Model B lets the population mean differ by market size; the large-market increment's credible interval shows whether big metros have a systematically different price elasticity.
• OLS vs A vs B: the comparison panels show OLS's heavy tails tamed by A, with B further structuring the means by market — the classic partial-pooling story.