时间依存Cox回归模型(Time-Dependent Cox Regression Model)——Python软件实现

关键词：Python; Cox等比例回归; Cox回归; 生存分析; 等比例风险; 含时间依存协变量; 时依协变量

一、案例介绍

某肿瘤研究所收集了200例肺癌患者的生存数据：包括生存状态(status，1=删失，2=死亡)、生存时间(time，天)、性别(sex，1=男，2=女)、年龄 (age，岁)和卡氏评分(ph.karno)。现欲探究患者的性别、年龄、卡氏评分与生存结局的关系。部分数据见图1，本案例数据可从“附件下载”处下载。

二、问题分析

本案例的目的是探究肺癌患者的性别、年龄和卡氏评分与生存结局的关系，可以采用Cox比例风险回归模型进行分析，但需要满足5个条件：

条件1：因变量是含有时间信息的二分类变量。本案例中因变量是包含生存时间的二分类资料，time是生存时间(天)；status是生存结局，满足该条件。

条件2：各观测值之间相互独立，无互相干扰。由数据和研究设计可知，满足该条件。

条件3：一般要求结局事件的样本量为自变量个数的10~20倍(EPV原则)。该条件需要软件分析来判断。

条件4：自变量之间无严重多重共线性。该条件需要软件分析来判断。

条件5：等比例风险(Proportional hazards，PH)假设。该条件需要软件分析来判断。

三、软件操作及结果解读

(一) 导入数据

import pandas as pd
df = pd.read_csv(r'Cox比例风险模型.csv')     #导入CSV数据
df

在数据栏目(图1)中可以查看全部数据情况，数据集中共有6个变量和200例观察数据，6个变量分别代表患者编号(ID)、生存时间(time，天)、生存状态(status，1=删失，2=死亡)、年龄(age，岁)、性别(sex，1=男，2=女)和卡氏评分(ph.karno)。

(二) 适用条件判断

1. 条件3判断

(1) 软件操作

##查看status的分类频数##

df.groupby("status").count()

(2) 结果解读

删失与死亡的例数结果见图2，结果可知结局事件为158例。按照EPV为10~20的原则可满足多因素模型纳入7~15个变量的需求。因此，本案例数据满足条件3。

2. 条件4判断 (共线性检测)

(1) 软件操作

##共线性诊断##

import numpy as np
import pandas as pd
from statsmodels.stats.outliers_influence import variance_inflation_factor
##添加截距项，因为variance_inflation_factor函数拟合的OLS不包含截距项##
data_ls = df.copy()
data_ls.loc[:,'intercept'] = 1
##计算方差膨胀因子#
vif = [variance_inflation_factor(np.matrix(data_ls.drop(['ID','time','status'],axis = 1)), i) for i in range(data_ls.drop(['ID','time','status'],axis = 1).shape[1])]
pd.DataFrame(data = zip(data_ls.drop(['ID','time','status'],axis = 1),vif),columns = ['变量','vif'])

(2) 结果解读

自变量的方差膨胀因子结果见图3，可知3个自变量的VIF均远小于10，提示变量间不存在严重的多重共线性，满足条件4。

3. 条件5判断 (等比例风险判定)

PH假定指自变量对生存率[风险比值h(t)/h0(t)]的影响不会随时间的变化而变化。若不满足PH假定，则需要改用含时间依存协变量的Cox比例模型分析数据。本案例PH假定的检验详见本章后文。

 (三) 统计学描述及推断

1. 统计描述

(1) 软件操作

##统计描述##

df.groupby("sex").count()    #性别构成

df['time'].describe()     #对时间进行统计描述

df['age'].describe()     #对年龄进行统计描述

df['ph.karno'].describe()     #对评分进行统计描述

(2) 结果解读

    各变量描述统计结果见图4~7，本研究共纳入200例肺癌患者，平均年龄为(62.34±9.00)岁；其中男125例，女75例；随访过程中158例死亡，中位生存时间为285 d，；平均卡氏评分为(81.90±12.50)分。

2. 单因素分析

(1) 软件操作

##单因素Cox回归##

##对性别进行单因素分析##

df1=df[['time', 'status', 'sex']]
df1

from lifelines import CoxPHFitter
cph = CoxPHFitter()
cph.fit(df1, 'time', event_col='status')
cph.print_summary()

##对年龄进行单因素分析##

df2=df[['time', 'status', 'age']]
df2

cph = CoxPHFitter()
cph.fit(df2, 'time', event_col='status')
cph.print_summary()

##对评分进行单因素分析##

df3=df[['time', 'status', 'ph.karno']]
df3

cph = CoxPHFitter()
cph.fit(df3, 'time', event_col='status')
cph.print_summary()

(2) 结果解读

