2022年2022年广义线性模型_十_ .pdf

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

1、文章编号:10021566(2004)02007308广义线性模型(十)陈希孺(中国科学院研究生院,北京 100039)摘 要:本讲座是广义线性模型这个题目的一个比较系统的介绍。主要分3部分:建模、统计分析与模 型 选 择 和 诊 断。写 作 时 依 据 的 主 要 参 考 资 料 是L.Fahrmeir等 人 的MultivariateStatisticalModeling Basedon GeneralizedLinear Models。关键词:广义线性模型;建模;统计分析;模型选择和诊断中图分类号:O212文献标识码:AGeneralized linear modelsCHEN Xi2r

2、u(Graduate Schoolof ChineseAcademia of Science,Beijing 100039,China)Abstract:This setof articles gives an introduction to generalizedlinear models.They can be divided into three parts:Modelbuilding,Statistical inference and Model diagnostics.The presentation is mainly based on L.Fahrmeir et al.Multi

3、variateStatistical Modeling Based on GeneralizedLinear Models.Key words:generalizedlinear models;model building;statistical inference;model diagnostics3.3 诊断问题这是一个内涵不很确定,方法也不很规范的问题。一般在讲到这个题目时多要涉及以下几方面的内容:1.强影响点,即样本中对决定模型有较强影响的那种点。这种点之值得注意,是因为万一它们有较大的误差,将使模型产生较大的偏差。2.残差分析,残差是指因变量观测值Yi与其由模型估计的“理论值”Yi的

4、差。残差的整体状况对研判模型的正确性即数据中有无异常有用,这后一方面就是。3.异常值检测。所谓“异常值”是指那种样本点,它因受到某种系统性因素的影响(如记录错误)而带有很大的误差,表现为其取值远离其它样本之值。这种值起的作用总是负面的,异常值检测的目的就是发现这种值并将其从样本中剔除。即使在最简单的正态线性回归场合,也不能说对以上几个问题有了理论上坚实且(尤其重要的)在应用上行之有效的整套方法,更遑论远为复杂的广义线性模型了。因此,以下所讲的只宜被看作一些原则性的指点,其成功的使用在很大程度上依赖于经验以及对所论问题的专业知识和实际背景的了解。(一)帽子矩阵与高杠杆点先从最简单的线性回归来讲,

5、比较容易理解。设有模型Yi=0+xi0+ei,i=1,n(3.14)37广义线性模型(十)?1994-2007 China Academic Journal Electronic Publishing House.All rights reserved.http:/名师资料总结-精品资料欢迎下载-名师精心整理-第 1 页，共 8 页 -记 0=00,X3=1x11xn,Y=Y1Yn,e=e1en可将(3.14)写为 Y=X30+e(3.15)按 LS 估计,得 0的估计=S3-1X3Y,S3=X3X3于是由(3.15),得Y依模型得到的估计,即其“理论值”为Y=X3S3-1X3Y=PxY(3.

6、16)Px为X3的列向量空间的投影阵。因为它是由Y产生“Y加帽子”的矩阵,故常称为“帽子矩阵”并记为H:H=PX(3.17)假定(3.14)中的ei满足 G auss-Markov条件:Ee=0,Cov(e)=2In(In:n阶单位阵)(3.18)则残差 Y-Y=(I-H)Y(3.19)的期望为 0,因E=(I-H)EY=(I-H)X30=0(后一式由于X3=HX3)而协方差阵为COV()=2(I-H)(3.20)这里用到I-H为幂等阵。特别,记=(1,n),i=Yi-Yi,记H=(hij),有Var(i)=2(1-hii),i=1,n(3.21)(因H为投影阵,有 0 hii1,故 0 Va

7、r(i)2)从(3.21)式看出:若hii1,则Var(i)=Var(Yi-Yi)0。结合E(Yi-Yi)=0,有yi-yi0。这意思是:在(xi,yi)这个点,经验值yi与其理论值yi拟合得特别好,或者说,这样的点有把经验回归平面(即由数据估计得的回归平面)拖向自己的倾向。它被称为高杠杆点,名词的直观理由见下面的解释。往证(记号意义见下)hii=(1,xi)S3-11xi=1n+(xi-x)S-1(xi-x)(3.22)为此,记 1=(1,1,1),X=x1xn,x=1nn1xi=1X利用分块矩阵求逆公式:ABB C-1=A-1+A-1BGBA-1-A-1BG-G BA-1 G(3.23)其

