Hierarchical Linear Model -- Mixture of Normals Prior

rhierLinearMixture (bayesm)

Package: bayesm (Peter Rossi)
Function: rhierLinearMixture
Reference: Rossi, Allenby & McCulloch (2005). Bayesian Statistics and Marketing, Ch. 5.

---

Model

Same unit-level equation as rhierLinearModel: $$y_i = \mathbf{X}_i \boldsymbol{\beta}_i + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0, \tau_i)$$

Heterogeneity distribution replaced by a K-component mixture: $$\boldsymbol{\beta}_i \mid \mathbf{Z}_i \sim \sum_{k=1}^K \pi_k \mathcal{N}(\boldsymbol{\mu}_k + \boldsymbol{\Delta}' \mathbf{z}_i,\ \boldsymbol{\Sigma}_k)$$

Object	Meaning
$K$	Number of mixture components (chosen by analyst)
$\pi_k$	Mixing weights, $\sum_k \pi_k = 1$, prior Dirichlet$(a)$
$\boldsymbol{\mu}_k$	Component-specific mean of $\boldsymbol{\beta}$
$\boldsymbol{\Sigma}_k$	Component-specific covariance, prior IW$(\nu, V)$
$\boldsymbol{\Delta}$	Effect of (centered) unit covariates Z on beta mean

Key difference from single-normal HLM:
Stores are softly assigned to components -- the posterior membership probability $p(\text{component}_k \mid \text{store}_i, \text{data})$ gives automatic store segmentation. Components can have different means and different covariances.

Note on Z: unlike rhierLinearModel, Z should not include an intercept column and should be centered. The mixture means $\boldsymbol{\mu}_k$ serve as component-specific intercepts.

---

Notebooks

Notebook	Heterogeneity	Key output
hlm_bayesm.ipynb	Single normal	Store betas, Delta
hlm_mixture_bayesm.ipynb (this)	Mixture of normals	Store segments, component means

The fitted model

Data. bayesm::cheese — 88 retailers (stores). For store $i$ and week $t$: $y_{it}=\log(\text{VOLUME})$, regressed on $X_{it}=[\,1,\ \log(\text{PRICE}),\ \text{DISP}\,]$ (DISP = feature/display intensity).

Likelihood — a regression per store: $$y_i = X_i\,\beta_i + \varepsilon_i,\qquad \varepsilon_i\sim N(0,\ \sigma_i^2 I).$$

Hierarchy — a mixture of normals on the store coefficients (the new ingredient over the single-normal HLM): each store belongs to one of $K$ latent segments, $$\beta_i\mid s_i=k \sim N(\mu_k + \Delta' z_i,\ \Sigma_k),\qquad s_i\mid\pi\sim\mathrm{Cat}(\pi),$$ with optional unit covariates $z_i$ (e.g. a large-market dummy, §B). $K=1$ is the ordinary single-normal hierarchy; $K>1$ asks whether stores fall into discrete types or are just one (possibly non-normal) population.

Priors: $\mu_k\sim N(\bar\mu{=}0,\ A_\mu^{-1})$, $A_\mu{=}0.1$; $\Sigma_k\sim \mathrm{IW}(\nu{=}6,\ V{=}5I)$; mixing weights $\pi\sim\mathrm{Dir}(a{=}5)$; store error variances $\sigma_i^2\sim \mathrm{IG}(\nu_e{=}3,\ \cdot)$.

1. Setup

userLib <- file.path(Sys.getenv('USERPROFILE'), 'R', 'win-library', '4.6')
.libPaths(c(userLib, .libPaths()))
library(bayesm)

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

bayesm version: 3.1.7 

2. Load Cheese Data

Same data preparation as hlm_bayesm.ipynb Section 4. 88 retailers, log(VOLUME) ~ intercept + log(PRICE) + DISP.

data(cheese)
retailers <- levels(cheese$RETAILER)
nreg_c    <- length(retailers)   # 88
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], ]
    regdata_c[[i]] <- list(
        y = log(sub$VOLUME),
        X = cbind(1, log(sub$PRICE), sub$DISP)
    )
}

