Cox比例风险模型是由英国统计学家Cox^[1]于1972年提出的一种半参数模型。该模型用风险函数反映协变量对生存期的影响，能够解决截尾数据的问题，并且可以同时分析多因素对生存期的影响，在疾病预后分析中得到了广泛的应用^[2-3]。应用Cox比例风险模型要求资料满足两个前提假设，一是个体间的风险比（hazard ratio，HR）恒定，即模型中协变量的效应不随时间改变而改变；二是对数风险或对数累积风险与协变量间的关系为线性。在实际研究中，传统Cox比例风险模型的两个前提假设很难同时满足。当不满足线性关系时，研究者常对连续型变量进行分段转化为分类变量，但分类变量的类别数目以及节点位置的选择一般会带有主观性并且会损失部分信息^[4]。此外，研究者也可以通过构建多项式回归或样条回归来直接拟合自变量和因变量之间的非线性关系。然而，多项式回归使每个段的内部效应被强制统一，在节点位置跳跃，“瞬时变化”不合理，存在过度拟合和共线性的问题。样条回归本质上也是多项式回归，但一般要求在每个节点上连续且二阶可导^[5]，以保证曲线的平滑性。

限制性立方样条（restricted cubic spline，RCS）^[6]是在样条回归的基础上再加一个约束条件，即样条函数在自变量数据范围两端的两个区间内为线性函数，使得两端的预测更为准确。当连续型变量无法满足传统Cox比例风险模型的线性假设时，限制性立方样条Cox比例风险模型是分析非线性关系最常见的方法之一。

既往许多研究探讨了限制性立方样条在非线性回归中的应用^[2-3]，但鲜少进行模型评估。本研究将从基本原理、前提假设等方面介绍限制性立方样条Cox比例风险模型与传统Cox比例风险模型，并结合乳腺癌预后实例评估其预测效果，探讨限制性立方样条Cox比例风险模型在肿瘤预后分析中的应用。

材料和方法

基本原理  假定有$ n $个观测，每个观测$ i $有3个观测值，分别是$ {t}_{i}、{x}_{i}、{c}_{i} $，其中$ {t}_{i} $为生存时间，$ {x}_{i} $为预测变量$ , {c}_{i} $为截尾指示变量（非截尾观测：$ {c}_{i} $=1；截尾观测：$ {c}_{i} $=0）。Cox比例风险模型的一般形式为：

$ h\left(t, {x}_{i}\right)={h}_{0}\left(t\right){\mathrm{e}}{\mathrm{x}}{\mathrm{p}}\left(\beta {x}_{i}\right) $

式中：$ h\left(t, {x}_{i}\right) $是指在危险因素$ x $的影响下，$ t $时刻的风险函数（hazard function），$ {h}_{0}\left(t\right) $是指当所有自变量$ {x}_{i} $都为0时，$ t $时刻的基准风险函数，它是与时间有关的任意函数，函数形式无任何限定，$ \beta $是回归系数。

在应用Cox比例风险模型进行统计推断和预测前，必须考察生存资料是否满足两个基本假定：（1）比例风险假定：任意两个风险函数之比不随时间的改变而改变，在协变量不同状态下，个体的风险比在不同时间点为常数，即为“等比例风险（proportional hazards）”。其中任意两个风险函数之比为相对危险度RR或风险比（hazard ratio，HR）：

$ HR=\frac{{h}_{i}(t, {x}_{i})}{{h}_{j}(t, {x}_{j})}=\frac{{h}_{0}\left(t\right){\mathrm{e}}{\mathrm{x}}{\mathrm{p}}{\mathrm{ }}({\beta }_{1}{x}_{i1}+{\beta }_{2}{x}_{i2}+\cdots +{\beta }_{p}{x}_{ip})}{{h}_{0}\left(t\right){\mathrm{e}}{\mathrm{x}}{\mathrm{p}}{\mathrm{ }}({\beta }_{1}{x}_{j1}+{\beta }_{2}{x}_{j2}+\cdots +{\beta }_{p}{x}_{jp})} $

（2）对数线性假定：对数风险比应与模型中的连续型协变量呈线性关系：

$ log{h}_{i}\left(t\right)-log{h}_{0}\left(t\right)=\beta {x}_{i}, i={\mathrm{1, 2}}, \cdots n $

