Python 数据分析实例——生存分析

什么是生存？生存的意义很广泛，可以指人或动物的存活（相对于死亡），可以是患者的病情正处于缓解状态（相对于再次复发或恶化），还可以是某个系统或产品正常工作（相对于失效或故障），甚至可以是客户的流失与否，等等。在生存分析中，研究的主要对象是寿命超过某一时间的概率。还可以描述其他一些事情发生的概率，例如产品的失效、出狱犯人第一次犯罪、失业人员第一次找到工作等。在某些领域的分析中，经常用追踪的方式来研究事物的发展规律，比如研究某种药物的疗效、手术后的存活时间、某件机器的使用寿命等。

1．概念

生存分析是对一个或多个非负随机变量进行统计推断，研究生存现象和响应时间数据及其统计规律的一门学科。

生存分析是既考虑结果又考虑生存时间，并可充分利用截尾数据所提供的不完全信息对生存时间的分布特征进行描述，对影响生存时间的主要因素进行分析。

2．生存分析研究的内容

（1）描述生存过程

研究生存时间的分布特点，估计生存率及平均存活时间，绘制生存曲线等，根据生存时间的长短可以估算出各个时点的生存率，并根据生存率来估计中位生存时间，也可以根据生存曲线分析其生存特点，一般使用Kaplan-Meier法和寿命表法。

（2）比较生存过程

可通过生存率及其标准误对各样本的生存率进行比较，以探讨各组间的生存过程是否存在差异，一般使用Log-rank检验和Breslow检验。

（3）分析危险因素

通过生存分析模型来探讨影响生存时间和终点事件的保护因素和不利因素、因素作用的大小及方向、相对危险度的大小，基本使用Cox回归模型。

（4）建立数学模型

建立最终的数学模型通过Cox回归模型完成。

3．生存分析对资料的基本要求

·　样本由随机抽样方法获得，要有一定的数量，死亡例数和比例不能太少。

·　完整数据所占的比例不能太少，即截尾值不宜太多。

·　截尾值出现的原因无偏性，为防止偏性，经常对被截尾的研究对象的年龄、职业、地区、病情轻重等情况进行分析。

·　生存时间尽可能精确。

·　缺项要尽量补齐。

4．生存资料的共同特点

·　蕴含结局和时间两个方面的信息。

·　结局为两分类事件。

·　一般通过随访收集得到，随访观察往往是从某统一时间点（如入院或实施手术等某种处理措施后）开始，观察到某规定时间点截止。

·　常因失访等原因造成研究对象的生存时间数据不完整，分布类型复杂，不能简单地套用以前的方法。

5．一些相关的基本概念

·　起始事件：反映研究对象开始生存过程的起始特征事件，如研究某一治疗对病人生存的影响的起始时间是“开始接受该治疗”。

·　终点事件（死亡事件）：出现研究者所关心的特定结局，如“病人因该疾病死亡”。

·　观察时间：从研究开始观察到研究观察结束的时间。由于研究时长无法无限延伸下去，因此研究一定会在某个特定时刻截止，而研究截止时，所有观察对象并不一定全都出现终点事件。换言之，有的研究对象在观察结束之前出现终点事件，有的直到观察结束时也没有出现终点事件，还有一些特例中途因为某些原因（如失访、意外死亡等）被迫提前结束了观察研究。

·　生存时间：观察到的存活时间，用符号t表示。

·　完全数据：从观察起点到死亡事件所经历的时间，生存时间是完整的。

·　截尾数据（删失值）：观察时间不是由于终点事件而结束的，而是由于失访、死于非研究因素、观察结束而对象仍存活3种原因结束的。常在截尾数据的右上角放一个“+”表示其实该对象可能活得更久。

6．生存分析的主要方法

（1）非参数法

这类方法的特点是，无论分布形式如何，只根据样本的顺序统计量对生存率进行估计。对于两个及多个生存率进行比较，其无效假设只是假定两组或多组生存时间分布相同，而不对其具体的分布形式和参数进行推断。log-rank乘法极限法和寿命表法都是非参数法。

（2）参数法

特点是假定生存时间服从特定的参数分布，然后根据已知的分布特点对生存时间进行分析，如指数分布法、Weibull分布法、对数正态回归分布法和Logistic回归法。

（3）半参数法

Cox比例风险回归模型就是半参数法，具体介绍它时再说明为什么叫半参数法。

7．生存分析公式模型

（1）概率密度函数f(t)

