竞争风险下纵列持续数据随机效应模型的估计与模拟研究，本文主要内容关键词为：效应论文,模型论文,风险论文,竞争论文,数据论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
1 研究背景
持续数据分析是计量经济学近年来发展非常迅速的一个新研究领域。某事件从开始到结束或终止所经历的时间长度T称持续期，如从失业到再就业、从发放抵押贷款到违约，等等。持续数据分析在生物医学领域也称生存分析，着重研究持续时间的概率分布或者事件终止的概率（也称危险率）。如果分析以一组协变量(x)为条件，条件危险率的典型模型是Cox[1]比例危险模型：
其中，
给出了个体i在单位时间内结束当前状态的概率，协变量(covariate)
是解释个体i结束当前状态原因的外生变量；
(t)是基本危险率。比例危险的含义是对具有不同协变量和
的个体，危险率的比率是一个常数。
实际问题中，引起事件终止的原因通常不止一个。比如失业者找到工作或者退出劳动力市场都会导致失业状态结束，引起贷款抵押终止的原因可能是贷款者提前还贷(prepayment)也可能是违约(defaults)。此时对事件终止概率的研究需要考虑不同的原因，也即不同的风险，这就是竞争风险(competing risks)模型。Han和Hausman[2]对竞争风险模型进行了开拓性研究，目前这一模型已成为处理具有多重原因的经济持续数据的重要方法。竞争风险下，研究者只能观测到个体最小的事件发生时间T，通常看作连续的随机变量；导致事件发生的原因用C表示(C=1,…,L)，C是离散随机变量。竞争风险下的持续期研究需要识别C和T的联合概率分布，通过设定部分分布函数(Sub-Distribution Function)来解决这一问题。假设第j个风险下的部分分布函数为
其中f(j,t)是部分密度函数(Sub-Density Function)。
持续数据建模往往遇到不可观测的异质性(unobserved heterogeneity)问题，遗漏重要解释变量或变量存在测量误差，都可能导致个体之间的异质性。尤其是纵列持续数据的异质性，缘于对同一个体进行重复观测或者将数据按照不同来源分组所造成的数据之间的相关性，导致组内事件更容易（或更不容易）发生。这种异质性是内在的，不能通过模型设定予以消除，很多文献中将这种不可观测的异质性称作“脆弱性”(failty)。忽略了个体的异质性将导致估计的参数有偏和不一致，如低估协变量的效应、高估相对危险率（Keyfitz和Littman[3]、Hougarrd[4]与Aalen[5]、Lancaster[6]）。Pickles和Crouchley[7]对脆弱性问题的综述表明，忽略异质性的影响会使协变量系数的估计趋近于0。
单个风险下，为了描述个体不可忽略的异质性，往往在危险率函数上乘一个随机效应因子，也即将Cox比例危险模型扩展为Cox脆弱性比例危险模型。由于随机效应的分布以及基本危险率均未知，随机效应模型的估计有两种思路，第一种思路是对基本危险率或随机效应分布函数施加参数结构假设，如Lan-caster[8]等，一般采用边际极大似然估计，计算中需采用EM算法。但EM算法的收敛对初始值的选择和停止计算的规则非常敏感，而且，其中的E步期望形式涉及贝叶斯后验分布，除了少数几种特定的分布，一般不易写出解析形式的后验期望。出于数学上处理的方便，Klein[9]建立了伽马脆弱性模型的EM算法。Hougarrd[10]提出三种脆弱性因子的分布：伽马(Gamma)分布、逆高斯(Inverse Gaussian)分布和正稳定(Positive Stable)分布。贝叶斯抽样技术以及MCMC方法在这一类随机效应模型的估计中也有广泛的应用，虽然这些方法较为稳健，但计算相当复杂。第二种思路是应用混合广义线性模型的理论。Lee和Nelder[11]提出了估计混合广义线性模型的等级似然估计(Hierarchical Likelihood)，Lee等[12]进一步完善该理论，扩展了其应用范围。Aitkin和Clayton[3]、Ha和Lee[14]、Ma等[15]以及Ha和Lee[16]的研究表明，Cox比例危险模型以及Cox脆弱性比例危险模型本质上可以纳入广义线性模型和混合广义线性模型框架；Ha等[17]将等级似然估计方法用于单风险的Cox随机效应模型，随后Ha等[18]又应用于存在删失数据的混合线型模型，Ha和Lee[19]将其应用于多重混合线性模型(Mulitilevel Mixed Linear Model)。
竞争风险下考虑随机效应，情形更为复杂，目前还没有一个规范的框架和可操作的通用程序。Vaupel et al.[20]考虑对整体危险率加入脆弱性因子；假定在给定随机因子V的条件下，联合生存函数为：
其中
为累积部分危险率，无条件的联合生存函数通过积分得到：
这种积分往往形式复杂难于计算，更重要的是它将随机效应乘在整体危险率上，不符合实际情况。Malay等[21]对危险率函数和不同风险的选择分别建模，假定不同风险下的随机效应来自荻立克莱(Dirichlet)过程，推导了ECM估计过程；Robert等[22]对重复测量数据和竞争风险建立联合模型，假定正态随机效应，推导了估计的EM算法。这一类估计方法沿袭了单风险下随机效应模型的第一种思路，但计算难度更大，发展缓慢。
本文针对包含右删失的纵列持续数据，对每种风险建立一个随机效应Cox比例危险模型，将等级似然估计的运用由单风险扩展到多风险，解决了竞争风险下纵列持续数据随机效应模型的估计问题。本文后面的内容安排如下：第二部分建立竞争风险下纵列持续数据的随机效应模型，论证它属于混合广义线性模型范畴，可以采用等级似然方法进行估计；第三部分提出该模型等级似然估计的步骤；第四部分对所建立的模型进行估计和模拟研究；第五部分对全文进行总结。
2 模型的建立
在广义混合线型模型中，如果随机效应和协变量在线性预测部分中可加，那么这种随机效应尺度就是弱规范尺度，弱规范尺度的随机效应总可以找到；而且，此时协变量系数β的估计和随机效应v的实现可以通过最大化等级似然获得，等级似然函数二阶导数的相反数之逆仍然给出了协变量系数和随机效应实现估计值的方差。上面的这些性质为等级似然的统计推断提供了可行的理论基础，关于这些性质的进一步讨论可以参阅Lee等[12]。
在上述分析的基础上，我们定义竞争风险下比例危险模型(1)和(2)的等级似然函数：
3 模型估计
上面指出，竞争风险下的脆弱性比例危险模型能够纳入混合广义线性模型的框架，下面具体讨论如何运用等级似然方法来估计模型参数。
3.1 随机效应的参数化处理和基本危险率的估计
等级似然函数中所包含的随机效应事实上不可观测，但它可看作关于β和v联合似然函数，其中v可看作不可观测的随机效应的实现，因而可将v处理成固定参数。假定数据的生成过程为：随机效应v来自某一参数分布族，给定v的实现值和协变量x，感兴趣事件的持续期的生成过程符合上述设定的模型(1)、(2)，此时，对于观测到的持续时间，随机效应的实现值就是一组固定的常数。将
一起作为未知参数，最大化式(5)，这就避免了高维积分的计算，不仅比EM方法简便，而且能够同时得到随机效应的最优无偏预测和参数估计的渐近方差阵[11]。我们采用Newton-Rapson迭代方法最大化等级似然函数式(5)，大量的模拟研究表明，对于混合广义线性模型以及脆弱性模型的等级似然函数求解问题，Newton-Rapson算法能够很快收敛（见Ha和Lee[16]的讨论），本文后面的模拟研究也证实了这一点。
模型中累积基本危险率和基本危险率的估计，采用类似于Breslow[23]的方法，有
实际中，我们观察到的经济数据常常是区间删失的，比如失业持续时间是以周为单位的观测，知道的确切时间只能是某周之内；抵押贷款的持续时间是以月为单位的观测，抵押终止的发生是在某个月之内。区间删失数据同时考虑竞争风险和随机效应时更为复杂，虽然从理论上说可通过积分解决，但实际上往往难于计算。近年来经济学的实证研究中对于竞争风险下的区间观测数据均采用McCall[25]的设定，将区间观测数据的竞争危险模型的基本危险率设为分段线性函数，包括Deng等[26]、Ciochetti[27]等。Mark Yuying An[28]证明，对于竞争风险下区间观测数据的比例危险模型，非参数的基本危险率函数是不可识别的，为了得到可以求解的似然函数，需对基本危险率函数做参数设定。由于本文考虑更加灵活的半参数模型(1)和(2)，所以不需要基本危险率的参数形式。对于区间观测，我们考虑以区间的中点作为实际观测的近似。我们认为，当观测的区间个数比较多时，这种近似是合理的，模拟研究中同时给出了这种区间近似的估计结果。这模拟结果列在表1至表5，需要说明的是，表1、表2和表5模拟次数为200次，模拟估计过程中，由于涉及大量矩阵运算，程序运行时间相当长，对模拟结果的试分析发现模拟100次的结果和模拟200次的结果差异非常小，因此出于程序运行时间的考虑，表3和表4列出了模拟160次的结果。
根据上述设定，我们分别按照不同的组（个体）数目以及不同的组内（重复）观测，对随机效应下协变量参数以及随机效应方差进行了模拟估计，数据的存储和模拟估计过程采用SAS(IML)语言编程实现。为估计方便，对随机效应的方差采用了再参数化的方法，估计的是方差的倒数。模拟结果列于表1至表4。各表中列出了参数模拟估计的平均值、标准差，以及平均标准误（标准误由等级似然的二阶导数矩阵之逆的右下角对角元给出）；最后两列给出了参数模拟估计值95%的置信区间包含参数真值的比例。从表中看出，协变量的系数估计以及随机效应方差的估计在各种情形下都相当精确，协变量参数模拟估计结果的标准差和平均标准误非常接近。这说明对竞争风险下纵列持续数据的模型设定是合适的，同时也说明，等级似然函数的二阶导数矩阵的逆作为参数估计的方差阵的估计是正确的。
表1 组数目为10和组内个体数目为10的200次模拟估计结果
注：表中h_lik指等级似然估计的结果，Cox指标准Cox模型（不考虑随机效应）的估计结果；表中“区间近似”是文中提到的区间观测模拟的数据；表中b1指
的估计，b2指
的估计，p指两个随机效应相关系数的估计，sd1和sd2分别指两个随机效应标准差的倒数的估计。下面的表中符号相同。
表2 组数目为20和组内个体数目为10的200次模拟估计结果
表3 组数目为10和组内个体数目为20的160次模拟估计结果
注：由于这一程序运行的时间太长，表3和下面的表4各进行了160次模拟。
表4 组数目为10和组内个体数目为30的160次模拟估计结果
表5 组数目为20和组内个体数目为20的200次模拟估计结果
注：表5中列出了两个随机效应的相关系数以及它们的标准差的倒数的估计。
表1至表4同时给出了利用区间数据模拟估计的结果。对于区间观测的模拟，当协变量系数真值设为-1，-2时，模拟的事件持续时间主要集中在0～10之间，大于6的数据不多，个别数据在10以外，因此假设区间长度为0.05，将观测数据划分为120个区间，大于6的数据认为是删失的，这样真正的删失率不仅仅是0.1。可以看出，在各种情形下，区间观测数据的估计和真实观测数据的估计结果非常接近。这说明，当观测的区间比较密集时，将其当作连续数据是合理的，台湾的林左裕[29]的研究采用了这一形式。这一处理方式很有实际意义，经济持续数据观测的区间往往是等长的，也较为密集，如抵押贷款的持续数据，贷款期限可能是10年以上，观测是月度数据，这样，观测的区间有100个以上。
对比表1至表4的结果可以看出：①表1到表2显示组数目增加而带来的样本容量的增加并没有显著改善协变量的系数估计。究其原因，我们认为这是由于组数目的增加同时也带来了待估参数的增加，比如从表1到表2，需估计的随机效应实现由20个变为40个。②由表1、表2、表3的结果看出，随着组内个体数目的增加，协变量系数的估计变得更加精确，包括区间观测数据的估计也是如此。③若不考虑随机效应的影响，仅用Cox比例危险模型进行估计，协变量系数的估计的绝对值严重下偏。这是显然的，因为设定的数据生成过程带有随机效应，而Cox比例危险模型未能反映随机效应。这说明，如果我们实际获得数据的生成过程有随机效应，估计结果将会有偏，而且偏误不会随着样本容量的增加而消失。
应该指出的是，模拟中我们发现，两个随机效应的相关系数的估计不如协变量系数和随机效应标准差的估计准确，偏差较大。我们设定不同的删失率（0.1和0.3）和相关系数（0.1和0.3），按不同组合作了进一步的模拟，仍然设协变量系数分别为-1，-2，数据生成过程和上面的模拟过程一样，设两个随机效应的标准差均为1.5，个体（组）数目设为20，组内观测设为20，每种情形均模拟200次。模拟结果列于表5，可以看出，随机效应相关系数的估计不够稳定，标准差的估计虽较稳定但有偏小的趋势。之所以出现这样的问题，我们认为，原因是对随机效应分布参数的估计采用边际似然的一阶近似，如果能够采用更高阶的近似，应该能够得到较好的结果，但是高阶近似会使得模型估计的复杂度大大增加。考虑到我们主要关注的是存在随机效应的情况下协变量参数的估计，其精确度在本文的估计中是有保证的，而估计方法相对于目前文献上的估计方法要简便得多，所以结果是可以接受的。
5 结论
持续数据分析中经常遇到不可观测的异质性问题，尤其竞争风险下的纵列数据带来的异质性问题比较复杂，目前文献上还没有一个通用的分析框架和估计程序。本文针对这一问题展开研究，并考虑数据右删失情形下，论证竞争风险下的脆弱性比例危险模型属于混合广义线性模型的范畴，推导出用于模型估计的等级似然函数，给出了类似于Breslow[23]的基本危险率估计。采用等级似然方法，我们将随机效应的取值看作随机效应的实现值，当作未知参数和模型协变量参数一起进行估计以避免高维积分的计算，而且给出了随机效应参数的最佳无偏预测。通过模拟研究我们发现，协变量系数的估计是相当精确的，能够有效克服异质性产生的偏差。这说明，竞争风险下持续数据建模忽略异质性（随机效应）是相当危险的，这和文献中单风险下随机效应模型的研究结论一致。
标签：随机效应模型论文; 大数据论文; 似然函数论文; 协变量论文; 参数估计论文; 风险系数论文; 风险模型论文; 异质性论文;