用lifelines进行生存分析【Python版】

最新推荐文章于 2023-12-28 15:36:22 发布
最新推荐文章于 2023-12-28 15:36:22 发布
阅读量9.5k 收藏 174
点赞数 48
分类专栏： 生成分析 文章标签： python 机器学习
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_42990106/article/details/123119503
版权

开篇语

生存分析在医学研究中占有很大的比例，而且进行生存分析时，多用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男是1723存活
65女是2024存活
36女否1619死亡
38男是2521存活
69男否1623存活
70男否2820死亡

1.2 生存分析中常见专业术语

生存分析有一些常见的专业术语。例如什么是生存率，什么是生存曲线等。

1.3 生存分析的主要目的

  1. 估计研究对象的生存率
  2. 比较2不同组(影响因素)的生存率
  3. 分析影响研究对象生存期长短的因素有哪些和贡献度

1.4 生存分析的主要方法

  1. 寿命表法和Kaplan-Meier法（属于统计描述)
  2. Log-rank检验（属于统计推断)
  3. 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']
AgeRaceMarital_StatusT_StageN_StageSixth_StageGradeA_StageTumor_SizeEstrogen_StatusProgesterone_StatusRegional_Node_ExaminedReginol_Node_PositiveSurvival_MonthsStatus
04300224104000191110
1470021210450025920
267102012025004121
3461110010190026121
4631021210350021531

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~110100
1~21099
2~33198
3~41194
4~52092
5~38790

可以通过寿命表计算死亡概率和生存概率和生存率

