2.3 生存分析在二手车定价案例中的应用
前文详细阐述了生存分析的理论，本节将基于2.1节提到的二手车定价案例，带领大家通过Python代码实操完成生存分析及结果解读。在开始实操之前，我们先对二手车定价案例涉及的具体问题及数据进行界定，并对操作流程做简要介绍。
1. 问题聚焦
我们以帕萨特二手车定价为例，以90天为一个销售周期，如果二手车超过90天仍未销售出去，则车辆的收益为0，成本为90天消耗的仓储及维护总成本。在2.1节推导出的利润公式中，所有车辆购入成本Cp固定不变，为常数，为了简化问题，此处不考虑购入成本，我们只需要求得如下公式的最大值：
2. 数据获取
假设我们从平台上获得的数据信息如表2-1所示。
表2-1 二手车定价案例可获取的数据信息
3. 操作流程
下面围绕上述信息展开案例实操，操作流程分为以下几步。
- 了解软件包及数据要求：介绍软件包和数据格式，并对样例数据进行必要的处理。
- 掌握生存分析基础操作：绘制二手车销售生存曲线，使用对数秩检验判断不同生存曲线是否存在显著差异。
- 掌握Cox比例风险模型操作及解读：通过Cox比例风险模型对变量如何影响二手车销售生存曲线展开分析及解读，并对个体维度的生存曲线进行预测。
- Cox比例风险模型应用于最优价格求解：基于前三步的输出结果，结合利润公式求解不同定价下的利润水平，寻找最优价格。
2.3.1 软件包、数据格式和数据读入
1. 软件包
目前，Python中有支持生存分析的软件包Liflines和scikit-survival可供使用，其中Liflines对分析友好，可以支持生存函数的绘制、对数秩检验及Cox模型拟合，本节将通过Python代码实现基于Liflines包的生存分析。
2. 数据格式
生存分析的数据格式以研究对象为单位，每行包括生存时间（删失前观测生存时间），是否发生终点事件的标记及生存概率影响因子。在二手车定价案例中，原始数据以每一个待出售车辆为单位，每一行代表一个车辆的数据信息，信息可以分为以下两类。
- 车辆出售情况：包括车辆是否发生终点事件（是否售出）和生存时间（在售时长）。
- 车辆公开信息（影响因子）：包括车辆颜色、行驶里程、车辆照片、用户访问次数和车辆出售价格。
原始数据中各变量的类型及描述如表2-2所示。
表2-2 二手车定价案例数据各变量说明
3. 数据读入及处理
数据读入及处理如代码清单2-1、代码清单2-2所示。
代码清单2-1 引入软件包
##**软件包引入**重要的软件包的介绍在代码备注中 import pandas as pd ##引入基础分析软件包pandas import numpy as np ##引入基础分析软件包numpy import matplotlib.pyplot as plt ##引入绘图软件包 from lifelines import KaplanMeierFitter ##引入生存分析包-KM生存曲线 from lifelines.statistics import logrank_test ##引入生存分析包-logrank检验 from lifelines import NelsonAalenFitter ##引入生存分析包-风险曲线 from lifelines import CoxPHFitter ##引入生存分析包-Cox模型
代码清单2-2 读入原始数据
##**数据读入**
data_survival = pd.read_csv('input_survival_v2.csv', sep=',', encoding = 'GBK', index_col=False)
data_survival = data_survival[(data_survival['Publish_period']<=90)].reset_index(drop=True)
data_survival = data_survival.drop(['Departure_Date','End_Date','Car_id'],axis=1) ##删除不需要的字段
data_survival.head()
原始数据展示如图2-7所示。
图2-7 二手车定价案例原始数据格式
其中，Color为字符串，属于离散变量，需要对该字段进行处理，如代码清单2-3所示。
代码清单2-3 对离散变量进行处理
##**对离散变量进行处理** df = pd.get_dummies(data_survival, drop_first=True, columns=['Color']) df.head()
处理后的数据如图2-8所示。
图2-8 二手车定价案例处理后的数据格式
2.3.2 绘制二手车销售生存曲线及差异对比
1. 绘制整体生存曲线
数据处理完毕后，我们绘制KM曲线。首先基于全量数据绘制二手车销售情况的生存曲线，如代码清单2-4所示。
代码清单2-4 绘制二手车销售情况随时间变化的整体生存曲线
##**绘制二手车销售情况随时间变化的整体生存曲线
kmf = KaplanMeierFitter()
T = df["Publish_period"] ##T代表生存时长
E = df["is_sold"] ##E代表关注事件（终点事件）
kmf.fit(T, event_observed=E,label='Survival Curve')
#kmf.plot(show_censors=True,ci_show=False) ##绘图，展示删失数据
kmf.plot()
plt.xlim()
plt.ylabel("est. probability of survival $\hat{S}(t)$")
plt.xlabel("day $t$")
绘制结果如图2-9所示。
图2-9 二手车销售生存曲线
如图2-9可知，50%生存概率对应天数在20天左右。
2. 绘制不同访问次数车辆的生存曲线并对比
前面绘制的是整体生存曲线，但在分析中，更需要对不同生存曲线进行绘制及对比，接下来我们绘制两条生存曲线，实现不同访问次数的车辆销售生存曲线对比，如代码清单2-5所示。
代码清单2-5 绘制不同访问次数的二手车销售情况随时间变化的生存曲线
##绘制不同访问次数的二手车销售情况随时间变化的生存曲线
avg_inquires=np.mean(df['N_Inquires'])
df_less_inquires=df[(df['N_Inquires']<avg_inquires)]
df_more_inquires=df[(df['N_Inquires']>avg_inquires)]
ax = plt.subplot(111)
kmf=KaplanMeierFitter()
kmf.fit(df_less_inquires['Publish_period'],
event_observed=df_less_inquires['is_sold'],
label="less_inquires")
ax = kmf.plot(ax=ax)
kmf.fit(df_more_inquires['Publish_period'],
event_observed=df_more_inquires['is_sold'],
label="more_inquires")
ax = kmf.plot(ax=ax)
plt.title("survival curves of two different type");
绘制结果如图2-10所示。
图2-10 不同访问次数二手车生存曲线
对两条曲线进行观察后发现，访问次数更少的车辆生存曲线更靠左，说明在相同时间下，访问次数少的车辆出售的概率更高，接下来我们通过对数秩检验对曲线差异显著性进行检验，如代码清单2-6所示。
代码清单2-6 验证不同访问次数二手车销售生存曲线是否有显著差异
##**验证不同访问次数的二手车销售生存曲线是否有显著差异
T1 = df_less_inquires["Publish_period"]
E1 = df_less_inquires["is_sold"]
T2 = df_more_inquires["Publish_period"]
E2 = df_more_inquires["is_sold"]
results = logrank_test(T1, T2,
event_observed_A=E1,
event_observed_B=E2)
results.print_summary()
检验结果如图2-11所示。
图2-11 不同访问次数二手车生存曲线差异检验
图2-11所示的结果显示P值小于0.05，可以认为两条曲线存在显著差异，那么我们可以直接得出访问越少，销售速度越快的结论吗？答案是否定的，我们在分析过程中需要警惕一点，那就是相关并不代表因果，出现这一现象有可能是因为我们只考虑访问数量一个因素，而这些访问少的车辆恰好价格更低或者行驶里程更短。因此，对于没有控制其他因素下的单一因素对比需要谨慎对待。
2.3.3 二手车销售生存概率影响因素分析及个体预测
1. 影响因素分析
接下来我们引入Cox比例风险回归模型，对车辆颜色、行驶里程、车辆照片、用户访问次数、车辆出售价格这5个变量进行分析，找出对生存函数有显著影响变量，如代码清单2-7所示。
代码清单2-7 引入Cox模型对特征进行解释
##**引入Cox模型对特征进行解释 cph = CoxPHFitter() cph.fit(df, duration_col='Publish_period', event_col='is_sold',show_progress=True,step_size=0.5) cph.print_summary() cph.plot() ##图形对特征进行解释
模型结果如图2-12所示。
图2-12 模型拟合结果（各特征系数及显著性）
由模型结果可见，5个变量对结果均有显著影响（P值均小于0.05），对于系数的解读如下。
- 连续变量：行驶里程、照片数量、访问次数及价格的系数为负，代表均为保护因素，对销售速度有负向影响。控制其他因素后可以看到，访问次数依然有显著影响（一般情况下需要先对变量的相关性进行检验并剔除高度相关变量，此处在数据采样阶段已进行过相关性检验，因此分析过程中不再重复操作）。
- 分类变量：以黑色为基准，棕色的系数为负，表示保护因素，说明棕色车对销售速度相对黑色车更差，而白色、银色车的销售速度优于黑色，其中白色最优。
模型中Concordance是解释模型整体效果的参数，与一般模型中的AUC含义类似，结果显示模型的Concordance为0.69，表现尚可。
2. 个体生存概率预测
在模型拟合之后，可以实现对个体生存概率的预测，我们采用Cox模型拟合结果对个体进行预测并随机抽取两个样本（车辆ID为24及11166）展示生存概率预测曲线，如代码清单2-8所示。
代码清单2-8 对个体生存函数进行预测
##**对个体生存函数进行预测 X = df.drop(['Publish_period','is_sold'],axis=1) ##筛选特征集合，剔除销售时间与事件结局数据 surv_hat=cph.predict_survival_function(X) ## 抽取任意两个样本，观察预测结果 surv_hat[24].plot(label='24') surv_hat[11166].plot(label='11166') plt.legend()
两个样本的预测生存曲线绘制如图2-13所示。
图2-13 两个样本的预测生存曲线
以50%生存概率对应生存时间（中位生存时间）为预测生存时间，图2-13所示二者的预测生存时间分别为62天和41天。
2.3.4 基于Cox风险比例模型的最优价格求解
最优价格分析
在操作之前，我们首先来看一下最优价格的分析流程，如图2-14所示。
图2-14 二手车最优价格分析流程
基于分析流程，我们在目前发布的最低价与最高价之间创建等差序列，选定10个价格作为预选价格，并在其中选出利润最大时对应价格，如代码清单2-9所示。
代码清单2-9 基于预测生存函数寻找最优价格策略
##**基于预测生存函数寻找最优价格策略
## 第1步：创建获取预测在售天数函数
def predict_day(surv_hat):
days = np.zeros(surv_hat.shape[1])
prob = np.zeros(surv_hat.shape[1])
j = surv_hat.shape[1]
for i in range(1,surv_hat.shape[1]):
prob[i-1] = surv_hat[surv_hat[i-1] >= 0.5][i-1].min()
prob[j-1] = surv_hat[surv_hat[j-1] >= 0.5][j-1].min()
days[i-1] = surv_hat[surv_hat[i-1] == prob[i-1]].index.values.min()
days[j-1] = surv_hat[surv_hat[j-1] == prob[j-1]].index.values.min()
##以预测半衰期作为预测在售天数
return prob,days
## 第2步：创建90天是否卖出函数
def is_sold(data):
y = np.zeros(data.shape[0])
for i in range(1,data.shape[0]):
if data[i-1]>=0.6:
y[i-1]=0
else:
y[i-1]=1
return y
## 第3步：创建利润函数
def profit(data, predict_days,sold_tag):
d = list(predict_days)
y = list(sold_tag)
revenue = np.sum(data['Price']*y)
cost = np.sum(1000*d)
profit = revenue - cost
return profit
## 第4步：计算不同价格下的利润
min_price = df['Price'].values.min()
max_price = df['Price'].values.max()
sp_price = np.linspace(min_price, max_price, 10) ##选定十个预选价格
X = df.drop(['Publish_period','is_sold'],axis=1)
profit_list = []
price_list = list(sp_price)
for p in price_list:
X['Price'] = p
surv_hat = cph.predict_survival_function(X)
prob_result,days_result = predict_day(surv_hat)
sold_result = is_sold(prob_result)
profit_result = profit(X,days_result,sold_result)
profit_list.append(profit_result)
profit_res = pd.DataFrame({'price': price_list, 'profit': profit_list})
得出结果如图2-15所示。
图2-15 价格与利润对应关系图
由此可见，在目前价格水平下，收益处于单调递增阶段，即售价越高利润越高，目前最优价格处于价格上限21万。
至此，整个分析实操就结束了，我们通过生存分析实现了最优价格的筛选，在固定分析框架后，我们将代码整理打包，就能达到自动化定价的目的。需要说明的是，这一分析数据仅为样例数据，不代表真实情况，因此结果与现实世界相比的合理性不是最重要的，掌握分析方法并能够有效使用实操才是本节学习的关键。