Hierarchical Linear Model with Normal Prior -- bayesm (rhierLinearModel)¶

Package: bayesm (Peter Rossi)
Function: rhierLinearModel -- Gibbs sampler for hierarchical linear models
Reference: Rossi, P. E., Allenby, G. M., & McCulloch, R. (2005). Bayesian Statistics and Marketing, Ch. 3. Wiley.
Python companion: hlm_python.ipynb — the same analysis and graphs via the from-scratch hlm_gibbs.py Gibbs sampler.

Model¶

Unit-level (one regression per unit $i = 1,\ldots,n_{reg}$): $$y_i = \mathbf{X}_i \boldsymbol{\beta}_i + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0,\ \tau_i \mathbf{I})$$

Population (hierarchy): $$\boldsymbol{\beta}_i \sim \mathcal{N}(\mathbf{Z}_i \boldsymbol{\Delta},\ \mathbf{V}_{\beta})$$

Object Dim Meaning
$y_i$ $n_i \times 1$ Observations for unit $i$
$\mathbf{X}_i$ $n_i \times p$ Regressors for unit $i$ (can differ across units)
$\boldsymbol{\beta}_i$ $p \times 1$ Unit-specific coefficients
$\mathbf{Z}_i$ $1 \times n_z$ Unit characteristics / covariates
$\boldsymbol{\Delta}$ $n_z \times p$ How unit chars shift the prior mean of $\boldsymbol{\beta}$
$\mathbf{V}_{\beta}$ $p \times p$ Cross-unit heterogeneity covariance

Priors¶

$$\mathrm{vec}(\boldsymbol{\Delta}) \mid \mathbf{V}_{\beta} \sim \mathcal{N}(\mathrm{vec}(\bar{\boldsymbol{\Delta}}),\ \mathbf{V}_{\beta} \otimes \mathbf{A}^{-1})$$ $$\mathbf{V}_{\beta} \sim \mathcal{IW}(\nu,\ \mathbf{V}) \qquad \tau_i \sim \frac{\nu_e \cdot s_i^2}{\chi^2_{\nu_e}}$$

Hyperparameter Default Interpretation
Deltabar $0$ Prior mean of $\Delta$
A $0.01 I$ Prior precision on $\Delta$
nu $p + 3$ $\mathbf{V}_{\beta}$ IW degrees of freedom
V $\nu I$ $\mathbf{V}_{\beta}$ IW scale
nu.e $3$ Error variance df
ssq $\mathrm{var}(y_i)$ Error variance scale (per unit)

Returns: betadraw ($n_{reg} \times p \times R$), Deltadraw, Vbetadraw, taudraw

Notebooks in this series¶

File Language Content
hlm_bayesm.ipynb (this) R HLM with rhierLinearModel -- synthetic data + cheese
hlm_python.ipynb Python Python replication using hlm_gibbs.py
hlm_gibbs.py Python Gibbs sampler module (replicates bayesm)

1. Setup¶

In [1]:
userLib <- file.path(Sys.getenv('USERPROFILE'), 'R', 'win-library', '4.6')
.libPaths(c(userLib, .libPaths()))
library(bayesm)
library(MASS)      # mvrnorm

set.seed(42L)
cat('bayesm version:', as.character(packageVersion('bayesm')), '\n')

bayesm version: 3.1.7

2. Synthetic Data¶

Simulate from the exact model so we can verify posterior recovery.

Design: $n_{reg}=60$ units, $n_{obs}=80$ observations per unit, $p=3$ regressors (intercept, $x_1$, $x_2$), $n_z=2$ unit covariates (constant + one continuous characteristic, centered).

True $\boldsymbol{\Delta}$ ($n_z \times p$):

intercept $x_1$ $x_2$
const 2.0 -1.5 0.8
unit_char 0.5 0.3 -0.2

The base $x_1$ slope is strongly negative ($-1.5$); the base $x_2$ slope is moderately positive ($0.8$). The continuous unit characteristic slightly offsets both: higher unit_char raises the $x_1$ slope by $0.3$ and lowers the $x_2$ slope by $0.2$.

True $\mathbf{V}_{\beta}$: moderate heterogeneity — diagonal $(0.5, 0.4, 0.3)$ with small positive off-diagonal correlations. This produces unit-level betas ranging roughly:

param min max note
intercept 0.24 3.52 driven by unit_char spread
$x_1$ -3.61 0.13 almost entirely negative
$x_2$ -0.17 2.01 mostly positive

Error variance: $\tau = 0.5$ (common across units in the simulation).

In [2]:
nreg <- 60L   # number of units (e.g. stores)
nobs <- 80L   # observations per unit
nvar <- 3L    # regressors per unit: intercept, x1, x2
nz   <- 2L    # unit covariates: constant + one continuous char

# Unit covariates Z (nreg x nz): constant + centered continuous
Z         <- cbind(1, scale(runif(nreg, 0, 5)))
colnames(Z) <- c('const', 'unit_char')

# True Delta (nz x nvar): how Z shifts mean of beta
Delta_true <- matrix(c(
    2.0, -1.5,  0.8,    # const -> (intercept, x1-slope, x2-slope)
    0.5,  0.3, -0.2     # unit_char effect on each beta
), nrow=nz, ncol=nvar, byrow=TRUE)
rownames(Delta_true) <- colnames(Z)
colnames(Delta_true) <- c('intercept','x1','x2')

# True Vbeta (pxp): cross-unit heterogeneity
Vbeta_true <- matrix(c(0.5, 0.2, 0.0,
                        0.2, 0.4, 0.1,
                        0.0, 0.1, 0.3), nvar, nvar)

# Draw unit-level betas from N(Z*Delta, Vbeta)
Beta_true <- Z %*% Delta_true +
             mvrnorm(nreg, mu=rep(0, nvar), Sigma=Vbeta_true)

# Build regdata: list of lists with X and y per unit
tau_true  <- 0.5   # common error variance for simulation
regdata_s <- vector('list', nreg)
for (i in seq_len(nreg)) {
    Xi <- cbind(1, matrix(rnorm(nobs * (nvar-1L)), nobs, nvar-1L))
    yi <- Xi %*% Beta_true[i, ] + rnorm(nobs, sd=sqrt(tau_true))
    regdata_s[[i]] <- list(y=yi, X=Xi)
}

