#######################
# Logistic Regression #
#######################

# This script reproduces the full Chestnut Ridge logistic regression example in R:
#   Step 1. Logistic regression model (fit, pseudo R-square, odds ratios, predicted probabilities)
#   Step 2. Breakeven rate for the catalog campaign
#   Step 3. Predicted vs. actual purchases
#   Step 4. Confusion matrix and hit rate
#   Step 5. Profitability and ROMI compared to random and RFM targeting
# 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("ggplot2")
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) install.packages(new_packages)

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

## Read in Logistic Regression data
logit <- read.csv(file.choose()) ## Choose retail_logit.csv file


############################################
# Step 1: Logistic Regression Model        #
############################################

## Transform and Create Data
logit$number_of_orders2 <- logit$number_of_orders^2
logit$lnrevenue <- log(logit$revenue + 1)

## Run Logistic Regression using GLM
logit_result <- glm(formula = purchase ~ lnrevenue + number_of_orders + number_of_orders2 + recency_days +
	loyalty_card + married + income, data = logit, family = "binomial")
summary(logit_result)

## Pseudo R-square - McFadden
null_result <- glm(formula = purchase ~ 1, data = logit, family = "binomial")
1 - logLik(logit_result)/logLik(null_result)

## Odds Ratio
exp(logit_result$coefficients)

## Predicted Probability
logit$predict <- predict(logit_result, logit, type = "response")


############################################
# Step 2: Breakeven Rate                   #
############################################

## Campaign Economics
avg_purchase <- 40     ## Average purchase amount
cogs_rate <- 0.48      ## Cost of goods sold as a share of revenue ($19.20 / $40)
ship_cost <- 6         ## Average cost to ship the product
mail_cost <- 2         ## Average cost of the marketing campaign (catalog) per customer

## Profit per Purchasing Customer
profit_per_buyer <- avg_purchase - avg_purchase * cogs_rate - ship_cost - mail_cost
profit_per_buyer

## Breakeven Rate
breakeven <- mail_cost / profit_per_buyer
breakeven


############################################
# Step 3: Predicted vs. Actual             #
############################################

## Target Customers at or above the Breakeven Rate
logit$predict_purchase <- ifelse(logit$predict >= breakeven, 1, 0)

## Predicted vs. Actual Purchase for Each Customer
head(logit[, c("customer_id", "predict", "predict_purchase", "purchase")], 20)

## Predicted vs. Actual Purchase, Sorted by Predicted Probability
logit_sorted <- logit[order(logit$predict, decreasing = TRUE), ]
head(logit_sorted[, c("customer_id", "predict", "predict_purchase", "purchase")], 20)

## Distribution of Predicted Probability for Purchasers and Non-Purchasers
ggplot(logit, aes(x = predict,
	fill = factor(purchase, levels = c(0, 1), labels = c("Did Not Purchase", "Purchased")))) +
	geom_histogram(binwidth = 0.025, boundary = 0, position = "identity",
		alpha = 0.6, color = "white") +
	geom_vline(xintercept = breakeven, linetype = "dashed", linewidth = 1) +
	ggplot2::annotate("text", x = breakeven, y = Inf, vjust = 1.5, hjust = -0.05,
		label = paste0("Breakeven Rate = ", format(floor(breakeven * 10000 + 0.5) / 100, nsmall = 2), "%")) +
	scale_fill_manual(values = c("Did Not Purchase" = "#3366CC", "Purchased" = "#E68A1A")) +
	labs(title = "Predicted Probability of Purchase by Actual Purchase",
		x = "Predicted Probability of Purchase", y = "Number of Customers", fill = NULL) +
	theme_minimal() +
	theme(legend.position = "top")


############################################
# Step 4: Confusion Matrix                 #
############################################

## Confusion Matrix (Hit Rate Table)
confusion <- table(Predicted = logit$predict_purchase, Actual = logit$purchase)
confusion

true_negative <- confusion["0", "0"]
true_positive <- confusion["1", "1"]
false_negative <- confusion["0", "1"]
false_positive <- confusion["1", "0"]

## Hit Rate
hit_rate <- (true_negative + true_positive) / sum(confusion)
hit_rate


############################################
# Step 5: Profitability and ROMI           #
############################################

## Function to Compute Campaign Profitability
campaign_results <- function(n_customers, n_targeted, n_buyers) {
	revenue <- n_buyers * avg_purchase
	cogs <- revenue * cogs_rate
	shipping <- n_buyers * ship_cost
	marketing <- n_targeted * mail_cost
	profit <- revenue - cogs - shipping - marketing
	c(number_of_customers = n_customers,
		number_targeted = n_targeted,
		pct_targeted = n_targeted / n_customers,
		number_of_buyers = n_buyers,
		response_rate = n_buyers / n_targeted,
		revenue = revenue,
		cogs = cogs,
		shipping = shipping,
		targeted_marketing = marketing,
		profit = profit,
		romi = profit / marketing)
}

## Random Targeting (send the catalog to every customer)
random <- campaign_results(nrow(logit), nrow(logit), sum(logit$purchase))

## Targeted: Independent Sort RFM (results from Chapter 7)
rfm_independent <- campaign_results(10000, 4459, 1314)

## Targeted: Sequential Sort RFM (results from Chapter 7)
rfm_sequential <- campaign_results(10000, 4480, 1306)

## Targeted: Logistic Regression
logistic <- campaign_results(nrow(logit), true_positive + false_positive, true_positive)

## Compare Results
results <- data.frame(random, rfm_independent, rfm_sequential, logistic)
round(results, 2)

## Export Logistic Regression Results (used as the input data for the Chapter 9 example)
write.csv(logit, file = file.choose(new=TRUE), row.names = FALSE) ## Name file logit_pred.csv