开篇语
生存分析在医学研究中占有很大的比例，而且进行生存分析时，多用R语言、SPSS等工具进行生存分析，用python进行生存分析不多。因为发现一个python版的生存分析工具—lifelines ，这个库已经提供比较完善的生存分析相关的工具。自己又最近学习生存分析，然后结合lifelines开始编写这个项目。写代码的同时，也对一些生存分析中概念性的名词，根据自己的理解一起展示出来。因为是边学边写，有错误的地方请指正 。
#安装生存分析用的python库----lifelines
#lifelines相关文档地址https://lifelines.readthedocs.io/en/latest/
!pip install lifelines
!pip install numpy==1.19.2
!pip install pyzmq==18.1.1
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
1.什么是生存分析
在医学中，生存分析是非常常见和重要的分析方法之一,是一种将事件的结果和出现这种结果所经历的时间结合起来一起分析的统计方法。
案例一：
为了研究什么因素会影响直肠癌患者生存时间。会收集患者的资料。例如有年龄、性别、婚姻状态、TMN分期、T分期、肿瘤直径、肿瘤分化程度、治疗手段等。还有会收集患者在一定时间内发生死亡事件的时间，即生存时间。采用Kaplan–Meier 曲线单变量进行分析，那个因素是影响直肠癌患者的生存时间，或者采用Cox 比例风险回归进行多因素分析。并对新患有直肠癌患者进行生存时间预测。
案例二：为了对某种抗癌药物做临床试验，将癌症患者随机分成两组，一组服用该抗癌药物，另一组服用对照药。并记录两组患者开始服药直到死亡的生存时间。并分析两组中死亡事件发生与时间关系是否存在差异，来判断药物是否有效。
通过两个案例可以看出生存分析就是研究在不同条件下，事件发生与时间的关系 ，不仅考虑事件是否出现，还有考虑事件出现的时间长短。
1.1生存分析的数据是怎样的
生存分析的数据和分类数据有点类似。但是生存分析多了一列关于时间数据，这列时间数据是从特定事件发生了，记录整个观察期间研究对象多久才发生失效事件(例如死亡事件）。生存分析就是研究生存时间和结局与众多影响因素间的关系。
例如下面表格记录了6位癌症患者在2年观测期间中的存活时间和感兴趣的影响因素。这里的失效事件是死亡事件，是我们关心的。生存时间对应表格的生存时间。研究的影响因素有：年龄、性别、是否吸烟、肿瘤大小。
在表格中有个患者生存时间只有21个月，结局却是存活，可能由于是患者中途退出观察或者换了电话号码，不能联系患者，无法继续进行研究，这种数据属于删失数据，可见在现实生活中生存分析数据通常会含有删失数据。再例如在表格中有个患者生存时间是19个月，结局是死亡，这种数据属于完整数据。不过这些删失数据在进行生存分析时，也是会一起进行分析。
| 年龄 | 性别 | 是否吸烟 | 肿瘤大小(mm) | 生存时间(月) | 事件 |
|---|---|---|---|---|---|
| 55 | 男 | 是 | 17 | 23 | 存活 |
| 65 | 女 | 是 | 20 | 24 | 存活 |
| 36 | 女 | 否 | 16 | 19 | 死亡 |
| 38 | 男 | 是 | 25 | 21 | 存活 |
| 69 | 男 | 否 | 16 | 23 | 存活 |
| 70 | 男 | 否 | 28 | 20 | 死亡 |
1.2 生存分析中常见专业术语
生存分析有一些常见的专业术语。例如什么是生存率，什么是生存曲线等。
1.3 生存分析的主要目的
- 估计研究对象的生存率
- 比较2不同组(影响因素)的生存率
- 分析影响研究对象生存期长短的因素有哪些和贡献度
1.4 生存分析的主要方法
- 寿命表法和Kaplan-Meier法（属于统计描述)
- Log-rank检验（属于统计推断)
- Cox比例风险回归模型（属于统计推断)
【1】寿命表和KM法，都属于非参数法，都是用来对生存数据进行统计描述，例如估计各个时间段的生存函数并绘制生存曲线。寿命法多用于，样本庞大，随访时间跨间较大的数据，例如间隔1年才对研究对象随访一次，KM法和寿命法计算生存函数都是用到乘积极限法，但在实际情况中，寿命表分析不是太常用，多用KM法进行分析。KM法可以计算中位生存时间，对数据进行分组分别绘制生存曲线进行单因素分析（单因素分析只能对分类变量进行分析，如果对连续变量分析，需要对连续变量进行分组，变成分类)。
【2】Log-rank检验，属于非参数检验，用于比较两组或多组生存曲线或生存时间是否相同，可以比较KM法中分组后绘制生存曲线。
【3】Cox比例风险回归模型 ，属于半参数法，可以估计生存函数，可以比较两组或多组生存函数，可以单因素分析也可以多因素分析，单因素分析包括分类变量和连续变量，可以建立生存时间与影响因素之间的关系模型，计算影响因素的分险比。注意Cox回归要求满足比例风险假定，就是影响因素是不随时间变化而变化。
1.5 生存分析的评价指标
评价模型的好坏一方面有模型的拟合优度，例如R方，AIC等。一方面有模型预测精度，例如均方误差，AUC等。在临床上更关注的模型预测精度，在生存分析建模中常用的模型预测精度评价指标有（这里只讲讲C-index指数。。。。）:
C-index(一致性指数)
c-index反映的是模型预测结果与实际情况的一致程度。
计算方法是把所以研究对象随机两两组成一对，模型预测的生存时间更长的对象，他的实际生存时间也更长，或者预测的生存概率的高的生存时间，实际生存时间也更长。
2.数据处理和分析
2.1 查看本项目中乳腺癌数据的基线资料表
这个乳腺癌数据集研究的失效事件是"Status"中的"Dead",生存时间是"Survival_Months"单位是月。
这个乳腺癌数据集中具有的特征如下：
- Age(年龄)
- Race(种族）
- Marital_Status(婚姻状态)
- T_Stage(原发肿瘤的大小)
- N_Stage(癌细胞扩散的淋巴结情况)
- Sixth_Stage
- Grade
- A_Stage
- Tumor_Size(肿瘤大小)
- Estrogen_Status(雌激素受体是否阳性)
- Progesterone_Status(孕激素受体是否阳性)
- Regional_Node_Examined
- Reginol_Node_Positive
- Survival_Months（生存时间，单位:月)
通过基线资料表中根据“Alive”和“Dead”分成两组，可以看出每个特征的人数分布和比例。
2.2 读取、处理数据
通过pandas读取csv文件获取数据。通过head()可以看到某些列是由字符串组成的，例如“Race”列中每个值由“Black”、“White”等值组成，所以进行生存分析时需要编码成数字0、1、2等。
"""
读取数据并处理分类变量数据，不然使用cox模型会报错
对列值时字符串进行编码，将一种符号编成一种数字码
因为有些特征不具有大小之分，例如Race，那个是0，那个是1，都没问题，可以通过pd.factorize进行简单处理。
例如T_Stage明显有渐进的，可以处理成T1是0，T2是1。
例如['T2', 'T1', 'T3', 'T4'] ->[1,0,2,3]
这里我们关心的失效事件是Dead，删失数据Alive，所以“Status”需要处理成Alive为0，Dead为1
"""
df = pd.read_csv('SEER Breast Cancer Dataset .csv')
for name in df.columns.values:
if df[name].dtype == object:
print(f'Column name:{name}')
print(f"Value:{df[name].unique().tolist()}\n")
df['Race'] = pd.factorize(df['Race'])[0].astype(np.uint16)
df['Marital_Status'] = pd.factorize(df['Marital_Status'])[0].astype(np.uint16)
df['A_Stage'] = pd.factorize(df['A_Stage'])[0].astype(np.uint16)
df['T_Stage'] = df['T_Stage'].map({'T1':1, 'T2':2,'T3':3,'T4':4})
df['N_Stage'] = df['N_Stage'].map({'N1':0, 'N2':1,'N3':2})
df['Sixth_Stage'] = df['Sixth_Stage'].map({'IIA':0, 'IIB':1,'IIIA':2,'IIIB':3,'IIIC':4})
df['Grade'] = df['Grade'].map({'Well differentiated; Grade I':0,
'Moderately differentiated; Grade II':1,
'Poorly differentiated; Grade III':2,
'Undifferentiated; anaplastic; Grade IV':3})
df['Estrogen_Status'] = df['Estrogen_Status'].map({'Positive':0, 'Negative':1})
df['Progesterone_Status'] = df['Progesterone_Status'].map({'Positive':0, 'Negative':1})
df['Status'] = df['Status'].map({'Alive':0, 'Dead':1})
df.head()
Column name:Race
Value:['Other (American Indian/AK Native, Asian/Pacific Islander)', 'White', 'Black']
Column name:Marital_Status
Value:['Married (including common law)', 'Divorced', 'Single (never married)', 'Widowed', 'Separated']
Column name:T_Stage
Value:['T2', 'T1', 'T3', 'T4']
Column name:N_Stage
Value:['N3', 'N2', 'N1']
Column name:Sixth_Stage
Value:['IIIC', 'IIIA', 'IIB', 'IIA', 'IIIB']
Column name:Grade
Value:['Moderately differentiated; Grade II', 'Poorly differentiated; Grade III', 'Well differentiated; Grade I', 'Undifferentiated; anaplastic; Grade IV']
Column name:A_Stage
Value:['Regional', 'Distant']
Column name:Estrogen_Status
Value:['Positive', 'Negative']
Column name:Progesterone_Status
Value:['Positive', 'Negative']
Column name:Status
Value:['Alive', 'Dead']
| Age | Race | Marital_Status | T_Stage | N_Stage | Sixth_Stage | Grade | A_Stage | Tumor_Size | Estrogen_Status | Progesterone_Status | Regional_Node_Examined | Reginol_Node_Positive | Survival_Months | Status | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 43 | 0 | 0 | 2 | 2 | 4 | 1 | 0 | 40 | 0 | 0 | 19 | 11 | 1 | 0 |
| 1 | 47 | 0 | 0 | 2 | 1 | 2 | 1 | 0 | 45 | 0 | 0 | 25 | 9 | 2 | 0 |
| 2 | 67 | 1 | 0 | 2 | 0 | 1 | 2 | 0 | 25 | 0 | 0 | 4 | 1 | 2 | 1 |
| 3 | 46 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 19 | 0 | 0 | 26 | 1 | 2 | 1 |
| 4 | 63 | 1 | 0 | 2 | 1 | 2 | 1 | 0 | 35 | 0 | 0 | 21 | 5 | 3 | 1 |
2.3 简单数据分析
对数据绘制柱状图、箱线图等，展示某些数据与终末事件的关系
"""
一共有4024个研究对象，发生删失3408例，发生死亡事件616例
"""
print(len(df))
print(df['Status'].value_counts(),'\n')
4024
0 3408
1 616
Name: Status, dtype: int64
#不同T(肿瘤大小)分期，删失人数和死亡人数的柱状分布图
sns.countplot(x="T_Stage",hue="Status",data=df,)
<AxesSubplot:xlabel='T_Stage', ylabel='count'>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mJ86lnlJ-1645703648858)(output_15_1.png)]
#不同淋巴结转移情况，删失人数和死亡人数的柱状分布图
sns.countplot(x="N_Stage",hue="Status",data=df,)
<AxesSubplot:xlabel='N_Stage', ylabel='count'>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MVpYq555-1645703648859)(output_16_1.png)]
"""
绘制发生删失组和发生死亡组的年龄箱线图
发生死亡事件组的年龄中位数比删失组高一点
"""
ax = df.boxplot(column='Age',by="Status",figsize=(6,6))
ax.set_xticklabels(["censor","death"])
[Text(1, 0, 'censor'), Text(2, 0, 'death')]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-52x0zxw0-1645703648860)(output_17_1.png)]
2.4 绘制生命周期图
在数据集中随机挑选30例患者，绘制生命周期图。假设观测时间定位96个月，在96个月中，有5位患者发生失效事件(死亡),其余患者大多数还没到96个月就失去了联系/无法随访到患者(这种数据属于删失数据）。例如图中1347号患者，在随访到第62个月，患者就失去了联系，并不知道患者后面的真实存活时间。
通过生命周期图，可以看到出现很多数据都是删失数据，无论不结合删失数据，单独计算发生失效事件的患者的平均存活时间，还是结合删失数据一起一起计算患者的平均存活时间，这样计算都是要低估了真实的存活时间，可以认识到了右删失数据对生存时间估计的影响很大，因此生存分析就是解决此类问题。
from lifelines.plotting import plot_lifetimes
df_sample = df.sample(n=30,random_state=10)
df_sample['Status'] = np.where(df_sample['Status']== 1,True,False)
plt.figure(figsize=[8,6])
ax = plot_lifetimes(df_sample['Survival_Months'], event_observed=df_sample['Status'],sort_by_duration=False)
ax.vlines(96, 0, 30, lw=2, linestyles='--',colors='black')
plt.xlabel("Months")
3.开始—生存分析
3.1 寿命表分析法
寿命表将事件分成时间段（月、年），按时间段统计事件的发生情况，从而计算生存率。一般用于样本较大的，时间跨度区间较大的（例如1年、2年）的随访资料。因为某些原因，对于研究对象并不是每天进行联系，而是在某个时间区间，例如一年后获取对象的资料。因此并不是准确记录结局的发生时间，只知道某个区间中发生了终点事件。这种情况可以绘制寿命表计算生存率进行分析。
假设有寿命表(类似频数表)如下：
| 确诊年数 | 死亡人数 | 删失人数 | 初期人数 |
|---|---|---|---|
| 0~1 | 1 | 0 | 100 |
| 1~2 | 1 | 0 | 99 |
| 2~3 | 3 | 1 | 98 |
| 3~4 | 1 | 1 | 94 |
| 4~5 | 2 | 0 | 92 |
| 5~ | 3 | 87 | 90 |
可以通过寿命表计算死亡概率和生存概率和生存率
死亡概率= 死亡人数/(初期人数-删失人数/2)
生存概率=1-死亡概率
生存率= 生存概率1 X 生存概率2 X 生存概率3。。
例如2-3区间的死亡概率 = 3/(98-1/2) = 0.0307。
生存概率= 1-0.0307 =0.969
生存率= 0.99 * 0.9899 * 0.969=0.949
就可以通过计算出来的生存率绘制生存曲线
| 确诊年数 | 死亡人数 | 删失人数 | 初期人数 | 死亡概率 | 生存概率 | 生存率 |
|---|---|---|---|---|---|---|
| 0~1 | 1 | 0 | 100 | 0.01 | 0.99 | 0.99 |
| 1~2 | 1 | 0 | 99 | 0.0101 | 0.9899 | 0.9799 |
| 2~3 | 3 | 1 | 98 | 0.0307 | 0.969 | 0.949 |
| 3~4 | 1 | 1 | 94 | 0.017 | 0.989 | 0.938 |
| 4~5 | 2 | 0 | 92 | 0.022 | 0.978 | 0.918 |
| 5~ | 3 | 87 | 90 | 0.0645 | 0.9355 | 0.859 |
进行单因素分析时，可以根据影响因素进行分组，分别制作寿命表，进行计算各个指标并绘制生存曲线进行比较。
3.2 对乳腺癌数据集生成寿命表
计算生存率并绘制生存曲线
from lifelines.utils import survival_table_from_events
#原本数据是生存时间是以月位单位的，现在特意计算成以年位单位，然后制作寿命表
def months_to_year(x):
if x % 12 ==0:
return x /12
else:
return int(x/12)+1
df['year'] = df['Survival_Months'].map(months_to_year)
T = df['year']
E = df['Status']
survival_table_from_events(T, E,columns=['removed', 'observed', 'censored', 'entrance', 'at_risk'])
| removed | observed | censored | entrance | at_risk | |
|---|---|---|---|---|---|
| event_at | |||||
| 0.0 | 0 | 0 | 0 | 4024 | 4024 |
| 1.0 | 71 | 47 | 24 | 0 | 4024 |
| 2.0 | 114 | 88 | 26 | 0 | 3953 |
| 3.0 | 117 | 100 | 17 | 0 | 3839 |
| 4.0 | 225 | 118 | 107 | 0 | 3722 |
| 5.0 | 739 | 105 | 634 | 0 | 3497 |
| 6.0 | 725 | 60 | 665 | 0 | 2758 |
| 7.0 | 731 | 54 | 677 | 0 | 2033 |
| 8.0 | 674 | 33 | 641 | 0 | 1302 |
| 9.0 | 628 | 11 | 617 | 0 | 628 |
| event_at | removed | observed | censored | entrance | at_risk |
|---|---|---|---|---|---|
| 0.0 | 0 | 0 | 0 | 4024 | 4024 |
| 1.0 | 71 | 47 | 24 | 0 | 4024 |
| 2.0 | 114 | 88 | 26 | 0 | 3953 |
| 3.0 | 117 | 100 | 17 | 0 | 3839 |
| 4.0 | 225 | 118 | 107 | 0 | 3722 |
| 5.0 | 739 | 105 | 634 | 0 | 3497 |
| 6.0 | 725 | 60 | 665 | 0 | 2758 |
| 7.0 | 731 | 54 | 677 | 0 | 2033 |
| 8.0 | 674 | 33 | 641 | 0 | 1302 |
| 9.0 | 628 | 11 | 617 | 0 | 628 |
通过以上的寿命表计算出死亡概率和生存概率和生存率
| event_at | removed | observed | censored | entrance | at_risk | probability of death | probability of surviving | survival rate |
|---|---|---|---|---|---|---|---|---|
| 0.0 | 0 | 0 | 0 | 4024 | 4024 | 0 | 1 | 1 |
| 1.0 | 71 | 47 | 24 | 0 | 4024 | 0.0117 | 0.9883 | 0.9883 |
| 2.0 | 114 | 88 | 26 | 0 | 3953 | 0.0223 | 0.9777 | 0.9663 |
| 3.0 | 117 | 100 | 17 | 0 | 3839 | 0.0261 | 0.9739 | 0.941 |
| 4.0 | 225 | 118 | 107 | 0 | 3722 | 0.03217 | 0.9678 | 0.911 |
| 5.0 | 739 | 105 | 634 | 0 | 3497 | 0.033 | 0.967 | 0.881 |
| 6.0 | 725 | 60 | 665 | 0 | 2758 | 0.0247 | 0.975 | 0.858 |
| 7.0 | 731 | 54 | 677 | 0 | 2033 | 0.0319 | 0.9681 | 0.8313 |
| 8.0 | 674 | 33 | 641 | 0 | 1302 | 0.0336 | 0.9664 | 0.803 |
| 9.0 | 628 | 11 | 617 | 0 | 628 | 0.0344 | 0.9656 | 0.776 |
#然后通过计算出来的生存率绘制生存曲线
x = [0,1,2,3,4,5,6,7,8,9]
y = [1, 0.9883, 0.9663, 0.941,0.911,0.881,0.858,0.8313,0.803,0.776]
plt.figure(figsize=(8,8))
plt.title("survival curve")
plt.step(x,y,color="#8dd3c7", where="pre",lw=2)
plt.xlim(0,9)
plt.ylim(0.7,1)
plt.xlabel("year")
plt.ylabel('S(t)')
plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-m5vjpVYP-1645703648861)(output_25_0.png)]
3.3 Kaplan-Meier法
上面尝试通过寿命表来计算生存率并绘制生存曲线，在估算生存率的时候，其实用到的是乘积极限法，也就是Kaplan-Meier法，简称K-M法，且K-M法是一种无参数方法
生存率S(t)，当t=3时，意思研究对象活过3年时的生存率，是对象活过2年时的前提下计算活过3年的生存率，而知道活过2年的生存率，是研究对象活过1年的前提下计算出来。
所以S(t=3) = P1 * P2 * P3 ,P是每一年的生存概率。
原始的乳腺癌数据集的生存时间是以月位单位的，样本量有4000多个，如果通过寿命表来绘制，工作量很大。现在使用lifelines库中的KaplanMeierFitter工具就可以使用K-M法绘制生存曲线。
通过分组绘制生存曲线，展示不同因素下的研究对象的生存曲线，从而进行单因素分析。
通过KM法还能计算中位生存时间，中位生存时间就是纵坐标=0.5（即累积生存率=0.5的时候所对应的时间t），当绘制多组生存曲线的时候可以通过中位生存时间简单比较组与组之间的差异。
例如这个图，蓝色曲线的生存时间大概1200，黄色曲线的2100，黄色曲线的存活时间比蓝色的长。
#使用KM分析，并绘制生存曲线和计算中位生存时间
from lifelines import KaplanMeierFitter
from lifelines.utils import median_survival_times
kmf = KaplanMeierFitter()
kmf.fit(durations=df['Survival_Months'],
event_observed=df['Status'])
print("生成概率：")
print(kmf.survival_function_)
print("中间存活概率=0.5时的，95%区间")
median_confidence_interval = median_survival_times(kmf.confidence_interval_)
print(median_confidence_interval) #因为在研究期间生存概率没有低于0.5，所以没有中间生存概率对应的时间
plt.figure(figsize=[10,10])
kmf.plot_survival_function(at_risk_counts=True)#show_censors =True 是否显示删失 ，at_risk_counts=True 是否显示风险计数
生成概率：
KM_estimate
timeline
0.0 1.000000
1.0 1.000000
2.0 0.999503
3.0 0.999006
4.0 0.996767
... ...
103.0 0.785709
104.0 0.785709
105.0 0.785709
106.0 0.785709
107.0 0.785709
[108 rows x 1 columns]
中间存活概率=0.5时的，95%区间
KM_estimate_lower_0.95 KM_estimate_upper_0.95
0.5 inf inf
<AxesSubplot:xlabel='timeline'>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qkNumVLL-1645703648863)(output_27_2.png)]
#打印生存曲线的简短统计内容(寿命表)。包括观察值，事件数，删失数，生存率和其置信区间。
pd.concat([kmf.event_table,kmf.survival_function_,kmf.confidence_interval_],axis=1)
| removed | observed | censored | entrance | at_risk | KM_estimate | KM_estimate_lower_0.95 | KM_estimate_upper_0.95 | |
|---|---|---|---|---|---|---|---|---|
| 0.0 | 0 | 0 | 0 | 4024 | 4024 | 1.000000 | 1.000000 | 1.000000 |
| 1.0 | 1 | 0 | 1 | 0 | 4024 | 1.000000 | 1.000000 | 1.000000 |
| 2.0 | 3 | 2 | 1 | 0 | 4023 | 0.999503 | 0.998014 | 0.999876 |
| 3.0 | 4 | 2 | 2 | 0 | 4020 | 0.999006 | 0.997353 | 0.999627 |
| 4.0 | 10 | 9 | 1 | 0 | 4016 | 0.996767 | 0.994438 | 0.998121 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 103.0 | 50 | 0 | 50 | 0 | 251 | 0.785709 | 0.765872 | 0.804086 |
| 104.0 | 48 | 0 | 48 | 0 | 201 | 0.785709 | 0.765872 | 0.804086 |
| 105.0 | 45 | 0 | 45 | 0 | 153 | 0.785709 | 0.765872 | 0.804086 |
| 106.0 | 47 | 0 | 47 | 0 | 108 | 0.785709 | 0.765872 | 0.804086 |
| 107.0 | 61 | 0 | 61 | 0 | 61 | 0.785709 | 0.765872 | 0.804086 |
108 rows × 8 columns
3.4 对数据进行分组，用KM法进行单因素生存分析
使用KM法可以分析整个数据的研究对象的生存率。有的时候为了研究药物是否对疾病带来治疗功效，会设置实验组和对照组，使用KM分析法时可以分别画出两组的生存曲线来展示两组的差别，或者计算中位生存时间，假设实验组的中位生存时间大于对照组，就说明药物可以延迟患者的寿命。
但是KM法只能针对分类变量进行分组，对应连续值需要把连续值区间化，再进行分组.
比较两组之间的生存差异可以通过中位生存时间来比较。也可以指定时间点来计算P值两组曲线是否有显著性差异。也可以计算T时间内生存曲线下的面积，不同生存曲线之间的差异可以通过曲线下面积来比较（限制性平均生存时间/RMST）。也可以通过计算第p百分位生存时间来比较两组之间的生存差异
"""
对T——Stage 肿瘤大小进行分组并绘制生存曲线
从生存曲线中明显看到T4的患者（肿瘤大)的生存率下降比较快
"""
from lifelines import KaplanMeierFitter
plt.figure(figsize=(8,8))
kmf = KaplanMeierFitter()
for name, grouped_df in df.groupby('T_Stage'):
kmf.fit(grouped_df["Survival_Months"], grouped_df["Status"], label=name)
kmf.plot_survival_function()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8mLgXVID-1645703648864)(output_30_0.png)]
"""
因为年龄是连续值，km法只适用分类变量，因为需要把连续值的年龄进行区间分组
对年龄进行分组，并绘制生存曲线
"""
from lifelines import KaplanMeierFitter
df['age_'] = pd.cut(df['Age'],[0,40,50,60])
plt.figure(figsize=(8,8))
kmf = KaplanMeierFitter()
for name, grouped_df in df.groupby('age_'):
kmf.fit(grouped_df["Survival_Months"], grouped_df["Status"], label=name)
kmf.plot_survival_function()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-NVKc0lxQ-1645703648865)(output_31_0.png)]
"""
可以通过中位生存时间来比较两组之间的生存差异
（因为这个数据有些组的生存率最后都是大于0.5，无法计算中卫生存时间，所以显示inf）
"""
from lifelines.utils import median_survival_times
kmf = KaplanMeierFitter()
for name, grouped_df in df.groupby('T_Stage'):
kmf_temp = KaplanMeierFitter().fit(grouped_df["Survival_Months"], grouped_df["Status"])
median = median_survival_times(kmf_temp)
print(f"{name}组的中位生存时间:{median}")
T1组的中位生存时间:inf
T2组的中位生存时间:inf
T3组的中位生存时间:inf
T4组的中位生存时间:inf
"""
有些数据到随访结束时生存率都是大于0.5，即无法计算中位生存时间。
但是可以指定第qth分位数的生存时间，来求得对应的生存时间
"""
from lifelines.utils import qth_survival_time
qth = 0.9
kmf = KaplanMeierFitter()
for name, grouped_df in df.groupby('T_Stage'):
kmf_temp = KaplanMeierFitter().fit(grouped_df["Survival_Months"], grouped_df["Status"])
time = qth_survival_time(qth,kmf_temp)
print(f"{name}组的第{qth*100}分位数的生存时间:{time}")
T1组的第90.0分位数的生存时间:81.0
T2组的第90.0分位数的生存时间:49.0
T3组的第90.0分位数的生存时间:40.0
T4组的第90.0分位数的生存时间:25.0
/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/lifelines/fitters/__init__.py:273: ApproximationWarning: Approximating using `survival_function_`. To increase accuracy, try using or increasing the resolution of the timeline kwarg in `.fit(..., timeline=timeline)`.
exceptions.ApproximationWarning,
/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/lifelines/fitters/__init__.py:273: ApproximationWarning: Approximating using `survival_function_`. To increase accuracy, try using or increasing the resolution of the timeline kwarg in `.fit(..., timeline=timeline)`.
exceptions.ApproximationWarning,
/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/lifelines/fitters/__init__.py:273: ApproximationWarning: Approximating using `survival_function_`. To increase accuracy, try using or increasing the resolution of the timeline kwarg in `.fit(..., timeline=timeline)`.
exceptions.ApproximationWarning,
/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/lifelines/fitters/__init__.py:273: ApproximationWarning: Approximating using `survival_function_`. To increase accuracy, try using or increasing the resolution of the timeline kwarg in `.fit(..., timeline=timeline)`.
exceptions.ApproximationWarning,
"""
可以通过指定时间来比较这个时间对应的生存差异（例如我关心的时第60个月的生存差异)
"""
from lifelines.statistics import survival_difference_at_fixed_point_in_time_test
from lifelines.datasets import load_waltons
T = df['Survival_Months']
E = df['Status']
ix = df['A_Stage']== 0
kmf_1 = KaplanMeierFitter().fit(T[ix], E[ix])
kmf_2 = KaplanMeierFitter().fit(T[~ix], E[~ix])
point_in_time = 60. #指定比较的时间
results = survival_difference_at_fixed_point_in_time_test(point_in_time, kmf_1, kmf_2)
results.print_summary()
| null_distribution | chi squared |
|---|---|
| degrees_of_freedom | 1 |
| point_in_time | 60 |
| fitterA | <lifelines.KaplanMeierFitter:"KM_estimate", fi... |
| fitterB | <lifelines.KaplanMeierFitter:"KM_estimate", fi... |
| test_name | survival_difference_at_fixed_point_in_time_test |
| test_statistic | p | -log2(p) | |
|---|---|---|---|
| 0 | 40.84 | <0.005 | 32.49 |
"""
限制性平均生存时间/RMST
通过计算不同生存曲线下面积来比较生存差异
"""
from lifelines.utils import restricted_mean_survival_time
from lifelines.plotting import rmst_plot
T = df['Survival_Months']
E = df['Status']
ix = df['A_Stage']== 0 #0是'Regional'
time_limit = 60#指定T时间内
kmf_1 = KaplanMeierFitter().fit(T[ix], E[ix])
rmst_1 = restricted_mean_survival_time(kmf_1, t=time_limit)
print(f"曲线下面积:{rmst_1}")
kmf_2 = KaplanMeierFitter().fit(T[~ix], E[~ix])
rmst_2 = restricted_mean_survival_time(kmf_2, t=time_limit)
print(f"曲线下面积:{rmst_2}")
#把曲线面积绘制出来
plt.figure(figsize=(10,12))
ax = plt.subplot(311)
rmst_plot(kmf_1, t=time_limit, ax=ax)
ax = plt.subplot(312)
rmst_plot(kmf_2, t=time_limit, ax=ax)
ax = plt.subplot(313)
rmst_plot(kmf_1, model2=kmf_2, t=time_limit, ax=ax)
曲线下面积:57.22662022019032
曲线下面积:49.82516339891577
<AxesSubplot:xlabel='timeline'>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-lqgOVtUY-1645703648866)(output_35_2.png)]
3.5 Log-rank检验
在使用KM方法，根据分组绘制多条生存曲线后，只通过直接的观察来确定多条曲线之间是否具有显著性差异是不充分的。Log-rank检验弥补这方面。
Log-rank检验的零假设组与组的生存时间分布没有显著性差异。当计算出来的P值小于0.05,推翻零假设，认为两组的生存时间分布存在显著性差异。
"""
比较两个组的生存曲线是否有显著性差异
结果A_Stage组的P值小于0.005，是有显著性差异(直接观察图形也看出明显差异)
"""
from lifelines.statistics import logrank_test
from lifelines import KaplanMeierFitter
plt.figure(figsize=(8,8))
kmf = KaplanMeierFitter()
for name, grouped_df in df.groupby('A_Stage'):
kmf.fit(grouped_df["Survival_Months"], grouped_df["Status"], label=name)
kmf.plot_survival_function()
T = df['Survival_Months']
E = df['Status']
ix = df['A_Stage']== 'Distant'
result = logrank_test(T[ix],T[~ix], E[ix],E[~ix])
result.print_summary()
| t_0 | -1 |
|---|---|
| null_distribution | chi squared |
| degrees_of_freedom | 0 |
| test_name | logrank_test |
| test_statistic | p | -log2(p) | |
|---|---|---|---|
| 0 | 0.00 | nan | nan |
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-R66tN956-1645703648868)(output_37_1.png)]
"""
比较两个组以上的生存曲线是否有显著性差异
结果T_Stage组的P值小于0.005，是有显著性差异(直接观察图形也看出明显差异)
"""
from lifelines.statistics import multivariate_logrank_test
kmf = KaplanMeierFitter()
for name, grouped_df in df.groupby('N_Stage'):
kmf.fit(grouped_df["Survival_Months"], grouped_df["Status"], label=name)
kmf.plot_survival_function()
results = multivariate_logrank_test(df['Survival_Months'], df['N_Stage'], df['Status'])
results.print_summary()
print(results.p_value)
| t_0 | -1 |
|---|---|
| null_distribution | chi squared |
| degrees_of_freedom | 2 |
| test_name | multivariate_logrank_test |
| test_statistic | p | -log2(p) | |
|---|---|---|---|
| 0 | 298.86 | <0.005 | 215.58 |
1.2701885843852175e-65
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rNjOuauR-1645703648869)(output_38_2.png)]
3.6 参数估计
如果事先知道生存时间的分布，可以使用参数法：根据样本观测值来估计假定的分布模型中的参数，获得生存时间的概率分布模型。
生存时间经常服从的分布有：指数分布、Weibull分布、对数正态分布、对数Logistic分布、Gamma分布。分析过程类似KM法，由于往往对真实数据的生存时间分布并不知道，所有参数估计并不是常用的生存分析方法。
"""
简单展示如果用linelifes来使用参数法进行生存分析并绘制生存曲线。
"""
from lifelines.datasets import load_waltons
from lifelines.utils import median_survival_times
import matplotlib.pyplot as plt
import numpy as np
from lifelines import *
# 数据载入
T = df['Survival_Months']
E = df['Status']
fig, axes = plt.subplots(3, 3, figsize=(13.5, 7.5))
# 多种参数模型
kmf = KaplanMeierFitter().fit(T, E, label='KaplanMeierFitter')
wbf = WeibullFitter().fit(T, E, label='WeibullFitter')
exf = ExponentialFitter().fit(T, E, label='ExponentialFitter')
lnf = LogNormalFitter().fit(T, E, label='LogNormalFitter')
llf = LogLogisticFitter().fit(T, E, label='LogLogisticFitter')
pwf = PiecewiseExponentialFitter([40, 60]).fit(T, E, label='PiecewiseExponentialFitter')
ggf = GeneralizedGammaFitter().fit(T, E, label='GeneralizedGammaFitter')
sf = SplineFitter(np.percentile(T.loc[E.astype(bool)], [0, 50, 100])).fit(T, E, label='SplineFitter')
wbf.plot_survival_function(ax=axes[0][0])
exf.plot_survival_function(ax=axes[0][1])
lnf.plot_survival_function(ax=axes[0][2])
kmf.plot_survival_function(ax=axes[1][0])
llf.plot_survival_function(ax=axes[1][1])
pwf.plot_survival_function(ax=axes[1][2])
ggf.plot_survival_function(ax=axes[2][0])
sf.plot_survival_function(ax=axes[2][1])
<AxesSubplot:>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-RJkcN4V3-1645703648869)(output_40_1.png)]
3.7 比例风险回归模型–Cox模型
Cox比例风险回归模型是英国统计学家D.R.Cox（很遗憾，2022年1月，D.R.Cox去世了）在1972年提出的，用于肿瘤基本的预后分析。
Cox比例回归的公式如下：
其中，t表示生存时间，x为协变量，b为协变量的回归系数，h0为基准风险，表示所有协变量为0时的个体在t时间发生死亡的风险率。而h(t)是研究对象的风险函数，可以理解成时间t死亡额风险。exp(bi)为协变量的xi的风险比，风险比的值大于1的协变量为危险因素，大于1的是保护因素。
Cox比例风险回归模型 ，属于半参数法，可以估计生存函数，可以比较两组或多组生存分布函数，可以单因素分析也可以多因素分析。单因素分析包括分类变量和连续变量(KM法只能分析分类变量)，可以建立生存时间与影响因素之间的关系模型，计算影响因素的分险比。注意Cox回归要求满足比例风险假定，就是影响因素是不随时间变化而变化。假如变量不满足，则应对变量进行分层再进行cox回归。
3.8 对数据进行划分成训练集和验证集
先把数据中有缺失的行删除，然后8:2对数据进行分割成训练集和验证集
#加载数据
df = pd.read_csv('SEER Breast Cancer Dataset .csv')
df['Race'] = pd.factorize(df['Race'])[0].astype(np.uint16)
df['Marital_Status'] = pd.factorize(df['Marital_Status'])[0].astype(np.uint16)
df['A_Stage'] = pd.factorize(df['A_Stage'])[0].astype(np.uint16)
df['T_Stage'] = df['T_Stage'].map({'T1':1, 'T2':2,'T3':3,'T4':4})
df['N_Stage'] = df['N_Stage'].map({'N1':0, 'N2':1,'N3':2})
df['Sixth_Stage'] = df['Sixth_Stage'].map({'IIA':0, 'IIB':1,'IIIA':2,'IIIB':3,'IIIC':4})
df['Grade'] = df['Grade'].map({'Well differentiated; Grade I':0,
'Moderately differentiated; Grade II':1,
'Poorly differentiated; Grade III':2,
'Undifferentiated; anaplastic; Grade IV':3})
df['Estrogen_Status'] = df['Estrogen_Status'].map({'Positive':0, 'Negative':1})
df['Progesterone_Status'] = df['Progesterone_Status'].map({'Positive':0, 'Negative':1})
df['Status'] = df['Status'].map({'Alive':0, 'Dead':1})
#删除具有缺失值的行
df = df.dropna(axis=0)
#分割训练集和验证集
df_train = df.sample(frac=0.8,random_state=2022)
df_val = df.drop(index = df_train.index)
print(f"训练集大小：{len(df_train)}")
print(f"验证集大小：{len(df_val)}")
print("*"*20)
#统计失效事件和删失数据的数量
print(df_train['Status'].value_counts())
print(df_val['Status'].value_counts())
训练集大小：3219
验证集大小：805
********************
0 2740
1 479
Name: Status, dtype: int64
0 668
1 137
Name: Status, dtype: int64
3.9 建立Cox模型并开始训练
通过print_summary()模型会输出三个表格
表格一：duration col是时间列名，event col是事件列名， breslow是模型参数估计方法，3219是总数据量，479是发生失效事件的数量
| model | lifelines.CoxPHFitter |
|---|---|
| duration col | ‘Survival_Months’ |
| event col | ‘Status’ |
| baseline estimation | breslow |
| number of observations | 3219 |
| number of events observed | 479 |
| partial log-likelihood | -3516.57 |
| time fit was run | 2022-02-06 03:27:20 UTC |
表格二：coef是协变量的回归系数，exp(coef)是风险比，se(coef)是系数的标准误差,标准误差的95%置信区，分险比的95%置信区，Wald统计值(z),Wald统计值的P值是否显著。
如何解读风险比？当风险比是大于1属于危险因素，值增大会增加患者的死亡风险，小于1属于保护因素，降低患者的死亡风险。
例如Age的p值是小于0.005，风险比HR = 1.02，95%置信区为1.01-1.03，表明Age值与死亡风险增加之间的有关系。保持其他协变量不变，Age增加一个单位，患者的风险增加2%。
| coef | exp(coef) | se(coef) | coef lower 95% | coef upper 95% | exp(coef) lower 95% | exp(coef) upper 95% | z | p | -log2§ | |
|---|---|---|---|---|---|---|---|---|---|---|
| Age | 0.02 | 1.02 | 0.01 | 0.01 | 0.03 | 1.01 | 1.03 | 3.09 | <0.005 | 8.95 |
| Race | 0.36 | 1.43 | 0.12 | 0.13 | 0.59 | 1.14 | 1.80 | 3.09 | <0.005 | 8.96 |
| Marital_Status | 0.07 | 1.08 | 0.04 | -0.01 | 0.16 | 0.99 | 1.17 | 1.76 | 0.08 | 3.69 |
表格三：Concordance 是一致性指数，用了全部协变量训练cox模型，一致性指数为0.75
| Concordance | 0.75 |
|---|---|
| Partial AIC | 7059.14 |
| log-likelihood ratio test | 398.31 on 13 df |
| -log2§ of ll-ratio test | 253.44 |
"""
单变量分析
"""
from lifelines import CoxPHFitter
varList = df_train.columns.to_list()
varList.pop()#去掉时间列
varList.pop()#去掉状态列
#创建Cox回归模型
cph = CoxPHFitter()
for var in varList:
cph.fit(df=df_train,duration_col='Survival_Months',event_col='Status',formula=var)
print(f"{var} 的P值：{cph.summary.p[var]},<0.05 {cph.summary.p[var]<0.05}")
Age 的P值：0.0248649940872201,<0.05 True
Race 的P值：1.3550638342845673e-06,<0.05 True
Marital_Status 的P值：8.68912013912815e-05,<0.05 True
T_Stage 的P值：2.2304661752517098e-20,<0.05 True
N_Stage 的P值：1.6617003098137667e-45,<0.05 True
Sixth_Stage 的P值：1.0290156463874445e-47,<0.05 True
Grade 的P值：6.022575977628071e-23,<0.05 True
A_Stage 的P值：1.0633163182382273e-08,<0.05 True
Tumor_Size 的P值：5.79647504244188e-16,<0.05 True
Estrogen_Status 的P值：1.0415072293170037e-27,<0.05 True
Progesterone_Status 的P值：5.007000090288116e-25,<0.05 True
Regional_Node_Examined 的P值：0.1065243965415506,<0.05 False
Reginol_Node_Positive 的P值：1.5269471157224244e-46,<0.05 True
"""
多变量分析
"""
from lifelines import CoxPHFitter
#创建Cox回归模型
cph = CoxPHFitter()
#df传入DataFrame格式数据，duration_col传入时间的列名，event_col传入事件的列名，默认是使用全部协变量
#可以通过参数formula='Age+Race'，传入部分协变量
cph.fit(df=df_train,duration_col='Survival_Months',event_col='Status',show_progress=True)
#打印模型情况
cph.print_summary()
Iteration 1: norm_delta = 1.07058, step_size = 0.9000, log_lik = -3715.72192, newton_decrement = 255.86131, seconds_since_start = 0.2
Iteration 2: norm_delta = 0.47681, step_size = 0.9000, log_lik = -3575.51461, newton_decrement = 52.18325, seconds_since_start = 0.4
Iteration 3: norm_delta = 0.12309, step_size = 0.9000, log_lik = -3520.88106, newton_decrement = 3.96398, seconds_since_start = 0.6
Iteration 4: norm_delta = 0.01909, step_size = 1.0000, log_lik = -3516.63264, newton_decrement = 0.06359, seconds_since_start = 0.9
Iteration 5: norm_delta = 0.00045, step_size = 1.0000, log_lik = -3516.56817, newton_decrement = 0.00003, seconds_since_start = 1.1
Iteration 6: norm_delta = 0.00000, step_size = 1.0000, log_lik = -3516.56814, newton_decrement = 0.00000, seconds_since_start = 1.3
Convergence success after 6 iterations.
| model | lifelines.CoxPHFitter |
|---|---|
| duration col | 'Survival_Months' |
| event col | 'Status' |
| baseline estimation | breslow |
| number of observations | 3219 |
| number of events observed | 479 |
| partial log-likelihood | -3516.57 |
| time fit was run | 2022-02-23 17:30:27 UTC |
| coef | exp(coef) | se(coef) | coef lower 95% | coef upper 95% | exp(coef) lower 95% | exp(coef) upper 95% | z | p | -log2(p) | |
|---|---|---|---|---|---|---|---|---|---|---|
| Age | 0.02 | 1.02 | 0.01 | 0.01 | 0.03 | 1.01 | 1.03 | 3.09 | <0.005 | 8.95 |
| Race | 0.36 | 1.43 | 0.12 | 0.13 | 0.59 | 1.14 | 1.80 | 3.09 | <0.005 | 8.96 |
| Marital_Status | 0.07 | 1.08 | 0.04 | -0.01 | 0.16 | 0.99 | 1.17 | 1.76 | 0.08 | 3.69 |
| T_Stage | 0.37 | 1.45 | 0.11 | 0.16 | 0.58 | 1.17 | 1.79 | 3.42 | <0.005 | 10.63 |
| N_Stage | 0.36 | 1.43 | 0.18 | 0.02 | 0.70 | 1.02 | 2.02 | 2.05 | 0.04 | 4.63 |
| Sixth_Stage | -0.00 | 1.00 | 0.11 | -0.23 | 0.22 | 0.80 | 1.25 | -0.01 | 0.99 | 0.01 |
| Grade | 0.47 | 1.60 | 0.08 | 0.32 | 0.63 | 1.37 | 1.87 | 6.02 | <0.005 | 29.06 |
| A_Stage | 0.00 | 1.00 | 0.21 | -0.41 | 0.41 | 0.66 | 1.51 | 0.01 | 0.99 | 0.01 |
| Tumor_Size | -0.00 | 1.00 | 0.00 | -0.01 | 0.00 | 0.99 | 1.00 | -0.64 | 0.52 | 0.93 |
| Estrogen_Status | 0.54 | 1.72 | 0.15 | 0.24 | 0.84 | 1.28 | 2.31 | 3.58 | <0.005 | 11.50 |
| Progesterone_Status | 0.55 | 1.73 | 0.12 | 0.31 | 0.78 | 1.37 | 2.18 | 4.58 | <0.005 | 17.71 |
| Regional_Node_Examined | -0.03 | 0.97 | 0.01 | -0.05 | -0.02 | 0.95 | 0.98 | -4.69 | <0.005 | 18.49 |
| Reginol_Node_Positive | 0.05 | 1.05 | 0.01 | 0.03 | 0.08 | 1.03 | 1.08 | 4.10 | <0.005 | 14.58 |
| Concordance | 0.75 |
|---|---|
| Partial AIC | 7059.14 |
| log-likelihood ratio test | 398.31 on 13 df |
| -log2(p) of ll-ratio test | 253.44 |
#画出每个协变量的HR值和95%置信度区间（森林图)
plt.figure(figsize=(8,8))
cph.plot(hazard_ratios=True)
"""
结合上面的print_summary()表格看出，协变量的95%置信度区间跨过横坐标轴上1的协变量的p值都是大于0.05
"""
'\n结合上面的print_summary()表格看出，协变量的95%置信度区间跨过横坐标轴上1的协变量的p值都是大于0.05\n'
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-hxYOl86A-1645703648871)(output_47_1.png)]
#去掉不显著的变量重新拟合cox模型
FORMULA = 'Age + Race + T_Stage + N_Stage+Grade+Estrogen_Status+Progesterone_Status+Regional_Node_Examined+Reginol_Node_Positive'
cph.fit(df=df_train,duration_col='Survival_Months',event_col='Status',formula=FORMULA)
#打印模型情况
cph.print_summary()
| model | lifelines.CoxPHFitter |
|---|---|
| duration col | 'Survival_Months' |
| event col | 'Status' |
| baseline estimation | breslow |
| number of observations | 3219 |
| number of events observed | 479 |
| partial log-likelihood | -3518.31 |
| time fit was run | 2022-02-23 17:30:49 UTC |
| coef | exp(coef) | se(coef) | coef lower 95% | coef upper 95% | exp(coef) lower 95% | exp(coef) upper 95% | z | p | -log2(p) | |
|---|---|---|---|---|---|---|---|---|---|---|
| Age | 0.02 | 1.02 | 0.01 | 0.01 | 0.03 | 1.01 | 1.03 | 3.32 | <0.005 | 10.13 |
| Estrogen_Status | 0.55 | 1.73 | 0.15 | 0.25 | 0.84 | 1.29 | 2.33 | 3.65 | <0.005 | 11.87 |
| Grade | 0.47 | 1.61 | 0.08 | 0.32 | 0.63 | 1.38 | 1.87 | 6.06 | <0.005 | 29.43 |
| N_Stage | 0.35 | 1.42 | 0.09 | 0.16 | 0.53 | 1.18 | 1.71 | 3.68 | <0.005 | 12.05 |
| Progesterone_Status | 0.55 | 1.73 | 0.12 | 0.31 | 0.78 | 1.37 | 2.18 | 4.57 | <0.005 | 17.68 |
| Race | 0.40 | 1.49 | 0.12 | 0.17 | 0.62 | 1.19 | 1.87 | 3.45 | <0.005 | 10.78 |
| Reginol_Node_Positive | 0.05 | 1.05 | 0.01 | 0.03 | 0.08 | 1.03 | 1.08 | 4.31 | <0.005 | 15.89 |
| Regional_Node_Examined | -0.03 | 0.97 | 0.01 | -0.05 | -0.02 | 0.95 | 0.98 | -4.68 | <0.005 | 18.42 |
| T_Stage | 0.33 | 1.39 | 0.06 | 0.22 | 0.44 | 1.24 | 1.56 | 5.80 | <0.005 | 27.13 |
| Concordance | 0.75 |
|---|---|
| Partial AIC | 7054.62 |
| log-likelihood ratio test | 394.83 on 9 df |
| -log2(p) of ll-ratio test | 261.63 |
#绘制生存曲线S(t)
cph.baseline_survival_.plot()
<AxesSubplot:>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AXbYqi6u-1645703648872)(output_49_1.png)]
#累积风险曲线H(t)
cph.baseline_cumulative_hazard_.plot()
<AxesSubplot:>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Ai7eyoCy-1645703648873)(output_50_1.png)]
#因为累积风险函数 H(t) = -lnS(t)
#现在计算H(t) 再绘制曲线，结果和上面的图一致。
(-np.log(cph.baseline_survival_)).plot()
<AxesSubplot:>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GlExCSVJ-1645703648874)(output_51_1.png)]
#绘制每个时间点对应的风险函数h(t)
cph.baseline_hazard_.plot()
<AxesSubplot:>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-u17GVfeN-1645703648874)(output_52_1.png)]
print(cph.predict_partial_hazard(pd.DataFrame(df_val.iloc[0]).T))#预测个体的风险
print(cph.predict_percentile(pd.DataFrame(df_val.iloc[10]).T,p=0.9))#预测个体第90.0分位数的生存时间
print(cph.predict_median(pd.DataFrame(df_val.iloc[10]).T))#预测个体的中位生存时间
print(cph.predict_cumulative_hazard(df_val.iloc[0]))#预测个体的累积风险函数
print(cph.predict_survival_function(df_val.iloc[0]))#预测个体的生存函数
1 0.573792
dtype: float64
96.0
inf
1
1.0 0.000000
2.0 0.000118
3.0 0.000353
4.0 0.001059
5.0 0.001414
... ...
103.0 0.103245
104.0 0.103245
105.0 0.103245
106.0 0.103245
107.0 0.103245
[107 rows x 1 columns]
1
1.0 1.000000
2.0 0.999882
3.0 0.999647
4.0 0.998941
5.0 0.998587
... ...
103.0 0.901906
104.0 0.901906
105.0 0.901906
106.0 0.901906
107.0 0.901906
[107 rows x 1 columns]
from lifelines.utils import concordance_index
#计算c-index
#方法一
print(f"训练集的c-index:{cph.score(df_train, scoring_method='concordance_index')}")
print(f"验证集的c-index:{cph.score(df_val, scoring_method='concordance_index')}")
#方法二
print(concordance_index(df_train['Survival_Months'], -cph.predict_partial_hazard(df_train), df_train['Status']))
#预测验证集中某个研究对象的生存率
cph.predict_survival_function(df_val.iloc[0])
训练集的c-index:0.7460170801920486
验证集的c-index:0.7108462839198003
0.7460166658656356
| 1 | |
|---|---|
| 1.0 | 1.000000 |
| 2.0 | 0.999882 |
| 3.0 | 0.999647 |
| 4.0 | 0.998941 |
| 5.0 | 0.998587 |
| ... | ... |
| 103.0 | 0.901906 |
| 104.0 | 0.901906 |
| 105.0 | 0.901906 |
| 106.0 | 0.901906 |
| 107.0 | 0.901906 |
107 rows × 1 columns
"""
可以通过k折交叉验证，把所有协变量的组合都尝试一遍，找到最好的组合
这里的组合有8190多个，运行时间会很长很长，可以不必运行，不影响整个项目
"""
from lifelines.utils import k_fold_cross_validation
from itertools import combinations
def combine(temp_list, n):
'''根据n获得列表中的所有可能组合（n个元素为一组）'''
temp_list2 = []
for c in combinations(temp_list, n):
temp_list2.append(c)
return temp_list2
varList = df_train.columns.to_list()
varList.pop()#去掉"Status"
varList.pop()#去掉"Survival_Months"
end_list = []
for i in range(len(varList)):
if combine(varList, i)!=[()]:
end_list.extend(combine(varList, i))
#打印组合数量
print(len(end_list))
mean_score = []#保存每个组合的k折交叉的c-index得分
for i in set(end_list):
columns = list(i)+['Survival_Months','Status']
df_temp = df_train.loc[:,columns]
scores = k_fold_cross_validation(CoxPHFitter(), df_temp, 'Survival_Months', event_col='Status', k=10,scoring_method="concordance_index")
mean_score.append(np.mean(scores))
#打印最高得分
print(np.max(mean_score))
#打印最高得分的协变量组合
print(end_list[mean_score.index(np.max(mean_score))])
3.10 协变量值如何影响生存曲线
"""
将模型的基线生存曲线与一组中协变量值发生变化时发生的情况进行比较。
当我们改变协变量(s)时，在其他条件相同的情况下，这有助于比较受试者的生存期。
"""
cph.plot_partial_effects_on_outcome('T_Stage',values=[1,2,3,4])
<AxesSubplot:>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Q7TFBpGb-1645703648875)(output_57_1.png)]
"""
观测多个协变量的值是如何影响生存曲线
"""
cph.plot_partial_effects_on_outcome(['T_Stage', 'N_Stage'], values=[[0, 0], [1, 0], [1,1], [1,2]], cmap='coolwarm')
<AxesSubplot:>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-oRVpOhzx-1645703648876)(output_58_1.png)]
#绘制平滑校准曲线
from lifelines.calibration import survival_probability_calibration
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
#训练集的校准曲线，t0 ( float ) – 评估事件发生概率的时间
survival_probability_calibration(cph,df_train,t0=36,ax=axes[0])
#验证集的校准曲线
survival_probability_calibration(cph,df_val,t0=36,ax=axes[1])
"""
ICI – 预测和观察到的平均绝对差
E50 – 预测和观察到的绝对差中位数
"""
ICI = 0.0033427008106428234
E50 = 0.0022028810234795415
ICI = 0.008988198843408044
E50 = 0.008933834566536292
'\nICI – 预测和观察到的平均绝对差\nE50 – 预测和观察到的绝对差中位数\n'
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qDMleuTf-1645703648877)(output_59_2.png)]
3.11 检验比例风险(PH)假定
之前说过Cox回归模型需要满足比例风险假定，即协变量不是随时间变化而变化，例如某个协变量是是否接受某种药剂的治疗，但是这种治疗不是只发生一次，而是随着时间的增加，治疗的药量会越来越多。不像婚姻状态、人种、性别等协变量是不随时间变化的。
【如何检验PH假定】
满足PH假定和不满足PH假定，建立的cox回归模型需要特别的处理，这里可以通过check_assumptions()工具来检验所有协变量是否满足PH假定。这个工具通过统计检验法，将给出统计学检验的P值，通过P值判断。一般P小于0.05的协变量是不满足比例分险假定。
check_assumptions（）工具输出有三部分
第一部分：
| test_statistic | p | -log2§ | ||
|---|---|---|---|---|
| A_Stage | km | 4.91 | 0.03 | 5.23 |
| rank | 4.43 | 0.04 | 4.83 | |
| Estrogen_Status | km | 4.99 | 0.03 | 5.29 |
| rank | 6.11 | 0.01 | 6.22 |
零假设变量是满足PH假定，P值小于0.05时，拒绝零假设，即此变量不满足PH假定，例如表格的Estrogen_Status变量。
第二部分：
2. Variable 'Estrogen_Status' failed the non-proportional test: p-value is 0.0134.
Advice: with so few unique values (only 2), you can include `strata=['Estrogen_Status', ...]` in
the call in `.fit`. See documentation in link [E] below.
Bootstrapping lowess lines. May take a moment...
直接指出那些协变量是不满足PH假定，还给出了解决方法。
第三部分：
Estrogen_Status变量的 Schoenfeld residuals 图
假如变量是满足PH假定，再图中黑色的实线（变量影响系数随时间的变化）应该平行于蓝色的水平虚线。反之，如图那样，随着时间的变化，黑色的实线也发生变化。
【不满足PH假定如何解决】
-
方法一：【分层】对这个不满足PH假定的变量进行分层，对剩余的协变量进行cox建模
-
方法二：【时依变量协变量】
假如有不满足PH假定，应该采用含依时协变量的Cox回归进行拟合模型。
cph.check_assumptions(df_train, p_value_threshold=0.05, show_plots=True)
The ``p_value_threshold`` is set at 0.05. Even under the null hypothesis of no violations, some
covariates will be below the threshold by chance. This is compounded when there are many covariates.
Similarly, when there are lots of observations, even minor deviances from the proportional hazard
assumption will be flagged.
With that in mind, it's best to use a combination of statistical tests and visual tests to determine
the most serious violations. Produce visual plots using ``check_assumptions(..., show_plots=True)``
and looking for non-constant lines. See link [A] below for a full example.
| null_distribution | chi squared |
|---|---|
| degrees_of_freedom | 1 |
| model | <lifelines.CoxPHFitter: fitted with 3219 total... |
| test_name | proportional_hazard_test |
| test_statistic | p | -log2(p) | ||
|---|---|---|---|---|
| Age | km | 3.22 | 0.07 | 3.78 |
| rank | 4.69 | 0.03 | 5.04 | |
| Estrogen_Status | km | 4.99 | 0.03 | 5.29 |
| rank | 6.11 | 0.01 | 6.22 | |
| Grade | km | 0.03 | 0.86 | 0.22 |
| rank | 0.01 | 0.91 | 0.13 | |
| N_Stage | km | 1.31 | 0.25 | 1.99 |
| rank | 1.26 | 0.26 | 1.93 | |
| Progesterone_Status | km | 10.13 | <0.005 | 9.42 |
| rank | 9.77 | <0.005 | 9.14 | |
| Race | km | 0.22 | 0.64 | 0.65 |
| rank | 0.07 | 0.79 | 0.34 | |
| Reginol_Node_Positive | km | 0.85 | 0.36 | 1.49 |
| rank | 0.84 | 0.36 | 1.48 | |
| Regional_Node_Examined | km | 0.00 | 0.95 | 0.07 |
| rank | 0.04 | 0.84 | 0.25 | |
| T_Stage | km | 0.00 | 0.95 | 0.07 |
| rank | 0.03 | 0.86 | 0.21 |
1. Variable 'Age' failed the non-proportional test: p-value is 0.0304.
Advice 1: the functional form of the variable 'Age' might be incorrect. That is, there may be
non-linear terms missing. The proportional hazard test used is very sensitive to incorrect
functional forms. See documentation in link [D] below on how to specify a functional form.
Advice 2: try binning the variable 'Age' using pd.cut, and then specify it in `strata=['Age',
...]` in the call in `.fit`. See documentation in link [B] below.
Advice 3: try adding an interaction term with your time variable. See documentation in link [C]
below.
Bootstrapping lowess lines. May take a moment...
2. Variable 'Estrogen_Status' failed the non-proportional test: p-value is 0.0134.
Advice: with so few unique values (only 2), you can include `strata=['Estrogen_Status', ...]` in
the call in `.fit`. See documentation in link [E] below.
Bootstrapping lowess lines. May take a moment...
3. Variable 'Progesterone_Status' failed the non-proportional test: p-value is 0.0015.
Advice: with so few unique values (only 2), you can include `strata=['Progesterone_Status', ...]`
in the call in `.fit`. See documentation in link [E] below.
Bootstrapping lowess lines. May take a moment...
---
[A] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html
[B] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Bin-variable-and-stratify-on-it
[C] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Introduce-time-varying-covariates
[D] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Modify-the-functional-form
[E] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Stratification
[[<AxesSubplot:xlabel='rank-transformed time\n(p=0.0304)'>,
<AxesSubplot:xlabel='km-transformed time\n(p=0.0729)'>],
[<AxesSubplot:xlabel='rank-transformed time\n(p=0.0134)'>,
<AxesSubplot:xlabel='km-transformed time\n(p=0.0256)'>],
[<AxesSubplot:xlabel='rank-transformed time\n(p=0.0018)'>,
<AxesSubplot:xlabel='km-transformed time\n(p=0.0015)'>]]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GxkuF4CL-1645703648878)(output_61_4.png)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-guK9Oao6-1645703648879)(output_61_5.png)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mXXEf9qq-1645703648880)(output_61_6.png)]
通过check_assumptions()工具PH检验，发现一共有三个协变量是不满足PH假定，分别是Age(连续变量) 、Estrogen_Status(分类变量)、Progesterone_Status(分类变量)。对于这些变量可以通过"分层"来处理。例如对总样本根据Estrogen_Status分类变量对数据进行划分子样本，再对所有子样本进行cox建模。连续变量Age也可以通过分层来处理，先把连续变量通过区间划分成多个等级(连续变量->分类变量)，再进行对所有子样本进行cox建模。
"""
先处理分类变量Estrogen_Status和Progesterone_Status
"""
FORMULA = 'Age + Race + T_Stage + N_Stage+Grade+Regional_Node_Examined+Reginol_Node_Positive'
cph_wexp = CoxPHFitter()
cph_wexp.fit(df_train, 'Survival_Months','Status',formula =FORMULA,strata=['Estrogen_Status','Progesterone_Status'])
#现在剩下 Age是不满足PH假定
cph_wexp.check_assumptions(df_train, show_plots=True)
The ``p_value_threshold`` is set at 0.01. Even under the null hypothesis of no violations, some
covariates will be below the threshold by chance. This is compounded when there are many covariates.
Similarly, when there are lots of observations, even minor deviances from the proportional hazard
assumption will be flagged.
With that in mind, it's best to use a combination of statistical tests and visual tests to determine
the most serious violations. Produce visual plots using ``check_assumptions(..., show_plots=True)``
and looking for non-constant lines. See link [A] below for a full example.
| null_distribution | chi squared |
|---|---|
| degrees_of_freedom | 1 |
| model | <lifelines.CoxPHFitter: fitted with 3219 total... |
| test_name | proportional_hazard_test |
| test_statistic | p | -log2(p) | ||
|---|---|---|---|---|
| Age | km | 3.27 | 0.07 | 3.82 |
| rank | 6.04 | 0.01 | 6.16 | |
| Grade | km | 0.02 | 0.89 | 0.18 |
| rank | 0.89 | 0.35 | 1.54 | |
| N_Stage | km | 1.41 | 0.24 | 2.09 |
| rank | 0.05 | 0.82 | 0.28 | |
| Race | km | 0.39 | 0.53 | 0.91 |
| rank | 0.72 | 0.40 | 1.34 | |
| Reginol_Node_Positive | km | 1.08 | 0.30 | 1.74 |
| rank | 1.06 | 0.30 | 1.72 | |
| Regional_Node_Examined | km | 0.00 | 0.97 | 0.04 |
| rank | 0.34 | 0.56 | 0.84 | |
| T_Stage | km | 0.01 | 0.91 | 0.14 |
| rank | 0.01 | 0.94 | 0.09 |
1. Variable 'Age' failed the non-proportional test: p-value is 0.0140.
Advice 1: the functional form of the variable 'Age' might be incorrect. That is, there may be
non-linear terms missing. The proportional hazard test used is very sensitive to incorrect
functional forms. See documentation in link [D] below on how to specify a functional form.
Advice 2: try binning the variable 'Age' using pd.cut, and then specify it in `strata=['Age',
...]` in the call in `.fit`. See documentation in link [B] below.
Advice 3: try adding an interaction term with your time variable. See documentation in link [C]
below.
Bootstrapping lowess lines. May take a moment...
---
[A] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html
[B] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Bin-variable-and-stratify-on-it
[C] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Introduce-time-varying-covariates
[D] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Modify-the-functional-form
[E] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Stratification
[[<AxesSubplot:xlabel='rank-transformed time\n(p=0.0140)'>,
<AxesSubplot:xlabel='km-transformed time\n(p=0.0707)'>]]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-HoIfvY9C-1645703648880)(output_63_4.png)]
"""
现在处理连续变量 age，用区间对age进行划分，使变成分类变量.
用箱线图展示age的年龄分布，最小30岁，最大约68岁，集中分布在47到62岁之间
"""
df_train.boxplot(column='Age')
df_train2 = df_train.copy()
df_train2['age_strata'] = pd.cut(df_train2['Age'], np.array([30, 45, 60,75]))
df_train2['age_strata'] = pd.factorize(df_train2['age_strata'])[0].astype(np.uint16)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-FQq8tKTO-1645703648881)(output_64_0.png)]
FORMULA = 'Race + T_Stage + N_Stage+Grade+Regional_Node_Examined+Reginol_Node_Positive'
cph_strata = CoxPHFitter()
cph_strata.fit(df_train2, 'Survival_Months','Status',formula =FORMULA,strata=['age_strata','Estrogen_Status','Progesterone_Status'])
cph_strata.check_assumptions(df_train2, show_plots=True)
#所有不符合的PH假定的变量都处理完了，提示：Proportional hazard assumption looks okay
Proportional hazard assumption looks okay.
[]
#最后打印一下处理完不符合PH假定的协变量的cox模型
#发现c-index是变低了。
cph_strata.print_summary()
()
cph_strata.plot()
| model | lifelines.CoxPHFitter |
|---|---|
| duration col | 'Survival_Months' |
| event col | 'Status' |
| strata | [age_strata, Estrogen_Status, Progesterone_Sta... |
| baseline estimation | breslow |
| number of observations | 3219 |
| number of events observed | 479 |
| partial log-likelihood | -2590.82 |
| time fit was run | 2022-02-23 17:46:56 UTC |
| coef | exp(coef) | se(coef) | coef lower 95% | coef upper 95% | exp(coef) lower 95% | exp(coef) upper 95% | z | p | -log2(p) | |
|---|---|---|---|---|---|---|---|---|---|---|
| Grade | 0.45 | 1.57 | 0.08 | 0.29 | 0.60 | 1.34 | 1.83 | 5.73 | <0.005 | 26.54 |
| N_Stage | 0.36 | 1.44 | 0.10 | 0.18 | 0.55 | 1.19 | 1.73 | 3.83 | <0.005 | 12.95 |
| Race | 0.42 | 1.53 | 0.12 | 0.20 | 0.65 | 1.22 | 1.92 | 3.64 | <0.005 | 11.84 |
| Reginol_Node_Positive | 0.05 | 1.05 | 0.01 | 0.03 | 0.08 | 1.03 | 1.08 | 4.14 | <0.005 | 14.79 |
| Regional_Node_Examined | -0.04 | 0.97 | 0.01 | -0.05 | -0.02 | 0.95 | 0.98 | -4.70 | <0.005 | 18.59 |
| T_Stage | 0.32 | 1.38 | 0.06 | 0.21 | 0.43 | 1.23 | 1.54 | 5.59 | <0.005 | 25.38 |
| Concordance | 0.70 |
|---|---|
| Partial AIC | 5193.64 |
| log-likelihood ratio test | 257.89 on 6 df |
| -log2(p) of ll-ratio test | 172.99 |
<AxesSubplot:xlabel='log(HR) (95% CI)'>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Kb7FI84k-1645703648882)(output_66_2.png)]
3.12 真的需要关心比例风险假设吗
lifelines 作者认为：
具体的详细可以查看原文链接哈
结语
对于熟悉python语言的人来说，Lifelines这生存分析工具的确很方便。但是对于熟悉生存分析的人来说，lifelines又有很多的不足。因为这个库缺少生成各种评价指标的工具，例如没有提供绘制ROC曲线，绘制诺模图，绘制决策曲线等工具。而R语言却提供很多包来绘制这些图。
还有感觉LifeLines这个库对多分类变量的处理不是很好，在R语言中，多分类变量一般会处理成因子型，会对每个子变量都以一个“基础”变量作对比，各生成风险比。而lifelines感觉会把多分类的变量处理成连续值变量。
毕竟做生存分析绝大部分都是学医的，图的就是现成的工具去解决科研中的问题。Lifelines没有提供的工具，用的人也一般不会特意造*。不过还是希望Lifelines以后的更新里，可以提供更多生存分析的工具。
参考资料
1.生存分析简明教程
2.生存分析（Survival Analysis）、Cox风险比例回归模型（Cox proportional hazards model）及C-index
5.生存分析
8.谭令,刘紫麟,唐翎瀚,肖江卫.直肠癌长期生存预测模型的构建与验证
—Nomogram预测模型[J].中国普外基础与临床杂志,2021,28(08):1031-1038.