在前面文章中介绍了Cox比例风险回归模型(Cox Proportional Hazard Regression Model)的假设检验理论，本篇文章将实例演示在R软件中实现时间依存Cox回归模型的操作步骤。
关键词：R语言; R软件; Cox等比例回归; Cox回归; 生存分析; 等比例风险; 含时间依存协变量; 时依协变量
一、案例介绍
某肿瘤研究所收集了200例肺癌患者的生存数据：包括生存状态(status，1=“删失”，2=“死亡”)、生存时间(time，天)、性别(sex，1=“男”，2=“女”)、年龄 (age，岁)和卡氏评分(ph.karno)，部分数据见图1。现欲探究患者的性别、年龄、卡氏评分与生存结局的关系。部分数据见图1。本文案例可从“附件下载”处下载。
二、问题分析
本案例的目的是探究肺癌患者的性别、年龄和卡氏评分与生存结局的关系，可以采用Cox比例风险回归模型进行分析，但需要满足5个条件：
条件1：因变量是含有时间信息的二分类变量。本案例中因变量是包含生存时间的二分类资料，time是生存时间(天)；status是生存结局。本案例数据满足该条件。
条件2：各观测值之间相互独立，无互相干扰。由数据和研究设计可知，该条件满足。
条件3：一般要求结局事件的样本量为自变量个数的10~20倍(EPV原则)。该条件需要软件分析来判断。
条件4：自变量之间无严重多重共线性。该条件需要软件分析来判断。
条件5：等比例风险(Proportional hazards，PH)假设，该条件需要软件分析来判断。
三、软件操作及结果解读
(一) 导入数据
mydata <- read.csv("Cox比例风险模型.csv") #导入CSV数据
View(mydata) #查看数据
在数据栏目中可以查看全部数据情况，数据集中共有6个变量和200个观察数据，6个变量分别代表患者编号(ID)、生存时间(time，天)、生存状态(status，1=“删失”，2=“死亡”)、年龄 (age，岁)、性别(sex，1=“男”，2=“女”)和卡氏评分(ph.karno)。
如果数据集较大也可使用如下命令查看数据框结构：
str(mydata) #查看数据框结构
(二) 适用条件判断
1. 条件3判断
(1) 软件操作
##查看status的分类频数##
table(mydata$status) #计数 prop.table(table(mydata$status)) #计算百分比
(2) 结果解读
图3计算了删失与死亡的例数及其所占的百分比，结果可知结局事件为158例。按照EPV为10~20的原则可满足多因素模型纳入7~15个变量的需求。因此，本案例数据满足条件3。
2. 条件4判断 (共线性检测)
(1) 软件操作
##共线性诊断##
library(car) lm.reg<-lm(status~age+sex+ph.karno,data=mydata) #线性回归分析 vif(lm.reg) #计算vif
(2) 结果解读
图4结果中列出了自变量的方差膨胀因子，可见，三个自变量的VIF均远小于10，提示变量间不存在严重的多重共线性，可知满足条件4。
3. 条件5判断 (等比例风险判定)
PH假定指自变量对生存率[风险比值h(t)/h0(t)]的影响不会随时间的变化而变化。若不满足PH假定，则需要改用含时间依存协变量的Cox比例模型分析数据。本案例PH假定的检验详见本章后文。
(三) 统计学描述及推断
1.统计描述
(1) 软件操作
##统计描述##
table(mydata$sex) #性别构成 prop.table(table(mydata$sex)) # 男女占比
library(psych) #调用“psych”包 describe(mydata$time) describe(mydata$age) describe(mydata$ph.karno)
(2) 结果解读
如图5所示，本研究共纳入200例肺癌患者，其中男性125例(62.50%)，女性75例(37.50%)；如图6所示，所有研究对象的中位生存时间为285天，随访过程中158例(79.00%)死亡；所有研究对象的年龄为(62.34±9.00)岁，卡氏评分为(81.90±12.50)分。
2. 单因素分析
(1) 软件操作
##单因素cox回归##
library(survival) #调用包“survival” fit.sex<-coxph(Surv(time,status)~factor(sex),data=mydata) #性别的单因素分析 summary(fit.sex) #查看结果
fit.age<-coxph(Surv(time,status)~age,data=mydata) #年龄的单因素分析 summary(fit.age) #查看结果
fit.karno<-coxph(Surv(time,status)~ph.karno,data=mydata) # karno评分的单因素分析 summary(fit.karno) #查看结果
(2) 结果解读
单因素Cox比例风险模型分析显示(图7-1—图7-3)，性别、年龄和卡氏评分与生存结局之间均存在统计学关联(P<0.1)，女性死亡风险比男性低38.8% (HR=0.612，95%CI：0.438~0.856；P=0.004)；年龄每增加一岁肺癌患者的死亡风险增加1.7% (HR=1.017，95%CI：0.999~1.036；P=0.069)；卡氏评分每上升一个单位肺癌患者死亡的风险下降1.5% (HR=0.985，95%CI：0.974~0.997；P=0.013)。
3. 多因素分析
将单因素Cox比例风险模型分析中有统计学关联的变量(年龄、性别、卡氏评分)纳入到多因素Cox比例风险模型。
(1) 软件操作
##多因素cox回归##
fit<-coxph(Surv(time,status)~factor(sex)+age+ph.karno,data=mydata) summary (fit)
##PH假定检验##
ph.test<-cox.zph(fit) ph.test
(2) 结果解读
多因素Cox比例风险模型分析显示(图8)，性别和卡氏评分与生存结局之间存在统计学关联(P<0.05)，女性的死亡风险比男性低37.5% (HR=0.625，95%CI：0.447~0.874；P=0.006)；卡氏评分每增加一个单位，肺癌患者死亡风险下降1.2%(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.220)，但关联无统计学意义。
对多因素Cox比例风险模型中的自变量进行PH假定检验(图9)，结果显示年龄(age)和性别(sex)的P值>0.1，而卡氏评分(phkarno)的P值<0.05，提示年龄和性别满足PH假定，卡氏评分不满足PH假定；整体检验P值<0.05，提示不满足条件5。
4. 含时间依存协变量的Cox比例风险模型
由于卡氏评分(ph.karno)不满足PH假定，因此需要设立卡氏评分与时间的交互项，与年龄、性别一起建立含时间依存协变量的Cox比例风险模型。
(1) 软件操作
fit.ph<-coxph(Surv(time,status)~sex+age+ph.karno+tt(ph.karno),
tt=function(x,t,...) x*log(t+20),data=mydata)
summary(fit.ph)
(2) 结果解读
由图10可知所构造的时间依存协变量tt(ph.karno)的P=0.038<0.05，提示自变量ph.karno具有时间依存性，进一步证实了其不满足等比例风险Cox回归模型的PH假定要求，故此处应采用时间依存协变量Cox回归模型。
分析结果显示，性别的效应值HR为0.6159，表示女性死亡的风险比男性低38.4% (HR=0.616，95%CI：0.441~0.861；P=0.005)，年龄与生存结局的关联无统计学意义 (HR=1.012，95%CI：0.994~1.032；P=0.197)。卡氏评分的时依系数β(t)= -0.094+0.015×ln(t+20)，效应值HR=exp(-0.094+0.015×ln(t+20))。例如，当时间为200天时，卡氏评分对应的HR=exp(-0.094+0.015×ln(200+20))= 0.987。
四、结论
本研究采用Cox比例风险模型探究年龄、性别、卡氏评分与肺癌患者生存结局的关系，数据样本量满足要求，自变量之间无严重多重共线性；PH假定检验发现卡氏评分不满足PH假定(P=0.008)，需要建立含时间依存协变量的Cox比例风险模型。
含时间依存协变量的Cox比例风险模型分析结果表明，性别和卡氏评分均是肺癌患者生存结局的影响因素，其中女性死亡的风险比男性低38.4% (HR=0.616，95%CI：0.441~0.861；P=0.005)。年龄与生存结局的关联无统计学意义 (HR=1.013，95%CI：0.994~1.032；P=0.197)。卡氏评分的时依系数β(t)= -0.094+0.015×ln(t+20)，效应值HR=exp(-0.094+0.015×ln(t+20))。
五、分析小技巧
- 对于生存数据，如果自变量是分类变量，可以使用Kaplan-Meier生存曲线和log-rank检验进行单因素分析，也可以使用单因素Cox比例风险模型；如果自变量是定量变量，一般单因素分析采用Cox比例风险模型。
六、知识小贴士
PH假定检验可以通过以下四种方法判断：
- 分类协变量的Kaplan-Meier生存曲线间无交叉。
- 利用Schoenfeld残差法查看连续性变量对应的P值和残差图。
- 绘制log cumulative hazard curve(对数累积风险曲线)，曲线平行是满足PH假定的充分且必要条件。
- 可将每个协变量与对数生存时间的交互作用项放入模型中，如果交互项统计学上不显著，则满足等比例风险条件。