二项Logistic 回归参数最大似然估计的计算(共16页)

1、二项Logistic回归(hugu)参数(cnsh)最大似然估计的计算1 Logistic回归的基本(jbn)思想在地震风险评估中，研究者往往关心地震发生时，地表发生破裂的概率，地表破裂是由哪些因素引起的等问题。利用以往的相关数据找出统计规律性来解决这些问题，实质上可以转化为一个多元回归分析问题，其中，为随机变量。由于因变量的取值只有两个状态：破裂()和不破裂()，因此直接寻找因变量和自变量的关系非常困难。于是，可以把研究问题转换一个角度，不去直接分析和的关系，而是分析条件概率和的关系，这等价于寻找一个取值在0到1之间的连续函数。数学上满足这种条件的函数存在且不唯一，Logistic回归就是满

2、足这种要求的函数之一。和线性回归分析类似，Logistic 回归基本原理就是利用一组观测数据拟合一个Logistic模型，然后借助这个模型来揭示总体中若干个自变量与一个因变量取每个值的概率之间的依存关系，并评估用这一模型模拟相关事物变化规律的准确性。具体地说，Logistic回归分析可以从统计意义上确定在消除了其它变量的影响后，每一个自变量的变化是否引起因变量取某个值的概率的变化，并估计出在其它自变量固定不变的情况下，每个自变量对因变量取某个值的概率的数值影响大小。在使用Logistic回归分析前，需要明确模型的使用条件：1、要求因变量是分类变量，包括顺序变量和名义变量，不管哪种变量，都要用数

3、字表示它，如可以令表示地震发生时地表破裂，令表示地震发生时地表未破裂；2、自变量可以是(i)数值型连续变量，如地震的震级，(ii) 顺序变量，如覆盖层的厚度分组(10-20m，20-40m等)，(iii)名义变量，如地震类型，可令走滑型地震为1，正断型地震为2，逆冲型地震为3。2 多元二项Logistic回归模型的定义由于地震发生时地表是否破裂受到多个因素的影响，故引入多元Logistic回归模型。假设因变量是一个取值为1和0的二值变量，是影响的个因素，回归系数，则关于的元Logistic回归模型定义为 (1)由式(1)可得 (2)3 Logistic回归(hugu)参数估计我们用最大似然估计

4、方法(fngf)去求模型的参数。再假设(jish)从总体中抽取一个容量为的随机样本，其中，则有似然函数为 (3)两边取对数，整理可得 (4)写成向量形式为 (4)为求式(4)的驻点，可求对数似然函数关于的似然方程组为 (5)写成向量形式为式(5) 为非线性方程组，一般情况下没有解析解，可以用Newton-Raphson迭代(di di)方法求其数值解，令 (6)则关于(guny)的Jacobian矩阵(j zhn)为 (7)向量形式为 (7)根究Newton-Raphson方法的原理，可得参数迭代公式为 (8)算法(sun f)如下：Step 1: 给定(i dn)参数的初值参数(cnsh)和

5、误差容许精度，令；Step 2：计算；Step 3: 若或，即满足容许的精度，则结束，否则更新参数，转至Step2.当给定地震发生时，地表破裂是否发生的数据时，根据上面的算法，可以求出参数的最大似然估计。我们按照上述算法用MATLAB编写了多元Logistic回归参数估计的程序，可以给出参数估计值，详见附录。附录1用Newton-Raphson方法求解参数，附录2用直接优化对数似然函数方法求解参数，附录3用MATLAB自带的广义回归模型估计参数。附录4是别人做的Logistic回归的例子，用的软件是SAS（一种经过美国FDA认证的昂贵的商业统计软件）得到的结果。附录5是SPSS操作的过程和得到

6、的结果。附录1:Matlab Files for Logistic Regression% 假设我们有一个数据，45个观测值，四个变量，包括：% 1. age（年龄，数值型）；% 2. vision（视力状况，分类型，1表示好，0表示有问题）；% 3. drive（驾车教育，分类型，1表示参加过驾车教育，0表示没有） 和% 4. 一个分类型输出变量accident（去年是否出过事故，1表示出过事故，0表示没有）。% 我们的目的就是要考察前三个变量与发生事故的关系。 % 第1至4列分别为 accident age vision drive ; clear,clc,close alldata =