cat('nreg:', nreg, '  nobs per unit:', nobs, '  nvar:', nvar, '  nz:', nz, '\n')
cat('\nTrue Delta (how unit chars shift beta means):\n')
print(round(Delta_true, 3))
cat('\nTrue Vbeta (cross-unit heterogeneity):\n')
print(round(Vbeta_true, 3))
cat('\nBeta_true range (', nreg, 'units):\n')
print(round(apply(Beta_true, 2, range), 3))

nreg: 60   nobs per unit: 80   nvar: 3   nz: 2 

True Delta (how unit chars shift beta means):
          intercept   x1   x2
const           2.0 -1.5  0.8
unit_char       0.5  0.3 -0.2

True Vbeta (cross-unit heterogeneity):
     [,1] [,2] [,3]
[1,]  0.5  0.2  0.0
[2,]  0.2  0.4  0.1
[3,]  0.0  0.1  0.3

Beta_true range ( 60 units):
     intercept     x1     x2
[1,]     0.242 -3.608 -0.165
[2,]     3.523  0.134  2.006

3. rhierLinearModel — Synthetic Data¶

Pass Z so the sampler jointly estimates $\boldsymbol{\Delta}$, $\mathbf{V}_{\beta}$, unit-level betas $\{\boldsymbol{\beta}_i\}$, and error variances $\{\tau_i\}$.

Prior hyperparameters used:

Param Value Meaning
Deltabar $\mathbf{0}_{2\times 3}$ Diffuse: no prior knowledge of $\Delta$
A $0.01\,I_2$ Very loose prior precision on $\Delta$
nu.e 3 Minimal df for error variance prior
nu 6 ($= p+3$) Minimal proper IW prior on $\mathbf{V}_{\beta}$
V $5\,I_3$ ($= (p+2)\,I$) Prior mean $\mathbf{V}_{\beta} \approx 2.5\,I$

With $N = 60 \times 80 = 4{,}800$ total observations, the likelihood dominates these diffuse priors and the posterior should recover the true parameters.

20,000 draws, 5,000 burn-in. Runtime: ~0.9s.

In [3]:
R_s    <- 20000L
BURN_s <- 5000L

Data_s  <- list(regdata=regdata_s, Z=Z)
Prior_s <- list(Deltabar = matrix(0, nz, nvar),
                A        = 0.01 * diag(nz),
                nu.e     = 3L,
                nu       = nvar + 3L,
                V        = (nvar + 2L) * diag(nvar))
Mcmc_s  <- list(R=R_s, keep=1L, nprint=5000L)

t0    <- proc.time()[[3]]
set.seed(42L)
out_s <- rhierLinearModel(Data=Data_s, Prior=Prior_s, Mcmc=Mcmc_s)
cat(sprintf('\n%d draws in %.1fs\n', R_s, proc.time()[[3]] - t0))

 
Starting Gibbs Sampler for Linear Hierarchical Model
    60  Regressions
    2  Variables in Z (if 1, then only intercept)
 
Prior Parms: 
Deltabar
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
A
     [,1] [,2]
[1,] 0.01 0.00
[2,] 0.00 0.01
nu.e (d.f. parm for regression error variances)=  3
Vbeta ~ IW(nu,V)
nu =  6
V 
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]    0    5    0
[3,]    0    0    5
 
MCMC parms: 
R=  20000  keep=  1  nprint=  5000
 
 MCMC Iteration (est time to end - min) 
 5000 (0.0)
 10000 (0.0)
 15000 (0.0)
 20000 (0.0)
 Total Time Elapsed: 0.02 

20000 draws in 0.9s
In [4]:
keep_s <- seq(BURN_s + 1L, R_s)   # post-burn draws

# Posterior mean of Delta
Delta_post <- matrix(colMeans(out_s$Deltadraw[keep_s, ]), nz, nvar)
rownames(Delta_post) <- rownames(Delta_true)
colnames(Delta_post) <- colnames(Delta_true)

cat('Posterior mean Delta:\n'); print(round(Delta_post, 3))
cat('\nTrue Delta:\n');          print(round(Delta_true, 3))
cat('\nMax |error|:', round(max(abs(Delta_post - Delta_true)), 4), '\n')

# Posterior mean of Vbeta
Vbeta_post <- matrix(colMeans(out_s$Vbetadraw[keep_s, ]), nvar, nvar)
cat('\nPosterior mean Vbeta:\n'); print(round(Vbeta_post, 3))
cat('\nTrue Vbeta:\n');           print(round(Vbeta_true, 3))

Posterior mean Delta:
          intercept     x1     x2
const         2.056 -1.479  0.729
unit_char     0.519  0.366 -0.128

True Delta:
          intercept   x1   x2
const           2.0 -1.5  0.8
unit_char       0.5  0.3 -0.2

Max |error|: 0.0715 

Posterior mean Vbeta:
       [,1]  [,2]   [,3]
[1,]  0.486 0.159 -0.011
[2,]  0.159 0.397  0.073
[3,] -0.011 0.073  0.321

True Vbeta:
     [,1] [,2] [,3]
[1,]  0.5  0.2  0.0
[2,]  0.2  0.4  0.1
[3,]  0.0  0.1  0.3

Recovery Assessment¶

Delta (how unit characteristics shift beta means):

intercept $x_1$ $x_2$
const — posterior 2.056 -1.480 0.729
const — truth 2.000 -1.500 0.800
unit_char — posterior 0.519 0.367 -0.128
unit_char — truth 0.500 0.300 -0.200

Max absolute error across all 6 entries: 0.072. All estimates are within one posterior SD of the truth — very good recovery from 60 units — 80 obs.

$\mathbf{V}_{\beta}$ (cross-unit heterogeneity covariance):

Entry Posterior Truth
diag(1,1) 0.486 0.500
diag(2,2) 0.397 0.400
diag(3,3) 0.321 0.300
off(1,2) 0.159 0.200
off(2,3) 0.073 0.100
off(1,3) -0.011 0.000

Diagonal entries are recovered closely. Off-diagonal covariances are slightly underestimated — typical for IW posteriors with moderate sample size; the IW prior shrinks toward a diagonal matrix. All signs are correct.

In [5]:
# Compare unit betas: posterior mean vs truth
beta_pm_s <- apply(out_s$betadraw[, , keep_s], c(1,2), mean)  # nreg x nvar