# Large-market flag (centered -- required by rhierLinearMixture)
cities     <- sub(' - .*$', '', retailers)
large_pat  <- paste(c('LOS ANGELES','CHICAGO','NEW YORK','NEW ENGLAND',
                       'ATLANTA','PHILADELPHIA','SAN FRANCISCO',
                       'DETROIT','DALLAS','HOUSTON','MIAMI'), collapse='|')
large_raw  <- as.integer(grepl(large_pat, cities))
large_cent <- large_raw - mean(large_raw)   # centered
Z_c        <- matrix(large_cent, ncol=1L)   # nreg x 1, NO intercept

cat('Stores:', nreg_c, '  Large market (raw):', sum(large_raw), '\n')
cat('Z mean after centering:', round(mean(Z_c), 6), '\n')

Stores: 88   Large market (raw): 24 
Z mean after centering: 0 

3. Diagnostic -- Is a Mixture Warranted?

Load the Model A (single-normal) posterior beta means from hlm_bayesm.ipynb and examine the distribution of log-price elasticities across stores. Bimodality or heavy tails would justify a mixture prior.

Result: The 88 store-level price elasticities span [-4.01, -0.63] with mean -2.14 and SD 0.75. The Shapiro-Wilk test gives W=0.9759, p=0.100 -- borderline, just above the 5% threshold. A single normal is not strongly rejected.

Two reasons to try a mixture anyway:

1. Power: With n=88 stores, the Shapiro-Wilk test has limited power to detect mild bimodality. A visual inspection of the density plot may reveal more than the p-value alone.
2. Segmentation: Even when the marginal distribution is unimodal, a mixture model provides automatic store segmentation as a by-product -- each store receives a posterior probability of belonging to each component. This can be substantively useful regardless of the normality test outcome.

betas_hlm <- read.csv('cheese_hlm_betas.csv')
betas_A   <- betas_hlm[betas_hlm$model == 'A_no_Z', ]

price_A <- betas_A$log_price
cat('Model A log-price elasticity (88 stores):\n')
cat(sprintf('  mean=%.3f  sd=%.3f  min=%.3f  max=%.3f\n',
             mean(price_A), sd(price_A), min(price_A), max(price_A)))

par(mfrow=c(1,2), mar=c(4,4,2,1))

# Density plot
d <- density(price_A, bw='SJ')
plot(d, main='Density of price elasticities (Model A)',
     xlab='log-price posterior mean', lwd=2, col='steelblue')
rug(price_A, col='#33333355')

# Normal Q-Q plot
qqnorm(price_A, main='Normal Q-Q plot', pch=16, cex=0.7, col='steelblue')
qqline(price_A, col='firebrick', lwd=2)

# Shapiro-Wilk test
sw <- shapiro.test(price_A)
cat(sprintf('\nShapiro-Wilk test: W=%.4f  p=%.4f\n', sw$statistic, sw$p.value))
cat(if(sw$p.value < 0.05) 'Reject normality -- mixture may help\n'
    else 'Cannot reject normality -- single normal may suffice\n')

Model A log-price elasticity (88 stores):
  mean=-2.145  sd=0.752  min=-4.012  max=-0.628

Shapiro-Wilk test: W=0.9759  p=0.1001
Cannot reject normality -- single normal may suffice

4. rhierLinearMixture -- K=2, No Z

Start with two components and no covariates. This is the direct analogue of Model A from hlm_bayesm.ipynb.

Prior specification:

Parameter	Value	Role
ncomp=2	K=2	Number of mixture components
nu.e=3	3	Error variance prior df
nu=6, V=5*I	--	IW prior on each component covariance $\Sigma_k$
mubar=0	0	Prior mean for each component mean $\mu_k$
Amu=0.1	0.1	Precision for $\mu_k$ prior (default 0.01 is too diffuse)
a=[5,5]	symmetric	Dirichlet prior on mixing weights $\pi$; centered at 50/50

Note on Amu: the default Amu=0.01 implies a very flat prior on component means (SD ~ 10 per parameter). We use Amu=0.1 (SD ~ 3) as a mild tightening that is still quite diffuse relative to the data variation.

R_m    <- 25000L
BURN_m <- 5000L