7、1 17 1 11 44 0 01 48 1 01 55 0 01 75 1 10 35 0 10 42 1 10 57 0 00 28 0 10 20 0 10 38 1 00 45 0 10 47 1 10 52 0 00 55 0 11 68 1 01 18 1 01 68 0 01 48 1 11 17 0 01 70 1 11 72 1 01 35 0 11 19 1 01 62 1 00 39 1 10 40 1 10 55 0 00 68 0 10 25 1 00 17 0 00 45 0 10 44 0 10 67 0 00 55 0 11 61 1 01 19 1 01 69

8、 0 01 23 1 11 19 0 01 72 1 11 74 1 01 31 0 11 16 1 01 61 1 0;Y = data(:,1);X= data(:,3:4);beta0 = 0.110 ;1.7137 ; -1.5000+1*rand(3,1);%rand(4,1); %猜测(cic)的初始值%自带的fsolve 算法求解没有(mi yu)问题，核心原理也是Newton-Raphson 算法%options = optimset(Display,iter,TolFun,1e-8)beta = fsolve(be)LogisticF(be,Y,X),beta0)%,opti

9、ons) % 自编的Newton-Raphson 算法(sun f)，对初值比较敏感beta = LogisticRegressNR(Y,X,beta0)可以看到，自编函数与MATLAB自带函数Solve得到(d do)的结果相同，自编函数的缺点是对初值敏感，没有编写对应的策略，可用GA、PSO等算法定初始值。自编函数(hnsh)中用到的子函数：%子函数 极大(j d)似然方程组function F = LogisticF(beta0,Y,X)n = length(Y); %样本个数 n = n1+n2X1 = ones(n,1) X;dims = size(X1,2); %待估参数个数 in

