一篇项目走进生存分析(Survival Analysis)的世界【Python版

本文介绍了使用Python的lifelines库进行生存分析的项目。首先阐述生存分析的概念、数据特点、术语、目的及方法，接着以乳腺ai数据集为例，展示数据处理过程，然后详细演示寿命表法、Kaplan-Meier法、Log-rank检验及Cox比例风险回归模型的应用，包括模型构建、检验、评估及协变量处理等，为医学生存分析提供参考。

开篇语

生存分析在医学研究中占有很大的比例，而且进行生存分析时，多用R语言、SPSS等工具进行生存分析，用python进行生存分析不多。因为发现一个python版的生存分析工具---lifelines ，这个库已经提供比较完善的生存分析相关的工具。自己又最近学习生存分析，然后结合lifelines开始编写这个项目。写代码的同时，也对一些生存分析中概念性的名词，根据自己的理解一起展示出来。因为是边学边写，有错误的地方请指正 。

In [ ]

#安装生存分析用的python库----lifelines#lifelines相关文档地址https://lifelines.readthedocs.io/en/latest/!pip install lifelines!pip install  numpy==1.19.2!pip install pyzmq==18.1.1

登录后复制In [2]

import pandas as pdimport numpy as npimport seaborn as snsimport matplotlib.pyplot as plt

登录后复制

1.什么是生存分析

在医学中，生存分析是非常常见和重要的分析方法之一,是一种将事件的结果和出现这种结果所经历的时间结合起来一起分析的统计方法。

案例一： 为了研究什么因素会影响直肠ai患者生存时间。会收集患者的资料。例如有年龄、性别、婚姻状态、TMN分期、T分期、肿瘤直径、肿瘤分化程度、治疗手段等。还有会收集患者在一定时间内发生死亡事件的时间，即生存时间。采用Kaplan–Meier 曲线单变量进行分析，那个因素是影响直肠ai患者的生存时间，或者采用Cox 比例风险回归进行多因素分析。并对新患有直肠ai患者进行生存时间预测。

立即学习“Python免费学习笔记（深入）”；

案例二：为了对某种抗ai药物做临床试验，将ai症患者随机分成两组，一组服用该抗ai药物，另一组服用对照药。并记录两组患者开始服药直到死亡的生存时间。并分析两组中死亡事件发生与时间关系是否存在差异，来判断药物是否有效。

通过两个案例可以看出生存分析就是研究在不同条件下，事件发生与时间的关系 ，不仅考虑事件是否出现，还有考虑事件出现的时间长短。

1.1生存分析的数据是怎样的

