---
title: "The CLVTools Package"
output: 
  pdf_document: 
    latex_engine: xelatex
    number_sections: yes
vignette: >
  %\VignetteIndexEntry{The CLVTools Package}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
bibliography: bibliography.bib
---


```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  #fig.path = "figures/WALKTHROUGH-",
  out.width = "100%"
)
```

# Walkthrough for the CLVTools package


## Setup the R environment

Install the stable version from CRAN:
     
```{r install-package-CRAN, eval = FALSE}
install.packages("CLVTools")
```

Install the development version from GitHub (using the `devtools` [@devtools] package):

```{r install-package-GITHUB, eval = FALSE}
install.packages("devtools")
devtools::install_github("bachmannpatrick/CLVTools", ref = "master")
```

Load the package
```{r load-library}
library("CLVTools")
```

# Load sample data provided in the package
As Input data `CLVTools` requires customers' transaction history. Every transaction record consists of a purchase date and customer ID. Optionally, the price of the transaction may be included to allow for prediction of future customer spending using an additional Gamma/Gamma model[@Fader2005b; @Colombo1999]. Using the full history of transaction data allows for comprehensive plots and summary statistics, which allow the identification of possible issues prior to model estimation. Data may be provided as `data.frame` or `data.table` [@data.table].  

It is common practice to split time series data into two parts, an estimation and a holdout period. The model is estimated based on the data from the estimation period while the data from the holdout period allows to rigorously assess model performance. Once model performance is checked on known data one can proceed to predict data without a holdout period. The length of the estimation period is heavily dependent on the characteristics of the analyzed dataset. We recommend to choose an estimation period that contains in minimum the length of the average inter-purchase time. Note that all customers in the dataset need to purchase at least once during the estimation period, i.e. these models do not account for prospects who have not yet a purchase record.

Some models included in `CLVTools` allow to model the impact of covariates. These covariates may explain heterogeneity among the customers and therefore increase the predictive accuracy of the model. At the same time we may also identify and quantify the effects of these covariates on customer purchase and customer attrition. `CLVTools` distinguishes between time-invariant and time-varying covariates. Time-invariant covariates include customer characteristics such as demographics that do not change over time. Time-varying covariates are allowed to change over time. They include for example direct marketing information or seasonal patterns.

For the following example, we use simulated data comparable to data from a retailer in the apparel industry. The dataset contains transactional detail records for every customer consisting of customer id, date of purchase and the total monetary value of the transaction.The apparel dataset is available in the `CLVTools` package. Use the `data(apparelTrans)` to load it:


```{r load-data}
data("apparelTrans")
apparelTrans
```

## Initialize the CLV-Object

Before we estimate a model, we are required to initialize a data object using the `clvdata()` command. The data object contains the prepared transactional data and is later used as input for model fitting. Make sure to store the generated object in a variable, e.g. in our example `clv.apparel`.

Through the argument `data.transactions` a `data.frame` or `data.table` which contains the transaction records, is specified. In our example this is `data.transactions=apparelTrans`. The argument `date.format` is used to indicate the format of the date variable in the data used. The date format in the apparel dataset is given as "year-month-day" (i.e., "2005-01-03"), therefore we set `date.format="ymd"`. Other combinations such as `date.format="dmy"` are possible. See the documentation of `lubridate` [@lubridate] for all details. `time.unit` is the scale used to measure time between two dates. For this dataset and in most other cases The argument `time.unit="week"` is the preferred choice. Abbreviations may be used (i.e. "w"). `estimation.split` indicates the length of the estimation period. Either the length of the estimation period (in previous specified time units) or the date at which the estimation period ends can be specified. If no value is provided, the whole dataset is used as estimation period (i.e. no holdout period). In this example, we use an estimation period of 40 weeks. Finally, the three name arguments indicate the column names for customer ID, date and price in the supplied dataset. Note that the price column is optional.


```{r load-CreateObj}
clv.apparel <- clvdata(apparelTrans,  
                       date.format="ymd", 
                       time.unit = "week",
                       estimation.split = 40,
                       name.id = "Id",
                       name.date = "Date",
                       name.price = "Price")
```


## Check the `clvdata` Object


To get details on the `clvdata` object, print it to the console.

```{r print-CLVObject}
clv.apparel
```

Alternatively the `summary()` command provides full detailed summary statistics for the provided transactional detail. `summary()` is available at any step in the process of estimating a probabilistic customer attrition model with `CLVTools`. The result output is updated accordingly and additional information is added to the summary statistics.`nobs()` extracts the number of observations.  For the this particular dataset we observe a total of 250 customers who made in total 2257 repeat purchases. Approximately 26% of the customers are zero repeaters, which means that the only a minority of the customers do not return to the store after their first purchase.


