【原】生存分析之Cox回归

转自个人微信公众号【Memo_Cleon】的统计学习笔记：生存分析之Cox回归。

随访资料的生存分析是一个很大的题目。

从分析的因素上看，有单因素分析和多因素分析。正如“连续资料的单因素分析常用t检验、方差分析，对应的多因素分析是多重线性回归”、“分类资料的单因素分析方法卡方分析，对应的多因素分析有logistic回归”一样，生存分析的常用单因素（或少数因素）的分析有Life Tables法、Kaplan-Meier法，对应的多因素模型则常用Cox回归模型（Cox风险比例模型）。

从采取的分析方法上看，生存分析有非参数法（如Wilcoxon法、Log-rank法）、参数法（如Weibull回归、lognormal回归等）和半参数分析（Cox回归）。

Cox回归要求满足比例风险假定(proportional-hazards assumption)的前提条件。所谓比例风险假定，就是假定风险比(HR，Hazard Ratio)不随时间t变化而变化。

在进行生存分析前，你最好对以下的一些概念及其意义有所了解：起始事件、失效事件(Failure Event)/终点事件(Endpoint Event)、生存时间(Survival Time)/失效时间(Failure Time)、中位生存时间(Median Survival Time)、平均生存时间(Mean Survival Time)、删失值/截尾值(censored values)、生存概率(Survival Probability)、生存率(Survival Rate)/积累生存概率/生存函数/积累生存函数(Cumulative Survival Function)、风险函数(Hazard Function)、累积风险函数……风险函数h(t)=概率密度函数f(t)/生存函数S(t)，概率密度函数f(t)为累积分布函数F(t)的导数，而F(t)=1-S(t)。可参见《生存分析》。

模型结构与参数释义可参见颜虹等主编的《医学统计学》，如下。对此不感兴趣而只关心操作和结果解读的，可直接越过。

当前笔记用STATA演示Cox回归操作。STATA在进行Cox回归分析前首先需要声明生存时间变量，另外比例风险假定是进行Cox回归的前提条件，需要进行考察和检验。

示例（陈启光等.医学统计学第3版）：探讨某肿瘤的预后，某研究机构收集了41例患者的生存时间（月份）、生存结构及影响因素。影响因素包括性别、年龄、病理分级、是否复发、PD-L1分子。变量赋值与资料如下：

X1:gender; X2:age; X3:grade; X4:recur; X5:PD_L1

【1】数据录入：略

【2】声明时间变量：统计>>生存分析>>模型设定和实用工具>>声明数据集为生存时间数据，将在[时间变量]和[失效变量]中选择相应的变量即可。具体步骤如下：

相应命令如下：

stset time, failure(status==1)

也可以不指定失效值，默认不等于0的为失效值。另外scale选项可以重新定义生存时间，如本例生存时间单位是月，可以使用scale(12)后生存时间代表的就是年啦，命令为：stset time, failure(status) scale(12)

【3】Cox回归：统计>>生存分析>>回归模型>>Cox比例风险模型，选入相应的自变量即可。主要添加自变量时应选择合适的变量类型即参照水平（默认低水平为参照），如果未经步骤【2】的生存时间变量声明，也可以在此对话框中的[生存设置]中声明。需要说明的是，本例有序变量按连续变量处理，性别、病例分级、复发即PD-L1虽然是分类变量但都是二分类，直接按连续变量分析也不影响结果。命令为：stcox i.gender age i.grade i.recur i.PD_L1 或者stcox gender age grade recur PD_L1。

结果显示相比空白模型，纳入五个变量的整体模型是有统计学意义的（LR Chi2=22.3，P=0.0005）。

主要结果中默认显示的是风险比（HR），如果想显示模型系数β，需要在Cox比例风险模型对话框中的[报告]选项卡中选中复选框[报告系数，而不是风险比]，HR=exp(β)。根据系数可得出模型为：

h(t)=h₀(t)exp(-0.113·gender+0.365·age-5.44·grade+1.985·recur+0.190·PD-L1)

结果显示校正其他因素，女性相比男性、病例IV期相比III期死亡风险更低，但没有统计学意义；而年龄每增加一个等级、复发相对不复发、PD-L1阳性相比阴性死亡风险更高（分别是1.44、7.28、1.21倍），但只有复发与否具有统计学意义。

很明显，模型给出是强制纳入了所有的变量的结果，但实际上有很多因素并不具有统计学意义，为了精简模型，我们可能需要对模型的变量进行筛选，可以采用逐步回归，SPSS里面可直接在[method]中进行选择相应的方法。STATA里面则可借助stepwise命令，Cox逐步回归菜单操作：

统计>>其他>>逐步估计，具体操作如下。本例采用向前逐步回归，默认wald检验，变量剔除标准P值为0.1，纳入标准P值为0.05。如果直接进行逐步回归，分析前不要忘记声明生存时间变量。

相应命令：stepwise, pr(0.1) pe(0.05) forward: stcox gender age grade recur PD_L1

说明：

