SMCRM_Ch04.Rmd

---
title: "Statistical Methods in Customer Relationship Management"
subtitle: "Chapter 4. Customer Retention"
author: "Alexander Rodionov"
date: "6 июня 2018 г."
output: 
  html_document:
    highlight: tango 
    theme: readable
    code_folding: hide
    fig_caption: yes
    fig_height: 4
    fig_width: 6.3
    toc: yes
    toc_float: no
    number_sections: no
---

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

Мы продолжаем рассматриваем в **R** примеры **SAS** на наборах данных из книги Kumar V., Petersen Andrew J. "Statistical Methods in Customer Relationship Management".

# Глава 4. Удержание клиентов

После того, как клиенты были приобретены, важно обсудить второй шаг CRM - *сохранить клиентов*. Стратегии удержания клиентов используются как в договорных условиях (когда клиенты привязаны к контрактам, таким как догвор на мобильный телефон с пакетом связи или подписка на журнал) и внедоговорные условия (где клиенты не связаны контрактами, такими как покупки в продуктах или покупки одежды). Рейхельд и Сассер (Reichheld and Sasser, 1990) заявили, что 5% -ное улучшение удержания клиентов может привести к увеличению рентабельности в пределах от 25 до 85% с точки зрения чистой текущей стоимости в зависимости от отрасли. С тех пор компании постоянно выделяют ресурсы на управление удержанием клиентов, и исследователи уделяют большое внимание изучению сохранению клиентов.

Исследования по сохранению клиентов имеют два основных направления. Одно направление заинтересовано в исследовании влияния различных маркетинговых переменных на удержание клиентов, что, в свою очередь, влияет на производительность фирмы. Другое напарвление заинтересовано в создании эконометрических и статистических моделей для оценки или прогнозирования решений о сохранении клиентов как со стороны клиента, так и со стороны компании.

```{r Picture 1, echo=FALSE, out.width=1155, fig.cap = 'The Picture 1'}

knitr::include_graphics("C:/Soft/R/Examples/CRM/SMCRM_Ch04-1.png")

```

На рисунке 1 показана интегрированная структура, которая описывает различные отношения, рассмотренные в разных исследованиях во многих отраслях. Эти отрасли включают в себя, в частности, телекоммуникации, финансовые услуги, сферу услуг, рестораны и розничную торговлю. Общей темой этих отношений является мышление, которое направлено на:

* a) повышение качества продукции и услуг, приводящее к повышению удовлетворенности клиентов; 

* б) всемерная удовлетворенность клиентов приводит к увеличению удержания клиентов (что обусловлено качеством отношений, так что более высокое качество отношений положительно улучшает связь между удовлетворением и удержанием) и 

* c) увеличение удержания клиентов приводит к повышению эффективности работы фирмы.

```{r Picture 2, echo=FALSE, out.width=1155, fig.cap = 'The Picture 2'}

knitr::include_graphics("C:/Soft/R/Examples/CRM/SMCRM_Ch04-2.png")

```

Авторам книги по понятным причинам импонирует второе направление исследований по удержанию клиентов (см. рисунок 2), которое влияет на решения, принимаемых менеджерами по текущим клиентам. Часто возникает несколько ключевых вопросов, в разрешении которых менеджеры заинтересованы при формировании ответа после привлечения клиента. К ним относятся:

* Будет ли недавно приобретенный покупатель купать или нет в дальнейшем?

* Какова будет продолжительность жизни клиента (т. е. когда клиент завершит сотрудничество с фирмой)?

* Учитывая, что клиент собирается выкупить:

    + Сколько предметов покупает клиент?
    
    + Сколько этот клиент способен потратить?
    
    + Будет ли этот клиент покупать / заказывать в нескольких категориях продуктов?

* Покупает ли потребитель в основном у одной фирмы (значительная доля в кошельке) или у многих разных фирм (низкая доля в кошельке)?

* Какое долгосрочное влияние покупательского поведения клиента на стоимость фирмы?

Чтобы дать полное понимание того, как моделировать процесс удержание клиентов, авторы рассматривают вышеперечисленные вопросы в исследованиях один за другим вместе с соответствующими методами моделирования. Они также представляют эмпирические примеры в конце каждого подраздела, чтобы продемонстрировать, как применять эти знания к представительной выборке клиентов из фирмы **B2C**.

***

Подобно моделям приобретения клиентов, первый вопрос, на который необходимо ответить при выборе модели, заключается в том, встпуют ли клиенты договорные и внедоговорные отношения с фирмой. В большинстве случаев это определяет тип статистической модели, который необходимо использовать для получения информации из данных.

## Данные для эмпирических примеров

В этой главе авторы дают описание основных этапов моделирования, в ходе которого пытаются ответить на каждый ключевой вопрос исследования, поднятый в начале главы. Они также предоставляют по крайней мере один эмпирический пример в конце каждого подраздела, который покажет, какие данные могут использоваться для формирования ответа на эти ключевые вопросы исследования. Для всех эмпирических примеров в этой главе авторы предлагают набор данных под названием «Удержание клиентов», который разбит на две связанные таблицы данных. В этом наборе данных вы найдете две таблицы данных, которые включают репрезентативную выборку из 500 клиентов из типичной фирмы **B2C**, где все клиенты принадлежат к одной когорте. В этом случае когорта состоит из случайной выборки из 500 клиентов, которые совершили первую покупку у фирмы в четвертом квартале. В первой таблице данных предоставлена информацию о транзакциях для каждого клиента в течение 12 кварталов ("`customerRetentionTransactions`"). Таким образом, таблица данных состоит из 6000 строк (500 клиентов * 12 кварталов) и 8 столбцов. Во второй таблице данных предоставлена демографическая информация по каждому клиенту ("`customerRetentionDemographics`"). Таким образом, таблица данных состоит из 500 строк (500 клиентов) и 6 столбцов.

Первая таблица данных ("`customerRetentionTransactions`") включает нижеследующие переменные, которые будут использоваться в некоторой комбинации в ходе каждого последующего этапа анализа:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| Customer        | Номер клиента клиента (от 1 до 500) |
| Quarter         | Квартал (от 1 до 12), когда произошла транзакция |
| Purchase        | если **1**, когда покупатель приобрел в данном квартале, **0**, если в этом квартале не произошло покупок |
| Order_quantity  | Долларовая стоимость покупок в данном квартале |
| Crossbuy        | Количество различных категорий товаров / услуг, приобретенных в данном квартале |
| Ret_Expense     | Доллары потраченные на маркетинговые усилия по удержанию этого клиента в данном квартале |
| Ret_Expense_SQ  | Квадрат затрат на маркетинговые усилияпо удержанию этого клиента в данном квартале |

Вторая таблица данных ("`customerRetentionDemographics`") включает следующие переменные
которые будут использоваться в некоторой комбинации в ходе каждого последующего этапа анализа:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| Customer        | Номер клиента клиента (от 1 до 500) |
| Gender          | **1**, если клиент является мужчиной, **0**, если клиент является женщиной |
| Married         | **1**, если клиент женат, **0**, если клиент не состоял в браке |
| Income          | **1**, если доход < 30 000 долл. США, **2**; если 30 001 долл. США < доход < 45 000 долл. США; **3**, если 45 001 долл. США < доход < 60 000 долл. США; **4**, если 60 001 долл. США < доход < 75 000 долл. США; **5**, если 75 001 долл. США < доход < 90 000 долл. США; **6**, если доход > 90 001 долл. США |
| First_Purchase  | Стоимость первой покупки, сделанной клиентом в 1 квартале |
| Loyalty         | **1**, если клиент является членом программы лояльности, **0**, если нет |
| Share-of-Wallet (SOW) | Процент покупок клиента от данной фирмы, учитывая общий объем покупок по всем фирмам в этой категории |
| CLV             | Дисконтированная стоимость всех ожидаемых будущих прибылей или стоимости жизни клиента|

 Эти примеры будут охватывать темы "повторная покупка или нет"", "количество заказов", "размер заказа", "перекрестные покупки", "доля в кошельке (SOW)" и "доходность (CLV)".

```{r Ch04: Customer retention - Data}
library('tidyverse')
library('caret')
library('car')
utils::data('customerRetentionTransactions', package = 'SMCRM')
utils::data('customerRetentionDemographics', package = 'SMCRM')

```

Странно, но в авторском наборе данных отсутствует важнейший демографический показатель **возраст**.

## Эмпирический пример: Повторная покупка или нет (остаться или уйти)

Один из ключевых вопросов, на который авторы хотят ответить по удержанию клиентов, заключается в том, можем ли мы определить, какие клиенты имеют наибольшую вероятность повторной покупки ("`repurchase`"). Для этого сначала нужно узнать, какие текущие клиенты фактически совершили дополнительные покупки после их первоначальной первой покупки. В наборе данных, предоставленном для этой главы, имеется бинарная переменная, которая идентифицирует, покупает ли покупатель за данный период времени, в этом случае - квартале. Авторы также предоставляютм набор предикторов, которые могут помочь объяснить решение клиента о повторной покупке. В конце этого примера вы сможете сделать следующее:

1. Определите драйверы поведения клиентов для повторной покупки.

2. Интерпретируйте оценки параметров из модели повторной покупки.

3. Предскажите количество повторных покупок клиентами.

4. Определить предиктивную точность модели повторной покупки.

Авторы в обзоре литературы указали, что ряд исследователей включают в число предикторов в подобные модели показатели **RFM** (англ. **Recency** – How recently did the customer purchase? **Frequency** – How often do they purchase? **Monetary Value** – How much do they spend?) - CRM метрики, которая позволяет детально описать в трех переменных ранжированную совокупность клиентов. Однако сами авторы книги решили не прибегать к ее применению.

Компания **B2C** хочет увеличить доля клиентов, совершивших повторную покупку и сократить расходы на удержание клиентов, лучше понимая при этом, какие клиенты чаще всего покупают повторно за определенный период времени. Случайная выборка из 500 клиентов из одной когорты была взята из базы данных фирмы. Информация, необходимая для нашей модели, включает следующий список переменных:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| **Зависимая переменная** |  |
| Purchase        | если **1**, когда покупатель приобрел в данном квартале, **0**, если в этом квартале не произошло покупок |
| **Предикторы**          |  |
| Lag_Purchase    | **1**, если клиент приобрел в предыдущем квартале, **0**, если в предыдущем квартале покупка не произошла |
| Avg_Order_Quantity | Средняя долларовая стоимость покупок во всех предыдущих кварталах |
| Ret_Expense     | Доллары потраченные на маркетинговые усилия по удержанию этого клиента в данном квартале |
| Ret_Expense_SQ  | Квадрат затрат на маркетинговые усилияпо удержанию этого клиента в данном квартале |
| Gender          | **1**, если клиент является мужчиной, **0**, если клиент является женщиной |
| Married         | **1**, если клиент женат, **0**, если клиент не состоял в браке |
| Income          | **1**, если доход < 30 000 долл. США, **2**; если 30 001 долл. США < доход < 45 000 долл. США; **3**, если 45 001 долл. США < доход < 60 000 долл. США; **4**, если 60 001 долл. США < доход < 75 000 долл. США; **5**, если 75 001 долл. США < доход < 90 000 долл. США; **6**, если доход > 90 001 долл. США |
| First_Purchase  | Стоимость первой покупки, сделанной клиентом в 1 квартале |
| Loyalty         | **1**, если клиент является членом программы лояльности, **0**, если нет |

```{r Ch04 : Customer retention - Transformation data}
# Combine Data
library('sqldf') # Manipulate R Data Frames Using SQL

customerRetention <- sqldf("select distinct a.*, b.* 
                           from customerRetentionTransactions a left join customerRetentionDemographics
                           b on a.customer = b.customer order by customer;")
# Repurchase Probability
customerRetentionRepurchase <- customerRetention %>% 
  mutate(lcustomer = lag(customer, 1), cb = if_else(crossbuy > 1, 1, 0)) %>%  # Lag by customer
    mutate(lpurchase = if_else(customer == lcustomer, lag(purchase), 0),
           lcb = if_else(customer == lcustomer, lag(cb), 0),
           lcrossbuy = if_else(customer == lcustomer, lag(crossbuy), 0)) %>% 
      mutate(lpurchase = if_else((quarter == 2), 1, lpurchase), 
             lcrossbuy = if_else((quarter == 2), 1, lcrossbuy),
             lcb =       if_else((quarter == 2), 0, lcb)) %>% 
        group_by(customer) %>% 
          mutate(quantity_sum = cumsum(order_quantity)) %>% 
            mutate(avg_order_quantity = quantity_sum / quarter) %>%
              ungroup %>% 
                mutate(lavg_order_quantity = lag(avg_order_quantity, 1)) %>% # should `avg_order_quantity` by lag 1
                  mutate(ID = row_number()) %>% 
                    arrange(desc(ID)) %>%
                      select(-one_of("ID", "quantity_sum", "customer..8")) %>%
                        filter(quarter != 1)

# Check for Class Imbalances
customerRetentionRepurchase$purchase %>%
  factor() %>%
    table() %>%
      prop.test()

# ## Open workbook into temporary file
# openxlsx::addWorksheet(wb0 <- openxlsx::createWorkbook(), sheetName = "Output", gridLines = FALSE)
# openxlsx::writeData(wb0, sheet = 1, x = customerRetentionRepurchase, withFilter = TRUE); openxlsx::openXL(wb0)

```

Бинарный класс *"Повторной покупки"* получился сбалансированный - нулевую гипотезу о равном разбиении после теста пропорций долей следует признать верной.

Нужно смоделировать вероятность того, что клиент купит за данный период времени. Поскольку зависимая переменная (`purchase`) является бинарной, выбирают *логическую* регрессию для оценки модели. Также можно выбрать модель *пробита* и в целом достичь тех же результатов. В этом случае зависимая переменная представляет собой `Purchase`, а предикторы представляют девять независимых переменных.

### Построение и верификация модели *повторной покупки* {.tabset}

#### Logit