生存分析的数据和分类数据有点类似。但是生存分析多了一列关于时间数据，这列时间数据是从特定事件发生了，记录整个观察期间研究对象多久才发生失效事件(例如死亡事件）。生存分析就是研究生存时间和结局与众多影响因素间的关系。

例如下面表格记录了6位ai症患者在2年观测期间中的存活时间和感兴趣的影响因素。这里的失效事件是死亡事件，是我们关心的。生存时间对应表格的生存时间。研究的影响因素有：年龄、性别、是否吸烟、肿瘤大小。

在表格中有个患者生存时间只有21个月，结局却是存活，可能由于是患者中途退出观察或者换了电话号码，不能联系患者，无法继续进行研究，这种数据属于删失数据，可见在现实生活中生存分析数据通常会含有删失数据。再例如在表格中有个患者生存时间是19个月，结局是死亡，这种数据属于完整数据。不过这些删失数据在进行生存分析时，也是会一起进行分析。

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 查看本项目中乳腺ai数据的基线资料表

这个乳腺ai数据集研究的失效事件是"Status"中的"Dead",生存时间是"Survival_Months"单位是月。

这个乳腺ai数据集中具有的特征如下：

Age(年龄)Race(种族）Marital_Status(婚姻状态)T_Stage(原发肿瘤的大小)N_Stage(ai细胞扩散的淋巴结情况)Sixth_StageGradeA_StageTumor_Size(肿瘤大小)Estrogen_Status(雌激素受体是否阳性)Progesterone_Status(孕激素受体是否阳性)Regional_Node_ExaminedReginol_Node_PositiveSurvival_Months（生存时间，单位:月)

通过基线资料表中根据“Alive”和“Dead”分成两组，可以看出每个特征的人数分布和比例。 

2.2 读取、处理数据

通过pandas读取csv文件获取数据。通过head()可以看到某些列是由字符串组成的，例如“Race”列中每个值由“Black”、“White”等值组成，所以进行生存分析时需要编码成数字0、1、2等。

In [28]

"""读取数据并处理分类变量数据，不然使用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:RaceValue:['Other (American Indian/AK Native, Asian/Pacific Islander)', 'White', 'Black']Column name:Marital_StatusValue:['Married (including common law)', 'Divorced', 'Single (never married)', 'Widowed', 'Separated']Column name:T_StageValue:['T2', 'T1', 'T3', 'T4']Column name:N_StageValue:['N3', 'N2', 'N1']Column name:Sixth_StageValue:['IIIC', 'IIIA', 'IIB', 'IIA', 'IIIB']Column name:GradeValue:['Moderately differentiated; Grade II', 'Poorly differentiated; Grade III', 'Well differentiated; Grade I', 'Undifferentiated; anaplastic; Grade IV']Column name:A_StageValue:['Regional', 'Distant']Column name:Estrogen_StatusValue:['Positive', 'Negative']Column name:Progesterone_StatusValue:['Positive', 'Negative']Column name:StatusValue:['Alive', 'Dead']

登录后复制

   Age  Race  Marital_Status  T_Stage  N_Stage  Sixth_Stage  Grade  A_Stage  \0   43     0               0        2        2            4      1        0   1   47     0               0        2        1            2      1        0   2   67     1               0        2        0            1      2        0   3   46     1               1        1        0            0      1        0   4   63     1               0        2        1            2      1        0      Tumor_Size  Estrogen_Status  Progesterone_Status  Regional_Node_Examined  \0          40                0                    0                      19   1          45                0                    0                      25   2          25                0                    0                       4   3          19                0                    0                      26   4          35                0                    0                      21      Reginol_Node_Positive  Survival_Months  Status  0                     11                1       0  1                      9                2       0  2                      1                2       1  3                      1                2       1  4                      5                3       1

登录后复制

2.3 简单数据分析

对数据绘制柱状图、箱线图等，展示某些数据与终末事件的关系

In [4]

"""一共有4024个研究对象，发生删失3408例，发生死亡事件616例"""print(len(df))print(df['Status'].value_counts(),'\n')

登录后复制

40240    34081     616Name: Status, dtype: int64

登录后复制In [5]

#不同T(肿瘤大小)分期，删失人数和死亡人数的柱状分布图sns.countplot(x="T_Stage",hue="Status",data=df,)

登录后复制

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [6]

#不同淋巴结转移情况，删失人数和死亡人数的柱状分布图sns.countplot(x="N_Stage",hue="Status",data=df,)

登录后复制

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [7]

"""绘制发生删失组和发生死亡组的年龄箱线图发生死亡事件组的年龄中位数比删失组高一点"""ax = df.boxplot(column='Age',by="Status",figsize=(6,6))ax.set_xticklabels(["censor","death"])

登录后复制

[Text(1, 0, 'censor'), Text(2, 0, 'death')]

登录后复制

登录后复制

2.4 绘制生命周期图

在数据集中随机挑选30例患者，绘制生命周期图。假设观测时间定位96个月，在96个月中，有5位患者发生失效事件(死亡),其余患者大多数还没到96个月就失去了联系/无法随访到患者(这种数据属于删失数据）。例如图中1347号患者，在随访到第62个月，患者就失去了联系，并不知道患者后面的真实存活时间。

通过生命周期图，可以看到出现很多数据都是删失数据，无论不结合删失数据，单独计算发生失效事件的患者的平均存活时间，还是结合删失数据一起一起计算患者的平均存活时间，这样计算都是要低估了真实的存活时间，可以认识到了右删失数据对生存时间估计的影响很大，因此生存分析就是解决此类问题。 

In [ ]

from lifelines.plotting import plot_lifetimesdf_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年）的随访资料。因为某些原因，对于研究对象并不是每天进行联系，而是在某个时间区间，例如一年后获取对象的资料。因此并不是准确记录结局的发生时间，只知道某个区间中发生了终点事件。这种情况可以绘制寿命表计算生存率进行分析。

假设有寿命表(类似频数表)如下：

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

死亡概率= 死亡人数/(初期人数-删失人数/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

就可以通过计算出来的生存率绘制生存曲线

进行单因素分析时，可以根据影响因素进行分组，分别制作寿命表，进行计算各个指标并绘制生存曲线进行比较。

3.2 对乳腺ai数据集生成寿命表

计算生存率并绘制生存曲线

In [9]

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_riskevent_at                                                0.0             0         0         0      4024     40241.0            71        47        24         0     40242.0           114        88        26         0     39533.0           117       100        17         0     38394.0           225       118       107         0     37225.0           739       105       634         0     34976.0           725        60       665         0     27587.0           731        54       677         0     20338.0           674        33       641         0     13029.0           628        11       617         0      628

登录后复制

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

In [10]

#然后通过计算出来的生存率绘制生存曲线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()

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制

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是每一年的生存概率。

原始的乳腺ai数据集的生存时间是以月位单位的，样本量有4000多个，如果通过寿命表来绘制，工作量很大。现在使用lifelines库中的KaplanMeierFitter工具就可以使用K-M法绘制生存曲线。

通过分组绘制生存曲线，展示不同因素下的研究对象的生存曲线，从而进行单因素分析。

通过KM法还能计算中位生存时间，中位生存时间就是纵坐标=0.5（即累积生存率=0.5的时候所对应的时间t），当绘制多组生存曲线的时候可以通过中位生存时间简单比较组与组之间的差异。

 例如这个图，蓝色曲线的生存时间大概1200，黄色曲线的2100，黄色曲线的存活时间比蓝色的长。

In [11]

#使用KM分析，并绘制生存曲线和计算中位生存时间from lifelines import KaplanMeierFitterfrom lifelines.utils import median_survival_timeskmf = 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_estimatetimeline             0.0          1.0000001.0          1.0000002.0          0.9995033.0          0.9990064.0          0.996767...               ...103.0        0.785709104.0        0.785709105.0        0.785709106.0        0.785709107.0        0.785709[108 rows x 1 columns]中间存活概率=0.5时的，95%区间     KM_estimate_lower_0.95  KM_estimate_upper_0.950.5                     inf                     inf

登录后复制

登录后复制登录后复制

登录后复制In [12]

#打印生存曲线的简短统计内容(寿命表)。包括观察值，事件数，删失数，生存率和其置信区间。pd.concat([kmf.event_table,kmf.survival_function_,kmf.confidence_interval_],axis=1)

登录后复制

       removed  observed  censored  entrance  at_risk  KM_estimate  \0.0          0         0         0      4024     4024     1.000000   1.0          1         0         1         0     4024     1.000000   2.0          3         2         1         0     4023     0.999503   3.0          4         2         2         0     4020     0.999006   4.0         10         9         1         0     4016     0.996767   ...        ...       ...       ...       ...      ...          ...   103.0       50         0        50         0      251     0.785709   104.0       48         0        48         0      201     0.785709   105.0       45         0        45         0      153     0.785709   106.0       47         0        47         0      108     0.785709   107.0       61         0        61         0       61     0.785709          KM_estimate_lower_0.95  KM_estimate_upper_0.95  0.0                  1.000000                1.000000  1.0                  1.000000                1.000000  2.0                  0.998014                0.999876  3.0                  0.997353                0.999627  4.0                  0.994438                0.998121  ...                       ...                     ...  103.0                0.765872                0.804086  104.0                0.765872                0.804086  105.0                0.765872                0.804086  106.0                0.765872                0.804086  107.0                0.765872                0.804086  [108 rows x 8 columns]

登录后复制

3.4 对数据进行分组，用KM法进行单因素生存分析

使用KM法可以分析整个数据的研究对象的生存率。有的时候为了研究药物是否对疾病带来治疗功效，会设置实验组和对照组，使用KM分析法时可以分别画出两组的生存曲线来展示两组的差别，或者计算中位生存时间，假设实验组的中位生存时间大于对照组，就说明药物可以延迟患者的寿命。

但是KM法只能针对分类变量进行分组，对应连续值需要把连续值区间化，再进行分组.

比较两组之间的生存差异可以通过中位生存时间来比较。也可以指定时间点来计算P值两组曲线是否有显著性差异。也可以计算T时间内生存曲线下的面积，不同生存曲线之间的差异可以通过曲线下面积来比较（限制性平均生存时间/RMST）。也可以通过计算第p百分位生存时间来比较两组之间的生存差异

In [13]

"""对T——Stage 肿瘤大小进行分组并绘制生存曲线从生存曲线中明显看到T4的患者（肿瘤大)的生存率下降比较快"""from lifelines import KaplanMeierFitterplt.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()

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制In [14]

"""因为年龄是连续值，km法只适用分类变量，因为需要把连续值的年龄进行区间分组对年龄进行分组，并绘制生存曲线"""from lifelines import KaplanMeierFitterdf['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()

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制In [42]

"""可以通过中位生存时间来比较两组之间的生存差异（因为这个数据有些组的生存率最后都是大于0.5，无法计算中卫生存时间，所以显示inf）"""from lifelines.utils import median_survival_timeskmf = 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组的中位生存时间:infT2组的中位生存时间:infT3组的中位生存时间:infT4组的中位生存时间:inf

登录后复制In [43]

"""有些数据到随访结束时生存率都是大于0.5，即无法计算中位生存时间。但是可以指定第qth分位数的生存时间，来求得对应的生存时间"""from lifelines.utils import qth_survival_timeqth = 0.9kmf = 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.0T2组的第90.0分位数的生存时间:49.0T3组的第90.0分位数的生存时间:40.0T4组的第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,

登录后复制In [17]

"""可以通过指定时间来比较这个时间对应的生存差异（例如我关心的时第60个月的生存差异)"""from lifelines.statistics import survival_difference_at_fixed_point_in_time_testfrom lifelines.datasets import load_waltonsT = df['Survival_Months']E = df['Status']ix = df['A_Stage']== 0kmf_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 squareddegrees_of_freedom = 1     point_in_time = 60.0           fitterA =            fitterB =          test_name = survival_difference_at_fixed_point_in_time_test--- test_statistic      p  -log2(p)          40.84 <0.005     32.49

登录后复制抱歉，数学公式解析失败，可能是拼写错误或暂不支持。In [20]

"""限制性平均生存时间/RMST通过计算不同生存曲线下面积来比较生存差异"""from lifelines.utils import restricted_mean_survival_timefrom lifelines.plotting import rmst_plotT = 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.22662024019032曲线下面积:49.82516339891577

登录后复制

登录后复制登录后复制

登录后复制

3.5 Log-rank检验

在使用KM方法，根据分组绘制多条生存曲线后，只通过直接的观察来确定多条曲线之间是否具有显著性差异是不充分的。Log-rank检验弥补这方面。

Log-rank检验的零假设组与组的生存时间分布没有显著性差异。当计算出来的P值小于0.05,推翻零假设，认为两组的生存时间分布存在显著性差异。

In [21]

"""比较两个组的生存曲线是否有显著性差异结果A_Stage组的P值小于0.005，是有显著性差异(直接观察图形也看出明显差异)"""from lifelines.statistics import logrank_testfrom lifelines import KaplanMeierFitterplt.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 squareddegrees_of_freedom = 0         test_name = logrank_test--- test_statistic    p  -log2(p)           0.00  nan       nan

登录后复制抱歉，数学公式解析失败，可能是拼写错误或暂不支持。

登录后复制登录后复制登录后复制登录后复制登录后复制In [22]

"""比较两个组以上的生存曲线是否有显著性差异结果T_Stage组的P值小于0.005，是有显著性差异(直接观察图形也看出明显差异)"""from lifelines.statistics import multivariate_logrank_testkmf = 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 squareddegrees_of_freedom = 2         test_name = multivariate_logrank_test--- test_statistic      p  -log2(p)         298.86 <0.005    215.58

登录后复制抱歉，数学公式解析失败，可能是拼写错误或暂不支持。

1.2701885843852175e-65

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

3.6 参数估计

如果事先知道生存时间的分布，可以使用参数法：根据样本观测值来估计假定的分布模型中的参数，获得生存时间的概率分布模型。 生存时间经常服从的分布有：指数分布、Weibull分布、对数正态分布、对数Logistic分布、Gamma分布。分析过程类似KM法，由于往往对真实数据的生存时间分布并不知道，所有参数估计并不是常用的生存分析方法。

In [23]

"""简单展示如果用linelifes来使用参数法进行生存分析并绘制生存曲线。"""from lifelines.datasets import load_waltonsfrom lifelines.utils import median_survival_timesimport matplotlib.pyplot as pltimport numpy as npfrom 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])

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

登录后复制

3.7 比例风险回归模型--Cox模型

Cox比例风险回归模型是英国统计学家D.R.Cox（很遗憾，2024年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对数据进行分割成训练集和验证集

In [33]

#加载数据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=2024)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    27401     479Name: Status, dtype: int640    6681    137Name: Status, dtype: int64

登录后复制

3.9 建立Cox模型并开始训练

通过print_summary()模型会输出三个表格

表格一：duration col是时间列名，event col是事件列名， breslow是模型参数估计方法，3219是总数据量，479是发生失效事件的数量

表格二：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%。

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

In [34]

"""单变量分析"""from lifelines import CoxPHFittervarList = 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 TrueRace 的P值：1.3550638342845673e-06,<0.05 TrueMarital_Status 的P值：8.68912013912815e-05,<0.05 TrueT_Stage 的P值：2.2304661752517098e-20,<0.05 TrueN_Stage 的P值：1.6617003098137667e-45,<0.05 TrueSixth_Stage 的P值：1.0290156463874445e-47,<0.05 TrueGrade 的P值：6.022575977628071e-23,<0.05 TrueA_Stage 的P值：1.0633163182382273e-08,<0.05 TrueTumor_Size 的P值：5.79647504244188e-16,<0.05 TrueEstrogen_Status 的P值：1.0415072293170037e-27,<0.05 TrueProgesterone_Status 的P值：5.007000090288116e-25,<0.05 TrueRegional_Node_Examined 的P值：0.1065243965415506,<0.05 FalseReginol_Node_Positive 的P值：1.5269471157224244e-46,<0.05 True

登录后复制In [35]

"""多变量分析"""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.2Iteration 2: norm_delta = 0.47681, step_size = 0.9000, log_lik = -3575.51461, newton_decrement = 52.18325, seconds_since_start = 0.4Iteration 3: norm_delta = 0.12309, step_size = 0.9000, log_lik = -3520.88106, newton_decrement = 3.96398, seconds_since_start = 0.6Iteration 4: norm_delta = 0.01909, step_size = 1.0000, log_lik = -3516.63264, newton_decrement = 0.06359, seconds_since_start = 0.9Iteration 5: norm_delta = 0.00045, step_size = 1.0000, log_lik = -3516.56817, newton_decrement = 0.00003, seconds_since_start = 1.1Iteration 6: norm_delta = 0.00000, step_size = 1.0000, log_lik = -3516.56814, newton_decrement = 0.00000, seconds_since_start = 1.3Convergence success after 6 iterations.

登录后复制

             duration col = 'Survival_Months'                event col = 'Status'      baseline estimation = breslow   number of observations = 3219number of events observed = 479   partial log-likelihood = -3516.57         time fit was run = 2024-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%covariate                                                                                                                      Age                      0.02       1.02       0.01             0.01             0.03                 1.01                 1.03Race                     0.36       1.43       0.12             0.13             0.59                 1.14                 1.80Marital_Status           0.07       1.08       0.04            -0.01             0.16                 0.99                 1.17T_Stage                  0.37       1.45       0.11             0.16             0.58                 1.17                 1.79N_Stage                  0.36       1.43       0.18             0.02             0.70                 1.02                 2.02Sixth_Stage             -0.00       1.00       0.11            -0.23             0.22                 0.80                 1.25Grade                    0.47       1.60       0.08             0.32             0.63                 1.37                 1.87A_Stage                  0.00       1.00       0.21            -0.41             0.41                 0.66                 1.51Tumor_Size              -0.00       1.00       0.00            -0.01             0.00                 0.99                 1.00Estrogen_Status          0.54       1.72       0.15             0.24             0.84                 1.28                 2.31Progesterone_Status      0.55       1.73       0.12             0.31             0.78                 1.37                 2.18Regional_Node_Examined  -0.03       0.97       0.01            -0.05            -0.02                 0.95                 0.98Reginol_Node_Positive    0.05       1.05       0.01             0.03             0.08                 1.03                 1.08                           z      p   -log2(p)covariate                                     Age                     3.09 <0.005       8.95Race                    3.09 <0.005       8.96Marital_Status          1.76   0.08       3.69T_Stage                 3.42 <0.005      10.63N_Stage                 2.05   0.04       4.63Sixth_Stage            -0.01   0.99       0.01Grade                   6.02 <0.005      29.06A_Stage                 0.01   0.99       0.01Tumor_Size             -0.64   0.52       0.93Estrogen_Status         3.58 <0.005      11.50Progesterone_Status     4.58 <0.005      17.71Regional_Node_Examined -4.69 <0.005      18.49Reginol_Node_Positive   4.10 <0.005      14.58---Concordance = 0.75Partial AIC = 7059.14log-likelihood ratio test = 398.31 on 13 df-log2(p) of ll-ratio test = 253.44

登录后复制抱歉，数学公式解析失败，可能是拼写错误或暂不支持。In [36]

#画出每个协变量的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'

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制In [37]

#去掉不显著的变量重新拟合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()

登录后复制

             duration col = 'Survival_Months'                event col = 'Status'      baseline estimation = breslow   number of observations = 3219number of events observed = 479   partial log-likelihood = -3518.31         time fit was run = 2024-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%covariate                                                                                                                      Age                      0.02       1.02       0.01             0.01             0.03                 1.01                 1.03Estrogen_Status          0.55       1.73       0.15             0.25             0.84                 1.29                 2.33Grade                    0.47       1.61       0.08             0.32             0.63                 1.38                 1.87N_Stage                  0.35       1.42       0.09             0.16             0.53                 1.18                 1.71Progesterone_Status      0.55       1.73       0.12             0.31             0.78                 1.37                 2.18Race                     0.40       1.49       0.12             0.17             0.62                 1.19                 1.87Reginol_Node_Positive    0.05       1.05       0.01             0.03             0.08                 1.03                 1.08Regional_Node_Examined  -0.03       0.97       0.01            -0.05            -0.02                 0.95                 0.98T_Stage                  0.33       1.39       0.06             0.22             0.44                 1.24                 1.56                           z      p   -log2(p)covariate                                     Age                     3.32 <0.005      10.13Estrogen_Status         3.65 <0.005      11.87Grade                   6.06 <0.005      29.43N_Stage                 3.68 <0.005      12.05Progesterone_Status     4.57 <0.005      17.68Race                    3.45 <0.005      10.78Reginol_Node_Positive   4.31 <0.005      15.89Regional_Node_Examined -4.68 <0.005      18.42T_Stage                 5.80 <0.005      27.13---Concordance = 0.75Partial AIC = 7054.62log-likelihood ratio test = 394.83 on 9 df-log2(p) of ll-ratio test = 261.63

登录后复制抱歉，数学公式解析失败，可能是拼写错误或暂不支持。In [38]

#绘制生存曲线S(t)cph.baseline_survival_.plot()

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [39]

#累积风险曲线H(t)cph.baseline_cumulative_hazard_.plot()

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [40]

#因为累积风险函数 H(t) = -lnS(t)#现在计算H(t) 再绘制曲线，结果和上面的图一致。(-np.log(cph.baseline_survival_)).plot()

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [41]

#绘制每个时间点对应的风险函数h(t)cph.baseline_hazard_.plot()

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [45]

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.573792dtype: float6496.0inf              11.0    0.0000002.0    0.0001183.0    0.0003534.0    0.0010595.0    0.001414...         ...103.0  0.103245104.0  0.103245105.0  0.103245106.0  0.103245107.0  0.103245[107 rows x 1 columns]              11.0    1.0000002.0    0.9998823.0    0.9996474.0    0.9989415.0    0.998587...         ...103.0  0.901906104.0  0.901906105.0  0.901906106.0  0.901906107.0  0.901906[107 rows x 1 columns]

登录后复制In [46]

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.71084628391980030.7460166658656356

登录后复制

              11.0    1.0000002.0    0.9998823.0    0.9996474.0    0.9989415.0    0.998587...         ...103.0  0.901906104.0  0.901906105.0  0.901906106.0  0.901906107.0  0.901906[107 rows x 1 columns]

登录后复制In [ ]

"""可以通过k折交叉验证，把所有协变量的组合都尝试一遍，找到最好的组合这里的组合有8190多个，运行时间会很长很长，可以不必运行，不影响整个项目"""from lifelines.utils import k_fold_cross_validationfrom itertools import combinationsdef combine(temp_list, n):    '''根据n获得列表中的所有可能组合（n个元素为一组）'''    temp_list2 = []    for c in combinations(temp_list, n):        temp_list2.append(c)    return temp_list2varList = 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 协变量值如何影响生存曲线

In [47]

"""将模型的基线生存曲线与一组中协变量值发生变化时发生的情况进行比较。当我们改变协变量(s)时，在其他条件相同的情况下，这有助于比较受试者的生存期。"""cph.plot_partial_effects_on_outcome('T_Stage',values=[1,2,3,4])

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [48]

"""观测多个协变量的值是如何影响生存曲线"""cph.plot_partial_effects_on_outcome(['T_Stage', 'N_Stage'], values=[[0, 0], [1, 0], [1,1], [1,2]], cmap='coolwarm')

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [49]

#绘制平滑校准曲线from lifelines.calibration import survival_probability_calibrationfig, 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.0033427008106428234E50 =  0.0022028810234795415ICI =  0.008988198843408044E50 =  0.008933834566536292

登录后复制

'\nICI – 预测和观察到的平均绝对差\nE50 – 预测和观察到的绝对差中位数\n'

登录后复制

登录后复制

3.11 检验比例风险(PH)假定

之前说过Cox回归模型需要满足比例风险假定，即协变量不是随时间变化而变化，例如某个协变量是是否接受某种药剂的治疗，但是这种治疗不是只发生一次，而是随着时间的增加，治疗的药量会越来越多。不像婚姻状态、人种、性别等协变量是不随时间变化的。

【如何检验PH假定】

满足PH假定和不满足PH假定，建立的cox回归模型需要特别的处理，这里可以通过check_assumptions()工具来检验所有协变量是否满足PH假定。这个工具通过统计检验法，将给出统计学检验的P值，通过P值判断。一般P小于0.05的协变量是不满足比例分险假定。 check_assumptions（）工具输出有三部分

第一部分：

零假设变量是满足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', ...]` inthe 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回归进行拟合模型。

In [50]

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, somecovariates 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 hazardassumption will be flagged.With that in mind, it's best to use a combination of statistical tests and visual tests to determinethe 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 squareddegrees_of_freedom = 1             model =          test_name = proportional_hazard_test---                             test_statistic      p  -log2(p)Age                    km              3.22   0.07      3.78                       rank            4.69   0.03      5.04Estrogen_Status        km              4.99   0.03      5.29                       rank            6.11   0.01      6.22Grade                  km              0.03   0.86      0.22                       rank            0.01   0.91      0.13N_Stage                km              1.31   0.25      1.99                       rank            1.26   0.26      1.93Progesterone_Status    km             10.13 <0.005      9.42                       rank            9.77 <0.005      9.14Race                   km              0.22   0.64      0.65                       rank            0.07   0.79      0.34Reginol_Node_Positive  km              0.85   0.36      1.49                       rank            0.84   0.36      1.48Regional_Node_Examined km              0.00   0.95      0.07                       rank            0.04   0.84      0.25T_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 benon-linear terms missing. The proportional hazard test used is very sensitive to incorrectfunctional 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', ...]` inthe 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

登录后复制

[[,  ], [,  ], [,  ]]

登录后复制

登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制

登录后复制登录后复制登录后复制登录后复制

通过check_assumptions()工具PH检验，发现一共有三个协变量是不满足PH假定，分别是Age(连续变量) 、Estrogen_Status(分类变量)、Progesterone_Status(分类变量)。对于这些变量可以通过"分层"来处理。例如对总样本根据Estrogen_Status分类变量对数据进行划分子样本，再对所有子样本进行cox建模。连续变量Age也可以通过分层来处理，先把连续变量通过区间划分成多个等级(连续变量->分类变量)，再进行对所有子样本进行cox建模。

In [51]

"""先处理分类变量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, somecovariates 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 hazardassumption will be flagged.With that in mind, it's best to use a combination of statistical tests and visual tests to determinethe 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 squareddegrees_of_freedom = 1             model =          test_name = proportional_hazard_test---                             test_statistic    p  -log2(p)Age                    km              3.27 0.07      3.82                       rank            6.04 0.01      6.16Grade                  km              0.02 0.89      0.18                       rank            0.89 0.35      1.54N_Stage                km              1.41 0.24      2.09                       rank            0.05 0.82      0.28Race                   km              0.39 0.53      0.91                       rank            0.72 0.40      1.34Reginol_Node_Positive  km              1.08 0.30      1.74                       rank            1.06 0.30      1.72Regional_Node_Examined km              0.00 0.97      0.04                       rank            0.34 0.56      0.84T_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 benon-linear terms missing. The proportional hazard test used is very sensitive to incorrectfunctional 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

