【原】STATA：生存分析中的时变协变量与时依系数

Memo_Cleon 2022-12-05 发布于上海

展开全文

有读者让写个生存分析时依系数法的stata教程。

时依协变量（Time-Dependent Covariate），或称时变协变量（Time-Varying Covariate），顾名思义指的就是随时间变化而变化的解释变量。这里面存在两种情况。一种是随时间变化自变量本身发生了变化，比如有些受试者在随访过程戒烟了，由吸烟变为不抽烟，这种时依协变量也被称为内在时依协变量。还有一种情况，随着时间的变化，模型中自变量本身取值并未发生改变，但其效应却在发生变化，这种时依协变量被称为外在时依协变量。

目前时依协变量的教程多是针对外在时依协变量的，包括我们以前的笔记《非比例风险的Cox回归模型_时依系数法》和《非比例风险的Cox回归模型_分段模型》，而对自变量本身发生变化的介绍较少。

在stata里面，时变协变量指的是“狭义的”时变协变量，专指内在时依协变量，变化在变量本身；而在处理外在时依协变量时，变化在系数（时依系数）。如果你想考察时变协变量的临床效应，①对于外在时依协变量和连续性的内在时变协变量，stata的处理方法类似，不论是stcox（Cox比例风险模型）还是stcrreg（竞争风险模型），通过指定tvc()和texp()两个参数就可以搞定。tvc指定时依协变量，texp指定时依协变量中时间的纳入形式（比如常用经过ln变换），即时变协变量的时间乘数，默认是_t；②而对于分类的内在时变协变量，需要根据时变协变量变化前后将一条记录分裂成多条，以保证每条记录的时变协变量取值保持恒定，然后再利用时依数据集进行分析，有点类似我们在《非比例风险的Cox回归模型_分段模型》介绍过的将数据拆分成短期模型和长期模型，不同的是最终使用的这个时依数据集是短期模型和长期模型分析集的合集。

此次笔记内容：

【1】内在分类时变协变量效应考察

【2】内在连续时变协变量效应考察

【3】外在时变协变量效应考察（时依系数法）

示例1：内在分类时变协变量效应考察

斯坦福心脏移植项目(Crowley and Hu 1977) 的103例患者【id:患者编号；year：参加项目的年份；age：患者年龄；surgery：心脏手术史（0=No，1=Yes）；died：患者结局（0=alive；1=dead）；stime：生存时间（天）；waittime：等待供体时间（天）；transplant：心脏移植（0=未移植；1=移植）】，考察心脏移植手术是否会延长患者的生存时间。

这里的transplant就是一个时变协变量。在患者登记到该项目后，需要等待心脏供体，在等待期间，有的患者发生了死亡，有的患者则退出了该项目。最终有69例（66.99%）患者接受了心脏移植，这部分患者本身发生了变化（由等待移植变成移植）。

你可能觉得这很简单，同时考虑age、surgery，构造以一个时依协变量[transplant*(T_>waittime)]不就可以啦。但你很快就会发现当T_>waittime时，transplant全部取值为1，其效应显然没法计算。如果理清了变量效应与取值之间的关系就可以解决这个问题：transplant在模型里面不是纳入其主效应，而是transplant*(T_>waittime)。即最终模型应该包含age、surgery和transplant*(T_>waittime)：

h(t)=h₀(t)*exp(β_1*age+β₂*surgery+β₃*transplant+β₄*transplant*(T_>waittime))

这里的transplant*(T_>waittime)的效应实际上就是内在时变协变量transplant的效应，比较遗憾的是在stata中并不支持(T_>waittime)这种条件格式，在这点上笔者觉得SPSS做得更好。

我们先把SPSS过程和主要结果放上来：

STATA数据准备：重新构建分析数据集（时依数据集）

按接受心脏移植的时间构建新的变量posttran、time和status，分别标识心脏移植前后（0=移植前，1=移植后）、新的生存时间和新的生存结局（0=删失，1=死亡）。未接受心脏移植的患者posttran=transplant，生存时间和结局状态保持不变（time=stime=waittime，status=died）。接受受心脏移植的患者都被拆分成两条记录，一条记录保持原有信息不变（posttran=1，status=died，time=stime）；另外一条记录只包含移植前信息（posttran=0），生存时间定义为等待供体时间（time=waittime），生存状态为删失（status=0），相当于提供了生存时间为waittime的失访信息。

然后利用这个新的数据集进行Cox回归就可以啦！

COX回归

统计>>生存分析>>模型设定和实用工具>>声明数据集为生存时间数据

