##############################################################
# Chapter 16 - Using Marketing Experiments to Optimize the   #
#              Marketing Mix                                 #
# Chestnut Ridge Email Campaign Propensity Score Matching    #
# (R)                                                        #
##############################################################

# This script reproduces the full Chestnut Ridge email campaign propensity score matching example in R:
#   Step 1. Install packages (if needed), load packages, and set seed
#   Step 2. Read in the email marketing campaign data
#   Step 3. Look at pre-matched data
#   Step 4. Difference of covariates (pre-matching means)
#   Step 5. Determine optimal caliper size
#   Step 6. Execute the matching algorithm
#   Step 7. Tobit model with the matched data
#   Step 8. Predict the expected value of purchase_amt
#   Step 9. Compare the expected value for the email and no email groups (matched)
#   Step 10. Run the Tobit model with the unmatched data
#   Step 11. Look at the propensity scores before and after matching
#   Step 12. Create the revenue change variables
#   Step 13. Compare predicted and actual revenue gains (Figure 16.2)
#   Step 14. Estimated gain per customer from the email campaign
# Run it from top to bottom. Each step prints its result to the console or Plots pane.


##############################################################
# Step 1: Install packages (if needed), load packages, and   #
#         set seed                                           #
##############################################################

packages <- c("MatchIt", "dplyr", "AER", "ggplot2", "tidyr")
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) install.packages(new_packages)

library(MatchIt)
library(dplyr)
library(AER)
library(ggplot2)
library(tidyr)
set.seed(1)


##############################################################
# Step 2: Read in the email marketing campaign data          #
##############################################################

psmatch <- read.csv(file.choose())  ## Choose the file retail_psmatch.csv


##############################################################
# Step 3: Look at pre-matched data                           #
##############################################################

psmatch %>%
  group_by(email) %>%
  summarise(n_customers = n(),
            mean_purchase_amt = mean(purchase_amt),
            .groups = 'keep')


##############################################################
# Step 4: Difference of covariates (pre-matching means)      #
##############################################################

psmatch_keep <- c('email', 'revenue', 'number_of_orders', 'number_of_orders2',
                  'recency_days')
psmatch_cov <- subset(psmatch, select = psmatch_keep)
psmatch_cov %>%
  group_by(email) %>%
  summarise_all(mean)


##############################################################
# Step 5: Determine optimal caliper size                     #
##############################################################

ps_match <- glm(email ~ revenue + number_of_orders + number_of_orders2 +
                  recency_days, family = binomial(), data = psmatch)
ps_match_df <- data.frame(pr_score = predict(ps_match, type = "response"),
                          email = ps_match$model$email)
0.2*sd(ps_match_df$pr_score)


##############################################################
# Step 6: Execute the matching algorithm                     #
##############################################################

## Remove any rows with missing values
psmatch_nomiss <- psmatch %>%
  select(customer, purchase_amt, one_of(psmatch_keep)) %>%
  na.omit()

## Nearest neighbor 1-to-1 matching with a caliper
mod_match <- matchit(email ~ revenue + number_of_orders + number_of_orders2 +
                       recency_days,
                     method = "nearest", ratio = 1, caliper = 0.023,
                     replace = FALSE, data = psmatch_nomiss)
matched_data <- match.data(mod_match)
dim(matched_data)

## Look at post-matching means
matched_data_cov <- subset(matched_data, select = psmatch_keep)
matched_data_cov %>%
  group_by(email) %>%
  summarise_all(mean)


##############################################################
# Step 7: Tobit model with the matched data                  #
##############################################################

tobit_treat <- tobit(purchase_amt ~ email + revenue + number_of_orders +
                       number_of_orders2 + recency_days,
                     left = 0, data = matched_data)
summary(tobit_treat)


##############################################################
# Step 8: Predict the expected value of purchase_amt         #
##############################################################

mu <- fitted(tobit_treat)
sigma <- tobit_treat$scale
p0 <- pnorm(mu/sigma)
lambda <- function(x) dnorm(x)/pnorm(x)
ey0 <- mu + sigma * lambda(mu/sigma)
matched_data$pred_purchase_amt <- p0 * ey0


##############################################################
# Step 9: Compare the expected value for the email and no    #
#         email groups (matched)                             #
##############################################################

matched_data %>%
  group_by(email) %>%
  summarise(mean_pred_purchase_amt = mean(pred_purchase_amt),
            .groups = 'keep')


##############################################################
# Step 10: Run the Tobit model with the unmatched data       #
##############################################################

tobit_treat_nomatch <- tobit(purchase_amt ~ email + revenue + number_of_orders +
                               number_of_orders2 + recency_days,
                             left = 0, data = psmatch)
summary(tobit_treat_nomatch)

## Predict the expected value for the unmatched data
mu <- fitted(tobit_treat_nomatch)
sigma <- tobit_treat_nomatch$scale
p0 <- pnorm(mu/sigma)
ey0 <- mu + sigma * lambda(mu/sigma)
psmatch$pred_purchase_amt <- p0 * ey0