```{r summary-CLVObject}
summary(clv.apparel)
```


## Estimate Model Parameters

After initializing the object, we are able to estimate the first probabilistic latent attrition model. We start with the standard Pareto/NBD model [@Schmittlein1987] and therefore use the command `pnbd()` to fit the model and estimate model parameters. `clv.data` specifies the initialized object prepared in the last step. Optionally, starting values for the model parameters and control settings for the optimization algorithm may be provided: The argument `start.params.model` allows to assign a vector (e.g. `c(alpha=1, beta=2, s=1, beta=2)` in the case of the Pareto/NBD model) of starting values for the optimization. This is useful if prior knowledge on the parameters of the distributions are available. By default starting values are set to 1 for all parameters. The argument `optimx.args` provides an option to control settings for the optimization routine. It passes a list of arguments to the optimizer. All options known from the package `optimx` [@optimx1; @optimx2] may be used. This option enables users to specify specific optimization algorithms, set upper and/or lower limits or enable tracing information on the progress of the optimization. In the case of the standard Pareto/NBD model, `CLVTools` uses by default the optimization method `L-BFGS-G` [@byrd1995limited]. If the result of the optimization is in-feasible, the optimization automatically switches to the more robust but often slower `Nelder-Mead` method [@nelder1965simplex]. `verbose` shows additional output. 


```{r estimate-model}
est.pnbd <- pnbd(clv.data = clv.apparel)
est.pnbd
```

If we assign starting parameters and additional arguments for the optimizer we use: 
```{r estimate-model2, eval=FALSE}
est.pnbd <- pnbd(clv.data = clv.apparel, 
                     start.params.model = c(r=1, alpha = 2, s = 1, beta = 2), 
                     optimx.args = list(control=list(trace=5),
                                       method="Nelder-Mead" 
                                       ))
```

Parameter estimates may be reported by either printing the estimated object (i.e. `est.pnbd`) directly in the console or by calling `summary(est.pnbd)` to get a more detailed report including the likelihood value as well as AIC and BIC. Alternatively parameters may be directly extracted using `coef(est.pnbd)`. Also `loglik()`, `confint()` and `vcov()` are available to directly access the Loglikelihood value, confidence intervals for the parameters and to calculate the Variance-Covariance Matrix for the fitted model. For the standard Pareto/NBD model, we get 4 parameters $r, \alpha, s$ and $\beta$.  where $r,\alpha$ represent the shape and scale parameter of the gamma distribution that determines the purchase rate and $s,\beta$ of the attrition rate across individual customers. $r/\alpha$ can be interpreted as the mean purchase and $s/\beta$ as the mean attrition rate. A significance level is provided for each parameter estimates. In the case of the apparelTrans dataset we observe a an average purchase rate of $r/\alpha=0.147$ transactions and an average attrition rate of $s/\beta=0.031$ per customer per week. KKT 1 and 2 indicate the Karush-Kuhn-Tucker optimality conditions of the first and second order [@KKT]. If those criteria are not met, the optimizer has probably not arrived at an optimal solution. If this is the case it is usually a good idea to rerun the estimation using alternative starting values.


```{r param-summary}
#Full detailed summary of the parameter estimates
summary(est.pnbd)

#Extract the coefficients only
coef(est.pnbd)
#Alternative: oefficients(est.pnbd.obj)

```
To extract only the coefficients, we can use `coef()`. To access the confidence intervals for all parameters `confint()` is available.
```{r coef-model}
#Extract the coefficients only
coef(est.pnbd)
#Alternative: oefficients(est.pnbd.obj)

#Extract the confidence intervals
confint(est.pnbd)

```

In order to get the Likelihood value and the corresponding Variance-Covariance Matrix we use the following commands:
```{r ll-model}
# LogLikelihood at maximum
logLik(est.pnbd)

# Variance-Covariance Matrix at maximum
vcov(est.pnbd)

```

As an alternative to the Pareto/NBD model `CLVTools` features the BG/NBD model [@Fader2005c] and the GGomp/NBD [@Bemmaor2012a]. To use the alternative models replace `pnbd()` by the corresponding model-command. Note that he naming and number of model parameters is dependent on the model. Consult the manual for more details on the individual models.

|Command| Model | Covariates| 
|---|---|---|---|---|
| pnbd()   |Pareto/NBD|  time-invariant & time-varying  |
| bgnbd()  |BG/NBD   	| time-invariant  | 
| ggomnbd() |GGom/NBD   | time-invariant  |

