##########################################################
# Segmentation Analysis in R - Chestnut Ridge Example    #
# Cluster Analysis for Segmentation                      #
##########################################################

# This script reproduces the full Chestnut Ridge segmentation example in R:
#   Step 1. Visualize average order size by zip code (map)
#   Step 2. Hierarchical clustering and elbow plot to choose the number of segments
#   Step 3. K-means clustering with six segments
#   Step 4. Profile the segments on the bases (attitudes and behaviors)
#   Step 5. Profile the segments on the descriptors (demographics)
#   Step 6. Evaluate segment attractiveness (sales and profit)
# Run it from top to bottom. Each step prints its result to the console or Plots pane.

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

## Load Packages
library(dplyr)     # data manipulation
library(tidyr)     # reshaping summary tables
library(ggplot2)   # plots and maps
library(sf)        # working with map shapes
library(tigris)    # zip code and state boundaries from the US Census Bureau

## Set Seed
# Segmentation algorithms start at a random point. Setting the seed ensures
# your results are the same as the results in this example.
set.seed(1)

# Import Data
seg <- read.csv(file.choose()) ## Choose retail_segmentation.csv file

# Take a first look at the data
summary(seg)


#####################################################
# Step 1: Average Order Size by Zip Code (Map)      #
#####################################################

# Zip codes are stored as numbers, so zip codes that start with 0 (e.g., 08053)
# lost their leading zero. Add it back so every zip code has five digits.
seg$zip_code <- sprintf("%05d", seg$zip_code)

# Average order size for the customers in each zip code
zip_summary <- seg %>%
  group_by(zip_code) %>%
  summarise(avg_order_size = mean(avg_order_size),
            customers = n())

# Locations of the two brick-and-mortar stores
stores <- data.frame(store = c("Northern Virginia (22102)", "Western New Jersey (08053)"),
                     zip_code = c("22102", "08053"))

# Download zip code (ZCTA) and state boundaries from the US Census Bureau
# (the first download may take a few minutes; files are cached for later use)
options(tigris_use_cache = TRUE)
zip_shapes <- zctas(cb = TRUE, year = 2020,
                    starts_with = unique(c(zip_summary$zip_code, stores$zip_code)))
zip_map <- zip_shapes %>%
  inner_join(zip_summary, by = c("ZCTA5CE20" = "zip_code"))
state_map <- states(cb = TRUE, year = 2020) %>%
  filter(STUSPS %in% c("NJ", "PA", "DE", "MD", "DC", "VA", "WV", "NY"))

# Mark each store at the center of its zip code
store_shapes <- zip_shapes %>%
  inner_join(stores, by = c("ZCTA5CE20" = "zip_code"))
stores <- data.frame(store = store_shapes$store,
                     st_coordinates(st_centroid(st_geometry(store_shapes))))

# Map of average order size by zip code
ggplot() +
  geom_sf(data = state_map, fill = "grey95", color = "grey60") +
  geom_sf(data = zip_map, aes(fill = avg_order_size), color = NA) +
  geom_point(data = stores, aes(x = X, y = Y),
             shape = 24, size = 4, fill = "black", color = "white") +
  geom_text(data = stores, aes(x = X, y = Y, label = store),
            nudge_y = 0.25, size = 3.5, fontface = "bold") +
  scale_fill_distiller(palette = "RdYlBu", name = "Avg. Order Size ($)") +
  coord_sf(xlim = st_bbox(zip_map)[c(1, 3)], ylim = st_bbox(zip_map)[c(2, 4)]) +
  labs(title = "Average Order Size by Zip Code",
       subtitle = "Triangles mark the two stores") +
  theme_void()


#####################################################
# Step 2: Determine the Number of Customer Segments #
#####################################################

# The bases: Past Purchase Behaviors and Marketing Efforts
bases <- c("avg_order_size", "avg_order_freq", "crossbuy", "multichannel",
           "per_sale", "tenure", "avg_mktg_cnt", "return_rate")

