二项Logistic回归参数最大似然估计的计算资料.docx

1、二项Logistic回归参数最大似然估计的计算资料二项Logistic回归参数最大似然估计的计算1 Logistic回归的基本思想在地丧风险评估中，研究者往往关心地震发生时，地表发生破裂的概率，地表破裂是 由哪些因素引起的等问题。利用以往的相关数据找岀统讣规律性来解决这些问题，实质上可以转化为一个多元回 归分析问题y=，其中召,农，兀r，为随机变量。由于因变量丫的取值只有两个状态：破裂(y = i)和不破裂(r = o),因此直接寻找因变 量丫和自变量x的关系非常困难。于是，可以把研究问题转换一个角度，不去直接分析Y和 *的关系，而是分析条件概率p(r=i|x和x的关系，这等价于寻找一个取值在

2、o到1之 间的连续函数P(x) = PY = l|x o数学上满足这种条件的函数存在且不唯一,Logistic回归就是满足这种要求的函数之一。 和线性回归分析类似.Logistic回归基本原理就是利用一组观测数据拟合一个Logistic模型, 然后借助这个模型来揭示总体中若干个自变量与一个因变疑取每个值的概率之间的依存关 系，并评估用这一模型模拟相关事物变化规律的准确性。具体地说，Logistic回归分析可以 从统计意义上确左在消除了其它变量的影响后，每一个自变量的变化是否引起因变量取某个 值的概率的变化，并估计出在其它自变量固定不变的情况下，每个自变量对因变量取某个值 的概率的数值影响大小。

3、在使用Logistic回归分析前，需要明确模型的使用条件：1、要求因变量是分类变量， 包括顺序变量和划义变量，不管哪种变量，都要用数字表示它，如可以令丫 = 1表示地震发 生时地表破裂，令7 = 0表示地丧发生时地表未破裂：2、自变量可以是数值型连续变量, 如地丧的震级，(ii)顺序变量,如覆盖层的厚度分组(1020m, 20-40m等)，(iii)名义变量, 如地箴类型，可令走滑型地震为1,正断型地震为2,逆冲型地震为3。2多元二项Logistic回归模型的定义由于地震发生时地表是否破裂受到多个因素的影响，故引入多元Logistic回归模型。假 设因变量F是一个取值为1和0的二值变量，兀=為

4、宀，,xj是影响丫的R个因素，回 归系数“ =0。胡,0J ,则Y关于X的k元Logistic回归模型定义为 心)= PY = ix= cxp(0(严阳+0壬+. + 处丿=exp(L 八卩)(1)1 + exp(仇 + 0內 + /32x2 + +0內)1 + exp( 1, x1 fl)由式(1)可得PY = Ox= T1 l + exp（l, x7fi）3 Logistic回归参数估计我们用最大似然估计方法去求模型的参数。再假设从总体（Y,x）中抽取一个容量为耳+”2的随机样本（1,旺）,（1山2），（1，乙”（0心严1）,（0,%严2），（心严），其中=（兀，兀2,心），j = 12M

5、+s 则有似然函数为/1| 眄+丿|,(0) = 口川31|盯 P(eo|/-Ii 叫 +1=亓exp(炕+0厂0N厂+久$ ) 1亓 1 V1 + cxp(0o +朋| + 02看2 + 0皿)去占1 + exp(几+ 0丙+戸2齐2 +以心) =计 cxp(l.打 M)，幵 一 * 1 + cxp(L x,r0)七 1 + exp(L xj|0)两边取对数，整理可得In UP)=力(0()+ 0E + 际2 + + pkxik) 一工 lnl + exp(A, + xn + p2xi2 + + pkxik)写成向量形式为/tj +ru讥(0) = 1, x；p- lnl + exp(l,

6、xTfi) (4 J为求式（4）的驻点，可求对数似然函数InL（fl）关于“的似然方程组为=叫-/-I =01 + exp(-0()- 0 0)exp(0D + 0+ + A.q)1 + exp(久+ 0心+ 0丹)cln Lip) k十也= Z*Vh-Zi1 r-1心exp(0 + Q.D+ + 0H1 + exp( 0。+ Q兀+卩皿)1-! 1-11 +exp（-几一禹心 =00心cln Lip)n, n, +n2= Z-vrt-S il r-1.J exp(炕-0K：+ 0_土, l + exp(0|)+ 0.+ + 0血)i.i n i-i 1 + exp(-0仆一0.勺 0內)写成

7、向量形式为卄乞* = 0800 ti l + exp(-0)泌邑补 儿=。切 tr 台 l + exp(-xQ空込邑芝 “ =0氓 山台 l+exp（*）式（5）为非线性方程组，一般情况下没有解析解，可以用Ncwton-Raphson迭代方法求 其数值解，令F()=I l+CXp(-0-0內) -1! 坷.小 丫 J 山 I 1 +CXP（-仇一工0內）ME InA _ Y 匕 l+cxpi”)叫 叫： ry v - y 亠：zr11 幺 i+cxp(-x)=0m 円+心 rx-z I - l+CXp(-0_2A)j则F（“）关于“的Jacobian矩阵为5 1 + 心 vy v - y 2i

8、i 厶 A tr l + expc-x,#)i-iJ =exp(-几-几耳)-Z z l + exp(-0o-M.%)內杯环exp(-几一0声丿 ；z l + exp(几一士0几)F(屮“、心exp(-0“一Q.%)-2 rI 1 + exp(-几一儿)r4-11屮、時 exp(-A, - 內) z 11 + exp(-久-乞0 儿)FJ-!c计.总exp(-几一0心)V 】 / I l + exp(m-0r)F j-i1. 1 + exp(-鶴一乞0再厅 j-i入 exp(-0。-0%) /-|-Z x“ i+exp(-A.-ZM；)l2An.林乙exp（-炕一0人）“ 1+exp（-久-乞

9、0儿）F-j-i导 exp(-“0它心 exp(-“) 台l + exp(-x/)I2它心 exp(-x0 Z? 1 + exp(-x,0)Fm |l+exp(-x/)2exp(-兀0) 纟1+exp(-x(0F xJexp(-M0Zf l + exp(-x/)I2_V AnA；Aexp(-x/Z l + exp(-x.0)FW .和 exp(-x,0)旨 从环 exp(-x,“)它時 exp(-) 台l + exp(-x/)Fl + exp(-x,/?)rm l+expC-x/)2屮与叮p（/r）F（/r）M = oj2：Step 3:若V，即满足容许的精度，则结束，否则更新参数p=pn n

10、 = n + ,转至Stcp2当给左地後发生时，地表破裂是否发生的数据时，根据上而的算法，可以求出参数的最 大似然估计。我们按照上述算法用MATLAB编写了多元Logistic回归参数估计的程序，可以给出参数 估计值，详见附录。附录1用Newton-Raphson方法求解参数，附录2用直接优化对数似然函数方法求解参数, 附录3用MATLAB自带的广义回归模型估计参数。附录4是别人做的Logistic回归的例子, 用的软件是SAS （种经过美国FDA认证的昂贵的商业统汁软件）得到的结果。附录5是 SPSS操作的过程和得到的结果。附录 1: Mat lab Files for Logistic R

11、egression%假设我们有一个数据，45个观测值，四个变量，包括：% 1. age （年龄，数值型）：% 2. vision （视力状况，分类型，1表示好，0表示有问题）：% 3. drive （驾车教育，分类型，1表示参加过驾车教冇，0表示没有）和% 4. 一个分类型输出变M accident （去年是否岀过事故，1表示岀过事故，0表示没有）。%我们的目的就是要考察前三个变量与发生事故的关系。% 第 1 至 4 列分别为 accident age vision drive ;cleanclc,close alldata = 1 17 1 11 44 0 01 48 1 01 55 0 0

12、1 75 1 10 35 0 10 42 1 10 57000 28 0 10 20 0 10 38 1 0 0 45 0 1 0 47 1 10 52000 55 0 11 68 1 01 18 1 01 68 0 01 48 1 11 17001 70 1 11 72 1 01 35 0 11 19 1 01 62 1 00 39 1 10 40 1 10 55000 68 0 10 25 1 00 17000 45 0 10 44 0 10 67000 55 0 11 61 1 01 19 1 01 69 0 01 23 1 11 19001 72 1 11 74 1 01 31 0 1

13、1 16 1 01 61 1 0;Y = data(: J);X= data(:,3:4);beta0 = 0.110 ;17137 ; -1.5000+1 *rand(3,1 );%rand(4,1); %猜测的初始值%自带的fsolvc算法求解没有问题，核心原理也是Newton-Raphson算法%options = optimset(,DisplayiterTolFun 1 e-8)beta = fsoIve(be)LogisticF(be,Y,X),betaO)%,options)%自编的Newton-Raphson算法，对初值比较敏感beta = LogisticRcgressNR(

14、Y.X.betaO)可以看到，自编函数与MATLAB自带函数Solve得到的结果相同，自编函数的缺点是对初 值敏感，没有编写对应的策略，可用GA、PSO等算法泄初始值。自编函数中用到的子函数：%子函数极大似然方程组function F = LogisticF(betaO,Y,X)n = length(Y); %样本个数 n = nl+n2XI = ones(nj) X;dims = size(Xl,2); %待估参数个数indl = (Y=l); % Y=1 的样本个数coll =sum(XKindi,:);%似然方程组F的第一部分col2 = zeros(dims,l); %似然方程组F的第

15、二部分初值%以下的向量表达好像不对% col2 = sum(X 1./(1 +exp(X 1 betaO)1;%fori = l:ncol2 = col2 + (Xl(i,:)/(l+exp(-Xl(i,:)*betaO)y；endF = coll - co!2J;% Newton-Rapson 算法中的 Jacobian 矩阵function M = LogisticJM(betaO,Y,X)n = length(Y); % 样本个数XI = ones(nJ) X;% 变量个数dims = size(Xl,2);M = zeros(dims);fori = l:nM = M + Xl(i,:

16、),*Xl(i,:)*exp( - X1 (i,:)*betaO)/( 1 +exp(-X 1 (i,:)*betaO)A2); endM = - M;function beta = LogisticRcgressNR(YX.bctaO)%用牛顿一拉普生方法求极大似然估计%Y因变量样本观测值，取值为1,表示事件发生，取值为0,表示事件不发生% X多个自变量的样本观测值，默认X的第一列全为1% betaO猜测的beta的初始值% 王福昌 2015-6-15%主函数itermax = 10; %最多迭代次数errstol = le-4; % 容忍误差iters = 0; %迭代次数% n. k =

17、 size(X); 诙n样本容量，k自变量个数 deltabeta = ones(size(betaO);betal = betaO + deltabeta;%虚拟迭代误差向量while (iterserrstol)deltabeta = LogisticJM(betaO,Y.X)LogisticF(betaO,Y,X);betal = betaO + deltabeta;betaO = beta 1; iters = iters +1;endbeta = betaO;附录2:直接优化对数似然函数，也能得到同样结果function F = LogisticRegressOpt(betarYX)

18、%用最优化方法求极大似然估计%Y因变量样本观测值，取值为1,表示事件发生，取值为0,表示事件不发生% X多个自变量的样本观测值，默认X的第一列全为1% betaO猜测的beta的初始值% 王福昌 2015-6-15%主函数%迭代次数% 极大似然函数indl = (Y=l);n = length(Y);XI = ones(nj) X;F = sum(X 1 (ind 1 ,:)*beta) - sum(Iog( l+exp(Xl*beta);F = -F;%假设我们有一个数据，45个观测值，四个变量，包括:% 1.% 2.% 3.% 4. 有)。age （年龄,数值型）；vision （视力状况

19、，分类型，1表示好，0表示有问题）：drive （驾车教冇，分类型，1表示参加过驾车教育，0表示没有）和一个分类型输出变M accident （去年是否出过事故，1表示出过事故，0表示没%我们的目的就是要考察前三个变量与发生事故的关系。% 第 1 至 4 列分别为 accident age vision drive ;data = 1 17 1 11 44 0 0 数据与附录1中相同1 16 1 01 61 1 0;Y = data(:,l);X= data(:,3:4);beta0=l;l;0;% 0.1110 ;1.7137 ;-1.5000;%rand(4J); %猎测的初始值I bet

20、a fval = fminsearch(beta) LogisticRcgressOpt(beta,Y.X),bctaO)可见，参数估计结果相同附录3:b =glmfit(X,Y,binomial*, link*, logit*)p = glmval(beta,X, logit)可以看到，与前面的两种方法结果相同。附录4：对比的例子proc logistic（logistic回归的SAS实现一无哑变量）（2012-03-13 21:30:54）转拔原文mt： proc logistic （logistic H归的SAS实现一无呵变放）作者，贝搭数据统i|工作Logistic阿归主耍用来处別应变

21、肚为分类变煲的何題比如生存和死亡患病和未忠擒乐当然研宪者咲心的问题是哪牛因索导 致了患釣或不忠衲.哪些足生存和死亡.Logistic的sas语句很简单.其基本语句见下：PROC LOGISTIC DATA=SAS-data-set ;MODEL response = independents ;BY variables;OUTPUT :WEIGHT variable:Proc logistic语句戏认计味应变址值眾小阴性结果一股赋值为0）的概率.但足常密我们想要得到的是阳性給果的概率即赋的数值的概率（二分类变fit时一般赋値为1）.选项Mdescending解决了这一何题。Model语句用干定

22、义应变贵和自变址：实例：锻设我们有一个数据 45个观测值.四个变fit.包括:age （年龄.数值型：2.vision （视力状况.分类型.1表示好.0茨示有何题：3.drive （驾车教育，分类空.1表示参加过驾车教育 0表示没有）和4.一个分类製输出变 accident （去年足否出过审故.1表示岀过爭故.0茨示没有）。我们的目的就足要考察前：个变壇与发生事故的关系。data logistic;input accident age vision drive:datalines :1 17 1 11 44 0 01 48 1 01 55 0 01 75 1 10 35 0 10 42 1 1

23、0 57000 28 0 10 20 0 10 38 1 00 45 0 10 47 1 10 52000 550 11 68 1 01 18 1 01 68 0 01 48 1 11 17001 70 1 11 72 1 01 35 0 11 19 1 01 62 1 00 39 1 10 40 1 10 55000 68 0 10 25 1 00 17000 45 0 10 44 0 10 67000 55 0 11 61 1 01 19 1 01 69 0 01 23 1 11 19001 72 1 11 74 1 01 31 0 11 16 1 01 61 1 0 run;（1） pr

24、oc logistic data=logistic descending;model accident=age vision drive： run ;如果想要在选择适当的自变fit筛选方法則使用一下语句：（2） proc logistic data=logistic descending;model accident=age vision drive/selection=：stepwise sle=0.15 sls=0.15 stb: run :Selection用于选择筛选自变就的方法.有backward （向后法）、forward （向前法）、stepwise （逐步法）、score （最

25、优子集 法）、none （完全法五个选项.戏认为none：SLE=概率值.入选标准.規定变掀入选怏型的显著性水平.前进法的炊认足0.5逐步浓足0.15SLS=ftt率值.剔除标准.指定变扯保留在怏里的显著水平.后退法炊认为0.10.逐步浓是0.15标准化偏回归系数STB可用 来比较备个自变fit作用的人小还可以输出贾信区何.语句如“（3） proc logistic data=logistic descending;model accident=age vision drive/selection=stepwise sle=0.15 sls=0.15 stb cl; run ;结呆：The L

26、OGISTIC ProcedureModel InformationData Set WORK.LOGISTICResponse VariableacddentNumber of Response Levels2Number of Observations45Modelbinary logitOptimization TechniqueFishers scori ngResponse ProfileOrdered Value accidentTotal Frequency11 2520 20Probability modeled is accident=d给岀了木枳电的慕本信息总思人釦I明。缶

27、要注您的足Response Profile中.accident=1排在忤位。帕而我们说过.SAS 的Logistic回归方楼log (odds) IV：认的形式是处典那个变磧(ft比较小的加J: descending选项后accident=1就挣在闫位了.(2) Forward Selection ProcedureStep 0. Intercept entered:Model Convergence StatusConvergence criterion (GCONV=1E-8) satisfied.Residual Chi-Square TestChi-Square DFPr ChiSq1

28、0.7057 30.0134Step 1 Effect vision entered:Model Convergence StatusCon vergence criterion (GCONV=1E-8) satisfied.Model Fit StatisticsCriterion Intercept Only Intercept and CovariatesSC 65.633 62.857-2 Log L 61.827 55.244Testing Global Null Hypothesis: BETA=0TestChi-SquareDFPr ChiSqLikelihood Ratio6.583010.0103Score6.420910.0113Wald6.075610.0137Residual Chi-Square TestChi-SquareDF Pr ChiSq4.981