cat('Beta recovery (posterior mean vs truth) -- first 8 units:\n')
comp <- cbind(beta_pm_s[1:8, ], Beta_true[1:8, ])
colnames(comp) <- c('pm_int','pm_x1','pm_x2','tr_int','tr_x1','tr_x2')
print(round(comp, 3))

rmse <- sqrt(mean((beta_pm_s - Beta_true)^2))
cat(sprintf('\nRMSE (all %d units, %d params): %.4f\n', nreg, nvar, rmse))

Beta recovery (posterior mean vs truth) -- first 8 units:
     pm_int  pm_x1 pm_x2 tr_int  tr_x1  tr_x2
[1,]  2.473 -0.608 1.145  2.444 -0.602  1.270
[2,]  2.892 -0.391 0.087  2.895 -0.422  0.038
[3,]  2.308 -1.503 1.850  2.195 -1.448  1.734
[4,]  1.693 -1.195 1.166  1.603 -1.167  1.116
[5,]  2.525 -1.277 0.138  2.666 -1.277  0.189
[6,]  1.022 -2.148 0.051  0.873 -2.264 -0.123
[7,]  2.278 -1.965 0.120  2.119 -1.986  0.052
[8,]  0.988 -2.316 0.014  1.023 -2.462  0.089

RMSE (all 60 units, 3 params): 0.0797

Unit-Level Beta Recovery¶

Posterior means vs truth for the first 8 units (all 60 units, 3 params):

Unit pm int tr int pm $x_1$ tr $x_1$ pm $x_2$ tr $x_2$
1 2.473 2.444 -0.608 -0.602 1.145 1.270
2 2.892 2.895 -0.391 -0.422 0.087 0.038
3 2.308 2.195 -1.503 -1.448 1.850 1.734
4 1.693 1.603 -1.195 -1.167 1.166 1.116
5 2.525 2.666 -1.277 -1.277 0.138 0.189
6 1.022 0.873 -2.148 -2.264 0.051 -0.123
7 2.278 2.119 -1.965 -1.986 0.120 0.052
8 0.988 1.023 -2.316 -2.462 0.014 0.089

RMSE across all 60 units and 3 parameters: 0.0797 — excellent recovery.

The unit-level betas are pulled toward $\mathbf{Z}_i \boldsymbol{\Delta}$ by the hierarchy, but with $n_{obs}=80$ observations per unit the likelihood is strong enough that posterior means track the true values closely. Units with more extreme betas (e.g. unit 6: $x_1 = -2.26$) are recovered just as well as near-average units, because the IW/normal hierarchy is flexible enough to accommodate the true spread.

In [6]:
# Trace plots: Delta[1,1] and Delta[2,2] (two key entries)
par(mfrow=c(2,2), mar=c(3,3,2,1))

for (iz in 1:2) {
    for (iv in 1:2) {
        idx   <- (iv-1)*nz + iz
        draws <- out_s$Deltadraw[keep_s, idx]
        plot(draws, type='l', col='steelblue', lwd=0.5,
             main=sprintf('Delta[%s, %s]',
                          rownames(Delta_true)[iz],
                          colnames(Delta_true)[iv]),
             xlab='draw', ylab='')
        abline(h=Delta_true[iz, iv], col='firebrick', lwd=2)
    }
}
No description has been provided for this image

4. Cheese Data¶

Source: bayesm::cheese (Peter Rossi). IRI scanner panel, collected 1989-1992.

Product: Helvatia brand sliced cheese sold in US supermarkets.

Structure: 5,555 store-week observations across 88 grocery retailers in major US markets. The panel is highly balanced: observations per store range from 52 to 68 weeks (median 61), so each retailer contributes a similar amount of information to the hierarchy.

Variable Type Range Description
RETAILER Factor (88 levels) -- Store identifier: CITY - CHAIN format
VOLUME Integer -- Weekly unit sales
PRICE Numeric $1.32 - $4.64 Shelf price (3.5x range across stores/weeks)
DISP Numeric 0 - 1 ACV-weighted display/feature intensity

Markets represented (46 distinct cities extracted from retailer names; e.g. ATLANTA - KROGER CO, BALTI/WASH - SAFEWAY): Albany NY, Atlanta, Baltimore/Washington, Birmingham, Boston (New England), Chicago, Dallas, Denver, Detroit, Houston, Los Angeles, Miami, Minneapolis, New York, Philadelphia, Pittsburgh, and others.

Why hierarchical? Each store has its own price sensitivity ($\beta_{1i}$) and display response ($\beta_{2i}$). With only 52-68 observations per store, OLS estimates are noisy. The HLM borrows strength across all 88 stores: stores with noisier data are pulled more strongly toward the cross-store mean, while stores with cleaner price/display variation retain their own signal.

Unit-level model: $$\log(\text{VOLUME}_{it}) = \beta_{0i} + \beta_{1i}\log(\text{PRICE}_{it}) + \beta_{2i}\,\text{DISP}_{it} + \varepsilon_{it}$$

We expect $\beta_{1i} < 0$ (higher price reduces sales) and $\beta_{2i} > 0$ (display boosts sales). The wide price range ($1.32-$4.64) should identify price elasticity well within each store.

In [7]:
data(cheese)
retailers <- levels(cheese$RETAILER)
nreg_c    <- length(retailers)

cat('Retailers (stores):', nreg_c, '\n')
cat('Total observations:', nrow(cheese), '\n')

# Obs per retailer
cnt <- sapply(retailers, function(r) sum(cheese$RETAILER == r))
cat(sprintf('Obs per store: min=%d  max=%d  median=%.0f\n',
             min(cnt), max(cnt), median(cnt)))

cat('\nSample retailers:\n')
print(head(retailers, 8))
cat('\nPrice range: [', round(min(cheese$PRICE),2), ',',
    round(max(cheese$PRICE),2), ']\n')
cat('DISP range:  [', round(min(cheese$DISP), 2), ',',
    round(max(cheese$DISP), 2), ']\n')

Retailers (stores): 88 

Total observations: 5555 
Obs per store: min=52  max=68  median=61

