生存资料回归模型分析——基于MCMC过程构建生存资料Cox比例风险回归模型

刘媛媛 ，李长平 ，2，胡良平

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

频率学派和贝叶斯学派是统计学发展历史上两个重要的学派［1］。通常，前者认为随机事件的概率是客观存在并假设固定不变的；而后者则认为此“概率”是随机的而不是固定不变的，并服从于某种分布。也即，经典统计学认为“参数”是固定不变的“常量”，而贝叶斯统计学认为“参数”是随机变量，这是两者的根本分歧所在。1763年由Richard Price整理发表了贝叶斯的成果《An Essay towards solving a Problem in the Doctrine of Chances》提出了贝叶斯公式，并介绍了贝叶斯思想，其核心内容就是对参数的估计并不是单纯取决于客观数据，而是取决于客观数据（包括总体和样本信息）和先验信息的共同作用［2］。随着计算机技术的发展和贝叶斯方法的改进，特别是马尔科夫蒙特卡洛（Markov chain Monte Carlo，MCMC）方法的提出和应用，使得参数后验分布的模拟得以更方便地实现，从而体现出该法在处理小样本数据时的优势［3］。现在，越来越多的新的统计分析方法将经典统计分析方法和贝叶斯思想有机地结合起来，例如，基于贝叶斯理论和生存分析相结合的贝叶斯生存分析在近年来越来越多地被应用于不同的研究领域，尤其是医学科学研究中［4-5］。因此，本文将介绍基于PHREG过程和MCMC过程分别构建贝叶斯统计思想框架下生存资料的Cox比例风险回归模型的相关内容。

1 基于“贝叶斯统计思想”构建Cox比例风险回归模型

姚婷婷等［6］的文章已经介绍了Cox比例风险回归模型，见式（1）：

式（1）中，X1、X2、…、Xp为与生存时间可能有关的自变量（即影响因素）；h(t)为具有自变量X1、X2、…、Xp的个体在t时刻的风险函数；h0(t)为所有自变量为0时t时刻的基准风险函数；β1、β2、…、βp为各自变量的偏回归系数。偏回归系数的估计需借助偏似然函数，用最大似然估计方法得到。偏似然函数的计算公式见式（2）：

式（2）中，qi为第i个死亡时点的条件死亡概率，其分子部分为第i个个体在ti（t1≤t2≤…≤ti≤…≤tk）死亡时点的风险函数h(ti)，分母部分为生存时间T≥ti的所有个体（包括死亡和删失）的风险函数之和。

基于贝叶斯统计思想构建生存资料回归模型，即在原模型的基础上利用贝叶斯方法的基本原理对回归参数进行估计的过程。所以，需要先对这些参数指定适当的先验分布，如果先验分布选择不合适，则会对结果产生影响。故文献［7］建议将Cox回归模型系数βi的先验分布设定为正态分布，本研究也将按此进行先验分布的设置。

2 构建Cox比例风险回归模型

2.1 实例与数据

本文所用数据来自一项骨髓瘤研究，研究者用烷基化剂治疗65例患者，在随访期间内，死亡48例，存活17例。变量赋值情况见表1。

表1 变量赋值表

【说明】完整数据来自文献［8］。因篇幅所限，此处未呈现全部数据。

2.2 采用“PHREG过程”且基于贝叶斯统计思想构建Cox比例风险回归模型

可以利用PHREG过程的BAYES语句拟合Cox比例风险回归模型。

SAS程序如下：

【程序说明】MODEL语句是PHREG过程的必需语句，等号左边定义生存时间和生存结局变量（括号内为截尾数据标识），右边为各协变量（即自变量）。使用BAYES语句则要求回归模型的贝叶斯分析是采用Gibbs抽样，同时设定seed为随机数生成器种子，设为1；NMC为退火（退火是指为了使初始值对后验推断的影响最小化，需要在Markov Chain达到目标分布后弃掉先前的部分样本）后的迭代次数，设为10000；OUTPOST选项将后验分布样本保存在SAS数据集中以进行以后的处理。

【主要输出结果及解释】

