Bayesian Multivariate Regression — R bayesm companion

Companion to mvr_python.ipynb. Same tuna demand system and the same four graphs, here via bayesm::rmultireg — the reference conjugate Normal–Inverse-Wishart sampler ($Y=XB+E$, $E_t\sim N_M(0,\Sigma)$), drawing i.i.d. posterior samples (no MCMC). Same diffuse prior ($\bar B=0$, $A=0.01 I$, $\nu=M+3$, $V=(M+2)I$).

suppressMessages({library(bayesm); library(ggplot2)})
Y <- as.matrix(read.csv('tuna_lmove.csv')); X <- as.matrix(read.csv('tuna_X.csv'))
N <- nrow(Y); M <- ncol(Y); K <- ncol(X); brands <- c('StarKist','ChickSea','BB Solid','BB Chunk','Geisha','BB Large','HH Chunk')
Bbar <- matrix(0, K, M); A <- 0.01*diag(K); nu <- M+3L; V <- (M+2)*diag(M)
Rdraws <- 10000L; set.seed(42)
Bd <- array(0, c(Rdraws, K, M)); Cd <- array(0, c(Rdraws, M, M))
for (r in 1:Rdraws) {
  dr <- rmultireg(Y, X, Bbar, A, nu, V)
  Bd[r,,] <- dr$B
  S <- dr$Sigma; d <- sqrt(diag(S)); Cd[r,,] <- S / outer(d, d)
}
Bmean <- apply(Bd, c(2,3), mean); Cmean <- apply(Cd, c(2,3), mean)
Ecross <- t(Bmean[2:(M+1), ])                       # rows 2..8 = lp1..lp7 ; (sales eq, price brand)
cat('Tuna: N=', N, ' M=', M, ' K=', K, '\n', sep='')
cat('own-price elasticities (diagonal):', round(diag(Ecross), 2), '\n')

Tuna: N=338 M=7 K=15

own-price elasticities (diagonal): -4.29 -4.66 -5.37 -4.87 -4.02 3.09 -2.36 

options(repr.plot.width=6.5, repr.plot.height=5)
g <- expand.grid(rowb=factor(brands, levels=brands), colb=factor(brands, levels=brands)); g$val <- as.vector(Cmean)
ggplot(g, aes(colb, rowb, fill=val)) + geom_tile() + geom_text(aes(label=sprintf('%.2f', val)), size=2.6) +
  scale_fill_gradient2(low='blue', mid='white', high='red', midpoint=0, limits=c(-1,1), name='corr') +
  scale_y_discrete(limits=rev) + coord_equal() +
  labs(title='Sigma error correlation across brands (bayesm)', x=NULL, y=NULL) + theme_minimal(base_size=12) + theme(axis.text.x=element_text(angle=45, hjust=1))

The posterior-mean error correlation across the seven brands from bayesm::rmultireg. The strong off-diagonal block linking the Bumble Bee variants stands out, and the off-diagonals generally sit well away from zero — which is precisely why the seven demand equations are estimated jointly rather than one at a time. The same pattern the from-scratch conjugate sampler and PyMC produce.

options(repr.plot.width=6.8, repr.plot.height=5)
g2 <- expand.grid(sales=factor(brands, levels=brands), price=factor(brands, levels=brands)); g2$val <- as.vector(Ecross)
vmax <- max(abs(Ecross))
ggplot(g2, aes(price, sales, fill=val)) + geom_tile() + geom_text(aes(label=sprintf('%.1f', val)), size=2.6) +
  scale_fill_gradient2(low='blue', mid='white', high='red', midpoint=0, limits=c(-vmax, vmax), name='elast.') +
  scale_y_discrete(limits=rev) + coord_equal() +
  labs(title='Price-elasticity matrix (bayesm)\ndiagonal=own-price, off-diagonal=cross-price', x='price brand', y='sales brand') +
  theme_minimal(base_size=12) + theme(axis.text.x=element_text(angle=45, hjust=1))

The price-elasticity matrix: the diagonal is own-price (strongly negative — a brand's own price cuts its own sales) and the off-diagonals are cross-price (largely positive ⇒ the brands are substitutes, so a rival's price rise lifts your sales). Structurally identical to the Python conjugate and PyMC fits.

options(repr.plot.width=9, repr.plot.height=4)
own <- sapply(1:M, function(j) Bd[, j+1, j])        # own-price draws per brand
fo <- data.frame(brand=factor(brands, levels=rev(brands)),
                 mean=colMeans(own), lo=apply(own,2,quantile,.025), hi=apply(own,2,quantile,.975))
p1 <- ggplot(fo, aes(mean, brand)) + geom_pointrange(aes(xmin=lo, xmax=hi), color='firebrick') +
  geom_vline(xintercept=0, linetype='dotted') + geom_vline(xintercept=-1, linetype='dashed', color='steelblue') +
  labs(title='Own-price elasticities (95% CrI)', x='own-price elasticity', y=NULL) + theme_minimal(base_size=12)
# strongest off-diagonal correlation
off <- Cmean; diag(off) <- 0; idx <- which(abs(off)==max(abs(off)), arr.ind=TRUE)[1,]; i_ <- idx[1]; j_ <- idx[2]
cd <- Cd[, i_, j_]
p2 <- ggplot(data.frame(cd=cd), aes(cd)) + geom_histogram(bins=50, fill='steelblue', color='white') +
  geom_vline(xintercept=mean(cd), color='firebrick', linewidth=1) + geom_vline(xintercept=0, linetype='dotted') +
  labs(title=sprintf('Strongest cross-equation correlation: corr(%s,%s)\n95%% CrI [%.2f, %.2f]', brands[i_], brands[j_], quantile(cd,.025), quantile(cd,.975)),
       x='error correlation', y=NULL) + theme_minimal(base_size=12)
suppressMessages(print(p1)); suppressMessages(print(p2))

Top — own-price elasticities with 95% credible intervals. Six of the seven brands sit far below the unit-elastic line (−4.3, −4.7, −5.4, −4.9, −4.0, −2.4): canned-tuna buyers switch readily. The exception is BB Large at +3.1 — an economically impossible positive own-price elasticity, and the same collinearity artefact the other two engines produce (it is entangled with BB Solid, whose error correlation with it is the strongest in the system). Bottom — the posterior of that strongest cross-equation correlation; its 95% CrI excludes zero. Three engines (conjugate NumPy · PyMC/NUTS · bayesm) land on the same numbers, artefact included.

Results

• bayesm::rmultireg reproduces the from-scratch and PyMC results: the same $\Sigma$ correlation pattern (strong $B3$–$B6$ link), the same elasticity matrix (negative own-price for 6 of 7 brands, the $B6$ collinearity artifact), and the same strongest cross-equation correlation with a 95% CrI excluding 0.
• Three engines agree (conjugate NumPy · PyMC/NUTS · R bayesm), confirming the implementation and the economics.