#####################################################
# Conjoint Analysis                                 #
# Chapter 12 - Chestnut Ridge Smart Watch Example   #
#####################################################

# This script reproduces the full Chestnut Ridge smart watch conjoint analysis example in R:
#   Step 1. Conjoint study design (full and fractional factorial designs)
#   Step 2. Estimate the part-worths for each respondent
#   Step 3. Part-worths, attribute importance, willingness-to-pay, and market simulations
# Run it from top to bottom. Each step prints its result to the console or Plots pane.

#####################################################
# Step 1: Conjoint Study Design Setup               #
#####################################################

## Install Packages (if needed - only installs packages that are missing)
packages <- c("conjoint", "ggplot2")
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) install.packages(new_packages)

## Load Packages and Set Seed
library(conjoint)
library(ggplot2)
set.seed(1)

## Set up attributes and levels as a list
## The first level of each attribute is the baseline level
attrib.level <- list(brand = c("CR", "Apple", "Samsung", "FitBit"),
                     ship = c("$0", "$10", "$20"),
                     restock = c("0%", "5%", "10%", "15%"),
                     retdays = c("7 days", "14 days", "21 days"),
                     price = c("$150", "$200", "$250", "$300"))

## Create the full factorial design (576 product bundles)
experiment <- expand.grid(attrib.level)
nrow(experiment)

## Create the fractional factorial design (30 product profiles)
design <- caFactorialDesign(data = experiment, type = "fractional", cards = 30, seed = 1)
print(design)

## Check for correlation in fractional factorial design
print(cor(caEncodedDesign(design)))

## Export design for survey (the Python version of this example reads in this file)
write.csv(design, file.choose(new = TRUE), row.names = FALSE) ## Name the file conjoint_profiles.csv

#####################################################
# Step 2: Estimate the Part-Worths                  #
#####################################################

## Run the conjoint analysis study

## Read in the survey preference results
pref <- read.csv(file.choose()) ## Choose the file named conjoint_preferences.csv

## Set up the product profiles as factors
## Setting the levels makes the first level of each attribute the baseline case
## Base Case: Brand CR, Shipping $0, Restock 0%, Retdays 7 days, Price $150
profiles <- design
for (a in names(attrib.level)) {
  profiles[[a]] <- factor(profiles[[a]], levels = attrib.level[[a]])
}

## Estimate the part-worths for each respondent using OLS regression
## Each attribute enters the regression as dummy variables (one per non-baseline level)
part.worths <- NULL
for (i in 1:ncol(pref)) {
  fit <- lm(pref[, i] ~ brand + ship + restock + retdays + price, data = profiles)
  b <- coef(fit)
  ## Part-worths are 0 for each baseline level
  temp <- c(intercept = b[["(Intercept)"]],
            CR = 0, Apple = b[["brandApple"]], Samsung = b[["brandSamsung"]], FitBit = b[["brandFitBit"]],
            "$0" = 0, "$10" = b[["ship$10"]], "$20" = b[["ship$20"]],
            "0%" = 0, "5%" = b[["restock5%"]], "10%" = b[["restock10%"]], "15%" = b[["restock15%"]],
            "7 days" = 0, "14 days" = b[["retdays14 days"]], "21 days" = b[["retdays21 days"]],
            "$150" = 0, "$200" = b[["price$200"]], "$250" = b[["price$250"]], "$300" = b[["price$300"]])
  part.worths <- rbind(part.worths, temp)
}
part.worths <- data.frame(round(part.worths, 3), check.names = FALSE)
rownames(part.worths) <- colnames(pref)
print(part.worths)

#####################################################
# Step 3: Part-Worth Analysis, Willingness-to-Pay,  #
#         and Market Simulations                    #
#####################################################

## Function to create a horizontal bar chart with labels
## Bars are shown in the order they are listed
plot.bars <- function(values, title, labels = round(values, 3)) {
  df <- data.frame(name = factor(names(values), levels = rev(names(values))),
                   value = values, label = labels)
  ggplot(df, aes(x = value, y = name)) +
    geom_col(fill = "steelblue") +
    geom_text(aes(label = label, hjust = ifelse(value >= 0, -0.1, 1.1)), size = 3.5) +
    geom_vline(xintercept = 0) +
    scale_x_continuous(expand = expansion(mult = 0.2)) +
    labs(title = title, x = NULL, y = NULL) +
    theme_minimal()
}

## Average Attribute Part-Worths
## Exclude the baseline levels since their part-worths are always 0
baseline <- sapply(attrib.level, function(x) x[1])
avg.pw <- colMeans(part.worths)
avg.pw <- avg.pw[!(names(avg.pw) %in% baseline)]
print(round(avg.pw, 3))
plot.bars(avg.pw, "Average Attribute Part-Worths")