Data_K2  <- list(regdata=regdata_c)   # no Z
Prior_K2 <- list(
    ncomp  = 2L,
    nu.e   = 3L,
    nu     = nvar_c + 3L,
    V      = (nvar_c + 2L) * diag(nvar_c),
    mubar  = rep(0, nvar_c),
    Amu    = 0.1
)
Mcmc_K2  <- list(R=R_m, keep=1L, nprint=5000L)

set.seed(11L)
t0     <- proc.time()[[3]]
out_K2 <- rhierLinearMixture(Data=Data_K2, Prior=Prior_K2, Mcmc=Mcmc_K2)
cat(sprintf('\nK=2 (no Z): %d draws in %.1fs\n', R_m, proc.time()[[3]] - t0))

Z not specified
 
Starting MCMC Inference for Hierarchical Linear Model:
   Normal Mixture with 2 components for first stage prior
   for  88  cross-sectional units
 
Prior Parms: 
nu.e = 3
nu = 6
V 
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]    0    5    0
[3,]    0    0    5
mubar 
     [,1] [,2] [,3]
[1,]    0    0    0
Amu 
     [,1]
[1,]  0.1
a 
[1] 5 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 

K=2 (no Z): 25000 draws in 2.3s

keep_m <- seq(BURN_m + 1L, R_m)

# Posterior mixing weights
pi_post <- colMeans(out_K2$nmix$probdraw[keep_m, ])
cat('Posterior mixing weights:\n')
cat(sprintf('  Component 1: %.3f\n  Component 2: %.3f\n',
             pi_post[1], pi_post[2]))

# Component means (from compdraw)
mu1_draws <- t(sapply(out_K2$nmix$compdraw[keep_m],
                       function(d) d[[1]]$mu))
mu2_draws <- t(sapply(out_K2$nmix$compdraw[keep_m],
                       function(d) d[[2]]$mu))

mu1_pm <- colMeans(mu1_draws)
mu2_pm <- colMeans(mu2_draws)
names(mu1_pm) <- names(mu2_pm) <- c('intercept','log_price','disp')

cat('\nComponent 1 mean (posterior):\n'); print(round(mu1_pm, 3))
cat('\nComponent 2 mean (posterior):\n'); print(round(mu2_pm, 3))
cat('\nPrice elasticity gap between components:',
     round(abs(mu1_pm[2] - mu2_pm[2]), 3), '\n')

Posterior mixing weights:
  Component 1: 0.683
  Component 2: 0.317

Component 1 mean (posterior):
intercept log_price      disp 
    9.765    -2.056     1.144 

Component 2 mean (posterior):
intercept log_price      disp 
    8.326    -1.876     1.787 

Price elasticity gap between components: 0.18 

K=2 Component Summary

	Component 1	Component 2
Mixing weight	0.683 (68%)	0.317 (32%)
Intercept	9.765	8.326
log-price	-2.056	-1.876
Display	1.144	1.787

Interpretation:

The two components do not separate stores primarily by price sensitivity. The elasticity gap is only 0.18 (roughly 0.24 SD units of the overall distribution), consistent with the diagnostic in Section 3 showing the price elasticities are approximately normal.

Instead, the mixture captures two other dimensions of store heterogeneity:

• Base volume (intercept): Component 1 stores have ~1.4 log-units higher baseline volume than Component 2 stores -- a factor of roughly $e^{1.44} \approx 4\times$ higher unit sales at the same price and no display. Component 1 likely represents higher-volume retailers.
• Display sensitivity (disp): Component 2 stores respond more strongly to promotional displays (1.787 vs 1.144). Smaller-volume stores may rely more heavily on in-store promotion to drive sales.

Both components share a meaningfully negative price elasticity, confirming that price sensitivity is a feature common to all stores rather than a discriminating dimension between segments.

# Unit betas and component membership
beta_pm_K2 <- apply(out_K2$betadraw[, , keep_m], c(1,2), mean)
colnames(beta_pm_K2) <- c('intercept','log_price','disp')

# Posterior probability of component 1 for each store
# (from betadraw -- use the component draws indirectly via density comparison)
# Simpler: use hard classification based on which component mean is closer
dist1 <- abs(beta_pm_K2[,2] - mu1_pm[2])
dist2 <- abs(beta_pm_K2[,2] - mu2_pm[2])
seg   <- ifelse(dist1 < dist2, 1L, 2L)