Sample retailers:
[1] "ALBANY,NY - PRICE CHOPPER"   "ATLANTA - KROGER CO"        
[3] "ATLANTA - WINN DIXIE"        "BALTI/WASH - GIANT FOOD INC"
[5] "BALTI/WASH - SAFEWAY"        "BALTI/WASH - SUPER FRESH"   
[7] "BIRMINGHAM/MONTGOM - BRUNOS" "BIRMINGHAM/MONTGOM - KROGER"

Price range: [ 1.32 , 4.64 ]
DISP range:  [ 0 , 1 ]
In [8]:
# Exploratory data analysis: price and volume distributions
par(mfrow=c(2,3), mar=c(4,4,2,1))

hist(log(cheese$VOLUME), breaks=50, col='#4477aa88', border=NA,
     main='log(VOLUME)', xlab='log units sold')

hist(cheese$PRICE, breaks=50, col='#aa444488', border=NA,
     main='PRICE (shelf)', xlab='dollars')

hist(cheese$DISP, breaks=50, col='#44aa7788', border=NA,
     main='DISP (feature intensity)', xlab='ACV weight')

# log(VOLUME) vs log(PRICE) -- pooled
plot(log(cheese$PRICE), log(cheese$VOLUME), pch='.', col='#33333355',
     xlab='log(PRICE)', ylab='log(VOLUME)',
     main='Volume vs Price (pooled)')
abline(lm(log(VOLUME) ~ log(PRICE), data=cheese), col='firebrick', lwd=2)

# log(VOLUME) vs DISP
plot(cheese$DISP, log(cheese$VOLUME), pch='.', col='#33333355',
     xlab='DISP', ylab='log(VOLUME)',
     main='Volume vs Display')
abline(lm(log(VOLUME) ~ DISP, data=cheese), col='steelblue', lwd=2)

# Price variation across stores
price_sd <- sapply(retailers, function(r) sd(cheese$PRICE[cheese$RETAILER==r]))
hist(price_sd, breaks=30, col='#aa884488', border=NA,
     main='Within-store price SD', xlab='SD of PRICE')
No description has been provided for this image
In [9]:
# Build regdata: log(VOLUME) ~ intercept + log(PRICE) + DISP
nvar_c   <- 3L   # intercept, log(price), disp
regdata_c <- vector('list', nreg_c)
for (i in seq_len(nreg_c)) {
    sub <- cheese[cheese$RETAILER == retailers[i], ]
    y   <- log(sub$VOLUME)
    X   <- cbind(1, log(sub$PRICE), sub$DISP)
    colnames(X) <- c('intercept', 'log_price', 'disp')
    regdata_c[[i]] <- list(y=y, X=X)
}

# OLS per retailer (frequentist reference)
ols_c <- t(sapply(regdata_c, function(d) {
    coef(lm(d$y ~ d$X - 1))
}))
colnames(ols_c) <- c('intercept','log_price','disp')

cat('OLS log-price elasticity across stores:\n')
cat(sprintf('  mean=%.3f  sd=%.3f  range=[%.3f, %.3f]\n',
             mean(ols_c[,2]), sd(ols_c[,2]),
             min(ols_c[,2]), max(ols_c[,2])))

OLS log-price elasticity across stores:
  mean=-2.266  sd=1.775  range=[-12.651, 2.859]

OLS Baseline -- Why Pooling Matters¶

Per-store OLS results across 86 stores with clean fits (2 stores dropped -- insufficient price variation for OLS identification):

intercept log-price display
mean 10.365 -2.230 1.791
SD 1.755 1.780 3.591
p5 7.874 -4.269 -0.076
p95 12.432 -0.673 5.127

The cross-store SD of the log-price elasticity (1.78) is 80% as large as its mean (-2.23) -- OLS estimates are overwhelmed by noise. The 5th-95th percentile range [-4.27, -0.67] spans more than 3.5 units, and the full range includes values as extreme as -12.65 and +2.86 (economically implausible).

The display coefficient is even noisier (SD = 3.59), with some stores showing negative display effects -- again a noise artifact from limited within-store display variation.

This is exactly the setting where hierarchical pooling adds value: with only ~61 observations per store, OLS cannot reliably separate price/display effects from week-to-week volume noise. The HLM shrinks extreme estimates toward the cross-store mean, with shrinkage strength determined by each store's data quality.

In [10]:
# OLS coefficient distributions -- motivation for hierarchical model
# (some stores may have collinear X; drop NA rows)
ols_c_clean <- ols_c[complete.cases(ols_c), ]
cat(sprintf('Stores with clean OLS: %d / %d\n', nrow(ols_c_clean), nreg_c))

par(mfrow=c(1,3), mar=c(4,4,2,1))
param_labels <- c('Intercept', 'log-Price elasticity', 'Display effect')

for (j in 1:3) {
    hist(ols_c_clean[,j], breaks=30, col='#aaaaaa88', border=NA,
         main=param_labels[j], xlab='OLS estimate')
    abline(v=mean(ols_c_clean[,j]),   col='firebrick', lwd=2)
    abline(v=median(ols_c_clean[,j]), col='steelblue', lwd=2, lty=2)
}
legend('topright', c('mean','median'), col=c('firebrick','steelblue'),
       lwd=2, lty=c(1,2), bty='n')
mtext('OLS per-store estimates -- wide spread motivates pooling',
       outer=TRUE, line=-1.5, cex=0.85)

cat('OLS cross-store summary:\n')
print(round(apply(ols_c_clean, 2, function(x)
      c(mean=mean(x), sd=sd(x),
        p5=quantile(x,.05), p95=quantile(x,.95))), 3))

Stores with clean OLS: 86 / 88
OLS cross-store summary:
        intercept log_price   disp
mean       10.365    -2.230  1.791
sd          1.755     1.780  3.591
p5.5%       7.874    -4.269 -0.076
p95.95%    12.432    -0.673  5.127
No description has been provided for this image

5. Model A -- Hierarchy Without Unit Covariates¶

Omit Z: bayesm inserts a vector of ones automatically and prints Z not specified -- putting in iota. $\boldsymbol{\Delta}$ collapses to $1 \times 3$ -- a single pooled mean vector shared by all 88 stores. This is the simplest shrinkage model.

Prior hyperparameters:

Param Value Note
Deltabar $\mathbf{0}_{1\times 3}$ No prior knowledge of pooled mean
A $[0.01]$ Scalar -- diffuse prior on the 1-D Delta
nu.e 3 Minimal df for store error variances
nu 6 ($= p+3$) Minimal proper IW prior on $\mathbf{V}_{\beta}$
V $5\,I_3$ Prior mean $\mathbf{V}_{\beta} \approx 2.5\,I$