结果中首先显示相应的命令：stset time, id(id) failure(status==1) scale(1)

结果显示有一个错误，说是有（同一患者）多记录生存时间是一样的。由于接受心脏移植的患者会被拆成2条记录，一条是完整信息，一条是移植前信息，很显然完整信息的生存时间应该大于移植前信息的生存时间。我们筛选了一下，发现id=38的患者的这两个生存时间相等，该患者等待5天，在第5天接受了心脏移植，死于第5天，这在专业上是没问题的，但是在拆分数据时软件就会报错。这里我们只需要把移植前信息保持不变，而把完整信息的生存时间略微长于5就可以了，比如本例改为5.1。然后重新声明和进行Cox回归。

统计>>生存分析>>回归模型>>Cox比例风险（PH）模型

自变量选入age、surgery、posttran。由于surgery、posttran都是二分类变量，直接按连续变量处理不影响结果。另外在生存时间声明时，注意生存时间和失效时间应选入新生成的time和status。此外本例中每个接受心脏移植的患者都被拆分成了2条记录，在声明生存时间时，需要选择多记录的ID变量，即标识患者的编号。

数据集共有103个病例172条记录，存在75个失效事件，31938.1个分析人次，随访时间最长个体是1799天。

结果显示，年龄愈大死亡风险越大（HR=1.03，P=0.019），术前做过其他手术死亡风险会降低（HR=0.33,P=0.010）；在校正age和surgery后，心脏移植手术并没有延长患者的生存时间（HR=1.00，P=0.996）。

结果中默认的估计值是HR及相应结果，如想显示系数，可在[Cox比例风险（PH）模型]对话框的[报告]选显卡中选中“报告系数，而不是风险比”即可。

我们可以看到这个结果与前面用SPSS直接构建时依协变量transplant*(T_>waittime)的结果完全一致。相比之下，SPSS的时依协变量解决这种分类的内在协变量更为简便。但需要注意的是，不可以利用SPSS来计算stata最终使用的这个数据集来完成内在时变协变量的分析，因为SPSS的Cox回归中并没有标识个案的字段，当需要拆分的个案比较多（最终接受移植人数较多）时，结果会有较大的误差。

当然PH假定仍然是模型的的假设条件，Stata里面提供了足够多的方法，感兴趣的可以参照《生存分析之Cox回归》一文。

另外需要特别说明一下，时变协变量分析集的构建是整个分析中比较关键的一步，我是直接在WPS中进行的，只需要利用筛选功能选择transplant=1的记录，将其复制到数据集尾部，然后这部分记录的time=waittime，status=0即可。如果需要分割成多段，stata中的stsplit命令似乎更快捷一些。

生成新变量posttran、time和status：

数据>>创建或更改数据>>创建新变量

相应命令：generate posttran = transplant

同理生成time和status两个变量：

generate status = died

generate time = stime

声明生存时间：菜单操作[统计>>生存分析>>模型设定和实用工具>>声明数据集为生存时间数据]同前，相应命令为：

stset time, id(id) failure(status)

声明生存时间也可以在下面的拆分时间跨度记录的【生存设置】中进行。

拆分记录：统计>>生存分析>>模型设定和实用工具>>拆分时间跨度记录

相应命令为：

stsplit timecat if transplant==1, at(0) after(time=waittime)

因为未移植患者的生存时间等于等待时间，因此这个命令中去掉 if transplant==1结果是一样的。

接受心脏移植的患者被拆分成2条记录，一条是移植前，一条为移植后。移植前的posttran=0，status=0，移植后同原数据（posttran=transplant=1，status=died）。经过上述操作后，移植前的数据赋值存在错误，需要进行修改。可通过replace命令进行，菜单操作如下：

相应的命令如下：

replace status = _d

replace posttran = 0 if transplant==1 & timecat==-1

由于我们对数据集进行了生存数据声明，数据中自动生成几个变量：st（记录是否被使用，1使用）、_d（状态，0删失，1失效）、t0（进入风险集时间）、_t（退出风险集时间，t0=0时_t就是生存时间）。我们可以利用这些信息对数据进行快速地替换。

最后我们在运行stsplit 命令的时候，生成了一个timecat的分类变量，在我们后续分析中用不到，直接删除即可：

drop timecat

示例2：内在连续时变协变量效应考察

