back_end.Rmd

---
title: "MRR analysis tool"
date: "3/24/2021"
output: html_document
---

This tool has three major parts:   
1) Conversion rate analysis  
2) Mean MRR per NPC/PC analysis   
3) Expected MRR per visitor (or signup) analysis which merges the previous two parts.  

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

```{r warning=FALSE, message = FALSE}
library(tidyverse)
library(brms)
library(magrittr)
library(forcats)
library(tidyr)
library(modelr)
library(ggdist)
library(tidybayes)
library(cowplot)
library(ggrepel)
library(RColorBrewer)
library(bayestestR)
library(scales)
library(kableExtra)
library(patchwork)
```

## Conversion rate analysis

Here visitors and conversions are entered manually. You could actually use singups as the sample size and take the rate of NPC-s from that, which would give insight into MRR per signup (In this document we are getting MRR per visitor in the end).  

We are essentially using a binomial general linear model with uninformative priors that the brms() finds automatically. Since the model outputs logits we write a function to transform them back into probabilities
```{r warning=FALSE, message=FALSE, results='hide'}
##DATA PREP AND MODEL OF CONVERSION RATE

#Aggregate data entry for binomial test
conv_data2 <- tibble(variation = c("A", "B"),
                     n = c(1900, 2351),
                     conversions =c(68, 75))

#Binomial model
fit <- brm(family = binomial,
           conversions | trials(n) ~ variation,
           data = conv_data2,
           iter = 2000,
           warmup = 500,
           refresh = 0)

#Function for transforming logits to probability
logit2prob <- function(logit){
  odds <- exp(logit)
  prob <- odds / (1 + odds)
  return(prob)
}
```

The table shows model implied conversion rate statistics. Probability of outperforming is calculated by looking at the porportion of times the conversion rate is bigger in one variation as compared to the other in the posterior distribution. Uplift if the percentage difference in CR.
```{r}
## OVERVIEW TABLES BASED ON CONVERSION RATE MODEL

#Calculating the probability of out performing for each variation
outperforming <- posterior_samples(fit) %>%
  select(1:2) %>%
  mutate(varA_per = logit2prob(.[,1]) , VarB_per = logit2prob(.[,1]+.[,2])) %>%
  select(varA_per, VarB_per) %>% 
  summarise(a_better = paste0(round((sum(.[,1]>.[,2])/n())*100,2),"%"),
            b_better = paste0(round((sum(.[,2]>.[,1])/n())*100,2),"%")) %>% 
  gather(key = variation, value = `Prob. of outperforming`, a_better, b_better) %>% 
  select(`Prob. of outperforming`)

#Calculating uplift uring the back-transformed (with logit2prob) model coefficinets
uplift <- function(x,y){
   b_rate = logit2prob(x+y)
   dif = b_rate - logit2prob(x)
   uplift = dif/logit2prob(x)
   return(uplift)
}

#Creating the table from the numbers
conv_tab <- conv_data2 %>% 
  cbind(., outperforming) %>%
  mutate(CR = paste0(round((conversions/n)*100,2),"%")) %>%
  mutate(uplift = c("", paste0(round(uplift(posterior_summary(fit)[1,1],posterior_summary(fit)[2,1])*100,2
),"%"))) %>% 
  kbl(caption = "Conversion rate statistics") %>% 
  kable_styling()

conv_tab
```