样条函数是由具有某些连续性条件的子空间上的分段多项式构成，一般要求每个分段点上连续且二阶可导，以保证曲线的平滑性。假如给定$ n+1 $个点$ {t}_{0} $，$ {t}_{1} $，$ \cdots {t}_{n} $满足$ {t}_{0} < {t}_{1} < \cdots < {t}_{n} $，则这些点称为节点（knots）^[7]。具有$ n+1 $个节点$ {t}_{0} $，$ {t}_{1} $，$ \cdots {t}_{n} $的一个$ k $次样条函数（$ k\ge 0 $）满足条件：（1）在每一个区间$ [{t}_{i-1}, {t}_{i}] $上，S是一个次数$ \le k $的多项式；（2）在$ [{t}_{0}, {t}_{n}] $上S有（$ k $-1）阶连续导数。

由于样条函数在第一个节点前和最后一个节点后具有不稳定性，为减小“尾部效应”，限定函数在尾部为线性^[8]，即样条函数在自变量范围两端的两个区间$ [{t}_{0}, {t}_{1}) $和$ [{t}_{k-1}, {t}_{k}) $内是线性函数。具有$ k $个节点的限制性立方样条函数表达式为：

$ RCS\left(x, k\right)=\sum\limits _{i=1}^{k-1}{\beta }_{i}{S}_{i}\left(x\right) $

其中：

$ {S}_{1}\left(x\right)=x, $

$ {S}_{2}\left(x\right)=(x-{t}_{1}{)}_{+}^{3}-\frac{{\left(x-{t}_{k-1}\right)}_{+}^{3}\left({t}_{k}-{t}_{1}\right)}{{t}_{k}-{t}_{k-1}}+\frac{{\left(x-{t}_{k}\right)}_{+}^{3}\left({t}_{k-1}-{t}_{1}\right)}{{t}_{k}-{t}_{k-1}} $

$ {S}_{i}\left(x\right)={\left(x-{t}_{i-1}\right)}_{+}^{3}-\frac{{\left(x-{t}_{k-1}\right)}_{+}^{3}\left({t}_{k}-{t}_{i-1}\right)}{{t}_{k}-{t}_{k-1}}+\frac{{\left(x-{t}_{k}\right)}_{+}^{3}\left({t}_{k-1}-{t}_{i-1}\right)}{{t}_{k}-{t}_{k-1}}, i={\mathrm{1, 2}}, \cdots k-1 $

其中：

