Scotch Whisky — Multivariate Probit via bayesm (rmvpGibbs)¶
Reference: Edwards, Y. D. & Allenby, G. M. (2003). Multivariate Analysis of Multiple Response Data. Journal of Marketing Research, 40(3), 321–334.
Package: Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005).
Bayesian Statistics and Marketing. Wiley. R implementation: bayesm::rmvpGibbs.
This notebook replicates the Scotch whisky analysis using the bayesm R implementation
and compares results against the hand-coded Python Gibbs sampler in Scotch_gibbs.ipynb.
Key difference in data format: rmvpGibbs expects data stacked observation-by-observation:
y: vector of length $N \times J = 46{,}599$X: matrix of size $N \times J \times K$ rows, brand indicators as columns
1. Setup¶
In [50]:
# Point R at the user library where bayesm and readxl were installed
userLib <- file.path(Sys.getenv("USERPROFILE"), "R", "win-library", "4.6")
.libPaths(c(userLib, .libPaths()))
# Install readxl if not present
if (!requireNamespace("readxl", quietly = TRUE))
install.packages("readxl", lib = userLib, repos = "https://cloud.r-project.org")
library(bayesm)
library(readxl)
cat("bayesm version :", as.character(packageVersion("bayesm")), "\n")
cat("readxl version:", as.character(packageVersion("readxl")), "\n")
bayesm version : 3.1.7 readxl version: 1.5.0
In [51]:
# Load the Y matrix exported from Python (ensures identical brands/ordering)
Y_df <- read.csv("scotch_Y.csv", check.names = FALSE)
BRANDS <- colnames(Y_df)
Y <- as.matrix(Y_df)
mode(Y) <- "integer"
N <- nrow(Y)
J <- ncol(Y)
cat(sprintf("N = %d households, J = %d brands\n", N, J))
cat(sprintf("Avg brands / HH: %.2f\n", mean(rowSums(Y))))
rates <- colMeans(Y)
rate_df <- data.frame(brand = BRANDS, rate = round(rates, 3))
rate_df <- rate_df[order(-rate_df$rate), ]
print(rate_df, row.names = FALSE)
N = 2219 households, J = 21 brands
Avg brands / HH: 2.29
brand rate
Chivas.Regal 0.363
Dewar.s.White.Label 0.233
Johnnie.Walker.Black.Label 0.226
J...B 0.206
Johnnie.Walker.Red.Label 0.191
Other.Brands 0.187
Glenlivet 0.160
Cutty.Sark 0.153
Glenfiddich 0.151
Pinch..Haig. 0.053
Clan.MacGregor 0.046
Ballantine 0.045
Macallan 0.043
Passport 0.037
Black...White 0.037
Scoresby.Rare 0.036
Grants 0.033
Ushers 0.030
White.Horse 0.028
Knockando 0.021
the.Singleton 0.014
In [52]:
barplot(sort(rates, decreasing = TRUE),
names.arg = abbreviate(names(sort(rates, decreasing = TRUE)), minlength = 6),
col = "steelblue", las = 2, cex.names = 0.65,
ylab = "Purchase rate",
main = sprintf("Scotch whisky brand purchase rates (N=%d)", N))
3. Set up model¶
rmvpGibbs expects stacked data: for each of the $N$ subjects, the $J$ brand
outcomes are listed consecutively. Brand-specific intercepts are implemented
via $J$ indicator columns in X (one dummy per brand).
In [53]:
# Stack y: brand outcomes within each subject
y_vec <- as.vector(t(Y)) # length N*J, brand varies fastest within subject
# Stacked brand-indicator design matrix (N*J x J)
X_stack <- kronecker(matrix(1, nrow = N, ncol = 1), diag(J))
cat("y_vec length :", length(y_vec), " (should be", N * J, ")\n")
cat("X_stack dim :", nrow(X_stack), "x", ncol(X_stack), "\n")
cat("table(y_vec) :\n"); print(table(y_vec))
# ── Prior ──────────────────────────────────────────────────────────────────
# beta ~ N(0, 100*I)
# Sigma ~ IW(nu, V): need nu > J+1 for proper dist with finite mean.
# E[sigma_jj] = V_jj / (nu - J - 1)
# With V = (nu-J-1)*I_J => E[sigma_jj] = 1 (matches correlation-matrix scale)
betabar <- rep(0, J)
A <- (1/100) * diag(J) # prior precision on beta
nu <- J + 3L # = 24; gives finite IW mean
V <- (nu - J - 1) * diag(J) # = 2*I_J => E[sigma_jj] = 1
Data <- list(p = J, y = y_vec, X = X_stack)
Prior <- list(betabar = betabar, A = A, nu = nu, V = V)
cat(sprintf("\nPrior: nu=%d | V=%.0f*I_J | E[sigma_jj]=%.1f | beta_var=100\n",
nu, nu - J - 1, (nu - J - 1) / (nu - J - 1)))
y_vec length : 46599 (should be 46599 )
X_stack dim : 46599 x 21
table(y_vec) :
y_vec
0 1
41514 5085
Prior: nu=24 | V=2*I_J | E[sigma_jj]=1.0 | beta_var=100
4. Run rmvpGibbs¶
In [54]:
NDRAWS <- 30000L
BURN <- 3000L
KEEP <- 1L # store every draw
Mcmc <- list(R = NDRAWS, keep = KEEP, nprint = 2000L)
cat(sprintf("Running rmvpGibbs: %d draws, %d burn-in ...\n", NDRAWS, BURN))
t0 <- proc.time()["elapsed"]
out <- rmvpGibbs(Data = Data, Prior = Prior, Mcmc = Mcmc)
cat(sprintf("Done in %.1f seconds.\n", proc.time()["elapsed"] - t0))
Running rmvpGibbs: 30000 draws, 3000 burn-in ...
Table of y values
y
0 1
41514 5085
Starting Gibbs Sampler for MVP
2219 obs of 21 binary indicators; 21 indep vars (including intercepts)
Prior Parms:
betabar
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00
[10,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
[11,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00
[12,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00
[13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
[14,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[15,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[16,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[17,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[18,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[10,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[11,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[12,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[14,] 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[15,] 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00
[16,] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00
[17,] 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00
[18,] 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
[19,] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00
[20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00
[21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
nu
[1] 24
V
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 2 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 2 0 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 2 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 2 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 2 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 2 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 2 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 2 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0 2 0 0 0 0
[10,] 0 0 0 0 0 0 0 0 0 2 0 0 0
[11,] 0 0 0 0 0 0 0 0 0 0 2 0 0
[12,] 0 0 0 0 0 0 0 0 0 0 0 2 0
[13,] 0 0 0 0 0 0 0 0 0 0 0 0 2
[14,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[15,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[16,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[17,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[18,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[19,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[20,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[21,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 0
[11,] 0 0 0 0 0 0 0 0
[12,] 0 0 0 0 0 0 0 0
[13,] 0 0 0 0 0 0 0 0
[14,] 2 0 0 0 0 0 0 0
[15,] 0 2 0 0 0 0 0 0
[16,] 0 0 2 0 0 0 0 0
[17,] 0 0 0 2 0 0 0 0
[18,] 0 0 0 0 2 0 0 0
[19,] 0 0 0 0 0 2 0 0
[20,] 0 0 0 0 0 0 2 0
[21,] 0 0 0 0 0 0 0 2
MCMC Parms:
30000 reps; keeping every 1 th draw nprint= 2000
initial beta= 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
initial sigma=
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 1 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0 1 0 0 0 0
[10,] 0 0 0 0 0 0 0 0 0 1 0 0 0
[11,] 0 0 0 0 0 0 0 0 0 0 1 0 0
[12,] 0 0 0 0 0 0 0 0 0 0 0 1 0
[13,] 0 0 0 0 0 0 0 0 0 0 0 0 1
[14,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[15,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[16,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[17,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[18,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[19,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[20,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[21,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 0
[11,] 0 0 0 0 0 0 0 0
[12,] 0 0 0 0 0 0 0 0
[13,] 0 0 0 0 0 0 0 0
[14,] 1 0 0 0 0 0 0 0
[15,] 0 1 0 0 0 0 0 0
[16,] 0 0 1 0 0 0 0 0
[17,] 0 0 0 1 0 0 0 0
[18,] 0 0 0 0 1 0 0 0
[19,] 0 0 0 0 0 1 0 0
[20,] 0 0 0 0 0 0 1 0
[21,] 0 0 0 0 0 0 0 1
MCMC Iteration (est time to end - min)
2000 (10.7)
4000 (11.1)
6000 (10.5)
8000 (9.3)
10000 (8.3)
12000 (7.5)
14000 (6.6)
16000 (5.7)
18000 (4.8)
20000 (4.0)
22000 (3.2)
24000 (2.4)
26000 (1.6)
28000 (0.8)
30000 (0.0)
Total Time Elapsed: 11.88
Done in 712.6 seconds.
5. Extract posterior draws¶
In [55]:
keep_idx <- seq(BURN + 1L, NDRAWS, by = KEEP)
# Raw (unidentified) draws
beta_draws <- out$betadraw[keep_idx, ]
sigma_draws <- out$sigmadraw[keep_idx, ]
n_keep <- nrow(beta_draws)
# ── Normalize each draw to the identified parameterization ─────────────────
# rmvpGibbs does NOT fix sigma_jj = 1, so raw beta and Sigma are unidentified.
# The identified quantities are:
# beta*_j = beta_j / sqrt(sigma_jj)
# Sigma*_jk = sigma_jk / sqrt(sigma_jj * sigma_kk) (correlation matrix)
beta_norm_draws <- matrix(NA_real_, nrow = n_keep, ncol = J)
sigma_norm_draws <- matrix(NA_real_, nrow = n_keep, ncol = J * J)
for (r in seq_len(n_keep)) {
S <- matrix(sigma_draws[r, ], J, J)
sd <- sqrt(diag(S))
beta_norm_draws[r, ] <- beta_draws[r, ] / sd
sigma_norm_draws[r, ] <- as.vector(S / outer(sd, sd)) # correlation matrix
}
# Posterior means of identified parameters
B_mean_R <- colMeans(beta_norm_draws)
names(B_mean_R) <- BRANDS
R_draws_R <- array(sigma_norm_draws, dim = c(n_keep, J, J))
R_mean_R <- apply(R_draws_R, c(2, 3), mean)
dimnames(R_mean_R) <- list(BRANDS, BRANDS)
cat(sprintf("Stored draws : %d\n", n_keep))
cat(sprintf("Sigma* diagonal: min=%.4f max=%.4f (should be ~1)\n",
min(diag(R_mean_R)), max(diag(R_mean_R))))
# Model-implied rates — now pnorm(beta*) should match observed
impl_rates <- pnorm(B_mean_R)
cat("\nModel-implied vs observed rates:\n")
print(data.frame(brand = BRANDS,
observed = round(rates, 3),
implied = round(impl_rates, 3),
beta_star = round(B_mean_R, 3)),
row.names = FALSE)
Stored draws : 27000
Sigma* diagonal: min=1.0000 max=1.0000 (should be ~1)
Model-implied vs observed rates:
brand observed implied beta_star
Chivas.Regal 0.363 0.363 -0.351
Dewar.s.White.Label 0.233 0.234 -0.726
Johnnie.Walker.Black.Label 0.226 0.227 -0.749
J...B 0.206 0.207 -0.816
Johnnie.Walker.Red.Label 0.191 0.193 -0.867
Other.Brands 0.187 0.187 -0.890
Glenlivet 0.160 0.160 -0.994
Cutty.Sark 0.153 0.154 -1.019
Glenfiddich 0.151 0.151 -1.032
Pinch..Haig. 0.053 0.052 -1.622
Clan.MacGregor 0.046 0.047 -1.673
Ballantine 0.045 0.045 -1.695
Macallan 0.043 0.044 -1.709
Passport 0.037 0.038 -1.780
Black...White 0.037 0.037 -1.793
Scoresby.Rare 0.036 0.036 -1.798
Grants 0.033 0.034 -1.830
Ushers 0.030 0.031 -1.870
White.Horse 0.028 0.027 -1.919
Knockando 0.021 0.021 -2.034
the.Singleton 0.014 0.013 -2.228
6. Brand intercepts vs observed purchase rates¶
In [56]:
impl_rates <- pnorm(B_mean_R)
comp <- data.frame(
brand = BRANDS,
observed = round(rates, 3),
model_implied = round(impl_rates, 3),
beta = round(B_mean_R, 3)
)
comp <- comp[order(-comp$observed), ]
print(comp, row.names = FALSE)
# R² of observed vs model-implied
r2 <- cor(comp$observed, comp$model_implied)^2
cat(sprintf("\nR² (observed vs model-implied): %.6f\n", r2))
brand observed model_implied beta
Chivas.Regal 0.363 0.363 -0.351
Dewar.s.White.Label 0.233 0.234 -0.726
Johnnie.Walker.Black.Label 0.226 0.227 -0.749
J...B 0.206 0.207 -0.816
Johnnie.Walker.Red.Label 0.191 0.193 -0.867
Other.Brands 0.187 0.187 -0.890
Glenlivet 0.160 0.160 -0.994
Cutty.Sark 0.153 0.154 -1.019
Glenfiddich 0.151 0.151 -1.032
Pinch..Haig. 0.053 0.052 -1.622
Clan.MacGregor 0.046 0.047 -1.673
Ballantine 0.045 0.045 -1.695
Macallan 0.043 0.044 -1.709
Passport 0.037 0.038 -1.780
Black...White 0.037 0.037 -1.793
Scoresby.Rare 0.036 0.036 -1.798
Grants 0.033 0.034 -1.830
Ushers 0.030 0.031 -1.870
White.Horse 0.028 0.027 -1.919
Knockando 0.021 0.021 -2.034
the.Singleton 0.014 0.013 -2.228
R² (observed vs model-implied): 0.999937
7. Comparison with Python Gibbs sampler¶
In [57]:
# Load Python posterior means (exported from Scotch_gibbs.ipynb)
B_py <- read.csv("scotch_B_mean_py.csv")
R_py <- as.matrix(read.csv("scotch_R_mean_py.csv", row.names = 1))
# ── Beta comparison ───────────────────────────────────────────────────────────
beta_comp <- data.frame(
brand = BRANDS,
beta_R = round(B_mean_R, 4),
beta_Py = round(B_py$beta_py, 4),
diff = round(B_mean_R - B_py$beta_py, 4)
)
beta_comp <- beta_comp[order(-rates), ]
cat("Brand intercept comparison (R bayesm vs Python Gibbs):\n")
print(beta_comp, row.names = FALSE)
cat(sprintf("\nMax |diff| in beta: %.4f\n", max(abs(beta_comp$diff))))
Brand intercept comparison (R bayesm vs Python Gibbs):
brand beta_R beta_Py diff
Chivas.Regal -0.3507 -0.3503 -0.0005
Dewar.s.White.Label -0.7264 -0.7265 0.0001
Johnnie.Walker.Black.Label -0.7493 -0.7492 -0.0001
J...B -0.8162 -0.8162 0.0000
Johnnie.Walker.Red.Label -0.8666 -0.8669 0.0003
Other.Brands -0.8900 -0.8893 -0.0008
Glenlivet -0.9945 -0.9944 -0.0001
Cutty.Sark -1.0194 -1.0196 0.0002
Glenfiddich -1.0318 -1.0308 -0.0010
Pinch..Haig. -1.6221 -1.6221 0.0000
Clan.MacGregor -1.6729 -1.6732 0.0003
Ballantine -1.6952 -1.6955 0.0004
Macallan -1.7090 -1.7116 0.0026
Passport -1.7803 -1.7788 -0.0016
Black...White -1.7928 -1.7912 -0.0015
Scoresby.Rare -1.7975 -1.7961 -0.0014
Grants -1.8302 -1.8279 -0.0023
Ushers -1.8696 -1.8684 -0.0012
White.Horse -1.9192 -1.9178 -0.0015
Knockando -2.0338 -2.0333 -0.0005
the.Singleton -2.2284 -2.2214 -0.0070
Max |diff| in beta: 0.0070
In [58]:
# Lower-triangle correlation comparison
idx <- which(lower.tri(R_mean_R), arr.ind = TRUE)
b1 <- BRANDS[idx[,1]]
b2 <- BRANDS[idx[,2]]
corr_R <- R_mean_R[idx]
corr_Py <- R_py[idx]
corr_comp <- data.frame(
brand_1 = b1,
brand_2 = b2,
R_bayesm = round(corr_R, 3),
R_python = round(corr_Py, 3),
diff = round(corr_R - corr_Py, 3)
)
corr_comp <- corr_comp[order(-abs(corr_comp$R_bayesm)), ]
cat("Top 20 correlations — bayesm vs Python:\n")
print(head(corr_comp, 20), row.names = FALSE)
cat(sprintf("\nMax |diff| in correlations : %.4f\n", max(abs(corr_comp$diff))))
cat(sprintf("Mean |diff| in correlations: %.4f\n", mean(abs(corr_comp$diff))))
Top 20 correlations — bayesm vs Python:
brand_1 brand_2 R_bayesm R_python diff
Knockando Macallan 0.683 0.684 -0.001
Glenfiddich Glenlivet 0.598 0.598 0.001
the.Singleton Ushers 0.523 0.536 -0.013
the.Singleton Scoresby.Rare 0.509 0.508 0.001
the.Singleton Passport 0.499 0.496 0.003
Ushers Scoresby.Rare 0.497 0.493 0.004
the.Singleton Macallan 0.491 0.496 -0.005
Knockando Glenfiddich 0.489 0.487 0.002
Knockando Glenlivet 0.477 0.474 0.003
Grants Passport 0.472 0.469 0.003
the.Singleton White.Horse 0.450 0.449 0.000
Ushers Grants 0.449 0.447 0.001
Macallan Glenlivet 0.441 0.442 0.000
Black...White Ballantine 0.436 0.431 0.004
Ushers Passport 0.433 0.431 0.002
White.Horse Black...White 0.427 0.427 0.001
White.Horse Ballantine 0.426 0.428 -0.002
the.Singleton Knockando 0.419 0.415 0.005
White.Horse Passport 0.407 0.406 0.002
Passport Clan.MacGregor 0.404 0.409 -0.004
Max |diff| in correlations : 0.0150
Mean |diff| in correlations: 0.0028
8. Save results¶
In [59]:
write.csv(data.frame(brand=BRANDS, beta=B_mean_R), "scotch_B_mean_R.csv", row.names=FALSE)
write.csv(as.data.frame(R_mean_R), "scotch_R_mean_R.csv")
cat("Saved: scotch_B_mean_R.csv, scotch_R_mean_R.csv\n")
Saved: scotch_B_mean_R.csv, scotch_R_mean_R.csv