Plot of the model implied  probability distribution of possible conversion rates for the two variations and the plot of distribution of difference of CR between the variations.  The thicker interval beneath the curve is the 66% interval and the thinner one is the 95% one. They show the probability of the true conversion rate being in a certain interval.
```{r warning=FALSE, message=FALSE, results='hide'}
## DENSITY PLOTS FOR CONV.RATE MODEL
  
#Conversion rates distribution for both variations side by side
conv_r_dist <- posterior_samples(fit) %>% 
  select(1:2) %>%
  mutate(`Variation A` = logit2prob(.[,1]), `Variation B` = logit2prob(.[,1]+.[,2])) %>%
  select(`Variation A`, `Variation B`) %>% 
  gather(key = "variation", value = "CR", `Variation A`, `Variation B`) %>% 
  ggplot(aes(x = CR, fill = variation, color = fct_rev(variation))) +
  ggtitle("Conversion (visitor to pc) rates for the variations")+
  stat_slab(alpha = .5)+
  stat_pointinterval(position = position_dodge(width = .4, preserve = "single"))+
  scale_fill_brewer(palette="Dark2")+
  scale_x_continuous(name  = "Conversion Rate", labels = percent)+
  scale_y_continuous(NULL, breaks = NULL) +
  guides(color = FALSE, fill = guide_legend(title=NULL))+
  theme_bw()

#Distribution of difference of conversion rates
conv_dif_dist <- posterior_samples(fit) %>%
  select(1:2) %>%
  mutate(varA_per = logit2prob(.[,1]), VarB_per = logit2prob(.[,1]+.[,2])) %>%
  select(varA_per, VarB_per) %>% 
  transmute(dif = (VarB_per-varA_per)/mean(varA_per)) %>% 
  ggplot(aes(x = dif)) +
  ggtitle("Distribution of difference in conversion rate (visitor to pc)")+
  stat_slab(alpha = .5, fill = "#7570B3")+
  stat_pointinterval(position = position_dodge(width = .4, preserve = "single"))+
  scale_x_continuous(name  = "Difference in conversion Rate", labels = percent)+
  scale_y_continuous(NULL, breaks = NULL) +
  guides(color = FALSE, fill = guide_legend(title=NULL))+
  theme_bw() 

conv_dif_tab <- posterior_samples(fit) %>%
      select(1:2) %>%
      transmute(dif = uplift(.[[1]], .[[2]])) 

conv_r_dist
conv_dif_dist
```

## MRR per NPC analysis

From here on out everything pertains to a .csv datafile from zeppelin. This first code chunk here leads to a descriptive statistics about the MRR of the **raw data**. No models have been applied yet.

1) mean   
2) standard deviation  
3) max MRR contribution within a variation  
4) min MRR contribution within a variation  
5) 25% quantile: the value from which 25% of values (by count) are lower (and 75% higher)  
6) median: the value under which and over which 50% of values lie  
7) 75% quantile: the value from which 75% of value are lower (and 25% higher)

The graph displays the distribution of the actual data points (in purple) or in other words each company that became an NPC during the test with the MRR on the y axis.

The manta-ray looking things are violin-plots which are basically the regular probability distribution plot (like the ones above) but turned on its side and mirrored on both sides. So a more pronounced the peak means more values (PC-s) in that area of MRR. 

The box-plots in the middle show the median (the middle line of the box) and the edges of the boxes show the 25th and 75th quantiles so 50% of the data lie within the box. While also 50% of the data lie on both sides of the central line in the box.
```{r warning=FALSE, message=FALSE}
#EXPLORING RAW MRR DATA TABLES AND DISTRIBUTION PLOT
andmed <- read_csv("download (5).csv")

var_names <- andmed %>% 
  distinct(variation)

#Reading in the data from a .csv file
npc_dat <- andmed %>%
  select(1,2) %>% 
  mutate(variation = recode(variation, `hpttgswhlopta` = "A", `hpttgswhloptb`="B")) %>% #designate A and B
  rename(firstmrr = 2)  #rename second column to "firstmrr"

#Table for raw data with no transformations
raw_tab <- npc_dat %>%
  #filter(.[[2]] < 300) %>% #Possible to put in filter for removing outliers. Look into making optional? 
  group_by(variation) %>% 
  summarise(mean = mean(firstmrr),
            SD = sd(firstmrr),
            max = max(firstmrr),
            min= min(firstmrr),
            Q25 = quantile(firstmrr, prob = c(.25)),
            median = median(firstmrr),
            Q75 = quantile(firstmrr, prob = c(.75))) %>% 
  kbl(caption="Overview of raw MRR data") %>% 
  kable_styling()

#Plot for distribution of log-transformed raw data. (Look into making transformation optional?)
raw_dist <-  npc_dat %>% 
  ggplot(aes(x = as.factor(variation), y = log(firstmrr), fill = as.factor(variation))) +
  ggtitle("Distribution of log-transformed raw MRR data")+
  geom_violin(trim=FALSE, alpha =0.5)+
  geom_jitter(shape=10, position=position_jitter(0.05), color = "#7570B3")+
  geom_boxplot(width=0.15, fill = "white", alpha = 0.5, color = "black" )+
  scale_fill_brewer(palette="Dark2")+
  scale_x_discrete(name = "Variation")+
  scale_y_continuous(name = "MRR", n.breaks = 20)+
  theme_bw()

raw_tab
raw_dist
```