## Compare the expected value for the email and no email groups (no match)
psmatch %>%
  group_by(email) %>%
  summarise(mean_pred_purchase_amt = mean(pred_purchase_amt),
            .groups = 'keep')


##############################################################
# Step 11: Look at the propensity scores before and after    #
#          matching                                          #
##############################################################

## (distance is the propensity score that matchit() estimated)
ps_all <- data.frame(pr_score = mod_match$distance,
                     email = psmatch_nomiss$email,
                     sample = "Before Matching")
ps_matched <- data.frame(pr_score = matched_data$distance,
                         email = matched_data$email,
                         sample = "After Matching")
ps_compare <- rbind(ps_all, ps_matched)
ps_compare$sample <- factor(ps_compare$sample,
                            levels = c("Before Matching", "After Matching"))
ps_compare$email <- factor(ps_compare$email, levels = c(0, 1),
                           labels = c("No Email", "Email"))

ggplot(ps_compare, aes(x = pr_score, fill = email)) +
  geom_histogram(binwidth = 0.025, boundary = 0, position = "identity",
                 alpha = 0.6, color = "white") +
  facet_grid(sample ~ .) +
  scale_fill_manual(values = c("No Email" = "#2a78d6", "Email" = "#eb6834")) +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Propensity Scores Before and After Matching",
       x = "Propensity Score (Probability of Receiving the Email)",
       y = "Number of Customers", fill = NULL) +
  theme_minimal() +
  theme(legend.position = "top",
        panel.grid.minor = element_blank())


##############################################################
# Step 12: Create the revenue change variables               #
##############################################################

## Compare the purchase amount in the campaign month (January 2020) to the
## average monthly spending in the prior two years (revenue / 24 months)
matched_data$revenue_change <- matched_data$purchase_amt -
  (matched_data$revenue / 24)
matched_data$pred_revenue_change <- matched_data$pred_purchase_amt -
  (matched_data$revenue / 24)


##############################################################
# Step 13: Compare predicted and actual revenue gains        #
#          (Figure 16.2)                                     #
##############################################################

gains <- matched_data %>%
  group_by(email) %>%
  summarise(avg_pred_purchase_amt = mean(pred_purchase_amt),
            avg_pred_revenue_change = mean(pred_revenue_change),
            avg_purchase_amt = mean(purchase_amt),
            avg_revenue_change = mean(revenue_change))

## Put the measures in rows and the email groups in columns
gains_table <- gains %>%
  pivot_longer(cols = -email, names_to = "measure", values_to = "value") %>%
  pivot_wider(names_from = email, names_prefix = "email_",
              values_from = value) %>%
  mutate(difference = email_1 - email_0)
gains_table$measure <- c("Avg. Pred Purchase Amt", "Avg. Pred Revenue Change",
                         "Avg. Purchase Amt", "Avg. Revenue Change")
gains_table %>%
  mutate(across(where(is.numeric), ~ round(.x, 2)))

## Plot the comparison
gains_long <- gains %>%
  pivot_longer(cols = -email, names_to = "measure", values_to = "value")
gains_long$measure <- factor(gains_long$measure,
                             levels = c("avg_pred_purchase_amt",
                                        "avg_pred_revenue_change",
                                        "avg_purchase_amt",
                                        "avg_revenue_change"),
                             labels = c("Avg. Pred Purchase Amt",
                                        "Avg. Pred Revenue Change",
                                        "Avg. Purchase Amt",
                                        "Avg. Revenue Change"))
gains_long$email <- factor(gains_long$email, levels = c(0, 1),
                           labels = c("No Email", "Email"))

ggplot(gains_long, aes(x = measure, y = value, fill = email)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = scales::dollar(value, accuracy = 0.01)),
            position = position_dodge(width = 0.8), vjust = -0.4, size = 3.5) +
  scale_fill_manual(values = c("No Email" = "#2a78d6", "Email" = "#eb6834")) +
  scale_y_continuous(labels = scales::dollar,
                     expand = expansion(mult = c(0, 0.1))) +
  labs(title = "Comparison of Predicted and Actual Revenue Gains",
       x = NULL, y = "Average per Customer", fill = NULL) +
  theme_minimal() +
  theme(legend.position = "top",
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank())


##############################################################
# Step 14: Estimated gain per customer from the email        #
#          campaign                                          #
##############################################################

## Difference in average purchase amount (Email - No Email)
gains_table$difference[gains_table$measure == "Avg. Purchase Amt"]
gains_table$difference[gains_table$measure == "Avg. Pred Purchase Amt"]

## Difference-in-differences (Email - No Email change vs. prior monthly average)
gains_table$difference[gains_table$measure == "Avg. Revenue Change"]
gains_table$difference[gains_table$measure == "Avg. Pred Revenue Change"]