独家 | 规范性分析的实用介绍（附R语言案例研究&演示代码）

翻译：张睿毅

校对：吴金笛

文章来源：微信公众号 数据派THU

本文约4500字，建议阅读15分钟。

本文是作为Analytics Vidhya’s Internship Challenge的一部分提交的案例研究。

Analytics Vidhya’s Internship Challenge

https://datahack.analyticsvidhya.com/contest/av-internship-challenge/?utm_source=blog&utm_medium=practical-introduction-prescriptive-analytics

引言

“分析的不同分支是什么？“当我们开始分析之旅时，我们大多数人都被告知有两种类型——描述性分析和预测性分析。实际上还有第三个经常被忽视的分析——规范性分析。

处方分析是这三种方法中最强大的一种。让我用一个例子来演示一下。

最近，一场致命的飓风袭击了印度的奥迪沙，但幸运的是，大多数人已经撤离。奥迪沙气象部门已经预测到了这场可怕的气旋到来，并做出了拯救生命的决定，疏散可能有危险的地区。

与1999年相比，当时有10000多人死于类似的飓风。由于没有关于即将来临的暴风雨的预报，所以他们没有意识到将发生什么。那么，发生了什么变化？

奥迪沙政府是规范性分析的受益者。他们能够利用气象部门对气旋的准确预测的服务——它们的路径、强度和时间。他们用这个来决定何时以及需要做什么来防止生命的损失。

因此，在本文中，我们将首先了解“规定性分析”这一术语的含义。然后，我们将通过进行案例研究和实施分析分支（描述性、预测性和规定性）来巩固我们的学习。我们开始吧！

目录表

• 什么是规范性分析？
• 设置问题陈述
• 生成假设
• 制定我们的模型构建方法
• 数据可视化和数据准备–描述性分析
• 预测客户行为-预测分析
• 推荐提升性能—规范性分析

什么是规范性分析

我们可以将分析大致分为三个不同的部分——描述性分析、预测性分析和规范性分析。让我们来看看这些：

• 描述性分析是任何模型构建练习的第一部分。我们对历史数据进行分析，以确定因变量和自变量的模式和趋势。这一阶段也有助于假设生成、变量转换和特定行为模式的任何根本原因分析。
• 预测分析是分析的下一个阶段。这里，我们利用已清理和/或转换的数据，并在该数据上拟合一个模型，以预测因变量的未来行为。预测分析解决了可能发生的问题。
• 规规范性分析是最后一个阶段，预测用于规定（或建议）下一组要做的事情。这就是我们奥迪沙政府的例子来源。他们利用气象部门的预测，采取了一系列措施，如安置低洼地区的所有人员，提前安排食物、住所和医疗救助等，以确保损失有限。

下图很好地说明了规范性分析框架下的组件：

设置问题陈述

我发现学习一个主题的最好方法是实践它。所以，让我们通过一个案例研究和实现上面讨论的每个分析部分来理解规范性分析。

电信运营商组织的高级管理层担心客户流失率不断上升。此外，最近的一项独立调查表明，整个行业将面临不断增长的客户流失率和ARPU（每单位平均收入）的下降。

到目前为止，留住客户的努力是非常被动的。只有当客户要求关闭他们的帐户时，我们才会采取行动。这不是一个好策略，是吗？管理团队希望在这方面采取更积极的措施。

作为数据科学家，我们的任务是分析他们的数据，获得见解，预测客户的潜在行为，然后推荐改进性能的步骤。

为我们的问题准备数据集

您可以从这里下载数据集。我还提供了我的Github存储库上的完整代码。有三个R文件，您应该按以下顺序使用它们：

• DataPreparation.r
• Visualization.r
• ModelBuilding.r

假设生成

生成一个假设是解锁任何数据科学或分析项目的关键。我们应该首先列出我们试图通过我们的方法实现的目标，然后从那里着手。

客户流失是由以下因素驱动的（根据独立行业调查）：

• 成本和账单
• 网络与服务质量
• 数据使用连接问题

我们希望对我们的电信提供商进行同样的测试。通常，我们鼓励公司提出一套详尽的假设，以免遗漏任何变量或要点。但是，对于本文的范围，我们将把重点缩小到一个：