cat('Store count per component (hard assignment):\n')
print(table(seg))
cat('\nPrice elasticity by component:\n')
for (k in 1:2) {
    idx <- which(seg == k)
    cat(sprintf('  Component %d (%d stores): mean=%.3f  sd=%.3f  range=[%.3f, %.3f]\n',
                k, length(idx),
                mean(beta_pm_K2[idx, 2]), sd(beta_pm_K2[idx, 2]),
                min(beta_pm_K2[idx, 2]), max(beta_pm_K2[idx, 2])))
}

Store count per component (hard assignment):
seg
 1  2 
50 38 

Price elasticity by component:
  Component 1 (50 stores): mean=-2.582  sd=0.490  range=[-3.808, -1.969]
  Component 2 (38 stores): mean=-1.551  sd=0.262  range=[-1.938, -0.915]

Store Segment Membership

	Component 1	Component 2
Stores	50 (57%)	38 (43%)
Price elasticity mean	-2.582	-1.551
Price elasticity SD	0.490	0.262
Price elasticity range	[-3.81, -1.97]	[-1.94, -0.92]

Caution -- hard vs soft assignment: The component posterior means are close: -2.056 vs -1.876 (gap = 0.18). The hard assignment rule -- each store assigned to the nearer component mean -- places the threshold at the midpoint of the distribution (-1.97). Splitting a continuous, approximately normal distribution at its center automatically produces two groups whose store-level means are much further apart (gap = 1.03) than the component means themselves.

This is not evidence of bimodality. Rather, the mixture is providing a soft partition of a unimodal distribution into a more-elastic half and a less-elastic half. Stores near the threshold (-1.94 to -1.97) are genuinely ambiguous; only extreme stores are cleanly assigned.

Substantive reading: Conditional on belonging to each segment, Component 1 stores are substantially more price-sensitive (mean -2.58) and show tighter within-group variation (SD 0.26 vs 0.49 for Component 2). Component 2 stores are less elastic and more dispersed -- consistent with a heterogeneous group of smaller-volume, promotion-driven retailers identified in the component summary above.

par(mfrow=c(1,2), mar=c(4,4,2,1))
cols <- c('steelblue','firebrick')

# Density overlay: Model A vs K=2 components
d_all <- density(beta_pm_K2[,2], bw='SJ')
d1    <- density(beta_pm_K2[seg==1, 2], bw='SJ')
d2    <- density(beta_pm_K2[seg==2, 2], bw='SJ')
ylim  <- range(c(d_all$y, d1$y, d2$y))
plot(d_all, main='K=2 mixture -- price elasticity',
     xlab='log-price', lwd=2, col='black', ylim=ylim)
lines(d1, col=cols[1], lwd=1.8, lty=2)
lines(d2, col=cols[2], lwd=1.8, lty=2)
abline(v=mu1_pm[2], col=cols[1], lty=3)
abline(v=mu2_pm[2], col=cols[2], lty=3)
legend('topleft', c('All','Component 1','Component 2'),
       col=c('black',cols), lwd=2, lty=c(1,2,2), bty='n', cex=0.85)

# Component 1 vs 2 -- OLS vs Bayes K=2
ols_c <- t(sapply(regdata_c, function(d) coef(lm(d$y ~ d$X - 1))))
plot(ols_c[,2], beta_pm_K2[,2],
     col=cols[seg], pch=16, cex=0.8,
     xlab='OLS elasticity', ylab='Bayes K=2 elasticity',
     main='OLS vs Bayes K=2 (colour = component)')
abline(0, 1, lty=2, col='grey50')
legend('topleft', c('Component 1','Component 2'),
       col=cols, pch=16, bty='n', cex=0.85)

5. rhierLinearMixture -- K=2, With Centered Z

Add the centered large-market dummy as a covariate. Delta captures systematic shifts in beta means attributable to market size, layered on top of the mixture structure.

Data_K2Z  <- list(regdata=regdata_c, Z=Z_c)   # centered Z, no intercept
Prior_K2Z <- list(
    ncomp    = 2L,
    nu.e     = 3L,
    nu       = nvar_c + 3L,
    V        = (nvar_c + 2L) * diag(nvar_c),
    mubar    = rep(0, nvar_c),
    Amu      = 0.1,
    deltabar = rep(0, nvar_c),   # vec(Delta) prior mean
    Ad       = 0.01 * diag(nvar_c)
)
Mcmc_K2Z  <- list(R=R_m, keep=1L, nprint=5000L)

