---
title: 'Conjoint Analysis: What is the best laptop design?'
author: "Savaş Dayanık"
date: "3/9/2021"
output: 
  html_document:
    toc: yes
---

```{r setup, include=FALSE}
library(magrittr)
library(tidyverse)
library(readxl)
knitr::opts_chunk$set(echo = TRUE)
```

# Preference of one subject

```{r}
train_hande_df <- read_csv("laptops-hande.csv") %>%
  mutate_if(is.character, factor)
test_hande_df <- read_csv("laptops-hepsiburada-hande.csv", comment = "#") %>% 
  mutate_if(is.character, factor)

summary(train_hande_df)

train_hande_df %>% 
  arrange(RANK)
```

```{r}
options("contrasts")
options(contrasts = c("contr.sum", "contr.poly"))
```
# Estimate Hande's part worth and importance scores

```{r}
coef_df <- lm(-RANK ~ ., select(train_hande_df, - Alternative)) %>% 
  summary() %>% 
  coef() %>% 
  as_tibble(rownames = "Part") %>% 
  select(Part, Estimate) %>% 
  filter(Part != "(Intercept)") %>% 
  rename(Worth = Estimate)

parts_df <- coef_df %>% 
  mutate(ID = as.numeric(str_extract(Part, "\\d+")),
         Part = str_extract(Part, "\\D*")) %>% 
  nest(data=-Part) %>% 
  mutate(data2 = map(data, ~rbind(.x, tibble(Worth = -sum(.x$Worth),
                                              ID = max(.x$ID + 1))))) %>%
  select(-data) %>% 
  unnest(data2) %>% 
  mutate(PartLevel = map2_chr(Part, ID, ~levels(train_hande_df[[.x]])[.y])) %>% 
  group_by(Part) %>% 
  mutate(Importance = diff(range(Worth)))
```

# Plot Hande's part worth and importance scores

```{r, fig.asp=1}
parts_df %>% 
  ggplot(aes(PartLevel, Worth)) +
  geom_point() +
  geom_linerange(aes(ymin=0, ymax = Worth)) +
  coord_flip() +
  facet_grid(reorder(Part, -Importance) ~ ., space="free", scale="free") +
  theme(strip.text.y = element_text(angle=0))
```

# Estimate market participant part worths and and importance scores

Now repeat the same analysis for all subjects who filled up the questionnaire.

```{r}
# d2 <- read_csv("conjoint/all/laptops train survey results-all.csv") %>%
train_df <- read_excel("laptops train survey results-all-coded.xlsx") %>%  
  mutate_if(is.character, factor) %>% 
  mutate(id = as.character(id),
         Alternative = factor(Alternative, level = levels(Alternative)[order(parse_number(levels(Alternative)))]),
         Price = factor(Price, level = levels(Price)[order(parse_number(levels(Price)))]))

train_df %>% summary()

test_df <- read_excel("laptops test survey results-all-coded.xlsx") %>% 
    mutate_if(is.character, factor)

test_df %>% summary()
```


