广义线性模型_六_.pdf

《广义线性模型_六_.pdf》由会员分享，可在线阅读，更多相关《广义线性模型_六_.pdf（10页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、文章编号:1002-1566(2003)04-0055-10广义线性模型(六)陈希孺(中国科学研究生院,北京100039)摘 要:本讲座是广义线性模型这个题目的一个比较系统的介绍。主要分3部分:建模、统计分析与模型选择和诊断。写作时依据的主要参考资料是L.Fahrmeir等人的:Multivariate Statistical Model2ing Based on Generalized Linear Modles。关键词:广义线性模型;建模;统计分析;模型选择和诊断中图分类号:O212 文献标识码:AGeneralized Linear ModlesCHEN Xi2ru(Graduate S

2、chool of Chinese Academia of Science,Beijing 100039,China)Abstract:This set of articles gives an introduction to generalized linear models.They can be divided into three parts:Modelbuilding,Statistical inference and Model diagnostics.The presentation in mainly based onL.Fahrmeir et al.MultivariateSt

3、atisticar Modeling Based on Generalized Linear Modles.Key words:generalized linear models;model building;statistical inference;model diagnostics213 拟似然法到此为止所有的讨论都是在Y服从指数型分布的假定下进行的。这个假定的根据是,我们的目的在于离散数据统计分析,而在一些应用上很重要的情况下,这种数据的分布是二项分布、多项分布、Poisson分布等,它们都属于指数型。但是,在有些情况下,“指数型”这个假定不一定切合实际:当建模时,往往着眼在变量的均值