生存数据的随访过程就是一个重复测量过程，除了随访结局，研究者可能还会记录每一次随访时的协变量结果。有些协变量不会随时间的变化发生变化，但有些可能会发生变化（时变协变量），比如每次随访时的生化检验指标。我们在处理这些时变协变量时，一般会假设每个随访时间段内，时变协变量值保持不变。注意，这里的时变协变量跟复发事件不是一回事。

比如，梅奥诊所的312例原发性胆汁性肝硬化患者随访数据如下：

以第2例患者（id=2）为例，该患者有9次随访，每次随访的血液检验指标都会发生一些变换，这些指标就都是时变协变量，比如血清蛋白（albumin）在每次随访都会在具体值上有一些波动，在分析时我们会假设每两次临近的随访之间间albumin保持恒定，比如第0-182天，albumin是4.14mg/dl，而182-365天则是3.6mg/dl。这种数据虽然是计量资料，但实际上我们并不知晓两次测量间任意时刻的值，只是用一个值代替罢了，其分析方法同前面分类的内在时变协变量分析方法一样，重点在于时依数据集的构建。上面的数据来自R的survival包，在分析前还需要对变量status进行一些处理，即已经随访到终点事件的患者，在终点事件发生前的随访中status应界定为删失。

药物浓度的变化也是常见的一种时变协变量，在stata的帮助文件中提供了一个这样的例子。45例游走性肺炎【泛指所有由某种未知的病原体引起的肺炎，患者症状相对较轻仍可以行走活动】的康复时间（天），观察协变量有年龄、及drug1和drug2的使用。drug1和drug2为初始药物剂量水平。该研究在30天后终止，康复赋值1，30天后仍未恢复的患者设置为删失。

加载数据：

use https://www./data/r17/drugtr2

如果只是想了解初始给药水平对康复时间的影响，简单的Cox回归就可以了：

在Stata的帮助文件中，对结果的描述为：The output includes p-values for the tests of the null hypotheses that each regression coeffificient is 0 or, equivalently, that each hazard ratio is 1. That all hazard ratios are apparently close to 1 is a matter of scale; however, we can see that drug number 1 signifificantly increases the risk of being cured and so is an effective drug, whereas drug number 2 is ineffective (given the presence of age and drug number 1 in the model)。这种解释并不恰当。由于每种药物有5个剂量水平（0、50、100、125和150），并作为连续变量纳入模型，以drug1为例，这里的解释应为在考虑了年龄和drug2后，drug1给药剂量每增加1，游走肺炎治愈风险增加1.01倍（P=0.049）。这里的初始剂量没有给出单位，但相对50-150的给药水平，1显然并不大，所以即使drug2系数无统计学意义，也可能是因为给药剂量过低。如探讨药物对生存（治愈）时间的影响，这里设置为分类变量更为合适，感兴趣的可自行检验。

现在我们拟合另外一个模型，想了解下体内药物浓度对生存时间的影响。这里假设两种药物浓度在体内成指数下降趋势，半衰期均为2天，则体内药物浓度就是一个时变协变量：C_t=C₀*exp(-0.35t) 【由C_t=C₀*exp(-βt)和C₂=(1/2)*C₀=C₀*exp(-2β)可推导出β=0.35】。

操作与主要结果如下：

相应的命令为：

stcox age, tvc(drug1 drug2) texp(exp(-0.35*_t)) nohr

注意drug1和drug2是时变协变量，只在参数tvc()中出现而未出现在自变量中，texp是时变协变量的乘数，最终的模型为：

