# Load packages
library(MatchIt)
library(cobalt)

# Only include existing customers, since they are the only ones signing up to the email (for potential customers, email sign up is zero in our dataset)
existing_customers <- consumers[consumers$Existing_customer==1,]

# Log transform CLV, since it is highly skewed
existing_customers$Log_CLV <- log(existing_customers$CLV)

# Conduct matching of consumers who did and did not subscribe to the email (i.e. step 1 and 3 of the PSM process). In this example, the Email_group membership is predicted using the Relation_length and the log of CLV. The matching is automatically performed with the code.
m.out <- matchit(Email_group ~ Relation_length + Log_CLV, data=existing_customers, caliper=0.05)

# Investigate the PSM outcome
summary(m.out, standardize = TRUE)
love.plot(m.out)




# Plot balance diagnostics (i.e., step 4 of the PSM process)
plot(m.out, type = "jitter", interactive = FALSE)

plot(m.out, type = "qq", interactive = FALSE)
bal.plot(m.out, which = "both")
bal.plot(m.out, var.name = "Relation_length", which = "both")
bal.plot(m.out, var.name = "Log_CLV", which = "both")




# Get a dataframe with the matched data
matched_data <- match.data(m.out)

# Plot of email and purchase likelihood (unmatched)
myData <- aggregate(existing_customers$Purchase,
                    by = list(Email_group = existing_customers$Email_group),
                    FUN = function(x) c(mean = mean(x), sd = sd(x),
                                        n = length(x)))

myData <- do.call(data.frame, myData)

myData$se <- myData$x.sd / sqrt(myData$x.n)
colnames(myData) <- c("Email_group", "mean", "sd", "n", "se")
myData$names <- c(paste(myData$Email_group, "Email_group"))

p <- ggplot(data = myData, aes(x = factor(Email_group), y = mean))
p + geom_bar(stat = "identity",
             position = position_dodge(0.9), fill="steelblue") +
  geom_errorbar(aes(ymax = mean + 2*se,
                    ymin = mean - 2*se), position = position_dodge(0.9),
                width = 0.25) +
  labs(x = "Email_group", y = "Conversion rate") +
  ggtitle("Conversion rate by Email_group (full sample)") +
  scale_y_continuous(labels = function(x) paste0(x*100, "%")) +
  geom_text(position = position_dodge(width= .9), aes(y = mean, label = paste(format(mean*100, digits = 2, nsmall = 2), "%", sep = "")), size = 4, vjust = 5)


# Plot of email and purchase likelihood (matched)
myData <- aggregate(matched_data$Purchase,
                    by = list(Email_group = matched_data$Email_group),
                    FUN = function(x) c(mean = mean(x), sd = sd(x),
                                        n = length(x)))

myData <- do.call(data.frame, myData)

myData$se <- myData$x.sd / sqrt(myData$x.n)
colnames(myData) <- c("Email_group", "mean", "sd", "n", "se")
myData$names <- c(paste(myData$Email_group, "Email_group"))

p <- ggplot(data = myData, aes(x = factor(Email_group), y = mean))
p + geom_bar(stat = "identity",
             position = position_dodge(0.9), fill="steelblue") +
  geom_errorbar(aes(ymax = mean + 2*se,
                    ymin = mean - 2*se), position = position_dodge(0.9),
                width = 0.25) +
  labs(x = "Email_group", y = "Conversion rate") +
  ggtitle("Conversion rate by Email_group (matched sample)") +
  scale_y_continuous(labels = function(x) paste0(x*100, "%")) +
  geom_text(position = position_dodge(width= .9), aes(y = mean, label = paste(format(mean*100, digits = 2, nsmall = 2), "%", sep = "")), size = 4, vjust = 5)


# Load package to make an output table of all models
library(sjPlot)
library(sjmisc)
library(sjlabelled)

# Estimate logistic regression models
model1 <- glm(Purchase ~ Email_group, data=matched_data, family=binomial)
model2 <- glm(Purchase ~ Email_group, data=existing_customers, family=binomial)

model3 <- glm(Purchase ~ Email_group + Firm_banner + Email_group*Firm_banner, data=matched_data, family=binomial)
model4 <- glm(Purchase ~ Email_group + Firm_banner + Email_group*Firm_banner, data=existing_customers, family=binomial)

tab_model(model1, model2, model3, model4, transform = NULL, collapse.ci = TRUE, p.style = "stars")




# Add start channel
consumers_journey2 <- consumers_journey[consumers_journey$Channel_number==1,]
consumers_journey2$Channel_number <- 0
consumers_journey2$Channel_name <- "Start"
consumers_journey2$Channel <- 1

consumers_journey2 <- rbind(consumers_journey, consumers_journey2)

