Statistical Methods in Customer Relationship Management
Chapter 3. Customer acquisition
Alexander Rodionov
31 мая 2018 г.
- Глава 3. Привлечение потребителей
Мы начинаем рассматриваем в R примеры SAS на наборах данных из книги Kumar V., Petersen Andrew J. “Statistical Methods in Customer Relationship Management”.
Глава 3. Привлечение потребителей
В этой главе начинается обсуждение первого этапа CRM - приобретения новых клиентов. Цели моделирования привлечения клиентов включают в себя определение правильных клиентов, прогнозирование того, будут ли клиенты реагировать на рекламные кампании компании, планирование количества новых клиентов и анализирование краткосрочных и долгосрочных последствий маркетинговых и других бизнес-переменных при получении таких клиентов. Спецификации модели, которые рассматриваются в этой главе, включают логит, пробит, Тобит, линейную регрессию, лог-линейную модель, векторную авторегрессию, функцию риска (англ. “hazard function”) и исчисление решений (англ. “decision calculus”). Помимо изучения аспектов моделирования привлечения клиентов, в этой главе авторы также рассматривают большое количество исследований, в которых основное внимание уделяется влиянию маркетинговых переменных в качестве драйверов или предикторов приобретения клиентов.
На рисунке представлены вопросы, которые необходимо решить в эконометрических и статистических моделях удержания клиентов.
The Picture 1
Решение этих вопросов мы вместе с авторами книги рассмотрим ниже.
Данные для эмпирических примеров
Имеются набор данных “customerAcquisition” о 500 потенциальных потребителей типичной B2B компании, содержащие 16 признаков:
| Переменные | Описание |
|---|---|
| Customer | Номер потребителя (от 1 до 500) |
| Acquisition | 1, если клиент был приобретен, 0 в противном случае |
| First_Purchase | Долларовая стоимость первой покупки (0, если клиент не был приобретен) |
| CLV | Прогноз пожизненной ценности клиента. Это 0, если клиента не было или уже ушел в отток (’000) |
| Duration | Время в днях, когда компании была или продолжает быть клиентом, цензуированная до 730 дней |
| Censor | 1, если клиент оставался в конце окна наблюдения, 0 в противном случае |
| Acq_Expense | Долларовые расходы на маркетинг по привлечению клиента |
| Acq_Expense_SQ | Квадрат расходов на привлечение |
| Industry | 1, если клиент в секторе B2B, 0 в противном случае, т. е. B2C |
| Revenue | Годовой доход от продаж компании (млн долл.) |
| Employees | Количество сотрудников в компании |
| Ret_Expense | Долларовые расходы на маркетинг по удержание клиента |
| Ret_Expense_SQ | Квадрат расходов на удержание |
| Crossbuy | Количество категорий товаров / услуг, приобретенных клиентом |
| Frequency | Количество покупок клиента во время окна наблюдения |
| Frequency_SQ | Квадрат количества покупок |
Авторы настаивают на использовании как линейных, так и квадратичных величин расходов, поскольку ожидают, что для каждого дополнительного доллара, потраченного на усилия для приобретения и удержания для данной потенциального клианта, будет проявлятся закон убывающей доходности (англ. “Diminishing Returns”) выразающийся в том, что дополнительные затраты дают всё меньший объём дополнительной выручки.
library('tidyverse')## -- Attaching packages ----------------------------------------------------------------------------- tidyverse 1.2.1 --## v ggplot2 2.2.1 v purrr 0.2.4
## v tibble 1.4.2 v dplyr 0.7.5
## v tidyr 0.8.1 v stringr 1.3.1
## v readr 1.1.1 v forcats 0.3.0## -- Conflicts -------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library('caret')## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftlibrary('car')## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someutils::data('customerAcquisition', package = 'SMCRM')
# Check for Class Imbalances
writeLines("Distribution Variable 'Acquisition' by Response Levels")## Distribution Variable 'Acquisition' by Response LevelscustomerAcquisition$acquisition %>%
factor() %>%
table() # %>%## .
## 0 1
## 208 292 # prop.test()Эмпирический пример: вероятность отклика на привлечение клиента
Построим регрессионую модель с логистическим распределением отклика по набору предикторов согласно авторам книги. Обратим внимание, что имеется определенная несбалансированность классов, т.е. компаний, которые стали клиентами, заметно больше, чем отказавшихся ими стать.
Логистическая модель обычно представляется как:
\[ \displaystyle \large \pi(Y)=\frac{\exp(\beta_0+\beta_1X)}{1+\exp(\beta_0+\beta_1X)} \hspace{.5 in} [1]\] или переходя в обычной линейной модели:
\[ \displaystyle \large \ln\left(\frac{\pi(Y)}{1-\pi(Y)}\right)=\beta_0+\beta_1X \hspace{.5 in} [2]\]
Поэтому сам логит получается:
\[ \displaystyle \large \Pr(Y=1 \mid X) = [1 + e^{-X'\beta}]^{-1} \hspace{.5 in} [3]\]
Построение и верификация модели привлечения клиента (первый этап)
Logit
# Fit Logistic Regression Model for Customer Acquisition by Authors
# https://stats.idre.ucla.edu/sas/dae/logit-regression & https://stats.idre.ucla.edu/r/dae/logit-regression
Ch03.logit <- glm(acquisition ~ acq_expense + acq_expense_sq + industry + revenue + employees,
data = customerAcquisition, family = binomial(link = 'logit'))
summary(Ch03.logit)##
## Call:
## glm(formula = acquisition ~ acq_expense + acq_expense_sq + industry +
## revenue + employees, family = binomial(link = "logit"), data = customerAcquisition)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.48550 -0.06183 0.03275 0.26053 1.95191
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.621e+01 3.537e+00 -7.408 1.28e-13 ***
## acq_expense 6.720e-02 1.036e-02 6.488 8.69e-11 ***
## acq_expense_sq -4.138e-05 7.989e-06 -5.180 2.21e-07 ***
## industry 3.272e-02 3.868e-01 0.085 0.93259
## revenue 3.186e-02 1.182e-02 2.695 0.00703 **
## employees 4.959e-03 7.348e-04 6.749 1.48e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 678.97 on 499 degrees of freedom
## Residual deviance: 194.46 on 494 degrees of freedom
## AIC: 206.46
##
## Number of Fisher Scoring iterations: 8writeLines(sprintf("-2 Log L of Intercept and Only Covariates: %.3f", -2 * logLik(Ch03.logit)[1]))## -2 Log L of Intercept and Only Covariates: 194.464writeLines(sprintf(" AIC (smaller is better): %.3f", extractAIC(Ch03.logit)[2]))## AIC (smaller is better): 206.464# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")##
## Odds Ratio Estimates and 95% CIcar::Confint(Ch03.logit) %>%
exp() %>%
arm::pfround(digits = 3)## Estimate 2.5 % 97.5 %
## (Intercept) 0.000 0.000 0.000
## acq_expense 1.070 1.049 1.093
## acq_expense_sq 1.000 1.000 1.000
## industry 1.033 0.481 2.208
## revenue 1.032 1.009 1.057
## employees 1.005 1.004 1.007writeLines("\n Wald test of predictors")##
## Wald test of predictorscar::Anova(Ch03.logit, type="II", test="Wald")## Analysis of Deviance Table (Type II tests)
##
## Response: acquisition
## Df Chisq Pr(>Chisq)
## acq_expense 1 42.0957 8.692e-11 ***
## acq_expense_sq 1 26.8366 2.214e-07 ***
## industry 1 0.0072 0.932592
## revenue 1 7.2657 0.007028 **
## employees 1 45.5540 1.485e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1writeLines("\n Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with others")##
## Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with otherscar::vif(Ch03.logit) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity ## acq_expense acq_expense_sq industry revenue employees
## 60.292660 50.908495 1.030094 1.156517 2.565201writeLines("\n Logit Model: Association of Predicted Probabilities and Observed Responses \n")##
## Logit Model: Association of Predicted Probabilities and Observed Responsesprob <- predict(Ch03.logit, newdata = customerAcquisition, type = "response")
caret::confusionMatrix(data = ifelse(prob > 0.5, "1", "0") %>% factor,
reference = customerAcquisition$acquisition %>% factor,
positive = "1", mode = "everything")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 183 19
## 1 25 273
##
## Accuracy : 0.912
## 95% CI : (0.8837, 0.9353)
## No Information Rate : 0.584
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8181
## Mcnemar's Test P-Value : 0.451
##
## Sensitivity : 0.9349
## Specificity : 0.8798
## Pos Pred Value : 0.9161
## Neg Pred Value : 0.9059
## Precision : 0.9161
## Recall : 0.9349
## F1 : 0.9254
## Prevalence : 0.5840
## Detection Rate : 0.5460
## Detection Prevalence : 0.5960
## Balanced Accuracy : 0.9074
##
## 'Positive' Class : 1
## В качестве меры точности модели мы рекомендует коэффициент Cohen’s Kappa, так как в задаче наблюдается проблема несбалансированности классов. Коэффициент аккуратности из-за этой проблемы тут не применим. Коэффициент Kappa - мера того, насколько близко наблюдения, классифицированные входе машинного обучения, т.е. предсказанные метки классов, соответствуют истинному разбиению на классы, контролируя точность обученного классификатора. Более того, что показатель Kappa проливает свет на то, как обучился сам классификатор, коэффициент для одной модели напрямую сопоставим со коэффициентами для любой другой модели, применяемой для одной и той же задачи классификации, включая классификацию по многим классам.
В статистической практике решения классификационных задач считается , что Kappas > 0,80 характеризует отличную модель; интервал от 0,40 до 0,80, справедливо относить к хорошему классификатору, а < 0,40 - ко слабым в прогностическом отношении результатам, хотя их и можно применять для описательных целей.
Теперь обратимся к пробит-модели, которую авторы выбралди в качестве основной на первом этапе этого моделирования.
Probit
Подобно модели логита, модель пробит (англ. “Probit”), основанная на нормальном распределении, часто используется для моделирования двоичных откликов, особенно в тех случаях, когда есть желание оценить модель двух потоков. Модель пробита, еще более сложная для оценки из-за отсутствия решения с замкнутой формой, теоретически привлекательна в двухэтапной структуре моделирования из-за стандартного нормального распределения. Авторы пишут, что это делает его полезным, когда второй этап представляет собой линейную регрессию строемую методом наименьших квадратов с нормально распределенным остатками.
Probit: \(\Pr(Y=1 \mid X) = \Phi(X'\beta)\)
# Fit Probit Model for Customer Acquisition by Authors
# https://stats.idre.ucla.edu/sas/dae/probit-regression & https://stats.idre.ucla.edu/r/dae/probit-regression
Ch03.probit <- glm(acquisition ~ acq_expense + acq_expense_sq + industry + revenue + employees,
data = customerAcquisition, family = binomial(link = 'probit'))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(Ch03.probit)##
## Call:
## glm(formula = acquisition ~ acq_expense + acq_expense_sq + industry +
## revenue + employees, family = binomial(link = "probit"),
## data = customerAcquisition)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.47116 -0.01729 0.00416 0.23784 1.92274
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.513e+01 1.862e+00 -8.127 4.42e-16 ***
## acq_expense 3.891e-02 5.431e-03 7.166 7.74e-13 ***
## acq_expense_sq -2.409e-05 4.160e-06 -5.791 6.98e-09 ***
## industry 5.135e-02 2.154e-01 0.238 0.81161
## revenue 1.771e-02 6.608e-03 2.681 0.00735 **
## employees 2.853e-03 4.032e-04 7.077 1.48e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 678.97 on 499 degrees of freedom
## Residual deviance: 192.12 on 494 degrees of freedom
## AIC: 204.12
##
## Number of Fisher Scoring iterations: 9writeLines(sprintf("-2 Log L of Intercept and Only Covariates: %.3f", -2 * logLik(Ch03.probit)[1]))## -2 Log L of Intercept and Only Covariates: 192.122writeLines(sprintf(" AIC (smaller is better): %.3f", extractAIC(Ch03.probit)[2]))## AIC (smaller is better): 204.122writeLines("\n Wald test of predictors")##
## Wald test of predictorscar::Anova(Ch03.probit, type="II", test="Wald")## Analysis of Deviance Table (Type II tests)
##
## Response: acquisition
## Df Chisq Pr(>Chisq)
## acq_expense 1 51.3462 7.743e-13 ***
## acq_expense_sq 1 33.5400 6.981e-09 ***
## industry 1 0.0568 0.811607
## revenue 1 7.1860 0.007347 **
## employees 1 50.0807 1.476e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1writeLines("\n Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with others")##
## Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with otherscar::vif(Ch03.probit) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity ## acq_expense acq_expense_sq industry revenue employees
## 62.353531 53.432649 1.034644 1.137440 2.453772writeLines("\n Probit Model: Association of Predicted Probabilities and Observed Responses \n")##
## Probit Model: Association of Predicted Probabilities and Observed Responsesprob <- predict(Ch03.probit, newdata = customerAcquisition, type = "response")
caret::confusionMatrix(data = ifelse(prob > 0.5, "1", "0") %>% factor,
reference = customerAcquisition$acquisition %>% factor,
positive = "1", mode = "everything")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 183 20
## 1 25 272
##
## Accuracy : 0.91
## 95% CI : (0.8814, 0.9336)
## No Information Rate : 0.584
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8141
## Mcnemar's Test P-Value : 0.551
##
## Sensitivity : 0.9315
## Specificity : 0.8798
## Pos Pred Value : 0.9158
## Neg Pred Value : 0.9015
## Precision : 0.9158
## Recall : 0.9315
## F1 : 0.9236
## Prevalence : 0.5840
## Detection Rate : 0.5440
## Detection Prevalence : 0.5940
## Balanced Accuracy : 0.9057
##
## 'Positive' Class : 1
## Иначе выражаясь, логистическая функция имеет слегка более плоские хвосты у нуля и единицы, тогда как кривая пробита приближается к осям резче, чем логит.
Логистическую регрессию лучше интерпретитовать, чем пробит, поскольку она может рассматриваться как логарифмические коэффициенты отношений шансов.
Probit (AR version)
Вместе с тем очевидно, что предиктор industry не значим в обоих бинарных регрессиях. Хотя авторы в примечании на странице 41 упоминают об этом, но никаких шагов по устранению этого фактора не предпринимают.
Более того, проверка по VIF (англ. “Variance Inflation Factors”) демонстрирует большие значения по признакам acq_expense и acq_expense_sq. Однако, принимая во внимание проявление закона убывающей доходности мы должны оставить оба мультиколлинеарных предиктора в этой модели. Теперь бинарную модель мы проверяем более тщательно.
# Fit Probit Regression Model for Customer Acquisition (Improved)
uno = 'probit'
set.seed(2018) #From random.org
(Ch03pr.AR <- train(factor(acquisition) ~ acq_expense_sq + revenue + employees, metric = 'Kappa',
data = customerAcquisition, method = "glm", family = binomial(link = uno)))#,## Generalized Linear Model
##
## 500 samples
## 3 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 500, 500, 500, 500, 500, 500, ...
## Resampling results:
##
## Accuracy Kappa
## 0.8967953 0.7890784 # trControl = trainControl(method = "none", number = 1)))
summary(Ch03pr.AR)##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8985 -0.3397 0.0024 0.3150 1.8403
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.788e+00 3.877e-01 -9.769 < 2e-16 ***
## acq_expense_sq 9.542e-06 8.204e-07 11.632 < 2e-16 ***
## revenue 1.180e-02 5.332e-03 2.214 0.0268 *
## employees 1.686e-03 2.176e-04 7.748 9.34e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 678.97 on 499 degrees of freedom
## Residual deviance: 264.66 on 496 degrees of freedom
## AIC: 272.66
##
## Number of Fisher Scoring iterations: 8# Odds Ratio Estimates and 95% CI
writeLines("\n Odds Ratio Estimates and 95% CI")##
## Odds Ratio Estimates and 95% CIcar::Confint(Ch03pr.AR$finalModel) %>%
exp() %>%
arm::pfround(digits = 3)## Estimate 2.5 % 97.5 %
## (Intercept) 0.023 0.010 0.047
## acq_expense_sq 1.000 1.000 1.000
## revenue 1.012 1.001 1.023
## employees 1.002 1.001 1.002# Logistic regression diagnostics
writeLines("\n Wald test of predictors")##
## Wald test of predictorscar::Anova(Ch03pr.AR$finalModel, type="II", test="Wald")## Analysis of Deviance Table (Type II tests)
##
## Response: .outcome
## Df Chisq Pr(>Chisq)
## acq_expense_sq 1 135.3033 < 2.2e-16 ***
## revenue 1 4.9008 0.02684 *
## employees 1 60.0303 9.341e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1writeLines("\n Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with others")##
## Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with otherscar::vif(Ch03pr.AR$finalModel) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity## acq_expense_sq revenue employees
## 1.348118 1.030692 1.313025writeLines("\n Variable Importance for Model \n")##
## Variable Importance for Modelcaret::varImp(Ch03pr.AR) %>% .$importance %>% print.AsIs()## Overall
## acq_expense_sq 100.00000
## revenue 0.00000
## employees 58.75996# # Create the scatter plots Logit versus model predictors
# # prob <- predict(Ch06lg.AR, newdata = customerAcquisition, type = "prob")[, "1"]
# # mutate(logit = log(prob / (1 - prob))) %>%
# Create the scatter plots Probit versus model predictors
predictors <- Ch03pr.AR$coefnames
Ch03pr.AR$trainingData %>%
dplyr::select(one_of(predictors)) %>%
mutate(link = predict(Ch03pr.AR$finalModel, newdata = customerAcquisition, type = "link")) %>%
gather(key = "predictors", value = "predictor.value", -link) %>%
ggplot(aes(predictor.value, link))+
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(Ch03pr.AR$finalModel, title = paste(paste(formula(Ch03pr.AR)[c(2, 1, 3)], collapse = " ")))## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
# wide variety of diagnostic plots for checking the quality of regression fit
# https://bookdown.org/jefftemplewebb/IS-6489/logistic-regression.html
car::influenceIndexPlot(Ch03pr.AR$finalModel)
writeLines("\n Improved Probit Model: Association of Predicted Probabilities and Observed Responses \n")##
## Improved Probit Model: Association of Predicted Probabilities and Observed Responsescaret::confusionMatrix(data = predict(Ch03pr.AR, newdata = customerAcquisition),
reference = customerAcquisition$acquisition %>% factor,
positive = "1", mode = "everything")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 191 30
## 1 17 262
##
## Accuracy : 0.906
## 95% CI : (0.877, 0.9301)
## No Information Rate : 0.584
## P-Value [Acc > NIR] : < 2e-16
##
## Kappa : 0.8083
## Mcnemar's Test P-Value : 0.08005
##
## Sensitivity : 0.8973
## Specificity : 0.9183
## Pos Pred Value : 0.9391
## Neg Pred Value : 0.8643
## Precision : 0.9391
## Recall : 0.8973
## F1 : 0.9177
## Prevalence : 0.5840
## Detection Rate : 0.5240
## Detection Prevalence : 0.5580
## Balanced Accuracy : 0.9078
##
## 'Positive' Class : 1
## qplot(`Observed Classes`, `Predicted Classes`,
data=bind_cols(`Observed Classes`= factor(customerAcquisition$acquisition),
`Predicted Classes` = predict(Ch03pr.AR, newdata = customerAcquisition)),
colour= `Observed Classes`, geom = c("boxplot", "jitter"),
main = "Predicted Classes vs. Observed Classes", xlab = "Observed Classes", ylab = "Predicted Classes")
Полученая улучшенная регрессия без фактора industry (B2B либо B2C) мне представляется более надежной, так как избавилась от малозначимого предиктора.
Как это использовать?
Авторы полагают, что независимо от того, используется логит или пробит-структура для моделирования вероятности ответа, вывод модели весьма полезен для определения того, какой клиент может приобрести фирма. Кроме того, результаты модели бинарного выбора могут служить драйверами приобретения клиентов, которые полезны для менеджеров при принятии решений в будущих кампаниях по привлечению клиентов.
Теперь можно узнать, как изменения в расходах на приобретение и изменения характеристик потенциальных клиентов могут либо увеличить, либо уменьшить нашу вероятность заполучить потребителя. Эта информация может дать значительную информацию руководителям, которым поручено определить оптимальный объем ресурсов, которые необходимо потратить на привлечение клиентов.
Эмпирический пример: величина первоначального заказа вновь привлеченных клиентов
Многие фирмы поняли, что недостаточно просто сосредоточиться только на попытке приобрести как можно больше клиентов, не заботясь о той ценности, которую клиент может предоставить. Исследования в области маркетинга показали, что первоначальная стоимость заказа может быть ценным предиктором для будущей стоимости клиента или, по крайней мере, оправдать сумму денег, затрачиваемую на привлечение клиента. Таким образом, может быть полезно понять драйверы значения первоначальная заказа и, в свою очередь, иметь возможность прогнозировать ожидаемое значение первоначального заказа каждого клиента.
Таким образом, авторы предлагают двухэтапную модель прогноза:
\[ \displaystyle \large E(Initial \enspace Order \enspace Quantity) = P(Acquisition = 1) * E(First \enspace Purchase | Acquisition = 1) \hspace{.5 in} [4] \] И если с моделью, по которой рассчитывается количествм вновь привлеченных клиентов, авторы уже определились остается второй этап.
Поскольку первый этап этой модели авторами уже пройден и они остановились на пробит-модели. Когда мы оцениваем второй этап двухэтапной модели, мы получаем нижеследующие оценки параметров линейной регрессии уравнений. Обращаем внимание, что в построении этой модели участвуют только привлеченные клиенты (Acquisition == 1).
Построение и верификация величины первоначального заказа (второй этап)
Linear
# Fit Linear Regression Model for Initial Order Quantity by Authors
# https://stats.idre.ucla.edu/sas/webbooks/reg/
# SAS Code: imr_acquisition = (pdf(’Normal’, xb_probit))/(probnorm(xb_probit));
xbeta <- predict(Ch03.probit, newdata = customerAcquisition, type = "link")
customerAcquisition <- customerAcquisition %>%
mutate(imr_acquisition = dnorm(xbeta) / pnorm(xbeta)) # Cumulative normal pdf
(Ch03.linear <- lm(first_purchase ~ acq_expense + acq_expense_sq + industry + revenue + employees + imr_acquisition, data = filter(customerAcquisition, acquisition == 1) )) %>%
summary##
## Call:
## lm(formula = first_purchase ~ acq_expense + acq_expense_sq +
## industry + revenue + employees + imr_acquisition, data = filter(customerAcquisition,
## acquisition == 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -79.737 -19.651 -0.775 19.425 98.388
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.226e+01 4.896e+01 -1.272 0.2046
## acq_expense 7.048e-01 1.303e-01 5.409 1.34e-07 ***
## acq_expense_sq -7.674e-04 9.221e-05 -8.323 3.62e-15 ***
## industry -2.764e+00 3.706e+00 -0.746 0.4564
## revenue 3.026e+00 1.115e-01 27.149 < 2e-16 ***
## employees 2.541e-01 4.764e-03 53.331 < 2e-16 ***
## imr_acquisition 1.910e+01 9.657e+00 1.978 0.0489 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.45 on 285 degrees of freedom
## Multiple R-squared: 0.9623, Adjusted R-squared: 0.9615
## F-statistic: 1212 on 6 and 285 DF, p-value: < 2.2e-16writeLines("Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with others")## Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with otherscar::vif(Ch03.linear) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity ## acq_expense acq_expense_sq industry revenue employees imr_acquisition
## 120.587121 105.260319 1.022881 1.077557 1.741442 2.969913В этой линейной модели также наблюдается несколько малозначимых предикторов и явно обозначенная мультиколлинеарность признаков. При все этом качество линейной модели авторов очень высокое. Однако нам необходимо убедиться в ее устойчивости, ведь такое большое качество может быть порождением переобучения модели на конкретных данных.
Linear (AR version)
Теперь попробуем исключить из линейной модели незначимыq фактор industry и провести построение устойчивой регресии посредством бустинга, что позволит избежать проявления переобучения при построении модели.
# Fit Linear Regression Model for Initial Order Quantity (Improved)
uno = 'y_hat'
set.seed(2018)
(Ch03ln.AR <- train(first_purchase ~ acq_expense + acq_expense_sq + revenue + employees + imr_acquisition,
data = filter(customerAcquisition, acquisition == 1), method = "lm"))#,## Linear Regression
##
## 292 samples
## 5 predictor
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 292, 292, 292, 292, 292, 292, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 30.90976 0.9609925 24.50053
##
## Tuning parameter 'intercept' was held constant at a value of TRUE # trControl = trainControl(method = "none", number = 1)))
summary(Ch03ln.AR)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -80.752 -19.620 -1.614 19.205 97.491
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.789e+01 4.834e+01 -1.404 0.1613
## acq_expense 7.135e-01 1.297e-01 5.502 8.33e-08 ***
## acq_expense_sq -7.729e-04 9.185e-05 -8.415 1.91e-15 ***
## revenue 3.035e+00 1.107e-01 27.432 < 2e-16 ***
## employees 2.543e-01 4.747e-03 53.578 < 2e-16 ***
## imr_acquisition 1.975e+01 9.610e+00 2.055 0.0408 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.42 on 286 degrees of freedom
## Multiple R-squared: 0.9622, Adjusted R-squared: 0.9616
## F-statistic: 1457 on 5 and 286 DF, p-value: < 2.2e-16# Coefficient Estimates and 95% CI
writeLines("\n Coefficient Estimates and 95% CI")##
## Coefficient Estimates and 95% CIcar::Confint(Ch03ln.AR$finalModel) %>%
exp() %>%
arm::pfround(digits = 3)## Estimate 2.5 % 97.5 %
## (Intercept) 0.000000e+00 0.000000e+00 6.864880e+11
## acq_expense 2.041000e+00 1.581000e+00 2.635000e+00
## acq_expense_sq 9.990000e-01 9.990000e-01 9.990000e-01
## revenue 2.081000e+01 1.673800e+01 2.587400e+01
## employees 1.290000e+00 1.278000e+00 1.302000e+00
## imr_acquisition 3.781997e+08 2.307000e+00 6.198704e+16# Linear regression diagnostics
writeLines("\n Chisq test of predictors")##
## Chisq test of predictorscar::Anova(Ch03ln.AR$finalModel, type="II", test="Chisq")## Anova Table (Type II tests)
##
## Response: .outcome
## Sum Sq Df F value Pr(>F)
## acq_expense 28022 1 30.2742 8.334e-08 ***
## acq_expense_sq 65539 1 70.8053 1.911e-15 ***
## revenue 696523 1 752.4965 < 2.2e-16 ***
## employees 2657110 1 2870.6392 < 2.2e-16 ***
## imr_acquisition 3910 1 4.2243 0.04076 *
## Residuals 264726 286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1writeLines("\n Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with others")##
## Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with otherscar::vif(Ch03ln.AR$finalModel) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity## acq_expense acq_expense_sq revenue employees imr_acquisition
## 119.622521 104.600905 1.063681 1.731720 2.945703writeLines("\n Variable Importance for Model \n")##
## Variable Importance for Modelcaret::varImp(Ch03ln.AR) %>% .$importance %>% print.AsIs()## Overall
## acq_expense 6.690003
## acq_expense_sq 12.342602
## revenue 49.252460
## employees 100.000000
## imr_acquisition 0.000000# Create the scatter plots Linear Model versus model predictors
# https://stats.idre.ucla.edu/r/seminars/ggplot2_intro/
predictors <- Ch03ln.AR$trainingData %>% colnames(.) %>% .[-1]
Ch03ln.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(Ch03ln.AR$finalModel, title = paste(paste(formula(Ch03ln.AR)[c(2, 1, 3)], collapse = " ")))## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
# wide variety of diagnostic plots for checking the quality of regression fit
car::influenceIndexPlot(Ch03ln.AR$finalModel)
# Goodness-of-fit of normal distributions to residuals of Linear Model (if Kolmogorov-Smirnov statistic < 0.05)
residuals(Ch03ln.AR) %>% fitdistrplus::fitdist("norm") %>% fitdistrplus::gofstat()## Goodness-of-fit statistics
## 1-mle-norm
## Kolmogorov-Smirnov statistic 0.03277147
## Cramer-von Mises statistic 0.04864034
## Anderson-Darling statistic 0.31671278
##
## Goodness-of-fit criteria
## 1-mle-norm
## Akaike's Information Criterion 2821.092
## Bayesian Information Criterion 2828.445Улучшенная модель без признака industry стала значимее, ни сколько не потеряла в своей точности, но точно приобрела устойчивость после бустинга. Следует признать, что остатки линейной модели и в самом деле распределяются по нормальному распределению, поскольку значение статистики Колмогорова-Смирнова p < 0.05.
Error of Linear Model
Последний шаг - сравнить предсказанные значения First_Purchase по предложенной авторами модели, а также моей улучшенной модели с фактическими значениями. При этом не забываем, что уравнение у авторов задано в два этапа:
\[ \displaystyle \large E(Initial \enspace Order \enspace Quantity) - P(Acquisition = 1) * E(First \enspace Purchase | Acquisition = 1) \hspace{.5 in} [5] \]
Авторы делают это, вычисляя среднее абсолютное отклонение (MAD или MAE) и среднюю абсолютную процентную погрешность (MAPE).
# Computing the Mean Absolute Deviation (MAD) and Mean Absolute Percent Error (MAPE)
# SAS Code: pred_fp = probnorm(xb_probit) * (xbeta);
with(filter(customerAcquisition, acquisition == 1), {
pred_fp <- pnorm(xbeta[ifelse(customerAcquisition$acquisition == 1, TRUE, FALSE)]) *
predict(Ch03.linear) #, newdata = filter(customerAcquisition, acquisition == 1))
# mad = mean(abs(first_purchase - pred_fp));
writeLines(sprintf("Mean Absolute Deviation (MAD): %.0f долл.", mean(abs(first_purchase - pred_fp))))
# mape = mean(abs(first_purchase - pred_fp)/first_purchase);
writeLines(sprintf("Mean Absolute Percent Error (MAPE): %.3f %% \n", mean(abs(first_purchase - pred_fp) / first_purchase) * 100))
pred_fp <- pnorm(xbeta[ifelse(customerAcquisition$acquisition == 1, TRUE, FALSE)]) *
predict(Ch03ln.AR, newdata = filter(customerAcquisition, acquisition == 1))
# mad = mean(abs(first_purchase - pred_fp));
writeLines(sprintf("Improved Mean Absolute Deviation (MAD): %.0f долл.", mean(abs(first_purchase - pred_fp))))
# mape = mean(abs(first_purchase - pred_fp)/first_purchase);
writeLines(sprintf("Improved Mean Absolute Percent Error (MAPE): %.3f %% \n", mean(abs(first_purchase - pred_fp) / first_purchase) * 100))
# mad1 = mean(abs(first_purchase - mean(first_purchase));
mad1 <- mean(abs(first_purchase - mean(first_purchase)))
writeLines(sprintf("Naive Mean Absolute Deviation (MAD1): %.0f долл.", mad1))
# mape1 = mean(abs(first_purchase - mean(first_purchase))/first_purchase);
mape1 <- mean(abs(first_purchase - mean(first_purchase)) / first_purchase) * 100
writeLines(sprintf("Naive Mean Absolute Percent Error (MAPE1): %.3f %%", mape1))
})## Mean Absolute Deviation (MAD): 52 долл.
## Mean Absolute Percent Error (MAPE): 18.686 %
##
## Improved Mean Absolute Deviation (MAD): 52 долл.
## Improved Mean Absolute Percent Error (MAPE): 18.564 %
##
## Naive Mean Absolute Deviation (MAD1): 127 долл.
## Naive Mean Absolute Percent Error (MAPE1): 135.482 %Если бы авторы вместо построения этой двухэтапной модели просто использовали наивный подход - применяли для прогноза среднее значение first_purchase (372.47 долл.)то, они бы обнаружили, что наивный MAD1 и наивный MAPE1. Таким образом, авторская модель сформировала лучший прогноз значений начального стоимости заказа, чем наивный усредненный подход.
Как это использовать?
Благодаря этой двухэтапной модели менеджмент может точнее прогнозировать размеры первоначального заказа вновь привлеченных клиентов, что существенно улучшает аккуратность планирования доходов от поступлений новых потребителей. При этом руководство может четко представлять драйверы, которые влияют на величину первого заказа вновь привлеченных клиентов.
Эмпирический пример: Продолжительность / время
Приобретение клиентов имеет ключевое значение для многих фирм. Тем не менее, многие фирмы также хотят смотреть за пределы точки приобретения и отвечать на такие вопросы, как «Как долго будет вновь приобретенный клиент по-прежнему будет оставаться клиентом?» Итак, в этом примере авторы попытаются выявить драйверы новой продолжительности работы с клиентами. Мы ожидаем, что продолжительность обслуживания клиентов будет зависеть от характеристик взаимоотношений с клиентами и информации о фирме. Но мы также хотим узнать, влияет ли попытка привлечения фирмы на продолжительность сотрудничества каждого нового клиента с фирмой. В конце этого примера авторы объяснят:
Определить драйверы продолжительности жизни нового клиента.
Предсказывать продолжительность жизни каждого нового клиента.
Информация, необходимая для этой модели, включает следующий список переменных:
| Переменные | Описание |
|---|---|
| Зависимая переменная | |
| Duration | Время в днях, когда компании была или продолжает быть клиентом, цензуированная до 730 дней |
| Предикторы | |
| Acq_Expense | Долларовые расходы на маркетинг по привлечению клиента |
| Acq_Expense_SQ | Квадрат расходов на привлечение |
| Ret_Expense | Долларовые расходы на маркетинг по удержание клиента |
| Ret_Expense_SQ | Квадрат расходов на удержание |
| Crossbuy | Количество категорий товаров / услуг, приобретенных клиентом |
| Frequency | Количество покупок клиента во время окна наблюдения |
| Frequency_SQ | Квадрат количества покупок |
| Industry | 1, если клиент в секторе B2B, 0 в противном случае, т. е. B2C |
| Revenue | Годовой доход от продаж компании (млн долл.) |
| Employees | Количество сотрудников в компании |
| Censor | 1, если клиент оставался в конце окна наблюдения, 0 в противном случае |
В этом случае данные выборки взяты из фирмы, работающей на корпоративном рынке, которая фактически наблюдает, когда клиент уходит в отток до истечения срока договора. Из списка переменных мы видим, что мы имеем одну зависимую переменную (Duration), которая подвергается цензуированию справо. Окно наблюдения за клиентами составляет всего два года (730 дней) после того, как были заклюбчен договор. Таким образом, мы наблюдаем, как клиент покидает компанию, если это произойдет до конца второго года. В результате данные для Duration подвергаются правосторонней цензуированию на 730 день - это означает, что любой клиент, значение Duration которого равен 730 в таблице данных, еще не покинул фирму в конце окна наблюдения. У авторов также есть 10 независимых переменных, которые мы надеемся, что это объяснят изменения в продолжительности сотрудничества каждого нового клиента с фирмой.
Множительные оценки Каплана-Мейера
Чтобы дать оценку функции выживания графическим методом и определить разницу между различными претикторами и их ковариатами полезно построить кривые Каплана—Мейера (англ. “Kaplan–Meier Estimator”). Для построения подобных кривых необходимо обратиться к важнейшему пакету в R по тематике анализа выживания - survival, а также специализированным графикам из пакета survminer.
KM Curve by Duration
# survminer: Survival Analysis and Visualization
# http://www.sthda.com/english/wiki/survival-analysis-basics & https://habr.com/post/348754/
library('survival')##
## Attaching package: 'survival'## The following object is masked from 'package:caret':
##
## clusterlibrary('survminer')## Loading required package: ggpubr## Loading required package: magrittr##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extractfit <- survival::survfit(Surv(time = duration, event = (censor == 0), type = "right") ~ 1,
data = filter(customerAcquisition, acquisition == 1))
# Graph of Kaplan–Meier Estimator
ggsurvplot(fit, data = filter(customerAcquisition, acquisition == 1), # survfit object with calculated statistics.
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for
# point estimates of survival curves.
fontsize = 2,
xlim = c(0, max(customerAcquisition$duration)),
# present narrower X axis, but not affect survival estimates.
xlab = "Time in days", # customize X axis label.
break.time.by = 30, # break X axis in time intervals by 30 days.
ggtheme = theme_light(), # customize plot and risk table with a theme.
risk.table = "abs_pct", # to show both absolute number and percentage
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.height = 0.25, # the height of the risk table
risk.table.y.text = FALSE,# show bars instead of names in text annotations
# in legend of risk table.
surv.median.line = "hv"#, # add the median survival pointer.
# legend.labs =
# c("B2C", "B2B") # change legend labels.
) 
Очевидно, что половину клиентов завершили сотрудничество к 615 суткам.
KM Curve by Duration vs industry
Важно рассмотреть как ведут себя разные сегменты клиентов (страты) в ходе сотрудничества с фирмой. Например, страты B2C и B2B.
fit <- survival::survfit(Surv(time = duration, event = (censor == 0), type = "right") ~ industry, data = filter(customerAcquisition, acquisition == 1))
# Graph of Kaplan–Meier Estimator
ggsurvplot(fit, data = filter(customerAcquisition, acquisition == 1), # survfit object with calculated statistics.
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for
# point estimates of survival curves.
palette = c("coral", "#E7B800", "#2E9FDF"),
xlim = c(0, max(customerAcquisition$duration)),
xlab = "Time in days", # customize X axis label.
break.time.by = 30, # break X axis in time intervals by 30 days.
ggtheme = theme_light(), # customize plot and risk table with a theme.
fontsize = 2.5,
add.all = TRUE, # add the survival curve of pooled patients
risk.table = TRUE, # show risk table.
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.height = 0.25, # the height of the risk table
risk.table.y.text = FALSE,# show bars instead of names in text annotations
# in legend of risk table.
surv.median.line = "hv", # add the median survival pointer.
cumevents = TRUE,
legend.labs =
c("All", "B2C", "B2B") # change legend labels.
) 
На графике видно, что время сотрудничества у компаний “B2C” значимо ниже, чем у продолжительно сотрудничающих “B2B”.
KM Curve by Duration vs crossbuy
# survminer: Survival Analysis and Visualization
# http://www.sthda.com/english/wiki/survival-analysis-basics & https://habr.com/post/348754/
fit <- survival::survfit(Surv(time = duration, event = (censor == 0), type = "right") ~ crossbuy, data = filter(customerAcquisition, acquisition == 1))
# Graph of Kaplan–Meier Estimator
ggsurvplot(fit, data = filter(customerAcquisition, acquisition == 1), # survfit object with calculated statistics.
pval = TRUE, # show p-value of log-rank test.
conf.int = FALSE, # show confidence intervals for
# point estimates of survival curves.
xlim = c(0, max(customerAcquisition$duration)),
xlab = "Time in days", # customize X axis label.
break.time.by = 30, # break X axis in time intervals by 30 days.
ggtheme = theme_light(), # customize plot and risk table with a theme.
fontsize = 2.5,
add.all = TRUE, # add the survival curve of pooled patients
surv.median.line = "hv" # add the median survival pointer.
) 
# # Tests if there is a difference between two or more survival curves using the G-rho family of tests
# x <- survival::survdiff(Surv(time = duration, event = (censor == 0), type = "right") ~ crossbuy, data = filter(customerAcquisition, acquisition == 1))Изменение продолжительности в зависимости от количества категорий товаров / услуг, приобретенных клиентом, crossbuy уже не столь велико - p = 0.8.
Построение и верификация модели времени ускоренного отказа
Чтобы определить драйверы продолжительности работы клиента, мы используем модель Времени Ускореннного Отказа (англ. “Accelerated Failure Time (AFT) Model”). В этом случае имеется следующее уравнение:
\[ \displaystyle \large \ln(Duration) = X'\beta + \sigma \varepsilon. \hspace{.5 in} [6]\] где \(\ln(Duration)\) - натуральный логарифм времени, в течение которого клиент сотрудничал с фирмой;
\(X'\) - матрица из 10 независимых переменных;
\(\beta\) - вектор коэффициентов;
\(\sigma\) - параметр масштаба;
\(\varepsilon\) - ошибки модели.
Авторы оценивают модель, полагая, что \(\varepsilon\) следует распределению Вейбулла (англ. “Weibull distribution”) - одному из распространенных, наряду с экспоненциальным и распределение Гомпертца, среди моделей времени ускоренного отказа.
Функции плотности вероятности распределения Вейбулла -
\[ \displaystyle \large f(x; \alpha, \gamma) = \begin{cases} \frac{\gamma}{\alpha}\left(\frac{x}{\alpha}\right)^{\gamma - 1}e^{-(x/\alpha)^{\gamma}} & x\geq0 ,\\ 0 & x<0 \end{cases} \hspace{.5 in} [7] \]
Авторы оценивают модель только с 292 перспективными потребителями, которые были привлечены в качестве клиентов (Acquisition == 1).
Duration Models
# 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
UpLimit <- max(customerAcquisition$duration)
library('VGAM') # Vector Generalized Linear and Additive Models## Loading required package: stats4## Loading required package: splines##
## Attaching package: 'VGAM'## The following object is masked from 'package:car':
##
## logit## The following object is masked from 'package:caret':
##
## predictors## The following object is masked from 'package:tidyr':
##
## fill# Left- & Right-Censored Tobit (Normally distributed error term) AFT Model from 'VGAM' package
(Ch03.tb <- VGAM::vglm(duration ~ acq_expense + acq_expense_sq + ret_expense + ret_expense_sq + crossbuy +
frequency + frequency_sq + industry + revenue + employees,
family = tobit(Lower = 0, Upper = UpLimit, type.fitted = "censored"),
data = filter(customerAcquisition, acquisition == 1))) %>%
summary##
## Call:
## VGAM::vglm(formula = duration ~ acq_expense + acq_expense_sq +
## ret_expense + ret_expense_sq + crossbuy + frequency + frequency_sq +
## industry + revenue + employees, family = tobit(Lower = 0,
## Upper = UpLimit, type.fitted = "censored"), data = filter(customerAcquisition,
## acquisition == 1))
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## mu -1.9756 -0.5832 0.02388 0.486433 2394
## loge(sd) -0.9186 -0.4548 -0.14476 0.005865 2319
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 4.601e+02 2.070e+02 2.222 0.026252 *
## (Intercept):2 4.736e+00 5.331e-02 88.830 < 2e-16 ***
## acq_expense -5.288e-01 6.258e-01 -0.845 0.398147
## acq_expense_sq -1.015e-03 4.481e-04 -2.266 0.023436 *
## ret_expense 2.899e-01 6.479e-02 4.474 7.68e-06 ***
## ret_expense_sq -1.290e-05 2.128e-05 -0.606 0.544366
## crossbuy 1.907e+01 4.942e+00 3.858 0.000114 ***
## frequency 3.539e+01 7.521e+00 4.706 2.53e-06 ***
## frequency_sq -4.076e-01 4.262e-01 -0.956 0.338862
## industry 1.633e+02 1.741e+01 9.384 < 2e-16 ***
## revenue 2.515e+00 5.125e-01 4.908 9.19e-07 ***
## employees 2.099e-02 1.780e-02 1.179 0.238380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: mu, loge(sd)
##
## Log-likelihood: -1008.813 on 572 degrees of freedom
##
## Number of iterations: 8
##
## No Hauck-Donner effect found in any of the estimates# # Censored Normal Distribution of error term Or Censored Tobit
# Extra <- with(filter(customerAcquisition, acquisition == 1),
# list(leftcensored = 0, rightcensored = (duration == UpLimit)))
#
# (Ch03.nr <- VGAM::vglm(duration ~ acq_expense + acq_expense_sq + ret_expense + ret_expense_sq + crossbuy +
# frequency + frequency_sq + industry + revenue + employees,
# family = cens.normal(), extra = Extra,
# data = filter(customerAcquisition, acquisition == 1))) %>%
# summary
# Right-Censored Rayleigh Distribution (special case of Weibull distribution with the "shape" parameter gamma = 2)
Extra = with(filter(customerAcquisition, acquisition == 1), list(rightcensored = (duration == UpLimit) ))
(Ch03.ry <- VGAM::vglm(duration ~ acq_expense + acq_expense_sq + ret_expense + ret_expense_sq + crossbuy +
frequency + frequency_sq + industry + revenue + employees,
family = cens.rayleigh, extra = Extra,
data = filter(customerAcquisition, acquisition == 1))) %>%
summary##
## Call:
## VGAM::vglm(formula = duration ~ acq_expense + acq_expense_sq +
## ret_expense + ret_expense_sq + crossbuy + frequency + frequency_sq +
## industry + revenue + employees, family = cens.rayleigh, data = filter(customerAcquisition,
## acquisition == 1), extra = Extra)
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## loge(scale) -0.7192 -0.2283 0.1041 0.3123 1.06
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.601e+00 9.164e-01 3.929 8.52e-05 ***
## acq_expense 3.547e-03 2.776e-03 1.278 0.201
## acq_expense_sq -8.364e-06 1.948e-06 -4.293 1.76e-05 ***
## ret_expense 1.152e-03 2.876e-04 4.007 6.15e-05 ***
## ret_expense_sq -2.486e-08 9.665e-08 -0.257 0.797
## crossbuy 1.016e-01 2.266e-02 4.486 7.26e-06 ***
## frequency 1.497e-01 3.539e-02 4.231 2.33e-05 ***
## frequency_sq -1.959e-03 2.023e-03 -0.968 0.333
## industry 6.573e-01 8.029e-02 8.186 2.70e-16 ***
## revenue 1.298e-02 2.387e-03 5.440 5.34e-08 ***
## employees 1.204e-04 8.096e-05 1.487 0.137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 1
##
## Name of linear predictor: loge(scale)
##
## Log-likelihood: -950.746 on 281 degrees of freedom
##
## Number of iterations: 9
##
## No Hauck-Donner effect found in any of the estimates# Weibull distribution do not to handle other censoring situations yet (except gamma = 2)
# (Ch03.ln <- VGAM::vglm(duration ~ acq_expense + acq_expense_sq + ret_expense + ret_expense_sq + crossbuy + frequency + frequency_sq + industry + revenue + employees, family = weibullR(), data = filter(customerAcquisition, acquisition == 1))) %>% summaryК сожалению, в моем распоряжении оказалась только один пакет, который способен построить на этих данных специальный случай распределения Вейбулла с \({\gamma = 2, \alpha = {\sqrt 2} \sigma}\), т. е. распределением Рэлея (англ. “Rayleigh distribution”). Ошибка с более привычным пакетом survival происходит из-за включения в модель предикторов acq_expense_sq и ret_expense_sq, обладающих чрезвычайно большими значениями.
В принципе, Экспоненциальное распределение тоже специальный случай распределения Вейбулла с \({\gamma = 1}\). Поскольку разработчики пакета VGAM с 2007 г. так и не смогли подготовить обещанную ими полноценную цензуированную версию регрессионную модель распределения Вейбулла, то нам придется ограничиться частным случаем распределения Вейбулла.
# https://stats.stackexchange.com/questions/19866/how-to-fit-a-weibull-distribution-to-input-data-containing-zeroes
foo <- function(theta, x)
{
if (theta <= -min(x)) return(Inf);
f <- MASS::fitdistr(x+theta, 'weibull')
-2*f$loglik
}
x <- (filter(customerAcquisition, acquisition == 1)$duration -
predict(Ch03.tb, newdata = filter(customerAcquisition, acquisition == 1), type = "response")) %>%
as.vector
set.seed(2018)
# Cullen & Fley Graph: Empirical Weibull Distribution
fitdistrplus::descdist(x, boot = 500, discrete = FALSE)## summary statistics
## ------
## min: -269.0988 max: 133.9924
## median: 0
## mean: -22.0155
## estimated sd: 65.52328
## estimated skewness: -0.6274331
## estimated kurtosis: 3.643327mtext("Weibull distribution of `filter(customerAcquisition, acquisition == 1)$duration - predict(Ch03.tb)`",
line = 0.5, cex =0.75)
bar <- optimize(foo, lower=-min(x)+0.001, upper=-min(x)+10, x=x)
library('fitdistrplus')## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select(fitW <- fitdistrplus::fitdist(x+bar$minimum, 'weibull') )## Fitting of the distribution ' weibull ' by maximum likelihood
## Parameters:
## estimate Std. Error
## shape 4.584192 0.2128151
## scale 280.836008 3.7492335plot(fitW)
fitdistrplus::gofstat(fitW)## Goodness-of-fit statistics
## 1-mle-weibull
## Kolmogorov-Smirnov statistic 0.2415202
## Cramer-von Mises statistic 2.4557373
## Anderson-Darling statistic 11.0484376
##
## Goodness-of-fit criteria
## 1-mle-weibull
## Akaike's Information Criterion 3271.692
## Bayesian Information Criterion 3279.046# https://stats.stackexchange.com/questions/19866/how-to-fit-a-weibull-distribution-to-input-data-containing-zeroes
# # R-code for the Blischke-Scheuer method
# fit_weibull <- function(x)
# {
# xbar <- mean(x)
# varx <- var(x)
# f <- function(b){return(gamma(1+2/b)/gamma(1+1/b)^2 - 1 - varx/xbar^2)}
# bhat <- uniroot(f,c(0.02,50))$root
# ahat <- xbar/gamma(1+1/bhat)
# return(c(ahat,bhat))
# }
#
# (z <- fit_weibull(x + bar$minimum))Наряду с модификацией авторской модели была построена классическая Тобит-модель (англ. “Tobit model”), которая представляет собой полноценно цензуированную модель с нормальным распределением. Кроме того, проверялась и “наивная” модель по средней продолжительности нецензуриванных фирм, сотрудничество в которыми длилось менее двух лет, которая равна 327 дней. Заметим, что на странице 46 авторы описывают эту величину как 333 дней, а в листинге программы SAS на странице 56 уже как 327.
# Accuracy measures for Censored Tobit
y <- filter(customerAcquisition, acquisition == 1)$duration
y_hat <- predict(Ch03.tb, newdata = filter(customerAcquisition, acquisition == 1), type = "response") %>% as.vector # if_else(fitted(Ch03.tb) > UpLimit, UpLimit, if_else(fitted(Ch03.tb) < 0, 0, fitted(Ch03.tb)))
m <- forecast::accuracy(rep(mean(filter(customerAcquisition, acquisition == 1 & censor == 0)$duration), length(y)), y); attributes(m)$ dimnames[[1]] <- " Naïve Model (Mean) Accuracy: "; m## ME RMSE MAE MPE MAPE
## Naive Model (Mean) Accuracy: 186.1359 309.3805 275.1028 -49.74756 116.4551m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- " Censored Tobit's Accuracy: "; m## ME RMSE MAE MPE MAPE
## Censored Tobit's Accuracy: -22.0155 69.01651 44.92601 -6.518997 19.9415# # Accuracy measures for Censored Normal Distribution AFT Model
# y_hat <- if_else(fitted(Ch03.nr) > UpLimit, UpLimit, if_else(fitted(Ch03.nr) < 0, 0, fitted(Ch03.nr)))
#
# m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "Censored Normal Model's Accuracy: "; m
# Accuracy measures for Right-Censoring Rayleigh AFT Model
y_hat <- if_else(fitted(Ch03.ry) > UpLimit, UpLimit, if_else(fitted(Ch03.ry) < 0, 0, fitted(Ch03.ry)))
m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "Right-censored Model's Accuracy: "; m## ME RMSE MAE MPE MAPE
## Right-censored Model's Accuracy: -7.790355 50.34878 25.98158 -0.1785197 8.393484Средняя абсолютная процентная ошибка (‘MAPE’) по наивной модели составила 116.5%,
что существенно выше, чем по тобит-модели (19.9%)
и особенно по модели с остатками по цензуированному Вейбуллу (8.4%).
Weibull (AR version)
Поскольку ряд предикторов в полученной модели с остатками по цензуированному Вейбуллу незначимы, то можно построить более надежную модели выживания (продолжительности сотрудничества с клиентами).
# Right-Censored Rayleigh Distribution (special case of Weibull distribution with the "shape" parameter gamma = 2)
Extra = with(filter(customerAcquisition, acquisition == 1), list(rightcensored = (duration == UpLimit) ))
(Ch03.RY <- VGAM::vglm(duration ~ acq_expense_sq + ret_expense + crossbuy + frequency + industry + revenue,
family = cens.rayleigh, extra = Extra,
data = filter(customerAcquisition, acquisition == 1))) %>%
summary##
## Call:
## VGAM::vglm(formula = duration ~ acq_expense_sq + ret_expense +
## crossbuy + frequency + industry + revenue, family = cens.rayleigh,
## data = filter(customerAcquisition, acquisition == 1), extra = Extra)
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## loge(scale) -0.7943 -0.2269 0.1114 0.3159 1.386
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.092e+00 1.985e-01 25.657 < 2e-16 ***
## acq_expense_sq -5.905e-06 2.229e-07 -26.489 < 2e-16 ***
## ret_expense 1.075e-03 7.183e-05 14.971 < 2e-16 ***
## crossbuy 9.987e-02 2.233e-02 4.473 7.73e-06 ***
## frequency 1.159e-01 8.678e-03 13.355 < 2e-16 ***
## industry 6.463e-01 7.902e-02 8.178 2.88e-16 ***
## revenue 1.309e-02 2.335e-03 5.603 2.10e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 1
##
## Name of linear predictor: loge(scale)
##
## Log-likelihood: -953.2191 on 285 degrees of freedom
##
## Number of iterations: 8
##
## No Hauck-Donner effect found in any of the estimates# Accuracy measures for Right-Censoring Rayleigh AFT Model by AR
y_hat <- y_hat <- if_else(fitted(Ch03.RY) > UpLimit, UpLimit, if_else(fitted(Ch03.RY) < 0, 0, fitted(Ch03.RY)))
m <- forecast::accuracy(y_hat, y); attributes(m)$ dimnames[[1]] <- "Right-censored Model's Accuracy: "; m## ME RMSE MAE MPE MAPE
## Right-censored Model's Accuracy: -8.402617 56.74402 28.64044 -0.6893516 9.271426И хотя у этой модели немного ниже точность прогнозирования, тем не менее она более надежна, так как вместо ранее примененных 10 предикторов оставлено только 6 значимых, а ошибка по MAPE увеличилась на менее, чем один процент.
Kaplan–Meier Graphs
library('survival') # Survival Analysis routines (Surv objects, Kaplan-Meier curves, Cox and parametric AFT models)
# library('AER') # Applied Econometrics with R for the book Christian Kleiber and Achim Zeileis (2008)
# # Censored Tobit Model from 'AER' package weibull
# (Ch03.ln <- AER::tobit(duration ~ ret_expense + crossbuy + frequency + industry + revenue,
# right = UpLimit, dist = 'weibull', data = filter(customerAcquisition, acquisition == 1))) %>%
# summary
y_hat <- if_else(fitted(Ch03.ry) > UpLimit, UpLimit, fitted(Ch03.ry))
# Kaplan–Meier Estimator of Survival Curves
Ch03.km <- survfit(Surv(time = duration, event = (censor == 0), type = "right") ~ 1,
data = filter(customerAcquisition, acquisition == 1))
plot(Ch03.km, xlab="Time in days", ylab="Функция выживания")
lines(survfit(Surv(time = fitted(Ch03.tb), event = (fitted(Ch03.tb) != UpLimit), type = "right") ~ 1,
data = filter(customerAcquisition, acquisition == 1)), conf.int = FALSE, col = "green")
lines(survfit(Surv(time = y_hat, event = (y_hat != UpLimit), type = "right") ~ 1,
data = filter(customerAcquisition, acquisition == 1)), conf.int = FALSE, col = "red")
## Add legends
legend(x = "bottomleft",
legend = c("Kaplan-Meier Estimator",
sprintf("Censored Tobit (MAPE= %.1f %%)", forecast::accuracy(fitted(Ch03.tb)[, 1], y)[, "MAPE"]),
sprintf("Right-Censored Weibull [gamma=2] (MAPE= %.1f %%)", forecast::accuracy(y_hat, y)[, "MAPE"])),
lwd = 2, bty = "n",
col = c("black", "green", "red"))
Survival Probability Curves
Сравнение на графике выживания показывает, то кривая, построенная моделью, проходит вблизи фактической кусочно-линейной функции для аппроксимации функции выживания Каплана-Мейера мало отклоняясь от этой эмпирической кривой.
В некоторых случаях исследователи могут захотеть моделировать одновременно несколько показателей привлечения клиентов, например, количество вновь приобретенных клиентов и начальное количество заказа, а системы линейной регрессии могут выполнить такую задачу. Кроме того, возможно, что метрики привлечения, такие как количество клиентов, приобретенных посредством промоутеров или взаимодействие в интернете и прибыльность компании, коррелированы и разнонаправлены во времени, а ’VAR -а это подходящий метод моделирования для учета корреляционных и временных эффектов. Кроме того, исследователи могут захотеть смоделировать продолжительность сотрудничества со вновь приобретенными клиентами и определить последствия, которые могут иметь различные объясняющие переменные во до момента прекращения сотрудничества. Поскольку большинство данных о продолжительности содержат цензуированые наблюдения, регрессия методом наименьших квадратов обеспечивала бы предвзятые оценки. Классическим параметрическим методом для данных о выживании является модель времени ускоренного отказа (AFT), которая может быть оценена методом максимального правдоподобия.
В моделировании длительности сотрудничества часто используются вариации модели времени ускоренного отказа, в зависимости от предположения о её теоретическом законе распределении остатков, таких как Вейбулла (в частности, экспоненциальное, Рэлея), Гумбела, гамма, лог-логистическое и лог-нормальное распределения.
Как это использовать?
Авторы полагают, что вывод модели весьма полезен для определения того, какого клиента как именно можно удержать фирма. Кроме того, результаты модели времени ускоренного отказа могут служить драйверами удержаний клиентов, какие факторы воздействуют на увеличение или уменьшение шансов на продолжение сотрудничества, которые полезны для менеджеров при принятии решений в будущих кампаниях по сохранению клиентов.
Теперь можно узнать, как изменения в расходах на удержание, дополнительных услуг и частоты обращений,а также изменения характеристик имеющихся клиентов могут либо увеличить, либо уменьшить вероятность завершения сотрудничества с потребителем. Эта информация может дать значительную информацию руководителям, которым поручено определить оптимальный объем ресурсов, которые необходимо потратить на удержание клиентов.
Эмпирический пример: доходность фирмы
Последним этапом анализа привлечения потребителей является определение того, выгодны ли покупатели, которые были приобретены. Авторы также хотели знать, можем ли определить драйверы рентабельности клиента, чтобы увидеть, могут ли будущие усилия по приобретению помочь получить большее количество прибыльных клиентов. Поэтому в этом примере можно научиться:
Определять драйверы доходности клиента.
Рассчитывать предиктивную точность модели доходности клиентов.
Информация, которая нам понадобится для этой модели, включает следующий список переменных для 292 потенциальных клиентов, которые заключили договора с фирмой (Acquisition == 1):
| Переменные | Описание |
|---|---|
| Зависимая переменная | |
| Censor | 1, если клиент оставался в конце окна наблюдения, 0 в противном случае |
| CLV | Прогноз пожизненной ценности клиента. Это 0, если клиента не было или уже ушел в отток (’000) |
| Предикторы | |
| Acq_Expense | Долларовые расходы на маркетинг по привлечению клиента |
| Acq_Expense_SQ | Квадрат расходов на привлечение |
| imr_censor или Lambda (\(\lambda\)) | Рассчитанное обратное отношение Миллса из модели привлечения клиента |
| Ret_Expense | Долларовые расходы на маркетинг по удержание клиента |
| Ret_Expense_SQ | Квадрат расходов на удержание |
| First_Purchase | Долларовая стоимость первой покупки (0, если клиент не был приобретен) |
| Crossbuy | Количество категорий товаров / услуг, приобретенных клиентом |
| Frequency | Количество покупок клиента во время окна наблюдения |
| Frequency_SQ | Квадрат количества покупок |
| Industry | 1, если клиент в секторе B2B, 0 в противном случае, т. е. B2C |
| Revenue | Годовой доход от продаж компании (млн долл.) |
| Employees | Количество сотрудников в компании |
В этом случае мы имеем две зависимые переменные. Одна зависимая переменная, называемая Censor, помогает нам отделять тех клиентов, которые все еще сотрудничают, и теми клиентами, которые уже покинули фирму. Вторая зависимая переменная называется CLV (англ. “Customer Lifetime Value”) или Пожизненная финансовая ценность клиента. Для этого примера для вас был предсказан показатель CLV (в тысячах) для каждого из приобретенных клиентов. Поскольку CLV является перспективной мерой (т. Е. Измеряет ожидаемую будущую прибыльность клиента), единственными клиентами с оценкой CLV являются те, кто по-прежнему остается клиентами фирмы. Это клиенты, которые отвечают следующим двум критериям: они имеют значение 1 для приобретения и значение 1 для цензора.
Как и в случае с первоначальным объемом покупки, мы должны учитывать как вероятность того, что клиент по-прежнему является клиентом (Censor), так и ожидаемой будущей доходностью клиента (CLV). Для этого получим следующее уравнение:
\[ \displaystyle \large E(Future \enspace Profitability) = P(Censor = 1) * E(CLV \enspace | \enspace Censor = 1) \hspace{.5 in} [8] \] В этом случае мы также признаем, что нам придется учитывать только часть компанию, котоыре удасться привлечь из выборки и затем оценить модель в двухэтапной структуре моделирования. Опять же, подобно процессу моделирования величины первоначального заказа, мы сначала моделируем Censor как пробит, используя все независимые переменные. Затем мы вычисляем коэффициенты обратного отношения Миллса для тех клиентов, у которых есть Censor == 1. Наконец, мы запускаем построение линейной регрессию с CLV методом наименьших квадратов в качестве зависимой переменной только для 135 компаний, которые все еще являются клиентами. Мы используем все перечисленные независимые переменные и отношение обратных Миллсов (\(\lambda\)), которые объявлены зависимой переменной imr_censor.Таким образом, авторы получают следующие результаты:
Первый этап расчета CLV - вероятность дожития клиента
Probit
# Fit Probit Model for Censor by Authors
Ch03.censor <- glm(factor(censor) ~ acq_expense + acq_expense_sq + ret_expense + ret_expense_sq + first_purchase +
crossbuy + frequency + frequency_sq + industry + revenue + employees,
data = filter(customerAcquisition, acquisition == 1), family = binomial(link = 'probit'))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(Ch03.censor)##
## Call:
## glm(formula = factor(censor) ~ acq_expense + acq_expense_sq +
## ret_expense + ret_expense_sq + first_purchase + crossbuy +
## frequency + frequency_sq + industry + revenue + employees,
## family = binomial(link = "probit"), data = filter(customerAcquisition,
## acquisition == 1))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.05384 -0.00001 0.00000 0.00000 2.40735
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.215e+01 1.355e+01 -2.374 0.017619 *
## acq_expense 6.955e-02 4.205e-02 1.654 0.098093 .
## acq_expense_sq -9.213e-05 3.646e-05 -2.527 0.011516 *
## ret_expense 5.975e-03 2.794e-03 2.138 0.032483 *
## ret_expense_sq 4.557e-07 9.229e-07 0.494 0.621485
## first_purchase -1.706e-03 7.843e-03 -0.217 0.827831
## crossbuy 7.579e-01 2.159e-01 3.510 0.000447 ***
## frequency 8.515e-01 2.529e-01 3.368 0.000758 ***
## frequency_sq -3.533e-03 1.147e-02 -0.308 0.758145
## industry 4.588e+00 1.049e+00 4.372 1.23e-05 ***
## revenue 7.185e-02 2.983e-02 2.408 0.016022 *
## employees 1.832e-03 2.077e-03 0.882 0.377885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 403.139 on 291 degrees of freedom
## Residual deviance: 42.965 on 280 degrees of freedom
## AIC: 66.965
##
## Number of Fisher Scoring iterations: 12writeLines(sprintf("-2 Log L of Intercept and Only Covariates: %.3f", -2 * logLik(Ch03.censor)[1]))## -2 Log L of Intercept and Only Covariates: 42.965writeLines(sprintf(" AIC (smaller is better): %.3f", extractAIC(Ch03.censor)[2]))## AIC (smaller is better): 66.965writeLines("\n Wald test of predictors")##
## Wald test of predictorscar::Anova(Ch03.censor, type="II", test="Wald")## Analysis of Deviance Table (Type II tests)
##
## Response: factor(censor)
## Df Chisq Pr(>Chisq)
## acq_expense 1 2.7363 0.0980933 .
## acq_expense_sq 1 6.3839 0.0115160 *
## ret_expense 1 4.5728 0.0324826 *
## ret_expense_sq 1 0.2438 0.6214847
## first_purchase 1 0.0473 0.8278310
## crossbuy 1 12.3228 0.0004475 ***
## frequency 1 11.3408 0.0007582 ***
## frequency_sq 1 0.0948 0.7581454
## industry 1 19.1153 1.231e-05 ***
## revenue 1 5.8005 0.0160220 *
## employees 1 0.7776 0.3778848
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1writeLines("\n Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with others")##
## Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with otherscar::vif(Ch03.censor) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity ## acq_expense acq_expense_sq ret_expense ret_expense_sq first_purchase crossbuy frequency frequency_sq
## 297.567692 387.043589 38.565400 40.789388 22.546020 2.916964 37.415927 24.858065
## industry revenue employees
## 5.366119 5.952052 18.949890writeLines("\n Censor Probit Model: Association of Predicted Probabilities and Observed Responses \n")##
## Censor Probit Model: Association of Predicted Probabilities and Observed Responsesprob <- predict(Ch03.censor, newdata = filter(customerAcquisition, acquisition == 1), type = "response")
caret::confusionMatrix(data = ifelse(prob > 0.5, "1", "0") %>% factor,
reference = filter(customerAcquisition, acquisition == 1)$censor %>% factor,
positive = "1", mode = "everything")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 153 3
## 1 4 132
##
## Accuracy : 0.976
## 95% CI : (0.9512, 0.9903)
## No Information Rate : 0.5377
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9518
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9778
## Specificity : 0.9745
## Pos Pred Value : 0.9706
## Neg Pred Value : 0.9808
## Precision : 0.9706
## Recall : 0.9778
## F1 : 0.9742
## Prevalence : 0.4623
## Detection Rate : 0.4521
## Detection Prevalence : 0.4658
## Balanced Accuracy : 0.9762
##
## 'Positive' Class : 1
## Как и в случае с пробит-моделью из начала главы, эту пробит-модель не вполне можно назвать удачной, ведь она содержит ряд незначимых предикторов, а значимые являются мультиколлинеарными между собой. Более того, при построении бинарного фактора Censor получено предупреждение “glm.fit: возникли подогнанные вероятности 0 или 1”, поскольку наблюдается зависимая переменная, которая идеально разделяет на нули и единицы, т. е. обучение ведется на «идеальном или квази совершенном разбиении».
Elastic-net
Решением этой и других проблем является использование пенализированной регрессии (англ. “Penalized Regression”), которая обеспечивает регуляризация, позволяющая находить приближённое решение некорректно поставленных задач. Загрузите пакет glmnet в R, способного устранить эту проблему и мультиколлиниарность нескольких предикторов посредством регуляризацией по Lasso и Ridge .
# Fit Lasso and Ridge (Elastic-net) regularized Generalized Linear Model for Censor (Improved)
uno = 'logit'
set.seed(2018)
Ch03cs.AR <- train(factor(censor) ~ acq_expense_sq + ret_expense + crossbuy + frequency + industry + revenue,
data = filter(customerAcquisition, acquisition == 1), metric = 'Kappa',
tuneGrid = expand.grid(.alpha = seq(0, 1, length=11),
## The elasticnet mixing parameter: alpha=1 is the lasso penalty, and alpha=0 the ridge penalty
.lambda = seq(0.001, 2, length=11)),
## LASSO (Least Absolute Shrinkage and Selection Operator)
method = "glmnet", family = "binomial", type.logistic="modified.Newton")#,
# trControl = trainControl(method = "none", number = 1))
predict(Ch03cs.AR$finalModel, newx = NULL, type = "coefficients", s = Ch03cs.AR$finalModel$tuneValue$lambda)## 7 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) -9.107204e+00
## acq_expense_sq -3.655406e-05
## ret_expense 6.847033e-03
## crossbuy 6.906295e-01
## frequency 7.378043e-01
## industry 4.143457e+00
## revenue 6.644739e-02writeLines("\n Variable Importance for Model \n")##
## Variable Importance for Modelcaret::varImp(Ch03cs.AR) %>% .$importance %>% print.AsIs()## Overall
## acq_expense_sq 0.0000000
## ret_expense 0.1643685
## crossbuy 16.6672196
## frequency 17.8057687
## industry 100.0000000
## revenue 1.6028024# Create the scatter plots Logit versus model predictors
prob <- predict(Ch03cs.AR, newdata = filter(customerAcquisition, acquisition == 1), type = "prob") %>% .[, "1"]
predictors <- Ch03cs.AR$coefnames # Ch03cs.AR$finalModel$xNames
Ch03cs.AR$trainingData %>%
dplyr::select(one_of(predictors)) %>%
mutate(logit = log(prob / (1 - prob))) %>%
gather(key = "predictors", value = "predictor.value", -logit) %>%
ggplot(aes(predictor.value, logit))+
geom_point(size = 0.5, alpha = 0.5) +
geom_smooth(method = "loess") +
facet_wrap(~ predictors, scales = "free_x") +
ylab(stringr::str_to_title(uno))## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : pseudoinverse used at -0.005## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : neighborhood radius 1.005## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : reciprocal condition number
## 7.9984e-031## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : There are other near singularities as
## well. 1.01## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at -0.005## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 1.005## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition number 7.9984e-031## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near singularities as well. 1.01
writeLines("\n Elastic-net Model: Association of Predicted Probabilities and Observed Responses \n")##
## Elastic-net Model: Association of Predicted Probabilities and Observed Responsescaret::confusionMatrix(data = predict(Ch03cs.AR, newdata = filter(customerAcquisition, acquisition == 1)),
reference = filter(customerAcquisition, acquisition == 1)$censor %>% factor,
positive = "1", mode = "everything")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 152 7
## 1 5 128
##
## Accuracy : 0.9589
## 95% CI : (0.9293, 0.9786)
## No Information Rate : 0.5377
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9173
## Mcnemar's Test P-Value : 0.7728
##
## Sensitivity : 0.9481
## Specificity : 0.9682
## Pos Pred Value : 0.9624
## Neg Pred Value : 0.9560
## Precision : 0.9624
## Recall : 0.9481
## F1 : 0.9552
## Prevalence : 0.4623
## Detection Rate : 0.4384
## Detection Prevalence : 0.4555
## Balanced Accuracy : 0.9582
##
## 'Positive' Class : 1
## qplot(`Observed Classes`, `Predicted Classes`,
data=bind_cols(`Observed Classes`= filter(customerAcquisition, acquisition == 1)$censor %>% factor,
`Predicted Classes` = predict(Ch03cs.AR, filter(customerAcquisition, acquisition == 1))),
colour= `Observed Classes`, geom = c("boxplot", "jitter"),
main = "Predicted Classes vs. Observed Classes", xlab = "Observed Classes", ylab = "Predicted Classes")
Качество классификационной модели после отбрасывания нескольких незначимых и малозначащих предикторов (acq_expense, ret_expense_sq, first_purchase, frequency_sq & employees) осталось весьма высоким, но после регуляризации можно ожидать, что это идеальное разбиение - это не просто побочный продукт данной выборки, но может реально содержаться в генеральной совокупности.
Второй этап расчета CLV - Пожизненная финансовая ценность клиента
Тем не менее авторы используют предсказания своей пробит-модели как отдельный предиктор imr_censor (\(\lambda\)) на втором этапе моделирования Пожизненная финансовая ценность клиента (CLV), когда обращаются в линейной регрессии. Эту регрессию строят методом наименьших квадратов на усеченной выборке, содержащей оставшихся (цензуированных) клиентов, т. е. clv > 0.
# Fit Linear Regression Model for Customer Lifetime Value by Authors
# SAS Code: imr_censor = pdf(’Normal’,xb_censor)/probnorm(xb_censor);
# customerAcquisition <- customerAcquisition %>%
# mutate(imr_censor = predict(Ch03cs.AR, newdata = customerAcquisition, type = "prob") %>% .[, "1"]) # Logit
xbeta <- predict(Ch03.censor, newdata = customerAcquisition, type = "link")
customerAcquisition <- customerAcquisition %>%
mutate(imr_censor = dnorm(xbeta) / pnorm(xbeta)) # Cumulative normal pdf
(Ch03.clv <- lm(clv ~ acq_expense + acq_expense_sq + ret_expense + ret_expense_sq + first_purchase + crossbuy +
frequency + frequency_sq + industry + revenue + employees + imr_censor,
data = filter(customerAcquisition, clv > 0) )) %>%
summary##
## Call:
## lm(formula = clv ~ acq_expense + acq_expense_sq + ret_expense +
## ret_expense_sq + first_purchase + crossbuy + frequency +
## frequency_sq + industry + revenue + employees + imr_censor,
## data = filter(customerAcquisition, clv > 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08356 -0.27824 0.01656 0.23935 1.12707
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.613e+00 9.317e-01 -2.805 0.005862 **
## acq_expense 8.700e-03 2.863e-03 3.039 0.002909 **
## acq_expense_sq -1.041e-05 2.649e-06 -3.929 0.000142 ***
## ret_expense 3.559e-03 4.035e-04 8.819 9.81e-15 ***
## ret_expense_sq -6.768e-07 1.117e-07 -6.058 1.58e-08 ***
## first_purchase 2.591e-03 1.082e-03 2.395 0.018130 *
## crossbuy 2.034e-01 2.013e-02 10.106 < 2e-16 ***
## frequency 1.552e-01 3.415e-02 4.545 1.30e-05 ***
## frequency_sq -4.854e-03 1.807e-03 -2.686 0.008247 **
## industry 5.664e-01 7.845e-02 7.220 4.86e-11 ***
## revenue 8.172e-03 3.696e-03 2.211 0.028885 *
## employees 1.693e-04 2.811e-04 0.602 0.548116
## imr_censor 1.738e-01 1.350e-01 1.287 0.200491
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.385 on 122 degrees of freedom
## Multiple R-squared: 0.8554, Adjusted R-squared: 0.8412
## F-statistic: 60.14 on 12 and 122 DF, p-value: < 2.2e-16writeLines("Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with others")## Variance Inflation Factors - if vif() > 2 - feature has multicollinearity with otherscar::vif(Ch03.clv) # Variance Inflation Factors - if vif() > 2 - feature has multicollinearity ## acq_expense acq_expense_sq ret_expense ret_expense_sq first_purchase crossbuy frequency frequency_sq
## 78.474401 81.385983 40.664391 40.003803 21.842561 1.064404 20.156080 19.876452
## industry revenue employees imr_censor
## 1.151873 3.431655 17.728522 1.300894Какие выводы могут получить авторы из этих результатов? Мы видим, что коэффициент imr_censor (\(\lambda\)) является положительным, но незначимым предиктором в модели CLV. Таким образом, авторы не могут интерпретировать это как означающее, что проблема выбора не существует, только маловероятно, что остатки уравнения бинарного выбора Censor не коррелируют с остатками уравнения линейной регрессии. Мы также видим, что за исключением Employees, остальные предикторы в модели CLV значимы при p < 0.05. Это означает, что мы обнаружили множество ключевых драйверов CLV для этого набора клиентов.
Как только мы предсказали значение CLV для каждой из перспективных компаний, следует сравнить это с фактическим значением из базы данных. Мы сравниваем значения CLV, где Censor == 1, то есть клиенты, которые все еще активные. вычислим среднее абсолютное отклонение (MAD или MAE) и среднюю абсолютную процентную погрешность (MAPE).
# Computing the Mean Absolute Deviation (MAD) and Mean Absolute Percent Error (MAPE)
with(filter(customerAcquisition, censor == 1), {
pred_clv <- predict(Ch03.clv) #, newdata = filter(customerAcquisition, censor == 1))
# mad = mean(abs(pred_clv - clv));
writeLines(sprintf("Mean Absolute Deviation (MAD): %.3f тыс. долл.", mean(abs(clv - pred_clv))))
# mape = mean(abs(clv - pred_clv)/clv);
writeLines(sprintf("Mean Absolute Percent Error (MAPE): %.3f %%", mean(abs(clv - pred_clv) / clv) * 100))
# mad1 = mean(abs(clv - mean(clv));
mad1 <- mean(abs(clv - mean(clv)))
writeLines(sprintf("Naive Mean Absolute Deviation (MAD1): %.3f тыс. долл.", mad1))
# mape1 = mean(abs(clv - mean(clv))/clv);
mape1 <- mean(abs(clv - mean(clv)) / clv) * 100
writeLines(sprintf("Naive Mean Absolute Percent Error (MAPE1): %.3f %%", mape1))
})## Mean Absolute Deviation (MAD): 0.290 тыс. долл.
## Mean Absolute Percent Error (MAPE): 4.517 %
## Naive Mean Absolute Deviation (MAD1): 0.758 тыс. долл.
## Naive Mean Absolute Percent Error (MAPE1): 11.857 %На этих данных авторы нашли, что MAD составляет 0.29 (или 290 долл. США), а MAPE - 4.52%. Для наивной тестовой модели авторы используют среднее значение CLV по цензуированным клиентам 6.579 тыс. долл. как прогнозируемое значение для каждого остающегося клиента. В этом случае получаем MAD1 - 758 долл. США и MAPE1 11.52%. Мы видим, что модель CLV от авторов обеспечивает лучшее предсказание CLV, чем наивная усредненная модель. Это показывает, что определение драйверов CLV поможет фирме глубже понять, какие клиенты, скорее всего, станут прибыльными в будущем.
Как это использовать?
Авторы находят, что Acq_Expense и Ret_Expense положительны связаны с уменьшающейся доходностью, что отмечено положительным коэффициентом на этих двух переменных и и отрицательными коэффициентами на их квадраты. Это говорит о том, что чем больше тратятся усилий, связанные с приобретением и удержанием, до предела, тем выше ожидаемая финансовая ценность жизни клиента, т. е. CLV. Для First_Purchase видно, что существует положительный коэффициент, предполагающий, что чем выше первоначальная сумма покупки, тем выше ожидаемая стоимость жизни клиента. Для Crossbuy наблюдается, что существует положительный коэффициент, предполагающий, что чем больше категорий продуктов приобретает покупатель, тем вполне логично выше ожидаемая стоимость жизни. Авторы замечают, что частота покупок Frequency положительна с уменьшающейся доходностью, что отмечено положительным коэффициентом на этом признаке и отрицательным коэффициентом на квадрат.
Это говорит о том, что клиенты, которые покупают последовательно с умеренной скоростью в течение срока их сотрудничества, скорее всего, имеют наивысшую ценность жизни. Затем через термины описания клиентов обнаруживается, что клиенты, которые являются B2B (Industry) и клиентами с более высоким уровнем дохода (Revenue), с большей вероятностью имеют более высокую жизненную ценность (CLV), чем клиенты, которые не входят в отрасль B2B и имеют более низкий доход.
Должен признаться, что крупнейшие вендоры уже прошли эти этапы анализа клиентов и подготовили преднастроенные модели для них. Примером в области телекоммуникации служит Oracle Communications Data Model Data Mining Models и связанная с ним отчетность Oracle Communications Data Model Sample Reports.
devtools::session_info()## Session info ----------------------------------------------------------------------------------------------------------## setting value
## version R version 3.5.0 (2018-04-23)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate Russian_Russia.1251
## tz Asia/Dhaka
## date 2018-07-20## Packages --------------------------------------------------------------------------------------------------------------## package * version date source
## abind 1.4-5 2016-07-21 CRAN (R 3.5.0)
## arm 1.10-1 2018-04-13 CRAN (R 3.5.0)
## assertthat 0.2.0 2017-04-11 CRAN (R 3.5.0)
## backports 1.1.2 2017-12-13 CRAN (R 3.5.0)
## base * 3.5.0 2018-04-23 local
## bindr 0.1.1 2018-03-13 CRAN (R 3.5.0)
## bindrcpp * 0.2.2 2018-03-29 CRAN (R 3.5.0)
## broom 0.4.4 2018-03-29 CRAN (R 3.5.0)
## car * 3.0-0 2018-04-02 CRAN (R 3.5.0)
## carData * 3.0-1 2018-03-28 CRAN (R 3.5.0)
## caret * 6.0-79 2018-03-29 CRAN (R 3.5.0)
## cellranger 1.1.0 2016-07-27 CRAN (R 3.5.0)
## class 7.3-14 2015-08-30 CRAN (R 3.5.0)
## cli 1.0.0 2017-11-05 CRAN (R 3.5.0)
## cmprsk 2.2-7 2014-06-17 CRAN (R 3.5.0)
## coda 0.19-1 2016-12-08 CRAN (R 3.5.0)
## codetools 0.2-15 2016-10-05 CRAN (R 3.5.0)
## colorspace 1.3-2 2016-12-14 CRAN (R 3.5.0)
## compiler 3.5.0 2018-04-23 local
## crayon 1.3.4 2017-09-16 CRAN (R 3.5.0)
## curl 3.2 2018-03-28 CRAN (R 3.5.0)
## CVST 0.2-1 2013-12-10 CRAN (R 3.5.0)
## data.table 1.11.2 2018-05-08 CRAN (R 3.5.0)
## datasets * 3.5.0 2018-04-23 local
## ddalpha 1.3.3 2018-04-30 CRAN (R 3.5.0)
## DEoptimR 1.0-8 2016-11-19 CRAN (R 3.5.0)
## devtools 1.13.5 2018-02-18 CRAN (R 3.5.0)
## digest 0.6.15 2018-01-28 CRAN (R 3.5.0)
## dimRed 0.1.0 2017-05-04 CRAN (R 3.5.0)
## dplyr * 0.7.5 2018-05-19 CRAN (R 3.5.0)
## DRR 0.0.3 2018-01-06 CRAN (R 3.5.0)
## e1071 1.6-8 2017-02-02 CRAN (R 3.5.0)
## evaluate 0.10.1 2017-06-24 CRAN (R 3.5.0)
## fitdistrplus * 1.0-9 2017-03-24 CRAN (R 3.5.0)
## forcats * 0.3.0 2018-02-19 CRAN (R 3.5.0)
## foreach 1.4.4 2017-12-12 CRAN (R 3.5.0)
## forecast 8.3 2018-04-11 CRAN (R 3.5.0)
## foreign 0.8-70 2017-11-28 CRAN (R 3.5.0)
## fracdiff 1.4-2 2012-12-02 CRAN (R 3.5.0)
## geometry 0.3-6 2015-09-09 CRAN (R 3.5.0)
## GGally 1.4.0 2018-05-17 CRAN (R 3.5.0)
## ggplot2 * 2.2.1 2016-12-30 CRAN (R 3.5.0)
## ggpubr * 0.1.6 2017-11-14 CRAN (R 3.5.0)
## glmnet 2.0-16 2018-04-02 CRAN (R 3.5.0)
## glue 1.2.0 2017-10-29 CRAN (R 3.5.0)
## gower 0.1.2 2017-02-23 CRAN (R 3.5.0)
## graphics * 3.5.0 2018-04-23 local
## grDevices * 3.5.0 2018-04-23 local
## grid 3.5.0 2018-04-23 local
## gridExtra 2.3 2017-09-09 CRAN (R 3.5.0)
## gtable 0.2.0 2016-02-26 CRAN (R 3.5.0)
## haven 1.1.1 2018-01-18 CRAN (R 3.5.0)
## highr 0.6 2016-05-09 CRAN (R 3.5.0)
## hms 0.4.2 2018-03-10 CRAN (R 3.5.0)
## htmltools 0.3.6 2017-04-28 CRAN (R 3.5.0)
## httr 1.3.1 2017-08-20 CRAN (R 3.5.0)
## ipred 0.9-6 2017-03-01 CRAN (R 3.5.0)
## iterators 1.0.9 2017-12-12 CRAN (R 3.5.0)
## jsonlite 1.5 2017-06-01 CRAN (R 3.5.0)
## kernlab 0.9-26 2018-04-30 CRAN (R 3.5.0)
## km.ci 0.5-2 2009-08-30 CRAN (R 3.5.0)
## KMsurv 0.1-5 2012-12-03 CRAN (R 3.5.0)
## knitr 1.20 2018-02-20 CRAN (R 3.5.0)
## labeling 0.3 2014-08-23 CRAN (R 3.5.0)
## lattice * 0.20-35 2017-03-25 CRAN (R 3.5.0)
## lava 1.6.1 2018-03-28 CRAN (R 3.5.0)
## lazyeval 0.2.1 2017-10-29 CRAN (R 3.5.0)
## lme4 1.1-17 2018-04-03 CRAN (R 3.5.0)
## lmtest 0.9-36 2018-04-04 CRAN (R 3.5.0)
## lubridate 1.7.4 2018-04-11 CRAN (R 3.5.0)
## magic 1.5-8 2018-01-26 CRAN (R 3.5.0)
## magrittr * 1.5 2014-11-22 CRAN (R 3.5.0)
## MASS * 7.3-49 2018-02-23 CRAN (R 3.5.0)
## Matrix 1.2-14 2018-04-13 CRAN (R 3.5.0)
## memoise 1.1.0 2017-04-21 CRAN (R 3.5.0)
## methods * 3.5.0 2018-04-23 local
## minqa 1.2.4 2014-10-09 CRAN (R 3.5.0)
## mnormt 1.5-5 2016-10-15 CRAN (R 3.5.0)
## ModelMetrics 1.1.0 2016-08-26 CRAN (R 3.5.0)
## modelr 0.1.2 2018-05-11 CRAN (R 3.5.0)
## munsell 0.4.3 2016-02-13 CRAN (R 3.5.0)
## nlme 3.1-137 2018-04-07 CRAN (R 3.5.0)
## nloptr 1.0.4 2017-08-22 CRAN (R 3.5.0)
## nnet 7.3-12 2016-02-02 CRAN (R 3.5.0)
## openxlsx 4.0.17 2017-03-23 CRAN (R 3.5.0)
## parallel 3.5.0 2018-04-23 local
## pillar 1.2.2 2018-04-26 CRAN (R 3.5.0)
## pkgconfig 2.0.1 2017-03-21 CRAN (R 3.5.0)
## plyr 1.8.4 2016-06-08 CRAN (R 3.5.0)
## prodlim 2018.04.18 2018-04-18 CRAN (R 3.5.0)
## psych 1.8.4 2018-05-06 CRAN (R 3.5.0)
## purrr * 0.2.4 2017-10-18 CRAN (R 3.5.0)
## quadprog 1.5-5 2013-04-17 CRAN (R 3.5.0)
## quantmod 0.4-13 2018-04-13 CRAN (R 3.5.0)
## R6 2.2.2 2017-06-17 CRAN (R 3.5.0)
## RColorBrewer 1.1-2 2014-12-07 CRAN (R 3.5.0)
## Rcpp 0.12.17 2018-05-18 CRAN (R 3.5.0)
## RcppRoll 0.2.2 2015-04-05 CRAN (R 3.5.0)
## readr * 1.1.1 2017-05-16 CRAN (R 3.5.0)
## readxl 1.1.0 2018-04-20 CRAN (R 3.5.0)
## recipes 0.1.2 2018-01-11 CRAN (R 3.5.0)
## reshape 0.8.7 2017-08-06 CRAN (R 3.5.0)
## reshape2 1.4.3 2017-12-11 CRAN (R 3.5.0)
## rio 0.5.10 2018-03-29 CRAN (R 3.5.0)
## rlang 0.2.0 2018-02-20 CRAN (R 3.5.0)
## rmarkdown 1.9 2018-03-01 CRAN (R 3.5.0)
## robustbase 0.93-0 2018-04-24 CRAN (R 3.5.0)
## rpart 4.1-13 2018-02-23 CRAN (R 3.5.0)
## rprojroot 1.3-2 2018-01-03 CRAN (R 3.5.0)
## rstudioapi 0.7 2017-09-07 CRAN (R 3.5.0)
## rvest 0.3.2 2016-06-17 CRAN (R 3.5.0)
## scales 0.5.0 2017-08-24 CRAN (R 3.5.0)
## sfsmisc 1.1-2 2018-03-05 CRAN (R 3.5.0)
## splines * 3.5.0 2018-04-23 local
## stats * 3.5.0 2018-04-23 local
## stats4 * 3.5.0 2018-04-23 local
## stringi 1.1.7 2018-03-12 CRAN (R 3.5.0)
## stringr * 1.3.1 2018-05-10 CRAN (R 3.5.0)
## survival * 2.41-3 2017-04-04 CRAN (R 3.5.0)
## survminer * 0.4.2 2018-01-31 CRAN (R 3.5.0)
## survMisc 0.5.4 2016-11-23 CRAN (R 3.5.0)
## tibble * 1.4.2 2018-01-22 CRAN (R 3.5.0)
## tidyr * 0.8.1 2018-05-18 CRAN (R 3.5.0)
## tidyselect 0.2.4 2018-02-26 CRAN (R 3.5.0)
## tidyverse * 1.2.1 2017-11-14 CRAN (R 3.5.0)
## timeDate 3043.102 2018-02-21 CRAN (R 3.5.0)
## tools 3.5.0 2018-04-23 local
## tseries 0.10-44 2018-04-15 CRAN (R 3.5.0)
## TTR 0.23-3 2018-01-24 CRAN (R 3.5.0)
## urca 1.3-0 2016-09-06 CRAN (R 3.5.0)
## uroot 2.0-9 2017-01-29 CRAN (R 3.5.0)
## utils * 3.5.0 2018-04-23 local
## VGAM * 1.0-5 2018-02-07 CRAN (R 3.5.0)
## withr 2.1.2 2018-03-15 CRAN (R 3.5.0)
## xml2 1.2.0 2018-01-24 CRAN (R 3.5.0)
## xtable 1.8-2 2016-02-05 CRAN (R 3.5.0)
## xts 0.10-2 2018-03-14 CRAN (R 3.5.0)
## yaml 2.1.19 2018-05-01 CRAN (R 3.5.0)
## zoo 1.8-1 2018-01-08 CRAN (R 3.5.0)Выводы по главе 3
Цель этой главы состояла в том, чтобы изучить текущие модели привлечения клиентов и представить некоторые эмпирические примеры того, как фирмы могут применять эти знания. Приобретение клиентов является первым ключевым шагом в процессе CRM. Авторы также показали, что когда фирмы способны участвовать в оптимальном выборе перспективных потребителей, понимая движущие силы ведущие к клиентам и «правильного» объема усилий по их подключению к сотрудничеству, осознавая взаимосвязь между маркетинговыми расходами и стоимостью, которую может принести клиент. Получаемый в ходе этого процесса результат может обеспечить существенную рентабельность клиентов и как следствие значительную прибыльность фирмы.