表示每时刻死亡的概率：

分布函数F(t)：

（2）生存函数S(t)

生存函数S(t)：P(T≥t)个体生存时间大于t的概率。

（3）危险函数h(t)

（4）累计危险函数H(t)

8．生存分析的目的

估计：根据样本生存资料估计总体生存率及其他有关指标（如中位生存期等），估计不同时间的生存率、生存曲线以及中位生存期等。

比较：对不同处理组的生存率进行比较，以了解哪种方案较优。

影响因素分析：目的是为了探索和了解影响生存时间长短的因素，或平衡某些因素的影响后，研究某个或某些因素对生存率的影响。

预测：具有不同因素水平的个体生存预测。

9．模拟实验代码实现

【例1】

问题描述：泰坦尼克号的沉没是历史上非常有名的沉船之一。1912年4月15日，泰坦尼克号在处女航时与冰山相撞沉没，2224名乘客和船员中有1502人遇难。这一耸人听闻的悲剧震惊了国际社会，并导致制定了更好的船舶安全法规。船难造成如此巨大的人员伤亡的原因之一是船上没有足够的救生艇供乘客和船员使用。虽然在沉船事件中幸存下来有运气因素，但有些人比其他人更有可能存活下来，比如妇女、儿童和上层阶级。

模拟实验分析中需要预测哪一类人更有可能存活下来。需要用机器学习的工具去预测哪些乘客在这次灾难中幸存。本例采用Anaconda 3下的Jupyter Notebook开发环境。

内容：

·　导入必要的库。

·　阅读并研究数据。

·　数据分析。

·　数据可视化。

·　清理数据。

·　选择最佳模式。

·　创建提交文件。

#导入必要的库
In[1]:
#data analysis libraries
import numpy as np
import pandas as pd

#visualization libraries
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

#ignore warnings
import warnings
warnings.filterwarnings('ignore')
#读取并探索数据
In[2]:
#import train and test CSV files
train = pd.read_csv("D:\Anaconda3\workspace/train.csv")
test = pd.read_csv("D:\Anaconda3\workspace/test.csv")

#take a look at the training data
train.describe(include="all")           

输出：

数据分析，将考虑数据集中的特性以及它们的完整性：

In[3]:
#get a list of the features within the dataset
print(train.columns)           

输出训练数据的列属性：

参阅数据集示例了解变量：

训练数据集的摘要：

检查其他不可用的值：

画一个以性别做对比条件的生存几率图，这是因为全球都说lady first，女士优先：

性别生存率对比图如图1所示。

图1　性别生存率对比图

按阶层生存率对比，坐船有等级之分，像高铁、飞机都有。这个属性会对生产率产生影响。因为一般有钱人、权贵才会住头等舱。

按阶层生存率对比图如图2所示。

图2　按阶层生存率对比图

兄弟姐妹生存率对比，一些人是和兄弟姐妹一起上船的。这个会有影响，因为有可能因为救兄弟姐妹而导致自己没有上救生船。

兄弟姐妹生存率对比图如图3所示。

图3　兄弟姐妹生存率对比图

父母和小孩生存率对比，有些人是带着父母、小孩上船的，有可能因为要救父母、小孩耽误上救生船。​

In[10]:
#draw a bar plot for Parch vs. survival
sns.barplot(x="Parch", y="Survived", data=train)
plt.show()           

父母和小孩生存率对比图如图4所示。

图4　父母和小孩生存率对比图

按年龄逻辑分类：

按年龄逻辑分类图如图5所示。

图5　按年龄逻辑分类图

计算CabinBool与存活的百分比：

CabinBool与存活的百分比图如图6所示。

图6　CabinBool与存活的百分比

填充登录地点缺失的值：

创建两个数据集（train & test）组合：

用常用名称替换标题：

将每个标题组映射到一个数值：

为每一个标题填充缺失的年龄与模式年龄组：

将每个年龄值映射到一个数值：​

In[18]:
#map each Age value to a numerical value
age_mapping = {'Baby': 1, 'Child': 2, 'Teenager': 3, 'Student': 4, 'Young Adult':
5, 'Adult': 6, 'Senior': 7}
train['AgeGroup'] = train['AgeGroup'].map(age_mapping)
test['AgeGroup'] = test['AgeGroup'].map(age_mapping)

train.head()

#dropping the Age feature for now, might change
train = train.drop(['Age'], axis = 1)
test = test.drop(['Age'], axis = 1)           

删除name特性，因为它不包含更多有用的信息：

