（概率论与数理统计专业论文）广义加速失效时间frailty模型.pdf

m a s t e rd i s s e r t a t i o no f 抬a r2 01 0 u n i v e r s i t yi d ：1 0 2 6 9 s t u d e n ti d ：510 7 0 6 0 5 0 0 4 g e n e r a l i z e da c c e l e r a t e df a i l u r et i m e d e p a r t m e n t m a j o r f r a i l t ym o d e l s c h o o lo ff i n a n c ea n ds t a t i s t i c s p r o b a b i l i t ya n dm a t h e m a t i c a ls t a t i s t i c s r e s e a r c hdi r e c t i o nsu r v i v a la n a l y s i s s u p e r v i s o r a u t h o r d a t e 删黼 华东师范大学学位论文原创性声明 郑重声明：本人呈交的学位论文广义加速失效时间f r a i l t y 模型，是在华 东师范大学攻读石页蚜博士( 请勾选) 学位期间，在导师的指导下进行的研究工作 v 及取得的研究成果。除文中已经注明引用的内容外，本论文不包含其他个人已经发 表或撰写过的研究成果。对本文的研究做出重要贡献的个人和集体，均已在文中作 了明确说明并表示谢意。 作者签名巧懈日期：蜱 华东师范大学学位论文著作权使用声明 广义加速失效时间f r a i l t y 模型系本人在华东师范大学攻读学位期间在导 师指导下完成的砸切博士( 请勾选) 学位论文，本论文的研究成果归华东师范大 学所有。本人同意华东师范大学根据相关规定保留和使用此学位论文，并向主管部 门和相关机构如国家图书馆、中信所和“知网”送交学位论文的印刷版和电子版； 允许学位论文进入华东师范大学图书馆及数据库被查阅、借阅：同意学校将学位论 文加入全国博士、硕士学位论文共建单位数据库进行检索，将学位论文的标题和摘 要汇编出版，采用影印、缩印或者其它方式合理复制学位论文。 本学位论文属于( 请勾选) ( ) 1 经华东师范大学相关部门审查核定的“内部”或“涉密学位论文木， 于年月日解密，解密后适用上述授权。 ( 乃2 不保密，适用上述授权。 作者签名，譬车触 日期：趔兰么12 导师签名： ，“涉密”学位论文应是已经华东师范大学学位评定委员会办公室或保密委员会审定过的学位论文( 需 附获批的华东师范大学研究生申请学位论文“涉密”审批表方为有效) ，未经上述部门审定的学位 论文均为公开学位论文。此声明栏不填写的，默认为公开学位论文，均适用上述授权) 。 陈鹏程硕士学位论文答辩委员会成员名单 姓名职称单位备注 汤银才教授华东师范大学主席 曾林蕊副教授华东师范大学 程依明副教授华东师范大学 目录 摘要 a b s t r a c t ( 英文摘要) 目录 第一章引言 1 1 问题的提出 1 2 研究现状 1 2 1 基本概念 1 2 2 比例风险模型 1 2 3 加速失效时间模型 1 2 4 比例风险f r a i l t y 模型 51 2 5 加速失效时间f r a i l t y 模型 1 3 论文安排 第二章带有广义伽马分布f r a i l t y 项的加速失效时间模型 2 1 模型分析 2 2 估计过程 2 3 数值模拟 2 4 渐进推导 第三章广义加速失效时间打a i l t y 模型的m l e m 算法 3 1 模型分析 3 2 估计过程 3 3 数值模拟 第四章实例分析 结论 参考文献 致谢 3 7 3 9 4 1 4 4 v 1 4 4 5 6 7 8 8 加 m n 鹅 勰 嚣驼 1 1 心绞痛观测时间的生存函数的k m 估计 2 2 1 偏似然函数的核光滑处理效果示意图1 5 表格目录 表格目录 1 1 心绞痛观测时间数据 3 2 1 模拟产生5 0 0 个数据集，每个数据集的病人分为2 0 0 组，每组2 名病人，在 删失率为7 0 的状态下考察模型参数的估计1 9 2 2 模拟产生5 0 0 个数据集，每个数据集的病人分为2 0 0 组，每组2 名病人，在 删失率为8 5 的状态下考察模型参数的估计2 0 2 3 模拟产生5 0 0 个数据集，每个数据集的病人分为2 0 0 组，每组4 名病人，在 删失率为7 0 的状态下考察模型参数的估计2 1 2 4 模拟产生5 0 0 个数据集，每个数据集的病人分为2 0 0 组，每组4 名病人，在 删失率为8 5 的状态下考察模型参数的估计2 2 2 - 5 模拟产生5 0 0 个数据集，每个数据集的病人分为4 0 0 组，每组2 名病人，在 删失率为7 0 的状态下考察模型参数的估计2 3 2 6 模拟产生5 0 0 个数据集，每个数据集的病人分为4 0 0 组，每组2 名病人，在 删失率为8 5 的状态下考察模型参数的估计2 4 譬1 模拟产生5 0 0 个数据集，每个数据集的病人分为2 0 0 组，每组2 名病人时考 察模型参数的估计3 4 3 - 2 模拟产生5 0 0 个数据集，每个数据集的病人分为2 0 0 组，每组4 名病人时考 察模型参数的估计3 5 3 - 3 模拟产生5 0 0 个数据集，每个数据集的病人分为4 0 0 组，每组2 名病人时考 察模型参数的估计3 6 4 1 心绞痛临床试验参数估计值3 7 中文摘要 摘要 生存分析是一种处理生存数据的统计方法。传统的生存分析模型都是建立在生存 数据相互独立的基础上的。但是在实际应用过程中，生存数据不独立的情况处处可 见。经典的比例风险模型和加速失效时间模型都无法对此类数据进行合理有效的解 释。因此我们需要一种新的模型来处理这种相关的生存数据。f r a i l t y 模型是分析此类 数据的经典模型。我们在模型中用f r a i l t y 项来评估数据之间的相关性。然而到目前为 止，大多数研究都是集中在比例风险f r a i l t y 模型上面。加速失效时 f i f r a i l t y 模型由于其 计算过于复杂，并未引起人们的太多关注。但是随着统计的发展，越来越多的人发现 了加速失效时间f r a i l t y 模型的优越性。在这篇论文中，我们结合广义加速失效时间模 型，首先对加速失效时间f r a i l t y 模型的形式进行一些修改。让加速失效时间f r a i l t y 模 型与广义加速失效时间模型的形式更为接近。在此模型的假设下，我们给出两种不同 的模型估计方法。在第一种方法中，我们基于普通的e m 算法，把模型中f r a i l t y 项的 分布拓展到了广义伽马分布。这让我们的估计方法更具灵活性和实用性。随后本文引 入了另一种基于m l e m 算法的模型估计方法。在这种方法中，我们把极大似然方法 和e m 算法进行了有效的融合，从而大大提高了我们的模型估计速度。最后，本文还 用组真实的临床数据检验了我们的算法。 关键词：生存分析，加速失效时间模型，广义伽马分布，e m 算法，m 己一e m 算法 a b s t r a c t s u r v i v a la n a l y s i si sak i n do fs t a t i s t i c a lm e t h o df o rs u r v i v a ld a t a t h et r a d i t i o n a l s u r v i v a lm o d e li sb a s e do nt h ei n d e p e n d e n ts u r v i v a lt i m e b u tt h ed e p e n d e n tt i m ei sv e r y c o m m o ni np r a c t i c e t h o s ec l a s s i c a ls u r v i v a lm o d e l ss u c ha sp r o p o r t i o nh a z a r dm o d e lo r a c c e l e r a t e df a i l u r et i m em o d e lc o u l dn o tb eu s e dt of i tt h ed a t a s ow en e e dd e v e l o pa n e wm o d e lf o rt h i ss i t u a t i o n t h ef r a i l t ym o d e li so n eo ft h em o s tp o p u l a rm o d e l su s e d t oa n a l y z ec l u s t e r e df a i l u r et i m ed a t a a n dt h ef r a i l t yt e r mi nt h em o d e li su s e dt oa s s e s s a s s o c i a t i o n si ne a c hc l u s t e r t h ef r a i l t ym o d e lb a s e do nt h es e m i p a r a m e t r i ca c c e l e r a t e d f a i l u r et i m em o d e la t t r a c t sl e s sa t t e n t i o nt h a nt h eo n eb a s e do nt h ep r o p o r t i o n a lh a z a r d s m o d e ld u et oi t sc o m p u t a t i o n a ld i f f i c u l t i e s i nt h i sp a p e r ，w ef i r s tm o d i f i e dt h em o d e lt o m a k es u r eo u rm o d e li sm o r ec l o s e dt ot h eg e n e r a l i z e da c c e l e r a t e df a i l u r et i m em o d e lt h a n t h eo t h e r s u n d e rt h i sa s s u m p t i o n ，w ed e v e l o pt w od i f f e r e n tn e w e s t i m a t i o nm e t h o d sf o r t h es e m i p a r a m e t r i ca c c e l e r a t e df a i l u r et i m ef r a i l t ym o d e lb a s e do nt h ee m l i k ea l g o r i t h m a n dt h em l e ml i k ee s t i m a t i o nm e t h o d t h ef i r s tp r o p o s e dm e t h o di sc o m p a r e dw i t h t h ee x i s t i n ge ma l g o r i t h m ，w h i c hi n c o r p o r a t e st h ek e r n e ls m o o t h i n gi nt h em - s t e p w e e x t e n dt h ef r a i l t yd i s t r i b u t i o nt ot h eg e n e r a l i z e dg a m m ad i s t r i b u t i o nw h i c hm a k e so u r m e t h o di sm o r ef l e x i b l e t h e nw ed e v e l o pan e wm e t h o db a s e do nt h em l e ma l g o r i t h m w ei n c o r p o r a t et h em a x i m u ml i k e l i h o o dm e t h o da n dt h en o r m a le ma l g o r i t h mi nt h i s m e t h o d w ei m p r o v eo u re s t i m a t i o ne f f i c i e n c yt h r o u g hi t f o ri l l u s t r a t i o n ，w ea p p l yt h e p r o p o s e dm e t h o dt ot h ed a t as e to fs u b l i n g u a ln i t r o g l y c e r i na n d o r a li s o s o r b i d ed i n i t r a t e o na n g i n ap e c t o r i so fc o r o n a r yh e a r td i s e a s ep a t i e n t s k e yw o r d s ： s u r v i v a la n a l y s i s ，a c c e l e r a t e df a i l u r et i m ef r a i l t ym o d e l ，g e n e r a l i z e d g a m m ad i s t r i b u t i o n ，e m ，m l e m 一v 一 第一章引言 i 一- 1 1 - 立 弟一覃 芦i 已 ji 苗 生存分析是近二三十年发展起来的数理统计新分支。它是针对生存时间提出的一 种统计分析方法，主要对随机删失型数据进行研究。随机删失数据是生命科学、生物 制药、可靠性实验以及其他一些实际问题中常常碰到的一种重要统计数据。其主要特 征就是数据观测不完全，我们只能大概了解事件发生在某段时间，而无法得知具体的 时间点。这种特殊的数据要求有新的统计方法来对其进行分析。生存分析就是在这种 背景下被提出并发展的。生存分析理论及方法在生命科学、可靠性工程、保险精算、 社会学、市场学等各个领域都有广泛的应用。 1 1 问题的提出 生存分析领域最重要的一类分析方法就是时间到事件的分析方法。而据此提出的 模型也被称之为生存模型。经典的生存分析模型是比例风险模型( p r o p o r t i o n a lh a z a r d s m o d e l ) 和加速失效时间模型( a c c e l e r a t e df a i l u r et i m em o d e l ) 。当生存数据相互独立的时 候，这两种模型可以有效的解释和拟合数据。然而，在现实环境下，生存数据之间并 不一定存在独立关系。这就驱使我们寻找这些经典模型的替代品，以便可以合理的解 释相关的生存数据。 在临床试验的过程中，我们经常可以遇到生存时间非独立的情况。而这些数 据又可以被细分为相关失效时间数据( c o r r e l a t e df a i l u r et i m ed a t a ) $ 口群组失效时间数 据( c l u s t e r e df a i l u r et i m ed a t a ) 。这些数据非独立的原因往往是由于受测个体来自相同的 受测环境。进步细化这种相关性，当测试数据来自同一个病人的时候，例如测量癫 痫病发病周期的数据，此时的数据被称为相关失效时间数据；当我们的数据来自于同 一个医院的病人的时候，由于这家医院的设备和医生的共享性，导致我们的数据在同 一群组内产生相关性，此时的数据就是群组失效时间数据。 我们来看组真实的临床数据。下面的这组数据来自于一个冠心病的临床试 验( d a n a h ye ta 1 ，1 9 7 7 ) 。在这个临床试验中，我们共招募t 2 1 名冠心病患者，采用 交叉设计的试验方法，对口服二硝酸异山梨醇的疗效进行评估。在这个临床试验 中，病人要参加两组不同的治疗。起初，病人被先后要求舌下压硝化甘油( s l t ) 和 舌下压安慰剂( s l p ) 进行脚踏车运动。我们记录病人的心绞痛复发时间。在随后的 试验中病人将先后口服二硝酸异山梨醇( 0 i ) 和口服安慰剂( o p ) 。每服用一种药物之 前( o i oa n do p o ) ，服用药物之后的l 小时( 0 i la n do p l ) ，3 d , 时( 0 i 3a n do p 3 ) 署u 5 d , 时( 0 i 5a n do p 5 ) ，病人将进行脚踏车运动，运动会一直持续到病人的心绞痛复发。 一1 一 1 1 问题的提出 受测个体的心绞痛复发时间将作为我们的观测数据。需要注意的是，在个别情况下， 由于病人的体力原因，我们的观测数据产生了删失，这种情况在此试验中一共出现 了1 4 例。观测数据从1 1 0 秒至i j 7 4 3 秒，这表明病人在两次不同的测试间隔内获得了足够的 休息时间，从而保证了数据的真实有效性。 我们令聪为第i 组的第j 个病人真实的心绞痛复发时间，j = 1 ，1 0a n di = 1 ，2 1 。从而我们的观测值t i j = m i n ( t i ；，c i j ) ，其中g j 被称作随机删失变量， 它和受测个体是一一对应的。而当正j = 砭时，删失指标以f = 1 ，相对应的， 当t i j = g f 时，删失指标瓯f = 0 ，此时我们的数据产生了删失。需要明确的是，我们 在此假设删失时间并不包含参数相关的信息。我们用x ；表示协变量，在此临床试验中 它被用来表示受测药物在固定时间点的疗效。为了更加简单的表示数据，我们用0 i ，来 表示我们最终得到的数据( t 讲瓯f ，x t j ) 。我们在表1 - 1 中列出了我们的观测数据，表中 的木号表明是删失数据。 下面的图1 1 用k a p l a n m e i e r 生存曲线来估计相应药物在固定时间点的医疗效 果。我们可以看出，o i l 有着较高的生存概率，而相对应的o p o 的生存概率比较低。 o 誊啦 豸 。 巳 4 嘉6 , q - o o t i m et oa n a g i n ap e c t o r i s 图1 - 1 心绞痛观测时间的生存函数的k m 估计 一2 一 第一一章引言 表1 - 1 心绞痛观测时间数据 s l p s l to p 0o p l0 p 3 0 p 50 i o0 1 1 0 1 30 1 5 3 一 1 2 研究现状 由于我们是从同一病人身上反复测量得到的数据，这就导致了观测数据之间不是 相互独立的。此时，若我们依旧使用独立的生存模型对数据进行拟合，则我们对药物 疗效的估计就会产生一定的偏差。在某些特定的情况下，这种偏差是致命的。因此为 了更准确的评估口服二硝酸异山梨醇疗效，我们必须在生存模型中考虑数据的相关 性。进而对药物疗效进行准确的评估。 1 2 研究现状 本部分将着重介绍一些生存分析的基本概念和经典的生存模型。 1 2 1 基本概念 生存模型采取的是从时间到事件的分析方法。因此，事件的发生时间也就成了最 重要的观测变量。生存分析领域的绝大多数函数都是以时间为随机变量进行定义的。 在某些情况下，我们无法观测到事件发生的具体时间，此时我们就称数据发生了删 失。具体在本文中，我们主要研究右删失数据。这种数据是由于事件发生的时间超过 了我们记录的时间节点。而对删失数据的处理也为生存分析的发展提出了挑战。下面 我们就介绍一些生存分析领域的基本函数： 生存函数： s ( t ) = p ( t t ) = p ( t t ) ( 1 一1 ) 密度函数： 风险函数： ，( ) = l i t m 一。1 j 尸( t t + ) = 出m 0 壶 s ( ) 一s ( 汁) d s ( t ) = 一一 d t ， 1 州) _ l i m o 面 一 1 = l i r ai 一 0 t f ( t ) s ( t ) p ( t t t + a t l t t ) p ( t t t + a t ) p ( t t ) ，4 ( 1 - 2 ) ( 1 - 3 ) 第一章引言 累计风险函数： 右删失数据的似然函数： 何( ) = o 。庇( z ) d x = - - 1 0 9 s ( ) ( 1 4 ) l = f ( t t ) 正s ( t t ) 1 也= h ( t i ) 正s ( t ) ( 1 5 ) i=1扛：1 1 2 2 比例风险模型 d r c o x ( 1 9 7 2 ) 在1 9 7 2 年提出了如下模型 h ( t x ) = h o ( t ) e x p ( 3 x ) ( 1 - 6 ) 其中h ( t l x ) 表示在协变量为x 条件下的风险函数，而相对应的h o ( t ) 被称为基准风险函 数。当x = 0 的时候，风险函数h ( t l x ) 和基准风险函数h o ( ) 相等。 我们现在考察一下在两个不同协变量的条件下，风险函数之间的关系。我们 用x 1 来表示第一个协变量，用x 2 来表示第二个协变量。那么两个风险函数的比可以 写成： 坐盟：一一尥)(1-7)h(_ h l x 2 ) 2 e x p ( ( x l 一五2 j ) 我们通过上式可以看出，两个不同的风险函数之间的比值是一个由x l 和恐决定的常 数。所以，我们把此模型称为比例风险模型。而e x p ( 卢( x 1 一x 2 ) ) 被称为x 1 和磁的相 对风险。令h ( t i x ) 和i - l o ( ) 为h ( t l x ) 和h o ( ) 相对应的累计风险函数。则累计风险函数 之间满足下式： l o g ( h ( t t i x ) ) = p x + l o g ( h o ( t ) ) ( 1 - 8 ) 因此，两个累计风险函数的对数值曲线应该是一种平行关系，这也被称为比例风险模 型的模型假设。通常我们在使用比例风险模型之前就必须检验此假设。 比例风险模型最吸引人的地方在于，对参数p 的估计并不依赖于基准风险函数的 取值。而这种估计将通过部分似然的方法来实现。下面我们将简要介绍针对比例风险 模型的部分似然估计方法。在此之前我们先对一些符号进行定义。遵照前文，我们 用( t i ，瓯，z t ) ，i = l ，n 来表示已知数据，其中的n 表示样本总量。除此之外我们还需 要假设，删失时间不包含任何与参数相关的信息，以及，在协变量x 给定后，事件发 生的时间和删失时间之间相互独立。 5 一 1 2 研究现状 在满足以上前提条件下，为了简化似然函数的表达式，我们假设数据没有结。 令t a t 2 0 ) 是形状参数，a ( 0 ) 是尺度参数，且一。 q 跳= 喜喜w 射训嘞以矿一el ( z a q 、础扩 根据全似然函数，我们可以推导出z 的后验分布： 泡，。) n , 5 - a - - - 2 - 1 唧( - t z q 姜础扩觑卜g t ) ( 2 - 5 ) ( 2 - 6 ) 上式中的p ，厅，方是最新一次估计的参数值，而0 t j 是我们得到的真实观测值。我们希望 可以求出z 的后验分布的期望表达式。但是在实际操作的过程中，由于期望的表达式 众多，且由于此后验分布不能被化为经典分布族，所以我们不可以用简单的参数组合 来表示z 的一系列期望。因此在估计过程中，我们考虑利用m c m c 方法。在估计过 程中，我们基于z 的后验分布对z 进行抽样，然后用它们的矩估计来代替似然函数里 面的期望部分。从而，我们在e - s t e p 消除了潜在变量的不利影响，使我们的似然函数 变成了可估的形式。 一1 2 一 第二章带有广义伽马分布厂，a i t q 项的加速失效时间蟮型 m s t e p 在消除完潜在变量对似然函数的影响之后，我们要在m - s t e p 对似然函数进行极 大化，从而求出p ，盯，q p 2 极大似然估计值。然而由于似然函数e f 2 包含另外的未知 量h o 他j e 一卢锄1 和风( t i i e 一眠，) ，我们依旧无法对似然函数进行任何处理。所以，我们 此时基于z e n g 年a l i n ( 2 0 0 7 ) 里面针对普通加速失效时间模型的分析方法，对我们的似然 函数进行核光滑处理，把我们的似然函数由半参数形式转化成参数形式。在下面的叙 述中，我们用e ( w t ) 来代替e ( z 7 q a ) ，用f 2 ( p ，吒q ) 来代替e 1 2 ( p ，口，q ) 。 我们把含有所有t i j e p 的最小区间等分成厶份，其中0 三t o t l ( 2 - s ) 已经 的偏 参数 体光 我们 们不 极大 。它 样本 整体 体的 程是 证明 进是 量足 求出 第二章带有广义伽马分布厂，f ，项的加速失效时问模型 l _ l _ - l i i - - l _ - l _ l i _ - _ - - i - - _ _ _ _ - _ _ _ _ _ _ _ _ - _ _ - _ _ - _ _ - _ - _ _ - l - - _ _ - - _ _ _ _ _ - _ - _ 一i i i - - - l o ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! o “ 0 23 r e g r e s s i o np a r a m e t e r i c ) 023 r e g r e s s i o np a r a m e t e r o n 0123 r e g r e s s i o np a r a m e t e r 0 23 r e g r e s s i o np a r a m e t e r 图2 - 1 偏似然函数的核光滑处理效果示意图 1 5 一 o n n i o 1 n o o o 可ooc=m卫=jo一卫#o_i正 n n o n m 1 _ o p o oo吨o 口ooc=学暑，o一卫ij星 崛i o r 竹o 可oo点一心兰|-西oi o 嚣芎 d 峤p o 1 _ n o 口ooc=mj=口oi口一譬口i正 2 2 估计过程 其中 经过对似然函数的核光滑处理，我们的似然函数被转化成了如下形式： f 8 ( p ，盯，q ) = l l ( a ，q ) + 呓( p ，盯，口) 州叫，= 妻i = 1 就j = l ( o s z 汕g 入) + q - 2 l o gq - 2 + 扩2 - 1 ) 劫 一q - 2 e ( z i ) 一l o gr ( q 一 以p q ，= 妻i = 1 妻j = l 嘣。g 一嘣。g 去喜喜，e ( 半) k c 州s ) ) 上式的r i j = l o g t i j p z 巧，而k ( s ) 为选定的核函数，a 。为核函数的窗宽。 此时，似然函数被转化成了可估形式。我们可以极大化似然函数f 8 ( p ，盯，q ) ，以求 得参数p ，盯，q 的极大似然估计。 第二章带有广义伽马分布f l ( 1 i l q 顶的加速失效时间蠼型 第三步( m s t e p ) ：同时对p ，盯，q 极大化似然函数1 8 ，求出参数的极大似然估计 值，并依据此估计值更新h i j ( o o ) 的估计。 第四步：迭代第二步和第三步，直到参数口，t 7 ，q 同时收敛。 2 3 数值模拟 我们用冗语言进行编程，用模拟的方法来检验我们算法的有效性。我们采用参数 形式的加速失效时间f r a i l t y 模型来产生数据。在参数模型中，假设风险函数h ( t ) 满足 某种经典分布。那么我们可以通过确定风险函数和f r a i l t y 项的分布函数推导出个体实 际的生存时间。当然生存数据的特殊性在于它是一种随机删失类型的数据，对于删失 的处理，我们利用f 0 ，a 1 上面的均匀分布来对删失时间c 进行模拟，然后通过与真实生 存时间之间取小，实现对数据的删失。这里的a 是一个变量，我们可以通过控制a 的大 小，来实现对数据删失率的控制。 具体的模拟过程中，我们分别模拟风险函数服从对数正态分布和韦伯分布。 通过这两种不同的分布函数，我们可以考察风险函数的变化是否会对我们方法的 估计效果产生影响。对于f r a i l t y 项的分布，我们假设其服从对数正态分布，韦伯 分布和伽马分布。这三种分布都是广义伽马分布的特殊情况，它们也是生存数 据最常见的三种分布。也就是说，我们在估计过程中要考虑六种不同的参数模 型，以考察我们的算法在假设分布出现波动的时候，是否依旧可以给出精确的模 型参数估计值。当然，对参数p 的估计不仅仅依靠参数模型的假设分布，数据的 删失率对参数估计的影响也很大，所以，我们通过调整a 的值，分别考察删失率 为7 0 和8 5 两种情况下本章方法的估计效果。最后，我们需要考虑不同的群组数 量是否会影响估值方法的稳定性。基于这样的考虑，我们考察三种群组情况，分别 为n ：2 0 0 ，m = 2 ：r t = 2 0 0 ，n t = 4 和t t = 4 0 0 ，n t = 2 ，这里的n 为组数，n i 为每组个 体数。通过对以上一系列的情况的考察，来评估我们算法的稳定性和精确性。 除此之外，我们假设模型参数p 的真实值为( 0 ，1 ) t 。把协变量x 【，l 】设为常 量1 ，x | 2 1 为服从均值为2 ，方差为0 5 的正态分布。由于f r a i l t y 项的分布参数决定 了该分布的期望值，为了让其满足我们模型的假设前提条件，我们需要对其分布函数 里面的参数进行一些假设。当f r a i l t y 项服从伽马分布的时候，我们假设伽马分布的形状 参数为昙，而它的尺度参数为2 。在这种假设前提下，f r a i l t y 分布的期望值为1 ，方差 为2 。当f r a i l t y 项服从韦伯分布的时候，我们假设它的形状参数为1 ，尺度参数也为1 。 在此假设条件下，f r a i l t y 分布的期望为1 ，方差也为1 。在对数正态的情况下，我们假设 它的参数m e a n l o g = 一1 ，而它的另外一个参数s d l o g = 2 ，在这种前提下，f r a i l t y 分 布的期望为1 ，而方差为6 3 9 。通过我们的一系列假设，我们严格的把f r a i l t y 分布的期 一1 7 2 3 数值模拟 望限定为1 ，这是符合背景和逻辑的。而三种分布的不同方差也可以让我们从侧面考 察f r a i l t y 分布的方差对估计有何种影响。 在对似然函数的处理过程中，我们会用核光滑的方法光滑似然函数。因此我们 也要对核光滑方面的参数进行选择。在模拟中，我们采用标准正态分布的密度函数 作为核函数。另外一个重要的参数就是窗宽。在核光滑的处理过程中，窗宽的选择 对估计的会产生很大的影响。j o n e s ( 1 9 9 0 ) 和s h e a t h e r ( 1 9 9 1 ) 曾经给出过两个最优窗 宽：( s v s 3 ) 1 5 口1 n - 1 5 和4 1 3 a 2 n 。3 。其中o r l 为以j = 1 的数据计算出的l o g t z x 的 样本方差。盯2 为全体数据计算出的l o g t p x 的样本方差。我们在实际的操作中取两者 中的小者，作为我们的最终窗宽。 当然，为了评估方法的优越性，我们需要找寻其他方法作为比较。但是由于模型 形式之间的差异，我们把对其他模型提出的估计方法作为评估对象显然是有失公允而 且效率低下的。所以我们考虑r 中的固有模型估计方法作为对比对象。我们分别把独 立的加速失效时间模型和可加的加速失效时间f r a i l t y 模型( l a m b e r te ta 1 ，2 0 0 4 ；c o l l e t t ， 2 0 0 3 ) 的估计值作为对比。在下面的章节中，我们用g e m 来代表本章建议的估计方 法，a f t 来代表独立的加速失效时间模型，a f t f 来代表可加的加速失效时间f r a i l t y 模 型。我们对每种情况模拟生成5 0 0 个不同的数据集，用这5 0 0 次的估计值来计算最终 的m s e 和b i a s 。 表2 1 到表2 6 分别列出了各种情况下，上述三种方法的参数估计值。通过估计结 果我们不难看出，独立的加速失效时间模型的估计精度基本上落后于另外两种算法。 但是可加的加速失效时间厂r o i l t y 模型虽然在大多数情况下估计效果较好，但是它的稳 定性很不理想。在真值为1 的情况下，估计偏差有时候会超过土o 5 ，甚至出现了偏差达 到- 1 2 7 5 2 的情况。而纵观整个估计过程，本章所建议的方法在不同情况下的参数估计 值远远优于另外两种方法。m s e 和b i a s 在三种方法中基本都是最小值。这不仅说明我 们的估计方法有着极强的稳定性，也说明了本章方法的有效性和合理性。虽然估计的 效果非常理想，但是由于本方法中涉及到了对后验分布的抽样，以及要求e m 算法收 敛。这两个步骤严重影响了我们估计的效率，我们的方法要运行相当长的时间才能使 参数达到收敛。这也促使我们改进算法。因此在下个章节，我们会阐述一种全新的参 数估计方法，它的估计效率远远高于本章的方法。 1 8 第二章带有广义伽马分布f