基于Cox模型的虚拟响应算法，本文主要内容关键词为：算法论文,模型论文,Cox论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
引言
在设备可靠性验证过程中，通常需要综合各种不同环境下的实验数据，从而尽可能多的利用实验数据中获得的信息。传统的可靠性综合验证方法主要利用单一环境下的实验数据，因此如直接套用传统方法，需要对于各种不同环境下的实验数据分别进行分析，这样不但费时费力，而且分析结果不具有直接可比性，为此我们希望能够将各种不同环境下的实验数据放在一个统一的框架内进行研究。一个自然的想法是将环境变量作为协变量引入模型中，从而可以在统一的模型中处理不同环境下的数据。可以引入协变量的模型有很多，比如Cox模型[1]，Aalen可加模型[2]等等，其中应用最为广泛的是Cox模型。
Cox模型自提出以来，受到了研究者的广泛关注，并且由于模型的灵活性在实际应用中取得了巨大的成功，经过数十年的研究，得到了许多针对Cox模型的行之有效的处理方法，Cox模型也从最基本常数回归系数的情况推广到了具有时变效应的情况，更进一步，推广到了Cox-Aalen[2]模型，加入了可加模型的部分。可以说，Cox模型及其推广模型已经成为生存分析中最重要的模型。
Cox模型从本质上来讲是一种半参数模型，传统的处理办法是基于偏似然函数得到回归系数的估计。偏似然方法的好处是不必确定基准生存函数就可以直接得到协变量回归系数的估计。该方法的结果具有较好的性质，但是由于似然函数本身的特点，在求极值的过程中需要用到Newton法或者EM算法等等数值方法，求解过程不直观，当协变量维数较多的时候这样的高维优化问题解决起来比较困难，在现实应用中可操作性较差。此外，传统的偏似然方法是假定了协变量作用的对数线性形式，很多时候这个假定能够满足要求。但是有时候对于我们的特定模型以及能够观测到的协变量，该假定未必成立。当协变量作用的参数形式未知，传统的方法无从谈起。针对这种情况，文献[3]中提出了局部部分似然的方法，在局部利用多项式来逼近，将多项式的系数引入作为参数来最优化。该方法能够处理协变量作用形式为非参的情形，但是由于似然函数本身的特点仍然无法避免用数值方法最优化的问题。
一个直观的想法是，根据Cox模型的结构特征，我们可以根据样本数据得到一个以协变量作为自变量的响应变量，以下称为虚拟响应，然后根据协变量和相应的虚拟响应来拟合出协变量对于基准生存函数影响的形式。该算法的好处就是简单直观，对于协变量的作用形式是参数形式还是非参没有限制。
一、模型与算法
（一）基本Cox模型
Cox在1972年提出了如下模型[1,4]
至于关于协变量的大样本性质，假设检验以及时变性检验的研究，在文献[2]中有较完整的总结。
Efron[8]给出了无节点情况下的精确似然函数，从中可以看出Cox偏似然只是完全似然的一部分，因此又叫作部分似然函数。事实上，偏似然函数是我们的样本的秩统计量的似然函数[10]。虽然偏似然并不是完全的似然，但是许多学者的研究证明[8，11]，偏似然基本上包含原模型绝大部分的信息，而且在一定的正则条件下基于偏似然的估计结果具有相合性和渐近正态性。虽然偏似然的结果具有前述的优良性，但是求解过程本身却不容易。从形式上看，乘积的每一项都是0到1之间的正数，乘积的项数越多，则偏似然越趋向于0。而且我们看到分母上是一系列的加和，导致我们的最优化过程比较困难，当β维数较高时尤其如此。在实际应用中，求解过程不直观。
2.虚拟响应方法
我们把(1)式变形得到：
在上述迭代过程中，我们以两步迭代中回归系数
的变化作为迭代的终止规则：当两步迭代得到的回归系数没有变化（实际中是两步回归系数之差的范数小于某个指定值），则迭代终止。迭代终止时step 2所得结果即为回归系数的估计，step 3结果即为基准生存函数的估计。
（二）推广的Cox棋型
注意到我们的迭代过程中对于
的形式并未作任何的限制，所以我们并没有降低偏似然方法的灵活性。更进一步，我们的算法事实上并未限制协变量对于基准生存函数作用的具体形式。也就是说，我们可以将模型推广到所谓的比例风险模型：
二、数值模拟
算例1 我们以基准生存函数为指数型为例，采用随机截断的方法得到观测数据，同时利用偏似然方法和虚拟响应算法进行基准生存函数的估计以及回归系数的估计。
（一）模拟条件
基准生存函数取为均值为100的指数分布，随机截断c=30+t，t抽自均值30的指数分布。回归系数的真值为
=(0.30，0.21，0.53，0.8)，基准样本量n=[70，62，41，50，37]分别取样本量N=n，5n，20n，40n进行模拟，其中协变量水平数为5，具体为
基准样本量n=(70，62，41，50，37)。当两部迭代之间回归系数变化值的范数小于
时，迭代终止。下文中图形是50次模拟的结果，MSE是500次模拟的结果。
（二）模拟结果
偏似然方法和虚拟响应算法对于回归系数估计的结果比较。
表1
下面的图形分别为基准生存函数的拟合情况（图1），各样本量下回归系数各分量以及响应变量各分量随迭代进行的变化趋势（图2）。
算例2 仍然以指数型基准生存函数为例，利用随机截断产生的随机样本，对于非参数的模型估计基准生存函数和g(x)。
（三）模拟条件 基准生存函数为均值150的指数分布，
，我们将x的取值范围设为[2π，4π]，在此区间上等间隔分别取60，120，300个点，作为协变量水平。在每个协变量水平下，各取20，80个样本，分别进行估计。当两步迭代中向量V（如前定义）的差的范数小于
时，迭代终止。
（四）定义
图9～图14为基准生存函数在各情况下模拟50次的拟合情况。
（五）结果分析
在算例1中以MSE和回归系数均值为标准来比较，可以看到虚拟响应算法和偏似然方法的结果差距并不大，而且随着样本量的增大结果相差也越来越小。考虑到Cox偏似然方法的良好理论性质，我们可以认为虚拟响应算法也有这些良好的性质。而从实际应用的角度来讲，虚拟响应算法的直观性和求解的简单性优势更为重要。从图形上看起来随着样本量的逐渐增大，很明显的我们的算法最终结果越来越趋向于真值。因此可以预期我们的算法有相合性。
在算例2中计算结果对于协变量的组数有比较高的要求，如果协变量水平数较少，对于g(x)的估计有时候偏差较大，极端情形下会非常大，导致出现基准生存函数的估计结果近似恒等于1的不合理情况。分析的结论是当协变量水平数较少，则在作局部多项式估计不能保证有足够的点落在每个邻域内，从而导致非参估计的结果较差所致。至于每组协变量下的样本个数，模拟中可以看到，只要达到一定数量即可，提高样本个数对于结果的影响并不如增加协变量水平数来的明显。
上述两个算例在尾部的拟合不是很好，主要因为计算过程中为了避免ln0导致的数值上的麻烦，我们在定义乘积限估计的时候通过控制最后一个观测时刻的截断指示变量使得估计结果永远不会下降到0。模拟时在尾部没有产生随机样本，而我们最终的估计是一个阶梯跳跃函数，因此在此区间内估计结果不再跳跃。此问题可以通过加大每个协变量水平下样本量来解决。
此外，在以上两个模拟中，我们发现无论是虚拟响应变量还是回归系数，在具体一次模拟过程中都有单调性，正如我们图形中显示的那样。这也是算法收敛性的一个很好的佐证和线索。当然由于我们在拟合虚拟响应变量过程中所用的具体方法根据不同的模型而不尽相同，目前还没有统一的算法收敛性证明。这一点工作留待日后进一步研究。
三、问题与展望
总的来说，对于每个协变量下有多个样本的模型，我们的处理方法取得了较好的数值结果。但是对于每个协变量水平下只有一个样本的情况，我们的算法则不能处理，每个协变量下样本量较少的情况，可以通过加大观测的协变量水平数来达到很好的结果。
当采用最小二乘迭代算法的时候，依照最小二乘的理论要求协变量之间不能有很强的复共线性，否则数值结果不理想。由于试验条件往往是人为设定，这一要求较容易实现。如果不能满足上述要求，我们可以参考文献[12]中的内容，采用岭回归方法，取得能够获得较小MSE的结果，当然这一结果往往牺牲了无偏性。
当协变量作用形式为非参数的情形，我们也可能遇到协变量为多维情形，在文献[3]中提到了一些降维的技术，将问题转化为一维的情况。
我们的算法收敛性以及结果的大样本性质目前还没有得到理论证明。但是从数值模拟的情况来看，我们有足够的理由相信在一定条件下两者的都有较好的结果。这一工作留待继续研究。
本文最基本的思想是根据模型的特点构造虚拟响应，然后通过对拟合的得到协变量作用形式的估计。考虑到Cox模型的广泛应用，本文是基于Cox模型做的研究。事实上这一思想的应用可以不仅仅局限于Cox模型，只要模型的特点能够构造合适的虚拟响应，我们都可以采用此思路。一个潜在的例子是加速寿命试验，在此模型中不同的外加应力也会影响生存函数的形式，并直接的反应在寿命上。此时基准寿命和受外加应力影响的情况的下的寿命存在一个比例关系，这个比值是协变量的函数，我们可以以此作为虚拟响应。我们可以利用统计方法得到各情况下平均寿命的估计，并构造比值，然后采用类似Cox模型中的算法进行迭代。
标签：回归系数论文; 协变量论文; 似然函数论文; 样本量论文; 虚拟变量论文; 生存分析论文; 迭代模型论文; 迭代计算论文; 算法论文;