8、中A,C为对称方阵,而G=(C-BA-1B)-1。利用此公式,注意到S3=n1XX 1XX(3.24)以及 XX-X11Xn ni=1(xi-x)(xi-x)(3.25)有 S3-1=1n+xS-1x-xS-1-S-1xS-1(3.26)因此 hii=1n+xSx-2xiS-1x+xiS-1xi=1n+(xi-x)S-1(xi-x),即(3.22)x是数据x1,xn的中心。(3.22)显示:某点xi距此中心越远,则hii之值愈高 如前所述,47中文核心期刊 数理统计与管理23卷 2期 2004年 3 月?1994-2007 China Academic Journal Electronic P

9、ublishing House.All rights reserved.http:/名师资料总结-精品资料欢迎下载-名师精心整理-第 2 页，共 8 页 -它愈有把经验回归拉向自己的倾向。这一点在x为 1 维时,不难从图形上得到理解:看图 3.1,xi远离x,其hii接近 1,是一个高杠杆点。如果(xi,yi),在 点处(属正常状态),则尚不致太影响经验回归线的基本走向。反之,若(xi,yi)在B处,则它把经验回归线由l至l(拖向自己),而大为影响了经验回归线的走向。看图上,好像有一股力把l的右端抬上去,好像以(x,y)为心的一个杠杆作用。这正是高杠杆点这一名称的由来。反观图 3.2,xi离x

10、近,非高杠杆点。这时,如果(xi,yi)处在离开样本群体较远的位置B,其作用大略是把l平行地带上去一点点而方向变化不大,意味着这不大会影响对回归系数的估计而只是影响了其常数项。从实用上看,回归系数比模型中的常数项更重要,因它反映了所考虑的因素对目标变量的影响,而常数项则只反映原点位置,因此,高杠杆点的影响值得注意。这一点到后面再谈。对高杠杆点要避免一种误解,即认为它是不好的,情况不然。一则样本xi离中心x有远有近,在样本量大时,高杠杆点的存在无可避免。二则xi离中心远,除非yi异乎寻常地“出格”很多,它也不足以使经验回归线转动一个太大的角度。三则如果一切样本xi都云集在其中心x附近,情况反而不

