数据分析：使用PYTHON做生存分析（生存曲线、COX回归）

逻辑回归不用多说了，到处都在用。网上的生存分析大多是使用R或者SPSS。本文主要使用python进行生存分析。

生存分析的理论很多，本文不涉及，推荐写的很好的生存分析系列文章生存分析简介。写的非常清晰，本人也是在长时间不使用后，看着这个系列的文章才慢慢回想起来。

在python中，生存分析的包是lifelines，可以利用lifelines进行累计生存曲线的绘制、Log Rank test、Cox回归等。

一、绘制S(T)

先使用自带的数据集

from lifelines.datasets import load_waltons
from lifelines import KaplanMeierFitter
from lifelines.utils import median_survival_times

df = load_waltons()
print(df.head(),'n')
print(df['T'].min(), df['T'].max(),'n')
print(df['E'].value_counts(),'n')
print(df['group'].value_counts(),'n')

可以看到数据有三列，其中T代表min(T, C)，其中T为死亡时间，C为观测截止时间。E代表是否观察到“死亡”，1代表观测到了，0代表未观测到，即生存分析中的**“删失”**数据，共7个。 group代表是否存在病毒， miR-137代表存在病毒，control代表为用药，求中存在miR-137病毒人数34人，不存在129人。

利用此数据取拟合拟生存分析中的Kaplan Meier模型（专用于估计生存函数的模型），并绘制全体人群的生存曲线。

kmf = KaplanMeierFitter()
kmf.fit(df['T'], event_observed=df['E'])

kmf.plot_survival_function()

median_ = kmf.median_survival_time_
median_confidence_interval_ = median_survival_times(kmf.confidence_interval_)
print(median_confidence_interval_)

图中蓝色实线为生存曲线，浅蓝色带代表了95%置信区间。随着时间增加，存活概率S(t)越来越小，这是一定的，同时S(t)=0.5时，t的95%置信区间为[53, 58]。这并不是我们关注的重点，我们真正要关注的实验组（存在病毒）和对照组（未存在病毒）的生存曲线差异。因此我们要按照group等于“miR-137”和“control”分组，分别观察对应的生存曲线:

groups = df['group']
ix = (groups == 'miR-137')

kmf.fit(df['T'][ix], df['E'][ix], label='miR-137')
ax = kmf.plot()
treatment_median_confidence_interval_ = median_survival_times(kmf.confidence_interval_)
print("带有miR-137病毒存活50%对应的存活时间95%置信区间：'n'", treatment_median_confidence_interval_, 'n')

kmf.fit(df['T'][~ix], df['E'][~ix], label='control')
#共享一个画布
ax = kmf.plot(ax=ax)

control_median_confidence_interval_ = median_survival_times(kmf.confidence_interval_)
print("未带有miR-137病毒存活50%对应的存活时间95%置信区间：'n'", control_median_confidence_interval_)

可以看到，带有miR-137病毒的生存曲线在control组下方。说明其平均存活时间明显小于control组。同时带有miR-137病毒存活50%对应的存活时间95%置信区间为[19,29]，对应的control组为[56,60]。差异较大，这个方法可以应用在分析用户流失等场景，比如我们对一组人群实行了一些防止流行活动，我们可以通过此种方式分析我们活动是否有效。

二、COX回归

通常存活时间与多种因素都存在关联，因此我们的面临的数据是多维的。下面使用一个更复杂的数据集。

from lifelines.datasets import load_regression_dataset
from lifelines import CoxPHFitter

regression_dataset = load_regression_dataset()

print(regression_dataset.head())
print(regression_dataset['E'].value_counts())

其中T代表min(T, C)，其中T为死亡时间，C为观测截止时间。E代表是否观察到“死亡”，1代表观测到了，0代表未观测到，即生存分析中的**“删失”**数据，删失数据共11个。var1,var2,var3代表了我们关系的变量，可以是是否为实验组的虚拟变量，可以是一个用户的渠道路径，也可以是用户自身的属性

我们利用此数据进行Cox回归

cph = CoxPHFitter()
cph.fit(regression_dataset, 'T', event_col='E')
cph.print_summary()

从结果来看，我们认为var1和var2在5%的显著性水平下是显著的。我们看到结果中第二列为exp(coef)，我们可以利用这个数字带入PAF的计算公式进行归因分析。