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)kxxepxp 1011ln显然 ,故上述模型表明 是 的线性函数。此时我们称满足上面条件yE)( )(lnyEkx,21的回归方程为 Logistic 线性回归。Logist
3、ic 线性回归的主要问题是不能用普通的回归方式来分析模型,一方面离散变量的误差形式服从伯努利分布而非正态分布,即没有正态性假设前提;二是二值变量方差不是常数,有异方差性。不同于多元线性回归的最小二乘估计法则(残差平方和最小) ,Logistic 变换的非线性特征采用极大似然估计的方法寻求最佳的回归系数。因此评价模型的拟合度的标准变为似然值而非离差平方和。定义 1 称事件发生与不发生的概率比为 优势比(比数比 odds ratio 简称 OR),形式上表示为OR= (4)kxep101定义 2 Logistic 回归模型是通过极大似然估计法得到的,故模型好坏的评价准则有似然值来表征,称-2 为估
4、计值 的拟合似然度 ,该值越小越好,如果模型完全拟合,则似然值 为 1,而拟合ln()L ()L似然度达到最小,值为 0。其中 表示 的对数似然函数值。()lnL定义 3 记 为估计值 的方差-协方差矩阵, 为 的标准差矩阵,则称)(Var 21)()(VarS(5)kiSwi ,21,为 的 Wald 统计量,在大样本时, 近似服从 分布,通过它实现对系数的显著性检验。i i )(定义 4 假定方程中只有常数项 ,即各变量的系数均为 0,此时称0(6)20ln()l()L为方程的显著性似然统计量,在大样本时, 近似服从 分布。22k1.2 Logistic 模型的分类及主要问题根据研究设计的
5、不同,Logistic 回归通常分为成组资料的非条件 Logistic 回归和配对资料的条件Logistic 回归两种大类。还兼具两分类和多分类之分,分组与未分组之分,有序与无序变量之分。具体如下:两分类非条件 Logistic 回归:分组数据的 Logistic 回归,未分组数据的 Logistic 回归;多分类非条件 Logistic 回归:无序变量 Logistic 回归,无序变量 Logistic 回归;条件 Logistic 回归:1:1 型、1:M 型和 M:N 型 Logistic 回归。关于 Logistic 回归,主要研究的内容包括:1 模型参数的估计及检验2 变量模型化及自
6、变量的选择3 模型评价和预测问题4 模型应用2 Logistic 模型的参数估计及算法实现2.1 两分类分组数据非条件 Logistic 回归因变量(反应变量)分为两类,取值有两种,设事件发生记为 y=1,不发生记为 y=0,设自变量是分组数据,取有限的几个值;研究事件发生的概率 与自变量 的关TkxX),21 )|1(XyP系,其 Logistic 回归方程为:或 kxXyP10)|(ln kxxey10)|(例 2.1.1 分组数据 1 在一次住房展销会上,与房地产商签订初步购房意向书的有 n=325 人,在随后的 3个月时间内,只有一部分顾客购买了房屋。购买房屋的顾客记为 1,否则记为
7、0。以顾客的年家庭收入(万元)作为自变量 ,对数据统计后如表 2.1.1 所示,建立 Logistic 回归模型。表 2.1.1 购房分组数据序号 年家庭收入 X(万元) 签订意向人数 实际购买人数1 1.5 25 82 2.5 32 133 3.5 58 264 4.5 52 225 5.5 43 206 6.5 39 227 7.5 28 168 8.5 21 129 9.5 15 10例 2.1.2 药物疗效数据 2 为考察某药物疗效,随机抽取 220 例病人并分配到治疗组和对照组,治疗组采用治疗药物,对照组采用安慰剂。治疗一段时间后观察病人的疗效,得到表 2.1.2 数据。设 y 为疗
8、效指标(y=1 有效,y=0 无效), 为治疗组指标(1 为治疗组,0 为对照组) , 为年龄组指标(1 为45 岁,0 为其1x 2x他)。表 2.1.2 药物疗效数据序号 治疗分组 1x年龄分组 2x有疗效 无效 合计1 1 1 32 18 502 1 0 40 20 603 0 1 21 31 524 0 0 18 40 58上述两个例子数据都是经过统计加工后的分组数据,对此类数据进行 Logistic 回归,首先要明确应变量对应事件的发生概率如何确定和进行 Logit 变换,其次才能建立 Logistic 回归。为便于数据处理,我们将此类数据的格式作个约定,排列格式为(组序号,自变量
9、,该组事件发生数,该组总例数) 。X表 2.1.3 分组数据的标准格式表 2.1.1 改造表序号年家庭收入X(万元)实际购买人数 im签订意向总人数 in1 1.5 8 252 2.5 13 323 3.5 26 584 4.5 22 525 5.5 20 436 6.5 22 397 7.5 16 288 8.5 12 219 9.5 10 15表 2.1.2 改造表序号治疗分组 1x年龄分 组 2有效例数 im观察例 数 in1 1 1 32 502 1 0 40 603 0 1 21 524 0 0 18 58经过改造后,可得我们关心的事件的发生的频率为 。其中ninmpi ,21,i该
10、 组 总 例 数该 组 发 生 事 件 数为分组数,然后作 Logit 变换,即 。变换后的数据,形式上已经可以采用一n iiiLogtp1ln)(般的线性回归的处理方式来估计回归参数了。此时方程变为: kjiji nxp10 ,2,当然这样处理并没有解决异方差性,当 较大时, 的近似方差为:ini(7))(,)1()( iiiii yEpD所以选择权重 ,最后采用加权最小二乘法估计参数。nipnii ,2),1(注意,分组数据的 Logistic 回归只适用于大样本分组数据,对小样本的为分组数据不适用,并且以组数 为回归拟合的样本量,明显降低了拟合精度,在实际应用中必须谨慎。n求解算法及步骤
11、:1依据分组数据的标准格式,计算频率 、Logit 变换 和权重ipipi2构建加权最小二乘估计:(8) ni kjijiiini kjiji xyxy1 1201120 )(m)(m 令 , ,ii* Tikiii xX),(1* Tk),(10则方程又变成一般的线性回归模型: (9ni iTiXy12*)()3构造增广矩阵 利用消去法得 矩阵,得到估计21*kTTYX)(VarI 其中 为残差平方和 , 回归方差2,1KISE1knSE各系数检验采用 )1(kntItii总平方和 ,回归平方和niniiiiyST1122)(SETR总平方和求解相当于拟合 方程的残差平方和,故得上式 STi
12、iy*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.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.6
13、92308 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 方
14、程检验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 1XT0.086479 -0.014517 -0.014517 0.002876 本例 Logistic 模型的回归方程: xepi 14932
15、.08.42.01对于多分类无序自变量的 Logistic 回归,即某个自变量为 m 个水平的名义变量(如治疗方法 A,B,C) ,只需要引入 m-1(2 个)个哑变量,然后采用上述方法进行分析。例 2.1.3 研究三种治疗方法对不同性别病人的治疗效果 2,数据如表 2.1.4表 2.1.4 性别和治疗法对某病治愈情况的影响性别 治疗方法 有效 im无效 总例数 inA 78 28 106B 101 11 112男C 68 46 114A 40 5 45B 54 5 59女C 34 6 40由于治疗方法有三种,没有等级关系,所以属于无序的名义变量,故引入两个哑变量 分别代表32,xA 和 B
16、疗法,其中 表示方法 A, 表示方法 B, 表示方法 C,将0,132x1,032x0,32x上述数据转化成标准格式,得表 2.1.5。表 2.1.5 性别和治疗法对某病治愈情况的影响性别 1x23x有效 im总例数 in1 1 0 78 1061 0 1 101 1121 0 0 68 1140 1 0 40 450 0 1 54 590 0 0 34 40对于分类数据,也可以采用极大似然法进行参数估计,具体见 2.2 节最后部分内容。2.2 两分类未分组(连续)非条件 Logistic 回归应变量 取值为 0 和 1,设事件发生记为 y=1,否则为 0,设自变量 ,n 组观测y Tkxx)
17、,(21数据记为 , 。记 , ,则 与),(21ikiyx n,2TikiixX),1(20iiy的 Logistic 回归模型是:ikix,21(10)iXTiikxixikiii eexxfyE 11)()( 1010 易知, 是均值为 的 0-1 型分布,其分布律为 ii,iyiiyif 1)()(nii ,21;,0则 的似然函数和对数似然函数分别为: ny,21 ni iyiyL11)( ni iiini iiii yyyL11 )ln(l)l()(ll 代入 ,得ikxixie10(11)ni iXTiTii ikxiikii eyexL1 1010)ln( )ln(l 记 ,选
18、取 的估计 使得 达到极大,)(l)(LTk),(10 Tk),(10)(L这就是 Logistic 回归模型的极大似然估计,该过程的求解需要采用牛顿迭代法。构造得分函数 ,共 k+1 个非线性方程组,令其=0 求解 ,其中gFg ,2,)()( (12 )kgexyni iXTiigig ,210,1)(0 构造信息矩阵 ,即 二阶导矩阵的负矩阵,其中khLIgh ,2,)()(2)(L(13 )khgexIniiXTihiggh ,210,)1()(02很明显 ,故 是一个对称矩阵。)()(hgghII)I求解算法及步骤:1 根据公式(12 ) 计算得分函数 ,公式(13 )计算信息矩阵)
19、(gF)(ghI给定初值 , k =1 和精度 ,可取 0.0000010,(02 采用牛顿迭代式 , ,通过以下方式求解。1k )()(11kkFI构造增广矩阵 = ,通过对 IF 矩阵作 k+1 次 ij 消去变换求解)(I)(kI 若 或者 或者 ,则转 3kg02| kg0| |max0gkg否则 k = k +1,继续执行第 2 步3 此时 就是回归系数 的数值估计 ,k 就是迭代次数,消去变换后的 矩阵的前k IF子阵就是 方差-协方差矩阵的估计阵 =V ,下面给出检验有关计算:1k 1)(kghVar计算 Wald 统计量 ,近似服从 分布,检验 p 值 gW2)1(2)(2gg
20、WP标准误 , , ggVES).(geOR)(k,0例 2.2.1 公共交通调查数据 1 在一次关于公共交通的社会调查中,调查项目为“是乘坐公共汽车上下班,还是骑自行车上下班” 。因变量 y=1 表示乘坐公共汽车,y=0 表示骑自行车。自变量 是年龄,作为连续1x变量; 是月收入(元); 是性别, =1 表示男性, =0 表示女性。调查对象为工薪族群体,数据如2x3x33x表 2.2.1 所示。表2.2.1 公共交通社会调查 序号 年龄 1x月收入 2性别 3交通 y1 18 850 0 02 21 1200 0 03 23 850 0 14 23 950 0 15 28 1200 0 16
21、 31 850 0 07 36 1500 0 18 42 1000 0 19 46 950 0 110 48 1200 0 011 55 1800 0 112 56 2100 0 113 58 1800 0 114 18 850 1 015 20 1000 1 016 25 1200 1 017 27 1300 1 018 28 1500 1 019 30 950 1 120 32 1000 1 021 33 1800 1 022 33 1000 1 023 38 1200 1 024 41 1500 1 025 45 1800 1 126 48 1000 1 027 52 1500 1 12
22、8 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
23、 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 分组数据的极大似然估计法,主要过程如下:ni imniimiCL1)( ni iiiiimni iiiii n
24、Cii 11 )1l(1ll)ln()(lll 代入 ,得 ikxixie10 i iXTiiTimneXLi1 )l(ll 则有 ggLF)()( kgexnmni iXTiigig ,20,1 ;hIgh)()(2hniiXTihig ,1,)(12其中 , 分别表示分组 i 中事件发生次数和总观察数,如表 2.1.4 和 2.1.5 所示。然后可采用imnNewton-Raphson 迭代法进行求解。由 LLLStat 计算得到如下结果。表 2.2.4 性别和疗法对某病治愈的影响(未分组 Logistic 似然估计法)序号 均值 回归系数 系数标准误 wald统计量 自由度df 检验P值
25、常数项 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
26、.000078 0.000128-0.029931 -0.000078 0.069753 0.029993-0.030097 0.000128 0.029993 0.0998312.3 条件 Logistic 回归 2,3条件 Logistic 回归是配对设计(病例- 对照)中常用的一种统计分析方法,通过配对方法收集资料:每一配对组可包括一个病例和一个或多个对照,有 1:1 型、1:m 型配对。假设收集了如下数据:表 2.3.1 n 个 1:m 配对组,k 个协变量的比例资料 配对组号 病例组 0X第 1 对照组 1X 第 m 个对照组 mX1 1201,kx 2,kx kx121,2 0 1
27、1 n0201,nknx 112,nknx mnkmnx,21配对资料用配对的方法来控制影响因素的干扰,并且每个配对组都可以建立一个 Logistic 回归方程:ipLogit ki ,)(10 为此需要估计的参数有 n 个常数项 和 k 个回归系数 ,配对数越多估计的参数就n1, k1越多,但是一般的数据量难以支撑这样的估计,故一般的 Logistic 回归不适合配对资料。不过在参数估计时,常数项会被消去,所以方程组减少了 n 个常数项 的估计,复杂度大大降低。对于回归参数n01,的估计采用条件似然函数替代一般的似然函数进行。对于第 i 个配对组而言,共有 m+1 个观察对象,记为 ,其中仅有一例发病,且正好mBA,21是病例组 A 发病,而对照组均没有发病的条件概率 (类似 Bayes 概率)可以表示成:ip(14)mj mji BAPBAPp1212)()()( 其中 = ,而)(21mBAP |0|0(|00 iiiiii XyXyy, , (15 )jiXTjijiji e1| jiXTjiji eP1)| ,21