《陈希孺-广义线性模型_五_.pdf》由会员分享，可在线阅读，更多相关《陈希孺-广义线性模型_五_.pdf（8页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、文章编号:10021566(2003)03005608广义线性模型(五)陈希孺(中国科学院研究生院,北京 100039)摘 要:本讲座是广义线性模型这个题目的一个比较系统的介绍。主要分3部分:建模、统计分析与模型选择和诊断。写作时依据的主要参考资料是L.Fahrmeir等人的 Multivariate StatisticalModeling Based on Generalized Linear Models。关键词:广义线性模型;建模;统计分析;模型选择和诊断中图分类号:O212文献标识码:AGeneralized Linear ModelsCHEN Xi2ru(Graduate Schoo
2、l 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:Model building,Statistical inference and Model diagnostics.The presentation is mainly based on L.Fahrmeir et al.Multivariate Sta
3、tistical Modeling Based on Generalized Linear Models.Key Words:generalized linear models.;model building;statistical inference;model diagnostics第二部分 统计推断2.1极大似然估计(MLE)(一)单参数情况设有独立样本(Xi,Yi),i=1,n。Yi的分布写成指数标准形c(yi)exp(yii-b(i)d(yi),i=1,n(2.1)i与Xi有关,也与参数有关。似然函数为L=ni=1c(yi)exp(yii-b(i)要找=n,使L达到最大。这就是的ML
4、E。取对数logL=6ni=1logc(yi)+6ni=1(yii-b(i)其中c(yi)不依赖于,对估计无影响,故可略去。记li(i)=yii-b(i),i=1,n(2.2)称 l()=6ni=1li(i()=6ni=1(yii()-b(i()(2.3)为对数似然函数。有9l()9=6ni=1 yi-?b(i()9 i()9(?b()=db()d)(2.4)似然方程 9l()9=6ni=1 yi-?b(i()9 i()9=0(2.5)65中文核心期刊 数理统计与管理 22卷 3期 2003年5月 1994-2008 China Academic Journal Electronic Publ
5、ishing House.All rights reserved.http:/的解就是的MLEn,因为 i=EYt=?b(i()(2.6)且为突出i与的关系,常写为i(),可以把(2.5)写为6ni=1(yi-i()9 i()/9=0(2.8)至于 9 i/9,那就取决于联系函数g的形式。以h记g的反函数,有i=h(zi),这里zi=z(0)(xi)是由xi决定的一个向量,也可以就是xi。有9 i9=didi9 i9,9 i9=9h(zi)9zizi(2.9)9h(zi)9zi=dh(t)dtt=zi注意到i=?b(i),有di/di=b(i)=V ar(yi)=2i()(2.10)由(2.9
6、),(2.10)得9 i/9=Di()-2i()zi,Di()=dh(t)/dtt=zi(2.11)代入(2.8),得似然方程为(i()=h(zi)s()=6si()=0,si()=Di()-2i()(yi-i()zi(2.12)在自然联系的情况,i=zi,有 9 i/9=zi。于是由(2.8)直接得到似然方程为6ni=1si()=0,si()=(yi-i()zi(2.13)形式大为简化(在自然联系的情况,h=?b,故Di()=b(zi)=b(i)=2i(),而(2,12)式中的si()为(yi-i)()zi)MLE的极限分布(启发式推导,非严格证明)以0记的真值,以与流动值相区别。则 Jn=
7、6ni=1si(0)(2.14)为独立之和,有期望0,协方差阵为(COV在0之下计算)n(0)=n=COV(Jn)=6ni=1COV(si(0)=6D2i(0)-2i(0)zizi(2.15)于是按中心极限定理,在一定的条件下,当n充分大时,有 JnN(0,n)(2.16)按n的定义,有 6ni=1si(n)=0,故令ci()=Di()-2i(),有Jn=Jn-6ni=1si(n)=6ni=1(si(0)-si(n)=6ni=1zi ci(0)(yi-i(0)-ci(n)(yi-i(n)=6ni=1zici(0)(i(n)-i(0)+6zi(ci(0)-ci(n)(yi-i(n)(2.17)在
8、一定的条件下n是0的相合估计,故n0,当n充分大。因为ci()对连续(这要求函数h连续可微),故ci(0)-ci(n)为无穷小。因此,把(2.17)右边最后一项与Jn=6ni=1zici(0)(yi-i(0)比较,知它相对于Jn为无穷小。由(2.17)知,Jn的极限分布应与75广义线性模型(五)1994-2008 China Academic Journal Electronic Publishing House.All rights reserved.http:/?Jn=6ni=1zici(0)(i(n)-i(0)(2.18)的极限分布相同。因n0,近似地有i(n)-i(0)9 i(0)9
9、(n-0)=Di(0)zi(n-0)(2.19)以此代入(2.18)有?Jn=6ci(0)Di(0)zizi(n-0)=n(n-0)(2.20)由此可知,当n充分大时(n-0)的分布近似于-1n?Jn的分布近似于-1nJn的分布近似于-1nN(0,n)的分布近似于多元正态N,(0,n)(2.21)或者说 1/2n(n-0)dN(0,Ip)(2.22)Ip为p阶单位矩阵p为0的维数。严格的证明即是加上一些条件,把上面的heuristic(启发式的)步骤变严格,得分两步走:第一步要证明n为0的相合估计,因为上面在作近似Taylor展开时(2.19)式)用到这一点。这里就存在问题:似然方程(2.12
10、)是否有解?确可举出解不存在的例子:设Yi取0,1两值,i=P(Yi=1),为1维。用自然联系函数,有 i=xi1+exi方程(2.13)成为 6ni=1xi(yi-xi1+exi)=0(2.23)设xi 0,i=1,n 又Yi=1,i=1,n则上式成为6ni=1xiexi1+exi=0(2.24)对任何,上式左边总 0,故方程(2.23)无解。即使在方程有解的场合,还有一个解的唯一性的问题 解可以不唯一。不过,在一定的条件下(在应用上这些条件都满足),以下两点成立)1.limnP(方程(2.12)有解)=1(2.25)2.有一个解n0,.s.(2.26)如果是自然联系函数,则解必唯一。这时,
11、(2.25)、(2.26)可总和成一条:“似然方程以概率1有唯一根n0,.s.”n的估计为了利用(2.21)或(2.22)对0作统计推断,需要知道 n(见(2.15)式)。n的公式中包含真参数0,可以用n替代它,得到 n估计量n=6ni=1D2i(n)-2i(n)zizi(2.27)考虑到n0及Di,i的连续性,在一定的条件下易证n-n=nos(1)os(1)是一个p阶方阵,.s.0。这样,在(2.21)和(2.22)中,把n改为n仍成立。而这可用于统计推断。记 G()=-6ni=19si()/9 =-9l()/9 9 (2.28)往证 E(G(0)=n(2.29)85中文核心期刊 数理统计与
12、管理 22卷 3期 2003年5月 1994-2008 China Academic Journal Electronic Publishing House.All rights reserved.http:/因为E(Yi-i()=0,由(2.12)式看出,在计算 9si()/9 时所得的3项中,只有一项,即ziDi()-2i()(-9 i()/9 )这一项,其期望值不为0,这一项为Di()zi。故E(G(0)=6ni=1D2i(0)-2i(0)zizi=n(2.30)按这个结果,也可以用G(n)去估计n。在自然联系的情况,h=?b,而Di()=b(i)=2i(),故n=6ni=12i(0)z
13、izi(2.31)形式得到简化。如在Yi取0,1两值的情形:n=6ni=1i(1-i)zizi i=P(Yi=1)=ezi0/(1+ezi0)MLEn的迭代计算设从某一初始值出发,已经k步算到(k),用牛顿法求下一步的(k+1):令(k+1)=(k)+,H()=6ni=1D2i()-2i()zizi0=s(k+1)=s(k)+)s(k)+9s(k)/9 ,以2H(k)近似取代 9s(k)/9 ,得 H-1(k)s(k)(2.32)而把(k+1)=(k)+H-1(k)s(k)(2.33)作为k+1步迭代值。用一种修正的牛顿法,即用(k+1)=(k)+依次取=1,0.9,0.5,(12)2,(12
14、)3,直到一个(k+1)满足 s(k+1)s(k)为止。停止迭代的条件可取为(k+1)-(k)/(k)0是一个设定的数。初始条件0可取为(g(yi),zi),1in下线性回归系数的LS估计,即0=(6n1zizi)-16n1ziyi(2.35)引进记号?yi()=zi+D-1()(yi-i(),i=1,n(2.36)i()=D2i()-2i(),i=1,n(2.37)有(k+1)=(6ni=1wi(k)zizi)-16ni=1wi(k)yi(k)zi(2.38)这是一个加权最小二乘估计的形式,即min6ni=1(yi(k)-zi)2wi(k)(2.39)的解。这个形式使我们可以在计算(k)时,
15、使用在线性回归中已有的程序。(2.38)易证。事实上,根据(2.36),(2.37)及(2.33),有(2.38)右边=(6ni=1wi(k)zizi)-16ni=1wi(k)yi(k)zizi(k)+(6ni=1wi(k)zizi)-16ni=1wi(k)D-1i(k)(yi-i(k)zi=(k)+H-1(k)s(k)=(k+1)95广义线性模型(五)1994-2008 China Academic Journal Electronic Publishing House.All rights reserved.http:/(二)多参数情况在原则上与一维情况完全相似。仿前面的记法,有li(i)
16、=yii-b(i),(,y:q维;:p维)l()=6ni=1li(i()=6ni=1(yi i()-b(i()似然方程为9l()9 =6ni=1 yi-?b(i()9 i()9 =0(?b()=9b()9)(2.40)以g记联系函数,g-1=h,则(g,h:q1)i=E(yi)=h(Zi),(Zi:qp,由xi决定,Zi=Z(xi)注意到 9 i/9 i=9b()/9 9 =b()=V ar(yi)=6i(),有9 i/9 =6i()9 i/9 9 i/9 =6i()-19 i/9 =6i()-1DiZi(2.41)这里Di()=9 i/9tt=Zi=9h(t)/9tt=Zi=9h(Zi)/9
17、(Zi),(Di()=9h(Zi)/9(Zi),代入(2.40),有9l()9 =6ni=1ZiDi()6ni=1()(yi-i()=0,(i()=h(Zi)(2.42)上式可写为s()=9l()9=6ni=1si()=6ni=1Ziwi()9g(i)9 (yi-i()=0(2.43)si()=ZiDi()6-1i()(yi-i()=Ziwi()9g(i)9 (yi-i()(2.44)wi()=Di()6-1i()Di()=(9g(i)9 6i()9g(i)9)-1(2.45)此处利用了 g(h(t)=t9g(h(t)9h(t)9h(t)9t=Iq(令t=Zi 并注意h(Zi)=i)g(i)9
18、 Di()=Iqg(i)9 =(Di()-1(2.46)有(COV表示协方差在参数为时计算)n()=COV(s()=6ni=1COV(si()=6ni=1ZiDi()6-1i()6i()Di()Zi=6ni=1Ziwi()Zi(2.47)=ZW()Z(2.48)Z=Z1ZnW()=wi()00wn()(2.49)若记 Y=(Y1Yn)Y=Y1Yn,()=(1()n()06中文核心期刊 数理统计与管理 22卷 3期 2003年5月 1994-2008 China Academic Journal Electronic Publishing House.All rights reserved.ht
19、tp:/D()=D1()00Dn(),6()=61()006n()(2.50)可将对数似然方程写为 s()=ZD()6-1()(Y-()=0(2.51)在一定条件下可以证明:1.P(s()=0有解)1,当n。2.存在s()=0之一解n,是0的相合估计。3.1/2n(0)(n-0)dN(0,Ip),或(ZW(0)Z)1/2(n-0)dN(0,Ig),或以n代替0。4.(ZW(n)Z)1/2(n-0)dN(0,Ip)(2.52)n-0N(0,(ZW(n)Z)-1/2)(2.53)迭代求n的方法也与1维相似:(k+1)=k+(ZW(k)Z)-1s(k)(2.54)s见(2.69)初始值(0)可以由(
20、g(yi),Zi),i=1,n.用最小二乘法:min6ni=1g(yi)-Zi2,解为(0)(2.55)2.2 假设检验考虑y为1维的情况,多维的情况相似。在广义线性模型,实用上最常见的假设检验问题还是线性假设(0为参数真值)H0:C0=H1:C0(2.56)一个重要的特例是的一个或几个分量为0。例如 0(1)=0 0(1)0(2.57)在(2.56)中,0为p维,C是已知的rp常数矩阵,为行满秩(因而rp)1.Wald检验.检验统计量为n=(Cn-)(C-1nC)-1(Cn-)(2.58)这里n为0的MLE,n为COV(s(0)的估计,见(2.27)式。该式中的Di由(2.11)式给出 (2
21、.11)式中的h是连续函数g的反函数。若H0成立,则=C0,令n=(1/2n(n-0)(C-1/2n)(C-1/2n)(C-1/2n)-1(C-1/2n)(1/2n)(n-0)=n Bnn(2.59)按(2.22)(在有关假定(使(2.22)成立)之下),有ndN(0,Ip),又注意到Bnr为投影阵,因而为幂矩阵,秩为r,因此有 nd2(r)(2.60)注意到n=(C(n-0)(C-1nC)-1C(n-0)=(Cn-)(C-1nC)-1(Cn-)与,(2.58)比较,看出它与n的区别,只在于用 n取代了n,但n为n的估计(n是在n=ni=1D2i(0)-2i(0)zizi中,用n取代0)。故(
22、2.60)用n取代n仍成立。于是有nd2(r)(当原假设H0成立时)(2.61)16广义线性模型(五)1994-2008 China Academic Journal Electronic Publishing House.All rights reserved.http:/利用(2.61)可得到H0的检验:n2(r)时,否定H0(2.62)其渐进水平为。其所以在n取大值时否定H0,是因为当H0不成立时,E(Cn-)C0-0,n倾向于取大值。2.约束检验 利用约束下的MLE以?n记在原假设C=这个约束下0的MLE。以un=s(?n)-1n(?n)s(?n)?n:在约束下的MLE(2.63)作为
23、检验统计量。当u 某个常数时,否定原假设。此检验的直观背景如下:因s(n)=0,若原假设成立,则?n和n为同一0的估计,理应比较接近;因而s(?n)s(n)=0,这时u将取小值。反之,若原假设不成立,则?n与n不接近,s(?n)与0有距离,而u将取较大之值。可以证明:当原假设成立,且一定的条件满足时,有 und2(r),n(2.64)于是上文提到的常数可取为2(r),(0,1)为给定的水平。此处r的意义同前:C为rp的行满秩矩阵,1rp。3.拟似然比检验仍以ln()记对数似然函数。n和?n分别为真参数0的不受任何约束的MLE以及受到原假设约束的MLE。检验统计量为n=2(ln(n)-ln(?n
24、)(2.65)因ln(n)为ln()的最大值,总有n0。若原假设成立,则n和?n都是0的相合估计,n和?n理应接近而n倾向于小。反之若0不满足原假设,则n和?n有较大差距而n倾向于大。因此取 n const(2.66)为检验的否定域。可以证明:当原假设成立时,有 nd2(r)(2.67)这里r的意义同前,依(2.67),(2.66)中的const可取为2(r)以上3个检验出发点不同,但其(在原假设成立时的)极限分布一致。这有一个解释。设对数似然函数ln()为的2次函数。因其在MLEn处达到最大值,必有形状ln()=const-12(-n)A(-n)(2.68)对某个A 0,为求在约束C=之下的
25、极值点?n,用L agrange乘数法,作ln()+(C-),对求导命之为0,得C=A(?n-n)?n=n+A-1C Cn+CA-1C =(CA-1C)-1(-Cn)?n=n+A-1C(CA-1C)-1(-Cn)于是得到拟似然比统计量为212(?n-n)A(?n-n)=(Cn-)(CA-1C)-1(Cn-)2.69 由(2.68)有-52ln()/5 5=A,故A=n。因A为已知的常数方阵,n=A。把(2.69)右边与(2.58)比较,知拟似然比统计量与Wald统计量一致。至于约束检验统计量,按(2.68),s(?n)=A(n-?n),n(?n)=n=A,故(2.63)式的un为(n-?n)A
26、-1(n-?n),与(2.69)比较,知它与拟似然比统计量和Wald统计量一致。当原假设成立时,?n与n很接近。因此在计算拟似然比统计量及约束检验统计量时,只26中文核心期刊 数理统计与管理 22卷 3期 2003年5月 1994-2008 China Academic Journal Electronic Publishing House.All rights reserved.http:/涉及对数似然函数ln()在n很小邻域内之值,而在这样的邻域内,ln()可以很好的用一个二次型去逼近之。这解释了为何二者有统一的极限分布(在原假设成立下)。参考文献1L.Fahrmeir.Multivari
27、ate Statistical Modeling Based on Generalized Linear Models M.New York,Springer-Verlag,1994.2McCullagh.Generalized Linear ModelsM.London/New York,Chapman&Hill,1989 2ndedition.3L.Fahrmeir.Consistency and asymptotic normality of the maximum likelihood estimator in generalized linearmodelsJ.Ann.Statist
28、.,1985,342-368.上接第46页参考文献1 高旅端,陈志.不完全数据下Weibull分布的参数估计J.数理统计与应用概率,1995,10(3):832892 王松桂,程维虎,高旅端.概率论与数理统计M.北京:科学出版社,2000.3DEMPSTER A P,LAIRD N M,RUBIN D B.Maximum likelighood from incomplete data via the EM algo2rithmJ.J R Statist Soc,1977,B39:1238.4 王松桂.EM算法J.应用数学与计算数学,1983(6):43248.36广义线性模型(五)1994-2008 China Academic Journal Electronic Publishing House.All rights reserved.http:/