与成本、账单、网络和服务质量相关的变量是否对客户留下或离开服务提供商的决定做出了重大贡献？

制定我们的模型构建方法

现在我们有了数据集、问题陈述和假设要测试，是时候让我们开始做事情了。让我们仔细研究一下数据，看看能得出什么样的见解。

我在下面的插图中总结了我的方法。通常，任何模型构建练习都将经历类似的步骤。请注意，这是我的方法——您可以改变内容，并在您的终端上利用数据。例如，我们正在删除丢失值超过30%的变量，但您可以对此进行自己的调用。

下面是查找丢失值超过30%的变量的代码：

mydata=read.csv("Telecom_Sampled.csv")

mydata$churn=as.factor(mydata$churn)

anyNA(mydata)

Percentage_missing=round(colSums(is.na(mydata[,colnames(mydata)[colSums(is.na(mydata))>0]]))/nrow(mydata)*100,2)

data.frame(MissingProportion=sort(Percentage_missing, decreasing = TRUE))

 

#Finding variable names with more than 30% missing values

Variables_with_High_NAs=colnames(mydata)[colSums(is.na(mydata))/nrow(mydata)>0.3]

 

#Removing the variables with more than 30% missing values from the original dataset

mydata=mydata[,!names(mydata)%in%Variables_with_High_NAs]

 

#13 variables removed

如上图所示，我们删除了所有丢失值超过30%的变量。以下是我们数据集的总结：

我们把变量的数量从82减少到69。

数据可视化和数据分析准备——规范性分析

让我们做一个单变量、双变量和多变量分析，分析各种自变量和目标变量。这应该能让我们了解客户流失的影响。我在下面分享了一些可视化效果。您可以在Github存储库中找到整个探索性分析。

让我们从绘制三个图开始（输出在代码块下面）：

#Univariate Analysis & Multivariate Analysis

a=ggplot(Telecom_Winsor, aes(x=mou_Mean, y=..density.., fill=1))

a+geom_histogram(stat = "bin", bins = 15)+geom_density(alpha=0.5)+

guides(fill=FALSE)+labs(y="Density", title="Density graph - MOU_Mean Vs Churn")+

theme_bw()+facet_grid(~churn)+theme(plot.title = element_text(size = 10, hjust = 0.5))



a=ggplot(Telecom_Winsor, aes(x=totmrc_Mean, y=..density.., fill=1))

a+geom_histogram(stat = "bin", bins = 15)+geom_density(alpha=0.5)+

guides(fill=FALSE)+labs(y="Density", title="Density graph - totmrc_Mean Vs Churn")+

theme_bw()+facet_grid(~churn)+theme(plot.title = element_text(size = 10, hjust = 0.5))

 

a=ggplot(Telecom_Winsor,aes(x=F_eqpdays), alpha=0.5)

a+geom_bar(stat = "count", aes(fill=models), position = "dodge")+

facet_grid(~churn)+labs(x="", fill="Models",y="Count", title="F_eqpdays Impact Churn?")+

theme(legend.position = c(0.8,0.8), plot.title = element_text(size = 10, hjust = 0.5),

legend.key.size = unit(0.5,"cm"), legend.title = element_text(size = 8),

legend.text = element_text(size = 8), axis.text.x = element_text( angle=45,size = 8, vjust = 1,

hjust = 1),

legend.direction = "horizontal")

首先，我们将根据目标变量分析平均使用分钟数、收入范围、平均每月总经常性费用和平均中断或阻塞呼叫数-客户流失：

同样，我们将分析掉话（失败）的平均数、客户生命周期内的总呼叫数、出站无线到无线语音呼叫数的范围以及针对客户流失变量等待的平均呼叫数：

让我们改变一下。我们将使用Awesome ggplot2软件包中的刻面功能，根据客户流失变量绘制使用月份、信用等级代码、呼叫中断和当前设备的天数：

我们将分别分析数值变量，看看是否有共线度高的特征。这是因为共线变量的存在总是会降低模型的性能，因为它们在模型中引入了偏差。