25,000 draws, 5,000 burn-in. Runtime: ~1.5s on 88 stores -- 61 obs (5,555 total).

With diffuse priors and 88 stores, Delta will be essentially the cross-store OLS average. The Bayesian value is the full posterior over all 88 store-level betas, regularised toward that pooled mean.

In [11]:
R_c    <- 25000L
BURN_c <- 5000L
nvar_c <- 3L

Data_A  <- list(regdata=regdata_c)   # no Z => intercept-only hierarchy
Prior_A <- list(Deltabar = matrix(0, 1L, nvar_c),
                A        = 0.01 * diag(1L),
                nu.e     = 3L,
                nu       = nvar_c + 3L,
                V        = (nvar_c + 2L) * diag(nvar_c))
Mcmc_A  <- list(R=R_c, keep=1L, nprint=5000L)

t0     <- proc.time()[[3]]
out_A  <- rhierLinearModel(Data=Data_A, Prior=Prior_A, Mcmc=Mcmc_A)
cat(sprintf('\nModel A: %d draws in %.1fs\n', R_c, proc.time()[[3]] - t0))

Z not specified -- putting in iota
 
Starting Gibbs Sampler for Linear Hierarchical Model
    88  Regressions
    1  Variables in Z (if 1, then only intercept)
 
Prior Parms: 
Deltabar
     [,1] [,2] [,3]
[1,]    0    0    0
A
     [,1]
[1,] 0.01
nu.e (d.f. parm for regression error variances)=  3
Vbeta ~ IW(nu,V)
nu =  6
V 
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]    0    5    0
[3,]    0    0    5
 
MCMC parms: 
R=  25000  keep=  1  nprint=  5000
 
 MCMC Iteration (est time to end - min) 
 5000 (0.1)
 10000 (0.0)
 15000 (0.0)
 20000 (0.0)
 25000 (0.0)
 Total Time Elapsed: 0.03 

Model A: 25000 draws in 1.6s
In [12]:
keep_c <- seq(BURN_c + 1L, R_c)

# Pooled mean = Delta (1 x nvar)
Delta_A <- matrix(colMeans(out_A$Deltadraw[keep_c, ]), 1L, nvar_c)
colnames(Delta_A) <- c('intercept','log_price','disp')
cat('Posterior mean Delta (= pooled beta mean):\n')
print(round(Delta_A, 3))

# Vbeta
Vbeta_A <- matrix(colMeans(out_A$Vbetadraw[keep_c, ]), nvar_c, nvar_c)
cat('\nPosterior mean Vbeta (cross-store heterogeneity):\n')
print(round(Vbeta_A, 3))

# Unit betas
beta_pm_A <- apply(out_A$betadraw[, , keep_c], c(1,2), mean)
colnames(beta_pm_A) <- c('intercept','log_price','disp')
cat('\nPosterior mean log-price elasticity per store (first 10):\n')
print(round(beta_pm_A[1:10, 2], 3))
cat(sprintf('\nShrinkage check: pooled=%.3f  OLS mean=%.3f  Bayes mean=%.3f\n',
    Delta_A[1,2], mean(ols_c[,2]), mean(beta_pm_A[,2])))

Posterior mean Delta (= pooled beta mean):
     intercept log_price  disp
[1,]    10.289    -2.145 0.987

Posterior mean Vbeta (cross-store heterogeneity):
       [,1]   [,2]   [,3]
[1,]  1.308 -0.695  0.125
[2,] -0.695  0.727 -0.078
[3,]  0.125 -0.078  0.579

Posterior mean log-price elasticity per store (first 10):
 [1] -3.683 -1.672 -1.985 -3.455 -2.196 -2.914 -1.487 -1.084 -2.499 -2.691

Shrinkage check: pooled=-2.145  OLS mean=-2.266  Bayes mean=-2.145

Model A -- Results¶

Pooled mean Delta ($1 \times 3$):

intercept log-price display
Posterior mean 10.289 -2.145 0.987
OLS mean (86 stores) 10.365 -2.230 1.791

The pooled log-price elasticity of -2.145 means a 1% price increase reduces weekly sales by about 2.1% on average across stores -- a moderately elastic product. The display coefficient of 0.987 represents a large positive effect: moving from zero to full display intensity increases log(sales) by ~1.0 (roughly a 170% increase in units).

Posterior mean $\mathbf{V}_{\beta}$ (cross-store heterogeneity):

intercept log-price display
intercept 1.307 -0.695 0.125
log-price -0.695 0.727 -0.078
display 0.125 -0.078 0.579

The diagonal entries confirm substantial genuine heterogeneity across stores: SD of intercepts $\approx 1.14$, SD of price elasticities $\approx 0.85$, SD of display effects $\approx 0.76$. These are not just OLS noise ? the hierarchy has filtered out noise and the remaining spread is real.

The strong negative off-diagonal Vbeta[1,2] = -0.695 (correlation $\approx -0.71$) is a meaningful economic finding: stores with higher baseline sales volume (high intercept) tend to be less price-sensitive. This is consistent with high-traffic stores attracting less price-conscious shoppers, or high-volume stores having more stable pricing that limits within-store elasticity identification.

Shrinkage summary: OLS pooled mean = -2.266, Bayes pooled mean = -2.145 (close, as expected with diffuse priors and 88 stores). But the store-level estimates are heavily regularised: the first 10 Bayes elasticities range from -3.68 to -1.08, compared to the OLS range of -12.65 to +2.86. All Bayes estimates are negative and economically plausible.

In [13]:
# Shrinkage plot: OLS vs Bayes posterior mean, per store
par(mfrow=c(1,3), mar=c(4,4,2,1))
param_names <- c('intercept','log_price','disp')

for (j in 1:3) {
    plot(ols_c[, j], beta_pm_A[, j],
         pch=16, cex=0.7, col='steelblue',
         xlab='OLS', ylab='Posterior mean',
         main=param_names[j])
    abline(0, 1, col='grey50', lty=2)
    abline(v=Delta_A[1, j], col='firebrick', lty=3, lwd=1.5)
    abline(h=Delta_A[1, j], col='firebrick', lty=3, lwd=1.5)
}
mtext('OLS vs Bayes posterior mean (dashed = pooled prior mean)',
       outer=TRUE, line=-1.5, cex=0.9)