10、d1 = (Y=1); % Y=1的样本个数 col1 = sum(X1(ind1,:); % 似然方程组 F 的第一部分 col2 = zeros(dims,1); % 似然方程组 F 的第二部分初值% 以下的向量表达好像不对% col2 = sum(X1./(1+exp(X1*beta0);% for i = 1:n col2 = col2 + (X1(i,:)/(1+exp(-X1(i,:)*beta0);endF = col1 - col2;% Newton-Rapson算法中的 Jacobian 矩阵function M = LogisticJM(beta0,Y,X)n = leng

11、th(Y); % 样本个数X1 = ones(n,1) X; % 变量个数dims = size(X1,2);M = zeros(dims);for i = 1:n M = M + X1(i,:)*X1(i,:)*exp( - X1(i,:)*beta0)/(1+exp(-X1(i,:)*beta0)2);endM = - M;function beta = LogisticRegressNR(Y,X,beta0)% 用牛顿(ni dn)拉普生方法(fngf)求极大似然估计% Y 因变量样本观测(gunc)值，取值为1，表示事件发生，取值为0，表示事件不发生% X 多个自变量的样本观测值, 默

12、认X的第一列全为1% beta0 猜测的beta 的初始值% % 王福昌 2015 - 6 - 15%主函数itermax = 10; % 最多迭代次数errstol = 1e-4; % 容忍误差iters = 0; % 迭代次数% n, k = size(X); % n 样本容量， k 自变量个数deltabeta = ones(size(beta0); beta1 = beta0 + deltabeta;% 虚拟迭代误差向量while (iterserrstol) deltabeta = -LogisticJM(beta0,Y,X)LogisticF(beta0,Y,X); beta1 =

13、 beta0 + deltabeta; beta0 = beta1; iters = iters +1;endbeta = beta0;附录2: 直接优化对数似然函数，也能得到同样结果function F = LogisticRegressOpt(beta,Y,X)% 用最优化方法求极大似然估计% Y 因变量样本观测值，取值为1，表示事件发生，取值为0，表示事件不发生% X 多个自变量的样本观测值, 默认X的第一列全为1% beta0 猜测的beta 的初始值% % 王福昌 2015 - 6 - 15%主函数% 迭代次数% 极大似然函数ind1 = (Y=1);n = length(Y);X1

14、 = ones(n,1) X;F = sum(X1(ind1,:)*beta) - sum(log( 1+exp(X1*beta);F = -F;% 假设我们有一个数据，45个观测值，四个变量，包括：% 1. age（年龄，数值型）；% 2. vision（视力状况，分类型，1表示好，0表示有问题）；% 3. drive（驾车教育，分类型，1表示参加过驾车教育，0表示没有） 和% 4. 一个(y )分类型输出变量accident（去年是否出过事故，1表示出过事故，0表示没有）。% 我们的目的就是(jish)要考察前三个变量与发生事故的关系。 % 第1至4列分别(fnbi)为 accident

15、age vision drive ;data = 1 17 1 11 44 0 0数据与附录1中相同1 16 1 01 61 1 0;Y = data(:,1);X= data(:,3:4);beta0 = 1;1;0;% 0.1110 ;1.7137 ; -1.5000;%rand(4,1); %猜测的初始值beta fval = fminsearch(beta) LogisticRegressOpt(beta,Y,X),beta0)可见，参数估计结果相同附录3:b =glmfit(X,Y,binomial, link, logit)p = glmval(beta,X, logit)可以看到

16、，与前面的两种方法结果相同。附录(fl)4：对比的例子proclogistic(logistic回归(hugu)的SAS实现-无哑变量)(2012-03-13 21:30:54) HYPERLINK javascript:; 转载(zhunzi)原文地址： HYPERLINK /s/blog_8baaaca501010cj9.html o proclogistic(logistic回归的SAS实现-无哑变量) t _blank proclogistic(logistic回归的SAS实现-无哑变量)作者： HYPERLINK /u/2343218341 o 贝塔数据统计工作 t _blank 贝

17、塔数据统计工作Logistic回归主要用来处理应变量为分类变量的问题，比如生存和死亡，患病和未患病等。当然研究者关心的问题是哪些因素导致了患病或不患病，哪些是生存和死亡。Logistic的sas语句很简单，其基本语句见下：PROC LOGISTIC DATA=SAS-data-set ;MODEL response = independents ;BY variables;OUTPUT / ;WEIGHT variable;Proc logistic语句默认计算应变量值最小（阴性结果-一般赋值为0）的概率，但是常常我们想要得到的是阳性结果的概率，即赋值最大的数值的概率（二分类变量时一般赋值为1

18、），选项“descending”解决了这一问题。Model语句用于定义应变量和自变量；实例(shl)：假设我们有一个数据，45个观测(gunc)值，四个变量，包括：1.age（年龄(ninlng)，数值型）；2.vision（视力状况，分类型，1表示好，0表示有问题）；3.drive（驾车教育，分类型，1表示参加过驾车教育，0表示没有）和4.一个分类型输出变量accident（去年是否出过事故，1表示出过事故，0表示没有）。我们的目的就是要考察前三个变量与发生事故的关系。datalogistic;inputaccident age vision drive;datalines;1 17 1 1

19、1 44 0 01 48 1 01 55 0 01 75 1 10 35 0 10 42 1 10 57 0 00 28 0 10 20 0 10 38 1 00 45 0 10 47 1 10 52 0 00 55 0 11 68 1 01 18 1 01 68 0 01 48 1 11 17 0 01 70 1 11 72 1 01 35 0 11 19 1 01 62 1 00 39 1 10 40 1 10 55 0 00 68 0 10 25 1 00 17 0 00 45 0 10 44 0 10 67 0 00 55 0 11 61 1 01 19 1 01 69 0 01 23

20、 1 11 19 0 01 72 1 11 74 1 01 31 0 11 16 1 01 61 1 0;run;(1)proclogisticdata=logisticdescending;modelaccident=age vision drive;run;如果想要在选择适当的自变量筛选方法则使用(shyng)一下语句：(2)proclogisticdata=logisticdescending;modelaccident=age vision drive/selection=stepwise sle=0.15 sls=0.15 stb;run;Selection用于选择筛选(shixun

21、)自变量的方法，有backward（向后法）、forward（向前(xin qin)法）、stepwise（逐步法）、score（最优子集法）、none（完全法）五个选项，默认为none；SLE=概率值，入选标准，规定变量入选模型的显著性水平，前进法的默认是0.5，逐步法是0.15SLS=概率值，剔除标准，指定变量保留在模型的显著水平，后退法默认为0.10，逐步法是0.15标准化偏回归系数STB可用来比较各个自变量作用的大小还可以输出(shch)置信区间，语句如下：(3)proclogisticdata=logisticdescending;modelaccident=age vision d

22、rive/selection=stepwise sle=0.15 sls=0.15 stbcl;run;结果(ji gu)：The LOGISTIC ProcedureModel InformationData SetWORK.LOGISTICResponse VariableaccidentNumber of Response Levels2Number of Observations45Modelbinary logitOptimization TechniqueFishers scoringResponse ProfileOrderedValueaccidentTotal Frequen

23、cy11252020Probability modeled is accident=1.(1)给出了本模型的基本信息，意思大多自明。需要(xyo)注意的是Response Profile中，accident=1排在首位。前面我们说过，SAS的Logistic回归方程log(odds)默认的形式是处理那个变量值比较小的，加上descending选项后，accident=1就排在首位了。(2)Forward Selection ProcedureStep 0. Intercept entered:Model Convergence StatusConvergence criterion (GCON

24、V=1E-8) satisfied.Residual Chi-Square TestChi-SquareDFPr ChiSq10.705730.0134Step 1. Effectvisionentered:Model Convergence StatusConvergence criterion (GCONV=1E-8) satisfied.Model Fit StatisticsCriterionInterceptOnlyInterceptand CovariatesAIC63.82759.244SC65.63362.857-2 Log L61.82755.244Testing Globa

25、l Null Hypothesis: BETA=0TestChi-SquareDFPr ChiSqLikelihood Ratio6.583010.0103Score6.420910.0113Wald6.075610.0137Residual Chi-Square TestChi-SquareDFPr ChiSq4.981820.0828Step 2. Effectdriveentered:Model Convergence StatusConvergence criterion (GCONV=1E-8) satisfied.Model Fit StatisticsCriterionInter

26、ceptOnlyIntercept and CovariatesAIC63.82756.287SC65.63361.707-2 Log L61.82750.287Testing Global Null Hypothesis: BETA=0TestChi-SquareDFPr ChiSqLikelihood Ratio11.539120.0031Score10.597620.0050Wald8.594920.0136Residual Chi-Square TestChi-SquareDFPr ChiSq0.129310.7191NOTE: No (additional) effects met

27、the 0.05 significance level for entry into the model.(2)给出了自变量进入模型的次序(cx)。先是截距项Step 0了，不管(bgun)它。Step 1vision第一个进入模型，附带了很多评估(pn )它对因变量预测能力的指标。-2 Log L和Score用来检测自变量是否显著。-2 Log L中的L就是Likelihood Ratio，它的p值是0.0103，Score的p值是0.0113，都小于0.05，故vision是一个很显著的解释变量。还有，AIC(Akaike Information Criterion)和SC(Schwarz

28、 Criterion)两个信息量标准用来比较不同的模型，它们数值越小，模型变现就越好，这个接下来我们看step2 drive变量进入模型后的情况。我们可以看到模型的表现变好了，因为这是AIC和SC的值变小了，-2 Log L和Score对应的p值也更小。(3)Summary of Forward SelectionStepEffect EnteredDFNumber In Score Chi-SquarePr ChiSq1vision116.42090.01132drive124.86800.0274(3)总结了我们模型使用的前向选择方法，包括自变量进入模型的次序，以及每个自变量的卡方值和p值

29、。(4)Analysis of Maximum Likelihood EstimatesParameter DFEstimateStandardError Wald Chi-SquarePr ChiSqIntercept10.11100.54570.04140.8389vision11.71370.70495.91130.0150drive1-1.50000.70374.54400.0330(4)给出了模型参数的估计，据此可以写出改回归方程的形式是log(odds of having an accident)=log(p/(1-p)=0.1110+1.7137*vision-1.5000*drive。我们知道，odds=p/(1- p),有p=odds/(1+odds)。假设有个哥们，视力没问题但没有受过驾车教育（vision=0,drive=0），代入方程，有log(odds)=0.1110，再odds=exp(0.110)=1.1174，p=1.1174/2.1174=0.5277，即我们说这人发生事故的概率为0.5277；又另一个，视力有问题同样没受过驾车教育