# Run hierarchical clustering with bases variables
# scale() standardizes each basis so all are measured on the same scale
seg_hclust <- hclust(dist(scale(seg[, bases])), method = "complete")

# Elbow plot for first 10 segments
elbow <- data.frame(segments = 1:10,
                    height = sort(seg_hclust$height, decreasing = TRUE)[1:10])

ggplot(elbow, aes(x = segments, y = height)) +
  geom_line(color = "blue") +
  geom_point(size = 2) +
  scale_x_continuous(breaks = 1:10) +
  labs(title = "Elbow Plot from the Hierarchical Clustering Algorithm",
       x = "Number of Segments", y = "Distance (Height)") +
  theme_minimal()

# The kink in the curve around six segments suggests a six-segment solution.


#####################################################
# Step 3: K-Means Clustering with Six Segments      #
#####################################################

# Run k-means with 6 segments
seg_kmeans <- kmeans(x = seg[, bases], centers = 6)

# Add segment number back to original data
segmentation_result <- seg %>%
  mutate(segment = seg_kmeans$cluster)

# Number of customers in each segment
table(segmentation_result$segment)

# Export data to a CSV file (used as the input data for the Chapter 4 example)
write.csv(segmentation_result, file = file.choose(new = TRUE), row.names = FALSE) ## Name file segmentation_result.csv


#####################################################
# Step 4: Profile the Segments on the Bases         #
#####################################################

# Sales per year = average order size x average order frequency x (1 - return rate)
# Profit per year = sales per year x margin (52%) - catalogs x cost per catalog ($0.75)
segmentation_result <- segmentation_result %>%
  mutate(sales_per_year = avg_order_size * avg_order_freq * (1 - return_rate),
         profit_per_year = sales_per_year * 0.52 - avg_mktg_cnt * 0.75)

# Average values of the bases by segment, plus segment size, total sales, and total profit
bases_profile <- segmentation_result %>%
  group_by(segment) %>%
  summarise(across(all_of(bases), mean),
            segment_size = n(),
            total_sales = sum(sales_per_year),
            total_profit = sum(profit_per_year))

# Show segments as columns, as in the book
bases_table <- bases_profile %>%
  pivot_longer(-segment, names_to = "measure") %>%
  pivot_wider(names_from = segment, values_from = value)

cat("\nAverage Values of Bases by Segment\n")
print(as.data.frame(mutate(bases_table, across(-measure, ~ round(.x, 2)))), row.names = FALSE)


#####################################################
# Step 5: Profile the Segments on the Descriptors   #
#####################################################

descriptors <- c("age", "household_size", "income", "loyalty_card", "married", "own_home")

# Average values of the descriptors by segment, plus segment size
descriptors_profile <- segmentation_result %>%
  group_by(segment) %>%
  summarise(across(all_of(descriptors), mean),
            segment_size = n())

descriptors_table <- descriptors_profile %>%
  pivot_longer(-segment, names_to = "measure") %>%
  pivot_wider(names_from = segment, values_from = value)

cat("\nAverage Values of Descriptors by Segment\n")
print(as.data.frame(mutate(descriptors_table, across(-measure, ~ round(.x, 2)))), row.names = FALSE)


#####################################################
# Step 6: Segment Attractiveness                    #
#####################################################

# Revenue and profit from each segment per year
attractiveness <- bases_profile %>%
  mutate(sales_per_customer = total_sales / segment_size,
         profit_per_customer = total_profit / segment_size,
         share_of_customers = segment_size / sum(segment_size)) %>%
  select(segment, segment_size, share_of_customers, total_sales,
         sales_per_customer, total_profit, profit_per_customer)

cat("\nRevenue and Profit from Each Segment per Year\n")
print(as.data.frame(mutate(attractiveness, across(-segment, ~ round(.x, 2)))), row.names = FALSE)