set.seed(22L)
t0      <- proc.time()[[3]]
out_K2Z <- rhierLinearMixture(Data=Data_K2Z, Prior=Prior_K2Z, Mcmc=Mcmc_K2Z)
cat(sprintf('\nK=2 (with Z): %d draws in %.1fs\n', R_m, proc.time()[[3]] - t0))

 
Starting MCMC Inference for Hierarchical Linear Model:
   Normal Mixture with 2 components for first stage prior
   for  88  cross-sectional units
 
Prior Parms: 
nu.e = 3
nu = 6
V 
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]    0    5    0
[3,]    0    0    5
mubar 
     [,1] [,2] [,3]
[1,]    0    0    0
Amu 
     [,1]
[1,]  0.1
a 
[1] 5 5
deltabar
[1] 0 0 0
Ad
     [,1] [,2] [,3]
[1,] 0.01 0.00 0.00
[2,] 0.00 0.01 0.00
[3,] 0.00 0.00 0.01
 
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 

K=2 (with Z): 25000 draws in 2.4s

# Delta: effect of large-market on beta means
Delta_K2Z <- matrix(colMeans(out_K2Z$Deltadraw[keep_m, ]),
                    nrow=1L, ncol=nvar_c)
colnames(Delta_K2Z) <- c('intercept','log_price','disp')
cat('Posterior mean Delta (large-market shift, K=2 model):\n')
print(round(Delta_K2Z, 3))

# Component means
mu1z_pm <- colMeans(t(sapply(out_K2Z$nmix$compdraw[keep_m],
                              function(d) d[[1]]$mu)))
mu2z_pm <- colMeans(t(sapply(out_K2Z$nmix$compdraw[keep_m],
                              function(d) d[[2]]$mu)))
names(mu1z_pm) <- names(mu2z_pm) <- c('intercept','log_price','disp')

pi_K2Z <- colMeans(out_K2Z$nmix$probdraw[keep_m, ])
cat('\nMixing weights (K=2 with Z):\n')
cat(sprintf('  Component 1: %.3f   Component 2: %.3f\n', pi_K2Z[1], pi_K2Z[2]))
cat('\nComponent 1 mean:\n'); print(round(mu1z_pm, 3))
cat('\nComponent 2 mean:\n'); print(round(mu2z_pm, 3))

Posterior mean Delta (large-market shift, K=2 model):
     intercept log_price disp
[1,]     0.808     -0.03 0.28

Mixing weights (K=2 with Z):
  Component 1: 0.102   Component 2: 0.898

Component 1 mean:
intercept log_price      disp 
    3.310    -0.799     0.983 

Component 2 mean:
intercept log_price      disp 
   10.188    -2.088     0.995 

K=2 + Z Results

Delta -- large-market effect:

Parameter	Posterior mean	Interpretation
Intercept shift	+0.808	Large markets sell ~$e^{0.81} \approx 2.2\times$ more at baseline
log-price shift	-0.030	Virtually no difference in price sensitivity by market size
Display shift	+0.280	Slightly stronger display response in large markets

Mixture components:

	Component 1	Component 2
Weight	0.102 (~9 stores)	0.898 (~79 stores)
Intercept	3.310	10.188
log-price	-0.799	-2.088
Display	0.983	0.995

Interpretation:

The structure changes dramatically once Z is included. With market size accounted for by Delta, the mixture now separates:

• Component 2 (89.8%): The dominant segment. Price elasticity -2.09, closely matching the no-Z component means and the Model A posterior mean. This is the "normal" cheese retailer behavior.

• Component 1 (10.2%): A small outlier cluster of ~9 stores with very weak price sensitivity (-0.80, not even unit-elastic). These stores may include retailers with little price variation in the data (the near-constant price stores flagged in hlm_bayesm.ipynb OLS diagnostics), or niche channels where cheese is not a price-driven purchase.

Key contrast with no-Z model:

Without Z, the two components split the distribution of price elasticities roughly in half (68/32). Once Z absorbs the market-size dimension, the residual mixture finds a qualitatively different structure: a large homogeneous majority and a small heterodox minority. The display sensitivity (disp) is nearly identical across components (0.983 vs 0.995) -- this dimension is no longer driving the segmentation once the market-size covariate is included.

Delta vs Model B in hlm_bayesm.ipynb: The Delta log-price coefficient (-0.030) replicates the Model B finding that large markets are no more price-sensitive than small ones after correcting the Z covariate definition.

6. rhierLinearMixture -- K=3

Try three components to see if a third segment emerges. With only 88 stores, K=3 may be over-parameterised -- watch for a near-zero mixing weight (degenerate component).

Result: All three components have non-trivial posterior weight:

Component	Weight	~Stores
1	0.374	33
2	0.443	39
3	0.182	16

No component collapses toward zero -- all three carry at least 18% of the mass. The split is more balanced than K=2 (68/32), but this does not imply three real segments: the model-comparison Conclusion shows WAIC prefers K=2 only mildly over a single normal, K=3 does not improve on K=2, and no store commits to a component. The added components are not substantively distinct.

Data_K3  <- list(regdata=regdata_c)
Prior_K3 <- list(
    ncomp = 3L,
    nu.e  = 3L,
    nu    = nvar_c + 3L,
    V     = (nvar_c + 2L) * diag(nvar_c),
    mubar = rep(0, nvar_c),
    Amu   = 0.1
)

set.seed(33L)
t0     <- proc.time()[[3]]
out_K3 <- rhierLinearMixture(Data=Data_K3, Prior=Prior_K3, Mcmc=Mcmc_K2)
cat(sprintf('\nK=3 (no Z): %d draws in %.1fs\n', R_m, proc.time()[[3]] - t0))

pi_K3 <- colMeans(out_K3$nmix$probdraw[keep_m, ])
cat('\nMixing weights (K=3):\n')
cat(sprintf('  Component 1: %.3f\n  Component 2: %.3f\n  Component 3: %.3f\n',
             pi_K3[1], pi_K3[2], pi_K3[3]))
cat(if(min(pi_K3) < 0.05) '\nWarning: near-empty component -- K=3 may be over-specified\n'
    else '\nAll three components have non-trivial weight.\n')

Z not specified
 
Starting MCMC Inference for Hierarchical Linear Model:
   Normal Mixture with 3 components for first stage prior
   for  88  cross-sectional units
 
Prior Parms: 
nu.e = 3
nu = 6
V 
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]    0    5    0
[3,]    0    0    5
mubar 
     [,1] [,2] [,3]
[1,]    0    0    0
Amu 
     [,1]
[1,]  0.1
a 
[1] 5 5 5
 
MCMC Parms: 
R=  25000  keep=  1  nprint=  5000

 MCMC Iteration (est time to end - min) 
 5000 (0.0)
 10000 (0.0)
 15000 (0.0)
 20000 (0.0)
 25000 (0.0)
 Total Time Elapsed: 0.03 

K=3 (no Z): 25000 draws in 2.9s

Mixing weights (K=3):
  Component 1: 0.374
  Component 2: 0.443
  Component 3: 0.182

All three components have non-trivial weight.

7. Comparison -- Single Normal vs Mixture

Compare the distribution of store-level price elasticities across estimation approaches.

Cross-store SD of log-price elasticity:

Model	SD	Change vs Model A
Model A -- single normal	0.752	--
K=2 mixture	0.654	-13%
K=3 mixture	0.672	-11%

All three models agree closely on the broad picture: stores are heterogeneous, with price elasticities spanning roughly -4 to -1. The mixture models show less cross-store variation than Model A, not more. This may seem counterintuitive but has a clear explanation:

• Model A shrinks every store toward the single global mean.
• Mixture models shrink each store toward its local component mean. With fewer stores per component, within-component shrinkage is stronger, and stores end up closer to their respective component mean than they would be to the global mean. The net effect on total SD depends on whether this tighter local pooling dominates the between-component spread -- here it does.

Bottom line: the choice between K=1, K=2, and K=3 has modest impact on the estimated store betas. The primary benefit of the mixture is not that it changes the point estimates dramatically, but that it provides automatic store segmentation -- each store is softly assigned to a component, enabling downstream analyses (targeted pricing, promotional strategy) by segment.