我们应该处理共线性问题。目前，有许多方法可以解决这一问题，例如使用主成分分析（PCA）进行变量变换和简化。我删除了高度相关的变量：

###Finding highly correlated variables

Numeric=Filter(is.numeric,mydata_Rev)

library(corrplot)

M=cor(na.omit(Numeric))

corrplot(M, method = "circle", type = "lower",

tl.srt = 45, tl.col = "black", tl.cex = 0.75)

客户行为预测-预测分析

这是大多数人都熟悉的部分——根据训练数据构建模型。我们将构建一些模型，以便在整个频谱中比较它们的性能。

从简单的线性模型到复杂的非参数和非线性模型，对多个模型进行训练是一种很好的实践。模型的性能取决于依赖变量和独立变量之间的关系。如果关系是线性的，那么简单的模型会产生好的结果（而且它们更容易解释）。

或者，如果关系是非线性的，复杂的模型通常会给出更好的结果。随着模型复杂度的增加，模型引入的偏差减小，方差增大。对于我们的问题，我们将在训练集上构建大约10个模型，并根据未公开的测试数据对其进行验证。

我们将建立的模型是：

• 具有不同分类阈值的逻辑回归和判别分析等简单模型
• 使用合成少数族过采样技术平衡数据集后的随机森林（SMOTE）
• 五个独立模型的集合，并通过平均单个输出概率预测输出
• XGBoost 算法

以下是逻辑回归模型的代码（您可以使用我的Github存储库中提供的代码来尝试其他代码）：

LGM1=glm(churn~., data = Telecom_Winsor, family = "binomial")

summary(LGM1)

#Remove hnd_wecap since it did not seem to be a significant and is introducing NAs in the model

names(Telecom_Winsor)

Telecom_Winsor_Lg=subset(Telecom_Winsor,select = -hnd_webcap)

Telecom_Winsor_Lg=droplevels(Telecom_Winsor_Lg)

 

#Data Splitting

library(caret)

set.seed(1234)

 

Index=createDataPartition(Telecom_Winsor_Lg$churn,times = 1,p=0.75,list = FALSE)

Train=Telecom_Winsor_Lg[Index,]

Test=Telecom_Winsor_Lg[-Index,]

 

prop.table(table(Train$churn))

 

LGM1=glm(churn~., data = Train, family = "binomial")

summary(LGM1)

step(LGM1, direction = "both")

 

library(car)

 

LGMF=glm(formula = churn ~ mou_Mean + totmrc_Mean + rev_Range + drop_blk_Mean + 

           drop_vce_Range + callwait_Mean + callwait_Range + ccrndmou_Range + 

           adjqty + rev_Mean + ovrmou_Mean + avgqty + age1 + age2 + 

           hnd_price + actvsubs + uniqsubs + datovr_Range + adjmou + 

           adjrev + plcd_dat_Mean + crclscod + asl_flag + mouR_Factor + 

           change_mF + F_months + F_eqpdays + F_iwylis_Vmean, family = "binomial", 

         data = Train)

car::vif(LGMF)

#adjmou, avgqty and adjqty have very high VIF but with Df considered GVIF is low enough. Hence no requirement for further collinesrity treatment 

 

summary(LGMF)

 

Pred=predict(LGMF, Test, type = "response")

options(scipen = 9999)

L=data.frame(LogOfOdds=round(exp(coef(LGMF)),3))

L$Variable=row.names(L)

row.names(L)=NULL

L=L[,c(2,1)]

 

#The variables which if undergo a change of 1 unit there is more than 50% probability of the customer decision changing from terminating service to staying with the servicer

 

L%>%arrange(desc(LogOfOdds))%>%filter(LogOfOdds>=1)%>%mutate(Probability=round(LogOfOdds/(1+LogOfOdds),3))

 

Pred.class=ifelse(Pred>0.24,1,0)

CM=confusionMatrix(as.factor(Pred.class),Test$churn)

CM$table

fourfoldplot(CM$table)

Acc_Log24=CM$overall[[1]]

Sensitivity_Log24=CM$byClass[[1]]

Specificity_Log24=CM$byClass[[2]]

F1sq_Log24=CM$byClass[[7]]