11、好,因为在 G M 模型下有 COV()=2S3-1由此式及(3.26)看出:COV()=2S-1,而S=ni=1(xi-x)(xi-x)。当xi都与x接近时,S会偏小,而S-1偏大,即COV()偏大,这意味着的精度低。高杠杆点的意义只在于,在这个异常的y值会对回归系数的估值造成重大影响,因而值得注意。hii要多大才能算是高杠杆点?这并无公认的规定。根据H为幂等,有ni=1hii=tr(H)=tr(X3S3-1X3)=tr(S3-1X3X3),(p为的维数)=tr(S3-1S3)=tr(Ip+1)=p+1这里假定了n(p+1)矩阵X3的秩为p+1(否则S3-1不存在,H也就只能通过广义逆表示:

12、H=X3S3-1X3,这时仍可证明H为幂等,但其 tr 即其秩是等于X3的秩 2(p+1)/n时认为此点为高杠杆点。如果n相对于p很大,则甚至接近于0 的hii也可能被认为是高杠杆点,不甚合理。另一种判法根据如下的结果:若(3.14)中的(x1,Y1),(xn,Yn)为 iid 多维正态,则57广义线性模型(十)?1994-2007 China Academic Journal Electronic Publishing House.All rights reserved.http:/名师资料总结-精品资料欢迎下载-名师精心整理-第 3 页，共 8 页 -Fn-p-1phii-1/n1-hii

13、F(p,n-p-1)(3.27)选定 :(0,1)(如=0.05,0.01),当FF(p,n-p-1)时,认为F值因而hii值异常之高,即当 hii(n-p-1p+pF(p,n-p-1)/(n-p-1+pF(p,n-p-1)(3.28)时,认该点为高杠杆点。对广义线性模型,则是移植以上的想法。当自变量取离散值时,在x的一个值xj处有若干个(nj个)Y值Yj1,Yjnj,起作用的是这些值的平均Yj而不是单个的Yjk。这是因为Y有指数型分布,因而似然含数gj=1exp(jnjk=1Yjk)(略去与 无关的因子)只依赖于Yj,j=1,g,故的 MLE 也只依赖于Yj,所以,个别Yjk其的作用已综合在

14、Yj中,不必一一去考察,只需看Yj就行。按模型,得 jEy-jh(zi0)的估计为j=h(zin),均方误差为E(Y-j-j)2=E(Y-j-j)-(j-j)2=E(Y-j-j)2+E(j-j)2-2E(Y-j-j)(j-j)J1+J2-2J3(3.29)有 J1=2j/nj,(j2=Var(yjk)(3.30)当n较大时,n0,故近似的有j-j=h(zjn)-h(zj0)Djzj(n-0)Djzj-1nS(0)=Djzj-1ngt=1ztDt-2tnt(y-t-t)(3.31)这里用到了(2.21)(给出 n-0-1nS(0),n见(2.15)以及(2.12)(给出S(0),这里利用了 nt

15、k=1(Ykt-t)=nt(Y-t-t),又Dt=dh(u)duu=zt0)注意到当jt时,Y-j与Y-t独立因而E(Y-t-t)(Y-j-j)=0 由(3.31)有J3Djzj-1nzjDj(3.32)J2Djzj-1n(gt=1ztDtnt-2tztDt)-1nzjDj=Djzj-1nn-1nzjDj=Djzj-1nzjDj(3.33)由(3.32),(3.33),结合(3.29)和(3.30),得E(Y-j-j)21nj2j(1-njD2j-2jzj-1nzj)(3.34)令(g为x取不同值的个数)z=z1zg,D=D2100D2g,E=21002g,N=n100ng而作方阵 K=ND1

16、/26-1/2z-1nz 6-1/2D1/2N(3.35)则易见 njD2j-2jzj-1nzj=K的(j,j)元(3.36)而K为幂等阵,此因 K2=ND1/261/2zn-1z 61/2D1/2N?ND1/261/2z-?-1nz 61/2D1/2N67中文核心期刊 数理统计与管理23卷 2期 2004年 3 月?1994-2007 China Academic Journal Electronic Publishing House.All rights reserved.http:/名师资料总结-精品资料欢迎下载-名师精心整理-第 4 页，共 8 页 -=gt=1ntD2t-2tztzt

17、=n=ND1/261/2z-1nz 61/2D1/2N=K(3.37)在线性回归(且 G M 条件满足),有E(Yi-Y-i)=2(1-hii),而hii为一幂等阵H的(i,i)元。由(3.36)知,广义线性模型之下的公式(3.34)形式上与之完全相似,而我们可以把njD2j-2jzj-1nzj 1(3.38)的点(y-j,xj)称为高杠杆点(Dj,2j,n与0有关,可用n代0)不同的是:1.(3.36)左边不只依赖zj(zj相当于线性模型中的(1,xj),还依赖于未知的参数真值0,须用其估计值n代替;2.公式(3.34)只是近似而非严格成立;3.(3.34)右边还有因子2j,n-1j。因此,

18、除了高杠杆点(3.38)成立之点)有“把经验回归拖向自己”的效应外,2j很小或nj很大之点,也有这种效应,即后两种情况也是值得重视之点,或者看其综合效应:使(3.34)右边很小之点,由把经验回归拖向自己的效应。以上讲的是Y为一维的情况,若Y为q维,则仍以Y-j=(Y-j1,Y-jq)计自变量x取值xj时Y y的所有观察值的算术平均。zj=z(xj)有形状zj=zj1zjq每个zjk为p维行向量,p:的维数替代一维的E(Y-j-j)2,取COV(Y-j-j),此处 j=EY-j=(j1,jq),它依赖于参数真值0,而j=(j1,jq)则是以 n取代0的结果。有COV(Y-j-j)=COV(Y-j

19、-j)+COV(j-j)-E(Y-j-j)(j-j)-E(j-j)(Y-j-j)使用关系式(2.51),以及n-0-1nS(0),仿照(3.29)(3.34)的算法,易得COV(Y-j-j)=1/2j(Iq-nj-1/2jDjzj-1nzjDj-1/2j)1/2jnj(3.39)此处Dj见(2.41):Dj=9h(t)9tt=zj0(h=g-1,g为联系函数,6j=COV(Yjk)令z=z1zg,D=D100Dg,6=61006g,N=n100ng而作方阵 K=N6-1/2Dz-1nzD 6-1/2NK为gq阶方阵,将其分块为g2个q阶方阵,则主对角线第j块即为nj6-1/2jDjzj-1nz

20、jDj6-1/2j,且K为幂等阵。这个形式统一了q=1和q 1两种情况。由(3.39),仿照q=1时(3.34)的情况,可以把满足nj6-1/2jDjzj-1nzjDj6-1/2jIq(3.40)的点称为高杠杆点。由于nj6-1/2jDjzj-1/2nzjDj6-1/2jIq,也可以用tr(nj6-1/2jDjZj-1nzjDj6-1/2j)q(3.41)来取代(3.40)。与q=1 的情况一样,一个点对经验回归影响的大小,还与 6j与nj有关。77广义线性模型(十)?1994-2007 China Academic Journal Electronic Publishing House.Al

21、l rights reserved.http:/名师资料总结-精品资料欢迎下载-名师精心整理-第 5 页，共 8 页 -(二)强影响点一个样本点(yi,xi)对回归系数0的估计有多大的影响?考察此问题的一个方法是:把这样本点剔除,用剩下的样本点,对0作一估计,记为(i),与原来(使用全部样本)的估计 比较。若二者差距不大,则认为此样本点影响不大,否则就比较大。对简单线性模型(3.14),容易证明-(i)=(X(i)3X3(i)-1x3(i)Yi-S3-1x3(i)Yi/(1-hii)(3.42)这里X3(i)是在X3(见(3.14)中删去第i个行向量x3(i)所得的矩阵,Yi是(3.16)的Y

22、的第i分量,hii见(3.22)。在误差e1,en为 iid,且有公共分布N(0,2),则有-(i)N(0,2M),M是某个只与x1,xn有关的矩阵。故12(-(i)M-1(-(i)2(p+1)。用S2=1n-p-1Y-Y2估计 2,于是当n充分大时,有(-(i)M-1(-(i)/S2 2(p+1)(近似)(3.43)(p是的维数,(,)=的维数为p+1)注意此处我们不能断言(3.43)左边有F分布F(p+1,n-p-1),虽则S2和独立,但与(i)并非独立。有了(3.43),我们可以计算任一个样本点的“p值”。设(3.43)左边之值记为,则p(2(p+1)称为该样本点的p值。p值愈小,表示(

23、i)离愈远,而(yi,xi)对决定回归系数估计的影响愈强。若取=0.1,则大致可以说,样本点中大约只有10%,其值(i)能达到与如此距离的位置。可以约定取一个适当的,0 c)=。(三)残差样本Yi与其从模型中估得的“理论值”之差,称为(该样本点的)残差。所谓“理论值”,即Eyii的估计i。在简单线性模型且误差e1,en为 iid,N(0,1)的场合,残差的确切分布不难求得。对广义线性模型则只能满足于近似结果。残差分析的目的,一是根据全部样本的残差的整体表现,去研判模型是否存在问题和存在那种问题。一般的说,当模型没有严重的系统性错误(如,假定回归函数为线性,其实是二次的或周期的;假定误差为等方差

24、,其实误差方差是随x的值增长而增长等)时,残差在整体表现上应给人一种杂乱无章的印象,从中看不出什么明显的倾向性。从这也可领悟到,残差分析是一种 informal 的工具。专家告诫要避免“overinterpletation”(“过度解释”,“求甚解”,陶渊明:“好读书,不求甚解”),此语大体上可理解为:只有那种一眼望去就能感觉到的(残差图上的)倾向性,才是比较可靠的。如残差图成明显的弯曲形态判断线性假设可能有问题,残差沿某一方向系统的增大表明等方差假设可能有问题等。这都无明确的规律可循,使用者以自己对问题背景的了解和经验,结合某些统计知识,去进行研判。它说到底只是一种辅助性的工具,其作用是有限

25、的,不能对之期望值太高。另外,残差有多种类型(上面提到的y-y只是其一,且不是97广义线性模型(十)?1994-2007 China Academic Journal Electronic Publishing House.All rights reserved.http:/名师资料总结-精品资料欢迎下载-名师精心整理-第 7 页，共 8 页 -最常用的一种),有时需要进行适当的变换后才计算残差,这些问题的恰恰相反当处理无不依赖使用者的经验(以及有关的专业知识,当然也包括有关的统计理论和方法)。不列颠百科全书称统计学是一门“scienceandart”。Art 这个词用在残差分析上是再贴切不过

26、了。残差分析的另一个目的是根据个别点的残差去研判该点是否为异常点。这通常是指那种点,其残差(绝对值)与其它点的残差比较起来是突出的大。这个问题在上面(二)中提到过了。在广义线性模型中,最常用的残差有两种:1.学生化残差,又称 Pears on 残差 rip=(Yi-i)/(i)其中 i=EYi=h(zi0),i=h(zn),zi=z(xi),n为0的 MLE,2()为方差函数:EY=,Var(Y)=2()。这是一种中心化和标准化:(Yi-i)/(i)有均值 0 方差 1,而 i用其估计值i取代。这样做是为了有一个统一的比较标准,如果单看Yi-i,则这个值大不一定意味着(yi,xi)这个点可能有

27、问题,而有可能是由于Yi的方差大。如果在同一个xi处有多个(ni个)y值,则通常考察其平均Y-i,而定义 rip=ni(Yi-i)(i)当ni很大时,如果模型没有错误,将有ripN(0,1)。可在这一背景下去研判残差图的表现。当ni不很大时,rip不仅与正态有距离,且其分布也往往是偏斜(非对称)的。有的作者设计了一些变换以纠正之,就是以yi的某个函数取代yi,例如,在y为两点分布(=P(y=1)=1-P(y=0)时取 rip=nit(y-i)-t(i)+i(1-i)-1/3(2i-1)/6nii(1-i)1/6在y为 P oisson 分布时取 rip=32niy-i-(2/3i-1/3i)/

28、91/6i2.Deviance残差把deviance定义(3.19)和中每一项的平方根定义为该点的残差:riD=sgn(y-i-(i)2(li(i,i)-li(i,(i)当y-ii时取正值或0,当y-ii取负值。按导致(3.22)的推理,可知当ni甚大时,基本上有rip=riD。参考文献1L.Fahrmeir.Multivariate Statistical Modeling Basedon G eneralized Linear Models M.New York,Springer-Verlag,1994.2 McCullagh.Generalized Linear ModelsM.Lond

29、on/New York,Chapman&Hill,1989 2ndedition.3 L.Fahrmeir.Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear modelsJ.Ann.S tatist.,1985,342368.08中文核心期刊 数理统计与管理23卷 2期 2004年 3 月?1994-2007 China Academic Journal Electronic Publishing House.All rights reserved.http:/名师资料总结-精品资料欢迎下载-名师精心整理-第 8 页，共 8 页 -