python数据分析实战：生存分析与电信用户流失预测

1.背景

1.1 生存分析、KM曲线及Cox回归

常见的回归模型聚焦在事件结果与影响因素上，生存分析既关注结果又关注发生事件。既研究结果影响因素，又研究影响因素与结果出现事件长短的关系，是研究生存现象和发生事件关系及统计规律的一门学科。
什么是生存分析？这篇文章介绍了生存分析中的一些基本概念。风险函数$h(t)$指在$t$事件之前未发生任何事件而在$t$发生事件的概率，累积风险函数$H(t)={\int }_0^{t}h(u)du=-ln[S(t)]$，生存函数$S(t)=e^{-H(t)}=e^{-{\int }_0^{t}h(u)du}$表示研究对象到该时刻仍未发生事件的概率。
目前，对生存函数及风险函数的刻画方法分为非参数和参数两类。其中，最常用的是非参数方法中对生存函数进行刻画的Kaplan-Meier曲线法和对累积风险函数进行刻画的Nelson-Aalen曲线法。KM是一类非参数估计方法，不对数据分布做任何假设，而是直接用概率乘法定理估计生存率。
KM方法的优势在于能直观地观察生存曲线，便于不同生存曲线之间进行简单比较，但无法建立数学模型对多个影响因素进行分析。为了解决这个问题，Cox比例风险回归模型诞生了，实现了对包含分类变量及连续变量的多变量回归。

1.2 案例背景

数据来源于kaggle，中文名称为电信用户流失情况表，也被网上很多文章当作案例使用。下面为字段说明，事件为Churn流失，生存时间为tenure使用月数。

接下来我们对这个数据利用Cox比例风险回归模型进行分析。Cox模型使用需要满足比例风险假定Hazard Ratio（HR）：各危险因素的作用不随时间的变化而变化，即不随时间的变化而变化。但一般情况下是不用考虑的，官网文档给出了解释我们需不需要考虑PH假定？。
python中支持生存分析的包为lifelines ，可以支持生存函数的绘制、对数秩检测以及Cox模型拟合。此外还要介绍一个包pandas_profiling，这个包能快速对数据进行描述性统计，包含了数据的分布、均值最大值最小值、缺失值以及相关度。

2.AIC向前逐步回归法进行特征选择

在进行后续操作之前，依然要进行数据清洗及哑变量处理，在此不赘述了。这里介绍一种新的特征选择方式：AIC向前逐步回归法。

逐步回归分为三种，分别是向前逐步回归，向后逐步回归，逐步回归。向前逐步回归的特点是将自变量一个一个当如模型中，每当放入一个变量时，都利用相应的检验准则检验，当加入的变量不能使得模型变得更优良时，变量将会被剔除，如此不断迭代，直到没有适合的新变量加入为止。向后逐步回归的特点是，将所有变量都放入模型之后，一个一个的剔除变量，将某一变量拿出模型而使得模型更优良时，将会剔除此变量。如此反复迭代，直到没有合适的变量剔除为止。逐步回归则是结合了以上的向前和向后逐步回归的特点。

AIC即赤池值，是衡量模型拟合优良性和模型复杂性的一种标准，在建立多元线性回归模型时，变量过多，且有不显著的变量时，可以使用AIC准则结合逐步回归进行变量筛选。AIC数学表达式如下：
$$AIC=2p+n(log(SSE/n))AIC=2p+n(log(SSE/n))$$
AIC越小我们认为模型更优良，下面以AIC作为判定标准进行向前逐步回归，对特征进行筛选。但是向前逐步回归函数与每次选择的模型有关，并没有现成的包使用。上面的文章是以线性回归为例写的代码，所以我们要先依据Cox模型对其进行改写。

import pandas as pd
import pandas_profiling
import numpy as np
from sklearn import metrics
from sklearn.model_selection import train_test_split
from lifelines import CoxPHFitter
import matplotlib.pyplot as plt
#定义向前逐步回归函数
def forward_select(data):
    variate=set(data.columns)  #将字段名转换成字典类型
    variate.remove('tenure')  #去掉因变量的字段名
    variate.remove('Churn')
    selected=[]
    current_score,best_new_score=float('inf'),float('inf')  #目前的分数和最好分数初始值都为无穷大（因为AIC越小越好）
    #循环筛选变量
    while variate:
        aic_with_variate=[]
        for candidate in variate:  #逐个遍历自变量
            df=data[[candidate]+list(selected)+['tenure','Churn']]
            cph=CoxPHFitter(penalizer=0.01)
            cph.fit(df,duration_col='tenure',event_col='Churn',step_size=0.01)
            aic=cph.AIC_partial_
            aic_with_variate.append((aic,candidate))  #将第每一次的aic值放进空列表
        aic_with_variate.sort(reverse=True)  #降序排序aic值
        best_new_score,best_candidate=aic_with_variate.pop()  #最好的aic值等于删除列表的最后一个值，以及最好的自变量等于列表最后一个自变量
        if current_score>best_new_score:  #如果目前的aic值大于最好的aic值
            variate.remove(best_candidate)  #移除加进来的变量名，即第二次循环时，不考虑此自变量了
            selected.append(best_candidate)  #将此自变量作为加进模型中的自变量
            current_score=best_new_score  #最新的分数等于最好的分数
            print("aic is {},continuing!".format(current_score))  #输出最小的aic值
        else:
            print("for selection over!")
            break
    return(selected)
#拆分训练集和测试集
train, test = train_test_split(df,test_size=0.2)
#进行特征筛选
selected=forward_select(train)

筛选出的特征：