For the model we are using an alternative parametrization (no intercept) to make variation comparisons easier. Also using a shifted_lognormal model here because the data can't be under 0, is skewed to the right and is also shifted to the right because MRR values usually don't fall under a certain value thats above 0.
```{r warning=FALSE, message=FALSE, results='hide'}
#DATA PREP AND MODEL BUILD FOR MRR ANALYSIS

#preparing table with parametrization of variales for logit regression without intercept
npc_wide <- npc_dat %>% 
  select(1, 2) %>% 
  mutate(A = as.numeric(as.character(factor(variation, labels=c(1,0))))) %>% 
  mutate(B = 1-A)

#Bayesian logit model with log-transformed data.
b4.3 <- 
  brm(data = npc_wide, family = shifted_lognormal(link_sigma = "identity"),
      bf(firstmrr ~ 0 + A + B, sigma ~ 0 + A + B), #sigma is modeled as also depending on variation
      #prior = c(prior(student_t(3,0,1), class = b)), #Look into  more informative priors for accuracy managment
      control = list(adapt_delta = 0.95),
      iter = 3000, warmup = 600, chains = 4, cores = 4, #Look into this for managing computing time 
      seed = 4)
```

Here one sample of the model generated posterior predictive distribution can be compared to the original data to see if the model understands the data well enough to re generate it reasonably well just from knowing the variations. In the model the medians are quite unstable between samples so the box plot and the median implied by it is a very general view of where the model thinks the data lies and is just for general model performance assessment. For the  calculations of means later we are using a larger sample that is much more stable.
```{r warning=FALSE, message=FALSE}
#CHECKING MODEL CONVERGEANCE AND SIMILIARITY TO REAL WORLD

#write "plot(b4.3)" for convergence plot (Questionable if needed in the final output)

#Adding model predicted datapoints to raw data point table
post_pred <- as.data.frame(t(posterior_predict(object = b4.3, 
                                                       newdata = npc_wide[,3:4], nsamples = 1)))
npc_dat2 <- cbind(npc_wide[,1:2], post_pred)

#Comparing real world log-transformed raw data with model implied log-transformed data
model_dist <- npc_dat2 %>% 
  ggplot(aes(x = as.factor(variation), y = log(npc_dat2[[3]]), fill = as.factor(variation))) +
  ggtitle("Distribution of model predicted data \n (log-transformed)")+
  geom_violin(trim=FALSE, alpha =0.5)+
  geom_jitter(shape=10, position=position_jitter(0.05), color = "#7570B3")+
  geom_boxplot(width=0.15, fill = "white", alpha = 0.5, color = "black" )+
  scale_fill_brewer(palette="Dark2")+
  scale_x_discrete(name = "Variation")+
  scale_y_continuous(name = "MRR", n.breaks = 20, limits = ggplot_build(raw_dist)$layout$panel_scales_y[[1]]$range$range)+
  guides(fill=guide_legend(title = "Variations"))

#Plot of posterior predictive dist next to original data aka what the model predicts only knowing wether its variation A or B
compare_plot <- raw_dist | model_dist
compare_plot
```

This part is about model implied mean of MRR 

The table here shows the output of the model which simulates mean MRR per NPC.  

2) **Estimate of mean** column tell you what the models implies the most likley mean MRR per one NPC is for the variations.    
3) **Est.Error of mean** column tells you what the standard-deviation of that prediction of mean MRR per PC is in the model.
4) **Q2.5 and Q97.5** are basically the thin lines on the graph below or in other words the borders of MRR values within which the true MRR per NPC lies in with 95% probability according to the model.  
5) **Prob. of outperforming** tells you the probability that the mean MRR per NPC is bigger for one or the other variation

