在前面文章中介绍了Cox比例风险回归模型(Cox proportional hazard regression model)，本篇文章将实例演示在Stata软件中实现时间依存Cox回归模型的操作步骤。
关键词：Stata; Cox等比例回归; Cox回归; 生存分析; 等比例风险; 含时间依存协变量; 时依协变量
一、案例介绍
某肿瘤研究所收集了200例肺癌患者的生存数据：包括生存状态(status，0=“删失”，1=“死亡”)、生存时间(time，天)、性别(sex，1=“男”，2=“女”)、年龄 (age，岁)和卡氏评分(phkarno)，部分数据见图1。本文案例可从“附件下载”处下载。现欲探究患者的性别、年龄、卡氏评分与生存结局的关系。
二、问题分析
本案例的目的是探究肺癌患者的性别、年龄和卡氏评分与生存结局的关系，可以采用Cox比例风险回归模型进行分析，但需要满足5个条件：
条件1：因变量是含有时间信息的二分类变量。本案例中因变量是包含生存时间的二分类资料，time是生存时间(天)；status是生存结局。本案例数据满足该条件。
条件2：各观测值之间相互独立，无互相干扰。由数据和研究设计可知，该条件满足。
条件3：一般要求结局事件的样本量为自变量个数的10~20倍(EPV原则)。该条件需要软件分析来判断。
条件4：自变量之间无严重多重共线性。该条件需要软件分析来判断。
条件5：等比例风险(Proportional hazards，PH)假设，该条件需要软件分析来判断。
三、软件操作及结果解读
(一) 适用条件判断
1. 条件3判断
(1) 软件操作
*查看status的分类频数*
tab status
(2) 结果解读
由图2结果可知因变量中结局事件为158例。按照EPV为10~20的原则可满足多因素模型纳入7~15个变量的需求。因此，本案例数据满足条件2。
2. 条件4判断 (共线性检测)
(1) 软件操作
变量之间的共线性诊断可以采用容忍度法 (Tolerance)、方差膨胀因子法 (VIF)、相关系数法和条件数法。本案例采用容忍度法、方差膨胀因子法来查看变量之间的多重共线性，可在线性回归后计算容忍度和方差膨胀因子。结果如图3所示。
*拟合线性回归模型*
regress status age sex phkarno
*计算VIF和容忍度*
vif
(2) 结果解读
图3结果中列出了自变量的方差膨胀因子(VIF)和容忍度(1/VIF)。可见，三个自变量的VIF均远小于10，容忍度均远大于0.1，提示变量间不存在严重的多重共线性，可知满足条件4。
3. 条件5判断 (等比例风险判定)
PH假定指自变量对生存率[风险比值h(t)/h0(t)]的影响不会随时间的变化而变化。若不满足PH假定，则需要改用含时间依存协变量的Cox比例模型分析数据。本案例PH假定的检验详见本章后文。
(二) 统计描述及推断
1. 统计描述
(1) 软件操作
描述性统计分析，结果如图4-1—图4-3所示。
*描述sex变量*
tab sex
*描述time、age、phkarno变量*
tabstat time age phkarno, stat(count mean median sd min max) col(stat) long
(2) 结果解读
如图4-1—图4-2和图2所示，本研究共纳入200例肺癌患者，其中男性125例(62.50%)，女性75例(37.50%)；所有研究对象的中位生存时间为285天，随访过程中158例(79.00%)死亡；所有研究对象的年龄为(62.34±9.00)岁，卡氏评分为(81.90±12.50)分。
2. 单因素分析
(1) 软件操作
①将数据转化为生存数据格式，结果如图5所示。
stset time, id(id) failure(status==1)
②对年龄、性别和卡氏评分进行单因素Cox比例风险模型分析，结果如图6-1—图6-3所示。
*age单因素分析*
stcox age
*sex单因素分析*
stcox i.sex
* phkarno单因素分析*
stcox phkarno
(2) 结果解读
单因素Cox比例风险模型分析显示(图6-1—图6-3)，年龄、性别和卡氏评分与生存结局之间均存在统计学关联(P<0.1)，年龄每增加一岁肺癌患者的死亡风险增加1.7% (HR=1.017，95%CI：0.999~1.036；P=0.069)；女性死亡风险比男性低38.8% (HR=0.612，95%CI：0.438~0.856；P=0.004)；卡氏评分每上升一个单位肺癌患者死亡的风险下降1.5% (HR=0.985，95%CI：0.974~0.997；P=0.013)。
3. 多因素分析
将单因素Cox比例风险模型分析中有统计学关联的变量(年龄、性别、卡氏评分)纳入到多因素Cox比例风险模型。
(1) 软件操作
①模型拟合，结果如图7所示。
stcox age i.sex phkarno
②对模型中的变量进行PH假定检验，结果如图8所示。
estat phtest, detail
(2) 结果解读
多因素Cox比例风险模型分析显示(图7)，性别和卡氏评分与生存结局之间存在统计学关联(P<0.05)，女性的死亡风险比男性低37.5% (HR=0.625,95%CI：0.447~0.874；P=0.006)；卡氏评分每增加一个单位，肺癌患者死亡风险下降1.3%(HR=0.988；95%CI：0.977~1.000；P=0.047)。年龄每增加一岁肺癌患者的死亡风险增加1.2% (HR=1.012，95%CI：0.993~1.031；P=0.221)，但关联无统计学意义。
对多因素Cox比例风险模型中的自变量进行PH假定检验(图8)，结果显示年龄(age)和性别(sex)的P值>0.1，而卡氏评分(phkarno)的P值<0.05，提示年龄和性别满足PH假定，卡氏评分不满足PH假定；整体检验P值<0.05，提示不满足条件5。
4. 含时间依存协变量的Cox比例风险模型
由于卡氏评分(phkarno)不满足PH假定，因此需要设立卡氏评分与时间的交互项，与年龄、性别一起建立含时间依存协变量的Cox比例风险模型。
(1) 软件操作
拟合具有回归系数的模型，结果如图9所示。
stcox age i.sex phkarno, tvc(phkarno) texp(ln(_t+20)) nohr
拟合具有效应量HR值的模型，结果如图10所示。
stcox age i.sex phkarno, tvc(phkarno) texp(ln(_t+20))
(2) 结果解读
图9和图10中的结果分为两个部分，上方的结果(main)是对原始变量的相应估计。下方的结果(tvc)是对构造的时间依存协变量的估计[即phkarno*ln(time)]，P=0.038，<0.05，提示自变量tvc具有时间依存性，进一步证实了其不满足风险比例Cox回归模型的PH假定要求，故此处应采用时间依存协变量Cox回归模型。
分析结果显示，性别的效应值HR为0.6162506，表示女性死亡的风险比男性低38.4% (HR=0.616，95%CI：0.441~0.861；P=0.005)，年龄与生存结局的关联无统计学意义 (HR=1.012，95%CI：0.993~1.032；P=0.202)。卡氏评分的时依系数β(t)= -0.0936858+0.0146965×ln(t+20)，效应值HR=exp(-0.0936858+0.0146965×ln(t+20))。例如，当时间为200天时，卡氏评分对应的HR=exp(-0.0936858+0.0146965×ln(200+20))= 0.9856851。
四、结论
本研究采用Cox比例风险模型探究年龄、性别、卡氏评分与肺癌患者生存结局的关系，数据样本量满足要求，自变量之间无严重多重共线性；PH假定检验发现卡氏评分不满足PH假定(P=0.0059)，需要建立含时间依存协变量的Cox比例风险模型。
含时间依存协变量的Cox比例风险模型分析结果表明，性别和卡氏评分均是肺癌患者生存结局的影响因素，其中女性死亡的风险比男性低38.4% (HR=0.616，95%CI：0.441~0.861；P=0.005)。年龄与生存结局的关联无统计学意义 (HR=1.012，95%CI：0.993~1.032；P=0.202)。卡氏评分的时依系数β(t)= -0.0936858+0.0146965×ln(t+20)，效应值HR=exp(-0.0936858+0.0146965×ln(t+20))。
五、分析小技巧
- 对于生存数据，如果自变量是分类变量，可以使用Kaplan-Meier生存曲线和log-rank检验进行单因素分析，也可以使用单因素Cox比例风险模型；如果自变量是定量变量，一般单因素分析采用Cox比例风险模型。
- 使用Stata进行生存分析，必须将数据转化为生存数据格式后(stset之后)，才能使用任何其他的st开始的命令。命令stcox后只需要输入自变量，不需输入因变量，原因是在将数据转化为生存数据时(stset)，已指定过终点事件(Failure variable)和时间变量(Time variable)。
- 在使用estat phtest命令检验PH假定前，需要首先运行要检验的Cox模型，即estat phtest命令需要紧跟在Cox回归的命令之后，否则Stata不知道检验哪个回归的PH假定。
- 构建时依协变量时，常选择协变量与时间函数的乘积项。案例中对ln(t)设置为ln(t+20)，增加了一个常数，主要作用是使卡氏评分随着时间的死亡风险变化在一段时间范围内呈直线，常数项的大小一般对模型估计结果影响不大。