$ (x-{t}_{j}{)}_{+}^{3}=\left\{\begin{array}{l}(x-{t}_{j}{)}^{3}, \;if\;x\ge {t}_{j}\\ 0, \;\;\;\;\;{\mathrm{其}}{\mathrm{他}}\end{array}\right. $

对于每个样条变量x需要估计的参数个数是$ k $-1个，其中非线性部分$ k $-2个。限制性立方样条Cox比例风险模型若为时间三节点，则会产生一个新的时协变量$ {S}_{i}\left(t\right) $，其中$ t $表示时间，新的风险函数表达式为：

$ h\left(t, x\right)={h}_{0}\left(t\right){\mathrm{e}}{\mathrm{x}}{\mathrm{p}}{\mathrm{ }}[{\beta }_{1}{x}_{1}+{\beta }_{2}{x}_{1}t+{\beta }_{3}{x}_{1}{S}_{1}\left(t\right)] $

$ {x}_{1} $的单位变化，相应的对数风险比函数（log hazards ratio function，LHRF）为：

$ LHRF={\beta }_{1}+{\beta }_{2}t+{\beta }_{3}{S}_{1}\left(t\right) $

此时，对数线性假定检验为$ {\beta }_{3} $＝0是否成立，比例风险假定检验为$ {\beta }_{2}={\beta }_{3}=0 $是否成立。

鞅残差（martingale residuals）为研究中存在的但未被模型预测到的部分事件数的估计值，即个体观察事件数与模型给定的“期望”事件数之差，可用于检验Cox模型的对数线性假定。Bootstrap是非参数统计中一种重要的估计统计量方差的统计学方法，可检验预测模型的准确度。基本思想是：针对相应样本资料采用有放回的大量重复随机抽样，对每次随机抽样所获得一个Bootstrap样本计算某个特殊指标的统计量的值，并把大量的Bootstrap样本的统计量计算值视为新的数据进行统计描述和统计推断。

一致性指数（concordance index）也称为C-index，最早由范德堡大学生物统计教授Harrell^[9]在1996年提出，本质上是估计预测结果与实际观察的结果相一致的概率，主要反映预测模型的区分能力，常用于评价肿瘤预后模型的预测精度，其在训练数据集较小而在验证数据集较大，是模型过拟合的重要表现之一。净重分类改善指标（net reclassification improvement，NRI）评价的是在两种模型采用最优诊断截点进行预测时，使用新模型使得个体的预测结果得到改善的概率，包括事件发生组改善概率与不发生组改善概率两部分；若NRI > 0，则为正改善，说明新模型比旧模型的预测能力有所改善。综合判别指数（integrated discrimination improvement，IDI）反映的是两种模型预测概率差距上的变化，是基于疾病模型对每个个体的预测概率计算所得，可用来反映模型的整体改善状况；若IDI指数 > 0，提示新模型比旧模型的预测能力有所改善。校准度（calibration）^[10]指模型预测的事件发生概率与观察到的事件发生概率的一致程度，反映了模型正确估计绝对风险的程度。

实例分析  乳腺癌数据资料来源于Royston等^[11]用于验证模型的GBSG数据，包括686名德国原发性淋巴结阳性乳腺癌患者的预后变量信息。数据集为R4.0.2中内置数据集gbsg，包括的变量有：初次手术年龄（age，岁），绝经状态（meno，1=绝经后，0=绝经前），肿瘤大小（size，mm），肿瘤分级（grade，3=Ⅲ级，2=Ⅱ级，1=Ⅰ级），阳性淋巴结数（nodes），孕激素受体（pgr，fmol/L），雌激素受体（er，fmol/L），荷尔蒙疗法（hormon，1=有，0=无），无复发生存时间（rfstime，天），生存状态（status，1=存活且无复发，0=复发或死亡）。利用R4.0.2中的survival包拟合传统Cox比例风险模型；利用R4.0.2中的rms包拟合限制性立方样条Cox比例风险模型；利用R4.0.2中的pec包对比两种模型的C-index与校准度。利用R4.0.2中的survIDINRI包对比两种模型的NRI与IDI。P < 0.05为差异有统计学意义。

结果

传统Cox比例风险模型  将初次手术年龄、肿瘤大小、阳性淋巴结数，孕激素受体、雌激素受体等5个连续型变量转换为分类变量后进行分析。单因素分析结果显示绝经状态与生存状态无显著关联（P > 0.05），故纳入初次手术年龄、肿瘤大小、肿瘤分级、阳性淋巴结数，孕激素受体、雌激素受体、荷尔蒙疗法构建传统多因素Cox比例风险模型，回归分析结果见表 1。

限制性立方样条Cox比例风险模型  根据martingale residuals判断满足非线性关系的变量有初次手术年龄、阳性淋巴结数、孕激素受体与雌激素受体。R²与Dxy值越大，则拟合的模型越优。表 2列出了不同节点数的拟合效果，根据R²与Dxy判断初次手术年龄和孕激素受体最佳节点数为4，阳性淋巴结数和雌激素受体最佳节点数为5。纳入初次手术年龄、肿瘤大小、肿瘤分级等因素构建限制性立方样条Cox比例风险模型，多因素回归分析结果见表 3。

模型评估  采用Bootstrap重抽样法评估C-index（图 1），预测时间小于511天时，传统Cox比例风险模型的区分度优于限制性立方样条Cox比例风险模型；当预测时间大于511天时，限制性立方样条Cox比例风险模型的区分度优于传统Cox比例风险模型。进一步比较模型的NRI为0.179（95%CI：0.065~0.320），IDI为0.041（95%CI：0.018~0.083），表明限制性立方样条Cox比例风险模型比传统Cox比例风险模型的预测能力有所改善。采用bootstrap重抽样法评估校准度，绘制校准曲线直观展示预测风险值与实际风险值的关系。限制性立方样条Cox比例风险模型的校准度优于传统Cox比例风险模型（图 2）。传统Cox比例风险模型训练集与验证集的C-index分别为0.70和0.71，限制性立方样条Cox比例风险模型分别为0.71和0.74，两种模型均未出现过度拟合现象。

讨论

本研究采用限制性立方样条Cox比例风险模型分析了初次手术年龄、阳性淋巴结数、孕激素受体及雌激素受体与乳腺癌复发之间的非线性关联，进一步使用C-index、NRI、IDI以及校准度对限制性立方样条Cox比例风险模型与传统Cox比例风险模型进行评估，证明了限制性立方样条进行肿瘤预后分析的应用价值。由于乳腺癌预后资料缺少淋巴管浸润、肿瘤分期以及HER2表达数据，预后模型不包括以上3项必要的预测参数，可能会影响预后模型本身的预测效力。

建立模型时，如何处理模型中的连续型预测因子将影响模型的区分度与校准度。通常对连续型预测因子的处理方法有：（1）使用一个或多个切入点对连续型预测因子进行二分类或多分类；（2）假定预测因子与结果呈线性关系，对连续型预测因子直接进行建模；（3）预测因子保持连续，使用分数多项式或限制性立方样条对非线性关系进行建模。Collins等^[12]曾比较上述3种方法对预测模型性能的影响。结果表明，分数多项式或限制性立方样条模型的区分度和校准度明显优于其他方法；而对连续型预测因子进行分类会损失大量的方差信息，同时节点附近的跳跃使模型拟合变差，预测性能明显降低。Nieboer等^[13]的研究亦表明，当连续型预测因子与结局为非线性关系时，使用限制性立方样条可以更准确地描述其与结局发生的关系。既往研究基于限制性立方样条Cox比例风险模型分析BMI与全因死亡风险之间的非线性关联^[3]，弥补了传统Cox比例风险模型对变量的暴露-效应曲线应用条件限制。在本研究中，我们采用限制性立方样条同时探讨多个存在非线性关联的自变量，也取得了较理想的预测效果。

应用限制性立方样条拟合非线性关系时，需要设置样条函数的节点个数及位置进行分段回归。Stone等^[7]研究发现，大多数情况下节点的位置对限制性立方样条的拟合效果影响较小，而节点的个数会决定拟合曲线的形状。当节点个数等于样本量时，相当于将各个点用线段相连，得到的是完全拟合但不平滑的折线。由于节点个数的选择和自由度有关，所以当样本量较大时可以取较多的节点，但存在的问题是节点越多，自由度越大，模型越复杂，结果越难解释^[5]，且样条估计量的方差和过度拟合的风险会增加。限制性立方样条Cox比例风险模型中防止过度拟合的常见方法：通过增加训练数据更好地识别信号，避免噪声；基于交叉验证生成多个训练集，测试划分并调整模型；根据样本量控制节点个数，降低模型的复杂性。对于大部分数据集，4个节点可使模型较好拟合，既可以兼顾曲线的平滑程度，又可避免过度拟合造成的精确度降低。如果样本量较小，则建议使用3个节点，以便节点之间有足够的观测值能够拟合每个多项式。当样本量较大时，5个节点是更好的选择。本研究中，限制性立方样条Cox比例风险模型使用的最佳节点数为4和5，取得了较好的模型预测效果，且未出现过度拟合现象。

综上所述，限制性立方样条应用多样的平滑且合理的剂量-效应曲线来表达连续性暴露与结果之间的关联，可以解决非对数线性问题，拟合资料具有更大的灵活性。但限制性立方样条也存在一定的限制性：首先，限制性立方样条对节点的数量敏感，无法结合高次项参数背景进行解释^[5]；其次，全变量模型中协变量较多，选择变量的最佳拟合形式工作量很大，在统计软件中编程也较为复杂，容易引起结果偏差。所以当随访数据不满足传统Cox比例风险模型的对数线性假定且影响因素较少时，应用限制性立方样条方法进行肿瘤预后分析更加高效。

作者贡献声明  张彭燕  文献检索，统计分析，论文构思、撰写和修订。刘振球，樊虹，索晨，陈兴栋  结果解释，论文修改。张铁军  论文构思、指导和修订。

利益冲突声明  所有作者均声明不存在利益冲突。