（概率论与数理统计专业论文）广义线性模型中的参数估计问题.pdf

中文摘要广义线性模型在理论和实际应用中均有十分重要的意义，关于这方面的理论成果广泛应用于生物、医药等领域因此对广义线性模型问题的研究是很重要的，而在这类问题的研究中参数估计一直是大家关注的一个焦点问题之一参数估计在对模型作统计推断时是至关重要的方面，但是多重数值积分的问题带来了很大的困难，提出能避开数值积分困难的参数估计方法也就成了其中的研究重点本文主要讨论了广义线性模型中经常使用的一些参数估计方法，分别讨论了不带随机效应的广义线性模型和广义混合线性模型中参数估计问题，对这些方法给出综述和评价，并通过模拟计算对这些方法进行比较，模拟结果表明了m c e m 算法和m c n r 算法的优良性关键词：极大似然估计随机效应m c e mm c n rs m l 分层广义线性模型h - 似a b s t r a c tg e n e r a l i z e dl i n e a rm o d e l sh a v eu n i f i e dt h ea p p r o a c ht or e g r e s s i o nf o raw i d ev a r i e t yo fc o n t i n u o u sa n dd i s c r e t ed a t a t h el a t t e ri np a r t i c u l a rh a sav e r yi m p o r t a n ts i g n i f i c a n c ei np r a c -t i c a la p p l i c a t i o n t h e r e f o r et h es t u d yo ft h eg e n e r a l i z e dl i n e a rm o d e li sv e r ym e a n i n g f u l ，t h e nt h ep a r a m e t e r se s t i m a t i o ni nt h es t u d yo fs u c hp r o b l e m sh a sb e e naf o c a lp o i n to fc o n c e r nt oe v e r y o n e p a r a m e t e r se s t i m a t i o no ft h em o d e lw a sc r u c i a la s p e c to fs t a t i s t i c a li n f c r e n c e 。e x -t e n s i o n so fg e n e r a l i z e dl i n e a rm o d e l st oi n c l u d er a n d o me f f e c t sh a s ，t h u sf a r , b e e nh a m p e r e db yt h en e e df o rn u m e r i c a li n t e g r a t i o nt oe v a l u a t el i k e l i h o o d s b a s e do np r e v i o u sr e s e a r c h e st h i sa r t i c l eg i v eas u m m a r yo f s o m ep a r a m e t e r se s t i m a t i o na l g o r i t h m si ng e n e r a l i z e dl i n e a rm o d e l s ，a n de v a l u a t e d p a r t i c u l a r l y , t h o s ei ng e n e r a l i z e dl i n e a rm i x e dm o d e l s l a s t l y , c o m p a r i n gt h e s em e t h o d sa n dm a k ear e a s o n a b l ea s s e s s m e n tt h r o u g hs i m u l a t i o n s k e yw o r d s ：m a x i m u ml i k e l i h o o de s t i m a t i o n ；r a n d o me f f e c t s ；m c e m ；m c n r ；s m l ；h i e r a r c h i c a lg e n e r a l i z e dl i n e a rm o d e l s ；h - l i k e l i h o o d ；独创性声明本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的研究成果据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得东北师范大学或其他教育机构的学位或证书而使用过的材料与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意学雠文作者始j 鬈弛日期沙田学位论文版权使用授权书本学位论文作者完全了解东北师范大学有关保留、使用学位论文的规定，即- 东北师范大学有权保留并向国家有关部门或机构送交学位论文的复印件和磁盘，允许论文被查阅和借阅本人授权东北师范大学可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印，缩印或其它复制手段保存、汇编学位论文学位论文作者签名j 磊l ! 芝盟指导教师签名：垄! 筮至日期! = = 2 1 ：罗日期：之竺窆：! 妒学位论文作者毕业后去向：工作单位：通讯地址：电话：邮编：第一章引言自从1 9 7 2 年n e l d e r 和w e d d e r b u r n 提出广义线性模型这个概念到现在已经有3 0 多年历史了，广义线性模型的研究已经发展到了高峰时期，是统计研究当中一个经典的研究课题，研究文献，专著数以万计作为正态线性模型在形式上的直接推广，广义线性模型在很多方面都得到了广泛的应用，尤其是在离散型数据的问题处理中发挥了巨大的作用，有很大的实际意义如生物，医学，教育，经济和社会中的属性数据和记数数据等，广义线性模型为研究这些方面的问题提供了一个很重要的工具，理解和掌握广义线性模型的参数估计方法对实际工作中的问题解决是一个强有力的帮助另外，在广义线性模型中利用联系函数把不是线性关系的变量转变成线性模型条件下作统计研究，极大地方便了研究工作，从而使一些复杂的混合效应问题得到很好的解决而在广义线性模型的研究问题当中，参数的估计问题更是大家关注的焦点问题研究广义线性模型的参数估计问题有很高的学术价值和应用价值因为优良的参数估计是进行统计推断，对未来进行预测的主要依据，为决策提供更可靠的保证寻找更高效的算法和估计方法在实际工作中提高了效率和准确度在参数的估计方法中应用最多的就是极大似然估计。但很多时候似然方程的求解是一个很困难的问题，我们只能通过迭代算法来求其近似解，这时n e w t o n r a p h s o n 算法和e m 算法就是我们主要使用的基础工具现代计算机技术的快速发展，解决了许多很难的数值计算问题，为广义线性模型的参数估计问题研究带来了更加有利的条件，促进了这个领域的发展使得广义线性模型中很多复杂的参数估计值的迭代方程求解问题得以解决一开始，国内外这方面的研究工作者大量的工作都是在集中在这样的一类广义线性模型上t 对模型中的随机效应成分作正态性假设这当然使得研究工作变得更加方便，但也限制了它的发展，因为实际中很多的数据都不是满足正态性假设的学科的发展和实际需要促使人们把注意力集中在不是正态随机效应的广义线性模型上，广义线性混合模型就在这样的背景下产生的其中很有实际意义的这样一类广义线性混合模型就是l e e , z和n e l d e r 在1 9 9 6 年提出的分层广义线性模型( h i e r a r c h i c a lg e n e r a l i z e dl i n e a rm o d e l s ) 以及在2 0 0 6 年提出的双分层广义线性模型( d o u b l eh i e r a r c h i c a lg e n e r a l i z e dl i n e a r m o d e l s ) 由于在实际应用中的重要意义，对这样的广义线性模型的参数估计问题的研究也就成了热点问题，解决这类广义线性模型的参数估计方法是很有应用价值的研究方向广义线性模型的参数估计问题很重要的一个方面就是如何找到估计方程，寻找有效的算法求解理论方面，由于证明都比较复杂，这方面的文献都不涉及严格的数学推导，关于这方面的问题可以参考文献【1 】i 和【1 l 】文章共分为五部分，第一章为引言部分；第二章介绍不带随机效应的广义线性模型的参数估计问题；第三章主要讨论广义线性混合模型中的参数估计问题，重点介绍了m c e m ，m c n r ，s m l 这几种主要算法；第四章介绍了一类重要的广义线性混合模型一分层广义线性模型；第五章为计算机模拟部分，对算法进行模拟比较和评价2第二章不带随机效应的广义线性模型的参数估计l 模型简介n e l d e r 和w e d d e r b u r n ( 1 9 7 2 ) 提出了广义线性模型，这里从简单的模型说起假设y “，为相互独立的一维响应变量，期，而为自变量，而一般为多维变量，弘的分布属于指数型分布m ；劝= c f v j ) e x p f v f ) 一“b ) )( 2 1 )其中d ) 和6 ( ) 为已知的具体函数，岛称为自然参数，且有妇) 强= 昱( y d = ) 蜥= 踟r ( y )( b ) e t y t ) = 胁= “西，j i 为严格单调，充分光滑的函数，g = h - 1 称为联系函数引进记号协= 搬则在模型中，( 2 1 ) 式通过r i 而依赖届从而对卢进行统计推断2 参数极大似然估计早在1 9 8 9 年，m c c u l l a g h 和n e l d e r 就在专著广义线性模型中详细的讨论了参数的极大似然估计在这里主要参考陈希孺( 见参考文献【l 】) 给出的似然方程和迭代算法从( 2 1 ) 式出发得到似然函数t从而对数似然为工= hc y 1 ) e x p y i o l 一姗) l扛1i o g l = e l o g c o i ) + 一) )( 2 2 )因为c ( v i ) 对卢的估计无影响，且研与x i 和卢有关，为了突出0 和卢的关系，( 2 2 ) 式可写成托回= ：眦岛仰一6 ( 毋仰) j( 2 3 )百进一步得到似然方程，筹= 叫c 岛鬻- o方程的解就是参数口的m l e 接下来的问题就是方程的求解问题，但是通常情况下，方程( 2 4 ) 往往只能用迭代的方法来求解陈希孺( 见参考文献f 4 】) 提出了一种n e w t o n - r a p h s o n 迭代的算法算法基本步骤如下3l 。设卯，楞为矶，粕的当前值，算出：0 ) = g - 1 ( ) ，i = 1 棚的值2 。利用撵本y i ，k ( 它在迭代过程中当然是不变的) 及上述拶和的值，算出帮= ，7 0 ) + ( k p ：0 ) 珂0 ：0 ) ) i = 1 卅，的值3 。以而记自变量在第i 次观察中所取的值，找晟使表达式岫( z ：0 ) 一( 2 5 )括l的值达到最小，其中u 产( g ，：o ) ) ) 一2 ( - ：0 ) ) y 叫，( 2 5 ) 式的解，设为o ( 这就是一个以x 为自变量，z 为因变量，1 ，为权的加权最小二乘估计) 4 。由o 产生们，柏的更新值t毫”= 印( o ) f - 1 ，然后又回到第一步，直到口的值变化很小为止，最后得到的卢就是所要求的极大似然估计值直上述算法的理论依据见参考文献 4 】，这里就不加讨论了对于不带随机效应的广义线性模型，由于没有随机效应的影响，参数估计相对来说比较容易实现其中得分算法是最为重要的基本工具在第五章中我们将用一个简单的例子来模拟有关这类模型的参数估计，验证算法的可行性4第三章广义线性混合模型的参数估计l 广义线性混合模型在第二章中我们讨论了不带随机效应的广义线性模型的参数估计问题，但是在实际应用当中，我们经常遇到的是带随机效应的广义线性模型，而且随机效应不作正态性假设，就是这一章要研究的广义线性混合模型研究广义线性混合模型这这样的模型是很重要的，而且很有实际意义但是在这种情形下不可观测的随机效应给参数估计带来了很大的困难，使得广义线性模型的发展受到了限制在很多情况下广义线性混合模型的参数估计问题都要涉及到高维的积分问题，这是很难用一般的数值计算方法解决的如何解决这些计算的问题，避免直接计算高维积分的困难? 因此这方面的研究工作受到越来越多人的关注，成了广义线性模型这个领域研究的重点这一章我们就详细讨论广义线性混合模型中的参数估计问题假定y = ( y l ，k ) 为响应变量，巧之间相互独立，u 为不可观测的随机效应部分，为了使问题具体化，假定i i 来自参数为d 的某分布，i l l 一五( u i d ) ，在给定的情况下的，y的条件分布属于指数型分布( 便于考虑，这里联系函数取为自然联系函数) 氐( y f l u ，励= e x p t v i 聃一c ( 聃) + d ( y i ) l( 3 1 )其中，琅= p + ，分别为矩阵x 和z 的第i 行从模型( 3 1 ) 出发可得到似然工伊，o l y ) - - f 兀瓜协i u ，跳( l d 胁( 3 2 )v持l通常情况下，由于i l l 是多维的，( 3 2 ) 式这个积分很难求解，这就给参数估计带来了不可避免的麻烦如何用有效的方法来求解基于( 3 2 ) 式的参数极大似然估计就是这类广义线性模型中的关键问题如果能够解决这个问题将是很有意义的，因为随机效应带来的问题是一类很广泛的研究问题，而且在实际应用当中具有重要的实用价值针对这种情况。广义线性混合模型的参数估计问题在广义线性模型这领域中成了研究的焦点问题，下面将重点讨论这一类模型的参数估计问题2 参数估计为了解决( 3 2 ) 式中的参数估计问题，很多工作者都致力于寻找有效准确的方法得到参数的极大似然估计m c c u l l o c h ，c e ( 1 9 9 7 ) 提出了三种有效的算法很好的解决了模型( 3 1 ) 中的参数估计问题，分别是m o n t ec a r l oe m ( m c e m ) 算法，m o n t ec a r l on e w t o n 一5r a p h s o n ( m c n r ) 算法和s i m u l a t e dm a x i m u ml i k e l i h o o d ( s m l ) 算法，下面我们详细讨论这几种算法2 1m c e m 算法从e m 算法的思想出发，把随机效应u 当作缺失数据w 记为完全数据，则w =( ku ) 完全数据下的对数似然可由下式给出ti n l w = i n 矗。挑l u ，历+ i n 五( u l d )( 3 3 )7这样选择缺失数据具有两个优点。首先，在给定u 的情况下，嚣是相互独立的；其次，e m算法中的m 步相应于参数口和d 是分开的，即求最大值时可以分两部分单独进行计算根据这种情况，按以下步骤给出e m 算法t1 。选择初始值和d ( o ) ，令m = 0 2 。计算( 在”和d 【”下计算) ta 妒“) 使得e t n 名o l u ，国l y 达到最大b d ”使得e h a ( u l d ) l y 】达到最大c 令m = m + i 3 。如果达到收敛，则”1 和d c m “) 为极大似然估计，否则返回第二步一般情况下，上面步骤中的2 a 和2 b 中的期望值很难从模型( 3 1 ) 中求出，因为条件分布f ( u t v ) 的计算包含了石，面菇是未知的为了避免直接计算的困难，采用m o n t e c a r l o方法去近似我们想要知道的期望值在这里介绍的方法中，我们使用的是m e t r o p o l i s 算法m e t r o p o l i s 算法的本质想法是，目标分布( 又叫靶分布) 的样本不好产生，我们选定一个候选分布钆。从h 。中抽取样本，以接受概率函数接受接受从玩中抽取的样本，一直到产生完我们需要的样本个数结束。这样就得到了来自目标分布的一组样本，避免了从目标分布中直接抽样的困难在这里选择h u = 工作为候选分布，接受概率函数a k ( u ，u + ) 记u _ ( u l ，u 2 ，u k l ，啦，u q ) ，从而在算法中以概率a k ( u ，i i + ) 接受i i + ，否则保留u 这里4 c u ，u 。，= m i n t - 。； ：；：；言篇，c ，。，在这种选择下，( 3 4 ) 式中的后面那项可以简化为下面的形式t讹刊n u ，然瑞篇，6：圣! 望五! 坐竖：! 堕! ! ：! 旦堕尘! 旦!兀冬li n j ( y i i n ，卢埔( u l d 肌( u i d )：畦! ! 竺盘坐坠旦兀缸l u ，励( 3 5 )从( 3 5 ) 式可以看出，经过选择以= 工作为候选分布，接受概率函数变得容易计算，使我们的算法可以进行m e t r o p o l i s 算法用在这里的作用就是为了避开条件分布f ( u t v ) 的计算，这样就可以不用知道工的具体形式我们用m e t r o p o l i s 产生n 个u 的值，等价于从f 咖( u l y ，反。) d ( 帕) 中选取了n 个u 的值，然后用m o n t ec a r l o 近似代替2 a 和2 b 中的期望值就解决了e 1 n f 坤, ( y l u ，p ) i y l 和e b n a ( u l o ) l y 】的计算问题把m e t r o p o l i s 算法加入到e m 算法中，得到m c e m 算法的基本步骤如下，1 。选择初始值f 町和d ( 0 ) ，令m = 0 2 。使用m e t r o p o l i s 算法从f 咖( u l y ，鼬，d o , o ) 中产生n 个i l l 的值，n ( ”，i m ( 2 ) ，u ( 帅：a 选择卢1 ) 使得e 1 0 矗。o ，i ，历| y 】达到最大，也即是使l ：i n 加0 , 1 u 噼，卢)( 3 6 )雷达到最大b 选择d ( “1 使得l ，鉴ll n 工( u p ) 达到最大c 令m = m + 1 3 。如果达到收敛，则”和d ( m “为极大似然估计，否则返回第二步以上就是m c e m 算法的主要步骤，主要思想就是把m e t r o p o l i s 抽样方法应用到e m算法中，抽取不可观测随机效应成分u 的值使得e m 算法中极大化时很难求解的期望值变得容易计算，从而实现e m 算法的程序，求出我们想得到的参数的极大似然估计( m l e ) 从后面的模拟结果我们可以看出，m c e m 算法得收敛速度比较快，准确率也比较高，这说明这个算法对广义线性混合模型的参数估计还是比较好的一个方法，但是要达到很高的准确度的需要的运行次数还是很大的2 2m c n r 算法在不带随机效应的广义线性模型的参数估计问题中，对于固定效应的估计方程常常很难用一般的数值解法得到，所以n e w t o n r a p h s o n 方法是经常使用的迭代算法求解参数的极大似然估计自然的，n e w t o n r a p h s o n 方法在广义线性混合模型中同样是重要的工具，如何把n e w t o n r a p h s o n 方法应用到这一类模型的参数估计当中是很有意义的7研究问题m o n t ec a r l on e w t o n - r a p h s o n ( m c n r ) 算法就是在这种背景下产生的一种合成算法下面讨论算法的具体实现在模型( 3 1 ) 中，考虑到参数口只和a l u 这一部分似然有关，从n e w t o n r a p h s o n方法的基本思想出发，( 3 2 ) 式表示的似然是独立的两部分，矗。和工，因此，参数卢和d的极大似然估计方程可由以下形式给出；( 3 7 )和研等驴| y 】= 0( 3 8 )( 3 ，8 ) 式很容易由u 的分布求解，因为在的分布给定情况下，( 3 8 ) 式的期望值就只是关于d 的方程了但是对于( 3 7 ) 式的求解，一般情况下是很困难的，经验表明大多数情况下( 3 7 ) 式的求解都需要用到迭代算法，经常使用的方法就是n e w t o n 一只印h s o n 算法所以在这里我们提出的m o n t e c a r l o n e w t o n - r a p h s o n ( m c n r ) 算法也是以n e w t o n - r a p h s o n方法为基本工具，结合m o n t ec a r l o 方法实现程序的进行，这也找到了参数卢的极大似然估计的一种有效的估计方法n e w t o n r a p h s o n 方法的基本思想出发，把业学作为卢的一个函数在风处展开to l nj ；1 l t ( y l u ，固。, 9 i n a o ( y l u , 芦 ) lr2 筇塘铺+ 訾伊圳( 3 39 一)+ 鼍= - = 一l r ：m l 口一f h li 1根据m c c u l l a g h 和n e l d e r ( 1 9 8 9p 4 2 x 见参考文献【11 】) 中的结果，我们可以把上式写成，其中，0l n 厶历u ( y l u , 3 ) ! x w o ，u ) 享b 确【y 一乒愉，t o o t删搿w l 岛，u 瑚p 一廓)胁假u ) = e ( k i u )8( 3 1 0 )喇一眠t o = d i a g l ( a | l i l a ) 2 9 a , _ ( r d u ) l却，舡= d i a g & h a l u ，l用( 3 1 0 ) 式近似代替( 3 8 ) 式，得到迭代方程t芦( 肿1 ) = 哪+ e 【x ，w ( 秒) ，u ) x l y _ 1x x e 【w ，i i ) 宝b 妒，( ) ，一p b ，u ) ) l y 】( 3 1 0 )在通常情况下，( 3 1 0 ) 式中的期望同样是很难求解的，运用m e t r o p o l i s 算法，得到m o n t e c a r l o n e w t o n - r a p h s o n 算法，具体步骤如下l1 。选择初值郇和d p ，令m = 0 2 。运用m e t r o p o l i s 算法从f ( u t v ，呻，d 呻) 中产生( ”，n ( 2 ) ，柙个值，并用它们得到( 3 1 0 ) 式中期望的m e n t e c a r l o 估计( 记为自：1 ) ：a 计算“) = 户+ 营【x w 0 日( 埘) ，u ) x l y x x 懂 w ( 眇，u ) 象忙( y p ( 矽，u ”l y 】( 3 1 1 )b 选择d 1 使得1 篓ll n f ( u * j l d ) 达到最大c 令= m + 1 3 。如果达到收敛，则”，d ( ”) 为m l e ，否则返回第2 步n e w t o n 一见脚鲫l 方法是一种很有效的迭代算法，m c n r 方法同样对参数的估计是很有效的算法在5 中的模拟结果，我们也可以看到，m c n r 方法的收敛速度比m c e m方法还快，准确程度也差不多，但是要达到比较高的精确度同样需要运行次数非常大2 3s i m u l a t e dm a x i m u ml i k e l i h o o d 算法m c e m 和埘c m r 方法都是从t o g 似然函数出发进行的估计g e y e r 和t h o m p s o n ( 1 9 9 2 )以及g e t f a n d 和c a r t i n ( 1 9 9 3 ) 建议从直接从似然出发进行参数估计从( 3 2 ) 式开始工o t 。l y ) = f y 从y l u ，1 日，五( u i v )= 厂趔篇塑姒u 胁j丸( u )”“一9s 专姜趔铲h u8 忉( u ( ) )其中i l l 从重要抽样分布 。( u ) 中选取，为模拟值的个数，不管轧( u ) 的分布怎样选取，( 3 1 2 ) 式给出了似然的一个无偏估计对( 3 1 2 ) 式极大化可得参数的极大似然估计以上我们就主要讨论了m c e mm c n r ，s m l 这三种算法的具体实现步骤，但是一些理论上的根据就没有加以讨论对以上的方法作一些改进，我们能得到合成的算法，如m c n r + s m l ，m c n r + 2 s m l 等，这样的合成算法提高了估计的精确度限于篇幅这里就不作深入探讨了1 0第四章分层广义线性模型1 模型l e e 。y 和n e l d e r , j , ：l ( 1 9 9 6 ) 提出了一类重要的广义线性混合模型一一分层广义线性模型( h i e r a c h i c a l g e n e r a l i z e d l i n e a rm o d e l s ) 这类模型拓宽了广义线性模型的应用范围，更具实际意义，是很值得研究的问题记y 为响应变量，4 为不可观测的随机成分，满足：( 口) y 在给定的条件下的l o g - 似然具有以下广义线性模型的形式：i ( o ；y l 4 ) = 一6 ( 矿) + 枷( 4 1 )其中为自然参数，记，= e ( y l u ) ，= g m ) ，矿= q + v( 4 2 )其中”= 邪，v = v ( 4 ) 为“的严格单调函数( 6 ) “的分布是合理的假定的某一分布就称之为分层广义线性模型( h g l m ) 2h - 似然( h i e r a r c h i c a ll i k e l i h o o d )对于模型( 4 1 ) ，l e e , y 和n e l d e r ( 1 9 9 6 ) 定义了h 一似然，记为h ，由以下形式给出th = ，( 矿：卅v ) + 口；v )( 4 3 )( 4 3 ) 式把( 4 1 ) 中的似然拆成两部分，其中比口；力是v 在参数口下的对数密度函数，1 ( o ；y l v )为，i v 的对数密度函数由于v 的不可观测性，缸似然不是传统的似然极大化h 似然( 4 3 ) 得到参数卢和v 的估计，我们称之为极大h - l i k e l i h o o d 估计( m h l e ) ，它们可从以下的似然方程求得。o h 帮= 0( 4 4 )o h 加= 0( 4 5 )从( 4 4 ) ，( 4 5 ) 中得到m h l e 已经被证明了是渐近到m l e 的，在大样本情况下，这就提供了一种解决广义线性混合模型参数估计问题的方法3 共轭分层广义线性模型在分层广义线性模型中，有一类很重要的模型，就是共轭分层广义线性模型，这里我们用数据模型来研究它们记y i ( j = i 2 ，i = l ，2 ，t ”= m s ) 为响应变量的观测数据，嘶为不可观测的随机效应成分在自然联系函数下，2e u + v i其中吃= 舢；) ，劬= 毗u ) ，及v j = 战嘶) ，( y l u ) 的分布为“的共轭分布在共轭分层广义线性模型中，( 4 4 ) 式变为l a 艮h = i , jl 厂螨溉h - 似然中的( 口；具有以下的核形式ty l q m 一啦幢) “v f ) )厶_ _ ，“2 “，”、。，其中q i ( ) 和口2 ( ) 是。的函数，从而h 一似然可写成以下的核形式- 的一取州+ a ) v j 一啦( a ) 6 ( ) 由a 6 ( 眦) ) 1 0 0 = 卢可知o b ( v ) a v = “，故从( 4 8 ) 式出发有关= 叫十) + 州妒眈，舢i 卜。”。一j从方程o h a v i = 0 ，即( 4 1 0 ) 式等于零可解得随机效应u 的m h l e) + 一p ，+ + a l ( o )砷2 瓦丁( 4 6 )“7 )( 4 8 )( 4 9 )( 4 1 0 )( 4 1 1 )其中y i + = z 您，芦，+ = p ；，。这表明在共轭h g l m 中，随机效应的m l e 具有简单的表达形式接下来从几个具体的例子来介绍共轭h g l m 3 1p o i s s o n g a m m a 模型如果( y l u ) 的分布为p o i s s o n 分布，= p “，u 的分布为g a m m a 分布。则在这样的共轭h g l m 中有= l o g , u ：= e i i + m其中0 = 郦，= l o g u i ，从而有，( 口，“f 口u + o ! l o gc r - a , u f - 1 。gr ( a ) j( 4 1 2 )1 2所以a l ( = 啦( = 口，由( 4 1 1 ) 式可得瓤= o “一口，+ + a , = ( 片+ + ) ( 声“+ 口) ( 4 1 3 )我们就得到了e ( u j = t y + 叻，o “+ 口) 的一个估计另一方面，关于卢的极大h 似然估计方程变为嘏o h2 善一鼢+ 岫= 。( 4 1 4 )从( 4 1 4 ) 式得到1 3 的m h l e 与卢的m l e 是等价的3 2b i n o m i a l b e t a 模型( ，l “) 的条件分布为b i n o m i a l 分布，= m x ，“的分布假定为b ，啦) 这就得到了一个共轭h g l m吒= o , j + 坼其中矿= l o g l x o 一一) ，= l o g l u i ( 1 一蜥) 0 = l o g , r ( 1 一j r ) ) = 班e ( “) = 口i ( 口l + n )在这个模型下，有v ) = h 坼呻- + a 2 ) 1 。g ( 击) 一l o g 弛) 】这里a l ( 曲= 口2 ( = m + 0 1 2 ，由( 4 1 1 ) 式可得关于卢的极大h 一似然估计方程为 y h 一以 + 8撕2 1 i 五丁差2 善一m , j g j ) x , , j 2 。( 4 1 5 )( 4 1 6 )还有g a m m a i n v e r s eg a m m a 模型，i n v e r s eg a u s s i a n g a m m a 模型等也属于共轭分层广义线性模型的范畴，这里就不详细讨论了因为随机效应在实际中是我们很关心的一个方面，很多情况下随机效应在模型的建立中是不能随便忽略的，而且由于实际的需要，也不能对它作正态性假设这就给模型的1 3参数估计带来很大的困难分层广义线性模型的提出为解决这类问题提供了一种很重要的建模方法和思路，使得广义线性模型在实际中发挥出更大的作用在分层广义线性模型中，对固定效应和随机效应的估计都是从h 一似然出发，但是似然方程在大多数情况下是没有显式解的e m 算法和得分算法依然是解决估计方程的主要工具由于分层广义线性模型属于广义线性混合模型的范畴，所以可以利用第三章中给出的方法进行参数估计m h l e 的性质在文献【1 3 】中给出了证明，限于篇幅，这里就不讨论了第五章模拟比较在这一章，我们将通过计算机模拟数据的方法对前面给出的一些方法进行比较分析，验证这些方法对参数估计的好坏这里我们两个简单的例子来做模拟ln e w t o n r a p h s o n 算法对于不带随机效应的广义线性模型，以l o g i s t i c 回归模型为例；l o g ( 尚一堆l从以上模型得到，d g 似然如下 y i l o g ( p ) + ( 1 - p i ) l o g ( 1 一翩) 】= y t o t + 卢x i ) - l o g ( 1 + 矿咖) 】知l- - i在模拟中取，= 5 0 ，即产生5 0 个随机数，乃的取值为0 ，1 这里分别对口，卢取不同的值进行模拟计算，采用n e w t o n - r a p h s o n 迭代的方法进行求解模拟结果见表一，表二，表三口真值口真值口卢a 估计均方误卢估计均方误n r 估计l0 2o 8 1 7 2o 2 5 1 40 3 7o 3 l表二a 真值芦真值口卢口估计均方误卢估计均方误i n r 估计3o 32 6 1 0 50 2 5 7 4o 3 2o _ 3 5表三口真值卢真值。卢口估计均方误卢估计均方误i n r 估计40 0 43 7 6 4 90 0 3 4 9o 3 8o 3 62m c e mm c n r ，s m l 算法m c e mm c n r ，s m l 算法都是广义线性混合模型的参数估计方法，这里对它们进行模拟比较以一个l o g i t n o r m a l 模型为例进行模拟计算，模型如下，k | u i n d e pb e r n o u l l i ( p q )1 5i = 1 ，2 ，h ；，= 1 ，2 ，q ，t n ( p q o p q ) ) = 肛“+ u j ，蜥一i i d n ( o ，一)在模拟中取月= 1 5 ，9 = 1 0 ，从上述模型中产生数据在m c e m , m c n r 算法中每次用m e t r o p o l i s 抽样方法抽取“的个数为1 0 0 在s m l 算法中，重要抽样分布的循环更新次数为5 0 0 0 次具体过程见第三章中的算法步骤下面四个表格分别给出了在抽取随机数时对只，取不同值得到的模拟结果表四够= 1 ，一= o 4 )芦铲岔的均方误萨的均方误m c e m1 2 1 3 50 4 8 5 lo 0 50 4 5m c n r1 2 6 5 70 5 0 6 7o 0 8o 5 9s m l1 8 1 3 4o 5 9 1 4o 1 9o 9 3表五够= 2 ，一= o 4 )卢萨口的均方误护的均方误m c e m1 8 7 4 60 4 7 6 20 0 4o 3 6m c n r2 2 3 4 90 4 9 8 30 0 70 4 5s m l2 7 6 7 50 6 0 7 lo 1 51 2 2表六= 2 ，c r 2 = o 5 )卢铲声的均方误萨的均方误m c e m2 1 1 4 50 3 8 5 7o 0 50 4 1m c n r2 3 3 5 6o 6 0 2 20 0 9o t 8 7s m l2 6 1 3 40 7 0 9 4o 1 31 0 4表七啦= 5 ，o - 2 = o 5 )芦萨p 的均方误萨的均方误m c e m4 9 1 0 60 5 8 6 lo 0 3o 2 9m c n r5 2 7 1 5o 6 1 3 70 0 5o 3 7s m l5 7 5 4 10 ，6 9 5 8o 1 8l 。1 21 6从模拟结果我们可以发现n e w t o n r a p h s o n 算法的收敛速度很快。但是很依赖于初值的选取，尤其在估计值是多峰的情况下很难得到好的估计在对广义线性模型做参数估计时，它表现的比较稳健，估计的效果不错m c e m 算法继承了e m 算法的收敛性质，在足够大的模拟样本情况下，保证了算法的收敛性，收敛的速度很快，估计的准确度也比较高，但是有可能出现不是收敛到真值的情况m c n r 算法的收敛速度是最快的，但有可能出现不收敛的情况在不是大样本的情况下，算法的收敛性有可能得不到保证单一的s m l 算法表现不好，离真值可能会相差较大，收敛速度也没前两种快，但是和m c n r算法合成后的表现非常好，估计的准确性非常高1 7参考文献 1 】陈希孺广义线性模型( 一) 至( 十) 专题讲座数理统计与管理，2 0 0 2 2 0 0 4f 2 】王松桂，史建红等线性模型引论 m 】北京：科学出版社，2 0 0 4f 3 】3 张尧庭。方开泰多元统计分析引论【m 】北京一科学出版社，1 9 9 7【4 】陈希孺近代回归分析 m 】【5 】陈希孺高等数理统计学【m 】，合肥t 中国科技大学出版社，1 9 9 9 6 1 张金槐线性模型参数估计及其改进【m 】长沙t 国防科技大学出版社，1 9 9 9 7 】何晓群现代统计方法与应用【m 】北京t 中国人民大学出版社，1 9 9 8 8 】8 王松桂e m 算法【j 】，应用数学与计算数学【9 】茆诗松，王静龙，濮晓龙，高等数理统计，高等教育出版柱，施普林格出版社，2 0 0 0 1 0 】王沫然m a t l a b 与科学计算【m j 北京t 电子工业出版社，2 0 0 3【11 】m c c u l l o g h ，pa n dn e l d e r , j a ，g e n e r a l i z e dl i n e a rm o d e l s ，2 n de d n l o n -d o n ：c h a p m a na n dh a l l ，1 9 8 9【1 2 】n e l d e r , j a a n dw e d d e r b u m g e n e r a l i z e dl i n e a rm o d e l s ，j r s s b ，1 9 7 2 。3 7 0 3 8 4【1 3 】l e e ，ya n dn e l d e r , j a h i e r a r c h i c a lg e n e r a l i z e dl i n e a rm o d e l s ，j r s s b 1 9 9 6 ，5 8 ，6 1 9 - 6 7 8【1 4 f a h r m e i rl ，k a u f m a n nh c o n s i s t e n c ya n da s y m p t o t i cn o r m a l i t yo f t h em a x -i m u ml i k e l i h o o de s t i m a t o ri ng e n e r a l i z e dl i n e a rm o d e l a n ns t a t i s t ，1 9 8 5 ，1 3 ：3 1 2 - 3 6 0【1 5 】c h e nk ，e t s t r o n gc o n s i s t e n c yo f m a x i m u mq u a s i l i k e l i h o o de s t i m a t o r si ng e n e m l 呓e dl i n e a rm o d e l sw i t hf i x e da n da d a p t i v ed e s i g n a n ns t a t i s t , 1 9 9 9 ，2 7 1 1 5 5 1 1 6 31 8【1 6 】m c c u l l o c h ，c e m a x i n l n n l l i k e l i h o o d a l g o r i t h m s f o r g e n e r a l i z e d l i n e a r m o d -e i s j a s a ，1 9 9 7 ，9 2 ，1 6 2 1 7 0【1 7 】k a r i m ，m 1 la n dz e g e r , s l g e n e r a l i z e dl i n e a rm o d e l sw i t hr a n d o me l -f e e l s ：s a l a m a n d e r m a t i n gr e v i s i t e d b i o m e t r i c s ，1 9 9 2 ，4 8 ，6 3 1 - 6 4 4【1 8 】s c h a l l ，i le s t i m a t i o ni ng e n e r a l i z e dl i n e a rm o d e l sw i t hr a n d o me f f e c t sb i o m e t r i k a , 1 9 9 1 ，7 8 ，7 1 9 7 2 7【1 9 】l e e ，ya n dn e l d o r , j a d o u b l eh i e r a r c h i c a lg e n e r a l i z e dl i n e a rm o d e l s j a p p l s t a t i s t2 0 0 6 ，5 5 ，i - 2 9【2 1 】l e e , y a n dn e l d e r d a h i e r a r c h i c a lg e n e r a l i z e dl i n e a rm o d e l s ：as y n t h e s i so fg e n e r a l i z e dl i n e a rm o d e l s ，r a n d o me f f e c t sm o d e l sa n ds t r u c t u r e dd i s p e r s i o n s ,b i o m e t r i k a , 2 0 0 1 ，8 8 ，8 7 1 0 6f 2 2 】z e g e r , s l l o n g i t u d i n a ld a t aa n a l y s i su s i n gg e n e r a l i z e dl i n e a rm o d e l s ，b i o m e t r i k a , 1 9 8 6 ，7 3 ，1 3 2 2【2 3 】m c g i l c h r i s t , c a e s t i m a t i o ni ng e n e r a l i z e dm i x e dm o d e l s ，j r s s b 1 9 9 4 ，5 6 ，6 1 - 6 9【2 4 】g u t i e r r e z ，p ea n ds m i t h ，a e m c o n j u g a t ep a r a m e t r i z a t i o n sf o rn a t u r a le x p o n e n t i a lf a m i l i e s j a s a ，1 9 9 5 ，9 0 ，1 3 4 7 1 3 5 6 2 5 】l e e , y c a nw cr e c o v e ri n f o r m a t i o nf r o mc o n c o r d a n tp a i r