```{r, warning=FALSE}
market_df <- train_df %>% 
  nest(train_data = -id) %>% 
  mutate(
    # learn individual preferences of each respondent / subject
    model = map2(train_data, id,  ~{
      # cat ("ID=", .y, "\n")
      d <- .x
      lm(-RANK ~ ., select(d, - Alternative))
      }),
    cat("model     : ranks are learned with multiple regression.\n"),
    
    # compile the part worths  
    pworth = pmap(list(model, train_data, id),  ~{
      # cat ("ID=", ..3, "\n")
      d <- ..2
      linmod <- ..1
      
      coef_df <- linmod %>% 
      summary() %>% 
      coef() %>% 
      as_tibble(rownames = "Part") %>% 
      select(Part, Estimate) %>% 
      filter(Part != "(Intercept)") %>% 
      rename(Worth = Estimate)
  
      parts_df <- coef_df %>%
        mutate(ID = as.numeric(str_extract(Part, "\\d+")),
               Part = str_extract(Part, "\\D*")) %>%
        nest(data=-Part) %>%
        mutate(data2 = map(data, ~rbind(.x, tibble(Worth = -sum(.x$Worth),
                                                    ID = max(.x$ID + 1))))) %>%
        select(-data) %>%
        unnest(data2) %>%
        mutate(PartLevel = map2_chr(Part, ID, ~levels(d[[.x]])[.y])) %>%
        group_by(Part) %>%
        mutate(Importance = diff(range(Worth)),
               # below we pad from left numerical Part levels in order
               # to maintain the numerical order of levels in partworth plots
               # stored in pwplot.
               PartLevel = if (any(str_detect(PartLevel, "\\d"))) {
                 str_pad(PartLevel, width = max(str_length(PartLevel))) 
               } else {
                   PartLevel
                 }
               ) %>% 
        ungroup()
      parts_df
      }),
    cat("pworth    : part worth scores are calculated.\n"),
    
    # extract importance scores
    importance = map2(pworth, id, ~{
      # cat ("ID=", .y, "\n")
      .x %>% 
        distinct(Part, Importance) %>% 
        mutate(Importance = 100*Importance/sum(Importance))
    }),
    cat("importance: part importance scores are calculated.\n"),
    
    # visualize the part worths  
    pwplot = map2(pworth, id,  ~{
      # cat ("ID=", .y, "\n")
      parts_df <- .x 

    parts_df %>% 
        ggplot(aes(PartLevel, Worth)) +
        geom_point() +
        geom_linerange(aes(ymin=0, ymax = Worth)) +
        coord_flip() +
        facet_grid(reorder(Part, -Importance) ~ ., space="free", scale="free") +
        theme(strip.text.y = element_text(angle=0)) +
        labs(title = .y)
      
  }),
  cat("pwplot    : part worth plots are ready.\n"),
  ) %>% 
  
  # append test data
  left_join(nest(test_df, test_data = -id), by = "id") %>% 
  mutate(
    # rank predictions on test
    test_pred = pmap(list(model, test_data, id), ~{
      # cat(..3, " ")
      pred <- ..1 %>% 
        predict(newdata=..2) %>% 
        `*` (-1) 
      rank <- pred %>% 
        rank(ties.method = "first")
      
      ..2 %>% 
        select(Alternative, RANK) %>% 
        mutate(PREDraw = pred, PREDrank = rank)
      }),
    cat("test_pred : rankings on test data are predicted.\n")
  )

  #   test_observed_pred = map2(test_data, test_pred, ~ {
  #     .x %>% 
  #       select(Alternative, RANK) %>% 
  #       mutate(PRED = .y)
  #   }),
  #   cat("test_opred: observed and predicted rankings are combined.\n"),
  # )

```

# Plot a few market participant part worths and and importance scores

Show a few part worth plots 

```{r, fig.asp=.8}
# for (i in 1:10) plot(market_df$pwplot[[i]])

market_df %>% 
  filter(str_detect(id, "y22")) %>% 
  pull(pwplot) %>% 
  walk(plot)

```

# Validation of market model

```{r}
confusion <- market_df %>% 
  select(id, test_pred) %>% 
  unnest(test_pred) %>% 
  xtabs(~ RANK + PREDrank, .)

confusion
confnames <- dimnames(confusion) %>% names()
heatmap(confusion, Rowv = NA, Colv = NA, revC = TRUE, xlab = confnames[2], ylab = confnames[1])

# confusion %>% 
#   as_tibble() %>% 
#   rename(Observed = RANK, Predicted = PRED) %>% 
#   ggplot(aes(Predicted, Observed)) +
#   geom_tile(aes(fill=n)) +
#   scale_fill_gradient(low = "yellow", high = "red") +
#   coord_fixed()
```


```{r}
market_df %>% 
  select(id, test_pred) %>% 
  unnest(test_pred) %>% 
  group_by(id) %>% 
  summarize(agreement = sum(RANK == PREDrank), .groups = "drop") %>% 
  arrange(desc(agreement)) %>% 
  count(agreement) %>% 
  ggplot(aes(agreement, n)) +
  geom_col()
```

# Market shares

```{r}
mshares <- market_df %>% 
  select(test_data) %>% 
  unnest(test_data) %>% 
  filter(RANK==1) %>% 
  count(Brand) %>% 
  column_to_rownames(var = "Brand") %>% 
  as.matrix %>% 
  drop() %>% 
  prop.table() %>% 
  `*`(100)

mshares %>% 
  pie(main="Market shares")
```

# Substitution rates

