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。

​