library(ROCR)

Pred.Storage=prediction(Pred,Test$churn)

 

AUC=performance(Pred.Storage,"auc")

AUC_Log24=AUC@y.values[[1]]

perf=performance(Pred.Storage,"tpr","fpr")

 

###########################

 

cut_offs=data.frame(cut=perf@alpha.values[[1]], fpr=perf@x.values[[1]], tpr=perf@y.values[[1]])

cut_offs=cut_offs[order(cut_offs$tpr, decreasing = TRUE),]

library(dplyr)

cut_offs%>%filter(fpr<=0.42,tpr>0.59)

 

#cutoff of 0.2356 seems to give the highest tpr and relatively low fpr

 

Pred.class235=ifelse(Pred>0.235,1,0)

CM=confusionMatrix(as.factor(Pred.class235),Test$churn)

fourfoldplot(CM$table)

Acc_Log23.5=CM$overall[[1]]

Sensitivity_Log23.5=CM$byClass[[1]]

Specificity_Log23.5=CM$byClass[[2]]

 

AUC=performance(Pred.Storage,"auc")

AUC_Log23.5=AUC@y.values[[1]]

F1sq_Log23.5=CM$byClass[[7]]

############

Pred.class=ifelse(Pred>0.25,1,0)

CM=confusionMatrix(as.factor(Pred.class),Test$churn)

fourfoldplot(CM$table)

Acc_Log25=CM$overall[[1]]

Sensitivity_Log25=CM$byClass[[1]]

Specificity_Log25=CM$byClass[[2]]

F1sq_Log25=CM$byClass[[7]]

AUC=performance(Pred.Storage,"auc")

AUC_Log25=AUC@y.values[[1]]

################################

Pred.class=ifelse(Pred>0.26,1,0)

CM=confusionMatrix(as.factor(Pred.class),Test$churn)

fourfoldplot(CM$table)

Acc_Log26=CM$overall[[1]]

Sensitivity_Log26=CM$byClass[[1]]

Specificity_Log26=CM$byClass[[2]]

F1sq_Log26=CM$byClass[[7]]

#Choice of cutoff at 24,25, 26 results in increasing accuracy and sensitivity but decreasing Specificity

AUC=performance(Pred.Storage,"auc")

AUC_Log26=AUC@y.values[[1]]

以下是对我们模型的评估比较：

与其他模型相比，逻辑回归似乎给出了最好的结果。LG_26是一个逻辑回归模型，阈值为26%。如果你在这方面有所提高，请告诉我——我很高兴听到你对如何处理这个问题的想法。

改进绩效的建议—规范性分析

现在，我们一直在等待的部分——规范性分析！让我们看看我们可以提出哪些建议来改进模型的性能。

在下面的图片中，我们列出了在各自独立变量中，每1个单位变化就有超过50%的概率改变客户决策的变量。这一见解是由我们上面看到的逻辑回归模型产生的。这本质上是因变量的概率对数与自变量之间的关系。

因此，如果我们计算因变量系数的指数，我们就得到了概率，从中我们得到了自变量中一个单位变化的顾客行为变化的概率（使用公式概率=赔率/（1+赔率））。

下图将使您更好地了解我所说的内容：

还记得我们之前通过独立调查得出的假设吗？这也是事实。以下来自物流模型的汇总统计数据证明：

以下是我们从分析中得出的结论：

影响成本和计费的变量非常重要。

Adjmou是前5个优势比之一。

月平均总经常性收费（totmrc_Mean）、收入（charge amount）、范围（rev_Range）、adjmou（Billing adjustments）等被发现具有高度显著性。这表明成本和账单会影响客户行为。

同样，网络和服务质量变量，如drop_bkl_mean（平均掉话和阻塞的呼叫数）也非常显著。Datovr_Range（数据覆盖的收入范围）并不显著，但其优势比大于1，表明1个单位的价值变化有50%以上的机会将客户行为从一个级别改变到另一个级别。也许我们需要注意一下。

此外，截距也很重要。这构成了模型删除的分类变量水平的影响。

推荐

让我们根据我们了解的内容来写下我们的建议。

建议将速率计划迁移作为主动保留策略