The plots shows the probability distribution of  **mean MRR per PC** for the two variations and the probability distribution for the difference of **the mean MRR per NPC**.While the violin plots show the distribution of the data (the actual PC-s), this plot here shows the probability distribution of possible **mean per PC** and not the distribution of the individual PCs as data points. The plots describe: "given that a variation gets a NPC, what is their expected MRR" or with the distribution of difference what is the expected difference between an average NPC of two variations.

```{r warning=FALSE, message=FALSE}
#MRR MODEL TABLES AND PLOTS

#Model sample of means
post_dist <- as.data.frame(t(posterior_epred(b4.3, nsamples = 6000)))
mean_samp <- cbind(npc_dat$variation, post_dist) %>%
  distinct() %>%
  pivot_longer(-`npc_dat$variation`) %>%    
  pivot_wider(names_from=`npc_dat$variation`, values_from=value) %>% 
  select(A,B) 


#Calculations with model distribution of means
mrr_mean <- mean_samp %>%
  pivot_longer(cols = c(A, B), values_to = "means", names_to = "Variation")  %>% 
  group_by(Variation) %>% 
  summarise(`Estimate of mean` = round(mean(means),2),
            `Est. error of mean` = round(sd(means),2),
            Q2.5 = round(quantile(means, prob = c(.025)),2),
            Q97.5 = round(quantile(means, prob = c(.975)),2))

#Calculating probabikity of each variation having higher true mean MRR
mrr_better <- mean_samp %>%  
  summarise(a_better = paste0(round((sum(.[,1]>.[,2])/n())*100,2),"%"),
            b_better = paste0(round((sum(.[,2]>.[,1])/n())*100,2),"%")) %>% 
  pivot_longer(cols = c(a_better, b_better), values_to = "Prob. of outperforming (on mean)") %>% 
  select("Prob. of outperforming (on mean)")

#Table of model results
mrr_tab <- cbind(mrr_mean, mrr_better) %>%
  kbl(caption = "Table of model implied estimates of mean/median MRR per PC per variation") %>% 
  kable_styling()

mean_dif <- mean_samp %>% 
  transmute(dif = B - A)

#Distribution plot of mean MRR per NPC for both variations
dist_mean <- mean_samp %>%
  gather(key = "variation", value = "mean_mrr", A, B) %>% 
  ggplot(aes(x = mean_mrr, fill = variation , color = fct_rev(variation)))+
  ggtitle("Distributions of model implied mean MRR per PC per variation")+
  geom_density(alpha = 0.7)+
  #stat_slab(alpha = .5)+
  #stat_pointinterval(position = position_dodge(width = .4, preserve = "single"))+
  scale_fill_brewer(palette="Dark2")+
  scale_x_continuous(name  = "Mean MRR per PC", labels = dollar_format())+
  scale_y_continuous(NULL, breaks = NULL) +
  guides(color = FALSE, fill = guide_legend(title=NULL))+
  theme_bw() 

#Distribution plot of difference of model implied mean MRR per PC
dist_dif_mean <- mean_dif %>% 
  ggplot(aes(x = dif)) +
  ggtitle("Distribution of model implied difference of mean MRR per NPC between variations")+
  stat_slab(alpha = .5, fill = "#7570B3")+
  stat_pointinterval(position = position_dodge(width = .4, preserve = "single"))+
  scale_y_continuous(NULL, breaks = NULL) +
  scale_x_continuous(name  = "mean MRR per PC of A - mean MRR per PC of B", labels = dollar_format())+
  theme_bw()

mrr_tab
dist_mean
dist_dif_mean
```
###  MRR per visitor (or signup) analysis   


In the third and last part where both of the previous parts are combined. Conceptually this means that we take the probability of a single **visitor** being converted (conversion rate from visitor to PC) and multiply that probability by the value they are expected to produce in MRR. With that we get the expected value of a single **visitor**. The great thing here is that, if in the first part of conversion rate analysis you don't enter the **visitor** and "PC" numbers but you enter the **"signup"** and "PC" numbers, then you will get here the expected MRR per **signup** not per **visitor**. So it is possible to see how much MRR a variation is expected to produce per visitor OR per singup. Thats why it is good to be able to enter the proportions manually in the beginning. 

