Logistic回归模型

Logistic 回归模型
1 Logistic 回归模型的基本知识 1.1 Logistic 模型简介
主要应用在研究某些现象发生的概率p ，比如股票涨还是跌，公司成功或失败的概率，以及讨论概率
p 与那些因素有关。

显然作为概率值，一定有10≤≤p ，因此很难用线性模型描述概率p 与自变量的关
系，另外如果p 接近两个极端值，此时一般方法难以较好地反映p 的微小变化。

为此在构建p 与自变量关系的模型时，变换一下思路，不直接研究p ，而是研究p 的一个严格单调函数)(p G ，并要求)(p G 在p 接近两端值时对其微小变化很敏感。

于是Logit 变换被提出来：
p
p
p Logit -=1ln
)( （1）
其中当p 从10→时，)(p Logit 从+∞→∞-，这个变化范围在模型数据处理上带来很大的方便，
解决了上述面临的难题。

另外从函数的变形可得如下等价的公式：
X
T X
T T e
e p X
p
p
p Logit ββ
β+=
⇒=-=11ln )( （2）
模型(2)的基本要求是，因变量（y ）是个二元变量，仅取0或1两个值，而因变量取1的概率)
|1(X y P =就是模型要研究的对象。

而T
k x x x X ),,,,1(21 =，其中i x 表示影响y 的第i 个因素，它可以是定性变量也可以是定量变量，T
k ),,,(10ββββ =。

为此模型(2)可以表述成：
k
x k x k x k x k
k e
e
p x x p
p βββββββββ+++++++=
⇒+++=- 11011011011ln （3）
显然p y E =)(，故上述模型表明)
(1)
(ln
y E y E -是k x x x ,,,21 的线性函数。

此时我们称满足上面条件
的回归方程为Logistic 线性回归。

Logistic 线性回归的主要问题是不能用普通的回归方式来分析模型，一方面离散变量的误差形式服从伯努利分布而非正态分布，即没有正态性假设前提；二是二值变量方差不是常数，有异方差性。

不同于多元线性回归的最小二乘估计法则(残差平方和最小)，Logistic 变换的非线性特征采用极大似然估计的方法寻求最佳的回归系数。

因此评价模型的拟合度的标准变为似然值而非离差平方和。

定义1 称事件发生与不发生的概率比为 优势比(比数比 odds ratio 简称OR)，形式上表示为
OR=
k
x k x e p
p βββ+++=- 1101 （4） 定义2 Logistic 回归模型是通过极大似然估计法得到的，故模型好坏的评价准则有似然值来表征，称
-2ˆln ()L β
为估计值βˆ的拟合似然度，该值越小越好，如果模型完全拟合，则似然值ˆ()L β为1，而拟合似然度达到最小，值为0。

其中ˆ()lnL β
表示βˆ的对数似然函数值。

定义3 记)ˆ(β
Var 为估计值βˆ的方差-协方差矩阵，2
1
)]ˆ([)ˆ(ββVar S =为βˆ的标准差矩阵，则称 k i S w ii
i i ,,2,1,]ˆ[
2 ==β （5）
为i
βˆ的Wald 统计量，在大样本时，i w 近似服从)1(2
χ分布，通过它实现对系数的显著性检验。

定义4 假定方程中只有常数项0β，即各变量的系数均为0，此时称
20
ˆˆ2[ln ()ln ()]L L χββ=-- （6） 为方程的显著性似然统计量，在大样本时，2
χ近似服从)(2
k χ分布。

1.2 Logistic 模型的分类及主要问题
根据研究设计的不同，Logistic 回归通常分为成组资料的非条件Logistic 回归和配对资料的条件Logistic 回归两种大类。

还兼具两分类和多分类之分，分组与未分组之分，有序与无序变量之分。

具体如下： 两分类非条件Logistic 回归：分组数据的Logistic 回归，未分组数据的Logistic 回归； 多分类非条件Logistic 回归：无序变量Logistic 回归，无序变量Logistic 回归； 条件Logistic 回归：1:1型、1:M 型和M:N 型Logistic 回归。

关于Logistic 回归，主要研究的内容包括：
1． 模型参数的估计及检验 2． 变量模型化及自变量的选择 3． 模型评价和预测问题 4． 模型应用
2 Logistic 模型的参数估计及算法实现
2.1 两分类分组数据非条件Logistic 回归
因变量(反应变量)分为两类，取值有两种，设事件发生记为y=1，不发生记为 y=0，设自变量
T k x x x X ),,,(21 =是分组数据，取有限的几个值；研究事件发生的概率)|1(X y P =与自变量X 的关
系，其Logistic 回归方程为：
k k x x X y P X y P βββ+++=== 110)|0()|1(ln 或 k
x k x k
x
k x e
e X y P ββββββ+++++++== 1101101)|1( 例2.1.1 分组数据[1]
在一次住房展销会上，与房地产商签订初步购房意向书的有n=325人，在随后的3个月时间内，只有一部分顾客购买了房屋。

购买房屋的顾客记为1，否则记为0。

以顾客的年家庭收入(万元)作为自变量X ,对数据统计后如表2.1.1所示，建立Logistic 回归模型。

表2.1.1 购房分组数据
例2.1.2 药物疗效数据[2]
为考察某药物疗效，随机抽取220例病人并分配到治疗组和对照组，治疗组采用治疗药物，对照组采用安慰剂。

治疗一段时间后观察病人的疗效，得到表2.1.2数据。

设y 为疗效指标(y=1 有效，y=0无效)，1x 为治疗组指标(1为治疗组，0为对照组)，2x 为年龄组指标(1为>45岁，0为其他)。

上述两个例子数据都是经过统计加工后的分组数据，对此类数据进行Logistic 回归，首先要明确应变量对应事件的发生概率如何确定和进行Logit 变换，其次才能建立Logistic 回归。

为便于数据处理，我们将此类数据的格式作个约定，排列格式为（组序号，自变量X ，该组事件发生数，该组总例数）。

表2.1.3 分组数据的标准格式
表2.1.1 改造表
表2.1.2 改造表
经过改造后，可得我们关心的事件的发生的频率为 n i n m p i ,,2,1,i
i
==
该组总例数该组发生事件数。

其中n
为分组数，然后作Logit 变换，即i
i
i i p p p Logit p -==1ln )(~。

变换后的数据，形式上已经可以采用一般
的线性回归的处理方式来估计回归参数了。

此时方程变为：∑==+=k
j ij j i n i x p 1
0,,2,1,~
ββ 当然这样处理并没有解决异方差性，当i n 较大时，i p ~
的近似方差为： )(,)
1(1
)~(i i i i i i y E n p D =-≈
πππ （7）
所以选择权重 n i p p n i i i i ,,2,1),1( =-=ω，最后采用加权最小二乘法估计参数。

注意，分组数据的Logistic 回归只适用于大样本分组数据，对小样本的为分组数据不适用，并且以组数n 为回归拟合的样本量，明显降低了拟合精度，在实际应用中必须谨慎。

求解算法及步骤：
1．依据分组数据的标准格式，计算频率i p 、Logit 变换i p ~
和权重i ω 2．构建加权最小二乘估计：
∑∑∑∑====--=--n i k
j ij j i i i i n i k j ij j i i x y x y 1
1
201
1
2
0)(min )(min βωβωωββω （8）
令 i i i y y ω=
*
，T ik i i i i i x x X ),,,(1*ωωω =，T k ),,,(10ββββ =
则方程又变成一般的线性回归模型：∑=-n
i i T i X y
1
2**
)(min
β （9）
3．构造增广矩阵21***
*][+⨯+k k T T
Y X X X
利用消去法得]ˆ)ˆ([ββ
Var I =矩阵，得到估计βˆ
其中2,1++K K I 为残差平方和SE , 回归方差1
ˆ2
--=k n SE
σ
各系数检验采用 )1(~ˆˆ--=
k n t I t ii i i σ
β
总平方和∑∑∑===-
=
n
i n
i i
n
i i
i i
i y y
ST 1
1
2
12
2)
()
(ωωω，回归平方和SE ST SR -=
总平方和求解相当于拟合i i y ωβ*
0*=方程的残差平方和，故得上式ST
所以方程的检验为)1,(~)
1/(/----=
k n k F k n SE k
SR F
例2.1.1的求解过程如下（由LLLStat 统计软件计算）：
表2.1.4 数据Logit 变换及权重
家庭年收入x 实际购买mi 签订意向ni
比例pi 逻辑变换Logit 权重ni*pi(1-pi) 1.500000 8 25 0.320000 -0.753772 5.440000 2.500000 13 32 0.406250
-0.379490
7.718750 3.500000 26 58 0.448276 -0.207639 14.344828 4.500000 22 52 0.423077 -0.310155 12.692308 5.500000 20 43 0.465116 -0.139762 10.697674 6.500000 22 39 0.564103 0.257829 9.589744 7.500000 16 28 0.571429 0.287682 6.857143 8.500000 12 21 0.571429 0.287682 5.142857 9.500000
10
15
0.666667 0.693147
3.333333
表2.1.5 回归模型基本信息 总样本 9
求解方法 加权最小二乘 仅常数项beta0 -0.095029 方程F 统计量 51.982160 F 分布自由度 1,7 方程检验p 值 0.000176 总平方和 8.798294 回归平方和 7.754112 残差平方和
1.044181
表2.1.6 分组Logistic 回归系数检验
序号 均值
回归系数
系数标准误 t 统计量 自由度df
检验P 值 常数项 2.837815 -0.848882 0.113578 -7.473994 7 0.000056 家庭年收入x
14.901140 0.149323
0.020711
7.209865
7
0.000056
表2.1.7 1
][-X X T
0.086479 -0.014517 -0.014517 0.002876
本例Logistic 模型的回归方程：x
e x
e p
i 149323.0848882.0149323.0848882.01ˆ+-+-+=
对于多分类无序自变量的Logistic 回归，即某个自变量为m 个水平的名义变量（如治疗方法A,B,C ），
只需要引入m -1（2个）个哑变量，然后采用上述方法进行分析。

例2.1.3 研究三种治疗方法对不同性别病人的治疗效果[2]
，数据如表2.1.4
表2.1.4 性别和治疗法对某病治愈情况的影响
由于治疗方法有三种，没有等级关系，所以属于无序的名义变量，故引入两个哑变量32,x x 分别代表A 和B 疗法，其中0,132==x x 表示方法A, 1,032==x x 表示方法B, 0,032==x x 表示方法C ，将上述数据转化成标准格式，得表2.1.5。

表2.1.5 性别和治疗法对某病治愈情况的影响
对于分类数据，也可以采用极大似然法进行参数估计，具体见2.2节最后部分内容。

2.2 两分类未分组(连续)非条件Logistic 回归
应变量y 取值为0和1，设事件发生记为y=1，否则为0，设自变量T
k x x x x ),,,(21 =，n 组观测数据记为),,,,(21i ik i i y x x x ,n i ,,2,1 =。

记T
ik i i i x x x X ),,,,1(21 =，10=i x ，则i y 与ik i i x x x ,,,21 的
Logistic 回归模型是：
i
X T i X T ik
x k i x ik x k i x ik k i i i e
e
e
e
x x f y E βββββββββββπ+=
+=
+++==++++++11)()(110110110 （10）
易知，i y 是均值为i π的0-1型分布，其分布律为
i
y i i
y i
i y f --=1)
1()(ππ，n i y i ,,2,1;1,0 ==
则n y y y ,,,21 的似然函数和对数似然函数分别为： ∏=--=
n
i i
y i i
y i
L 1
1)
1(ππ
∑∑==-+-=--+=n
i i i
i
i n i i i i i y y y L 1
1
)]1ln(1ln
[)]1ln()1(ln [ln πππππ 代入ik
x k i x ik x k i x i e
e
ββββββπ+++++++=
1101101,得
∑∑==++++-=+-+++=n
i i
X T i T i n
i ik
x k i x ik k i i e
X y e
x x y L 1
1110110)]
1ln([)]
1ln()([ln ββββββββ （11）
记)(ln )(ββL LL =，选取T k ),,,(10ββββ =的估计T k
)ˆ,,ˆ,ˆ(ˆ10ββββ =使得)(βLL 达到极大，这就是Logistic 回归模型的极大似然估计，该过程的求解需要采用牛顿迭代法。

构造得分函数k g LL F g
g ,,2,1,0,)
()( =∂∂=
βββ，共k +1个非线性方程组，令其=0求解β，其中 k g e
e
x x y F n
i i
X T i X T ig ig i g ,,2,1,0,]1[)(0
=+-
=∑=βββ (12 )
构造信息矩阵k h g h
g LL I gh ,,2,1,0,,)
()(2 =∂∂∂-
=ββ，即)(βLL 二阶导矩阵的负矩阵，其中 k h g e
e x x I n
i i
X T i
X T ih ig gh ,,2,1,0,,)
1()(0
2
=+=∑
=βββ (13 ) 很明显)()(ββhg gh I I =，故)(βI 是一个对称矩阵。

求解算法及步骤：
1． 根据公式(12 ) 计算得分函数)(βg F ，公式（13）计算信息矩阵)(βgh I
给定初值)0,,0,0(0
==0β, k =1 和精度ε，可取0.000001
2． 采用牛顿迭代式 βββ
∆+=-1k k
, )()]([111---=∆k k F I βββ，通过以下方式求解。

构造增广矩阵)(1
-k IF β=))()((11--k k F I ββ，通过对IF 矩阵作k +1次ij 消去变换求解β∆
若εβ
β<∆=
∆∑=k
g g
2|||| 或者 εββ<∆=∆∑=k
g g 0
|||||| 或者 εβ<∆≤≤|}{|max 0g k
g ，则转3
否则k = k +1，继续执行第2步
3． 此时k
β就是回归系数β的数值估计βˆ，k 就是迭代次数，消去变换后的IF 矩阵的前1
1+⨯+k k 子阵就是β方差-协方差矩阵的估计阵11)()ˆ(+⨯+=k k gh V Var β=V ，下面给出检验有关计算：
计算Wald 统计量 gg
g g V W 2ˆβ=
,近似服从)1(2
χ分布，检验p 值 ))1((2g g W P p >=χ
标准误gg g V E S =).(.β, g
g e
OR ββˆ)(=, k g ,,1,0 =
例2.2.1 公共交通调查数据
[1]
在一次关于公共交通的社会调查中，调查项目为“是乘坐公共汽车上下班，
还是骑自行车上下班”。

因变量y=1表示乘坐公共汽车，y=0表示骑自行车。

自变量1x 是年龄，作为连续变量；2x 是月收入(元)；3x 是性别，3x =1表示男性，3x =0表示女性。

调查对象为工薪族群体，数据如表2.2.1所示。

表2.2.1 公共交通社会调查
序号 年龄1x 月收入2x 性别3x 交通 y 1 18 850 0 0 2 21 1200 0 0 3 23 850 0 1 4 23 950 0 1 5 28 1200 0 1 6 31 850 0 0 7 36 1500 0 1 8 42 1000 0 1 9 46 950 0 1 10 48 1200 0 0 11 55 1800 0 1 12 56 2100 0 1 13 58 1800 0 1 14 18 850 1 0 15 20 1000 1 0 16 25 1200 1 0 17 27 1300 1 0 18 28 1500 1 0 19 30 950 1 1 20 32 1000 1 0 21 33 1800 1 0 22
33
1000
1
23 38 1200 1 0 24 41 1500 1 0 25 45 1800 1 1 26 48 1000 1 0 27 52 1500 1 1 28
56
1800
1
1
以下计算结果采用LLLStat 1.0 软件得到：
表2.2.2 主要计算结果
序号 均值 回归系数 系数标准误 wald 统计量 自由度df 检验p 值 OR=Exp(B) 常数项 0.535714 -3.655016 2.091223 3.054766 1 0.080501 0.025861 年龄 1273.214286 0.082168
0.052119 2.485516 1 0.114899 1.085639 月收入 0.464286
0.001517
0.001865 0.661466 1 0.416043 1.001518 性别
36.107143 -2.501844
1.157818
4.669175
1
0.030709 0.081934
表2.2.3 Logistic 模型基本信息
总样本 28
求解方法
极大似然法 & Newton 迭代 迭代次数(仅beta0) 7(4) -2LogLikelihood(Beta) 25.970652 仅常数项beta0 -0.143101 -2LogLikelihood(beta0) 38.673263
方程Wald 值(相减) 12.702611 方程自由度 4 方程检验p 值
0.012824
对于例2.1.3分组数据的极大似然估计法，主要过程如下：
∏=--=n
i i
m i n i i m
i m n i
i
C L 1)1(ππ
∑∑==-+-+=--++=n
i i i i
i
i m n n
i i i i i i m n n m C m n m C L i i
i
i
1
1
)]1ln(1ln [ln )]1ln()(ln [ln ln πππππ代入ik
x k i x ik x k i x i e
e
ββββββπ+++++++=
1101101,得 ∑=+-+=
n
i i
X T i i T i m n e
n X m C L i i 1
)]1ln([ln ln ββ
则有 g g LL F βββ∂∂=
)
()(k g e
e x n x m n
i i X T i
X T ig i ig i ,,2,1,0,]1[1 =+-=∑=ββ
=
∂∂∂-
=h
g LL I gh )
()(2
ββk h g e
e x x n n
i i
X T i
X T ih ig i ,,2,1,0,,)
1(1
2
=+∑
=ββ； 其中i m ,i n 分别表示分组i 中事件发生次数和总观察数，如表2.1.4和2.1.5所示。

然后可采用Newton-Raphson 迭代法进行求解。

由LLLStat 计算得到如下结果。

表2.2.4 性别和疗法对某病治愈的影响(未分组Logistic 似然估计法)
序号 均值 回归系数 系数标准误 wald 统计量 自由度df 检验P 值 常数项 1.000000 1.418399 0.298690 22.550513 1 0.000002 性别 0.500000 -0.961618 0.299797 10.288472 1 0.001339 治疗A 0.333333 0.584745 0.264108 4.901966 1 0.026826 治疗B
0.333333
1.560763
0.315961
24.400993
1
0.000001
表2.2.5回归系数方差矩阵V(beta)(信息矩阵I(Beta)的逆矩阵)
0.089215 -0.072957 -0.029931 -0.030097 -0.072957 0.089878 -0.000078 0.000128 -0.029931 -0.000078
0.069753 0.029993 -0.030097
0.000128
0.029993
0.099831
2.3 条件Logistic 回归[2,3]
条件Logistic 回归是配对设计(病例-对照)中常用的一种统计分析方法，通过配对方法收集资料：每一配对组可包括一个病例和一个或多个对照，有1：1型、1：m 型配对。

假设收集了如下数据：
表2.3.1 n 个1:m 配对组,k 个协变量的比例资料
配对资料用配对的方法来控制影响因素的干扰，并且每个配对组都可以建立一个Logistic 回归方程：
n i x x p Logit k
k i ,,2,1,)(110 =+++=βββ
为此需要估计的参数有n 个常数项n
01
0,,ββ 和k 个回归系数k ββ,,1 ，配对数越多估计的参数就越多，但是一般的数据量难以支撑这样的估计，故一般的Logistic 回归不适合配对资料。

不过在参数估计时，常数项会被消去，所以方程组减少了n 个常数项n
01
0,,ββ 的估计，复杂度大大降低。

对于回归参数的估计采用条件似然函数替代一般的似然函数进行。

对于第i 个配对组而言，共有m +1个观察对象，记为m B B B A ,,,,21 ，其中仅有一例发病，且正好是病例组A 发病，而对照组均没有发病的条件概率i p （类似Bayes 概率）可以表示成：
∑=+=
m
j m j m m i B B B A P B B B A P B B B A P p 112121)
()()
( （14）
其中)(21m B B B A P = )|0()|0()|1(1100m
i m i i i i i X y P X y P X y P === ，而
j
i X T j i X T j i j i e
e
X y P ββ+=
=1)|1(，j
i X T j i j i e
X y P β+=
=11)|0(，m j ,,2,1 = (15 )
故n 个配对组的条件似然函数表示为：
∏∏
∏
∑
∏
∏
∏
=-=-==-==≠===∑+=∑+=++++++++=n
i m
j i X j i X T n
i m
j i X j i X T n
i m
k m
k
j j j
i X T k
i X T k i X T i X T m
j j
i X T i X T i X T m
j j
i X T i X T i X T e
e
e
e e
e
e
e e
e e
e
L 1
1
1
)0(11
)
0(1110
1
01
0]
1[11
11111
111111)(βββββββββββββ (16 )
则对数似然函数)()(ββLnL LL =为：
∑==-∑+-==n
i m
j i X j i X T e
LnL LL 1
1
)
0()1ln()()(βββ （17）
令 )(0i j i j i X X D -=，它是一个与第i 个样本点有关的k 维向量，j
ig D 表示向量中的第g 个元素，
则有如下得分函数和信息矩阵：
g g LL F βββ∂∂=)()(=∑===∑+∑-n i m j j
i
D T m
j j
i D T j
ig
e e D 111
1ββ h
g LL I gh ∂∂∂-=)
()(2ββ
=
k h g e
e
D e D e
e
D D m
j j
i D T m
j j
i D T j ih
m
j j i D T j ig
n
i m
j j
i D T m
j j
i D T j ih
j ig ,,2,1,])
1(1[
2
1
1
1
1
1
1
=∑+∑∑-
∑+∑======∑βββββ
注意此时的T
k ),,,(21ββββ =，没有常数项0β。

至此(17)式中的参数β可采用Newton-Raphson
迭代法求解了，β初值依然取为0向量。

不过该方程的求解已经相对复杂多了。

方程似然度检验和回归系数的wald 检验同非条件Logistic 回归。

例2.3.1 研究肥胖、口服避孕药雌激素与子宫内膜癌的关系，随机抽取20名患者，对于每名患者，
在随机抽取年龄相近的正常人作为对照。

检测患者与正常人的肥胖程度和雌激素服用情况[3]。

例2.3.1 求解的主要结果，由LLLStat软件计算得到：
表2.3.2 条件Logistic回归系数检验
序号均值(病例) 回归系数系数标准误wald统计量自由度df 检验P值肥胖0.650000 1.823914 0.547192 11.110390 1 0.000859 雌激素0.850000 1.589621 0.450544 12.448367 1 0.000419
表2.3.3条件Logistic回归模型基本信息
样本量20
求解方法极大似然+牛顿迭代
迭代次数45
-2LogLikelihood(Beta) 33.306763
-2LogLikelihood(0) 43.944492
方程Wald值(相减) 10.637728
方程自由度 2
方程检验p值0.004898
2.4 多分类有序反应变量Logistic 回归
在实际应用中，经常遇到反应变量为多分类有序变量的情况，例如评价指标分为差、中、良、优等，各等级之间是有序的。

这种资料的Logistic 回归分析通常称为比例比数模型(累积概率模型) [4]
，它需要拟合m -1 (m 为水平或等级个数)个Logistic 回归模型。

有序累积概率Logistic 模型：
1,,2,1;,,2,1,1)|(-==+=
≤++m j n i e
e
X j y P i
X T j a i X T j a i i ββ 或 （18）
1,,2,1,
)
|(1)
|(ln
1
1
-=+==-=∑∑==m j X X k y P X k y
P i T j j k i i j
k i i
βα （19）
有序累积概率模型参数的极大似然估计就是寻找参数使得联合概率实现最大化，由于观测之间相互独立，联合概率被分解成边缘概率之积。

而观测到j y i =的概率就是累积概率之差：
)|1()|()|(i i i i i i X j y P X j y P X j y P -≤-≤==
第i 个观测值对应似然值的贡献取决于观测到哪一个j 值，因此对于次序响应的每个j 值，取所有j y i =的观测之的乘积，有似然函数：
∏∏====n i m
j ij
d i i X j y P L 11
)
|(，其中若j y i =，则1=ij d ，否则0=ij d
并且对于任一个观测i X 而言，只有一个等级事件发生，即
∑===m
j i i
X j y
P 1
1)|(，故有(19)式。

其对
数似然函数如下（对于分组数据，似然函数变为：∏∏====
n
i m
j ij
d i n i i X j y P L 1
1
)
|(，i n 分组中各分类例数）。

]
)11ln(
)11ln(1ln
[)
|(ln ln 1
1
2
111111111∑∑∑∑=-=+-+-+++-+-++==+-
+++-
++===n
i m j i
X T j a i X T j a i
X T j a i X T j a ij i
X T m a i X T m a im i
X T a i X T a i n
i m
j i i ij e
e
e e
d e
e
d e
e
d X j y P d L ββββββββ （20）
其中：m
j m j j e
e e e e e
e e X j y P i X T m a i X T m a i
X T j a i X T j a i X T j a i X T j a i X T a i
X T a =-≤<=⎪⎪⎪⎪⎪⎩⎪⎪⎪
⎪⎪
⎨⎧+-+-
++==+-+-+-+-++++111
11111)|(111111ββββββββ （21）
然后就可以通过极大似然法，就上Newton-Raphson 方法加以求解参数β,,,11-m a a 了，注意的是
121-<<<m a a a 。

下面给出具体推导β,,,11-m a a 求解的详细过程。

对（20）式进行化简，可得
)}
1ln()]1ln([]
)1ln()1ln()ln([{ln 11111
1
2
11
i
X T m a im i
X T a i T
i n i m j i
X T j a i
X T j a j a j
a i T
ij e
d e
X a d e
e
e
e
X d L ββββββ+-+=-=+-+-+-+-+++-+--+=∑∑ （22）
∑=+++++--+=∂∂n
i i X T a i
X T a a a a i i X T a i e
e
e e e d e d a L 111121
2111)]1(11[ln βββ （23） ∑=+-+-+-+------+-+--=∂∂n i i X T m a i
X
T m a
im i X T m a i X
T m a
m a m a m a
im m e e d e e e e e d a L 11
1
1121111]1)1([ln ββββ （24） 2,,2,)]1()1([ln 1111-=++--+--=∂∂∑=++++++-m g e
e
e e e d e e e e e d a L n
i i X T g a i
X T g a g a g a g
a ig i X T g a i
X T g a g a g a g
a ig g ββββ （25） k
g e
e
e e d e e d e d x L n
i m j i X T j a i
X T j a i X T j a i
X T j a ij
i X T m a i
X T m a im i X T a i ig g ,,2,1,])111(111[ln 11
2111111 =+-+-++-+=∂∂∑∑=-=+-+-+++-+-+ββββββββ （26） ∑=+++++++-++-=∂∂∂n
i i X T a i
X T a a a a a i i X T a i
X T a i e e
e e e d e e d a a L 12112122
122111112)])1()(()1([ln ββββ （27） ∑=+-+-+-+----+----++++--=∂∂∂n
i i X T m a i
X T m a im i X T m a i X T m a m a m a m a m a im m m e e
d e e e e e d a a L 12112112212
11112])
1())1()(([ln ββββ （28） 2
,,2)])
1()
(())
1()(([ln 12
2
1
112
211
2
-=++
-+++--=∂∂∂∑=++++++++--+m g e
e
e e
e d e
e
e e e
d a a L n
i i
X T g a i X T g a g
a g a g
a g a ig i
X T g a i X T g a g a g a g a g a ig g g ββββ （29）
2,,2,1,)
(ln 1211
112-=---=∂∂∂∑=+++++m g e e e d a a L n
i g a g a g a
g a ig g g （30） k h m g e
e
d d x a L
n
i i
X T g a i X T g a ig ig ih h g ,,2,1;1,,1,)1()
(ln 1
2
12 =-=++-=∂∂∂∑=+++βββ （31） k
h g e e
e e d e e d e e d x x L n
i m j i X T j a i
X T j a i X T j a i
X T j a ij
i X T m a i
X T m a im i X T a i
X T a i ih ig h g ,,2,1,)])
1()1(()1()1([ln 11
2211221121112
=+++++++-=∂∂∂∑∑=-=+-+-+++-+-++ββββββββββ （32） 由此构建信息矩阵),(βa I 和),(βa F ，并可迭代求解了。

注：若为分组数据，上述每项乘以i n 。

例2.4.1研究性别和两种治疗方法对某种疾病疗效的影响[3],将疗效分成效果显、有效和无效三个等级，根据试验调查，得到如下资料。

表2.4.2 多分类有序反应变量数据格式
行号性别治疗方法频数疗效等级
1 1 1 16 1
2 1 1 5 2
3 1 1 6 3
4 1 0 6 1
5 1 0 7 2
6 1 0 19 3
7 0 1 5 1
8 0 1 2 2
9 0 1 7 3
10 0 0 1 1
11 0 0 0 2
12 0 0 10 3
计算结果，由LLLStat统计软件给出：
表2.4.3 回归系数方差矩阵V(beta)(信息矩阵I(Beta)的逆矩阵)
0.374733 0.324880 -0.257757 -0.192823
0.324880 0.323782 -0.244457 -0.169612
-0.257757 -0.244457 0.289488 0.069404
-0.192823 -0.169612 0.069404 0.236257
表2.4.4 有序分类因变量Logistic回归系数检验
序号回归系数系数标准误wald统计量自由度df 检验P值
常数项a1 -2.693576 0.612155 19.361377 1 0.000011
常数项a2 -1.812040 0.569018 10.141059 1 0.001450
性别 1.052352 0.538041 3.825528 1 0.050477
治疗方法 2.187272 0.486063 20.249800 1 0.000007
表2.4.5 有序分类因变量Logistic回归模型基本信息
样本分组数12
求解方法极大似然+牛顿迭代
迭代次数17
注意：该结果与 SAS, DPS不一致。

Poisson 回归模型
1 简介
一般情况下，单位容积水中的细菌数，单位时间内某些事件发生的次数，单位面积上降落的灰尘的颗粒数等，都可以用Poisson 分布来描述。

一般Poisson 分布描述成随机变量)(~λP Y ，概率分布律为：
,2,1,0,!
)(===-y y e
y Y P y
λλ
易知λ=EY ，通常λ可能受到众多因素的影响，不妨假设这些因素为k x x x ,,,21 (自变量，协变量），
令),,,,1(21k x x x X =，对于分组数据，Poisson 分布的期望发生数假设为[7]
：
i
X T i ik
x k i x i i i i e
n e
n X y E ββββλ===+++ 110)|( （1）
其中T
k ),,,(10ββββ =为回归参数，i n 为第i 组的总观测数。

回归模型的似然函数为Poisson 分
布条件下各个格子概率函数的乘积，因此Poisson 分布的极大似然函数和对数似然函数具体形式分别为：
∏
∏∏=∑=-=-====n
i i i
y i n i i n
i i i
y i
i
n i i y e
y e
p L 1
1
11
!
!
λλλλ
∑∑∑===-+-=n
i n
i i i i n
i i y y L 1
1
1
)!ln(ln ln λλ
代入i
X T i i e
n βλ=，得
∑∑∑∑∑∑======--+=--=n
i n
i i
y j i
X T i i T
i i i n
i n
i i
y j i
X T i i
X T i i j
e
n X y n y j
e
n e
n y L 1
11
111
ln ])ln([ln ]ln([ln ββββ （2）
令
∑=-=∂∂=n
i i
X T ig i ig i g g e X n X y L F 1][)(ln )(ββββ （3）
∑==∂∂∂-=n i i
X T ih ig i h g gh e
X X n L I 1
2)(ln )(βββββ,k h g ,,1,0, = （4） 则可采用Newton-Raphson 迭代法求解参数T
k ),,,(10ββββ =的极大似然估计了。

对于仅有常数项的Poisson 模型，其估计值为∑∑===n
i i
n
i i
n
y 11
ln ˆβ，用于计算对数似然比。

2．案例分析
例1[ 3]
Doll 和Hill(1966)研究英国男性医生患冠心病与抽烟、年龄关系。

由于死亡与追踪人数和追踪时间有关，故用追踪人数和追踪时间的乘积(人年)作为观察单位数。

假定其目标变量(死亡人数)近似服从Poisson 分布，其调查取样共74588调查单位，死亡598例。

主要研究因素有抽烟（1为抽烟,0为不抽烟）；调查对象年龄分成4组(35-44岁,45-54岁,55-64岁,65-74岁)，此为多分类变量，需要设置三个变量加以区分，可将其中一个年龄组作为参照组，不妨取35-44岁，计算时不考虑对照组信息。

表1 英国男性医生患冠心病与抽烟、年龄关系
分组 抽烟 34-44岁
45-54岁
55-64岁
65-74岁
死亡数 总例数 1 1 1 0 0 0 32 52307 2 1 0 1 0 0 104 43248 3 1 0 0 1 0 206 28612 4 1 0 0 0 1 186 12663 5 0 1 0 0 0 2 18790 6 0 0 1 0 0 12 10673 7 0 0 0 1 0 28 5710 8
1
28
2585
由LLLStat 软件计算的如下结果：
表2 回归系数方差矩阵V(beta)(信息矩阵I(Beta)的逆矩阵) 0.040354
-0.013325 -0.028763 -0.028467 -0.028496
-0.013325 0.016227
-0.000790 -0.001151 -0.001115
-0.028763 -0.000790 0.038071 0.029468 0.029466 -0.028467 -0.001151 0.029468 0.033767 0.029491 -0.028496 -0.001115 0.029466
0.029491
0.034161
表3 分组Poisson 回归系数检验
序号 均值 回归系数 系数标准误 wald 统计量 自由度df
检验P 值 常数项 1.000000 -8.036018
0.200882 1600.289462 1 0.000000 抽烟
0.500000
0.500463 0.127384 15.435202 1 0.000085 45-54岁 0.250000 1.475012 0.195118 57.147595 1 0.000000 55-64岁 0.250000 2.615085 0.183758 202.526084 1 0.000000 65-74岁 0.250000
3.338412
0.184828
326.246641
1
0.000000
表4 分组Poisson 回归模型基本信息
总组数 8
求解方法
极大似然+牛顿迭代 迭代次数(仅Beta0) 13(10) -2LogLikelihood(Beta) 7283.685428 仅常数项beta0
-5.676593
-2LogLikelihood(beta0) 7985.205751 方程Wald值(相减) 701.520324 方程自由度 4
方程检验p值0.000000
参考文献：
[1] 何晓群. 多元统计分析[M]. 北京：中国人民大学出版社,2008.9
[2] 金丕焕. 医用统计方法[M]. 上海：复旦大学出版社，2004.7
[3] 唐启义,冯明光. 实用统计分析及其DPS数据处理系统[M]. 北京：科学出版社, 2002.5
[4] Deniel A. Powers, 谢宇著,任强等译. 分类数据分析的统计方法[M].北京：社会科学文献出版社2009.7
[5] 高惠璇. 统计计算[M]. 北京：北京大学出版社,1995.7
[6] 关冶,陆金甫. 数值分析基础[M]. 北京：高等教育出版社,1998.5
[7] 沈其君. SAS统计分析[M]. 北京：高等教育出版社,2005.8。