```{r Ch04 : Repurchase - Logit, warning=FALSE}
# # Fit Logistic Regression Model for Customer Retention by Authors

Ch04.logit <- glm(purchase ~ lpurchase + avg_order_quantity + ret_expense + ret_expense_sq + 
                  gender + married + income + first_purchase + loyalty,
                  # lavg_order_quantity + factor(income) + factor(loyalty),
                  data = customerRetentionRepurchase, family = binomial(link = 'logit'))
summary(Ch04.logit)
writeLines(sprintf("-2 Log L of Intercept and Only Covariates: %.3f", -2 * logLik(Ch04.logit)[1]))
writeLines(sprintf("                  AIC (smaller is better): %.3f", extractAIC(Ch04.logit)[2]))

# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")
car::Confint(Ch04.logit) %>%
  exp() %>%
    arm::pfround(digits = 3)

writeLines("\n Wald test of predictors")
car::Anova(Ch04.logit, type="II", test="Wald")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04.logit) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity

writeLines("\n Repurchase Probability (Logit): Association of Predicted Probabilities and Observed Responses \n")
prob <- predict(Ch04.logit, newdata = customerRetentionRepurchase, type = "response")
caret::confusionMatrix(data = ifelse(prob > 0.5, "1", "0") %>% factor,
                       reference = customerRetentionRepurchase$purchase %>% factor,
                       positive = "1", mode = "everything")
```

Во-первых, это означает, что `Lag_Purchase` оказывает положительное влияние на текущую покупку, то есть клиенты, совершившие покупку в предыдущем квартале, с большей вероятностью совершают покупку в текущем квартале. Во-вторых, поскольку коэффициент по `Avg_Order_Quantity является` положительным и статистически значимым, это означает, что клиенты, которые в прошлом потратили больше в среднем, также чаще покупают в текущем периоде времени. В-третьих, авторы обнаружилм положительный, но уменьшающий доход от эффекта удерживающих расходов (Ret_Expense) при покупке в том же квартале, поскольку коэффициент на Ret_Expense положителен, а коэффициент Ret_Expense_SQ отрицателен. В-четвертых, нашли небольшой положительный эффект для женщин (отрицательный коэффициент по полу), что означает, что женщины, как правило, чаще приобретают, чем мужчины. В-пятых, обнаружен положительный эффект дохода, предполагающий, что клиенты, имеющие более высокий доход, с большей вероятностью будут покупать в текущем квартале. Наконец, поскольку коэффициент при лояльности положителен, это говорит о том, что клиенты, которые являются членами программы лояльности, с большей вероятностью будут покупать в текущем квартале.

Поскольку проблема *несбалансированности классов* практически отсутствует, то в качестве меры точности модели можно рекомендовать коэффициента **Cohen's Kappa** и **коэффициент аккуратности**.

#### Logit (Version AR)

Следует заметить, что расчет среднеквартального объема заказа `avg_order_quantity` при подсчете по приведенному в книге коду **SAS** фактически включает также *текущий квартал*. Хотя в описании авторы указали, что среднеквартальных объем заказа относиться только в предыдущим месяцам. В результате возникает смещение, учитывающее влияние объема заказа текущего квартала, которые и следует предсказывать, что особенно заметно проявляется в первых кварталах наблюдений. На мой взгляд было бы более точным брать предиктором `avg_order_quantity` с однопериодным лагом (кварталом).

Замечено, что если предикторы `income` и `loyalty` объявить порядковым и бинарным факторами, соответствено, а незначимые признаки `married`, `first_purchase` вовсе удалить, то и с меньшим числом предикторов можно получить устойчивую логистическую модель, но без смещения средних продаж на текущий квартал. Кроме того, в связи со значительной **мультиколлинеарностью** показателей затрат на удержание в текущем квартале и их квадрата - `ret_expense` и `ret_expense_sq` следует исключить.

```{r Ch04 : Repurchase Logit Improved, warning=FALSE}
# Fit Logit Regression Model for Customer Repurchase (Improved)
uno = 'logit'
set.seed(2018) #From random.org
(Ch04lg.AR <- train(factor(purchase) ~ lpurchase + lavg_order_quantity + ret_expense +
                  ordered(income) + first_purchase + factor(loyalty), metric = 'Kappa',
                  data = customerRetentionRepurchase, method = "glm", family = binomial(link = uno)))#,
                  # trControl = trainControl(method = "none", number = 1)))
summary(Ch04lg.AR)
# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")
car::Confint(Ch04lg.AR$finalModel) %>%
  exp() %>%
    arm::pfround(digits = 3)

# Logistic regression diagnostics
writeLines("\n Wald test of predictors")
car::Anova(Ch04lg.AR$finalModel, type="II", test="Wald")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04lg.AR$finalModel) # Variance Inflation Factors - if vif() > 2  - feature has multicollinearity
writeLines("\n Variable Importance for Model \n")
caret::varImp(Ch04lg.AR) %>% .$importance %>% print.AsIs()

# Create the scatter plots Logit versus model predictors
# prob <- predict(Mod04lg.AR, newdata = customerRetentionRepurchase, type = "prob")[, "1"]
#     mutate(logit = log(prob / (1 - prob))) %>%
set.seed(2018)
Ch04LG.AR <- train(factor(purchase) ~ lpurchase + lavg_order_quantity + ret_expense + 
                  income + first_purchase + loyalty, metric = 'Kappa',
                  data = customerRetentionRepurchase, method = "glm", family = binomial(link = uno))
predictors <- Ch04LG.AR$coefnames
Ch04LG.AR$trainingData %>%
  dplyr::select(one_of(predictors)) %>%
    mutate(link = predict(Ch04LG.AR$finalModel, newdata = customerRetentionRepurchase, type = "link")) %>%
      gather(key = "predictors", value = "predictor.value", -link) %>%
        ggplot(aes(predictor.value, link))+
          geom_point(size = 0.05, alpha = 0.05) +
          geom_smooth(method = "loess") +
          facet_wrap(~ predictors, scales = "free_x") +
          ylab(stringr::str_to_title(uno))
remove(Ch04LG.AR)

# Plot matrix of statistical model diagnostics
GGally::ggnostic(Ch04lg.AR$finalModel, title = paste(paste(formula(Ch04lg.AR)[c(2, 1, 3)], collapse = " ")))

# wide variety of diagnostic plots for checking the quality of regression fit
# https://bookdown.org/jefftemplewebb/IS-6489/logistic-regression.html
car::influenceIndexPlot(Ch04lg.AR$finalModel)

writeLines("\n Improved Logit Model: Association of Predicted Probabilities and Observed Responses \n")
caret::confusionMatrix(data = predict(Ch04lg.AR, newdata = customerRetentionRepurchase),
                       reference = customerRetentionRepurchase$purchase %>% factor,
                       positive = "1", mode = "everything")

qplot(`Observed Classes`, `Predicted Classes`, 
      data=bind_cols(`Observed Classes`= factor(customerRetentionRepurchase$purchase),
                     `Predicted Classes` = predict(Ch04lg.AR, newdata = customerRetentionRepurchase)),  
      colour= `Observed Classes`, geom = c("boxplot", "jitter"),
      main = "Predicted Classes vs. Observed Classes", xlab = "Observed Classes", ylab = "Predicted Classes")
```

Но авторы на такие шаги не пошли, вероятно, для большей дидактической ясности. Правда, аккуратность модели с лагом, т. е. `**l**avg_order_quantity` будет немного ниже, но зато устойчивость (модель получена бутстреп-методом - англ. `bootstrap` с 25-тью "псевдовыборок") станет существенной. 

### Как это использовать?

Авторы полагали, что поведение при проведении транзакций в прошлом, скорее всего, объяснит будущее поведение покупателей. В результате в этом примере они использовали несколько операций из предыдущего квартала по отношению текущему в качестве независимых переменных (предикторов). В дальнейшем менеджерам будет полезно знать какие  драйверы влияют на поведения последующих покупок клиентами. 