h(t)=h(t₀)exp(-0.149*age+0.266*drug1*exp(-0.35*_t)+0.183*exp(-0.35*_t)

但这里的系数解释应特别注意，以drug1为例，drug1的血压浓度（drug1*exp(−0.35t)）每增加1，治愈的风险是原来的exp(0.266)=1.30倍（P=0.002）。你可能会问为什么在前面一个模型中drug2没有统计学意义，在第二个中就有了呢，这是因为给药剂量增加1和血药浓度[drug1*exp(−0.35t)]增加1是不一样的。

这里的分析逻辑是先对每条记录在每个失效事件发生的时间点上进行拆分[stsplit, at(failures) ]，然后生成两个时变协变量[generate drug1emt = drug1*exp(-0.35*_t)，generate drug2emt = drug2*exp(-0.35*_t)]，最后进行Cox分析[stcox age drug1emt drug2emt]，感兴趣的可以手动拆分试一下。

当然这是一个理想的示例，只是为展示连续时变协变量的可能情况，真实世界中往往会多次服药，药物浓度也会达到一个相对稳态的情况。

示例3：外在时变协变量效应考察（时依系数法）

在时依系数模型中，[β₁+γ₁g(t)]是协变量x₁的时依系数，其中含有一个时间函数g(t)。该系数有一个不随时间变化的分量β₁，γ₁决定了与β₁的时间依赖性偏差的程度。因此，γ₁=0的检验可以看做是检验x₁的系数是否是不变的，而确认一个系数是否固定也是检验比例风险假设的一种方法。

美国退伍军人管理局进行的一项研究中，晚期肺癌的男性患者接受了标准治疗或试验化疗。数据记录了137名患者的死亡时间，9名患者在死亡前离开了研究，同时记录了每个患者的各种协变量。研究主要目标是评估试验化疗是否有益，次要目标包括分析作为预后变量的协变量【age（年龄，年）、KPS（Karnofsky performance score，卡诺斯基评分）、celltype（1=squamous, 2=smallcell, 3=adeno, 4=large）、diagtime（从诊断到研究开始的时间，月）、prior（前期治疗，0=No，1=Yes）、status（生存状态）、time（从注册开始的随访时间，天）】。

我们以age、KPS、celltype、prior、trt为解释变量建立Cox回归模型（celltype设置为因子），同时用时依系数法考察是否满足比例风险：

统计>>生存分析>>回归模型>>Cox比例风险（PH）模型

跟内在时变协变量的区别是，外在时变协变量不仅出现在[时变协变量]中，同时也在[自变量]中。

stcox age i.celltype KPS prior trt, tvc(age i.celltype KPS prior trt) texp(_t) nohr

结果显示celltype的时依系数不为0，也就是说celltype在时间维度上不是一个固定值，不满足风险比例假定。

重新进行Cox回归，时变协变量只保留celltype，[报告]恢复默认设置。

stcox age i.celltype KPS prior trt, tvc(i.celltype) texp(_t)

结果显示，在校正了其他协变量的影响后，与标准治疗相比，试验化疗的死亡风险更高但无统计学意义（HR=1 .33，P=0.163）。卡诺斯基评分是死亡风险的一个重要的保护因素（HR=0.96，P<0.001）。Celltype是一个外在时依协变量，其效应随着时间的变化而变化。

对于外在的时依协变量，除了直接使用这种时依系数法来描述其随时间变化的效应，分段模型也是常用的处理方法。分段模型其实可以看做是时依系数的一个特例，时依系数法在每一个随访时间都进行分段，而分段模型是只是在某个时间点进行分段。虽然风险比例假定在整个随访期间不成立，但在一个较短的时间段内则可能是成立的，分段模型的逻辑就是把整个生存时间拆分成多个时间段，然后为每一段随访时间拟合一个比例风险模型。示例同《非比例风险的Cox回归模型_分段模型》，90例胃癌患者随机接受单纯化疗（chemo）和联合治疗（化疗加放疗，combin）的生存时间研究，SPSS操作也可参见此文，寻找最佳切割点可也参考此文。

统计>>生存分析>>回归模型>>Cox比例风险（PH）模型

结果显示trt是一个时依协变量。经多次尝试，当在svdays=300时，-LL值最小，最终模型我们以生存时间为300天分段界值（参见《非比例风险的Cox回归模型_分段模型》）。

结果显示，当生存时间t≤300时，放化疗联合治疗发生死亡的风险高于单纯化疗（HR=3.668，P=0.002），即在短期内（≤300天）放化疗联合治疗发生死亡的风险是单纯化疗的3.688倍，具有统计学意义。当生存时间t＞300时，放化疗联合治疗引起死亡的作用在下降，生存时间仅为t≤300天的0.149倍，也就是说生存时间t＞300时，放化疗联合治疗死亡风险低于单纯化疗，风险比为0.547（3.668*0.149=0.547）。

需要注意的是这里的P<0.001是针对0.149是否等于1的检验，而不是针对的0.547。实际上我们只能从上述的结果中推断出生存时间t＞300时的HR等于0.547，而不能推断出其是否具有统计学意义。

当然，如果你需要，可以先构建两个模型。短期模型包含所有的患者，其中生存时间≤300天的患者生存时间和结局状态都不变，而生存时间＞300天的患者生存时间都会重新定义成300天，状态改为0（即删失），相当于生存时间＞300天的患者提供的是生存时间为300天的失访信息。用这样一份新的资料来进行COX回归；长期模型只选用生存时间＞300天的患者进行COX回归，其生存时间和结局状态不变，这样就可以得到长期模型HR的P值。两个模型的PH假定也可分别进行检验了（参见《非比例风险的Cox回归模型_分段模型》）。