1、Logistic 回归模型 1 Logistic 回归模型的基本知识 1.1 Logistic 模型简介 主要应用在研究某些现象发生的概率 ，比如股票涨还是跌，公司成功或失败的概率，以及讨论概p 率 与那些因素有关。显然作为概率值，一定有 ，因此很难用线性模型描述概率 与自变量的p 10p 关系，另外如果 接近两个极端值，此时一般方法难以较好地反映 p 的微小变化。为此在构建 与自变 量关系的模型时，变换一下思路，不直接研究 ，而是研究 的一个严格单调函数 ，并要求)(G 在 接近两端值时对其微小变化很敏感。于是 Logit 变换被提出来：)(G （1）pLogit1ln)( 其中当 从 时，
2、 从 ，这个变化范围在模型数据处理上带来很大的方便，p10)(pLogit 解决了上述面临的难题。另外从函数的变形可得如下等价的公式： （2）XTTepXpit 11ln)( 模型(2)的基本要求是，因变量（y）是个二元变量，仅取 0 或 1 两个值，而因变量取 1 的概率 就是模型要研究的对象。而 ，其中 表示影响 的第 个因素，它可)|1(XP TkxX),(21 ixyi 以是定性变量也可以是定量变量， 。为此模型(2)可以表述成：k0 （3）kx xepxp 1011ln 显然 ，故上述模型表明 是 的线性函数。此时我们称满足上面条件yE)( )(lnyEkx,21 的回归方程为 Lo
3、gistic 线性回归。 Logistic 线性回归的主要问题是不能用普通的回归方式来分析模型，一方面离散变量的误差形式服从 伯努利分布而非正态分布，即没有正态性假设前提；二是二值变量方差不是常数，有异方差性。不同于 多元线性回归的最小二乘估计法则(残差平方和最小) ，Logistic 变换的非线性特征采用极大似然估计的方 法寻求最佳的回归系数。因此评价模型的拟合度的标准变为似然值而非离差平方和。 定义 1 称事件发生与不发生的概率比为 优势比(比数比 odds ratio 简称 OR)，形式上表示为 OR= （4）kxep101 定义 2 Logistic 回归模型是通过极大似然估计法得到的
4、，故模型好坏的评价准则有似然值来表征， 称 -2 为估计值 的拟合似然度 ，该值越小越好，如果模型完全拟合，则似然值 为 1，而拟合ln()L ()L 似然度达到最小，值为 0。其中 表示 的对数似然函数值。()lnL 定义 3 记 为估计值 的方差-协方差矩阵， 为 的标准差矩阵，则称)(Var 21)()(VarS （5）kiSwi ,21, 为 的 Wald 统计量，在大样本时， 近似服从 分布，通过它实现对系数的显著性检验。i i )( 定义 4 假定方程中只有常数项 ，即各变量的系数均为 0，此时称0 （6）20ln()l()L 为方程的显著性似然统计量，在大样本时， 近似服从 分布
5、。22k 1.2 Logistic 模型的分类及主要问题 根据研究设计的不同，Logistic 回归通常分为成组资料的非条件 Logistic 回归和配对资料的条件 Logistic 回归两种大类。还兼具两分类和多分类之分，分组与未分组之分，有序与无序变量之分。具体如 下： 两分类非条件 Logistic 回归：分组数据的 Logistic 回归，未分组数据的 Logistic 回归； 多分类非条件 Logistic 回归：无序变量 Logistic 回归，无序变量 Logistic 回归； 条件 Logistic 回归：1:1 型、1:M 型和 M:N 型 Logistic 回归。 关于 L
6、ogistic 回归，主要研究的内容包括： 1 模型参数的估计及检验 2 变量模型化及自变量的选择 3 模型评价和预测问题 4 模型应用 2 Logistic 模型的参数估计及算法实现 2.1 两分类分组数据非条件 Logistic 回归 因变量(反应变量)分为两类，取值有两种，设事件发生记为 y=1，不发生记为 y=0，设自变量 是分组数据，取有限的几个值；研究事件发生的概率 与自变量 的关TkxX),21 )|1(XyP 系，其 Logistic 回归方程为： 或 kxXyP10)|(ln kxxey10)|( 例 2.1.1 分组数据 1 在一次住房展销会上，与房地产商签订初步购房意向书
7、的有 n=325 人，在随后的 3 个月时间内，只有一部分顾客购买了房屋。购买房屋的顾客记为 1，否则记为 0。以顾客的年家庭收入(万 元)作为自变量 ,对数据统计后如表 2.1.1 所示，建立 Logistic 回归模型。 表 2.1.1 购房分组数据 序号 年家庭收入 X(万元) 签订意向人数 实际购买人数 1 1.5 25 8 2 2.5 32 13 3 3.5 58 26 4 4.5 52 22 5 5.5 43 20 6 6.5 39 22 7 7.5 28 16 8 8.5 21 12 9 9.5 15 10 例 2.1.2 药物疗效数据 2 为考察某药物疗效，随机抽取 220 例
8、病人并分配到治疗组和对照组，治疗组采 用治疗药物，对照组采用安慰剂。治疗一段时间后观察病人的疗效，得到表 2.1.2 数据。设 y 为疗效指标 (y=1 有效，y=0 无效)， 为治疗组指标(1 为治疗组，0 为对照组) ， 为年龄组指标(1 为45 岁，0 为其1x 2x 他)。 表 2.1.2 药物疗效数据 序号 治疗分组 1x年龄分组 2x有疗效 无效 合计 1 1 1 32 18 50 2 1 0 40 20 60 3 0 1 21 31 52 4 0 0 18 40 58 上述两个例子数据都是经过统计加工后的分组数据，对此类数据进行 Logistic 回归，首先要明确应变 量对应事件
9、的发生概率如何确定和进行 Logit 变换，其次才能建立 Logistic 回归。为便于数据处理，我们 将此类数据的格式作个约定，排列格式为（组序号，自变量 ，该组事件发生数，该组总例数） 。X 表 2.1.3 分组数据的标准格式 表 2.1.1 改造表 序 号 年家庭收入 X(万元) 实际购买 人数 im签订意向总人数 in 1 1.5 8 25 2 2.5 13 32 3 3.5 26 58 4 4.5 22 52 5 5.5 20 43 6 6.5 22 39 7 7.5 16 28 8 8.5 12 21 9 9.5 10 15 表 2.1.2 改造表 序 号 治疗分 组 1x年龄分
10、组 2 有效例 数 im观察例 数 in 1 1 1 32 50 2 1 0 40 60 3 0 1 21 52 4 0 0 18 58 经过改造后，可得我们关心的事件的发生的频率为 。其中ninmpi ,21,i该 组 总 例 数该 组 发 生 事 件 数 为分组数，然后作 Logit 变换，即 。变换后的数据，形式上已经可以采用一n iiiLogtp1ln)( 般的线性回归的处理方式来估计回归参数了。此时方程变为： kjiji nxp10 ,2, 当然这样处理并没有解决异方差性，当 较大时， 的近似方差为：ini （7）)(,)1()( iiiii yEpD 所以选择权重 ，最后采用加权最
11、小二乘法估计参数。nipnii ,2),1( 注意，分组数据的 Logistic 回归只适用于大样本分组数据，对小样本的为分组数据不适用，并且以组 数 为回归拟合的样本量，明显降低了拟合精度，在实际应用中必须谨慎。n 求解算法及步骤： 1依据分组数据的标准格式，计算频率 、Logit 变换 和权重ipipi 2构建加权最小二乘估计： （8） ni kjijiiini kjiji xyxy1 1201120 )(m)(m 令 ， ，ii* Tikiii xX),(1* Tk),(10 则方程又变成一般的线性回归模型： （9ni iTiXy12*)( ） 3构造增广矩阵 利用消去法得 矩阵，得到估
12、计21*kTTYX)(VarI 其中 为残差平方和 , 回归方差2,1KISE1knSE 各系数检验采用 )1(kntItii 总平方和 ，回归平方和niniiiiyST1122)(SETR 总平方和求解相当于拟合 方程的残差平方和，故得上式 STiiy*0 所以方程的检验为 )1,()1/(knFknSERF 例 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.440
13、000 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.14
14、2857 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 家
15、庭年收入x 14.901140 0.149323 0.020711 7.209865 7 0.000056 表2.1.7 1XT 0.086479 -0.014517 -0.014517 0.002876 本例 Logistic 模型的回归方程： xepi 14932.08. .42.01 对于多分类无序自变量的 Logistic 回归，即某个自变量为 m 个水平的名义变量（如治疗方法 A,B,C） ， 只需要引入 m-1（2 个）个哑变量，然后采用上述方法进行分析。 例 2.1.3 研究三种治疗方法对不同性别病人的治疗效果 2，数据如表 2.1.4 表 2.1.4 性别和治疗法对某病治愈情况
16、的影响 性别 治疗方法 有效 im无效 总例数 in A 78 28 106 B 101 11 112男 C 68 46 114 A 40 5 45 B 54 5 59女 C 34 6 40 由于治疗方法有三种，没有等级关系，所以属于无序的名义变量，故引入两个哑变量 分别代表32,x A 和 B 疗法，其中 表示方法 A, 表示方法 B, 表示方法 C，将0,132x1,032x0,32x 上述数据转化成标准格式，得表 2.1.5。 表 2.1.5 性别和治疗法对某病治愈情况的影响 性别 1x23x有效 im总例数 in 1 1 0 78 106 1 0 1 101 112 1 0 0 68
17、114 0 1 0 40 45 0 0 1 54 59 0 0 0 34 40 对于分类数据，也可以采用极大似然法进行参数估计，具体见 2.2 节最后部分内容。 2.2 两分类未分组(连续)非条件 Logistic 回归 应变量 取值为 0 和 1，设事件发生记为 y=1，否则为 0，设自变量 ，n 组观测y Tkxx),(21 数据记为 , 。记 ， ，则 与),(21ikiyx n,2TikiixX),1(20iiy 的 Logistic 回归模型是：ikix,21 （10）iXT iikxixikiii eexxfyE 11)()( 1010 易知， 是均值为 的 0-1 型分布，其分布
18、律为 ii ，iyiiyif 1)()(nii ,21;,0 则 的似然函数和对数似然函数分别为： ny,21 ni iyiyL11)( ni iiini iiii yyyL11 )ln(l)l()(ll 代入 ,得ikxi xie10 （11） ni iXTiTii ikxiikii eyexL1 1010)ln( )ln(l 记 ，选取 的估计 使得 达到极大，)(l)(LTk),(10 Tk),(10)(L 这就是 Logistic 回归模型的极大似然估计，该过程的求解需要采用牛顿迭代法。 构造得分函数 ，共 k+1 个非线性方程组，令其=0 求解 ，其中gFg ,2,)()( (12
19、)kgexy ni iXTiigig ,210,1)(0 构造信息矩阵 ，即 二阶导矩阵的负矩阵，其中khLIgh ,2,)()(2)(L (13 )khgexIniiXTihiggh ,210,)1()(02 很明显 ，故 是一个对称矩阵。)()(hgghII)I 求解算法及步骤： 1 根据公式(12 ) 计算得分函数 ，公式（13 ）计算信息矩阵)(gF)(ghI 给定初值 , k =1 和精度 ，可取 0.0000010,(0 2 采用牛顿迭代式 , ，通过以下方式求解。1k )()(11kkFI 构造增广矩阵 = ，通过对 IF 矩阵作 k+1 次 ij 消去变换求解)(I)(kI 若
20、 或者 或者 ，则转 3kg02| kg0| |max0gkg 否则 k = k +1，继续执行第 2 步 3 此时 就是回归系数 的数值估计 ，k 就是迭代次数，消去变换后的 矩阵的前k IF 子阵就是 方差-协方差矩阵的估计阵 =V ，下面给出检验有关计算：1k 1)(kghVar 计算 Wald 统计量 ,近似服从 分布，检验 p 值 gW2)1(2)(2ggWP 标准误 , , ggVES).(geOR)(k,0 例 2.2.1 公共交通调查数据 1 在一次关于公共交通的社会调查中，调查项目为“是乘坐公共汽车上下班， 还是骑自行车上下班” 。因变量 y=1 表示乘坐公共汽车，y=0 表
21、示骑自行车。自变量 是年龄，作为连续1x 变量； 是月收入(元)； 是性别， =1 表示男性， =0 表示女性。调查对象为工薪族群体，数据如2x3x33x 表 2.2.1 所示。 表2.2.1 公共交通社会调查 序号 年龄 1x月收入 2性别 3交通 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
22、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 0 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统计量 自
23、由度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 求解方法 极大
24、似然法 ,1,1)|( mjnieXjyPXTjajii （19）,)|( |ln1 jkyPiTjjkiijkii 有序累积概率模型参数的极大似然估计就是寻找参数使得联合概率实现最大化，由于观测之间 相互独立，联合概率被分解成边缘概率之积。而观测到 的概率就是累积概率之差：jyi)|1()|()|( iiiiii XjPXjPXjyP 第 i 个观测值对应似然值的贡献取决于观测到哪一个 j 值，因此对于次序响应的每个 j 值，取所 有 的观测之的乘积，有似然函数：jyi ，其中若 ，则 ，否则 nimj ijdiXjyPL1)|( jyi1ijd0ijd 并且对于任一个观测 而言，只有一个等
25、级事件发生，即 ，故有(19)式。其iXmj iiXjP1)|( 对数似然函数如下（对于分组数据，似然函数变为： ， 分组中各分类例 nij ijdniijyL1)|(i 数） 。 （20)11ln()1ln(1ln )|(ll 21 i mj iXTjaijiXTjaijijiXTmaiimiXTaiinimj iiij eedededjyPL ） 其中： mjjjeeeXjyPiXTmaiiXTjaijiTjaiXjiTai 111)|( 1 （21） 然后就可以通过极大似然法，就上 Newton-Raphson 方法加以求解参数 了，注意的是,11ma 。下面给出具体推导 求解的详细过程
26、。对（20）式进行化简，可得121maa ,11ma )ln()ln( )ln(ln 111 iXTaimiXTaiTiij iXTjaijjjiij edeXd eL （22） ni iXTaiaiiXTai eeeaL1 1121 )(l （23） ni iXTmaiiiXTmaiamim edeed1 1121 )(l （24） 2,)1()1(ln1 1 mgeeeeaLi iXTgaigaigiXTgaigaigg （25） kg eedededx ni mj iXTjaijiXTjaijijiXTaimiiXTaig,21 ,)11(11l 2 （26） ni iXTaiaiiXT
27、aii eededaL1 21212112 )()()(l （27） ni iXTmaiiiXTmaimaimm edee1 21212112 )()()(l （28） 2, )1()()1()( ln1 221222 mg eedeedaLni iXTgaigaigiXTgaigaig （29） 2,1,)(ln1212 mgedaLni agig （30） khgedxani iXTgaiigihhg ,21;,1,)1()(l1 22 （31） khg eedededxLni mj iXTjaijiXTjaijijiXTmaiiiXTaiihg,21, )1()1()1()1(l 2 2
28、222 （32） 由此构建信息矩阵 和 ，并可迭代求解了。注：若为分组数据，上述每项乘以 。),(aI),(F in 例 2.4.1 研究性别和两种治疗方法对某种疾病疗效的影响 3,将疗效分成效果显、有效和无效三个等 级，根据试验调查，得到如下资料。 表 2.4.1 性别和两种治疗方法对某种疾病疗效的影响 性别 治疗方法 显著 有效 无效 合计 新药 16 5 6 27女 传统 6 7 19 32 新药 5 2 7 14男 传统 1 0 10 11 表2.4.2 多分类有序反应变量数据格式 行号 性别 治疗方法 频数 疗效等级 1 1 1 16 1 2 1 1 5 2 3 1 1 6 3 4
29、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
30、.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 求解方法 极大似然+牛顿迭
31、代 迭代次数 17 注意：该结果与 SAS, DPS 不一致。 Poisson 回归模型 1 简介 一般情况下，单位容积水中的细菌数，单位时间内某些事件发生的次数，单位面积上降落的灰尘的 颗粒数等，都可以用 Poisson 分布来描述。一般 Poisson 分布描述成随机变量 ，概率分布律为：)(PY,210,!)(yeyYP 易知 ，通常 可能受到众多因素的影响，不妨假设这些因素为 (自变量，协变EY kx,21 量） ，令 ，对于分组数据，Poisson 分布的期望发生数假设为 7：),1(2kxX （1）iX Tiikxiiiii enenXyE10)|( 其中 为回归参数， 为第 i
32、组的总观测数。回归模型的似然函数为 Poisson 分Tk),(10i 布条件下各个格子概率函数的乘积，因此 Poisson 分布的极大似然函数和对数似然函数具体形式分别为： niiiyiniiiyni epL111 !niniii yy11)l(ll 代入 ，得iX Tiien （2） ni niiyjiXTiiTiiini niiyjiiiXTii eyeL1 11 1l)l( lll 令 （3） ni iXTigiggg eyLF1)(l)( , （4）ni iXTihghgghI12)(l)( k,10 则可采用 Newton-Raphson 迭代法求解参数 的极大似然估计了。Tk),
33、(0 对于仅有常数项的 Poisson 模型，其估计值为 ，用于计算对数似然比。nii y10l 2案例分析 例 1 3 Doll 和 Hill(1966)研究英国男性医生患冠心病与抽烟、年龄关系。由于死亡与追踪人数和追踪 时间有关，故用追踪人数和追踪时间的乘积(人年) 作为观察单位数。假定其目标变量 (死亡人数)近似服从 Poisson 分布，其调查取样共 74588 调查单位，死亡 598 例。主要研究因素有抽烟（ 1 为抽烟,0 为不抽烟） ； 调查对象年龄分成 4 组(35-44 岁,45-54 岁,55-64 岁,65-74 岁) ，此为多分类变量，需要设置三个变量加以区 分，可将其
34、中一个年龄组作为参照组，不妨取 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 0 0 0 0 1 28 2585 由 LLLStat 软件计算的如下结果： 表2 回归系数方差
35、矩阵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回归系数检验 序号 均值 回归系数 系数标准
36、误 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 分组P
37、oisson回归模型基本信息 总组数 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