# Load Model A betas
betas_A_pm <- betas_hlm[betas_hlm$model == 'A_no_Z', 'log_price']
cols <- c('steelblue','firebrick')

par(mfrow=c(1,3), mar=c(4,4,2,1))
breaks <- seq(-5, 0.5, by=0.4)

hist(betas_A_pm,    breaks=breaks, col='#4477aa88', border=NA,
     main='Model A (single normal)', xlab='log-price elasticity',
     xlim=c(-5, 0.5))
abline(v=mean(betas_A_pm), col='firebrick', lwd=2)

hist(beta_pm_K2[,2], breaks=breaks, col='#44aa7788', border=NA,
     main='K=2 mixture', xlab='log-price elasticity',
     xlim=c(-5, 0.5))
abline(v=mu1_pm[2], col=cols[1], lwd=2, lty=2)
abline(v=mu2_pm[2], col=cols[2], lwd=2, lty=2)

beta_pm_K3 <- apply(out_K3$betadraw[, , keep_m], c(1,2), mean)
hist(beta_pm_K3[,2], breaks=breaks, col='#aa774488', border=NA,
     main='K=3 mixture', xlab='log-price elasticity',
     xlim=c(-5, 0.5))

cat('Cross-store SD of price elasticity:\n')
cat(sprintf('  Model A:  %.3f\n  K=2 mix:  %.3f\n  K=3 mix:  %.3f\n',
             sd(betas_A_pm), sd(beta_pm_K2[,2]), sd(beta_pm_K3[,2])))

Cross-store SD of price elasticity:
  Model A:  0.752
  K=2 mix:  0.654
  K=3 mix:  0.672

8. Save Results

df_K2 <- data.frame(
    retailer     = retailers,
    large_market = large_raw,
    intercept    = beta_pm_K2[,1],
    log_price    = beta_pm_K2[,2],
    disp         = beta_pm_K2[,3],
    component    = seg,
    model        = 'K2_mix'
)
write.csv(df_K2, 'cheese_hlm_mixture_K2.csv', row.names=FALSE)
cat('Saved: cheese_hlm_mixture_K2.csv\n')

# Component summary
comp_summ <- data.frame(
    component = c(1L, 2L),
    weight    = round(pi_post, 3),
    mu_int    = round(c(mu1_pm[1], mu2_pm[1]), 3),
    mu_price  = round(c(mu1_pm[2], mu2_pm[2]), 3),
    mu_disp   = round(c(mu1_pm[3], mu2_pm[3]), 3)
)
print(comp_summ)
write.csv(comp_summ, 'cheese_hlm_mixture_K2_components.csv', row.names=FALSE)

Saved: cheese_hlm_mixture_K2.csv
  component weight mu_int mu_price mu_disp
1         1  0.683  9.765   -2.056   1.144
2         2  0.317  8.326   -1.876   1.787

---

Conclusion — harmonised across engines

Cheese model comparison (from-scratch marginal-over-β WAIC; lower = better):

K	WAIC	Δ vs K=1
1 (single normal)	−838.9	—
2	−843.3	−4.4 (mildly better)
3	−840.0	−1.1

A 2-component mixture is mildly preferred by WAIC (Δ≈4) over a single normal — the heterogeneity is slightly non-normal — but K=3 does not improve on K=2, and no store commits to a component (every store's maximum posterior membership clusters near 0.5 — none above ~0.65). The components overlap heavily, and the three engines find different modes: bayesm a base-volume split (μ-int 8.3/9.8, weights 0.32/0.68), the from-scratch sampler ≈0.5/0.5 (label-switching), and PyMC a display split (≈10/10) — with NUTS failing to converge (r̂≈1.5) on this multimodal target.

Unified reading. The cheese coefficient distribution is continuous and only weakly non-normal, not a set of discrete segments. The mixture earns its keep as a mild flexible-density refinement (the same role as in the margarine / Rossi et al. result), not as a segmentation device. The robust, cross-engine-consistent outputs are the single-normal hierarchy (population means, shrinkage) and the Δ covariate effect; the per-component means/weights are mode-dependent and should not be over-interpreted.