```{r, fig.show = "hold", out.width="50%"}
subst_df <- market_df %>% 
  select(id, test_data) %>% 
  unnest(test_data) %>% 
  select(id, Brand, RANK) %>% 
  filter(RANK<=2) %>% 
  pivot_wider(names_from = "RANK", values_from = "Brand") %>% 
  count(`1`,`2`) %>% 
  rename(Brand = `1`, Substitute = `2`) %>% 
  nest(Subs = -Brand)

oldpar <- par(mar = c(0,0,1,0))
  
for (i in seq(nrow(subst_df))){
  subst_df$Subs[[i]] %>% 
      column_to_rownames(var = "Substitute") %>% 
      as.matrix %>% 
      drop() %>% 
      pie(main=paste("Substitutes of", subst_df$Brand[i] ))
}
par(mar = oldpar$mar)  
```

# Market simulation

All possible Casper configurations

```{r}
config_df <- test_df %>% 
  select(Processor, Speed, Memory, HardDisk, Screen, Price) %>% 
  expand(Processor, Speed, Memory, HardDisk, Screen, Price) %>% # expand to all permutations of the values of listed variables
  mutate(Brand = "Casper") %>% # our company is Casper
  relocate(Brand) %>%  
  anti_join(test_df, by = c("Brand", "Processor", "Speed", "Memory", "HardDisk", "Screen", "Price")) %>% # remove the Casper products already in the market
  mutate(config = as.character(seq(n()))) %>% 
  relocate(config)
```

Predicted rank utilities of all configurations

```{r}
market_config_df <- market_df %>% 
  mutate(config_pred = map2(model, id, ~{
    pred <- predict(.x, config_df) %>% 
      `*`(-1)
    config_df %>% 
      select(config) %>% 
      mutate(PREDraw = pred)
  })) %>% 
  select(id, test_pred, config_pred) %>% 
  unnest(test_pred) %>% 
  select(-RANK, -PREDrank) %>% 
  pivot_wider(names_from = Alternative, values_from = "PREDraw") %>% 
  unnest(config_pred) %>% 
  rename(new = PREDraw)

market_config_df %>% 
  head(10)
```

New markets obtained by adding each config to the existing hepsiburada market

```{r}
new_markets_df <- market_config_df %>% 
  pivot_longer(cols = c(-id, -config), names_to = "Alternative", values_to = "PREDraw") %>% 
  group_by(id, config) %>% 
  mutate(PREDrank = rank(PREDraw, ties.method = "first")) %>% 
  ungroup() %>% 
  left_join(mutate_all(select(test_hande_df, Alternative, Brand), as.character), by = "Alternative") %>% 
  replace_na(list(Brand = "Casper")) %>% 
  nest(new_market = -config) 

new_markets_df  %>% 
  head(10)
```

```{r}
new_markets_df$new_market[[1]] %>% 
  filter(PREDrank == 1) %>% 
  count(Brand) %>% 
  mutate(Perc = 100*n/sum(n))
```

Market shares after each new configuration is introduced to the existing hepsiburada market

```{r, fig.asp=0.8}
new_markets_stat_df <- new_markets_df %>% 
  mutate(market_shares = map(new_market, ~{
    .x %>% 
      filter(PREDrank == 1) %>% 
      count(Brand) %>% 
      mutate(Perc = 100*n/sum(n))
  }),
  Casper_mshare = map_dbl(market_shares, ~{.x %>% filter(Brand == "Casper") %>% pull("Perc")})) %>% 
  left_join(config_df, by = "config") %>% 
  mutate_at(vars(Memory, Price), ~{parse_number(as.character(.), locale = locale(grouping_mark = "."))}) %>% 
  arrange(desc(Casper_mshare), Price)

new_markets_stat_df %>%
  select(-new_market, -market_shares) %>% 
  filter(Casper_mshare >= mshares["Casper"])

new_markets_stat_df %>% 
  filter(Casper_mshare >= mshares["Casper"]) %>% 
  select(-new_market, -market_shares) %>% 
  ggplot(aes(Price, Casper_mshare)) +
  geom_hline(yintercept = mshares["Casper"], size=2, col = "yellow") +
  geom_jitter(aes(col = Screen), width=200, height=0, alpha = 0.8, size=3) +
  scale_x_continuous(labels = scales::comma) +
  labs(title = paste("Current Casper market share is", round(mshares["Casper"]), "% (yellow horizontal line)"),
       y = "Predicted Casper market share (%)") +
  facet_grid(Memory~HardDisk, as.table = FALSE, labeller = label_both) +
  theme(strip.text.y = element_text(angle=0))
```