单因素Cox比例风险模型分析结果(图8~13)显示，性别、年龄和卡氏评分与生存结局之间均存在统计学关联(P<0.1)，女性死亡风险比男性低39.0% (HR=0.61，95%CI为0.44~0.86，P<0.005)；年龄每增加1岁肺癌患者的死亡风险增加2.0% (HR=1.02，95%CI为1.00~1.04；P=0.07)；卡氏评分每上升1个单位肺癌患者死亡的风险下降1.0% (HR=0.99，95%CI为0.97~1.00；P=0.01)。

3. 多因素分析

将单因素Cox比例风险模型分析中有统计学关联的变量(年龄、性别、卡氏评分)纳入到多因素Cox比例风险模型。

(1) 软件操作

##多因素Cox回归##

cph = CoxPHFitter()
cph.fit(df.drop('ID',axis = 1), 'time', event_col='status')
cph.print_summary(decimals=3)

##PH假定检验##

cph.check_assumptions(df.drop('ID',axis = 1), p_value_threshold=0.05) 

(2) 结果解读

多因素Cox比例风险模型分析结果(图14)显示，性别和卡氏评分与生存结局之间存在统计学关联(P=0.006)，女性的死亡风险比男性低37.5% (HR=0.625，95%CI为0.447~0.874；P=0.006)；卡氏评分每增加1个单位，肺癌患者死亡风险下降1.2%(HR=0.988；95%CI为0.977~1.000；P=0.047)。年龄每增加1岁肺癌患者的死亡风险增加1.2% (HR=1.012，95%CI为0.993~1.031；P=0.220)，但关联无统计学意义。

多因素Cox比例风险模型自变量PH假定检验结果(图15)显示，年龄(age)和性别(sex)的P＞0.1，而卡氏评分(phkarno)的P＜0.05，提示年龄和性别满足PH假定，卡氏评分不满足PH假定，不满足条件5，所以采用含时间依存协变量的Cox比例风险模型。

4. 含时间依存协变量的Cox比例风险模型

设立卡氏评分与时间的交互项(图16)，与年龄、性别一起建立含时间依存协变量的Cox比例风险模型。

(1) 软件操作

# 数据整理

from lifelines.utils import to_episodic_format
df2 = to_episodic_format(df, duration_col='time', event_col='status', time_gaps=1.)
df2['time*ph.karno'] = df2['ph.karno'] * (np.log(df2['stop']+20))
df2

# 模型拟合

from lifelines import CoxTimeVaryingFitter
cph = CoxTimeVaryingFitter()
cph.fit(df2.drop('ID',axis = 1), event_col='status',id_col='id')
cph.print_summary()

(2) 结果解读

含时间依存协变量的Cox比例风险模型多因素检验结果(图17)显示，时间依存协变量time*ph.karno的P=0.04＜0.05，提示自变量ph.karno具有时间依存性，进一步证实了其不满足等比例风险Cox回归模型的PH假定要求，故此处应采用时间依存协变量Cox回归模型。

分析结果(图17)显示，女性死亡的风险比男性低38.0% (HR=0.62，95%CI为0.44~0.86；P＜0.005)；年龄与生存结局的关联无统计学意义 (HR=1.01，95%CI：0.99~1.03；P=0.2)。卡氏评分的时依系数β(t)= -0.09+0.01×ln(t+20)，效应值HR= exp(-0.09+0.01×ln(t+20))。例如，当时间为200d时，卡氏评分对应的HR= exp(-0.09+0.01×ln(200+20))= 0.965。

四、结论

本研究采用Cox比例风险模型探究年龄、性别、卡氏评分与肺癌患者生存结局的关系，数据样本量满足要求，自变量之间无严重多重共线性；PH假定检验发现卡氏评分不满足PH假定(P＜0.05)，需要建立含时间依存协变量的Cox比例风险模型。

含时间依存协变量的Cox比例风险模型分析结果表明，性别和卡氏评分均是肺癌患者生存结局的影响因素，其中女性死亡的风险比男性低38.0% (HR=0.62，95%CI：0.44~0.86；P<0.005)。年龄与生存结局的关联无统计学意义(HR=1.01，95%CI：0.99~1.03；P=0.2)。卡氏评分的时依系数β(t)= -0.09+0.01×ln(t+20)，效应值HR= exp(-0.09+0.01×ln(t+20))。

五、分析小技巧

对于生存数据，如果自变量是分类变量，可以使用Kaplan-Meier生存曲线和log-rank检验进行单因素分析，也可以使用单因素Cox比例风险模型；如果自变量是定量变量，一般单因素分析采用Cox比例风险模型。

六、知识小贴士

PH假定检验可以通过以下4种方法判断：

• 分类协变量的Kaplan-Meier生存曲线间无交叉
• 利用Schoenfeld残差法查看连续性变量对应的P值和残差图
• 绘制log cumulative hazard curve(对数累积风险曲线)，曲线平行是满足PH假定的充分且必要条件
• 将每个协变量与对数生存时间的交互作用项放入模型中，如果交互项统计学上不显著，则满足等比例风险条件。