```{r warning=FALSE, message=FALSE}
#SET UP FOR CONVERSION RATE AND MRR MERGED ANALYSIS

#Transforming logits to probability of conversion for variations samples
post_samp_transformed <- posterior_samples(fit) %>%
  select(1:2) %>%
  mutate(varA_per = logit2prob(.[,1]), VarB_per = logit2prob(.[,1]+.[,2])) %>%
  select(varA_per, VarB_per) 

#Calculating MRR per visitor based on conversion and mean MRR
overall_results <- mean_samp%>%
  select(A, B) %>% 
  mutate(A_overall = A*post_samp_transformed$varA_per,
         B_overall = B*post_samp_transformed$VarB_per) %>% 
  select(A_overall, B_overall)
```

In the table 

1) **Mean** gives you the expected MRR per each visitor.  
2) **SD** gives you the standard deviation of that mean MRR per visitor. 
3) **Quantiles** give the points between which with 95% probability the MRR per visitor lies.   
4) **Prob. of outperforming**  tells you the probability of per visitor MRR being bigger in one or the other variation.  

The first plot is no different from the previous similar ones. It shows the distribution of expected MRR per **visitor** for both variations. For example what for variation A is the most probable ammount of MRR that variation A will make for every person that visits the Pipedrive website in that variaton.

The last plot shows again the probability distribution of the difference of MRR per visitor to the site between the two variations. This is basically where the Prob. of outperforming comes from (how much of the density is under or over 0 EUR).
```{r warning=FALSE, message=FALSE}
#TABLES AND PLOTS BASED ON CR AND MRR MODEL MERGED ANALYSIS

#Calculating probability of outperforming based on CR and MRR
overall_better <- overall_results%>% 
  summarise(a_better = paste0(round((sum(.[,1]>.[,2])/n())*100,2),"%"),
            b_better = paste0(round((sum(.[,2]>.[,1])/n())*100,2),"%")) %>% 
  gather(key = variation, value = `Prob. of outperforming`, a_better, b_better) %>% 
  select(`Prob. of outperforming`)

#Table for performance based on CR and MRR
tab_overall <- overall_results %>%
  rename(A = "A_overall", B = "B_overall") %>% 
  gather(key="Variation", value = "Values", A, B)%>%
  group_by(Variation) %>% 
  summarise(Mean = round(mean(Values),2),
            SD = round(sd(Values),2),
            Q2.5 = round(quantile(Values, prob = c(.025)),2),
            Q97.5 = round(quantile(Values, prob = c(.975)),2)) %>% 
  cbind(overall_better) %>%
  mutate(`MRR per 1000 visitors` = 1000*Mean) %>% 
  kbl(caption="MRR per visitor/signup") %>% 
kable_styling()

#Most likely difference 
mean_difference <- overall_results %>% 
  summarise(mean(A_overall-B_overall))

#How much we would expect to make from one VISITOR (since most of them are zero)
dist_overall <- overall_results %>%
  gather(key = "variation", value = "overall", A_overall, B_overall) %>% 
  ggplot(aes(x = overall, fill = variation , color = fct_rev(variation)))+
  ggtitle("Model implied distributions of MRR per visitor")+
  stat_slab(alpha = .5)+
  stat_pointinterval(position = position_dodge(width = .4, preserve = "single"))+
  scale_fill_brewer(palette="Dark2")+
  scale_x_continuous(name  = "MRR per visitor", labels = dollar_format())+
  scale_y_continuous(NULL, breaks = NULL) +
  guides(color = FALSE, fill = guide_legend(title=NULL))+
  theme_bw() 

#MRR difference per visitor
dist_dif_overall <- overall_results %>% 
  transmute(dif2 = A_overall-B_overall) %>% 
  ggplot(aes(x = dif2)) +
  ggtitle("Model implied distribution of difference of overall MRR")+
  stat_slab(alpha = .5, fill = "#7570B3")+
  stat_pointinterval(position = position_dodge(width = .4, preserve = "single"))+
  scale_y_continuous(NULL, breaks = NULL) +
  scale_x_continuous(name  = "Difference in overall revenue", labels = dollar_format(), n.breaks = 20)+
  theme_bw()

tab_overall
dist_overall
dist_dif_overall
```