To estimate the GGom/NBD model we apply the `ggomnbd()`to the `clv.apparel` object. The GGom/NBD model is more flexible than the Pareto/NBD model, however it sometimes is challenging to optimize. Note that in this particular case providing start parameters is essential to arrive at an optimal solution (i.e. `kkt1: TRUE` and `kkt2: TRUE`).
```{r estimate-ggomnbd, eval=FALSE}
est.ggomnbd <- ggomnbd(clv.data = clv.apparel, 
                     start.params.model = c(r=0.7, alpha=5, b=0.005,  s=0.02, beta=0.001), 
                     optimx.args = list(control=list(trace=5),
                                       method="Nelder-Mead"))
```

## Predicting Customer Behavior

Once the model parameters are estimated, we are able to predict future customer behavior on an individual level. To do so, we use `predict()` on the object with the estimated parameters (i.e. `est.pnbd`). The prediction period may be varied by specifying `prediction.end`. It is possible to provide either an end-date or a duration using the same time unit as specified when initializing the object (i.e  `prediction.end = "2006-05-08"` or `prediction.end = 30`). By default, the prediction is made until the end of the dataset specified in the `clvdata()` command. The argument `continuous.discount.factor` allows to adjust the discount rate used to estimated the discounted expected transactions (DERT). The default value is `0.1` (=10%). Probabilistic customer attrition model predict in general three expected characteristics for every customer:

* "conditional expected transactions" (CET), which is the number of transactions to expect form a customer during the prediction period,
* "probability of a customer being alive" (PAlive) at the end of the estimation period and
* "discounted expected residual transactions" (DERT) for every customer, which is the total number of transactions for the residual lifetime of a customer discounted to the end of the estimation period. 

If spending information was provided when initializing the `clvdata`-object, `CLVTools` provides prediction for 
 
* predicted spending estimated by a Gamma/Gamma model [@Colombo1999; @Fader2005c] and 
*  the customer lifetime value (CLV).

If a holdout period is available additionally the true numbers of transactions ("actual.x") and true spending ("actual.spending") during the holdout period are reported.

To use the parameter estimates on new data (e.g., an other customer cohort), the argument `newdata` optionally allows to provide a new `clvdata` object.

```{r predict-model}
results <- predict(est.pnbd)
print(results)

```

To change the duration of the prediction time, we use the `predicton.end` argument. We can either provide a time period (30 weeks in this example):

```{r predict-model2, eval = FALSE}
predict(est.pnbd, prediction.end = 30)
```

or provide a date indication the end of the prediction period:

```{r plot-model3, eval = FALSE}
predict(est.pnbd, prediction.end = "2006-05-08")
```


## Model Plotting
`clvdata` objects may be plotted using the `plot()` command. Similar to `summary()`, the output of `plot()` adapts to the current modeling step. It provides a descriptive plot of the actual transactional data if the model has not yet been fitted. Once the model has been estimated, `plot()` provides an aggregated incremental tracking plot of the actual data and the model based on the estimated parameters. The time-span for the plot may be altered using the `prediction.end` argument by providing either a duration or an end date. By default the plot is generated for the entire time-span of the provided dataset specified in the `clvdata()` command. The dashed line indicates the end of the estimation period. Alternatively cumulative actual and expected transactions can be plotted by setting `cumulative` to `TRUE`. The argument `transactions` disable for plotting actual transactions (`transactions=FALSE`). For further plotting options see the documentation

plot(x, prediction.end = NULL, cumulative = FALSE,
  transactions = TRUE, label = NULL, plot = TRUE, verbose = TRUE,
  ...)

```{r plot-model, fig.height=4.40, fig.width=9}
plot(est.pnbd)

```



To plot the *cumulative* expected transactions 30 time units (30 weeks in this example) ahead of the end of the estimation plot, we use:

```{r plot-model2, eval = FALSE}
plot(est.pnbd, prediction.end = 30, cumulative = TRUE)
```

Alternatively, it is possible to specify a date for the `prediction.end`argument. Note that dates are rounded to the next full time unit (i.e. week):

```{r predict-model3, eval = FALSE}
plot(est.pnbd, prediction.end = "2006-05-08", cumulative = TRUE)
```


## Covariates
`CLVTools` provides the option to include covariates into probabilistic customer attrition models. Covariates may affect the purchase or the attrition process, or both. It is also possible to include different covariates for the two processes. However, support for covariates is dependent on the model. Not all implemented models provide the option for covariates. In general, `CLVTools` distinguishes between two types of covariates: time-invariant and time-varying. The former include factors that do not change over time such as customer demographics or customer acquisition information. The latter may change over time and include marketing activities or seasonal patterns.