# Add final channel (i.e., conversion or not)
consumers_journey3 <- consumers_journey[consumers_journey$Channel_number==1,]
touchpoints <- aggregate(consumers_journey2$Channel_number,by=list(consumers_journey2$Consumer_ID),FUN=max)
consumers_journey3$Channel_number <- touchpoints$x + 1
consumers_journey3$Channel_name <- ifelse(consumers$Purchase==1, "Purchase", "No Purchase")
consumers_journey3$Channel <- ifelse(consumers$Purchase==1, 10, 11)

consumers_journey2 <- rbind(consumers_journey2, consumers_journey3)

consumers_journey2 <- consumers_journey2[order(consumers_journey2$Consumer_ID, consumers_journey2$Channel_number),]

rm(consumers_journey3)
rm(touchpoints)


# Split data up in firm and charity banner group
consumers_journey2_firm <- consumers_journey2[consumers_journey2$Firm_banner==1,]
consumers_journey2_charity <- consumers_journey2[consumers_journey2$Firm_banner==0,]


# Estimate transition matrix
library(msm)

markov_chain <- statetable.msm(Channel, Consumer_ID, data=consumers_journey2_firm)
markov_chain <- prop.table(markov_chain, margin=1)

markov_chain <- as.data.frame(t(markov_chain))

markov_chain <- markov_chain[,3]
markov_chain[100:121] <- c(rep(0,9),1,0,rep(0,10),1)

library(markovchain)

mcjourney <- new("markovchain", states = c("Start","Banner_impression","Banner_click",
                                           "SEA_product_click","SEA_brand_click","Price_comp_click",
                                           "Email_received", "Email_click", "Direct_visit",
                                           "Purchase", "No purchase"), 
                 transitionMatrix = matrix(data=markov_chain,
                                           byrow= TRUE, nrow = 11))

show(mcjourney)

plot(mcjourney, edge.arrow.size=0.25)




# Estimate the stage consumers are in after 5 and 50 phases of the path to purchase, the initial state is the start state
initialstate <- c(1,0,0,0,0,0,0,0,0,0,0)

after5touchpoints <- initialstate * (mcjourney ^ 5)
after5touchpoints


after50touchpoints <- initialstate * (mcjourney ^ 50)
after50touchpoints




# Estimate transition matrix with banner impressions and clicks dropped out

markov_chain_banner_off <- statetable.msm(Channel, Consumer_ID, data=consumers_journey2_firm)

markov_chain_banner_off <- markov_chain_banner_off[c(1,4:9),c(1,4:11)]

markov_chain_banner_off <- prop.table(markov_chain_banner_off, margin=1)

markov_chain_banner_off <- as.data.frame(t(markov_chain_banner_off))

markov_chain_banner_off <- markov_chain_banner_off[,3]
markov_chain_banner_off[64:81] <- c(rep(0,7),1,0,rep(0,8),1)

mcjourney_banner_off <- new("markovchain", states = c("Start", "SEA_product_click","SEA_brand_click","Price_comp_click",
                                                      "Email_received", "Email_click", "Direct_visit",
                                                      "Purchase", "No purchase"), 
                            transitionMatrix = matrix(data=markov_chain_banner_off,
                                                      byrow= TRUE, nrow = 9))

show(mcjourney_banner_off)

plot(mcjourney_banner_off, edge.arrow.size=0.25)


# Estimate the stage consumers are in after 5 and 50 phases of the path to purchase, the initial state is the start state
initialstate <- c(1,0,0,0,0,0,0,0,0)

after5touchpoints <- initialstate * (mcjourney_banner_off ^ 5)
after5touchpoints


after50touchpoints <- initialstate * (mcjourney_banner_off ^ 50)
after50touchpoints




# Estimate transition matrix for the charity banner group
markov_chain <- statetable.msm(Channel, Consumer_ID, data=consumers_journey2_charity)
markov_chain <- prop.table(markov_chain, margin=1)

markov_chain <- as.data.frame(t(markov_chain))

markov_chain <- markov_chain[,3]
markov_chain[100:121] <- c(rep(0,9),1,0,rep(0,10),1)


mcjourney <- new("markovchain", states = c("Start","Banner_impression","Banner_click",
                                           "SEA_product_click","SEA_brand_click","Price_comp_click",
                                           "Email_received", "Email_click", "Direct_visit",
                                           "Purchase", "No purchase"), 
                 transitionMatrix = matrix(data=markov_chain,
                                           byrow= TRUE, nrow = 11))

show(mcjourney)

plot(mcjourney, edge.arrow.size=0.25)


# Estimate the stage consumers are in after 5 and 50 phases of the path to purchase, the initial state is the start state
initialstate <- c(1,0,0,0,0,0,0,0,0,0,0)

after5touchpoints <- initialstate * (mcjourney ^ 5)
after5touchpoints


after50touchpoints <- initialstate * (mcjourney ^ 50)
after50touchpoints