stepwise, pr(0.1) pe(0.05) forward : stcox (gender) (age) (grade) (recur) (PD_L1), nohr

pr(剔除标准)，pe(纳入标准)。逐步回归方法：forward、backward，逐步回归方法：lr，默认是wald法。括号内为同进同出的变量为一个回归项，类似于SPSS里面的Block，同一个回归项中变量同进同出，比如同一个分类变量设置为哑变量后可以放在同一个回归项中。nohr表示报告系数而不是风险比。

逐步回归结果如下：

最终引入的变量是recur和age，模型：h(t)=h₀(t)exp(1.756·recur+0.451·age)

recur和age的HR分别为5.791、1.569，即是否复发与年龄是该肿瘤的死亡风险因素。固定其他因素的影响，患者肿瘤复发的死亡风险是不复发的5.791倍；固定其他因素的影响，患者年龄每增加一个等级，死亡风险增加1.569倍。

有时候我们还想得到生存曲线。以指定协变量recur，绘制生存率曲线为例，操作如下：

[图形>>生存分析图>>生存，风险，累积风险或累积发生函数]或者[统计>>后验估计]，在打开的[后验估计选择器]中依次选择：设定，诊断和拟合优度分析>>生存函数，风险函数，累积风险等图形，点击[开始]进入[曲线-绘制生存函数，风险函数，累积风险函数或累积发病率函数]对话框，界面操作如下：

相应命令为：stcurve, survival at1( recur=0 ) at2( recur=1 )

结果如下（你要是觉得图片丑可以启用图形编辑器进行编辑）：

你也可以绘制Kaplan–Meier生存函数图、Kaplan–Meier失效函数图、Nelson-Aalen累积风险函数图即平滑危险估计图等，可在[统计>>生存分析>>回归模型>>图形] 或 [图形>>生存分析图] 中进行。 

【4】比例风险（PH）假定检验：比例风险假定是Cox回归的前提条件，可以通过计算检验，也可以通过图示法。

4.1计算检验命令：整体检验estat phtest；单独检验每个协变量estat phtest, detail

菜单可通过以下途径进入

①统计>>生存分析>>回归模型>>比例风险假设检验；

②图形>>生存分析图>>stcox后比例风险假设检验；

③统计>>后验估计，在打开的[后验估计选择器]中选择[比例风险假定的检验]；

点击[开始]进入对话框[estat-后验估计统计量]，选择[基于Schoenfeld残差的比例风险假设检验(phtest)]。如想进一步检验每个协变量的比例风险假定，可在继续点击对话框[estat-后验估计统计量]中的[选项]按钮，在打开的对话框中选择复选框[单独检验每个协变量的比例风险假设]。

结果显示P>0.05，满足比例风险假设。

4.2 图示法可以通过生存函数的双对数图(Log-log plot of survival)与Kaplan–Meier与Cox预测的生存曲线图(Kaplan–Meier and predicted survival plot)

stphplot plots sln{−ln(survival)} curves for each category of a nominal or ordinal covariate versus ln(analysis time). These are often referred to as “log-log” plots. Optionally, these estimates can be adjusted for covariates. The proportional-hazards assumption is not violated when the curves are parallel.

stcoxkm plots Kaplan–Meier observed survival curves and compares them with the Cox predicted curves for the same variable. The closer the observed values are to the predicted, the less likely it is that the proportional-hazards assumption has been violated. 

4.2.1 Log-log plot of survival可以拟合命令单独或者分层的Cox模型，并可根据调整变量进行调整。

命令1：stphplot, by(recur)

根据分类变量recur拟合单独的Cox模型，结果是下图左；

命令2：stphplot, strata(recur) adjust(age) zero

根据分类变量recur拟合分层的Cox模型，并根据变量age进行调整，调增方法是将age值调整为0（默认是均值），结果为下图右。

菜单可通过以下途径进入stphplot-生存函数的双对数图(Log-log plot of survival)

①统计>>生存分析>>回归模型>>比例风险假设的图形评估；

②图形>>生存分析图>>比例风险假设检验；

选入相应的自变量和分层变量即可。

结果显示变量recur两个水平基本“平行”，满足比例风险的假定。

4.2.2 Kaplan–Meier and predicted survival plot

命令：stcoxkm, by(recur)

当然你可以单独绘制recur两个水平的观测预测曲线：stcoxkm, by(recur) separate

菜单操作：

①统计>>生存分析>>回归模型>>Kaplan–Meier生存曲线和Cox预测曲线比较；

②图形>>生存分析图>>Kaplan–Meier和Cox生存曲线比较；

选入相应的自变量即可。

结果显示Kaplan–Meier观测曲线与Cox回归预测曲线重合性较好，满足比例风险假设。

如果运气不好，你的数据违背了比例风险的假定，可以考虑含时依协变量的Cox回归。含时依协变量的Cox回归也可以作为验证比例风险的假设的一种手段。个人理解，这个所谓时依协变量，其实就是在Cox回归里构建的一个交互作用项。这个关于含时依协变量的Cox回归我们放在以后分享。