Data for time-invariant covariates must contain a unique customer ID and a single value for each covariate. It should be supplied as a `data.frame` or `data.table`. In the example of the apparel retailer we use demographic information "gender" as time-invariant and information on the acquisition channel as covariate for both, the purchase and the attrition process. Use the `data("apparelStaticCov")` command to load the time-invariant covariates. In this example gender is coded as a dummy variable with `male=0` and `female=1` and channel with `online=0` and `offline=1`.

```{r Cov-staticData}
data("apparelStaticCov")
apparelStaticCov
```

Data for time-varying covariates requires a time-series of covariate values for every customer. I.e. if the time-varying covariates are allowed to change every week, a value for every customer for every week is required. Note that all contextual factors are required to use the same time intervals for the time-series. In the example of the apparel retailer we use information on direct marketing (`Marekting`) as time-varying covariate. Additionally, we add gender as time-invariant contextual factors. Note that the data structure of invariant covariates needs to be aligned with the structure of time-varying covariate. Use `data("apparelDynCov")` command to load 

```{r Cov-dynData}
data("apparelDynCov")
apparelDynCov
```

To add the covariates to an initialized `clvdata` object the commands `SetStaticCovariates()` and `SetDynamicCovariates()` are available. The two commands are mutually exclusive. The argument `clv.data` specifies the initialized object and the argument `data.cov.life` respectively `data.cov.trans` specifies the data source for the covariates for the attrition and the purchase process. Covariates are added separately for the purchase and the attrition process. Therefore if a covariate should affect both processes it has to be added in both arguments: `data.cov.life` *and* `data.cov.trans`. The arguments `names.cov.life` and `names.cov.trans` specify the column names of the covariates for the two processes. In our example, we use the same covariates for both processes. Accordingly, we specify the time-invariant covariates "Gender" and "Channel" as follows:

```{r Cov-setStatic}
clv.static<- SetStaticCovariates(clv.data = clv.apparel, 
                                      data.cov.life = apparelStaticCov, 
                                      data.cov.trans = apparelStaticCov,
                                      names.cov.life = c("Gender", "Channel"), 
                                      names.cov.trans =c("Gender", "Channel"), 
                                      name.id = "Id")
```

To specify the time-varying contextual factors for seasonal patterns and direct marketing, we use the following:

```{r Cov-setDynamic, eval=FALSE, message=FALSE, warning=TRUE}
clv.dyn <- SetDynamicCovariates(clv.data = clv.apparel, 
                                     data.cov.life = apparelDynCov,
                                     data.cov.trans = apparelDynCov, 
                                     names.cov.life = c("Marketing", "Gender", "Channel"), 
                                     names.cov.trans = c("Marketing", "Gender", "Channel"), 
                                     name.id = "Id",
                                     name.date = "Cov.Date")
```

In order to include time-invariant covariates in a time-varying model, they may be recoded as a time-varying covariate with a constant value in every time period.

Once the covariates are added to the model the estimation process is almost identical to the standard model without covariates. The only difference is that the provided object now data for contains either time-invariant or time-varying covariates and the option to define start parameters for the covariates of both processes using the arguments `start.params.life` and `start.params.trans`. If not set, the staring values are set to 1. To define starting parameters for the covariates, the name of the corresponding factor has to be used. For example in the case of time-invariant covariates:

```{r Static-cov-estimate, message=TRUE, warning=FALSE}
est.pnbd.static <- pnbd(clv.static, 
                         start.params.model = c(r=1, alpha = 2, s = 1, beta = 2),
                         start.params.life = c(Gender=0.6, Channel=0.4),
                         start.params.trans = c(Gender=0.6, Channel=0.4))
```

