%md ##Calculating the Probability of Future Customer Engagement

In non-subscription retail models, customers come and go with no long-term commitments, making it very difficult to determine whether a customer will return in the future. Determining the probability that a customer will re-engage is critical to the design of effective marketing campaigns. Different messaging and promotions may be required to incentivize customers who have likely dropped out to return to our stores. Engaged customers may be more responsive to marketing that encourages them to expand the breadth and scale of purchases with us. Understanding where our customers land with regard to the probability of future engagement is critical to tailoring our marketing efforts to them.  

The *Buy 'til You Die* (BTYD) models popularized by Peter Fader and others leverage two basic customer metrics, *i.e.* the recency of a customer's last engagement and the frequency of repeat transactions over a customer's lifetime, to derive a probability of future re-engagement. This is done by fitting customer history to curves describing the distribution of purchase frequencies and engagement drop-off following a prior purchase. The math behind these models is fairly complex but thankfully it's been encapsulated in the [lifetimes](https://pypi.org/project/Lifetimes/) library, making it much easier for traditional enterprises to employ. The purpose of this notebook is to examine how these models may be applied to customer transaction history and how they may be deployed for integration in marketing processes.

%md ###Step 1: Setup the Environment

To run this notebook, you need to attach to a **Databricks ML Runtime** cluster leveraging Databricks version 6.5+. This version of the Databricks runtime will provide access to many of the pre-configured libraries used here.  Still, there are additional Python libraries which you will need to install and attach to your cluster.  These are:</p>

* xlrd
* lifetimes==0.10.1
* nbconvert

To install these libraries in your Databricks workspace, please follow [these steps](https://docs.databricks.com/libraries.html#workspace-libraries) using the PyPI library source in combination with the bullet-pointed library names in the provided list.  Once installed, please be sure to [attach](https://docs.databricks.com/libraries.html#install-a-library-on-a-cluster) these libraries to the cluster with which you are running this notebook.

%md With the libraries installed, let's load a sample dataset with which we can examine the BTYD models. The dataset we will use is the [Online Retail Data Set](http://archive.ics.uci.edu/ml/datasets/Online+Retail) available from the UCI Machine Learning Repository.  This dataset is made available as a Microsoft Excel workbook (XLSX).  Having downloaded this XLSX file to our local system, we can load it into our Databricks environment by following the steps provided [here](https://docs.databricks.com/data/tables.html#create-table-ui). Please note when performing the file import, you don't need to select the *Create Table with UI* or the *Create Table in Notebook* options to complete the import process. Also, the name of the XLSX file will be modified upon import as it includes an unsupported space character.  As a result, we will need to programmatically locate the new name for the file assigned by the import process.

Assuming we've uploaded the XLSX to the */FileStore/tables/online_retail/*, we can access it as follows:

%md The data in the workbook are organized as a range in the Online Retail spreadsheet.  Each record represents a line item in a sales transaction. The fields included in the dataset are:

| Field | Description |
|-------------:|-----:|
|InvoiceNo|A 6-digit integral number uniquely assigned to each transaction|
|StockCode|A 5-digit integral number uniquely assigned to each distinct product|
|Description|The product (item) name|
|Quantity|The quantities of each product (item) per transaction|
|InvoiceDate|The invoice date and a time in mm/dd/yy hh:mm format|
|UnitPrice|The per-unit product price in pound sterling (£)|
|CustomerID| A 5-digit integral number uniquely assigned to each customer|
|Country|The name of the country where each customer resides|

Of these fields, the ones of particular interest for our work are InvoiceNo which identifies the transaction, InvoiceDate which identifies the date of that transaction, and CustomerID which uniquely identifies the customer across multiple transactions. (In a separate notebook, we will examine the monetary value of the transactions through the UnitPrice and Quantity fields.)

%md ###Step 2: Explore the Dataset

To enable the exploration of the data using SQL statements, let's flip the pandas DataFrame into a Spark DataFrame and persist it as a temporary view:

%md Examining the transaction activity in our dataset, we can see the first transaction occurs December 1, 2010 and the last is on December 9, 2011 making this a dataset that's a little more than 1 year in duration. The daily transaction count shows there is quite a bit of volatility in daily activity for this online retailer:

%md We can smooth this out a bit by summarizing activity by month. It's important to keep in mind that December 2011 only consists of 9 days so the sales decline graphed for the last month should most likely be ignored:

NOTE We will hide the SQL behind each of the following result sets for ease of viewing.  To view this code, simply click the **Show code** item above each of the following charts.

%sql -- unique transactions by month

SELECT 
  TRUNC(InvoiceDate, 'month') as InvoiceMonth,
  COUNT(DISTINCT InvoiceNo) as Transactions
FROM orders
GROUP BY TRUNC(InvoiceDate, 'month') 
ORDER BY InvoiceMonth;

%md For the little more than 1-year period for which we have data, we see over four-thousand unique customers.  These customers generated about twenty-two thousand unique transactions:

%sql -- unique customers and transactions

SELECT
 COUNT(DISTINCT CustomerID) as Customers,
 COUNT(DISTINCT InvoiceNo) as Transactions
FROM orders
WHERE CustomerID IS NOT NULL;

%md A little quick math may lead us to estimate that, on average, each customer is responsible for about 5 transactions, but this would not provide an accurate representation of customer activity.

Instead, if we count the unique transactions by customer and then examine the frequency of these values, we see that many of the customers have engaged in a single transaction. The distribution of the count of repeat purchases declines from there in a manner that we may describe as negative binomial distribution (which is the basis of the NBD acronym included in the name of most BTYD models):

%sql -- the distribution of per-customer transaction counts 

SELECT
  x.Transactions,
  COUNT(x.*) as Occurrences
FROM (
  SELECT
    CustomerID,
    COUNT(DISTINCT InvoiceNo) as Transactions 
  FROM orders
  WHERE CustomerID IS NOT NULL
  GROUP BY CustomerID
  ) x
GROUP BY 
  x.Transactions
ORDER BY
  x.Transactions;

%md If we alter our last analysis to group a customer's transactions that occur on the same date into a single transaction - a pattern that aligns with metrics we will calculate later - we see that a few more customers are identified as non-repeat customers but the overall pattern remains the same:

%sql -- the distribution of per-customer transaction counts
     -- with consideration of same-day transactions as a single transaction 

SELECT
  x.Transactions,
  COUNT(x.*) as Occurances
FROM (
  SELECT
    CustomerID,
    COUNT(DISTINCT TO_DATE(InvoiceDate)) as Transactions
  FROM orders
  WHERE CustomerID IS NOT NULL
  GROUP BY CustomerID
  ) x
GROUP BY 
  x.Transactions
ORDER BY
  x.Transactions;

%md Focusing on customers with repeat purchases, we can examine the distribution of the days between purchase events. What's important to note here is that most customers return to the site within 2 to 3 months of a prior purchase.  Longer gaps do occur but significantly fewer customers have longer gaps between returns.  This is important to understand in the context of our BYTD models in that the time since we last saw a customer is a critical factor to determining whether they will ever come back with the probability of return dropping as more and more time passes since a customer's last purchase event:

%sql -- distribution of per-customer average number of days between purchase events

WITH CustomerPurchaseDates
  AS (
    SELECT DISTINCT
      CustomerID,
      TO_DATE(InvoiceDate) as InvoiceDate
    FROM orders 
    WHERE CustomerId IS NOT NULL
    )
SELECT -- Per-Customer Average Days Between Purchase Events
  AVG(
    DATEDIFF(a.NextInvoiceDate, a.InvoiceDate)
    ) as AvgDaysBetween
FROM ( -- Purchase Event and Next Purchase Event by Customer
  SELECT 
    x.CustomerID,
    x.InvoiceDate,
    MIN(y.InvoiceDate) as NextInvoiceDate
  FROM CustomerPurchaseDates x
  INNER JOIN CustomerPurchaseDates y
    ON x.CustomerID=y.CustomerID AND x.InvoiceDate < y.InvoiceDate
  GROUP BY 
    x.CustomerID,
    x.InvoiceDate
    ) a
GROUP BY CustomerID

%md ###Step 3: Calculate Customer Metrics

The dataset with which we are working consists of raw transactional history.  To apply the BTYD models, we need to derive several per-customer metrics:</p>

* **Frequency** - the number of dates on which a customer made a purchase subsequent to the date of the customer's first purchase
* **Age (T)** - the number of time units, *e.g.* days, since the date of a customer's first purchase to the current date (or last date in the dataset)
* **Recency** - the age of the customer (as previously defined) at the time of their last purchase

It's important to note that when calculating metrics such as customer age that we need to consider when our dataset terminates.  Calculating these metrics relative to today's date can lead to erroneous results.  Given this, we will identify the last date in the dataset and define that as *today's date* for all calculations.

To get started with these calculations, let's take a look at how they are performed using the built-in functionality of the lifetimes library:

%md The lifetimes library, like many Python libraries, is single-threaded.  Using this library to derive customer metrics on larger transactional datasets may overwhelm your system or simply take too long to complete. For this reason, let's examine how these metrics can be calculated using the distributed capabilities of Apache Spark.

As SQL is frequency employed for complex data manipulation, we'll start with a Spark SQL statement.  In this statement, we first assemble each customer's order history consisting of the customer's ID, the date of their first purchase (first_at), the date on which a purchase was observed (transaction_at) and the current date (using the last date in the dataset for this value).  From this history, we can count the number of repeat transaction dates (frequency), the days between the last and first transaction dates (recency), and the days between the current date and first transaction (T) on a per-customer basis:

%md Of course, Spark SQL does not require the DataFrame to be accessed exclusively using a SQL statement.  We may derive this same result using the Programmatic SQL API which may align better with some Data Scientist's preferences.  The code in the next cell is purposely assembled to mirror the structure in the previous SQL statement for the purposes of comparison:

%md Let's take a moment to compare the data in these different metrics datasets, just to confirm the results are identical.  Instead of doing this record by record, let's calculate summary statistics across each dataset to verify their consistency:

NOTE You may notice means and standard deviations vary slightly in the hundred-thousandths and millionths decimal places.  This is a result of slight differences in data types between the pandas and Spark DataFrames but do not affect our results in a meaningful way. 

%md The metrics we've calculated represent summaries of a time series of data.  To support model validation and avoid overfitting, a common pattern with time series data is to train models on an earlier portion of the time series (known as the *calibration* period) and validate against a later portion of the time series (known as the *holdout* period). In the lifetimes library, the derivation of per customer metrics using calibration and holdout periods is done through a simple method call.  Because our dataset consists of a limited range for data, we will instruct this library method to use the last 90-days of data as the holdout period.  A simple parameter called a widget on the Databricks platform has been implemented to make the configuration of this setting easily changeable:

NOTE To change the number of days in the holdout period, look for the textbox widget by scrolling to the top of your Databricks notebook after running this next cell

%md As before, we may leverage Spark SQL to derive this same information.  Again, we'll examine this through both a SQL statement and the programmatic SQL API.

To understand the SQL statement, first recognize that it's divided into two main parts.  In the first, we calculate the core metrics, *i.e.* recency, frequency and age (T), per customer for the calibration period, much like we did in the previous query example. In the second part of the query, we calculate the number of purchase dates in the holdout customer for each customer.  This value (frequency_holdout) represents the incremental value to be added to the frequency for the calibration period (frequency_cal) when we examine a customer's entire transaction history across both calibration and holdout periods.

To simplify our logic, a common table expression (CTE) named CustomerHistory is defined at the top of the query.  This query extracts the relevant dates that make up a customer's transaction history and closely mirrors the logic at the center of the last SQL statement we examined.  The only difference is that we include the number of days in the holdout period (duration_holdout):

%md And here is the equivalent Programmatic SQL API logic:

%md Using summary stats, we can again verify these different units of logic are returning the same results:

%md Our data prep is nearly done.  The last thing we need to do is exclude customers for which we have no repeat purchases, *i.e.* frequency or frequency_cal is 0. The Pareto/NBD and BG/NBD models we will use focus exclusively on performing calculations on customers with repeat transactions.  A modified BG/NBD model, *i.e.* MBG/NBD, which allows for customers with no repeat transactions is supported by the lifetimes library.  However, to stick with the two most popular of the BYTD models in use today, we will limit our data to align with their requirements:

NOTE We are showing how both the pandas and Spark DataFrames are filtered simply to be consistent with side-by-side comparisons earlier in this section of the notebook.  In a real-world implementation, you would simply choose to work with pandas or Spark DataFrames for data preparation.

%md ###Step 4: Train the Model

To ease into the training of a model, let's start with a simple exercise using a [Pareto/NBD model](https://www.jstor.org/stable/2631608), the original BTYD model. We'll use the calibration-holdout dataset constructed in the last section of this notebook, fitting the model to the calibration data and later evaluating it using the holdout data:

%md With our model now fit, let's make some predictions for the holdout period.  We'll grab the actuals for that same period to enable comparison in a subsequent step:

%md With actual and predicted values in hand, we can calculate some standard evaluation metrics.  Let's wrap those calculations in a function call to make evaluation easier in future steps:

%md While the internals of the Pareto/NBD model may be quite complex.  In a nutshell, the model calculates a double integral of two curves, one which describes the frequency of customer purchases within a population and another which describes customer survivorship following a prior purchase event. All of the calculation logic is thankfully hidden behind a simple method call.

As simple as training a model may be, we have two models that we could use here: the Pareto/NBD model and the BG/NBD model.  The [BG/NBD model](http://brucehardie.com/papers/018/fader_et_al_mksc_05.pdf) simplifies the math involved in calculating customer lifetime and is the model that popularized the BTYD approach.  Both models work off the same customer features and employ the same constraints.  (The primary difference between the two models is that the BG/NBD model maps the survivorship curve to a beta-geometric distribution instead of a Pareto distribution.) To achieve the best fit possible, it is worthwhile to compare the results of both models with our dataset.

Each model leverages an L2-norm regularization parameter which we've arbitrarily set to 0 in the previous training cycle.  In addition to exploring which model works best, we should consider which value (between 0 and 1) works best for this parameter.  This gives us a pretty broad search space to explore with some hyperparameter tuning.

To assist us with this, we will make use of [hyperopt](http://hyperopt.github.io/hyperopt/).  Hyperopt allows us to parallelize the training and evaluation of models against a hyperparameter search space.  This can be done leveraging the multiprocessor resources of a single machine or across the broader resources provided by a Spark cluster.  With each model iteration, a loss function is calculated.  Using various optimization algorithms, hyperopt navigates the search space to locate the best available combination of parameter settings to minimize the value returned by the loss function.

To make use of hyperopt, lets define our search space and re-write our model training and evaluation logic to provide a single function call which will return a loss function measure:

%md Notice that the evaluate_model function retrieves its data from a variable named inputs.  Inputs is defined in the next cell as a [broadcast variable](https://spark.apache.org/docs/latest/rdd-programming-guide.html#broadcast-variables) containing the inputs_pd DataFrame used earlier. As a broadcast variable, a complete stand-alone copy of the dataset used by the model is replicated to each worker in the Spark cluster.  This limits the amount of data that must be sent from the cluster driver to the workers with each hyperopt iteration.  For more information on this and other hyperopt best practices, please refer to [this document](https://docs.databricks.com/applications/machine-learning/automl/hyperopt/hyperopt-best-practices.html).

With everything in place, let's perform our hyperparameter tuning over 100 iterations in order to identify the best model type and L2 settings for our dataset:

%md When used with the Databricks ML runtime, the individual runs that make up the search space evaluation are tracked in a built-in repository called mlflow. For more information on how to review the models generated by hyperopt using the Databricks mlflow interface, please check out [this document](https://docs.databricks.com/applications/machine-learning/automl/hyperopt/hyperopt-spark-mlflow-integration.html).

The optimal hyperparameter settings observed during the hyperopt iterations are captured in the argmin variable.  Using the space_eval function, we can obtain a friendly representation of which settings performed best:

%md Now that we know our *best* parameter settings, let's train the model with these to enable us to perform some more in-depth model evaluation:

NOTE Because of how search spaces are *searched*, different hyperopt runs may yield slightly different results. 

%md ###Step 5: Evaluate the Model

Using a method defined in the last section of this notebook, we can calculate the MSE for our newly trained model:

%md While important for comparing models, the MSE metric is a bit more challenging to interpret in terms of the overall goodness of fit of any individual model.  To provide more insight into how well our model fits our data, let's visualize the relationships between some actual and predicted values.

To get started, we can examine how purchase frequencies in the calibration period relates to actual (frequency_holdout) and predicted (model_predictions) frequencies in the holdout period:

%md What we see here is that a higher number of purchases in the calibration period predicts a higher average number of purchases in the holdout period but the actual values diverge sharply from model predictions when we consider customers with a large number of purchases (>60) in the calibration period.  Thinking back to the charts in the data exploration section of this notebook, you might recall that there are very few customers with such a large number of purchases so that this divergence may be a result of a very limited number of instances at the higher end of the frequency range. More data may bring the predicted and actuals back together at this higher end of the curve.  If this divergence persists, it may indicate a range of customer engagement frequency above which we cannot make reliable predictions.

Using the same method call, we can visualize time since last purchase relative to the average number of purchases in the holdout period. This visualization illustrates that as time since the last purchase increases, the number of purchases in the holdout period decreases.  In otherwords, those customers we haven't seen in a while aren't likely coming back anytime soon:

NOTE As before, we will hide the code in the following cells to focus on the visualizations.  Use **Show code** to see the associated Python logic.

plot_calibration_purchases_vs_holdout_purchases(
  model, 
  input_pd, 
  kind='time_since_last_purchase', 
  n=90, 
  **{'figsize':(8,8)}
  )

display()

%md Plugging the age of the customer at the time of the last purchase into the chart shows that the timing of the last purchase in a customer's lifecycle doesn't seem to have a strong influence on the number of purchases in the holdout period until a customer becomes quite old.  This would indicate that the customers that stick around a long while are likely to be more frequently engaged:

plot_calibration_purchases_vs_holdout_purchases(
  model, 
  input_pd, 
  kind='recency_cal', 
  n=300,
  **{'figsize':(8,8)}
  )

display()

%md From a quick visual inspection, it's fair to say our model isn't perfect but there are some useful patterns that it captures. Using these patterns, we might calculate the probability a customer remains engaged:

%md The prediction of the customer's probability of being alive could be very interesting for the application of the model to our marketing processes.  But before exploring model deployment, let's take a look at how this probability changes as customers re-engage by looking at the history of a single customer with modest activity in the dataset, CustomerID 12383:

%md From this chart, we can see this customer made his or her first purchase in January 2011 followed by a repeat purchase later that month.  There was about a 1-month lull in activity during which the probability of the customer being alive declined slightly but with purchases in March, April and June of that year, the customer sent repeated signals that he or she was engaged. Since that last June purchase, the customer hasn't been seen in our transaction history, and our belief that the customer remains engaged has been dropping though as a moderate pace given the signals previously sent.

How does the model arrive at these probabilities? The exact math is tricky but by plotting the probability of being alive as a heatmap relative to frequency and recency, we can understand the probabilities assigned to the intersections of these two values:

%md In addition to predicting the probability a customer is still alive, we can calculate the number of purchases expected from a customer over a given future time interval, such as over the next 30-days: