#######################################################################################
#R Code for Hierarchical Linear Model with random coefficients
#######################################################################################
# --- packages ---
library(readxl)
library(dplyr)
library(bayesm)
library(writexl)

# --- read & clean ---
# load
df <- read_excel("Car Conjoint Data.xlsx")
df$person  <- as.integer(df$Person)
df$rating <- as.numeric(df$Rating)
df$nissan <- as.numeric(df$Nissan)
df$toyota <- as.numeric(df$Toyota)
df$v4_170 <- as.numeric(df$V4_170HP)
df$v6_200 <- as.numeric(df$V6_200HP)
df$P27K   <- as.numeric(df$P27K)
df$P29K   <- as.numeric(df$P29K)

# --- build regdata: random intercept + random slope (ln_price) per store ---
store_ids <- df %>% distinct(person) %>% arrange(person) %>% pull(person)

regdata <- lapply(split(df, df$person), function(dfi) {
  X <- cbind(Intercept = 1, nissan = dfi$nissan, toyota=dfi$toyota, v4_170=dfi$v4_170, v6_200=dfi$v6_200, P27K=dfi$P27K, P29K=dfi$P29K)
  y <- dfi$rating
  list(X = X, y = y)
})

# --- fit model (progress every 1000 iters) ---

set.seed(123)
out <- bayesm::rhierLinearModel(
  Data = list(regdata = regdata, Z = NULL),
  Mcmc = list(R = 20000L, keep = 1L, nprint = 1000L)
)

# bayesm::rhierLinearModel storage (per the function’s source):
# betadraw   : [N x K x S]
# Deltadraw  : [S x nz x K]   (with Z=NULL => nz = 1)
# Vbetadraw  : [S x (K*K)]    (each row is vec(Sigma) in column-major order)

TOTAL <- 20000L
KEEP  <- 5000L
idx   <- (TOTAL - KEEP + 1):TOTAL  # last 5000 draws

###############################################
# Posterior mean of Bbar (the mean)
###############################################
Delta_means <- colMeans(out$Deltadraw[idx, , drop = FALSE])
# turn into a data frame: column vector with coefficient names
Delta_means_df <- data.frame(
  Coefficient = colnames(regdata[[1]]$X),
  Mean = as.numeric(Delta_means),
  row.names = NULL
)
Delta_means_df   # prints nicely

###############################################
# Posterior mean of V (the covariance matrix)
###############################################
# Average across rows
Vbeta_mean_vec <- colMeans(out$Vbetadraw[idx, , drop = FALSE])

# Reshape back to 7x7 covariance matrix
K <- 7
Vbeta_mean <- matrix(Vbeta_mean_vec, nrow = K, ncol = K, byrow = FALSE)

# Label rows and cols
colnames(Vbeta_mean) <- colnames(regdata[[1]]$X)
rownames(Vbeta_mean) <- colnames(regdata[[1]]$X)

Vbeta_mean_df <- data.frame(
  Coefficient = rownames(Vbeta_mean),
  as.data.frame(Vbeta_mean),
  row.names = NULL,
  check.names = FALSE
)

Vbeta_mean_df
###############################################
# Posterior mean of Bi (individual beta)
###############################################
beta_means <- apply(out$betadraw[, , idx], c(1, 2), mean)

# beta_means is N x K: rows = persons, cols = coefficients
N <- dim(beta_means)[1]
K <- dim(beta_means)[2]

# set names correctly
colnames(beta_means) <- colnames(regdata[[1]]$X)              # K names
if (exists("store_ids") && length(store_ids) == N) {
  rownames(beta_means) <- paste0("Person_", store_ids)         # N names
} else {
  rownames(beta_means) <- paste0("Person_", seq_len(N))        # fallback
}

# make a view with coefficient names in the first column and each person as a column
beta_means_df <- data.frame(
  Coefficient = colnames(beta_means),
  t(beta_means),                    # transpose so persons become columns
  row.names = NULL,
  check.names = FALSE
)

beta_means_df  # prints nicely

# Write all three to one workbook
write_xlsx(
  list(
    Delta_Means = Delta_means_df,
    Vbeta_Mean  = Vbeta_mean_df,
    Beta_Means  = beta_means_df
  ),
  path = "conjoint_posterior_summaries.xlsx"
)