登录后复制

[[,  ]]

登录后复制

登录后复制登录后复制登录后复制登录后复制In [58]

"""现在处理连续变量 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)

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制In [59]

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.

登录后复制

[]

登录后复制In [61]

#最后打印一下处理完不符合PH假定的协变量的cox模型#发现c-index是变低了。cph_strata.print_summary()cph_strata.plot()

登录后复制

             duration col = 'Survival_Months'                event col = 'Status'                   strata = ['age_strata', 'Estrogen_Status', 'Progesterone_Status']      baseline estimation = breslow   number of observations = 3219number of events observed = 479   partial log-likelihood = -2590.82         time fit was run = 2024-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%covariate                                                                                                                      Grade                    0.45       1.57       0.08             0.29             0.60                 1.34                 1.83N_Stage                  0.36       1.44       0.10             0.18             0.55                 1.19                 1.73Race                     0.42       1.53       0.12             0.20             0.65                 1.22                 1.92Reginol_Node_Positive    0.05       1.05       0.01             0.03             0.08                 1.03                 1.08Regional_Node_Examined  -0.04       0.97       0.01            -0.05            -0.02                 0.95                 0.98T_Stage                  0.32       1.38       0.06             0.21             0.43                 1.23                 1.54                           z      p   -log2(p)covariate                                     Grade                   5.73 <0.005      26.54N_Stage                 3.83 <0.005      12.95Race                    3.64 <0.005      11.84Reginol_Node_Positive   4.14 <0.005      14.79Regional_Node_Examined -4.70 <0.005      18.59T_Stage                 5.59 <0.005      25.38---Concordance = 0.70Partial AIC = 5193.64log-likelihood ratio test = 257.89 on 6 df-log2(p) of ll-ratio test = 172.99

登录后复制抱歉，数学公式解析失败，可能是拼写错误或暂不支持。

登录后复制

登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制登录后复制

3.12 真的需要关心比例风险假设吗

lifelines 作者认为：

免责声明

游乐网为非赢利性网站，所展示的游戏/软件/文章内容均来自于互联网或第三方用户上传分享，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系youleyoucom@outlook.com。

同类文章