No description has been provided for this image

Model A -- Interpretation¶

Pooled mean: with diffuse priors and 88 stores, Delta converges to essentially the cross-store average: intercept ~10.3, log-price ~-2.1, display ~1.0. All 88 stores are shrunk toward this common target.

Shrinkage pattern: the OLS vs Bayes scatter plot shows the classic shrinkage picture -- points lying on a line with slope < 1, rotated toward the pooled mean (dashed lines). Extreme OLS values (e.g. elasticities beyond -5 or above 0) are pulled sharply toward -2.1. Stores near the pooled mean are barely moved.

Vbeta: the posterior mean covariance has SD ~0.85 for log-price and ~0.76 for display across stores. This genuine heterogeneity justifies the unit-level model -- a pooled OLS would be misspecified. The negative intercept-price correlation (-0.71) is an important structural feature recovered only by the full Bayesian hierarchy.

6. Model B -- Hierarchy With Unit Covariates (Z)¶

Add a store-level covariate: large metro market (yes/no). Cities are extracted from the CITY - CHAIN retailer name format. The dataset covers 46 distinct markets; we classify 24 of 88 stores as large metros using grepl partial matching on 11 city patterns: Atlanta, Chicago, Dallas, Detroit, Houston, Los Angeles, Miami, New England (Boston metro), New York, Philadelphia, San Francisco.

Note: BOSTON appears as a separate smaller-market label in the data, distinct from NEW ENGLAND (NORTH) (the broader Boston metro). Using grepl rather than exact %in% matching correctly captures NEW YORK (NEW) and NEW ENGLAND (NORTH) which have suffixes in the data.

$\mathbf{Z}_i = (1,\ \text{large}_i)$, so $\boldsymbol{\Delta}$ is $2 \times 3$: row 1 = base mean of betas for small-market stores; row 2 = shift in each beta for large-market stores.

In [14]:
# Extract city from retailer name (format: 'CITY - CHAIN')
cities <- sub(' - .*$', '', retailers)

# 46 distinct markets in the data; classify top-population metros as 'large'.
# Use grepl (partial match) because some city labels have suffixes
# e.g. 'NEW YORK (NEW)', 'NEW ENGLAND (NORTH)' for Boston metro.
large_pat <- paste(c('LOS ANGELES', 'CHICAGO', 'NEW YORK',
                     'NEW ENGLAND', 'ATLANTA', 'PHILADELPHIA',
                     'SAN FRANCISCO', 'DETROIT', 'DALLAS',
                     'HOUSTON', 'MIAMI'),
                   collapse='|')
large_flag <- as.integer(grepl(large_pat, cities))

# Z: nreg x 2 (constant + large-market dummy)
Z_c <- cbind(1, large_flag)
colnames(Z_c) <- c('const', 'large_market')

cat('Large-market stores:', sum(large_flag), '/ 88\n')
cat('Large markets matched:\n')
print(sort(unique(cities[large_flag == 1])))
cat('\nSmall markets (first 10):\n')
print(head(sort(unique(cities[large_flag == 0])), 10))

Large-market stores: 24 / 88
Large markets matched:
 [1] "ATLANTA"             "CHICAGO"             "DALLAS/FT. WORTH"   
 [4] "DETROIT"             "HOUSTON"             "LOS ANGELES"        
 [7] "MIAMI"               "NEW ENGLAND (NORTH)" "NEW YORK (NEW)"     
[10] "PHILADELPHIA"        "SAN FRANCISCO"      

Small markets (first 10):
 [1] "ALBANY,NY"          "BALTI/WASH"         "BIRMINGHAM/MONTGOM"
 [4] "BOSTON"             "BUFFALO/ROCHESTER"  "CHARLOTTE"         
 [7] "CINCINNATI"         "CLEVELAND"          "COLUMBUS,OH"       
[10] "DENVER"
In [15]:
nz_c   <- 2L

Data_B  <- list(regdata=regdata_c, Z=Z_c)
Prior_B <- list(Deltabar = matrix(0, nz_c, nvar_c),
                A        = 0.01 * diag(nz_c),
                nu.e     = 3L,
                nu       = nvar_c + 3L,
                V        = (nvar_c + 2L) * diag(nvar_c))
Mcmc_B  <- list(R=R_c, keep=1L, nprint=5000L)

t0    <- proc.time()[[3]]
set.seed(456L)
out_B <- rhierLinearModel(Data=Data_B, Prior=Prior_B, Mcmc=Mcmc_B)
cat(sprintf('\nModel B: %d draws in %.1fs\n', R_c, proc.time()[[3]] - t0))

 
Starting Gibbs Sampler for Linear Hierarchical Model
    88  Regressions
    2  Variables in Z (if 1, then only intercept)
 
Prior Parms: 
Deltabar
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
A
     [,1] [,2]
[1,] 0.01 0.00
[2,] 0.00 0.01
nu.e (d.f. parm for regression error variances)=  3
Vbeta ~ IW(nu,V)
nu =  6
V 
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]    0    5    0
[3,]    0    0    5
 
MCMC parms: 
R=  25000  keep=  1  nprint=  5000
 
 MCMC Iteration (est time to end - min) 
 5000 (0.1)
 10000 (0.0)
 15000 (0.0)
 20000 (0.0)
 25000 (0.0)
 Total Time Elapsed: 0.03 

Model B: 25000 draws in 1.7s
In [16]:
# Posterior mean Delta (2 x 3)
Delta_B <- matrix(colMeans(out_B$Deltadraw[keep_c, ]), nz_c, nvar_c)
rownames(Delta_B) <- colnames(Z_c)
colnames(Delta_B) <- c('intercept','log_price','disp')
cat('Posterior mean Delta (Model B with Z):\n')
print(round(Delta_B, 3))

cat('\nInterpretation:\n')
cat(sprintf('  Base log-price elasticity (small market): %.3f\n', Delta_B[1,2]))
cat(sprintf('  Large-market shift in elasticity:         %.3f\n', Delta_B[2,2]))
cat(sprintf('  Large-market log-price elasticity:        %.3f\n', sum(Delta_B[,2])))