['Contract_Month-to-month',
 'TotalCharges',
 'MonthlyCharges',
 'Contract_Two year',
 'OnlineSecurity_No',
 'InternetService_Fiber optic',
 'DeviceProtection_No internet service',
 'PaymentMethod_Bank transfer (automatic)',
 'PaymentMethod_Credit card (automatic)',
 'Partner_No',
 'OnlineBackup_No',
 'TechSupport_No',
 'PaperlessBilling_No',
 'Dependents_No',
 'MultipleLines_No',
 'StreamingMovies_Yes',
 'MultipleLines_No phone service',
 'StreamingTV_Yes',
 'InternetService_No']

3.Cox模型搭建

3.1 特征重要性分析

data=train[selected+['tenure','Churn']]
#实例化模型
cph=CoxPHFitter(penalizer=0.01)
#数据拟合
cph.fit(data,duration_col='tenure',event_col='Churn',step_size=0.5)
#打印结果
cph.print_summary()
#绘图
fig,ax = plt.subplots(figsize=(12,9))
cph.plot(ax=ax)
plt.show()

拟合结果


Cox数学表达式：
$$h(t;x_1,x_2,...,x_p)=h_0(t)exp({\beta }_1{x}_1+{\beta }_2{x}_2+...+{\beta }_p{x}_p)$$
即准风险率$h_0(t)$为所有自变量取值均为0的风险率，$({\beta }_1,{\beta }_2,...,{\beta }_p)$是模型的偏自回归系数，反应了自变量对风险率的影响。与一般的回归分析不同，Cox模型不是直接用生存时间作为回归方程的因变量，自变量对生存时间的影响通过风险函数与基准函数的比值来表示。
我们对模型的解读为：

• 当自变量$x_p$为连续型变量时，${\beta }_p$为在除$x_p$之外的其他变量均不变的情况下，$x_p$每增加一个单位，风险率变化$exp({\beta }_p)$倍。
• 当自变量$x_p$为离散型变量时，如有0和1两个分类，${\beta }_p$为在除$x_p$之外的其他变量均不变的情况下，类别1相对类别0风险率是$exp({\beta }_p)$倍。

总体来说，对于影响因子的解读可以分为以下三类：

• $\beta >0$，风险率随$x$增大而增大，为危险因素。
• $\beta <0$，风险率随$x$增大而减小，为保护因素。
• $\beta =0$，风险率不随$x$变化，为无关因素。

这样一来我们便得到了因子的分类，危险因素会加速客户流失，保护因素会减缓客户流失：

3.2 模型校准

校准？总结起来就是模型结果预测概率和真实的经验概率保持一致。

模型校准就是要让模型结果预测概率和真实的经验概率保持一致。说人话也就是，在一个二分类任务中取出大量（M个）模型预测概率为0.6的样本，其中有0.6M个样本真实的标签是1。总结一下，就是模型在预测的时候说某一个样本的概率为0.6，这个样本就真的有0.6的概率是标签为1。
上面是一个正面的例子，下面我再来举一个反面的例子说明模型校准的重要性。还是在一个二分类任务中取出大量（M个）模型预测概率为0.6的样本，而这些样本的真实标签全部都是1。虽然从accuracy的角度来考察，模型预测样本概率为0.6最后输出时会被赋予的标签就是1，即accuracy是100%。但是从置信度的角度来考察，这个模型明显不够自信，本来这些全部都是标签为1的样本，我们肯定希望这个模型自信一点，输出预测概率的时候也是1。

下面我们以训练的模型与之前分割出的测试集做校准分析：

from lifelines.calibration import survival_probability_calibration
fig,ax = plt.subplots(figsize=(12,9))
survival_probability_calibration(cph, test[selected+['tenure','Churn']], t0=40, ax=ax)
ax.grid() 

横坐标为预测的事件发生率（Predicted risk），纵坐标是观察到的实际事件发生率（Observed
risk），范围均为0到1，可以理解为事件发生率（百分比）。对角线的虚线是参考线，即预测值=实际值的情况，红线是曲线拟合线。

• 如果预测值=观察值，则红线与参考线完全重合。
• 如果预测值＞观察值，即高估了风险，则红线在参考线下面。
• 如果预测值＜观察值，即低估了风险，则红线在参考线上面。

以时长<=40为例，我们的图形十分接近对角线，整体来说是一个不错的校准。前8个月风险估计比较准确，8个月后风险被低估。

3.3 对个体进行预测

对测试集的一些个体进行预测，并绘制生存线，数字代表样本序号。

surv_hat=cph.predict_survival_function(test[selected].reset_index())
fig,ax=plt.subplots(figsize=(12,9))
for i in (5,10,50,100,150,200,500,1000):
    surv_hat[i].plot(label=str(i))
plt.legend()
plt.show()

3.3 用户流失预测

半衰期也可以称为中位生存时间。由于删失的存在，无法以平均时间反应时间发生的时间水平，因此生存分析在描述生存时一般以半衰期，也就是存活概率50%时对应消耗的时间来表述。
下面以个体半衰期作为预测生存时间，得到各个时间结点流失的用户占比。

def count_existing_users(df,timecut):
    s=[]
    for i in timecut:
        l=df.iloc[i,:]
        s.append(len(l[l<0.5])/len(l))
    result=pd.DataFrame({'Month':timecut,'Lost Users Ratio':s})
    return result
count_existing_users(surv_hat,[10*i for i in range(8)])    

4.总结

我们分别利用Cox比例风险回归模型比较了不同特征对于客户流失的影响，并且预测了个体客户的生存时长，计算了各时间点用户流失的比例，结果具有一定的参考性。