4、、方差。像二项、多项、Poisson这些分布,都是单参型分布,其均值方差依赖于一个参数。因此方差2是均值的函数:2=(),()称为方差函数。例如对Poisson分布,()=。对二项分布B(m,),=m,2=m(1-)=(1-/m)。对负二项分布NB(r,),2=r+3+2/r等等。但实际数据有时不符合这个关系,如以前提过的所谓“超散布性”(见111末尾):如对二项分布,Ey=,但2可以大于1-/m(若YB(m,),则DY=m,而Var(Y)=m(1-)=(1-/m)。这时,可以证明,不存在一个取0,1,m为值(每个值的概率 0),服从指数分布,而对1满足Var(Y)=E(Y)(1-1mEY)的

5、变量。Y的分布称为负二项分布NB(r,)。55广义线性模型(六)每次试验成功的概率为,失败的概率为1-。给定自然数r,Y=第r次成功时已失败的次数,则P(Y=y)=r-1+yr-1r(1-)y,y=0,1,注:设有指数分布族P(Y=i)=c()d(i)ei,i=0,1,m,则c()0对一切 0,有c()=(mi=0d(i)ei)21E(Y)=c()=(mi=0id(i)eiE(Y2)=c()=(mi=0i2d(i)eiVar(Y)=c()=(mi=0i2d(i)ei-c2()(mi=0id(i)ei)2设Var(Y)=E(Y)(1-E(Y)/m),则mi=0d(i)ei mi=0i2d(i)e

6、i-(mi=0id(i)ei)=mi=0d(i)ei mi=0id(i)ei-1m(mi=0id(i)ei)2比较两边e的系数,有d(0)d(1)=d(0)d(1),得出 0,故d(i)0,i=0,1,m。因此d(0)d(1)0而必须有=1)对Poisson分布也有这个情况。当“超散布性”出现时,样本的均值方差不一致。这时就不能以Poisson分布为模型。以上讲的是这样一种情况:原来样本可以认为服从某种指数型分布,由于相依性及非齐次性(指在同一x之下多次观察的因变量Y并非同分布,因为还有重要因素未纳入x内。这些因素每次观察时取值可以不同(非齐性),因而使Y值有不同的分布)使指数假定不能成立。另

7、有一些情况,一开始就没有充分理由取指数型分布作为模型。这就说明:在实际问题中往往有必要在对变量的分布并未确切的情况下去建模,并发展出相应的统计推断方法。拟似然法就是为了这个目的。对我们此处的问题而言,拟似然法着眼在均值和方差,尤其是前者,即必须对均值有一个比较确切的描述:=EY=h(Z)(2170)如前,z是由自变量x产生的一个向量,是参数,而h就是一个充分光滑的已知函数。如果这一点也做不到,建模就无法进行了。其次,对方差与的关系有一个描述:2=Var(Y)=()(2171)这通常比确定(2170)要难。如在前面“超散布性”这种情况,相依性和非齐次性并不影响均值。所以,如果可能破坏指数性的原因

8、只在这些,它将不影响均值,而(2170)仍可按指数型分布去建立,但方差则不然。由时,有理由认为对不同的x值,“超散布性”只是使方差增加一个与x无关的倍数(已知),则可定()=(按指数型分布所定的方差),1已知(2172)例如对二项分布B(m,)有()=E(29s()9)=ni=1D2i()22i()zizi(2176)在方差设定正确时,有Fn()=COV(s()。可以证明:在一定条件下,当n 时,拟似然方程(2174)以任意接近于1的概率有一解n为0的相合估计,且nN()0,F21nVnF21n)或 V-12nFn(n-0)dN(0,1)(2177)此处Fn,Vn分别是Fn(0)及Vn(0)。

9、解n不一定唯一(实际上,即使方差设定正确,但h1=g非自然联系函数时,解也不一定唯一)。为用于统计推断,必须对(2177)中的Fn,Vn作估计。对Fn的估计用 Fn=Fn(n)(2178)对Vn的估计则不能用Vnn,因Vn()的表达式中包含真方差2i0(),而2i0()并不知道,可以用75广义线性模型(六)QMLE:Quasi MaximumLIkelihood Estimate.Vn=ni=1D2i(n)Yi-h(zin)24inzizi(2.79)估计Vn。在一定条件下可以证明:(2177)的后一极限关系可用 V-12n Fn(n-0)dN(0,1)(2180)所代替,而(2180)可用于

10、统计推断。至于假设检验,Wald检验与以前无异:对原假设C=a用检验统计量(Cn-a)1(C F21n Vn F21nC)21(Cn-a)当它取大值时否定原假设。可以证明:在一定条件下(既:使(2180)成立的那些条件),当原假设成立时,此统计量依分布收敛于x2(r)(C为rp行满秩)。至于其他两个检验,因涉及似然函数,而此处似然函数未知,情况就有所不同。在可能有“超散布性”存在的情况下,方差函数有形状Var(y=V0()(2181)而0()是在无超散分布(即=1)时正确的方差函数。这时有估计 的问题可用估计量n=1m-pmj=1nj(?yj-j)20(j)式中m=样本中x取不同值的个数(如在

11、例111,m=7)把这些值记为x(x),x(m)nj=样本中x=x(f)的个数(如在例111,若x(1)=(110),则n1=28+30=58)。?Yj=上述nj个样本中Y值的平均(如在例111,若x(1)=(110),则?y1=28/58)。j=x=x(f)时,=E(Y)的估计(如在例111,若x(1)=(110),则按前面的表,在logit模型下1=0148),即j=h(Z(j)n)。(Z(j)=Z(x(j)。0(j)=将j代入方差函数0()的结果。如在Y取0,1两值的情况,0(j)=j(1-j)。可以证明:在一定的条件(使极限定理成立)下,且当n增加时m保持有界,nj 对一切j,则n是的

12、相合估计。例111(续)对前章(一)的例111,采用自然联系,利用下表数据,剖腹产事先计划临时决定感 染有 无感 染有 无用抗生素有危险因子没 有11702118700不 用有危险因子没 有283083223309得到回归系数的估计为(此处用(j)记(j)的估计)(0)=21189,(1)=1107,(2)=2103,(3)=23125于是得到估计log感染概率不感染概率=21.89+1.07x1+2.03x2-3.25x385中文核心期刊 数理统计与管理 22卷 4期 2003年7月(回忆:临时决定x1=1;有危险因子x2=1;服用抗生素x3=1)暂把感染概率/不感染概率称为“危险比”,则由

13、上式危险比=e21.89e1.07x1e2.03x2e23.25x3可知最有利的组合是x1=x2=0,x3=1,它的危险比,比之“最不利组合”x1=x1=1,x3=0要小e6132倍,或572倍。有危险因子者其危险比增大(较之无危险因子但其他因素相同者)e2103=716倍。关系最大的是是否服用抗生素,服用者,其危险比缩小e3.2526倍。而临时仓促决定剖腹者,其危险比增大e1.073倍。如果采用Probit模型,即感染概率=(0+1x1+2x2+3x3),N(0,1)则将得0,3的估计分别为?0=21109,?1=0161,?2=1120,?3=21190而得 危险比=1-=(21.09+0

14、.61x1+1.2x2-1.9x3)12(1.09+0.61x1+1.2x2-1.9x3)现把由上述两种模型算出的危险比列表如下:(0 0 0)(0 0 1)(0 1 0)(1 0 0)(0 1 1)1 0 1(1 1 0)(1 1 1)Logit0115110100591115030144040104460601713135350.1300Probit011515010014111731014612010381010087312416011351可以看出,虽则回归系数估计相差甚多,但危险比的估计却很接近,说明采用哪一种模型不很重要。其理由在前面(111(三)结尾)已指出过了(表头上,例如(0

15、 0 1)表示x1=0,x2=x3=1)COV(n)用A21n去估计,An见(2127),算出n各分量的方差及其t值为:(l)(l)(V ar(l)1/2l值i=021189014124161i=1110701432149i=2210301464141i=323125014826167所有的t值都显著。下表给出由logit模型给出的感染概率估值,及由数据得出的经验值:(0 0 0)(0 0 1)(0 1 0)(1 0 0)(0 1 1)1 0 1(1 1 0)(1 1 1)Logit0100018801200111010001480106经验01300177011301110101015301

16、04例如,对x1=1,x2=x3=0一栏,总共有8+32=40个样本,其中受感染者8人,频率(即经验概率)为8/40=012095广义线性模型(六)考虑到样本总量251不算很大,上表的符合程度似乎还过得去。但(0 0 0)一栏很差,究竟如何,需要做拟合优度检验,见第3章。例116(续)使用logit模型,即logP(型感染)P(不感染)=10+11x1+12x2+13x3logP(型感染)P(不感染)=20+21x1+22x2+23x3使用例中的数据,估计结果为:10=221621;11=11174;12=11829;13=23152020=221560;21=01996;22=21195;2

17、3=23.087按这个估计结果,在服用抗生素降低危险比的作用上,型比 型大,分别为e3152=3318倍和e31087=2119倍。而对危险因子在升高危险比的作用上来说,则 型比 型大,分别为e21195=9倍和e11829=612倍。例117(续)x=(x1,x2,x3)都是哑变量:x1=1或21,视年龄 40或否x2=1或21,视从不吸烟或否(2183)x3=1或21,视“以前吸烟,现在不吸”或否z(x)=(x1,x2,x3,x1x2,x1x3)在z(x)中有x1x2及x1x3的项,表示考虑了“年龄”与“吸烟状态”之间的交互作用。使用所给的数据,得到门限值1,2以及x1,x1x3的系数(1

18、),(5)的估计如下表积累logistic分组Cox极大值分布121370(010)01872(010)21429(010)231844(010)11377(010)31843(010)(1)(x1)01114(0129)01068(0104)01095(0137)(2)(x2)01905(010)01318(010)01866(0119)(3)(x3)201364(0101)201110(0102)201359(0114)(4)(x1x2)201557(010)201211(0100)201529(0119)(5)(x1x3)01015(0191)01004(0192)01021(0114)

19、括号内的值是所谓“P值”。它的意义如下:例如,在logistic模型下(1)的P值为0129,表示即使(1)真值为0,那么纯粹由于随机原因,其“估计值的绝对值大于01114”这个事件的概率,也有0129。因此,01114这个值不足以否定0(1)=0这个假设,或更确切地说,现在数据没有支持“在模型中必须包含x1这项”的说法。当然,这不等于说证明了模型中确实不包含这项。这需要有实际背景的参考,或更多的数据。从表面上看,在3个模型之下,同一系数的估值相差甚多。其原因一定程度上可归于111(三)段中所说。重要的是由不同模型算出的概率差异如何(差异必然存在,因这3个模型彼此间并非就只差一个线性变换)。值

20、得注意的是各系数估计值的符号都相同。这表明:虽然06中文核心期刊 数理统计与管理 22卷 4期 2003年7月 使用的模型不同,其对各因素的作用的方向的估计是相同的(见下文的注)。注:在本例中,Y取3值:1(正常),2(边缘),3(不正常)按积累logistic模型有P(Yr)=exp(r+Z(x)1+exp(r+Z(x),r=1,2拿(1)这个分量来说,因为x1取 1两值,而+1相应于年龄 0,则对年龄 0。由上式知P(Yr)之值升高,即情况对 40一组有利。以为P(Yr)的升高表明Y取小值(Y取小值是好的情况)的机会大。相反,对 40那组而言,x1=21,x1(1)A(x)01218201

21、1450106701575201711201548概率 P(Y=rx)=P(Yrx)2=P(Yr-1 x)的计算 1+A(x)1115301727019391144701161013241+A(x)116581123211444119520166601829r=1019680187401922019860170101749r=2010270109401064010130115601150r=301005010320101401001011430110126中文核心期刊 数理统计与管理 22卷 4期 2003年7月 P(Y=1|x)=1-exp(2exp(1+A(x)P(Y=2|x)=1-exp

22、(2exp(2+A(x)-(1-exp(2exp(1+A(x)P(Y=3|x)=exp(2exp(2+A(x)从以上结果可以看出以下几点结论:一是|(1)|很小,显示年龄这个因素对呼吸情况好坏影响不大。二是(2)0且较大,显示吸烟对呼吸情况有负面影响,比年龄因素大。三是不论在年轻一组和年老一组中,“过去吸烟”比之“现在吸烟”,其呼吸情况更坏一些。这可能是由于,这中间一部分人是由于呼吸出了问题才戒烟的。至于计算的过程,对迭代计算求估计值,令=(1,2,),(x)=10z(x)01z(x)=(1,2)=(1,2)=g()=(g1(),g2()=loglog11-1,loglog11-1-2)(=(

23、x)按(2172)式,先要算出(2166)式右边,这可按(2163)式算wj()。在本例有()=COV(Y)=1(1-1)2 122 12 2(1-2)9g()9 =1(1-1)log11-101(1-1-2)log11-1-21(1-1-2)log11-1-2再算Hi9g()9 ()9g()9 211=1-exp(2exp(1+Z(xi)1i2=exp(2exp(1+Z(xi)-exp(2exp(2+Z(xi)2i(2185)样本(Yi,xi),i=1,n.Yi=(Y(1)i,(Y(2)i)Y(1)i=1若第i样本呼吸正常,否则为0Y(2)i=1若第i样本呼吸边缘,否则为0gi=9g()9

24、1、2以(21118)中之值代入s()=ni=1s1(),si()=(xi)HigiY(1)i-1iY(2)i-2i(1i,2i以(21185)中之值代入)迭代公式(k+1)=(k)+(ni=1u(x1)Hiu(x1)21s(k)。初始值由ni=1g(yi)-u(x1)(0)2=minni=1g(yi)-u(x1)2确定。(2186)36广义线性模型(六)在x所能取的不同值的个数不多时,一般在样本(x1,Y1)中会有一些样本有同一值。设全部样本分成m堆,每堆依次有n1,nm个样本,第j堆x值为xj,(x1,xm两两不同),Y值平均为?Yj,且?Yj各分量之和小于1。则g(?Yj)有意义,而(2

25、186)应以mj=1njg(?Yi)-u(x1)(0)2=minnj=1g(?Yi)-u(x1)2(2187)所取代。若?Yj各分量和为1,则(2187)中的?Yj改为nj+1n+2?Yj。这当然也适合于nj=1的情况。COV()的估计值由(ni=1u(xi)Hiu(xi)-1确定。Hi=在Hi中以1,2,分别取代1,2,。参考文献1L.Fahrmeir.Multivariate Statistical Modeling Based on Generalized Linear ModelsM.New Y ork,Springer2Ver2lag,1994.2McCullagh.Generali

26、zed Linear ModelsM.London/New Y ork,Chapman&Hill,1989 2ndedition.3L.Fahrmeir.Consistency and asymptotic normality of the maximum likelihood setimator in generalized linear modelsJ.Ann.Statist.,1985,342-368.(上接第54页)参考文献1 何晓群.现代统计分析方法与应用M.北京:中国人民大学出版社.1998.2美菲利普 科勒特等著,郭国庆等译.市场营销管理 亚洲版M.北京:中国人民大学出版社.1998.3 湖南统计局.湖南统计年鉴M.北京:中国统计出版社.1999.4 郭志刚.社会统计分析方法 SPSS软件应用M.北京:中国人民大学出版社.1999.编者按:将以下图片补在本刊上期第40页图白处。这是第3期“住户抵押贷款的比例提前偿付模型”一文的附图。在此,为印刷错误的出现,特向作者和广大读者致歉。46中文核心期刊 数理统计与管理 22卷 4期 2003年7月