In[19]:
#drop the name feature since it contains no more useful information.
train = train.drop(['Name'], axis = 1)
test = test.drop(['Name'], axis = 1)           

将每个性别值映射到一个数值：

将每个登船值映射到一个数值：

根据舱级的平均票价，在test集中填写缺失的票价值：

检测train集的数据：

检测test集的数据：

生存预测模型分析：​

In[25]:
from sklearn.model_selection import train_test_split

predictors = train.drop(['Survived', 'PassengerId'], axis=1)
target = train["Survived"]
x_train, x_val, y_train, y_val = train_test_split(predictors, target, test_size
= 0.22, random_state = 0)           

高斯朴素贝叶斯：

In[26]:
# Gaussian Naive Bayes
from sklearn.naive_bayes import GaussianNB
from sklearn.metrics import accuracy_score

gaussian = GaussianNB()
gaussian.fit(x_train, y_train)
y_pred = gaussian.predict(x_val)
acc_gaussian = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_gaussian)
Out[26]:
78.68           

逻辑回归：​

In[27]:
# Logistic Regression
from sklearn.linear_model import LogisticRegression

logreg = LogisticRegression()
logreg.fit(x_train, y_train)
y_pred = logreg.predict(x_val)
acc_logreg = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_logreg)

Out[27]:
79.19           

支持向量机：​

In[28]:
# Support Vector Machines
from sklearn.svm import SVC

svc = SVC()
svc.fit(x_train, y_train)
y_pred = svc.predict(x_val)
acc_svc = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_svc)

Out[28]:
82.74           

线性支持向量机：​

In[29]:
# Linear SVC
from sklearn.svm import LinearSVC

linear_svc = LinearSVC()
linear_svc.fit(x_train, y_train)
y_pred = linear_svc.predict(x_val)
acc_linear_svc = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_linear_svc)

Out[29]:
78.17           

感知机：​

In[30]:
# Perceptron
from sklearn.linear_model import Perceptron

perceptron = Perceptron()
perceptron.fit(x_train, y_train)
y_pred = perceptron.predict(x_val)
acc_perceptron = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_perceptron)

Out[30]:
79.19           

决策树：​

In[31]:
#Decision Tree
from sklearn.tree import DecisionTreeClassifier

decisiontree = DecisionTreeClassifier()
decisiontree.fit(x_train, y_train)
y_pred = decisiontree.predict(x_val)
acc_decisiontree = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_decisiontree)

Out[31]:
81.22           

随机森林：​

In[32]:
# Random Forest
from sklearn.ensemble import RandomForestClassifier

randomforest = RandomForestClassifier()
randomforest.fit(x_train, y_train)
y_pred = randomforest.predict(x_val)
acc_randomforest = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_randomforest)

Out[32]:
82.74           

K-紧邻：​

In[33]:
# KNN or k-Nearest Neighbors
from sklearn.neighbors import KNeighborsClassifier

knn = KNeighborsClassifier()
knn.fit(x_train, y_train)
y_pred = knn.predict(x_val)
acc_knn = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_knn)

Out[33]:
77.66           

随机梯度下降：​

In[34]:
# Stochastic Gradient Descent
from sklearn.linear_model import SGDClassifier

sgd = SGDClassifier()
sgd.fit(x_train, y_train)
y_pred = sgd.predict(x_val)
acc_sgd = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_sgd)

Out[34]:
80.71           

梯度提升分类器：​

In[35]:
# Gradient Boosting Classifier
from sklearn.ensemble import GradientBoostingClassifier

gbk = GradientBoostingClassifier()
gbk.fit(x_train, y_train)
y_pred = gbk.predict(x_val)
acc_gbk = round(accuracy_score(y_pred, y_val) * 100, 2)
print(acc_gbk)

Out[35]:
84.77           

预测模型算法得分统计：

创建生存预测文件，将id设置为乘客id并预测生存：​

In[37]:
#set ids as PassengerId and predict survival
ids = test['PassengerId']
predictions = gbk.predict(test.drop('PassengerId', axis=1))

#set the output as a dataframe and convert to csv file named submission.csv
output = pd.DataFrame({ 'PassengerId' : ids, 'Survived': predictions })
output.to_csv('submission.csv', index=False)

Out[37]:           

预测数据片段：

本次实验的数据集titanic.csv、test.csv、train.csv从Kaggle（https://www.kaggle.com/）平台下载，软件平台采用Jupyter Notebook。

​