Во-первых, авторы знают, приобрел ли покупатель в последнем квартале (`Lag_Purchase`). Эта переменная может быть получена путем принятия запаздывающего значения переменной индикатора покупки, отмечая, что одно наблюдение будет потеряно для каждого клиента. В этом случае мы используем только однопериодное отставание - квартальное. Во-вторых, мы имеем среднюю величину прошлых заказов (`Avg_Order_Quantity`). В этом случае значение для среднего количества заказа является средним значением переменной `Order_Quantity` во всех кварталах до текущего периода времени (как я показал выше, фактически это нет и текущий период авторами включается в расчет средней). В-третьих, у авторам известно, сколько долларов фирма потратила на каждого клиента (`Ret_Expense`) за каждый период времени и квадрат значения этой переменной (`Ret_Expense_SQ`). Авторы хотили использовать как линейные, так и квадратичные термины, так как ожидали, что для каждого дополнительного доллара, потраченного на усилия по удержанию для данного клиента, будет уменьшаться возврат к стоимости этого доллара (*"закон убывающей доходности"*). Наконец, поскольку фирма этого примера относится к сектору **B2C**, остальные пять переменных являются социально-демографическими характеристиками клиентов. К ним относятся пол клиента `gender`, является ли клиент женатым (`married`), порядковый ранг дохода клиента (`income`), стоимость первой покупки клиента (`First_Purchase)` и является ли клиент членом программы *лояльности* (`loyalty`).

В результате мы теперь знаем, как изменения в расходах на удержание, прошлые транзакции клиентов и характеристики клиентов могут либо увеличить, либо уменьшить вероятность последующих покупок. И мы также знаем, что эти драйверы споосбы быть полезными, помогая нам предсказать, собирается ли клиент купать или нет. Эта информация может дать значительную информацию руководителям, которым поручено определить оптимальный объем ресурсов для проведения усилий по удержанию.

## Эмпирический пример: Продолжительность сотрудничества с клиентом

Один из ключевых вопросов, на которые авторы хотят ответить в отношении продолжительности жизни, заключается в том, можем ли мы определить, какие клиенты имеют наибольшую вероятность быть активными в будущем. В условиях без контракта это означает оценку вероятности того, что клиент в настоящее время активен, учитывая его прошлую историю покупок. В случае контрактных установок это означает оценку ожидаемого срока службы клиентов, которые еще не получили недостатков, учитывая историческую информацию обо всех клиентах в прошлом (в том числе тех, кто уже отказался от сотрудничества). Это часто делается с использованием моделей Времени Ускореннного Отказа (англ. `Accelerated Failure Time (AFT) Models`) или Пропорциональной Опасности (англ. `Proportional Hazards (PH) models`). Поскольку данные, приведенные в этой главе, представляют собой неконтрактную установку (т. е. авторы не наблюдали непосредственно завершение сотрудничества с клиентами), их цель - определить вероятность того, что клиент активен, учитывая прошлую историю покупок клиента. Для этого нам необходимо иметь информацию о поведении транзакций каждого клиента, включая время первой покупки, время последней покупки и количество транзакций, которые произошли во время окна наблюдения. В наборе данных, приведенном в этой главе, авторы подробно описывали историю транзакций для каждого клиента. Просто нужно вычислить значения для каждой из трех требуемых переменных. В ходе этого примера можно определить следующее:

1. Вероятность того, что клиент активен в конце конца наблюдений (12-го квартала).

Компания B2C хочет улучшить свою способность идентифицировать клиентов, которые, вероятно, будут активно участвовать в отношениях с фирмой. Случайная выборка из 500 клиентов из одной когорты была взята из базы данных клиентов. Информация, необходимая для модели авторов, включает следующий список переменных:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| x               | Количество транзакций данного клиента за все периоды времени. Здесь мы предполагаем, что это сумма переменной Purchase, где клиенты в наибольшей степени совершили 1 покупку за квартал |
| tx              | Это время последней транзакции, то есть последний квартал, где наблюдалось `Purchase == 1` |
| T               | Общее время между первой покупкой и окончанием окна наблюдения, то есть 12 кварталов для всех клиентов |

В этом случае авторы не имели зависимую переменную, так как фактически не наблюдали отток клиентов. Вместо этого они использовали атрибуты информации о транзакции клиента для формирования вероятности того, что клиент активен. Во первых, в этом случае нам требуется количество транзакций, которые клиент имел в окне наблюдения (`x`), в данном примере - 12 кварталов. Чтобы упростить этот случай, предполагается, что клиенты покупают не более одного раза в любом квартале. Таким образом, при `Purchase == 1` наблюдается как бы одну транзакция. Во-вторых, требуется, чтобы последний раз, когда у клиента была транзакция с фирмой (`tx`). В этом случае это последний квартал, где наблюдалось `Purchase == 1`. Наконец, также требуется продолжительность времени, когда клиент мог быть активным. Поскольку это когорта клиентов, которые сделали первые покупки в первом квартале, все они получают значение **12** для `T`.

Так как в сегменте **B2C** фирмы зачастую не имеют долгосрочных контрактов, то для построения модели продолжительности сотрудничества с клиентом без контрактных установок применяют подход BG/NBD (от англ. `"Beta Gamma / Negative Binomial Distribution"`), [описанный](http://brucehardie.com/papers/018/fader_et_al_mksc_05.pdf) в 2005 г.:

$$ L(r, \alpha, a, b|X=x, t_x, T) = \dfrac{B(a, b + x + 1)}{B(a, b)} \dfrac{\Gamma(r + x) \alpha ^r }{\Gamma(r)(\alpha + T) ^{r + x} } + \dfrac{B(a + 1, x)}{B(a, b)} \dfrac{\Gamma(r + x) \alpha ^r }{\Gamma(r)(\alpha + t_x) ^{r + x}} , $$

где B(˚) - функция бета-распределения,

G(˚) - функция гамма-распределения,

a и b - параметры функции бета-распределения, 

r и $\alpha$ - параметры функции гамма-распределения, и

x, tx, и T - данные по клиентам без контрактных установок от фирмы.

Описание в *MS Excel* данного подхода к построению модели BG/NBD, благодаря Брюсу Харди, имеется на странице [Implementing the BG/NBD Model for Customer Base Analysis in Excel](http://www.brucehardie.com/notes/004/)

Несложно заметить, что показатель `tx` соответствует признаку **Recency** из упоминавшейся выше модели **RFM**, а `x` - **Frequency**. В стороне остается только **Monetary Value**, сумма сделок `order_quantity` по каждому клиенту за весь период наблюдений `T`, т.е. за **12** кварталов. Далее мы рассчитаем недостающий показатель и полноценно применим эти сведения для прогнозирования.

Бельгийский исследователь Тобиас Вербеке предусмотрительно подготовил обозначенный авторами набор данных в пакете [`SMCRM`](http://cran.rstudio.com/web/packages/SMCRM) - `customerRetentionLifetimeDuration`, содержащий сведения о продолжительности сотрудничества клиентов с фирмой. Однако мы можем его получить самостоятельно прямо исполнив команды **SQL** из приведенного кода **SAS** благодаря пакету [`sqldf`](http://cran.rstudio.com/web/packages/sqldf) в составе **R**.

```{r Ch04 : customerRetention - LifetimeDuration Data }
# P(Alive)
palive_x <- sqldf("select customer, sum(purchase) as x
                  from customerRetentionTransactions
                  group by customer order by customer;")

palive_tx <- sqldf("select customer, max(quarter) as tx
                  from customerRetentionTransactions
                  where purchase = 1 group by customer order by customer;")

palive_T <- sqldf("select customer, max(quarter) as T, sum(order_quantity) as M
                  from customerRetentionTransactions group by customer order by customer;")

palive_xtx <- sqldf("select a.*, b.tx
                  from palive_x a left join palive_tx b
                  on a.customer = b.customer;")

customerRetentionLTDuration <- sqldf("select distinct a.*, b.T, b.M
                  from palive_xtx a left join palive_T b
                  on a.customer = b.customer order by customer;")

remove(palive_x, palive_tx, palive_T, palive_xtx)

# Description of an empirical distribution for non-censored data using Cullen & Frey graph
library('fitdistrplus')
fitdistrplus::descdist(customerRetentionLTDuration$x, boot = 500, discrete = FALSE)
x <- fitdist((customerRetentionLTDuration$x / max(customerRetentionLTDuration$x)),
             "beta", method = "mme", discrete = FALSE)
plot(x)
writeLines("\n Beta distribution of `customerRetentionLTDuration$x`")
gofstat(x)

fitdistrplus::descdist(customerRetentionLTDuration$tx, boot = 500, discrete = FALSE)
(x <- fitdist((customerRetentionLTDuration$tx / max(customerRetentionLTDuration$tx)),
             "beta", method = "mme", discrete = FALSE))
plot(x)
writeLines("\n Beta distribution of `customerRetentionLTDuration$x`")
gofstat(x)

```

Независимые переменные `x`, `tx`, `M` из нового набора данных `customerRetentionLTDuration` распределены по непрерывному бета-распределению

Теперь у нас в распоряжении полный набор данных для построения модели *продолжительности сотрудничества с клиентом* и при помощи пакета [`Buy ’Til You Die - BTYD`](http://cran.rstudio.com/web/packages/BTYD). Должен заметить, что в пакете **SAS** нельзя напрямую рассчитать эти параметры.

Вероятность активности на момент времени T клиента без контрактных установок рассчитывается по формуле:

$$ P(Alive \enspace | \enspace r, \alpha, a, b, x, t_x, T) = 1 / \left\{  1 + \dfrac{a}{b + x} \left( \dfrac{\alpha + T }{\alpha + t_x }\right) ^{r + x} \right\} $$

```{r Ch04 : BG/NBD}
# Customers without contractual settings : Buy ’Til You Die - Beta Gamma / Negative Binomial Distribution
# http://srepho.github.io/CLV/CLV
library('BTYD') # Implementing Buy 'Til You Die Models

cal.cbs <- customerRetentionLTDuration %>% 
            dplyr::select(x, tx, `T`) %>% # calibration period CBS (customer by sufficient statistic)
              dplyr::rename(x = x, t.x = tx, T.cal = `T`)
params <- BTYD::bgnbd.EstimateParameters(cal.cbs)

LL <- BTYD::bgnbd.cbs.LL(params, cal.cbs)
params <- c(params, LL)
names(params) <- c("r", "alpha", "a", "b", "LL (Log-likelihood)")

# Parameters from the book 
params0 <- c(126.5368069,	159.8644711,	0.512114268,	3.328751451, -4676.1)
names(params0) <- c("r", "alpha", "a", "b", "LL (Log-likelihood)")

writeLines("\n BG/NBD Model: Parametres of AR Prediction \n")
params
writeLines("\n BG/NBD Model: Parametres of V. Kumar, J. Andrew Petersen Prediction \n")
params0


P_Alive <- customerRetentionLTDuration %>% 
             mutate(p = 1/(1+(params["a"]/(params["b"] + x))*
                             ((params["alpha"] + `T`)/(params["alpha"] + tx))^(params["r"] + x)))

utils::data(customerRetentionLifetimeDuration, package = 'SMCRM')
P_Alive0 <- customerRetentionLifetimeDuration %>% 
             mutate(p = 1/(1+(params0["a"]/(params0["b"] + x))*
                             ((params0["alpha"] + `T`)/(params0["alpha"] + tx))^(params0["r"] + x)))

writeLines("\n BG/NBD Model: Association of AR Prediction vs. V. Kumar, J. Andrew Petersen Prediction \n")
caret::confusionMatrix(data = ifelse(P_Alive$p > 0.5, "1", "0") %>% factor,
                       reference = ifelse(P_Alive0$p > 0.5, "1", "0") %>% factor,
                       dnn = c("AR Prediction", "V. Kumar, J. Andrew Petersen Prediction"),
                       positive = "1", mode = "everything")

customerRetentionLTDuration$p <- P_Alive0$p
remove(params, P_Alive, params0, P_Alive0)

```

Хотя в моих расчетах, которые производились не в *MS Excel*, а в специализированном пакете [BTYD](https://cran.r-project.org/web/packages/BTYD/) для **R**, получены неcколько иные параметры бета- и гамма- распределений. Возможно из-за разных алгоритмов оптимизации, во всяком случае логарифмическое правдоподобие (англ. "`Log-likelihood [LL]`") у них практически совпадают. Схоже и попадание в классы по авторскому разбиению. Разница лишь в том, что по авторской модели BG/NBD клиенты которые совершили лишь две покупки (в начале окна наблюдения - в 1 квартале и почти за год до его завершения - в 9 квартале) по-прежнему считаются активными, а по более строгому алгоритму оптимизации нужно совершить между этими двумя крайними трансакциями хотя бы еще одну покупку. Это представляется более разумным при построении модели *продолжительности сотрудничества с клиентом*, чем опора только на две точки взаимодействия с фирмой.

### Как это использовать?

Такой подход к моделированию *продолжительности сотрудничества с клиентом* в условия отсутстия контрактных установок представляется крайне полезным для менеджемента, давая ему возможность оптимально управлять расходами на удержание клиентов.

## Эмпирический пример: величина заказа

Многие фирмы поняли, что недостаточно просто сосредоточиться только на попытке заставить клиента выкупить. Фирма также должна обратить внимание на то, каким размером может оказать покупка. Исследования в области маркетинга показали, что стоимость заказа может быть ценным предиктором в будущей стоимости клиента для фирмы - или, по крайней мере, оправдать сумму денег, затрачиваемую на усилия по удержанию клиентов. Таким образом, может быть полезно понять драйверы величины заказов и, в свою очередь, иметь возможность прогнозировать ожидаемую величину заказов каждого потенциального клиента при условии, что заказ будет иметь место. В конце этого примера мы должны иметь возможность сделать следующее:

1. Определить драйверы величину заказа (по стоимости).

2. Предсказать ожидаемую величину заказа для каждого клиента.

3. Измерить предиктивную точность модели.

Информация, необходимая для этой модели, включает следующий список переменных:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| **Зависимая переменная** |  |
| Purchase        | если **1**, когда покупатель приобрел в данном квартале, **0**, если в этом квартале не произошло покупок |
| Order_quantity  | Долларовая стоимость покупок в данном квартале |
| **Предикторы**          |  |
| imr_purchase или *Lambda* ($\lambda$) | Рассчитанное *обратное отношение Миллса* из модели повторного заказа клиента |
| Lag_Purchase    | **1**, если клиент приобрел в предыдущем квартале, **0**, если в предыдущем квартале покупка не произошла |
| Avg_Order_Quantity | Средняя долларовая стоимость покупок во всех предыдущих кварталах |
| Ret_Expense     | Доллары потраченные на маркетинговые усилия по удержанию этого клиента в данном квартале |
| Ret_Expense_SQ  | Квадрат затрат на маркетинговые усилияпо удержанию этого клиента в данном квартале |
| Gender          | **1**, если клиент является мужчиной, **0**, если клиент является женщиной |
| Married         | **1**, если клиент женат, **0**, если клиент не состоял в браке |
| Income          | **1**, если доход < 30 000 долл. США, **2**; если 30 001 долл. США < доход < 45 000 долл. США; **3**, если 45 001 долл. США < доход < 60 000 долл. США; **4**, если 60 001 долл. США < доход < 75 000 долл. США; **5**, если 75 001 долл. США < доход < 90 000 долл. США; **6**, если доход > 90 001 долл. США |
| First_Purchase  | Стоимость первой покупки, сделанной клиентом в 1 квартале |
| Loyalty         | **1**, если клиент является членом программы лояльности, **0**, если нет |
 
Из вышеописанных перменных мы видим, что для определения драйверов величины заказов нам нужно иметь две зависимые переменные: `Purchase` and `Order_Quantity`. Это связано с тем, что ожидаемая величина заказа получается из следующего уравнения:

$$ E(Order \enspace Quantity) = P(Purchase = 1) * E(Order \text _ quantity | Purchase = 1) $$

Смещение выборки является проблемой, которая распространена во многих маркетинговых проблемах и должна быть статистически учтена во многих процессах моделирования. В этом случае у клиента есть выбор: покупать или не покупать, прежде чем принимать решение о покупке. Если бы мы проигнорировали этот выбор, мы бы смещали оценки от модели, и у нас были бы менее точные предсказания для значения `Order_Quantity`. Чтобы учесть эту проблему, мы должны иметь возможность предсказать величину как вероятности покупки (аналогично тому, что мы сделали для первого эмпирического примера в этой главе), так и ожидаемого значения `Order_Quantity`, учитывая, что клиент должен сделать покупку. Важно отметить, что мы не можем просто запускать две модели независимо, так как, вероятно, существует корреляция между условиями ошибки двух моделей. Таким образом, авторы указывают на необходимость использовать структуру моделирования, которая может одновременно оценивать коэффициенты двух моделей или, по крайней мере, учитывать соотношение между `Order_Quantity` и `Purchase`. Для этого мы используем двухэтапную структуру моделирования, аналогичную описанной ранее в [третьей главе](https://rpubs.com/A_Rodionoff/SMCRM_Ch03). Первым этапом идет формирование *пробит-модели*, дающей нам *обратное отношение Миллса* ($\lambda$).

### Первый этап построения и верификации *величины заказа* {.tabset}

#### Probit

```{r Ch04 : Repurchase Probit}
# Fit Probit Model for Customer Repurchase by Authors
Ch04.probit <- glm(purchase ~ lpurchase + avg_order_quantity + ret_expense + ret_expense_sq +
                              gender + married + income + first_purchase + loyalty, 
             data = customerRetentionRepurchase, family = binomial(link = 'probit'))
summary(Ch04.probit)
writeLines(sprintf("-2 Log L of Intercept and Only Covariates: %.3f", -2 * logLik(Ch04.probit)[1]))
writeLines(sprintf("                  AIC (smaller is better): %.3f", extractAIC(Ch04.probit)[2]))

writeLines("\n Wald test of predictors")
car::Anova(Ch04.probit, type="II", test="Wald")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04.probit) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity 

writeLines("\n Probit Model: Association of Predicted Probabilities and Observed Responses \n")
prob <- predict(Ch04.probit, newdata = customerRetentionRepurchase, type = "response") 
caret::confusionMatrix(data = ifelse(prob > 0.5, "1", "0") %>% factor,
                       reference = customerRetentionRepurchase$purchase %>% factor,
                       positive = "1", mode = "everything")
```

Качество решения классификационной задачи *повторной покупки* довольно высокое.

#### Probit (AR Version)

Вместе с тем очевидно, что предикторы `married` и `first_purchase` не значимы в пробит-регрессии, что отмечают даже авторы. Недалеко от них ушел по значимости фактор `gender`. Также вместо полученной со смещением переменной `avg_order_quantity` я буду использовать корректно рассчитанную `lavg_order_quantity`. Кроме того, проверка по **VIF** демонстрирует большие значения по признакам `acq_expense` и `acq_expense_sq`. Однако, принимая во внимание проявление *закона убывающей доходности* мы должны оставить оба *мультиколлинеарных*  предиктора в этой модели. Теперь бинарную модель мы проверяем более тщательно, оставля только значимые предикторы.

```{r Ch04 : Repurshase Probit Improved, warning=FALSE}
# Fit Probit Regression Model for Customer Repurchase (Improved)
uno = 'probit'
set.seed(2018)
(Ch04pr.AR <- train(factor(purchase) ~ lpurchase + lavg_order_quantity + ret_expense + ret_expense_sq +
                                       income + loyalty, metric = 'Kappa',
                          data = customerRetentionRepurchase, method = "glm", family = binomial(link = uno)))#,
                  # trControl = trainControl(method = "none", number = 1)))
summary(Ch04pr.AR)
# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")
car::Confint(Ch04pr.AR$finalModel) %>%
  exp() %>%
    arm::pfround(digits = 3)

# Logistic regression diagnostics
writeLines("\n Wald test of predictors")
car::Anova(Ch04pr.AR$finalModel, type="II", test="Wald")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04pr.AR$finalModel) # Variance Inflation Factors - if vif() > 2  - feature has multicollinearity
writeLines("\n Variable Importance for Model \n")
caret::varImp(Ch04pr.AR) %>% .$importance %>% print.AsIs()

# Create the scatter plots Probit versus model predictors
predictors <- Ch04pr.AR$coefnames
Ch04pr.AR$trainingData %>%
  dplyr::select(one_of(predictors)) %>%
    mutate(link = predict(Ch04pr.AR$finalModel, newdata = customerRetentionRepurchase, type = "link")) %>%
      gather(key = "predictors", value = "predictor.value", -link) %>%
        ggplot(aes(predictor.value, link))+
          geom_point(size = 0.05, alpha = 0.05) +
          geom_smooth(method = "loess") +
          facet_wrap(~ predictors, scales = "free_x") +
          ylab(stringr::str_to_title(uno))

# Plot matrix of statistical model diagnostics
GGally::ggnostic(Ch04pr.AR$finalModel, title = paste(paste(formula(Ch04pr.AR)[c(2, 1, 3)], collapse = " ")))

# wide variety of diagnostic plots for checking the quality of regression fit
# https://bookdown.org/jefftemplewebb/IS-6489/logistic-regression.html
car::influenceIndexPlot(Ch04pr.AR$finalModel)

writeLines("\n Improved Probit Model: Association of Predicted Probabilities and Observed Responses \n")
caret::confusionMatrix(data = predict(Ch04pr.AR, newdata = customerRetentionRepurchase),
                       reference = customerRetentionRepurchase$purchase %>% factor,
                       positive = "1", mode = "everything")

qplot(`Observed Classes`, `Predicted Classes`, 
      data=bind_cols(`Observed Classes`= factor(customerRetentionRepurchase$purchase),
                     `Predicted Classes` = predict(Ch04pr.AR, newdata = customerRetentionRepurchase)),  
      colour= `Observed Classes`, geom = c("boxplot", "jitter"),
      main = "Predicted Classes vs. Observed Classes", xlab = "Observed Classes", ylab = "Predicted Classes")
```

Улучшенная пробит-регрессия *вероятности повторной покупки* хотя и снизила  немного качества, получена без незначимых факторов `gender` и `married`, а также предиктора `first_purchase`, поэтому мне представляется более надежной, так как избавилась от этих малозначимых признаков и более устойчивой, так как модель построена бутстреп-методом - англ. `bootstrap` с 25-тью "псевдовыборок".

### Второй этап построения и верификации модели *величины заказа* {.tabset}

Переходим ко второму этапу модели *величины заказа*, которая использует в качестве предиктора результаты пробит-модели из первого этапа, обозначенное как *обратное отношение Миллса* ($\lambda$).

#### Linear

```{r Ch04 : Order quantity - Linear Regression}
# Fit Linear Regression Model for Order quantity by Authors

# SAS Code: imr_acquisition = (pdf(’Normal’, xb_probit))/(probnorm(xb_probit));
xbeta <- predict(Ch04.probit, newdata = customerRetentionRepurchase, type = "link")
customerRetentionRepurchase <- customerRetentionRepurchase %>%
  mutate(imr_purchase = dnorm(xbeta) / pnorm(xbeta)) #  Cumulative normal pdf

(Ch04.linear <- lm(order_quantity ~ lpurchase + avg_order_quantity + ret_expense +
ret_expense_sq + gender + married + income + first_purchase + loyalty + imr_purchase,
data = filter(customerRetentionRepurchase, order_quantity > 0) )) %>%
  summary

writeLines("Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04.linear) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity

```

Авторы получили результаты, говорящие что регрессионная задача решена удовлетворительно. Мы видим, что $\lambda$ является положительной и значимой. Можно её интерпретировать так, что существует потенциальная проблема смещения выбора, поскольку коэффициент ошибки нашего уравнения выбора положительно коррелирует с погрешностью нашего уравнения пробит-регрессии. Также видно, что все другие переменные модели *величины заказов* являются значимыми, за исключением `Married`, что означает, вероятно, обнаружение многие из драйверов величины заказа.

#### Linear (AR version)

Теперь попробуем исключить из линейной модели незначимый фактор `married`, а  вместо полученной со смещением переменной `avg_order_quantity` использовать корректно рассчитанную `lavg_order_quantity`.
 
```{r Ch04 : Order quantity - Linear Regression Improved, warning=FALSE}
# Fit Linear Regression Model for Customer Order quantity  (Improved)
uno = 'y_hat'

# Log Transformation of Order_Quantity
customerRetentionRepurchase <- customerRetentionRepurchase %>% 
    mutate(log_order_quantity = log(order_quantity))

set.seed(2018)
(Ch04ln.AR <- train(log_order_quantity ~ lpurchase + lavg_order_quantity + ret_expense + ret_expense_sq + 
                    gender + income + first_purchase + imr_purchase, 
                    data = filter(customerRetentionRepurchase, order_quantity > 0), method = "lm"))#,
                    # trControl = trainControl(method = "none", number = 1)))
summary(Ch04ln.AR)
#Coefficient Estimates and 95% CI
writeLines("\n Coefficient Estimates and 95% CI")
car::Confint(Ch04ln.AR$finalModel) %>%
  exp() %>%
    arm::pfround(digits = 3)

# Linear regression diagnostics
writeLines("\n Chisq test of predictors")
car::Anova(Ch04ln.AR$finalModel, type="II", test="Chisq")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04ln.AR$finalModel) # Variance Inflation Factors - if vif() > 2  - feature has multicollinearity
writeLines("\n Variable Importance for Model \n")
caret::varImp(Ch04ln.AR) %>% .$importance %>% print.AsIs()

# Create the scatter plots Linear Model versus model predictors
# https://stats.idre.ucla.edu/r/seminars/ggplot2_intro/
predictors <- Ch04ln.AR$trainingData %>% colnames(.) %>% .[-1]
Ch04ln.AR$trainingData %>% 
  rename(y_hat = `.outcome`) %>% 
    gather(key = "predictors", value = "predictor.value", -y_hat) %>%   
      ggplot(aes(predictor.value, y_hat))+
        geom_point(size = 0.5, alpha = 0.5) +
        geom_smooth(method = "loess") +
        facet_wrap(~ predictors, scales = "free_x") +
        ylab(stringr::str_to_title(uno))

# Plot matrix of statistical model diagnostics
GGally::ggnostic(Ch04ln.AR$finalModel, title = paste(paste(formula(Ch04ln.AR)[c(2, 1, 3)], collapse = " ")))

# wide variety of diagnostic plots for checking the quality of regression fit
# https://bookdown.org/jefftemplewebb/IS-6489/linear-regression.html
car::influenceIndexPlot(Ch04ln.AR$finalModel)

``` 

Исключение из модели факторов `married` и даже `loyalty` не нанесло какого-либо вреда качеству модели. Однако избавление от смещения по независимой переменной `avg_order_quantity` понизило качество линейной модели. Чтобы его поднять пришлось прибегнуть к логарифмированию зависимой переменной `order_quantity` и независимой `lavg_order_quantity`, что часто применяется в эконометрических моделях на стоимостных показателях, имеющих склонность к лог-нормальному распределению.


### Как это использовать?

Авторы нашли, что `Lag_Purchase` положительна, предполагая, что покупатели, купившие в предыдущем квартале, с большей вероятностью потратят также в текущем квартале. Также обнаружено, что `Avg_Order_Quantity` как и `Lag_Avg_Order_Quantity` также является положительным, предполагая, что чем выше средние значения прошлые величины заказов клиента, тем выше текущее значение значение. Очевидно, что `Ret_Expense` имеет положителеный коэфициент с уменьшающимся возвратом, как отмечено положительным коэффициентом на `Ret_Expense` и отрицательным коэффициент для `Ret_Expense_SQ`. Это означает, что маркетинговые усилия, направленные на сохранение и построение отношений с клиентом, заставляют клиента приобретать больше, лишь до определенной степени. Затем, после достижения некоторого порога, маркетинговые усилия фактически снижают стоимость покупки в среднем. Вероятно, это связано с тем, что чрезмерное общение с клиентами может часто напрягать клиентов и подвергать эрозии отношения между ними и фирмой. Авторы обнаружили, что несколько характеристик клиента являются положительными (`Income`, `First_Purchase`), давая менеджерам вполне разумную мысль о том, что клиенты, которые имеют более высокий доход, как впрочем имеющие более высокую стоимость первой покупки, как правило, имеют большую величину заказа.

Следующий шаг - предсказать значение `Order_Quantity`, чтобы увидеть, насколько авторская модель совпадает с фактическими значениями.

```{r Ch04 : Error of Linear Model}
# Computing the Mean Absolute Deviation (MAD) and Mean Absolute Percent Error (MAPE)

with(customerRetentionRepurchase, {
  # pred_oq <- predict(Ch04pr.AR$finalModel, newdata = customerRetentionRepurchase, type = "link") %>% pnorm() *
  #              exp(predict(Ch04ln.AR$finalModel, newdata = customerRetentionRepurchase))
  pred_oq <- predict(Ch04.probit, newdata = customerRetentionRepurchase, type = "link") %>% pnorm() *
               predict(Ch04.linear, newdata = customerRetentionRepurchase)
  # mean_order_quantity
  writeLines(sprintf("Mean of Order_Quantity: %.2f долл.", mean(customerRetentionRepurchase$order_quantity)))
  
  # mad = mean(abs(first_purchase - pred_oq));
  writeLines(sprintf("Mean Absolute Deviation (MAD): %.2f долл.", mean(abs(order_quantity - pred_oq))))
  
  # mad1 = mean(abs(order_quantity - mean(order_quantity));
  mad1 <- mean(abs(order_quantity - mean(order_quantity)))
  writeLines(sprintf("Naive Mean Absolute Deviation (MAD1): %.2f долл.", mad1))
  })
```

Поскольку значительное количество клиентов не заказывало ещеквартально повторные заказы, то метрику MAPE рассчитать авторам не удалось. Среднее значение по всей выборочной совокупности заказов (без первого квартала) составило `r sprintf("%.2f", mean(customerRetentionRepurchase$order_quantity))` долл. Этот показатель авторы применили в качестве *наивного* прогноза. Созданная ими модель *величины заказов* оказалось заметно аккуратнее *наивной* модели среднего, у которого MAD1 = `r sprintf("%.2f", mean(abs(customerRetentionRepurchase$order_quantity - mean(customerRetentionRepurchase$order_quantity))))` долл.

## Эмпирический пример: перекрестные продажи

Перекрестные продажи (англ. "`Cross-buying`") - это наиболее распространенная технология, используемая компаниями для увеличения числа случаев повторного приобретения, величины заказа и доходов компаний. После того, как компании создали определенный уровень лояльности с существующими клиентами, эти клиенты с большей вероятностью будут перекрестно покупать у компании. 

Еще один ключевой вопрос, на который авторы книги хотели ответить в отношении удержания клиентов, заключается в том, можно ли определить, какие клиенты имеют наивысшую вероятность перекрестых продаж в нескольких категориях. Для этого сначала нужно знать, какие текущие клиенты фактически приобретали в нескольких категориях, когда они совершили покупку.

В наборе данных, представленном в этой главе, имеется переменная `Crossbuy`, которая идентифицирует, сколько категорий продуктов покупает покупатель за определенный период времени. Также предоставляется набор драйверов, которые, вероятно, помогут объяснить решение клиента о перекрестной покупке. В конце этого примера читатели смогут сделать следующее:

1. Определить драйверы поведения при перекрестных покупок клиента.

2. Интерпретировать оценки параметров из модели перекрестной покупки.

3. Предсказать может ли клиент перекрестно купить или нет.

4. Определить предиктивную точность модели перекрестной покупки.

Компания B2C хочет понять, какие клиенты, скорее всего, перекрестно покупают за данный период времени. Это важно знать, поскольку многие исследования показали, что клиенты, покупающие по нескольким категориям, с большей вероятностью будут более прибыльными, чем покупатели, которые покупают меньшее количество категорий. Случайная выборка из 500 клиентов из одной когорты была взята из базы данных клиентов. Информация, необходимая для нашей модели, включает следующий список переменных:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| **Зависимая переменная** |  |
| Crossbuy        | Количество различных категорий товаро / услуг, приобретенных в данном квартале |
| **Предикторы**          |  |
| Lag_Purchase    | **1**, если клиент приобрел в предыдущем квартале, **0**, если в предыдущем квартале покупка не произошла |
| Lag_Crossbuy        | Количество различных категорий товаров / услуг, приобретенных в предыдущем квартале |
| Ret_Expense     | Доллары потраченные на маркетинговые усилия по удержанию этого клиента в данном квартале |
| Ret_Expense_SQ  | Квадрат затрат на маркетинговые усилияпо удержанию этого клиента в данном квартале |
| Gender          | **1**, если клиент является мужчиной, **0**, если клиент является женщиной |
| Married         | **1**, если клиент женат, **0**, если клиент не состоял в браке |
| Income          | **1**, если доход < 30 000 долл. США, **2**; если 30 001 долл. США < доход < 45 000 долл. США; **3**, если 45 001 долл. США < доход < 60 000 долл. США; **4**, если 60 001 долл. США < доход < 75 000 долл. США; **5**, если 75 001 долл. США < доход < 90 000 долл. США; **6**, если доход > 90 001 долл. США |
| First_Purchase  | Стоимость первой покупки, сделанной клиентом в 1 квартале |
| Loyalty         | **1**, если клиент является членом программы лояльности, **0**, если нет |

В этом случае автору упоминают о дискретной зависимлй переменной (`Crossbuy`), которая сообщает сколько категорий покупает покупатель в данном квартале. Чтобы понять вероятность перекрестного покупки, следует преобразовать эту дискретную переменную в двоичную переменную. Авторы делали это, устанавливая `CB` равным 1, когда `Crossbuy > 1` и `CB` равны 0 в противном случае. Также имеется девять независимых переменных, которые, по мнению авторов, станут драйверами поведения повторной покупки.

Авторы полагали, что поведение клиентов на транзакциях в прошлом, скорее всего, объяснит поведение покупателей с перекрестными покупками. В результате в этом примере авторы использовали несколько отложенных переменных, т.е. взятых с лагом, в качестве независимых переменных. Во-первых, имеется информация о том, приобрел ли покупатель в предыдущем квартале (`Lag_Purchase`). Также имеется переменная для количества категорий товаров, приобретенных покупателем в предыдущем квартале (`Lag_Crossbuy`). Эти две переменные могут быть получены путем взятия стоимости покупки / `cross-buy` с отставание в один месяц, отметив, что одно наблюдение будет потеряно для каждого клиента. В этом случае мы используем однопериодное (квартальное) отставание для обеих переменных. Во-вторых, имеется среднеквартальное количество прошлых заказов (`Avg_Order_Quantity`). В этом случае значение для среднего количества заказа является средним значением переменной `Order_Quantity` во всех кварталах до текущего периода времени. В-третьих, есть информация сколько долларов фирма потратила на каждого клиента (`Ret_Expense`) за каждый период времени и квадрат значения этой переменной (`Ret_Expense_SQ`). Авторы хотели использовать как линейные, так и квадратичные термины, так как ожидали, что для каждого дополнительного доллара, потраченного на усилия по удержанию для данного клиента, будет уменьшаться возврат к стоимости этого доллара. Наконец, поскольку фокусная фирма этого примера является фирмой **B2C**, остальные пять переменных являлись социально-демографическими характеристиками клиентов. К ним относятся `Gender` клиента, является ли клиент `Married`, `Income` клиента, стоимость первой покупки клиента (`First_Purchase`) и является ли клиент членом программы лояльности (`Loyalty`).

Во-первых, авторы смоделировали вероятность того, что клиент будет перекрестно покупать в данный период времени. Поскольку зависимая переменная (`CB`) является бинарной, авторы остановились на логистической регрессии для оценки модели. Авторы упоминали, что могли бы также выбрать пробит-модель и в целом добиться тех же результатов. В этом случае переменная y является `CB`, а предикторы представляют 10 независимых переменных в базе данных. Следует также выбирать только те случаи, когда покупка произошла, так как фирмы заинтересованы в том, чтобы клиенты перекрестно покупали в условиях покупки. Это дает нам `r sprintf("%.0f", sum(customerRetentionRepurchase$purchase))` наблюдений для построения этой модели. 

### Построение и верификация модели *перекрестных продаж* {.tabset}

#### Logit

```{r Ch04 : Crossbuy - Logit, warning=FALSE}
# # Fit Logistic Regression Model for Customer Crossbuy by Authors

Ch04.cb <- glm(cb ~ lpurchase + lcrossbuy + avg_order_quantity + ret_expense + ret_expense_sq + 
                    gender + married + income + first_purchase + loyalty,
                  data = filter(customerRetentionRepurchase, purchase == 1), family = binomial(link = 'logit'))
summary(Ch04.cb)
writeLines(sprintf("-2 Log L of Intercept and Only Covariates: %.3f", -2 * logLik(Ch04.cb)[1]))
writeLines(sprintf("                  AIC (smaller is better): %.3f", extractAIC(Ch04.cb)[2]))

# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")
car::Confint(Ch04.cb) %>%
  exp() %>%
    arm::pfround(digits = 3)

writeLines("\n Wald test of predictors")
car::Anova(Ch04.cb, type="II", test="Wald")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04.cb) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity

writeLines("\n Crossbuy Probability (Logit): Association of Predicted Probabilities and Observed Responses \n")
prob <- predict(Ch04.cb, newdata = filter(customerRetentionRepurchase, purchase == 1), type = "response")
caret::confusionMatrix(data = ifelse(prob > 0.5, "1", "0") %>% factor,
                       reference = filter(customerRetentionRepurchase, purchase == 1) %>% .$cb %>% factor,
                       positive = "1", mode = "everything")
```

Авторы получили логистическую регрессию довольно хорошего качества для описательных целей. Несмотря на то, что признак `married` является малозначимым предиктором в этой модели.

Во-первых, это означает, что `Lag_Crossbuy` положительно влияет на текущую кросс-покупку, то есть клиенты, которые приобрели больше категорий в предыдущем квартале, с большей вероятностью перекрестно покупают в текущем квартале. Во-вторых, поскольку коэффициент по `Avg_Order_Quantity` является положительным и статистически значимым, то клиенты, которые в прошлом тратили больше средств, чем в среднем, также чаще перекрестно покупают в текущем периоде времени. В-третьих, авторы обнаружили положительную, но уменьшающуюся отдачу от эффекта удерживающих расходов (Ret_Expense) при перекрестной покупке в том же квартале, поскольку коэффициент на Ret_Expense положителен, а коэффициент Ret_Expense_SQ отрицателен. В-четвертых, авторы отметили, что мужчины чаще совершают перекрестные покупку, чем женщины. В-пятых, авторы нашли положительный эффект от клиентского дохода, т.е. клиенты, имеющие более высокий доход, с большей вероятностью будут перекрестно покупать в текущем квартале. В-шестых, авторы определили, что клиенты, у которых была более высокая `First_Purchase`, с большей вероятностью закупят несколько категорий продуктов / услуг. Наконец, поскольку коэффициент для `Loyality` положителен, это говорит о том, что клиенты, которые являются членами программы лояльности, чаще всего перекрестно покупают в текущем квартале.

#### Logit (AR version)

Авторы по-прежнему прибегают в `Avg_Order_Quantity`, который в их примере охватывает текущий квартал и тем самым дает смещение в оценке модели в сторону увеличения качества модели. Я рассмотрю эту модель с более корректным предиктором с однопериодным лагом - `Lag_Avg_Order_Quantity`.

```{r Ch04 : Crossbuy Logit Improved, warning=FALSE}
# Fit Logit Regression Model for Customer Crossbuy (Improved)
uno = 'logit'
set.seed(2018) #From random.org
(Ch04cb.AR <- train(factor(cb) ~ lpurchase + lcrossbuy + lavg_order_quantity + ret_expense + ret_expense_sq + 
                    gender + married + income + first_purchase + loyalty, metric = 'Kappa',
      data = filter(customerRetentionRepurchase, purchase == 1), method = "glm", family = binomial(link = uno)))#,
                  # trControl = trainControl(method = "none", number = 1)))
summary(Ch04cb.AR)
# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")
car::Confint(Ch04cb.AR$finalModel) %>%
  exp() %>%
    arm::pfround(digits = 3)

# Logistic regression diagnostics
writeLines("\n Wald test of predictors")
car::Anova(Ch04cb.AR$finalModel, type="II", test="Wald")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04cb.AR$finalModel) # Variance Inflation Factors - if vif() > 2  - feature has multicollinearity
writeLines("\n Variable Importance for Model \n")
caret::varImp(Ch04cb.AR) %>% .$importance %>% print.AsIs()

# Create the scatter plots Logit versus model predictors
predictors <- Ch04cb.AR$coefnames
Ch04cb.AR$trainingData %>%
  dplyr::select(one_of(predictors)) %>%
    mutate(link = predict(Ch04cb.AR$finalModel, newdata = filter(customerRetentionRepurchase, purchase == 1), type = "link")) %>%
      gather(key = "predictors", value = "predictor.value", -link) %>%
        ggplot(aes(predictor.value, link))+
          geom_point(size = 0.05, alpha = 0.05) +
          geom_smooth(method = "loess") +
          facet_wrap(~ predictors, scales = "free_x") +
          ylab(stringr::str_to_title(uno))

# Plot matrix of statistical model diagnostics
GGally::ggnostic(Ch04cb.AR$finalModel, title = paste(paste(formula(Ch04cb.AR)[c(2, 1, 3)], collapse = " ")))

# wide variety of diagnostic plots for checking the quality of regression fit
# https://bookdown.org/jefftemplewebb/IS-6489/logistic-regression.html
car::influenceIndexPlot(Ch04cb.AR$finalModel)

writeLines("\n Improved Logit Model: Association of Predicted Probabilities and Observed Responses \n")
caret::confusionMatrix(data = predict(Ch04cb.AR, newdata = filter(customerRetentionRepurchase, purchase == 1)),
                       reference = filter(customerRetentionRepurchase, purchase == 1) %>% .$cb %>% factor,
                       positive = "1", mode = "everything")

qplot(`Observed Classes`, `Predicted Classes`, 
      data=bind_cols(`Observed Classes`= filter(customerRetentionRepurchase, purchase == 1) %>% .$cb %>% factor,
          `Predicted Classes` = predict(Ch04cb.AR, newdata = filter(customerRetentionRepurchase, purchase == 1))),
      colour= `Observed Classes`, geom = c("boxplot", "jitter"),
      main = "Predicted Classes vs. Observed Classes", xlab = "Observed Classes", ylab = "Predicted Classes")
```

Конечно, без смещения качество модели немного пострадает, зато она будет более правдоподобной и безусловно устойчивой. При этом даже значимость переменной `Married` повысилось.

### Как это использовать?

Полезно рассмотреть значения коэффициентов *отношения шансов* (англ. "`odds ratio`"). Что касается `Lag_Crossbuy`, мы видим, что каждая дополнительная категория, купленная покупателем в предыдущем квартале, делает клиента на 19,7% более вероятным для перекрестную покупку в текущем квартале. Что касается `Avg_Order_Quantity`, авторы отмечали, что при каждом увеличении на 1 доллар шанс перекрестной покупки в текущем квартале увеличивается на 1,0%. Что касается расходов на удержание клиентов `Ret_Expense`, то и *отношение шансов* зависит от уровня `Ret_Expense`. Это связано с тем, что мы включаем как собственно его уровень, так и его квадрат `Ret_Expense_SQ`. Например, если обычно в этой выборке тратили `r sprintf("%.0f", mean(customerRetentionRepurchase$ret_expense))` долл. ежеквартально на данного клиента, тратя больше долларов, мы должны увидеть увеличение вероятности перекрестной покупки.

И важно отметить, что это будет зависеть от начального уровня Ret_Expense. Что касается гендерных факторов, авторы указывают, что мужчины на 33,7% имеют шансов чаще покупают перекрестно, чем женщины. Что касается клиентского дохода, мы видим, что при каждом повышении уровня доходов на один уровень (`Income`) вероятность шанса кросс-покупки  увеличится на 28,2%. Что касается `First_Purchase`, авторы пишут, что при каждом увеличении на 1 долл. вероятность кросс-покупки увеличивается на 2,5%. Наконец, что касается лояльности (`Loyalty`), становится ясным, что, будучи членом программы лояльности, вероятность перекрестной покупки в данном квартале на 30,8% выше, чем у клиента, который не входит в программу лояльности.

Менеджерам будет полезно знать, как изменения в расходах на удержание, прошлые транзакции клиентов и характеристики клиентов могут либо увеличить, либо уменьшить вероятность шанса перекрестной покупки. Эта модель поможет предсказать, собирается ли клиент покупать в другой категории или нет. Эта информация может дать существенную информацию руководителям, которым поручено определить, какие клиенты, скорее всего, склонны к кросс-покупкам.

## Эмпирический пример: доля покупок в кошельке

Программы лояльности и прямые рассылки - это два способа, которыми компании используют для управления отношениями с клиентами. Цель состоит в том, чтобы установить тесные отношения с клиентами и таким образом улучшить восприятие отношений с клиентами. Помимо понимания моделей покупки клиента у данной фирмы, многие фирмы также хотят знать, как клиент распределяют покупки в данной категории среди всех фирм. Для этого надо знать доля покупок в кошельке клиента (англ. `"Share-of-Wallet, SOW"`). Исследования показали, что понимание `SOW` клиента в данной фирме может помочь понять вероятность того, что клиент собирается купить у данной фирмы и неизбежно долгосрочную ценность клиента для фирмы (`CVL`). Таким образом, может быть полезно понять драйверы `SOW` и, в свою очередь, иметь возможность прогнозировать ожидаемое `SOW` каждого клиента. В этом примере авторы дают возможность возможность:

1. Определите драйверы SOW.

2. Предскажите ожидаемое SOW для каждого клиента.

3. Определите предиктивную точность модели.

Информация, необходимая для этой модели, включает следующий список переменных:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| **Зависимая переменная** |  |
| Share-of-Wallet (SOW) | Процент покупок клиента от данной фирмы, учитывая общий объем покупок по всем фирмам в этой категории |
| **Предикторы**          |  |
| Purchase_Rate | Средняя доля кварталов с покупками во всех 12 кварталах |
| Avg_Order_Quantity | Средняя долларовая стоимость покупок за все 12 кварталов |
| Avg_Crossbuy | Среднее значение для кросс-покупки за все 12 кварталов |
| Avg_Ret_Expense | Средние расходы, потраченные на маркетинговые усилия по удержанию данного клиента за все 12 кварталов |
| Avg_Ret_Expense_SQ | Квадрат средних расходов, потраченные на маркетинговые усилия по удержанию данного клиента за все 12 кварталов |

Для получения этих агрегированных показателей необходимо получить массив данных, обобщенный за все 12 кварталов:

```{r Ch04 : Customer retention - Share-of-Wallet (SOW) data}
# Calculation Share-of-Wallet (SOW) Database

library('sqldf')
customerRetentionSOW <- sqldf("select customer, avg(purchase) as purchase_rate,
                              avg(order_quantity) as avg_order_quantity,
                              avg(crossbuy) as avg_crossbuy, avg(ret_expense) as avg_ret_exp,
                              (avg(ret_expense) * avg(ret_expense)) as avg_ret_exp_sq, gender,
                              married, income, first_purchase, sow, loyalty
                              from customerRetentionRepurchase group by customer, gender,
                              married, income, first_purchase, sow, loyalty order by customer;")

ggplot(customerRetentionSOW, aes(x = purchase_rate)) +
  geom_histogram(breaks = seq(0, 1, by = .1), col = "gray", aes(fill = ..count..)) +
    scale_fill_gradient("Count", low = "green", high = "red") +
      geom_density(adjust = 1/5, alpha = .2, fill = "#FF6666") +
         geom_vline(xintercept = mean(customerRetentionSOW$purchase_rate), linetype="dashed", color = "coral")

# Description of an empirical distribution for non-censored data using Cullen & Frey graph
library('fitdistrplus')
descdist(customerRetentionSOW$purchase_rate, boot = 500, discrete = FALSE)
x <- fitdist(customerRetentionSOW$purchase_rate, "unif", method = "mle", discrete = FALSE)
plot(x)
writeLines("\n Uniform distribution of `customerRetentionSOW$purchase_rate`")
gofstat(x)

```

Полученный непрерывный признак `Purchase_Rate` представлен на графике выше - очевидно, что он имеет равномерное распределение (англ. "`Uniform distribution`"). Это подтвержает тест при помощи графика Cullen & Frey.

В этом случае мы имеем цензуированную зависимую переменную (`SOW`), которая падает на континуум между 1 и 100. Минимальный размер этого случая составляет 1%, поскольку все клиенты в нашей базе данных совершили по крайней мере одну покупку у данной фирмы, а максимальный составляет 100%, поскольку все клиенты в базе данных могут потенциально покупать эти продукты только у этой фирмы. Таким образом, мы должны учитывать эту цензуированную зависимую переменную. В этом случае авторы предлагали  использовать вариацию тобит-модели, которую использовали в предыдущих примерах в [третьей главе](https://rpubs.com/A_Rodionoff/SMCRM_Ch03). В стандартном случае тобит-модель имеет ситуацию, когда нижняя граница зависимой переменной определяется, как правило, равной 0, а верхняя граница тобит-модель бесконечна. Однако в этом случае нам нужно учитывать цензуирование как с нижней границей, так и верхней границы, где нижняя граница равна 1, а верхняя граница равна 100. Таким образом, мы имеем следующее определение для SOW:

$$ \displaystyle \begin{array}{ll}
SOW_i = \begin{cases}
100  & {\text{if}}\ SOW ^* _i \geq 100 \\
SOW ^* _i & {\text{if}}\ 1 < SOW ^* _i < 100 \\ 
1 & {\text{if}}\ SOW ^* _i \leq 1 \\
\end{cases}
\end{array} $$

### Построение и верификация модели *доли покупок в кошельке* {.tabset}

#### Tobit

```{r Ch04: Share-of-Wallet (SOW) - Tobit}
# Censored Regression of Accelerated Failure Time (AFT) model on Time-to-Failure Data
# https://stats.idre.ucla.edu/sas/dae/tobit-analysis/ & https://stats.idre.ucla.edu/r/dae/tobit-models/  & http://rpubs.com/Joaquin_AR/381600


library('VGAM') # Vector Generalized Linear and Additive Models
# Left- & Right-Censored Tobit (Normally distributed error term) AFT Model from 'VGAM' package
(Ch04.tb <- VGAM::vglm(sow ~ purchase_rate + avg_order_quantity + avg_crossbuy + avg_ret_exp + avg_ret_exp_sq +
                         gender + married + income + first_purchase + loyalty,
                       family = tobit(Lower = 1, Upper = 100, type.fitted = "censored"),
                       data = customerRetentionSOW)) %>%
  summary

# # Censored Normal Distribution of error term Or Censored Tobit
# Extra <- with(customerRetentionSOW,
#               list(leftcensored = (sow < 1), rightcensored = (sow > 100)))
# 
# (Ch04.nr <- VGAM::vglm(sow ~ purchase_rate + avg_order_quantity + avg_crossbuy + avg_ret_exp + avg_ret_exp_sq +
#                        gender + married + income + first_purchase + loyalty,
#                        family = cens.normal(), extra = Extra,
#                        data = customerRetentionSOW)) %>%
#   summary

```

Задача авторов состояла в том, чтобы максимизировать функцию логарифмического правдоподобия путем оценки коэффициентов и стандартной ошибки уравнения. При оценке модели получаны вышеприведенные результаты. 

Авторы нашли, что все переменные, за исключением `Purchase_Rate`, `Gender` и `Married`, статистически значимы при `p < 0,05`. 

Следующий шаг - предсказать значение `SOW`, чтобы увидеть, насколько хорошо наша модель сравнивается с фактическими значениями. Авторы сделали это, сравненивая предсказания по двухстороннему цензуированному регрессионому уравнению с фактическими значениями `SOW`. 

#### Error of SOW

```{r Ch04: Error of SOW Models}
# mean_order_quantity
y <- customerRetentionSOW$sow
writeLines(sprintf("Mean of Share-of-Wallet (SOW): %.2f", mean(y)))

# Accuracy measures for Naïve AFT Model
m <- forecast::accuracy(rep(mean(customerRetentionSOW$sow), length(y)), y)
attributes(m)$ dimnames[[1]] <- "    Naïve Model (Mean) Accuracy: "; m

# Accuracy measures for Censored Tobit 
y_hat <- fitted(Ch04.tb)[, 1]
m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "      Censored Tobit's Accuracy: "; m

# # Accuracy measures for Censored Normal Distribution AFT Model
# y_hat <- if_else(fitted(Ch04.nr) > 100, 100, if_else(fitted(Ch04.nr) < 1, 1, fitted(Ch04.nr)))
# m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "Censored Normal Model's Accuracy: "; m

```

Полученная авторами модель дает для приобретенных клиентов MAD = `r sprintf("%.2f", forecast::accuracy(y_hat, y)[, "MAE"])` или в среднем `r sprintf("%.2f", forecast::accuracy(y_hat, y)[, "MAPE"])`% от фактического SOW. Если бы мы вместо этого использовали среднее значение `SOW` (`r sprintf("%.2f", mean(customerRetentionSOW$sow))`) для всех клиентов в качестве нашего прогноза для всех клиентов (это было бы *наивным* примером модели), мы бы обнаружили, что MAD `r sprintf("%.2f", forecast::accuracy(rep(mean(customerRetentionSOW$sow), length(y)), y)[, "MAE"])` или в среднем на `r sprintf("%.2f", forecast::accuracy(rep(mean(customerRetentionSOW$sow), length(y)), y)[, "MAE"])`% от фактического `SOW`. Следовательно, авторская модель значительно улучшает работу по прогнозированию значения `SOW`, чем прогнозирование по среднему значению `SOW`.

#### Tobit (AR version)

```{r Ch04: SOW - Tobit Improved}
# Censored Regression of Accelerated Failure Time (AFT) model on Time-to-Failure Data
# Left- & Right-Censored Tobit (Normally distributed error term) AFT Model from 'VGAM' package
(Ch04tb.AR <- VGAM::vglm(sow ~ avg_order_quantity + avg_crossbuy + avg_ret_exp_sq +
                        income + first_purchase + loyalty,
                        family = tobit(Lower = 1, Upper = 100, type.fitted = "censored"),
                        data = customerRetentionSOW)) %>%
  summary

# Calculate the upper and lower 95% confidence intervals for the coefficients.
b <- coef(Ch04tb.AR)
se <- sqrt(diag(vcov(Ch04tb.AR)))
writeLines("\n The Upper and Lower 95% confidence intervals for the coefficients of Left- & Right-Censored Tobit")
cbind(`Lower Intervals` = b - qnorm(0.975) * se, `Upper Intervals` = b + qnorm(0.975) * se)

dat <- customerRetentionSOW
dat$yhat <- fitted(Ch04tb.AR)[, 1]
dat$rr <- resid(Ch04tb.AR, type = "response")
dat$rp <- resid(Ch04tb.AR, type = "pearson")[, "mu"]

# Goodness-of-fit of normal distributions to residuals of Tobit-Model (if Kolmogorov-Smirnov statistic < 0.05)
resid(Ch04tb.AR, type = "response")[, 1] %>% fitdist("norm") %>% gofstat()

par(mfcol = c(2, 3))

with(dat, {
  plot(yhat, rr, main = "Fitted vs Residuals")
  qqnorm(rr, main = "Normal Q-Q Plot of Residuals")
  plot(yhat, rp, main = "Fitted vs Pearson Residuals")
  qqnorm(rp, main = "Normal Q-Q Plot of Pearson Residuals (mu)")
  plot(sow, rp, main = "Actual vs Pearson Residuals")
  plot(sow, yhat, main = "Actual vs Fitted")
})

par(mfcol = c(1, 1))

# Accuracy measures for Censored Tobit Improved
y <- customerRetentionSOW$sow
y_hat <- fitted(Ch04tb.AR)[, 1]
m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "Censored Tobit Improved Accuracy: "; m

remove(y, y_hat, m, b, se)

```

Если исключить три незначащие признаки `Purchase_Rate`, `Gender` и `Married` из числа предикторов, то получилась более устойчивую модель c незначительно увеличившейся ошибкой по MAE и MAPE. Первый свободный член `(Intercept):1` равный 10 надо интерепретировать как обычный сводобный член линейного уравнения, а второй сводобный член `(Intercept):2` равный 2.4 следует проэкспоненцировать и применить полученное число `11.0 = exp(2.4)` как стандартное отклонение, задаваемое при формировании нормально распределенных остатков sd в этом линейном уравнении. При этом следует признать, что что остатки улучшинной тобит-модели и в самом деле распределяются по нормальному распределению, поскольку значение статистики Колмогорова-Смирнова близко **p < 0.05**

### Как это использовать?

Авторы обнаружили, что коэффициент `Avg_Order_Quantity` положителен, что указывает на то, что чем выше средние значения заказа клиента в прошлом, тем выше значение `SOW`. Авторы также нашли, что `Avg_Crossbuy` позитивен, предполагая, что чем больше клиент закупил в нескольких категориях в прошлом, тем выше клиентский `SOW`. Мы находим, что `Ret_Expense` положителен с уменьшающимся возвратом, как отмечено положительным коэффициентом на `Ret_Expense` и отрицательным коэффициентом на `Ret_Expense_SQ`. Это означает, что маркетинговые усилия, направленные на сохранение и построение отношений с клиентом, заставляют клиента иметь более высокое значение `SOW`. Затем, после достижения порога, маркетинговые усилия фактически уменьшают `SOW` в среднем. Вероятно, это связано с тем, что чрезмерное общение с клиентами может часто напрягать отношения между клиентом и фирмой. Оказывается, что три из характеристики клиента являются положительными (`Income`, `First_Purchase` и `Loyalty`), предполагая, что клиенты с более высоким доходом, более высокая стоимость первой покупки и кто являются членами программы лояльности, вероятно, будут иметь более высокое значение`SOW`.

## Эмпирический пример: Рентабельность или Пожизненная финансовая ценность клиента (CLV)

Конечный вопрос для фирм, заинтересованных в удержании клиентов, связан с тем, насколько рентабельным может быть нынешний клиент в будущем. В этом примере авторы не фокусировались на прогнозе ожидаемой будущей прибыльности клиента или `CLV` (англ. "`Customer Lifetime Value`") или Пожизненная финансовая ценность клиента. Вместо этого попробуем сосредоточиться на драйверах `CLV`. Таким образом, авторы построили прогноз `CLV` для каждого из клиентов в выборке из 500. Затем авторы использовали это предсказание, чтобы понять, какая из переменных в данной базе данных поможет объяснить будущую ценность клиента. Если драйверы эффективно объясняют будущую ценность клиента, можно иметь возможность использовать результаты оценки в предсказании `CLV` для любого клиента, не входящего в текущую выборку. В конце этой главы авторы книги должны иметь возможность сделать следующее:

1. Определить драйверы `CLV`.

2. Предскажить ожидаемый `CLV` для каждого клиента.

3. Рассчитать предиктивную точность модели.

Информация, необходимая для этой модели, включает следующий список переменных:

| Переменные      | Описание                                                                             |
|:--------------- |:-------------------------------------------------------------------------------------|
| **Зависимая переменная** |  |
| CLV             | Дисконтированная стоимость всех ожидаемых будущих прибылей или стоимости жизни клиента|
| **Предикторы**          |  |
| Purchase_Rate | Средняя доля кварталов с покупками во всех 12 кварталах |
| Avg_Order_Quantity | Средняя долларовая стоимость покупок за все 12 кварталов |
| Avg_Crossbuy | Среднее значение для кросс-покупки за все 12 кварталов |
| Avg_Ret_Expense | Средние расходы, потраченные на маркетинговые усилия по удержанию данного клиента за все 12 кварталов |
| Avg_Ret_Expense_SQ | Квадрат средних расходов, потраченные на маркетинговые усилия по удержанию данного клиента за все 12 кварталов |
| Gender          | **1**, если клиент является мужчиной, **0**, если клиент является женщиной |
| Married         | **1**, если клиент женат, **0**, если клиент не состоял в браке |
| Income          | **1**, если доход < 30 000 долл. США, **2**; если 30 001 долл. США < доход < 45 000 долл. США; **3**, если 45 001 долл. США < доход < 60 000 долл. США; **4**, если 60 001 долл. США < доход < 75 000 долл. США; **5**, если 75 001 долл. США < доход < 90 000 долл. США; **6**, если доход > 90 001 долл. США |
| First_Purchase  | Стоимость первой покупки, сделанной клиентом в 1 квартале |
| Loyalty         | **1**, если клиент является членом программы лояльности, **0**, если нет |
| Share-of-Wallet (SOW) | Процент покупок клиента от данной фирмы, учитывая общий объем покупок по всем фирмам в этой категории |

```{r Ch04 : Customer retention - Customer Lifetime Value (CLV) data}
# Calculation Customer Lifetime Value (CLV) Database

library('sqldf')
customerRetentionCLV <- sqldf("select customer, avg(purchase) as
                              purchase_rate, avg(order_quantity) as avg_order_quantity,
                              avg(crossbuy) as avg_crossbuy, avg(ret_expense) as avg_ret_exp,
                              (avg(ret_expense) * avg(ret_expense)) as avg_ret_exp_sq, gender,
                              married, income, first_purchase, sow, clv, loyalty
                              from customerRetentionRepurchase group by customer, gender,
                              married, income, first_purchase, sow, clv, loyalty order by customer;")

ggplot(customerRetentionCLV, aes(x = clv)) +
  geom_histogram(bins = 20, col = "gray", aes(fill = ..count..)) +
    scale_fill_gradient("Count", low = "green", high = "red") +
      geom_density(adjust = 1/5, alpha = .2, fill = "#FF6666") +
         geom_vline(xintercept = mean(customerRetentionCLV$clv), linetype="dashed", color = "coral")

# Description of an empirical distribution for non-censored data using Cullen & Frey graph
library('fitdistrplus')
descdist(customerRetentionCLV$clv, boot = 500, discrete = FALSE)
min(customerRetentionCLV$clv) # clv < 0 !!!!!, but values must be in [0-1] to fit a beta distribution
(x <- fitdist(((customerRetentionCLV$clv - min(customerRetentionCLV$clv)) /
                (max(customerRetentionCLV$clv) - min(customerRetentionCLV$clv))),
             "beta", method = "mme", discrete = FALSE))
plot(x)
writeLines("\n Beta distribution of `customerRetentionCLV$clv`")
gofstat(x)

# (x <- fitdist(customerRetentionCLV$clv, "gamma", method = "mme", discrete = FALSE))
# plot(x)
# writeLines("\n Gamma distribution of `customerRetentionCLV$clv`")
# gofstat(x)
# 
# library('actuar') # 19 continuous tailed distributions (Burr, Gumbel, Inverse..., Paralogistic, Pareto, etc)
# (x <- fitdist(((customerRetentionCLV$clv - min(customerRetentionCLV$clv)) /
#                 (max(customerRetentionCLV$clv) - min(customerRetentionCLV$clv))),
#               "pareto", start = list(shape = 1, scale = 500)))
# plot(x)
# writeLines("\n Pareto distribution of `customerRetentionSOW$clv`")
# gofstat(x)

# # fitting a non-standard distribution to data and initial values
# set.seed(2018); x <- rweibull(500, shape = .5, scale = 1)
# x = (x - min(x) + 1) / (max(x) + 1)
# descdist(x, boot = 500, discrete = FALSE)
#  
# fitLog <- fitdist(x, "logis")
# fitln <- fitdist(x, "lnorm")
# fitW <- fitdist(x, "weibull")
# fitg <- fitdist(x, "gamma")
# fitn <- fitdist(x, "norm")
# fitexp <- fitdist(x, "exp")
# fitB <- fitdist(x, "beta", method = "mme")
# fitUni <- fitdist(x, "unif")
# fitP <- fitdist(x, "pareto", start = list(shape = 1, scale = 500))
# 
# (z <- gofstat(list(fitLog, fitln, fitW, fitg, fitn, fitexp, fitB, fitUni, fitP),
#               fitnames = c("Logistic", "Log-normal", "Weibull", "Gamma", "Normal", "Exponential", "Beta",
#                            "Uniform", "Pareto")))
# writeLines("\n Reference distribution of Sample is defined as:")
# which.min(z$cvm)

```

В этом случае имеется зависимая переменная `CLV` (ее гистограмма приведена на графике выше), которая представляет ожидаемую дисконтированную будущую прибыль от каждого клиента, распределенную по непрерывному бета-распределению. Переменная `CLV` непрерывна и не связана никаким порогом (и даже меньше нуля - AR). Таким образом, учитывая предположение, что `CLV` *нормально* распреден вокруг определенного среднего, авторы просто предложили построить линейную регрессию методом наименьших квадратов, чтобы выявить драйверы `CLV`. Поэтому авторы предлагают уравнение линейной регрессии в следующем формате:

$$ \bf CLV_i = X' \boldsymbol \beta + \varepsilon, \enspace \varepsilon = \bf \textit{ N } \boldsymbol (0, \enspace \sigma^2 ),  $$

где `CLV` - значение жизненного цикла клиента для данного клиента,

`X` - матрица независимых переменных,

$\beta$ - вектор коэффициентов независимых переменных,

$\varepsilon$ - случайная ошибка модели, который *нормально* распределен со средним значением 0 и дисперсией $\sigma^2$.

### Построение и верификация модели *Пожизненной финансовой ценности клиента* {.tabset}

При оценке модели `CLV` авторы получили следующие результаты:

#### Linear

```{r Ch04 : CLV - Linear Regression}
# Fit Linear Regression Model for Customer Lifetime Value (CLV) by Authors

(Ch04.clv <- lm(clv ~ purchase_rate + avg_order_quantity + avg_crossbuy + avg_ret_exp + avg_ret_exp_sq + 
                gender + married + income + first_purchase + loyalty + sow,
                data = customerRetentionCLV)) %>%
  summary

writeLines("Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04.clv) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity

```

Мы получаем следующие данные из этих результатов. Мы находим, что все переменные, за исключением `Purchase_Rate` и `Married`, статистически значимы при **p < 0,05**. Тем не менее сама модель имеет очень большую точность, но вот о ее устойчивости сказать сложно. Возможно имела место переподгонка под конкретные данные и требуется пройдти процедуры обеспечивающие устойчивость.

#### Linear (AR version)

Теперь попробуем исключить из линейной модели незначимые факторы `married` и `income` и провести построение устойчивой регресии посредством бустинга.
 
```{r Ch04 : CLV - Linear Regression Improved, warning=FALSE}
# Fit Linear Regression Model for Customer Lifetime Value (CLV)   (Improved)
uno = 'y_hat'

set.seed(2018)
(Ch04clv.AR <- train(clv ~ purchase_rate + avg_order_quantity + avg_ret_exp + avg_ret_exp_sq + 
                gender + first_purchase + loyalty + sow,
                data = customerRetentionCLV, method = "lm"))#,
                    # trControl = trainControl(method = "none", number = 1)))
summary(Ch04clv.AR)
# Coefficient Estimates and 95% CI
writeLines("\n Coefficient Estimates and 95% CI")
car::Confint(Ch04clv.AR$finalModel) %>%
  exp() %>%
    arm::pfround(digits = 3)

# Linear regression diagnostics
writeLines("\n Chisq test of predictors")
car::Anova(Ch04clv.AR$finalModel, type="II", test="Chisq")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04clv.AR$finalModel) # Variance Inflation Factors - if vif() > 2  - feature has multicollinearity
writeLines("\n Variable Importance for Model \n")
caret::varImp(Ch04clv.AR) %>% .$importance %>% print.AsIs()

# Create the scatter plots Linear Model versus model predictors
# https://stats.idre.ucla.edu/r/seminars/ggplot2_intro/
predictors <- Ch04clv.AR$trainingData %>% colnames(.) %>% .[-1]
Ch04clv.AR$trainingData %>% 
  rename(y_hat = `.outcome`) %>% 
    gather(key = "predictors", value = "predictor.value", -y_hat) %>%   
      ggplot(aes(predictor.value, y_hat))+
        geom_point(size = 0.5, alpha = 0.5) +
        geom_smooth(method = "loess") +
        facet_wrap(~ predictors, scales = "free_x") +
        ylab(stringr::str_to_title(uno))

# Plot matrix of statistical model diagnostics
GGally::ggnostic(Ch04clv.AR$finalModel, title = paste(paste(formula(Ch04clv.AR)[c(2, 1, 3)], collapse = " ")))

# wide variety of diagnostic plots for checking the quality of regression fit
car::influenceIndexPlot(Ch04clv.AR$finalModel)

# Goodness-of-fit of normal distributions to residuals of Linear Model (if Kolmogorov-Smirnov statistic < 0.05)
residuals(Ch04clv.AR) %>% fitdist("norm") %>% gofstat()

``` 

Улучшенная модель без признаков `avg_crossbuy`, `married` и `income` стала устойчивее и не практически не потеряла своей точности. Следует признать, что остатки линейной модели и в самом деле распределяются по нормальному распределению, поскольку значение статистики Колмогорова-Смирнова **p < 0.05**.

#### Error of CLV Models

Последний шаг - предсказать значение `Customer Lifetime Value (CLV)`, чтобы увидеть, насколько авторская и моя модели совпадают с фактическими значениями и в какой мере они превосходят *наивную* модель по среднему значению.

```{r Ch04 : Error of CVL Model}
# Computing the Mean Absolute Deviation (MAD) and Mean Absolute Percent Error (MAPE)

# mean_clv
y <- customerRetentionCLV$clv
writeLines(sprintf("Mean of Customer Lifetime Value (CLV): %.2f", mean(y)))

# Accuracy measures for Naïve AFT Model
m <- forecast::accuracy(rep(mean(customerRetentionCLV$clv), length(y)), y)
attributes(m)$ dimnames[[1]] <- "     Naïve Model (Mean) Accuracy: "; m

# Accuracy measures for Linear CLV Model 
y_hat <- predict(Ch04.clv, newdata = customerRetentionCLV)
m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "     Linear CLV Model's Accuracy: "; m

# Accuracy measures for Linear CLV Improved Model 
y_hat <- predict(Ch04clv.AR, newdata = customerRetentionCLV)
m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "Linear Improved Model's Accuracy: "; m

```

Среднее значение Пожизненная финансовая ценность клиента в данной выборке составляет `r sprintf("%.2f", mean(customerRetentionCLV$clv))` долл. Этот показатель авторы применили в качестве *наивного* прогноза. Созданная ими модель, а также моя модель *CLV* оказалось заметно аккуратнее *наивной* модели среднего: MAD = `r sprintf("%.2f", forecast::accuracy(y_hat, y)[, "MAE"])` долл. или в среднем `r sprintf("%.2f", forecast::accuracy(y_hat, y)[, "MAPE"])`% от фактического CLV.

### Как это использовать?

Надо отметить, что коэффициент `Avg_Order_Quantity` положителен, что указывает на то, что чем выше средние значения прошлых заказов клиента, тем выше CLV. Также оказалось, что в авторском модели коэффициент `Avg_Crossbuy` позитивен, поэтому следует предполагать, что чем больше клиент купил в нескольких категориях в прошлом, тем выше CLV клиента. Мы находим, что `Ret_Expense` положителен с уменьшающимся возвратом, как отмечено положительным коэффициентом на `Ret_Expense` и отрицательным коэффициентом на `Ret_Expense_SQ`. Это означает, что маркетинговые усилия, направленные на сохранение и установление отношений с клиентом, заставляют клиента иметь более высокий CLV. Затем, после достижения порога, маркетинговые усилия на самом деле уменьшают CLV в среднем. Вероятно, это связано с тем, что чрезмерное общение с клиентами может часто напрягать отношения между клиентом и фирмой, а поскольку, когда маркетинговые усилия теряют эффективность, затраты продолжают расти без соответствующего увеличения прибыли. Мы обнаруживаем, что характеристики клиента являются положительными (`Gender`, `First_Purshase` и `Loyality`), предполагая, что клиенты, которые являются мужчинами, с более высокой суммой первой покупки и являясь членами программы лояльности, скорее всего, будут иметь более высокую CLV.

Несколько сложнее понять чем вызван отрицательный коэффициент `Purchase_rate`, высокий уровень которого свидетельствует о высокой доли кварталов, в которые клиентом проводилась закупка. Однако невысокая значимость этого члена регрессионного уравнения в обоих моделях позволяет его проигнорировать.

## Дискуссия по RFM Модели *вероятности повторной покупки*

Наконец, рассмотрим модель *вероятности повторной покупки* для последнего - 12-го квартала с использованием **RFM** данных, рассчитанных за предыдущие 11 кварталов. Кроме трех предикторов: `recency_score`, `frequency_score` и `monetary_score` в пробит-модель включены расходы на удержание клиента в 12 квартале и их квадрат: `ret_expense` и `ret_expense_sq`.

```{r Ch04 : RFM Model - Purchase in Last Quarter, warning=FALSE}
# RFM: Recency, Frequency and Monetary Value Analysis - Customer Level Data
# https://rdrr.io/cran/rfm/f/vignettes/rfm-customer-level-data.Rmd
library('rfm') 

analysis_date <- lubridate::as_date((Sys.Date() - 91), tz = 'UTC')
CH04.rfm <- rfm_table_customer(data = customerRetentionTransactions %>%  
                              filter(quarter < 12 & purchase == 1) %>% 
                              group_by(customer) %>% 
                              summarise(n_transactions = sum(purchase), recency_days = (11-max(quarter)) * 91,
                                             total_revenue = sum(order_quantity) ), 
                              customer_id = customer, n_transactions = n_transactions,
                              recency_days = recency_days, total_revenue = total_revenue,
                              analysis_date = lubridate::as_date((Sys.Date() - 91), tz = 'UTC'))

# RFM Score Charts
rfm_heatmap(CH04.rfm)
rfm_histograms(CH04.rfm)

# Scatter Plots
rfm_rm_plot(CH04.rfm) # Recency vs Monetary Value
rfm_fm_plot(CH04.rfm) # Frequency vs Monetary Value
rfm_rf_plot(CH04.rfm) # Recency vs Frequency

# 11 Segments of Customers
segment <- c(
    "Champions", "Loyal Customers", "Potential Loyalist",
    "New Customers", "Promising", "Need Attention",
    "About To Sleep", "At Risk", "Can't Lose Them", "Hibernating",
    "Lost"
)
description <- c(
    "Bought recently, buy often and spend the most",
    "Spend good money. Responsive to promotions",
    "Recent customers, spent good amount, bought more than once",
    "Bought more recently, but not often",
    "Recent shoppers, but haven't spent much",
    "Above average recency, frequency & monetary values",
    "Below average recency, frequency & monetary values",
    "Spent big money, purchased often but long time ago",
    "Made big purchases and often, but long time ago",
    "Low spenders, low frequency, purchased long time ago",
    "Lowest recency, frequency & monetary scores"
)
recency <- c("4 - 5", "2 - 5", "3 - 5", "4 - 5", "3 - 4", "2 - 3", "2 - 3", "<= 2", "<= 1", "1 - 2", "<= 2")
frequency <- c("4 - 5", "3 - 5", "1 - 3", "<= 1", "<= 1", "2 - 3", "<= 2", "2 - 5", "4 - 5", "1 - 2", "<= 2")
monetary <- c("4 - 5", "3 - 5", "1 - 3", "<= 1", "<= 1", "2 - 3", "<= 2", "2 - 5", "4 - 5", "1 - 2", "<= 2")
segments <- tibble(
    Segment = segment, Description = description,
    R = recency, `F` = frequency, M = monetary
)


rfm_segments <- CH04.rfm %>%
    magrittr::use_series(rfm) %>%
    mutate(
        segment = case_when(
            (recency_score %>% between(4, 5)) & (frequency_score %>% between(4, 5)) &
                (monetary_score %>% between(4, 5))  ~ "Champions",
            (recency_score %>% between(2, 5)) & (frequency_score %>% between(3, 5)) &
                (monetary_score %>% between(3, 5)) ~ "Loyal Customers",
            (recency_score %>% between(3, 5)) & (frequency_score %>% between(1, 3)) &
                (monetary_score %>% between(1, 3)) ~ "Potential Loyalist",
            (recency_score %>% between(4, 5)) & (frequency_score == 1) &
                (monetary_score == 1) ~ "New Customers",
            (recency_score %>% between(3, 4)) & (frequency_score == 1) &
                (monetary_score == 1) ~ "Promising",
            (recency_score %>% between(2, 3)) & (frequency_score %>% between(2, 3)) &
                (monetary_score %>% between(2, 3)) ~ "Needs Attention",
            (recency_score %>% between(2, 3)) & (frequency_score <= 2) &
                (monetary_score <= 2) ~ "About To Sleep",
            (recency_score <= 2) & (frequency_score %>% between(2, 5)) &
                (monetary_score %>% between(2, 5)) ~ "At Risk",
            (recency_score == 1) & (frequency_score %>% between(4, 5)) &
                (monetary_score %>% between(4, 5)) ~ "Cant Lose Them",
            (recency_score %>% between(1, 2)) & (frequency_score %>% between(1, 2)) &
                (monetary_score %>% between(1, 2)) ~ "Hibernating",
            (recency_score <= 2) & (frequency_score <= 2) &
                (monetary_score <= 2) ~ "Lost",
            TRUE ~ "Others"
        )
    ) %>%
    dplyr::select(
        customer_id, segment, rfm_score, transaction_count, recency_days,
        amount
    )

# use datatable
rfm_segments %>%
    DT::datatable(
        filter = "top",
        options = list(pageLength = 5, autoWidth = TRUE),
        colnames = c(
            "Customer", "Segment", "RFM",
            "Orders", "Recency", "Total Spend"
        )
    )

rfm_top_segments <- rfm_segments %>%
  count(segment) %>%
    arrange(desc(n)) %>%
      rename(Segment = segment, Count = n)

z <- bind_cols(CH04.rfm$rfm,
      filter(customerRetentionTransactions, quarter == 12) %>%
        dplyr::select(purchase, order_quantity, ret_expense, ret_expense_sq))

# Fit Probit Regression Model for Customer Repurchase (Recency, Frequency and Monetary Value Analysis)
set.seed(2018)
uno = 'probit'
Ch04RFM.pr <- train(factor(purchase) ~ # recency_score + frequency_score + monetary_score + 
                      factor(rfm_score) + ret_expense + ret_expense_sq, 
                    metric = 'Kappa', data = z, method = "glm", family = binomial(link = uno))

summary(Ch04RFM.pr)
# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")
car::Confint(Ch04RFM.pr$finalModel) %>%
  exp() %>%
    arm::pfround(digits = 3)

# Logistic regression diagnostics
writeLines("\n Wald test of predictors")
car::Anova(Ch04RFM.pr$finalModel, type="II", test="Wald")
writeLines("\n Variance Inflation Factors - if vif() > 2  - feature has multicollinearity with others")
car::vif(Ch04RFM.pr$finalModel) # Variance Inflation Factors - if vif() > 2  - feature has multicollinearity
writeLines("\n Variable Importance for Model \n")
caret::varImp(Ch04RFM.pr) %>% .$importance %>% print.AsIs()

# # Create the scatter plots Probit versus model predictors
# predictors <- Ch04RFM.pr$coefnames
# Ch04RFM.pr$trainingData %>%
#   dplyr::select(one_of(predictors)) %>%
#     mutate(link = predict(Ch04RFM.pr$finalModel, newdata = z, type = "link")) %>%
#       gather(key = "predictors", value = "predictor.value", -link) %>%
#         ggplot(aes(predictor.value, link))+
#           geom_point(size = 0.05, alpha = 0.05) +
#           geom_smooth(method = "loess") +
#           facet_wrap(~ predictors, scales = "free_x") +
#           ylab(stringr::str_to_title(uno))
# 
#
# # Plot matrix of statistical model diagnostics
# GGally::ggnostic(Ch04RFM.pr$finalModel, title = paste(paste(formula(Ch04RFM.pr)[c(2, 1, 3)], collapse = " ")))
#
# # wide variety of diagnostic plots for checking the quality of regression fit
# # https://bookdown.org/jefftemplewebb/IS-6489/logistic-regression.html
# car::influenceIndexPlot(Ch04RFM.pr$finalModel)

caret::confusionMatrix(data = predict((Ch04RFM.pr), newdata = z),
                       reference = z$purchase %>% factor,
                       positive = "1", mode = "everything")

```

Ниже представлены гистограмма сегментов в полярных координатах и классический вариант графика "**Итогового RFM Отчета**".

```{r RFM Plotting, warning=FALSE}
d = data_frame(ID = c(11, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1),
               x1 = c(0, 2, 3, 4, 3, 2, 2, 0, 0, 1, 4), x2 = c(2, 5, 5, 5, 4, 3, 3, 2, 1, 2, 5), 
               y1 = c(0, 3, 1, 0, 0, 2, 0, 2, 4, 1, 4), y2 = c(2, 5, 3, 1, 1, 3, 2, 5, 5, 2, 5),
               t = c('Lost', 'Loyal Customers', 'Potential Loyalist', 'New Customers', 'Promising',
                 'Needs Attention', 'About To Sleep', 'At Risk', 'Can`t Lose Them', 'Hibernating', 'Champions'))

segmentcolors <- c('slategray3', 'lightgoldenrod1', 'violet', 'chartreuse', 'lightcyan2',                                'cadetblue3', 'darkseagreen1', 'turquoise3', 'pink', 'darkslategray1', 'lightskyblue2')
             
d <- d %>% 
        dplyr::full_join(rfm_top_segments, by = c("t" = "Segment")) %>%
          mutate(Count = coalesce(Count, 0L), a = Count / sum(Count),
                 l = paste0(t, if_else(Count == 0, sprintf(" %.0f", Count),
                                       sprintf("\n%.0f (%.1f%%)", Count, a * 100))))  %>% 
            arrange(ID)

## 1. Plotting of Correlation Matrix for Feartures
library("GGally")
RFMLabels <- c("Recency", "Frequency", "Monetary")
p <- z %>% 
  dplyr::select(recency_score:monetary_score) %>% 
  GGally::ggpairs(., xlab = "Scores", ylab = "Scores",
                  title = "Correlations Between Scores",
                  lower = list(continuous = wrap("smooth_loess", alpha = 0.3),
                               combo = wrap("dot_no_facet", alpha = 0.4)),
                  columnLabels = RFMLabels)

p[2, 1] <- z %>% # Recency vs. Frequency
    ggplot(aes(x = recency_score, y = frequency_score)) + 
      stat_smooth(method = 'loess', col = 'coral') + 
      geom_jitter(aes(color = rfm_segments$segment), width = 0.2, height = 0.2) +
      scale_color_manual(values = segmentcolors, breaks = d$l)

p[3, 1] <- z %>% # Recency vs. Monetary
    ggplot(aes(x = recency_score, y = monetary_score)) + 
      stat_smooth(method = 'loess', col = 'dodgerblue') + 
      geom_jitter(aes(color = rfm_segments$segment), width = 0.2, height = 0.2) +
      scale_color_manual(values = segmentcolors, breaks = d$l)

p[3, 2] <- z %>% # Frequency vs. Monetary
    ggplot(aes(x = frequency_score, y = monetary_score)) + 
      stat_smooth(method = 'loess', col = 'gold') + 
      geom_jitter(aes(color = rfm_segments$segment), width = 0.2, height = 0.2) +
      scale_color_manual(values = segmentcolors, breaks = d$l)

p[1, 1] <- z %>% 
    ggplot(aes(recency_score)) +
      stat_density(fill = 'coral')

p[2, 2] <- z %>% 
    ggplot(aes(frequency_score)) +
      stat_density(fill = 'dodgerblue')

p[3, 3] <- z %>% 
    ggplot(aes(monetary_score)) +
      stat_density(fill = 'gold')

p

# Compute hierarchical clustering of clients
K = 5
cln.hc <- z %>% filter(frequency_score > 0) %>%
  # scale() %>%                         # Scale the data
  dist(method = "euclidean") %>%        # Compute dissimilarity matrix
  hclust(method = "ward.D")             # Compute hierachical clustering
z <- z %>%  mutate(clusters = cutree(cln.hc, k = K))

# Print of The Hierarchical Clustering of Clients by RFM Model
library('kableExtra') # Construct Complex Table with 'kable' and Pipe Syntax 
library('formattable') # Create 'Formattable' Data Structures
bind_cols( z %>% group_by(clusters) %>% 
               summarise(count = n()), 
           z %>% group_by(clusters) %>% 
               summarise_at(.vars = vars(recency_score:monetary_score),
                            .funs = c(Mean = "mean", Sd = "sd")) %>%
             .[, -1]) %>% 
  mutate(clusters = cell_spec(clusters, "html",
       color = ifelse(count <= (arrange(., count)[1, 'count'] %>%  as.numeric), "red", "auto")),
       count = color_bar("thistle")(accounting(count, digits = 0)),
       recency_score_Mean = color_tile("ivory1", "coral")(digits(recency_score_Mean, digits = 2)),
       frequency_score_Mean = color_tile("thistle1", "dodgerblue")(digits(frequency_score_Mean, digits = 2)),
       monetary_score_Mean = color_tile("chartreuse", "gold")(digits(monetary_score_Mean, digits = 2)),
       recency_score_Sd, frequency_score_Sd, monetary_score_Sd
      ) %>% 
    kable(format = "html", digits = c(0, 0, 2, 2, 2, 2, 2, 2), longtable = TRUE, booktabs = TRUE, escape = F,
          col.names = c("No Clusters", "Counts of Clients", "Recency Scores", "Frequency Scores",
                        "Monetary Scores", "Recency Scores", "Frequency Scores", "Monetary Scores"),
          caption="The Hierarchical Clustering of Clients by RFM Model") %>%
      kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", full_width = FALSE)) %>%
        column_spec(2, width = "4 cm") %>%
        add_header_above(c(" ", " ", "Mean Of" = 3, "Standard Deviation Of" = 3)) %>% 
        group_rows(index = c("Low" = 2, "High" = 2, "Mid" = 1))

# 2. Plotting of Clients Dendrogram by RFM scores in K groups and color by Clusters
library('factoextra') # Extract and Visualize the Results of Multivariate Data Analyses
dend <- fviz_dend(cln.hc, k = K, # Cut in K Clusters
         cex = 0.5, # label size
         k_colors = c("deepskyblue", "deepskyblue3", "yellow3", "goldenrod", "gold"),
         show_labels = FALSE,
         color_labels_by_k = TRUE, # color labels by groups
         main = "Dendrogram of clusters", xlab = "Clients", ylab = "Distance",
         rect = TRUE) # Add rectangle around groups
dend + annotate("text", x = c(35, 140, 310, 430, 470), y = rep(-15000, 5),
                label = c("Five\ncluster", "Three\ncluster", "Second\ncluster", "Four\nclus.", "First\ncluster"))

# 3. Plotting Customer Clusters by Interactive Web Graphics using plotly 
library('plotly')
plot_ly(data = z, x = z$recency_score + rnorm(500) / 8, y = z$monetary_score + rnorm(500) / 8,
  z = z$frequency_score  + rnorm(500)/8, type = "scatter3d", mode = "markers", color = cutree(cln.hc, k = K)) %>% 
  layout(title = "Customer clustering",
             scene = list(xaxis = list(title = "Recency"),
                          yaxis = list(title = "Monetary"),
                          zaxis = list(title = "Frequency"),
                          camera = list(eye = list(x = 1, y = 2, z = -1),
                                        center = list(x = 0, y = 0, z = 0)))
)

# 4. Plotting The RFM Polar Report
ggplot(d, aes(x = ID, y = Count, fill = l)) +
  geom_bar(width = 1, stat = "identity", color = "white") + 
  scale_fill_manual(values = segmentcolors, name = "Segments", breaks = rev(d$l),
                    guide = guide_legend(label.theme =  element_text(size = 8, angle = 0))) +
  theme(axis.ticks = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.line = element_blank()) +
  geom_text(data = filter(d, Count > 0), aes(label = l), lineheight = 0.8, size = 3, show.legend = FALSE) +
  coord_polar() +
  ggtitle("RFM Polar Report")
  
# 5. Plotting The RFM Summary Report
ggplot() + 
  geom_rect(data = d, aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = t), color = "dark gray",
            alpha = if_else(d$a == 0, 1 / max(d$a) * min(d$a[d$a > 0]) - 0.01, 1 / max(d$a) * d$a) ) +
  geom_text(data = d, aes(x = x1 + (x2 - x1) / 2, y = y1 + (y2 - y1) / 2 + if_else(t == 'Lost', -0.3, 0),
                          label = l), size = 4) +
  scale_fill_manual(values = segmentcolors, name="Segments") + 
  labs(title = "RFM Summary Report", x = "Recency Score", y = "Frequency Score") +
  theme_classic() + theme(legend.position="none")

```

Как видно во многих категориях клиенты не представлены. И если отсутствие в группах 'At Risk', 'Can`t Lose Them' или 'Lost' не вызывает опасений, то нули в 'New Customers' либо 'Promising' должно настораживать.

### Как это использовать?

Основная задача менеджмента переводить клиентов из голубых и розовых зон в зеленную зону.

Что касается новых (розовых) клиентов, то важнейший пункт - *вернуть покупателей за второй покупкой* — обратиться к ним как можно быстрее. Но не слишком рано и весьма деликатно. Что включить в [содержание письма](https://retailrocket.ru/blog/cases/rfm-segmentatsiya-keys-tehnosila/):

* **Текст**: «Поздравляем с первой покупкой»

* **Подборка**: «Персональные товарные рекомендации»

* **Баннер**: «Знакомство с категориями»

* **Блок**: «Описание Программы лояльности»

Некоторыми шагами могут стать:

1. **Реактивация клиента** (для фиолетовых и красных клиентов, склонных к уходу в отток), направлять 1 раз в месяц. "Мы соскучились. Вот вам суперскидка на 3 дня".

2. **Клиентам, которым необходимо внимание** направить специальное предложение, котоырей действует 1-2 дня.

3. Для **перспективных клиентов** голубой и бледно-зеленой зоны направить сообщение "Вам может понравиться этот товар", а также периодически рассылать обзор "Знакомство с категориями", "Знакомство с новинками".

4. Для зеленых, самых лояльных клиентов, предложение дать рекомендациии фирме.

Триггеры для всех клиентов:

* **брошенная корзина**, положил товар в корзину, но через 30 минут так ничего  не купил;

* **брошенный просмотр**, смотрел товары, но ничего в корзину не положил;

* **прошенный поиск**, искал на сайте, придя по поисковику или каталогу, но ничего не положил в корзину.



```{r The End of session}

devtools::session_info()

```

## Выводы по главе 4

Цель этой главы состояла в том, чтобы изучить эконометрические модели удержания клиентов и предоставить некоторые эмпирические примеры того, как фирмы могут применять эти знания к своим собственным базам данных клиентов. И авторы показали, что когда фирмы могут понять каждый аспектов удержания клиентов, то они могут эффективно управлять своими текущими клиентами для получения прибыли.