这是输出结果的“第1部分”，第1列“参数”实际上是拟创建的回归模型中的“自变量”；第2列指随机重复抽样一万次；第3列“均值”实际上是各自变量前的回归系数的估计值，而且，其中每个估计值都是一万次随机重复抽样计算所得到结果的算术平均值；第4列为与“各均值”对应的“标准差”；最后两列为与“各均值”对应的95%HPD（highest posterior density，HPD）区间，即95%最高后验密度置信区间。因此，根据此置信区间是否包含“0”（包含0时表明该变量对结果的影响无统计学意义），可得以下回归方程：

h(t)=h0(t)exp(1.7610×LogBUN)

图1 变量LogBUN回归参数的马尔可夫链迭代轨迹图

图2 变量LogBUN回归参数的自相关函数图

图1显示马尔可夫链已经收敛（因篇幅所限，其他协变量的相关结果此处从略）。

2.3 采用“MCMC过程”且基于贝叶斯统计思想构建Cox比例风险回归模型

图3 变量LogBUN回归参数的后验密度核密度图

2.3.1 利用LAG函数拟合模型

SAS程序如下：

【程序说明】ARRAY语句用于将回归系数的名称与协变量、常量相关联。PARMS语句给出模型中的参数名称，并为其指定初始值。PRIOR语句指定模型参数的先验分布为正态分布。程序中的bZ为似然函数中的回归项即每个观测的风险集项。本例所用似然函数为Breslow似然函数。符号“l”为每个观测计算对数似然的公式。IF-ELSE语句将所有的观测分成三部分，并使用lag1函数来检验两个相邻的生存时间time是否不同。符号“v”为生存状态（Vstatus）的总和，因为删失数据不进入似然公式的计算，所以需要将其去掉。MODEL语句用于在给定参数的似然函数的情况下指定数据的条件分布。

【主要输出结果及解释】

这里的输出结果（特指第3列到第6列）与前面输出结果的第1部分（特指第3列到第6列）是基本相同的，各列的含义相同，此处从略。

其中beta1～beta9分别对应各协变量，实际上就是前面输出结果中第1部分的“第1列”。

2.3.2 利用JOINTMODEL选项拟合模型

若利用PROC语句中的JOINTMODEL选项，则可使对数似然函数适用于整个数据集，而不只是单个的观察值。但在使用此选项前，还要为数据风险集S指定包含其中的观察的数量，为此需先创建一个stop变量。因篇幅所限，这部分的SAS程序和输出结果从略。

3 讨论与小结

3.1 讨论

以Cox比例风险回归模型为例，对于满足PH假定的有删失数据的生存资料来说，该模型能够很好地利用偏似然（Partial Likelihood，PL）估计理论识别响应变量的影响因素，但其应用仍需一定的样本量。贝叶斯思想的提出则有效地弥补了小样本的不足。二者相结合的基于贝叶斯统计思想的Cox比例风险回归模型很好地融合了其各自优势。

本研究通过实例，基于贝叶斯统计思想并借助PHREG过程和MCMC过程分别构建了Cox比例风险回归模型，并且分别利用MCMC过程中的LAG函数、JOINTMODEL选项拟合模型。通过对程序及结果的解释比较，不难发现此三种方法均可得到后验样本统计描述指标（包括均值、标准差及95%最大后验密度置信区间）、有效样本大小、马尔可夫链迭代轨迹图等主要结果，而且后验样本的均值等指标在数值上相差不大，MCMC过程中两种方法的结果更是完全一致，只是在结果展示形式上稍有不同。由于PHREG过程本身是实现经典统计学中Cox比例风险模型回归分析的标准过程，这里只需加入BAYES语句用于指定进行贝叶斯估计。因此，通过比较不同过程的具体语句可以发现，PHREG过程相对MCMC过程要简洁，并且可以提供普通的Cox比例风险模型的结果——回归系数的最大似然估计的结果（包括回归系数估计值、标准误差及其95%置信区间），而MCMC过程只是提供了“平均意义”下的类似结果。

3.2 小结

在SAS中，可以借助PHREG过程或MCMC过程构建基于贝叶斯统计思想的Cox比例风险回归模型，两种做法的主要结果相似，但前者程序语句相对简洁。研究者可根据具体情况选择其一进行回归模型的构建。