# Vbeta
Vbeta_B <- matrix(colMeans(out_B$Vbetadraw[keep_c, ]), nvar_c, nvar_c)
cat('\nPosterior mean Vbeta (residual heterogeneity after Z):\n')
print(round(Vbeta_B, 3))
cat('\nModel A Vbeta (no Z):\n')
print(round(Vbeta_A, 3))
cat('(Vbeta should shrink in Model B if Z explains some heterogeneity)\n')

Posterior mean Delta (Model B with Z):
             intercept log_price  disp
const           10.074    -2.141 0.904
large_market     0.782    -0.004 0.306

Interpretation:
  Base log-price elasticity (small market): -2.141
  Large-market shift in elasticity:         -0.004
  Large-market log-price elasticity:        -2.145

Posterior mean Vbeta (residual heterogeneity after Z):
       [,1]   [,2]   [,3]
[1,]  1.182 -0.689  0.079
[2,] -0.689  0.723 -0.073
[3,]  0.079 -0.073  0.557

Model A Vbeta (no Z):
       [,1]   [,2]   [,3]
[1,]  1.308 -0.695  0.125
[2,] -0.695  0.727 -0.078
[3,]  0.125 -0.078  0.579
(Vbeta should shrink in Model B if Z explains some heterogeneity)

Model B -- Results¶

Posterior mean Delta ($2 \times 3$):

intercept log-price display
const (small market base) 10.074 -2.141 0.904
large_market shift +0.782 -0.004 +0.306
large market total 10.856 -2.145 1.210

The price elasticity shift is essentially zero (-0.004): large and small markets have virtually identical price sensitivity (-2.145 vs -2.141). This is a striking reversal from a previous (misclassified) run that found a -0.476 shift -- that result was an artifact of classifying only 3 cities (LA, Chicago, Atlanta) as large markets instead of all 11 major metros. With the correct 24/88 classification, the price-sensitivity difference disappears.

Large markets have credibly higher baseline sales (+0.782 intercept shift) -- consistent with higher population density and store traffic.

Display shift = +0.306: large-market stores may be more display-responsive, but the CrI below determines whether this is credible.

Vbeta -- residual heterogeneity after controlling for market size:

Entry Model A (no Z) Model B (with Z) Reduction
intercept var 1.308 1.182 9.6%
log-price var 0.727 0.723 0.6%
display var 0.579 0.557 3.8%
intercept-price cov -0.695 -0.689 0.9%

The intercept variance shrinks by 9.6% -- the large-market dummy explains a meaningful share of the baseline sales variation across stores. Price sensitivity variance barely moves (0.6%), consistent with the near-zero price shift in Delta. Most store-level heterogeneity remains idiosyncratic and unexplained by market size alone.

In [17]:
# 95% credible intervals for Delta rows
cat('95% CrI for large-market shift (row 2 of Delta):\n')
for (j in 1:nvar_c) {
    idx   <- (j-1L)*nz_c + 2L   # row 2 column j
    draws <- out_B$Deltadraw[keep_c, idx]
    cat(sprintf('  %s: mean=%.3f  [%.3f, %.3f]\n',
                colnames(Delta_B)[j],
                mean(draws),
                quantile(draws, 0.025),
                quantile(draws, 0.975)))
}

95% CrI for large-market shift (row 2 of Delta):
  intercept: mean=0.782  [0.240, 1.322]
  log_price: mean=-0.004  [-0.441, 0.436]
  disp: mean=0.306  [-0.118, 0.733]

Large-Market Shift -- Credible Intervals¶

Parameter Posterior mean 95% CrI Excludes zero?
intercept +0.782 [+0.240, +1.322] Yes
log-price -0.004 [-0.441, +0.436] No
display +0.306 [-0.118, +0.733] No

Intercept shift is credible: large-metro stores have genuinely higher baseline log-sales -- the 95% CrI [0.24, 1.32] excludes zero comfortably. This reflects higher population density and store traffic in major metros.

Price elasticity shift: no credible difference. Mean = -0.004, CrI [-0.44, +0.44] is nearly symmetric around zero. Once all 11 major metro markets are correctly included (24/88 stores), the price sensitivity of large and small markets is statistically indistinguishable.

Display shift: not credible. Mean = +0.306, but CrI [-0.12, +0.73] includes zero. There is a directional signal that large-market stores respond more to display, but the evidence is too weak to act on.

Summary: the only robust finding from Model B is the intercept shift -- large markets sell more. Market size does not predict how stores respond to price changes or promotional display. A richer Z matrix (chain identity, store format, local competition index) would be needed to explain the substantial residual heterogeneity remaining in Vbeta.

Model B -- Interpretation¶

Delta row 1 (const): base-level betas for small-market stores -- intercept 10.074, log-price -2.141, display 0.904.

Delta row 2 (large_market): the shift for large-metro stores. The intercept shift (+0.782) is the main finding: large markets have higher baseline sales. The price shift is negligible (-0.004) -- contrary to the competitive-intensity hypothesis, large-metro stores are not detectably more price-sensitive once all 11 major metros are included.

Vbeta shrinks mainly on the intercept (9.6%) but barely on price (0.6%). Market size explains baseline volume differences but not cross-store variation in price sensitivity -- that variation is store-specific.

Practical implication: market-size-based price differentiation is not supported by this model. Store-specific elasticity estimates from Model A (or a richer Z with chain/format covariates) would be more actionable for pricing decisions.

7. Model Comparison — OLS vs Model A vs Model B¶

Compare store-level log-price elasticities across estimation approaches. Bayesian shrinkage pulls extreme OLS estimates toward the pooled mean; Model B further conditions shrinkage on market size.

In [18]:
beta_pm_B <- apply(out_B$betadraw[, , keep_c], c(1,2), mean)

# Summary table: price elasticity by market type
cat('Log-price elasticity summary (beta[:,2])\n')
cat(sprintf('%-20s  %8s  %8s  %8s\n', 'Group', 'OLS', 'Bayes A', 'Bayes B'))
cat(strrep('-', 52), '\n')

groups <- list(
    'All stores'    = seq_len(nreg_c),
    'Large markets' = which(large_flag == 1),
    'Small markets' = which(large_flag == 0)
)
for (nm in names(groups)) {
    idx <- groups[[nm]]
    cat(sprintf('%-20s  %8.3f  %8.3f  %8.3f\n', nm,
                mean(ols_c[idx, 2]),
                mean(beta_pm_A[idx, 2]),
                mean(beta_pm_B[idx, 2])))
}