## Average Attribute Importance
## Importance is the range (max - min) of the part-worths within each attribute
importance <- NULL
for (a in names(attrib.level)) {
  pw.attrib <- part.worths[, attrib.level[[a]]]
  importance <- cbind(importance, apply(pw.attrib, 1, max) - apply(pw.attrib, 1, min))
}
colnames(importance) <- c("Brand", "Shipping", "Restock", "Return Days", "Price")
avg.importance <- colMeans(importance)
print(round(avg.importance, 3))
plot.bars(avg.importance, "Average Attribute Importance")

## Percentage Average Attribute Importance
pct.importance <- avg.importance / sum(avg.importance)
names(pct.importance) <- paste(names(pct.importance), "%")
print(round(pct.importance * 100, 2))
plot.bars(pct.importance, "Percentage Average Attribute Importance",
          labels = sprintf("%.2f%%", pct.importance * 100))

## Average Willingness to Pay for a Feature
## The change in price from $150 to $300 ($150) tells us how many dollars one util is worth
dollars.per.util <- 150 / (colMeans(part.worths)[["$150"]] - colMeans(part.worths)[["$300"]])
features <- c("Apple", "Samsung", "FitBit", "$10", "$20", "5%", "10%", "15%", "14 days", "21 days")
wtp <- colMeans(part.worths)[features] * dollars.per.util
names(wtp) <- c("WTP - Brand Apple", "WTP - Brand Samsung", "WTP - Brand FitBit",
                "WTP - Shipping $10", "WTP - Shipping $20",
                "WTP - Restock 5%", "WTP - Restock 10%", "WTP - Restock 15%",
                "WTP - Return Days 14", "WTP - Return Days 21")
print(round(wtp, 2))
plot.bars(wtp, "Average Willingness to Pay for a Feature",
          labels = ifelse(wtp < 0, sprintf("-$%.2f", abs(wtp)), sprintf("$%.2f", wtp)))

## Function to calculate each respondent's utility for a set of products
## Each row of products is one product profile
calc.util <- function(products) {
  util <- NULL
  for (p in 1:nrow(products)) {
    u <- part.worths$intercept
    for (a in names(attrib.level)) {
      u <- u + part.worths[[products[p, a]]]
    }
    util <- cbind(util, u)
  }
  colnames(util) <- products$product
  rownames(util) <- rownames(part.worths)
  util
}

## Function to estimate market shares using three choice rules
calc.shares <- function(util) {
  ## Maximum utility (first choice) rule: each respondent chooses the product with the highest utility
  first.choice <- t(apply(util, 1, function(u) as.numeric(u == max(u))))
  ## Share of preference rule: choice probability is proportional to utility
  ## Negative utilities are set to 0 so that no product receives a negative probability
  util.pos <- pmax(util, 0)
  share.pref <- util.pos / rowSums(util.pos)
  ## Logit share of preference rule: choice probability is proportional to exp(utility)
  logit.share <- exp(util) / rowSums(exp(util))
  data.frame(rule = rep(c("First Choice", "Share of Preference", "Logit Share"), each = ncol(util)),
             product = rep(colnames(util), 3),
             share = c(colMeans(first.choice), colMeans(share.pref), colMeans(logit.share)))
}

## Function to plot market shares by choice rule
plot.shares <- function(shares, title) {
  shares$rule <- factor(shares$rule, levels = c("First Choice", "Share of Preference", "Logit Share"))
  shares$product <- factor(shares$product, levels = rev(unique(shares$product)))
  ggplot(shares, aes(x = share, y = product)) +
    geom_col(fill = "steelblue") +
    geom_text(aes(label = sprintf("%.2f%%", share * 100)), hjust = -0.1, size = 3.5) +
    facet_wrap(~ rule, ncol = 1) +
    scale_x_continuous(labels = function(x) paste0(x * 100, "%"), expand = expansion(mult = c(0, 0.2))) +
    labs(title = title, x = NULL, y = NULL) +
    theme_minimal()
}

## Market Share - Current Products
current <- data.frame(product = c("Apple", "Samsung", "FitBit"),
                      brand = c("Apple", "Samsung", "FitBit"),
                      ship = c("$0", "$20", "$10"),
                      restock = c("15%", "0%", "10%"),
                      retdays = c("7 days", "14 days", "14 days"),
                      price = c("$200", "$300", "$250"))
util.current <- calc.util(current)
print(round(util.current, 3))
shares.current <- calc.shares(util.current)
print(shares.current)
plot.shares(shares.current, "Market Share - Current Products")

## Market Share - Proposed Products (add the CR smart watch)
proposed <- rbind(data.frame(product = "CR", brand = "CR", ship = "$10", restock = "10%",
                             retdays = "14 days", price = "$250"),
                  current)
util.proposed <- calc.util(proposed)
print(round(util.proposed, 3))
shares.proposed <- calc.shares(util.proposed)
print(shares.proposed)
plot.shares(shares.proposed, "Market Share - Proposed Products")