Analogously, we can estimate the model containing time-varying covariates. We recommend to enable the built-in support for multithreading when estimating more complex models like this. See section [Multithreading](#multithreading).

```{r Dyn-cov-estimate, eval=FALSE}
est.pnbd.dyn <- pnbd(clv.dyn, 
                     start.params.model = c(r=1, alpha = 2, s = 1, beta = 2),
                     start.params.life = c(Marketing=0.5, Gender=0.6, Channel=0.4),
                     start.params.trans = c(Marketing=0.5, Gender=0.6, Channel=0.4))
```

To inspect the estimated model we use `summary()`, however all other commands such as `print()`, `coef()`, `loglike()`, `confint()` and `vcov()` are also available. Now, output contains also parameters for the covariates for both processes. Since covariates are added separately for the purchase and the attrition process, there are also separate model parameters for the two processes. These parameters are directly interpretable as rate elasticity of the corresponding factors: A 1% change in a contextual factor $\bf{X}^{P}$ or $\bf{X}^{L}$ changes the purchase or the attrition rate by $\gamma_{purch}\bf{X}^{P}$ or $\gamma_{life}\bf{X}^{L}$ percent, respectively  [@Gupta1991]. In the example of the apparel retailer, we observe that female customer purchase significantly more (`trans.Gender=1.42576`). Note, that female customers are coded as 1, male customers as 0. Also customers acquired offline (coded as Channel=1), purchase more (`trans.Channel=0.40304`) and stay longer (`life.Channel=0.9343`). Make sure to check the Karush-Kuhn-Tucker optimality conditions of the first and second order [@KKT] (KKT1 and KKT1) before interpreting the parameters. If those criteria are not met, the optimizer has probably not arrived at an optimal solution. If this is the case it is usually a good idea to rerun the estimation using alternative starting values.


```{r Cov-summary}
summary(est.pnbd.static)
```


To predict future customer behavior we use `predict()`.  Note that dependent on the model, the predicted metrics may differ. For example, in the case of the Pareto/NBD model with time-varying covariates, instead of DERT, DECT is predicted. DECT only covers a finite time horizon in contrast to DERT. Time-varying covariates must be provided for the entire prediction period. If the data initially provided in the `SetDynamicCovariates()` command does not cover the complete prediction period, the argument `new.data` offers the ability to supply new data for the time-varying covariates in the from of a `clvdata` object.

## Add Correlation to the model
To relax the assumption of independence between the purchase and the attrition process, `CLVTools` provides the option to specify the argument `use.cor` in the command t fit the model (i.e. `pnbd`). In case of `use.cor=TRUE`, a Sarmanov approach is used to correlate the two processes. `start.param.cor` allows to optionally specify a starting value for the correlation parameter.

```{r Cov-cor, eval=FALSE}
est.pnbd.cor <- pnbd(clv.apparel, 
                     use.cor= TRUE)
summary(est.pnbd.cor)
```

The parameter `Cor(life,trans)` is added to the parameter estimates that may be directly interpreted as a correlation. In the example of the apparel retailer the correlation parameter is not significant and the correlation is very close to zero, indicating that the purchase and the attrition process are independent.

## Advanced Options for Contextual Factors
`CLVTools` provides two additional estimation options for models containing covariates (time-invariant or time-varying): regularization and constraints for the parameters of the covariates.  Both options are included in the command to fit the model (i.e., `pnbd()`. Support for this option is dependent on the model. They may be used simultaneously.

* The argument `reg.lambdas` provides the possibility to specify separate `\lambda_{reg}` for the two processes (i.e. `reg.lambdas = c(trans=100, life=100)`. The larger the `\lambda_{reg}` the stronger the effects of the regularization. Regularization only affects the parameters of the covariates.
* The argument `names.cov.constr` implements equality constraints for contextual factors with regards to the two processes. For example the variable "gender" is forced to have the same effect on the purchase as well as on the attrition process. To do so, the option `names.cov.constr` is available (i.e. `names.cov.constr=c("Gender")`). To provide starting parameters for the constrained variable use `start.params.constr`.

To enable regularization for the covariates, we use the following command:

```{r reg-advOptions,}
est.pnbd.reg <- pnbd(clv.static, 
                         start.params.model = c(r=1, alpha = 2, s = 1, beta = 2),
                         reg.lambdas = c(trans=100, life=100))
summary(est.pnbd.reg)
```
To constrain "Gender" to have the same effect on both processes we use the following command. Note, that the output now only contains one parameter for "Gender" as it is constrained to be the same for both processes.

```{r constr-advOptions,}
est.pnbd.constr <- pnbd(clv.static, 
                         start.params.model = c(r=1, alpha = 2, s = 1, beta = 2),
                         start.params.constr = c(Gender=0.6),
                         names.cov.constr=c("Gender"))
summary(est.pnbd.constr)
```

## Multithreading: Enable parallel processing for CLVTools {#multithreading}
`CLVTools` supports parallel processing when estimating models containing time-varying covariates using the package `future` [@future] and the corresponding `foreach` [@foreach] parallel adapter `doFuture` [@doFuture]. To enable the distribution of the workload across all available cores, use execute the following commands before estimating the model (i.e. calling `pnbd()):

```{r doFuture, eval=FALSE}
  # disable multithreading for data.table (to avoid nested parallelism)
  setDTthreads(1)

  library("doFuture")
  registerDoFuture()
  plan("multisession")
```
To limit the number of parallel processes simultaneously executed, we can specify the number of workers `plan(multisession, workers = 2)`.




### Literature