Mou_Mean（使用分钟数）是最重要的变量之一。因此，积极主动地与客户合作，以增加他们的MOU，使其保留更长时间是有意义的。

此外，mouR_Factor是非常显著的。记住，这是mou_Range的派生变量。

MOU的变化也非常显著。Change_mF是change_mou的派生变量。

为了补充上述内容，我们还发现ovrmou_Mean也是一个非常显著的变量，其优势比大于1。该变量对系数有一个正估计，表明过度流失的增加。

如果我们公司能与客户合作，这会有所帮助。根据它们的使用情况，我们可以将它们迁移到最佳计划费率，以避免过多的费用。

对顾客制定的主动保留策略

确定客户流失率最高的客户，并为他们制定主动的保留策略。可能问如果预算有限怎么办？然后，公司可以建立一个提升图，并通过接触目标客户来优化其保留工作：

在这里，该模型拥有30%的总客户群，准确地提供了33%的潜在客户流失候选人：

所取得的提升将有助于我们接触到客户流失的候选人，目标是公司的总客户群要少得多。还要注意前30个十分位数是如何给我们带来最高增益的。这会给我们大约33%的可能终止服务的客户。

简单地说，公司选择了整个客户数据库的30%，该数据库覆盖了33%可能离职的人员。这比随机给客户打电话要好得多，因为客户可能会对所有潜在的客户流失候选人造成15%的点击率。

您可以使用下面的代码来测试模型，方法是识别20%需要积极合作以防止客户流失的客户：

gains(as.numeric(Telecom_Winsor$churn),predict(LGMF,type="response",newdata=Telecom_Winsor[,-42])

,groups = 10)



Telecom_Winsor$Cust_ID=mydata$Customer_ID



Telecom_Winsor$prob<-predict(LGMF,type="response",newdata=Telecom_Winsor[,-42])



quantile(Telecom_Winsor$prob,prob=c(0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90,1))



targeted=Telecom_Winsor%>%filter(prob>0.3224491 & prob<=0.8470540)%>%dplyr::select(Cust_ID)

客户流失率大于32.24%，小于84.7%。ModelBuilding.r代码将帮助您处理上述代码块的逻辑流。

最后几点

如果公司能够恰当地利用规范性分析，那么它真的是一件了不起的事情。就分析的三个分支而言，这仍然是人们关注的焦点。

但是，随着我们在分析层级上的不断提升，规范性分析是最受欢迎的领域，因为它可以帮助组织以相当的信心规划和准备未来。

规定性分析旨在确定各种选择中的最佳解决方案或结果。请记住，我们不能分别地分析三个分支。我们需要先做描述性的和预测性的工作，然后再进行规范性的工作。

原文标题:

A Practical Introduction to Prescriptive Analytics (with Case Study in R)

原文链接：

https://www.analyticsvidhya.com/blog/2019/05/practical-introduction-prescriptive-analytics/

编辑：王菁

校对：杨学俊

译者简介

张睿毅，北京邮电大学大二物联网在读。我是一个爱自由的人。在邮电大学读第一年书我就四处跑去蹭课，折腾整一年惊觉，与其在当下焦虑，不如在前辈中沉淀。于是在大二以来，坚持读书，不敢稍歇。资本主义国家的科学观不断刷新我的认知框架，同时因为出国考试很早出分，也更早地感受到自己才是那个一直被束缚着的人。太多真英雄在社会上各自闪耀着光芒。这才开始，立志终身向遇到的每一个人学习。做一个纯粹的计算机科学里面的小学生。

翻译组招募信息

工作内容：将选取好的外文前沿文章准确地翻译成流畅的中文。如果你是数据科学/统计学/计算机专业的留学生，或在海外从事相关工作，或对自己外语水平有信心的朋友，数据派翻译组欢迎你们加入！

你能得到：提高对于数据科学前沿的认知，提高对外文新闻来源渠道的认知，海外的朋友可以和国内技术应用发展保持联系，数据派团队产学研的背景为志愿者带来好的发展机遇。

其他福利：和来自于名企的数据科学工作者，北大清华以及海外等名校学生共同合作、交流。