死亡概率= 死亡人数/(初期人数-删失人数/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~1101000.010.990.99
1~210990.01010.98990.9799
2~331980.03070.9690.949
3~411940.0170.9890.938
4~520920.0220.9780.918
5~387900.06450.93550.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'])
removedobservedcensoredentranceat_risk
event_at
0.000040244024
1.071472404024
2.0114882603953
3.01171001703839
4.022511810703722
5.073910563403497
6.07256066502758
7.07315467702033
8.06743364101302
9.0628116170628
event_atremovedobservedcensoredentranceat_risk
0.000040244024
1.071472404024
2.0114882603953
3.01171001703839
4.022511810703722
5.073910563403497
6.07256066502758
7.07315467702033
8.06743364101302
9.0628116170628

通过以上的寿命表计算出死亡概率和生存概率和生存率

event_atremovedobservedcensoredentranceat_riskprobability of deathprobability of survivingsurvival rate
0.000040244024011
1.0714724040240.01170.98830.9883
2.01148826039530.02230.97770.9663
3.011710017038390.02610.97390.941
4.0225118107037220.032170.96780.911
5.0739105634034970.0330.9670.881
6.072560665027580.02470.9750.858
7.073154677020330.03190.96810.8313
8.067433641013020.03360.96640.803
9.06281161706280.03440.96560.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)
removedobservedcensoredentranceat_riskKM_estimateKM_estimate_lower_0.95KM_estimate_upper_0.95
0.0000402440241.0000001.0000001.000000
1.0101040241.0000001.0000001.000000
2.0321040230.9995030.9980140.999876
3.0422040200.9990060.9973530.999627
4.01091040160.9967670.9944380.998121
...........................
103.05005002510.7857090.7658720.804086
104.04804802010.7857090.7658720.804086
105.04504501530.7857090.7658720.804086
106.04704701080.7857090.7658720.804086
107.0610610610.7857090.7658720.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_distributionchi squared
degrees_of_freedom1
point_in_time60
fitterA<lifelines.KaplanMeierFitter:"KM_estimate", fi...
fitterB<lifelines.KaplanMeierFitter:"KM_estimate", fi...
test_namesurvival_difference_at_fixed_point_in_time_test
test_statisticp-log2(p)
040.84<0.00532.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_distributionchi squared
degrees_of_freedom0
test_namelogrank_test
test_statisticp-log2(p)
00.00nannan

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(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_distributionchi squared
degrees_of_freedom2
test_namemultivariate_logrank_test
test_statisticp-log2(p)
0298.86<0.005215.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是发生失效事件的数量

modellifelines.CoxPHFitter
duration col‘Survival_Months’
event col‘Status’
baseline estimationbreslow
number of observations3219
number of events observed479
partial log-likelihood-3516.57
time fit was run2022-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%。

coefexp(coef)se(coef)coef lower 95%coef upper 95%exp(coef) lower 95%exp(coef) upper 95%zp-log2§
Age0.021.020.010.010.031.011.033.09<0.0058.95
Race0.361.430.120.130.591.141.803.09<0.0058.96
Marital_Status0.071.080.04-0.010.160.991.171.760.083.69

表格三：Concordance 是一致性指数，用了全部协变量训练cox模型，一致性指数为0.75

Concordance0.75
Partial AIC7059.14
log-likelihood ratio test398.31 on 13 df
-log2§ of ll-ratio test253.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.
modellifelines.CoxPHFitter
duration col'Survival_Months'
event col'Status'
baseline estimationbreslow
number of observations3219
number of events observed479
partial log-likelihood-3516.57
time fit was run2022-02-23 17:30:27 UTC
coefexp(coef)se(coef)coef lower 95%coef upper 95%exp(coef) lower 95%exp(coef) upper 95%zp-log2(p)
Age0.021.020.010.010.031.011.033.09<0.0058.95
Race0.361.430.120.130.591.141.803.09<0.0058.96
Marital_Status0.071.080.04-0.010.160.991.171.760.083.69
T_Stage0.371.450.110.160.581.171.793.42<0.00510.63
N_Stage0.361.430.180.020.701.022.022.050.044.63
Sixth_Stage-0.001.000.11-0.230.220.801.25-0.010.990.01
Grade0.471.600.080.320.631.371.876.02<0.00529.06
A_Stage0.001.000.21-0.410.410.661.510.010.990.01
Tumor_Size-0.001.000.00-0.010.000.991.00-0.640.520.93
Estrogen_Status0.541.720.150.240.841.282.313.58<0.00511.50
Progesterone_Status0.551.730.120.310.781.372.184.58<0.00517.71
Regional_Node_Examined-0.030.970.01-0.05-0.020.950.98-4.69<0.00518.49
Reginol_Node_Positive0.051.050.010.030.081.031.084.10<0.00514.58

Concordance0.75
Partial AIC7059.14
log-likelihood ratio test398.31 on 13 df
-log2(p) of ll-ratio test253.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()
modellifelines.CoxPHFitter
duration col'Survival_Months'
event col'Status'
baseline estimationbreslow
number of observations3219
number of events observed479
partial log-likelihood-3518.31
time fit was run2022-02-23 17:30:49 UTC
coefexp(coef)se(coef)coef lower 95%coef upper 95%exp(coef) lower 95%exp(coef) upper 95%zp-log2(p)
Age0.021.020.010.010.031.011.033.32<0.00510.13
Estrogen_Status0.551.730.150.250.841.292.333.65<0.00511.87
Grade0.471.610.080.320.631.381.876.06<0.00529.43
N_Stage0.351.420.090.160.531.181.713.68<0.00512.05
Progesterone_Status0.551.730.120.310.781.372.184.57<0.00517.68
Race0.401.490.120.170.621.191.873.45<0.00510.78
Reginol_Node_Positive0.051.050.010.030.081.031.084.31<0.00515.89
Regional_Node_Examined-0.030.970.01-0.05-0.020.950.98-4.68<0.00518.42
T_Stage0.331.390.060.220.441.241.565.80<0.00527.13

Concordance0.75
Partial AIC7054.62
log-likelihood ratio test394.83 on 9 df
-log2(p) of ll-ratio test261.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.01.000000
2.00.999882
3.00.999647
4.00.998941
5.00.998587
......
103.00.901906
104.00.901906
105.00.901906
106.00.901906
107.00.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_statisticp-log2§
A_Stagekm4.910.035.23
rank4.430.044.83
Estrogen_Statuskm4.990.035.29
rank6.110.016.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假定如何解决】

  1. 方法一：【分层】对这个不满足PH假定的变量进行分层，对剩余的协变量进行cox建模

  2. 方法二：【时依变量协变量】

假如有不满足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_distributionchi squared
degrees_of_freedom1
model<lifelines.CoxPHFitter: fitted with 3219 total...
test_nameproportional_hazard_test
test_statisticp-log2(p)
Agekm3.220.073.78
rank4.690.035.04
Estrogen_Statuskm4.990.035.29
rank6.110.016.22
Gradekm0.030.860.22
rank0.010.910.13
N_Stagekm1.310.251.99
rank1.260.261.93
Progesterone_Statuskm10.13<0.0059.42
rank9.77<0.0059.14
Racekm0.220.640.65
rank0.070.790.34
Reginol_Node_Positivekm0.850.361.49
rank0.840.361.48
Regional_Node_Examinedkm0.000.950.07
rank0.040.840.25
T_Stagekm0.000.950.07
rank0.030.860.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_distributionchi squared
degrees_of_freedom1
model<lifelines.CoxPHFitter: fitted with 3219 total...
test_nameproportional_hazard_test
test_statisticp-log2(p)
Agekm3.270.073.82
rank6.040.016.16
Gradekm0.020.890.18
rank0.890.351.54
N_Stagekm1.410.242.09
rank0.050.820.28
Racekm0.390.530.91
rank0.720.401.34
Reginol_Node_Positivekm1.080.301.74
rank1.060.301.72
Regional_Node_Examinedkm0.000.970.04
rank0.340.560.84
T_Stagekm0.010.910.14
rank0.010.940.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()
modellifelines.CoxPHFitter
duration col'Survival_Months'
event col'Status'
strata[age_strata, Estrogen_Status, Progesterone_Sta...
baseline estimationbreslow
number of observations3219
number of events observed479
partial log-likelihood-2590.82
time fit was run2022-02-23 17:46:56 UTC
coefexp(coef)se(coef)coef lower 95%coef upper 95%exp(coef) lower 95%exp(coef) upper 95%zp-log2(p)
Grade0.451.570.080.290.601.341.835.73<0.00526.54
N_Stage0.361.440.100.180.551.191.733.83<0.00512.95
Race0.421.530.120.200.651.221.923.64<0.00511.84
Reginol_Node_Positive0.051.050.010.030.081.031.084.14<0.00514.79
Regional_Node_Examined-0.040.970.01-0.05-0.020.950.98-4.70<0.00518.59
T_Stage0.321.380.060.210.431.231.545.59<0.00525.38

Concordance0.70
Partial AIC5193.64
log-likelihood ratio test257.89 on 6 df
-log2(p) of ll-ratio test172.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

3.生存分析不止中位生存时间

4.临床研究中常用的评价指标AUC和C-index

5.生存分析

6.怎样理解生存分析中的风险函数

7.生存分析的COX回归

8.教你三招 Cox回归比例风险(PH)假定的检验

8.谭令,刘紫麟,唐翎瀚,肖江卫.直肠癌长期生存预测模型的构建与验证
—Nomogram预测模型[J].中国普外基础与临床杂志,2021,28(08):1031-1038.

python生存分析入门
Andy_shenzl的博客
06-03 5431
引入 什么是生存分析 生存分析是对生存资料的分析。所谓生存资料是指描述寿命或者一个发生时间的数据。通过对某一具有相同特质的群体进行数据分析，我们可以得到这类人活过一定时间的概率。这就叫做生存分析。一个人的生存时间长短与许多因素有关，研究各个因素与生存时间有无关系以及关联程度大小，这也是生存分析。推广来说，疾病的复发可以看做“健康期”的生存，机器的故障可以看做“机器”的生存，甚至下岗职工再就业可以看做“失业期”的生存，生存分析可以被推广至很多问题的分析中去。 应用场景 生存分析最早在生物医学中使用的最多，用来
生命线：Python中的生存分析
02-05
生存分析最初是由精算和医学界开发并广泛应用的。 其目的是回答为什么现在发生事件而不是后来在不确定性下发生（事件可能指死亡，疾病缓解等）。 这对于对测量寿命感兴趣的研究人员非常有用：他们可以回答诸如哪些因素可能影响死亡的问题？ 但是，除了医学和精算科学外，生存分析还有许多其他有趣而激动人心的应用。 例如： SaaS提供商对衡量订户的生存期或采取某些第一行动的时间感兴趣 库存缺货是对商品真正“需求”的审查事件。 社会学家对衡量政党的一生，人际关系或婚姻感兴趣 A / B测试可确定不同组执行一项操作需要多长时间。 lifelines是生存分析的最佳部分的纯Python实现。 记录和生存分析简
python实现生存分析（Survival Analysis）
qq_39579290的博客
09-15 5568
python实现KM生存分析
生存分析序章2——生存分析之Python篇：lifelines库入门
theskylife的博客
12-26 3096
数据分析和挖掘在当今信息时代扮演着至关重要的角色。生存分析作为其中一项重要的统计学方法，被广泛应用于医学、工程、金融等领域，帮助我们理解事件发生与时间的关系。在本文中，我们将深入介绍Python中的lifelines库，探讨其在生存分析中的应用。
python做生存分析
qq_26108393的博客
10-08 5362
这里跟大家分享一下生存分析，并且介绍一个用python做生存分析的包。 生存分析 应用场景 生存分析最早在生物医学中使用的最多，用来预估某个群体的存活时间，后来被推广到了更广泛的领域。一个做生物的专家可能通常会关心这样一个问题：这个群体的样本能活多久？这个问题我们通常会使用生存分析来回答。这个群体可以理解为某个国家的人民，或者注射过某种药物的一些病人。 同样也可以推广到更一般的场景，公司...
生存分析随机森林实验与代码_python_生存分析_随机森林_
10-03
在Python中，我们可以使用`lifelines`库进行生存分析，这是一个强大的生存分析库，支持多种生存模型。对于随机森林，可以使用`sklearn.ensemble.RandomForestClassifier`或`sklearn.ensemble.RandomForestRegressor`...
cox回归模型python实现_Python用Cox比例风险模型进行生存分析-基于lifelines库
weixin_39588445的博客
12-06 5066
由于需要对数据进行生存分析，我之前都是在R语言上面进行的，这次想在python上实现，在网上找了好多教程，为了以后忘记，在此记录一下Cox模型的详细理论这里就不讨论了，网上有很多介绍的，自行搜索，这里直接上代码吧首先我们需要在cmd下安装一个lifelines库pip install lifelines接下来，在python中导入相关的模块，这里使用官方的数据集进行演示#导入相关的库from li...
python中e库是干什么的_用Python做生存分析--lifelines库简介
weixin_39687359的博客
11-24 892
Python提供了一个简单而强大的生存分析包——lifelines，可以非常方便的进行应用。这篇文章将为大家简单介绍这个包的安装和使用。安装lifelines支持用pip的方法进行安装，您可以使用以下命令进行一键安装：pip install lifelines使用简介在python中，可以利用lifelines进行累计生存曲线的绘制、Log Rank test、Cox回归等。下面以lifeline...
python数据分析案例-利用生存分析Kaplan-Meier法与COX比例风险回归模型进行客户流失分析与剩余价值预测
热门推荐
maiyida123的博客
08-16 1万+
利用Kaplan-Meier法与COX比例风险回归模型比较了不同特征对于客户流失的影响，并且预测了客户的生存时长，计算了其剩余价值以及潜在可提升的价值。
python 生存分析,利用python进行泰坦尼克生存预测——数据探索分析
weixin_42327291的博客
03-25 470
最近一直断断续续的做这个泰坦尼克生存预测模型的练习，这个kaggle的竞赛题，网上有很多人都分享过，而且都很成熟，也有些写的非常详细，我主要是在牛人们的基础上，按照数据挖掘流程梳理思路，然后通过练习每一步来熟悉应用python进行数据挖掘的方式。数据挖掘的一般过程是：数据预览——>数据预处理(缺失值、离散值等)——>变量转换(构造新的衍生变量)——>数据探索(提取特征)——&gt...
Python库 | lifelines-0.16.2.tar.gz
04-10
资源分类：Python库 所属语言：Python 资源全名：lifelines-0.16.2.tar.gz 资源来源：官方 安装方法：https://lanzao.blog.csdn.net/article/details/101784059
Python库 | lifelines-0.7.1.0.tar.gz
05-17
资源分类：Python库 所属语言：Python 资源全名：lifelines-0.7.1.0.tar.gz 资源来源：官方 安装方法：https://lanzao.blog.csdn.net/article/details/101784059
一篇项目走进生存分析(Survival Analysis)的世界【Python版
m0_63642362的博客
03-05 6113
使用Python版生存分析工具LifeLines，对乳腺癌数据进行生存分析。
Python绘图——生存曲线
aobulaien001的博客
09-09 2089
以上就是本次推送的全部内容，希望小伙伴们能够多提提建议！
生存分析序章3——生存分析之Python篇：介绍生存分析和scikit-survival库
最新发布
theskylife的博客
12-28 1744
在现代数据科学领域，生存分析是一项关键的统计方法，用于研究个体在一定时间内事件发生的概率。无论是医学、金融还是社会科学，生存分析都为我们提供了一种深入洞察数据的途径。在Python生态系统中，有一个强大的工具——scikit-survival库，专门用于处理生存分析问题。本篇博客将为您介绍生存分析的基本概念以及如何使用scikit-survival库进行生存分析，结合实际应用场景，揭示数据分析和数据挖掘中的潜在价值。
Survival Analysis with Python 生存分析及Python实现
百载文枢江左的博客
04-03 1098
survival analysis with Python
Python实战|利用生存分析预测用户流失周期（一）
weixin_41432734的博客
07-12 956
点击上方蓝字关注我们利用生存分析预测用户流失周期 上篇 ...
生存分析及其在Python中的应用
HackNebulaZ的博客
09-09 128
累积风险描述了在给定时间点之前发生事件的累积概率，而生存曲线是生存函数的补数，表示个体在给定时间点之前发生事件的概率。累积风险描述了在给定时间点之前发生事件的累积概率，而生存曲线是生存函数的补数，表示个体在给定时间点之前发生事件的概率。生存分析中的一个重要概念是生存函数（Survival Function），它描述了个体在给定时间段内存活下来的概率。生存分析中的一个重要概念是生存函数（Survival Function），它描述了个体在给定时间段内存活下来的概率。首先，我们创建一个示例数据集，其中。
数据分析学习总结笔记13：生存分析及Python实现
weixin_41961559的博客
04-18 7824
文章目录1 引言2 定义3 数学直观4 Kaplan-Meier 估计5 Cox比例风险模型6 总结 1 引言 生存分析是一套统计方法，用来解决诸如“多长时间后，某个特定事件发生”这样的问题; 换句话说，也可以称之为**“事件时间分析”**。 这种方法被称为生存分析，是由于主要是由医学研究人员开发的，他们更感兴趣的是寻找不同群体患者的预期寿命(例如: 用药物a治疗的组群1和用药物b治疗的组群2)。...
生存分析python
09-16
生存分析是一种统计方法，用于分析个体在研究期间的生存时间或事件发生的概率。在Python中，我们可以使用不同的库来进行生存分析，如lifelines、scikit-survival和survival等。 在lifelines库中，可以使用Kaplan-Meier估计方法来估计生存曲线，并使用Cox比例风险回归模型来评估预测因素的影响。在scikit-survival库中，可以使用Cox比例风险回归模型、Aalen's加性模型和Kaplan-Meier估计方法来进行生存分析。同时，在survival库中，也提供了各种生存分析的方法和模型。 下面是一个使用lifelines库进行生存分析的示例代码： ```python # 导入库 from lifelines import KaplanMeierFitter from lifelines.datasets import load_waltons # 加载数据集 data = load_waltons() # 创建Kaplan-Meier估计器对象 kmf = KaplanMeierFitter() # 根据不同的组别计算生存曲线 groups = data['group'] time = data['time'] event = data['event'] kmf.fit(time, event, label='All') # 绘制生存曲线 kmf.plot() # 显示图形 plt.show() ```
评论 12
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值