生存资料回归模型分析——Cox比例风险假设的拟合优度法

宋德胜 ，李长平 ，2，刘媛媛 ，崔 壮 *，胡良平

（1.天津医科大学公共卫生学院流行病与卫生统计学教研室，天津 300070；2.世界中医药学会联合会临床科研统计学专业委员会，北京 100029；3.军事科学院研究生院，北京 100850*通信作者：崔 壮，E-mail：cuizhuang@tmu.edu.cn）

在生存分析中，生存函数是指发生事件的时间超过指定时间的概率［1］，这个概率取决于两个方面：基线风险函数和效应参数。前者体现了在基线协变量下，事件发生的风险如何随时间变化；后者表示发生事件的风险如何随协变量变化。因此，若效应参数随时间发生变化，使用经典的Cox比例风险回归模型估计的生存函数会存在较大偏差。鉴于上述情况，对比例风险假设的检验也就成为了应用Cox比例风险模型的重要前提。目前，对于比例风险的检验主要包括图示法（前文已阐述，本文不再赘述）和拟合优度检验法。前者虽然简便直观，但是存在主观性。因此，本文从拟合优度的角度来介绍比例风险假设的检验方法。

1 原理阐述

Grambsch和Therneau证明了缩放Schoenfeld残差在诊断比例风险假设方面的作用［2］。在比例风险假设成立的条件下，每个模型中的协变量的缩放Schoenfeld残差与生存时间的函数斜率为0。基于此，可以构建检验统计量。

对于单个预测变量（即自变量）的检验统计量见（1）：

式（1）中，rs是缩放Schoenfeld残差，g(t)是检验之前预先指定的时间函数形式，δi是事件的指示变量，Δ是事件总数，是目标预测变量的参数估计值的方差。该统计量近似服从自由度为1的χ2分布。

对于p个预测变量的总体检验的统计量见式（2）：

式（2）中，ri是目标预测变量未缩放的Schoenfeld残差向量。该统计量服从自由度为p的渐近χ2分布［3］。

除了上述检验方法外，也可以在SAS软件RHREG过程的MODEL语句中加入协变量与时间函数的交互项，使用Wald检验来判断协变量是否符合比例风险假设。

2 SAS软件实现

2.1 示例数据说明

本文使用的数据集Myeloma来源于多发性骨髓瘤的研究数据［4］。原始数据、变量说明和SAS数据集等内容参见前文，此处从略。

2.2 比例风险假设检验

以下程序演示了如何通过SAS软件的PHREG过程进行Schoenfeld残差检验以及如何向模型中加入自变量与时间函数交互项检验比例风险假设。

【程序说明】①是基于加权Schoenfeld残差进行检验的SAS程序。Proc phreg表示调用PHREG过程，Data选项指定要分析的数据集为myeloma，后面的ZPH选项（注意：此选项在SAS 9.3中无效）要求进行比例风险假设检验，使用的是加权Schoenfeld残差；ZPH选项中的transform指定时间的变换形式，本例使用的是对数变换；Class语句指定platelet和frac是分类变量；Model语句指定要拟合的模型变量；等号左边time*vstatus（0）中的time是生存时间，vstatus指定截尾指示变量，括号中的0表示vstatus=0为截尾；等号右边指定要分析的自变量，Selection=s表示使用逐步选择进行变量筛选。②是在Cox回归模型中加入自变量与时间函数的交互项来进行检验，语句和①基本一致。不同之处是在MODEL语句中加入了自变量与时间函数的交互项，即Log-BUNt、HGBt、Aget、LogWBCt、LogPBMt、Proteint 和Scalct。MODEL语句后面指定交互项的具体计算式。

输出结果分别见表1、表2、表3。

表1 Cox比例风险回归模型中参数的最大似然估计结果

表2 ZPH选项给出的比例风险假设检验结果

表3 Cox非比例风险回归模型中参数的最大似然估计结果

表1是使用逐步回归进行变量筛选后的结果。可以看到最终模型中保留的变量为LogBUN和HGB（P＜0.05）；危险比分别为5.336、0.888，它们是基于Cox比例风险回归模型计算得到的。

表2是ZPH选项给出的比例风险假设检验结果。可以看出模型的全局检验结果：卡方统计量为6.6966（P＜0.05），因此，模型中至少存在一个变量不满足比例风险假设。而LogBUN的卡方统计量为6.6837（P＜0.05），即LogBUN不满足比例风险假设。

表3为最终的最大似然估计结果，可以看出，LogBUN和LogBUN与时间对数的交互项均保留在模型中（P＜0.05），因此，可以得出LogBUN不满足比例风险假设。

综合以上两种检验结果，可以看出，LogBUN不满足比例风险假设。

3 讨论与小结

3.1 讨论

残差可以体现构建模型的拟合程度。由于生存分析资料的特殊性，其残差的计算相比一般的参数回归分析方法更为复杂。在生存资料中，常见的残差包括Martingale残差、Deviance残差、Score残差以及Schoenfeld残差。其中，Schoenfeld残差在检验比例风险假设时更为常用。因此，本文介绍了一种基于Schoenfeld残差进行比例风险假设检验的方法。该方法克服了图示法中只能直观判断某变量是否满足比例风险假设的缺点，以构建统计量的形式给出了一种客观的数值度量方式。这种方法在本文中是通过使用SAS软件PHREG过程ZPH选项实现的，但这个选项在SAS/STAT 13.1（即SAS 9.4）以后版本才可用。使用SAS/STAT 13.1以前版本的读者可考虑使用Chen［3］编写的PH_score_test宏程序来实现，或者考虑构建协变量与时间函数的交互项，将之放入MODEL语句中进行拟合的方式来检验比例风险假设是否成立。但这种方法存在一个局限，即无法确定合理的时间函数形式。

在明确存在时间依赖型协变量以后，可以考虑扩展的Cox回归模型来进行后续的模型构建，也可以选用RMST回归［5］来进行拟合。

3.2 小结

本文主要介绍了使用残差检验法和构建协变量与时间函数的交互项、从客观角度进行检验的方法。这两种方法在本文示例中检验的结果一致。但考虑到时间函数形式无固定标准，读者需要根据实际资料灵活选择检验方法。