# Shrinkage: variance across stores
cat('\nCross-store SD of log-price elasticity:\n')
cat(sprintf('  OLS: %.3f   Bayes A: %.3f   Bayes B: %.3f\n',
             sd(ols_c[,2]),
             sd(beta_pm_A[,2]),
             sd(beta_pm_B[,2])))

Log-price elasticity summary (beta[:,2])
Group                      OLS   Bayes A   Bayes B
---------------------------------------------------- 
All stores              -2.266    -2.145    -2.143
Large markets           -2.176    -2.151    -2.146
Small markets           -2.300    -2.143    -2.142

Cross-store SD of log-price elasticity:
  OLS: 1.775   Bayes A: 0.752   Bayes B: 0.750
In [19]:
# Three-panel: OLS vs A vs B for log-price elasticity
par(mfrow=c(1,3), mar=c(4,4,2,1))
col_flag <- ifelse(large_flag==1, 'firebrick', 'steelblue')

# Panel 1: OLS distribution
hist(ols_c[,2], breaks=25, col='#aaaaaa88', main='OLS (per store)',
     xlab='log-price elasticity', border=NA)
abline(v=mean(ols_c[,2]), col='black', lwd=2)

# Panel 2: Bayes A
hist(beta_pm_A[,2], breaks=25, col='#4477aa88', main='Bayes A (no Z)',
     xlab='log-price elasticity', border=NA)
abline(v=Delta_A[1,2], col='firebrick', lwd=2, lty=2)

# Panel 3: Bayes B
hist(beta_pm_B[,2], breaks=25, col='#aa444488', main='Bayes B (with Z)',
     xlab='log-price elasticity', border=NA)
abline(v=Delta_B[1,2], col='steelblue', lwd=2, lty=2)
abline(v=sum(Delta_B[,2]), col='firebrick', lwd=2, lty=2)
legend('topright', c('small','large'), col=c('steelblue','firebrick'),
       lwd=2, lty=2, bty='n', cex=0.8)
No description has been provided for this image
In [20]:
# Side-by-side: posterior mean log-price elasticity, all 88 stores
# sorted by Model A estimate
ord <- order(beta_pm_A[,2])
xs  <- seq_len(nreg_c)

par(mar=c(5,4,3,1))
plot(xs, ols_c[ord,2], pch=16, cex=0.5, col='#aaaaaa',
     ylim=range(c(ols_c[,2], beta_pm_A[,2], beta_pm_B[,2])),
     xlab='Store (sorted by Bayes A estimate)',
     ylab='log-price elasticity',
     main='Per-store log-price elasticity: OLS vs Bayes A vs Bayes B')
points(xs, beta_pm_A[ord,2], pch=16, cex=0.6, col='steelblue')
points(xs, beta_pm_B[ord,2], pch=16, cex=0.6, col='firebrick')
abline(h=0, col='grey60', lty=2)
abline(h=Delta_A[1,2], col='steelblue', lty=3, lwd=1.5)
legend('topleft',
       c('OLS', 'Bayes A (no Z)', 'Bayes B (with Z)', 'Pool mean A'),
       col=c('grey50','steelblue','firebrick','steelblue'),
       pch=c(16,16,16,NA), lty=c(NA,NA,NA,3), lwd=c(NA,NA,NA,1.5),
       pt.cex=0.7, bty='n', cex=0.8)
No description has been provided for this image

Comparison -- Key Takeaways¶

Log-price elasticity by group and method:

Group OLS Bayes A Bayes B
All stores -2.266 -2.145 -2.143
Large markets (24) -2.176 -2.151 -2.146
Small markets (64) -2.300 -2.143 -2.142
Cross-store SD 1.775 0.752 0.750

  1. OLS is noisy -- shrinkage removes 58% of the spread. SD drops from 1.775 (OLS) to 0.752 (Bayes A), a 58% reduction. The remaining SD of 0.75 represents genuine store-level heterogeneity that the hierarchy has filtered from noise.

  2. OLS group means are misleading. OLS suggests large markets are less elastic (-2.18 vs -2.30). After shrinkage, that difference disappears entirely (-2.151 vs -2.143). The OLS group difference was driven by the uneven distribution of extreme noisy estimates across market types, not a real phenomenon.

  3. Models A and B are virtually identical for price sensitivity. SD = 0.752 (A) vs 0.750 (B). The large-market dummy adds no explanatory power for cross-store price elasticity variation. Model B earns its keep only through the credible intercept shift.

  4. Store-specific estimates are the main deliverable. The 88 posterior mean elasticities from Model A span roughly [-3.7, -1.1] -- a 2.6-unit range of genuine heterogeneity. A retailer or manufacturer should use these store-specific estimates rather than any pooled number for pricing or promotional planning.

  5. What Z covariate would help? Chain identity (Kroger vs Safeway vs independent) or store format (warehouse vs conventional) would likely explain more of the residual Vbeta than a simple metro/non-metro split.

8. Save Results¶

In [21]:
# Unit-level posterior means
df_A <- as.data.frame(beta_pm_A)
df_A$retailer    <- retailers
df_A$large_market <- large_flag
df_A$model       <- 'A_no_Z'

df_B <- as.data.frame(beta_pm_B)
colnames(df_B)[1:3] <- c('intercept','log_price','disp')
df_B$retailer     <- retailers
df_B$large_market <- large_flag
df_B$model        <- 'B_with_Z'

write.csv(rbind(df_A, df_B), 'cheese_hlm_betas.csv', row.names=FALSE)

# Delta draws for Python comparison
write.csv(out_A$Deltadraw[keep_c,], 'cheese_hlm_A_Delta_draws.csv', row.names=FALSE)
write.csv(out_B$Deltadraw[keep_c,], 'cheese_hlm_B_Delta_draws.csv', row.names=FALSE)

cat('Saved: cheese_hlm_betas.csv\n')
cat('       cheese_hlm_A_Delta_draws.csv\n')
cat('       cheese_hlm_B_Delta_draws.csv\n')

Saved: cheese_hlm_betas.csv
       cheese_hlm_A_Delta_draws.csv
       cheese_hlm_B_Delta_draws.csv