文章编号:1002-1566(2003)04-0055-10
陈希孺
(中国科学研究生院，北京100039)
摘 要：本讲座是广义线性模型这个题目的一个比较系统的介绍。主要分3部分：建模、统计分析
与模型选择和诊断。写作时依据的主要参考资料是L.Fahrmeir等人的:《Multivariate S tatistical M odel2
ing Based on G eneralized Linear M odles》。
关键词：广义线性模型；建模；统计分析；模型选择和诊断
中图分类号:O212 文献标识码:A
G eneralized Linear Modles
CHE N X i2ru
(G raduate School of Chinese Academia of Science,Beijing100039,China)
Abstract:This set of articles gives an introduction to generalized linear m odels.They can be divided into three parts:M odel building,S tatistical in ference and M odel diagnostics.The presentation in mainly based onL.Fahrmeir et al.《Multivariate S tatisticar M odeling Based on G eneralized Linear M odles》.
K ey w ords:generalized linear m odels;m odel building;statistical in ference;m odel diagnostics
213 拟似然法
到此为止所有的讨论都是在Y服从指数型分布的假定下进行的。这个假定的根据是，我们的目的在于离散数据统计分析，而在一些应用上很重要的情况下，这种数据的分布是二项分布、多项分布、P oiss on分布等，它们都属于指数型。 但是，在有些情况下,“指数型”这个假定不一定切合实际：当建模时，往往着眼在变量的均值、方差。像二项、多项、P oiss on这些分布，都是单参型分布，其均值方差依赖于一个参数。因此方差σ2是均值μ的函数:σ2=ν(μ),ν(μ)称为方差函数。例如对Poisson分布,ν(μ)=μ。对二项分布B(m,θ),μ=mθ,σ2=mθ(1-θ)=μ(1-μ/m)。对负二项分布NB(r,θ)①,σ2=r+3μ+μ2/r等等。但实际数据有时不 符合这个关系，如以前提过的所谓“超散布性”(见ξ111末尾)：如对二项分布,Ey=μ，但σ2可以大于μ1-μ/m(若Y～B(m,θ)，则DY=mθ>μ，而Var(Y)= mθ(1-θ)=μ(1-μ/m))。这时，可以证明，不存在一个取0,1,…,m为值(每个值的概率>0), 服从指数分布，而对Φ≠1满足Var(Y)=ΦE(Y)(1-1
m
EY)的变量。
Y的分布称为负二项分布NB(r,θ)。
①每次试验成功的概率为θ，失败的概率为1-θ。给定自然数r,Y=第r次成功时已失败的次数，则
P(Y=y)=r-1+y
r-1
θr(1-θ)y,y=0,1,…
注：设有指数分布族P (Y =i )=c (θ)d (i )e i
θ,i =0,1,…m ，则c (θ)>0对一切θ>0，有c (θ
)=(∑m
i =0d (i )e i θ)21
E (Y )=c (θ
)=(∑m
i =0id (i )e i θ E (Y 2
)=c (θ)=(∑m
i =0
i 2d (i )e i
θVar (Y )=c (θ)=(∑m i =0
i 2d (i )e i θ-c 2(θ
)
(∑m
i =0
id (i )e i θ)2设Var (Y )=<E (Y )(1-E (Y )/m )，则
∑m
i =0
d (i )
e i
θ・∑m
i =0
i 2d (i )e i
θ-(∑m
i =0
id (i )e i
θ)=<∑m
i =0weipu
d (i )
e i
θ・∑m
i =0
id (i )e i
θ-1
m
(∑m
i =0
id (i )e i θ)
2
比较两边e θ
的系数，有d ((0)d (1)=<d (0)d (1)，得出<=1(因Y 取每个值的概率>0，故d (i )>0,i =0,1,…m 。因此d (0)d (1)≠0而必须有<=1)
对P oiss on 分布也有这个情况。当“超散布性”出现时，样本的均值方差不一致。这时就不能以P oiss on 分布为模型。
以上讲的是这样一种情况：原来样本可以认为服从某种指数型分布，由于相依性及非齐次性(指在同一x 之下多次观察的因变量Y 并非同分布，因为还有重要因素未纳入x 内。这些因素每次观察时取值可以不同(非齐性)，因而使Y 值有不同的分布)使指数假定不能成立。另有一些情况，一开始就没有充分理由取指数型分布作为模型。这就说明：在实际问题中往往有必要在对变量的分布并未确切的情况下去建模，并发展出相应的统计推断方法。拟似然法就是为了这个目的。
对我们此处的问题而言，拟似然法着眼在均值和方差，尤其是前者，即必须对均值有一个比较确切的描述:
μ=EY =h (Z ′
β)(2170)
如前,z 是由自变量x 产生的一个向量,β是参数，而h 就是一个充分光滑的已知函数。如果这一点也做不到，建模就无法进行了。
其次，对方差与μ的关系有一个描述:
σ2=Var (Y )=ν(μ
)(2171)
这通常比确定(2170)要难。如在前面“超散布性”这种情况，相依性和非齐次性并不影响均值。
所以，如果可能破坏指数性的原因只在这些，它将不影响均值，而(2170)仍可按指数型分布去建立，但方差则不然。由时，有理由认为对不同的x 值,“超散布性”只是使方差增加一个与x
无关的倍数<(已知)，则可定ν(μ)=<・
(按指数型分布所定的方差),<>1已知(2172)例如对二项分布B (m ,θ
)有ν(μ)=<・μ(1-μ/m )(2173)
有时，与问题实际背景有关的考虑可以给出一个适当的ν(μ)。如果这些条件都不存在，但有充分
多的样本，也可以利用。具体说，在n 个不同的x 值x 1,…x n 处有若干个Y 观察值：在x 1处有Y n ,…Y m i ，算出
Y i =
∑
n i
j -1Y y
/m i ,s 2
i =1
n 1-
1∑n i
j -
1
(Y
y
- Y i )2
,i =1,…,n
在平面坐标上标出点( Y i ,s 2
i ),i =1,…,n 。根据其走势，配以适当的函数作为ν(μ)。这个散点图也可用来检验某种模型是否合适。如在上图，点
大致落在45°角线附近，说明用P oiss on 分布作为模型基本合适，若是下图的情况则不行。
好在有一点，用拟似然法建模并不绝对依靠方差函数规定的正确性。即使所设定的方差函数与实际有出入，仍可对未知参数(2170)中的β作
估计并证明(在一定条件下)估计量的相合性与渐进正态性。要不是这样，这个方法的可用性就大大降低了。当然，这不是说我们不须努力使规定
的方差函数尽量接近于实际。可以证明：接近的程度愈好，估计的效率也就愈高。
下面转到估计问题，它是仿照对数似然方程(2112);
夜尿停∑n
i =1
D i
(β)σ22i
(β)(y 1-μi (β))z i =0(2174)
在此,z i 已有,μi (β)由(2170)给出，为h (z i ’β)。D 1(β)由(2111)给出,σ2i (β)由(2171)给出:σ2
i (β)=ν(μ1)=ν
(h (z i ’β))。有了这些，可由方程(2174)求解得出β0
的(拟极大似然)估计值^βn (QML E )①。
仍把(2174)的左边记为s (β)，并记σ2
i 0(β)=Var (Y 1)。由于EY i =μi (β
)(因2170)为正确设定)，有
COV β(s (β))=∑n i =1D 2
i (β)σ2
i 0(β)σ4i (β)z i z i ’ƒV n (β
)(2175)如果方差函数设定正确:σ2i =σ2i 0，则(2175)转化为(2115)(改β0为β。(2115)是在真参数β0韩素音
处计算的)。又记
F n (β)>E β(29s (β)9β)=∑n
i =1D 2i (β)σ22
i (β)z i z i ’(2176)在方差设定正确时，有F n (β
)=COV β(s (β))。可以证明：在一定条件下，当n →∞时，拟似然方程(2174)以任意接近于1的概率有一解
^βn 为β0
的相合估计，且
^βn ～N (β)0,F 21n V n F 21n ) 或 V -1
2n F n (^βn -β0
)→d
N (0,1)
(2177)
此处F n ,V n 分别是F n (β0)及V n (β0)。解^βn 不一定唯一(实际上，即使方差设定正确，但h 1整合营销论文
=
g 非自然联系函数时，解也不一定唯一)。
为用于统计推断，必须对(2177)中的F n ,V n 作估计。对F n 的估计用
^F n =F n (^βn )
(2178)
对V n 的估计则不能用V n ^βn ，因V n (β)的表达式中包含真方差σ2i 0(β)，而σ2
i 0(β
)并不知道，可以用
①QM LE:Quasi M aximum LIkelihood Estimate.
^V n =
∑
n
i =1
D 2
i (^βn )Y i -h (z i ’^βn ))2
σ4
i ^βn
z i z i ’
(2.79)
估计V n 。在一定条件下可以证明:(2177)的后一极限关系可用
^
V -1
2n ^
F n (^βn -β0
)→d
N (0,1)
(2180)
所代替，而(2180)可用于统计推断。
至于假设检验,Wald 检验与以前无异：对原假设C
β=a 用检验统计量(C ^βn -a )1(C^F 21n ^V n ^F 21n C ’)21
(C ^βn -a )
当它取大值时否定原假设。可以证明：在一定条件下(既：使(2180)成立的那些条件)，当原假
设成立时，此统计量依分布收敛于x 2(r )(C 为r ×p 行满秩)。至于其他两个检验，因涉及似然函数，而此处似然函数未知，情况就有所不同。在可能有“超散布性”存在的情况下，方差函数有形状
Var (y =ΦV 0(μ
)(2181)而ν0(μ
)是在无超散分布(即Φ=1)时正确的方差函数。这时有估计Φ的问题可用估计量^Φn =
1
m -p
∑
m
j =1
n j ( y j -^μj )
2
ν0(^μ
j )
式中
m =样本中x 取不同值的个数(如在例111,m =7)
把这些值记为x (x ),…x (m )
n j =样本中x =x (f )的个数(如在例111，若x (1)
=(1 1 0)，则n 1=28+30=58)。
Y j =上述n j 个样本中Y 值的平均(如在例111，若x (1)=(1 1 0)，则 y 1=28/58)。^μj =x =x (f )时,μ=E (Y )的估计(如在例111，若x (1)=(1 1 0)，则按前面的表，在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 保持有界,n j →∞对一切j ，则^Φn 是Φ的相合估计。
例111(续) 对前章(一)的例111，采用自然联系，利用下表数据,
剖腹产事先计划
聚类分析案例
临时决定
感 染有 无
感 染有 无用抗生素有危险因子
没 有
1 170 211 870 0不 用有危险因子
没 有
28 308 32
23 30 9
得到回归系数的估计为(此处用^β(j )记β(j )的估计)
^β(0)=21189, ^β(1)=1107, ^β(2)=2103, ^β(3)=23125于是得到估计log 感染概率不感染概率
=21.89+1.07x 1+2.03x 2-3.25x 3
(回忆：临时决定x1=1；有危险因子x2=1；服用抗生素x3=1)
暂把感染概率/不感染概率称为“危险比”，则由上式
危险比=e21.89・e1.07x1・e2.03x2・e23.25x3
可知最有利的组合是x1=x2=0,x3=1，它的危险比，比之“最不利组合”x1=x1=1,x3=0要小e6132倍，或572倍。有危险因子者其危险比增大(较之无危险因子但其他因素相同者)e2103= 716倍。关系最大的是是否服用抗生素，服用者，其危险比缩小e3.25≈26倍。而临时仓促决定剖腹者，其危险比增大e1.07≈3倍。
如果采用Probit模型，即
感染概率=Φ(β0+β1x1+β2x2+β3x3),Φ～N(0,1)
则将得β0,…,β3的估计分别为
β0=21109, β1=0161, β2=1120, β3=21190
而得 危险比=
Φ
1-Φ
=
Φ(21.09+0.61x
1
+1.2x2-1.9x3) 12Φ(1.09+0.61x1+1.2x2-1.9x3)
现把由上述两种模型算出的危险比列表如下:
(000)(001)(010)(100)(011)101(110)(111) Logit0115110100591115030144040104460601713135350.1300
Probit011515010014111731014612010381010087312416011351
可以看出，虽则回归系数估计相差甚多，但危险比的估计却很接近，说明采用哪一种模型不很重要。其理由在前面(ξ111(三)结尾)已指出过了(表头上，例如(001)表示x1=0,x2=x3 =1)
COV(^βn)用A21n去估计,^A n见(2127)，算出^βn各分量的方差及其t值为:
β
(l)
^β(l)(V^a r(^β(l)))1/2l值
i=021189014124161
i=1110701432149
高苯乙烯橡胶i=2210301464141
i=323125014826167
所有的t值都显著。
下表给出由logit模型给出的感染概率估值，及由数据得出的经验值:
(000)(001)(010)(100)(011)101(110)(111) Logit0100018801200111010001480106
经验0130017701130111010101530104
例如，对x1=1,x2=x3=0一栏，总共有8+32=40个样本，其中受感染者8人，频率(即经验概率)为8/40=0120
考虑到样本总量251不算很大，上表的符合程度似乎还过得去。但(000)一栏很差，究竟如何，需要做拟合优度检验，见第3章。
例116(续) 使用logit模型，即
log P (Ⅰ型感染)
P(不感染)
=β10+β11x1+β12x2+β13x3
log P (Ⅱ型感染)
P(不感染)
=β20+β21x1+β22x2+β23x3
使用例中的数据，估计结果为:
β
10=221621; β11=11174; β12=11829; β13=231520
β
20
=221560; β21=01996; β22=21195; β23=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),…β(5)的估计如下表
积累logistic分组C ox极大值分布
θ
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)
括号内的值是所谓“P值”。它的意义如下：例如，在logistic模型下β(1)的P值为0129，表示即使β(1)真值为0，那么纯粹由于随机原因，其“估计值的绝对值大于01114”这个事件的概率，也有0129。因此,01114这个值不足以否定β0(1)=0这个假设，或更确切地说，现在数据没有支持“在模型中必须包含x1这项”的说法。当然，这不等于说证明了模型中确实不包含这项。这需要有实际背景的参考，或更多的数据。
从表面上看，在3个模型之下，同一系数的估值相差甚多。其原因一定程度上可归于ξ111(三)段中所说。重要的是由不同模型算出的概率差异如何(差异必然存在，因这3个模型彼此间并非就只差一个线性变换)。值得注意的是各系数估计值的符号都相同。这表明：虽然
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议