SAS学习系列26. Logistic回归.doc

1、26. Logistic 回归（一）Logistic 回归一、原理二元或多元线性回归的因变量都是连续型变量，若因变量是分类变量（例如：患病与不患病；不重要、重要、非常重要） ，就需要用 Logistic 回归。Logistic 回归分析可以从统计意义上估计出在其它自变量固定不变的情况下，每个自变量对因变量取某个值的概率的数值影响大小。Logistic 回归模型有 “条件”与“非条件”之分，前者适用于配对病例对照资料的分析，后者适用于队列研究或非配对的病例-对照研究成组资料的分析。对于二分类因变量，y=1 表示事件发生；y=0 表示事件不发生。事件发生的条件概率 P y=1 | xi 与 xi

2、之间是非线性关系，通常是单调的，即随着 xi 的增加/减少，P y=1 | x i 也增加 /减少。Logistic 函数 F(x)= ，图形如下图所示： 11+该函数值域在(0,1)之间，x 趋于-时，F(x)趋于 0；x 趋于+时，F(x)趋于 1. 正好适合描述概率 P y=1 | xi . 例如，某因素 x 导致患病与否：x 在某一水平段内变化时，对患病概率的影响较大；而在 x较低或较高时对患病概率影响都不大。记事件发生的条件概率 P y=1 | xi = pi，则pi = =11+(+) +1+ 记事件不发生的条件概率为1- pi = 11+则在条件 xi 下，事件发生概率与事件不发

3、生概率之比为= 1- +称为事件的发生比，简记为 odds. 对 odds 取自然对数得到ln( 1- )= +上式左边（对数发生比）记为 Logit(y), 称为 y 的 Logit 变换。可见变换之后的 Logit(y)就可以用线性回归，计算出回归系数 和 值。若分类因变量 y 与多个自变量 xi 有关，则变换后 Logit(y)可由多元线性回归： 1logit()ln)kpx或 11()|, kkxpyxe二、回归参数的解释1. 三个名词发生比（odds）= = 事件 发 生 频 数事件未 发 生 频 数 1例如，事件发生概率为 0.6，不发生概率为 0.4，则发生比为 1.5（发生比1

4、，表示事件更可能发生） 。发生比率（OR）= = = = 12 1/(11)2/(12) 11/1221/2211221221即主对角线乘积/副对角线乘积，也称为交叉积比率，优势比。例如，说明：大于 1（小于 1）的发生比率，表明事件发生的可能性会提高（降低） ，或自变量对事件概率有正（负）的作用；发生比率为 1 表示变量对事件概率无作用。相对风险（RR ）= = 12 11/(11+12)21/(21+22)用来进行两组概率之间的比较。当 p1= p2 时，相对风险为 1，表明两组在事件发生方面没有差别。2. 连续型自变量回归参数的解释截距 ： 基准发生比的对数，即当 Logistic 回归

5、模型中没有任何自变量时（除常量外，所有自变量都取 0 值）所产生的发生比。由于理解发生比，比理解对数发生比更容易，故将 Logistic 回归模型改写为：odds = = 1 +11+= 11若 k0（ k1（ ;CLASS 分类变量;FREQ 频数变量;);MODEL 因变量 = 自变量列表 ;variable ;注：CLASS, EFFECT 语句必须在 MODEL 语句之前；CONTRAST, EXACT, ROC 语句必须在 MODEL 语句之后。说明：（1）输入数据集可选项DESCENDING 指定因变量按降序排序（“y=1”放前面）；ORDER= 指定因变量的排序顺序；PLOT 绘

6、图选项；（2）EFFECT 语句用原变量数据创建某种效应设计矩阵做对比用，例如 LAG 效应等。（3）CLASS 语句对分类变量进行 0-1 化处理，变成虚拟变量；（4）MODEL 语句是必不可少的，用来指定因变量和自变量。可以用可选项指定“y=1”，例如：model remiss(event=1) = cell smear infil li blast temp;可选项：selection=stepwise /forward/ backwardsle 或 sls指定变量进入或剔除出模型的显著水平；Aggregate 和 scale=n|p|d 计算偏差和 pearson 卡方拟和优度统计量，

7、n 表示对离差参数不进行校正；p 规定离差参数的估计为 pearson 卡方统计量除以自由度；d 规定离差参数的估计为偏差除以自由度；alpha= 设置置信限；cl/clparm=WALD估计所有参数/WALD 参数的置信区间；plrl对自变量估计比数比的置信区间；influence做回归诊断；RSQUARE输出拟合的调整的 R2；EFFECTPLOT输出模型拟合统计量；（5）ESTIMATE 语句 用来估计效应变量的线性组合的值；（6）EXACT 语句用其它变量的充分统计量对变量的充分统计量做精确检验；（7）CONTRAST 语句用来检验均值的线性组合关系的原假设。有三个基本参数，一是标签，

8、二是分类变量名，三是效应均值线性组合的系数表（系数的次序是匹配分类变量按字母数字次序的水平值） 。示例：contrast US vs NON-U.S. brand 1 1 -2;检验 H0： 1+2-23=0（8）ROC 语句绘制 ROC 曲线；（9）SCORE 语句输出若干结果到数据集。（10）TEST 语句对系数关系式做检验，示例：（ai 为自变量名）test1: test intercept + .5 * a2 = 0;test2: test intercept + .5 * a2;test3: test a1=a2=a3;test4: test a1=a2, a2=a3;例 1 不同治

9、疗方法对某病疗效的影响研究：用 Logistic 回归模型 Peffect=1 | treat = 拟合，即+1+ Logit(p) = ln( ) = 1 +代码：data effects;input treat effect count;cards;1 1 161 0 482 1 402 0 20;proc logistic data = effects DESCENDING; freq count; model effect=treat; /*或者用 model effect(event=1) = treat; 前面就可以不用DESCENDING 选项了*/run;运行结果及说明：响应概

10、况有序值effect 总频数1 1 562 0 68建模的概率为 effect=1。模型收敛状态满足收敛准则 (GCONV=1E-8)。模型拟合统计量准则 仅截距 截距和协变量AIC 172.737 152.361SC 175.558 158.001-2 L 170.737 148.361检验全局零假设: BETA=0检验 卡方 自由度 Pr 卡方似然比 22.3768 1 卡方Intercept 1 -2.8904 0.6390 20.4594 卡方偏差 0.2141 1 0.2141 0.6436Pearson 0.2155 1 0.2155 0.6425剩余差 D 和 Pearson 拟

11、合优度检验的 P 值分别为 0.6436 和0.6425 都远大于 0.05，故拟合结果较好可以接受。唯一轮廓数: 4模型拟合统计量准则 仅截距 截距和协变量AIC 109.669 101.900SC 112.026 108.970-2 L 107.669 95.900检验全局零假设: BETA=0检验 卡方 自由度 Pr 卡方似然比 11.7694 2 0.0028评分 11.2410 2 0.0036Wald 10.0644 2 0.0065假设检验 H0： 1= k=0. 似然比检验的 P 值=0.0028 卡方Intercept 1 1.1568 0.4036 8.2167 0.004

12、2sex 1 -1.2770 0.4980 6.5750 0.0103degree 1 -1.0545 0.4980 4.4844 0.0342拟合回归方程为：Logit(p)= 1.1568-1.2770*Sex 1.0545*Degree优比估计值效应 点估计值 95% Wald置信限sex 0.279 0.105 0.740degree 0.348 0.131 0.924优势比（OR）: 男性（Sex=1）治愈与未治愈的比值为：p1/(1-p1) = exp(1.1568-1.2770*1 1.0545*Degree)女性（Sex=0 ）治愈与未治愈的比值为：p0/(1-p0) = ex

13、p(1.1568-1.2770*0 1.0545*Degree)两个比值的比即优势比为：OR = p1/(1-p1)/ p0/(1-p0) = =e-1.2770 = 0.279e1预测概率和观测响应的关联一致部分所占百分比 60.0 Somers D 0.419不一致部分所占百分比 18.1 Gamma 0.536结值百分比 21.9 Tau-a 0.211对 1512 c 0.709Obs sex degree effect count _LEVEL_ prob1 0 0 1 21 1 0.760752 0 0 0 6 1 0.760753 0 1 1 9 1 0.525554 0 1 0

14、 9 1 0.525555 1 0 1 8 1 0.469996 1 0 0 10 1 0.469997 1 1 1 4 1 0.236018 1 1 0 11 1 0.23601输出预测概率。最大似然估计值分析参数 自由度 估计值 标准误差Wald 卡方 Pr 卡方Intercept 1 1.2528 0.4629 7.3239 0.0068sex 1 -1.4759 0.6628 4.9587 0.0260degree 1 -1.2528 0.6607 3.5954 0.0579sex*degree 1 0.4643 1.0012 0.2151 0.6428Sex*Degree 的卡方检验

15、 P 值=0.6428 远大于其它参数的 P 值，故不用考虑该交互作用的影响。例 3 多分类自变量的虚拟变量处理（CLASS 语句将会自动完成自变量的“虚拟变量化”处理过程，因变量不需要用 CLASS 语句处理），以及对比检验（CONTRAST 语句）。代码：data uti; input diagnosis : $13. treatment $ response $ count ; datalines; complicated A cured 78 complicated A not 28 complicated B cured 101 complicated B not 11 compli

16、cated C cured 68 complicated C not 46 uncomplicated A cured 40 uncomplicated A not 5 uncomplicated B cured 54 uncomplicated B not 5 uncomplicated C cured 34 uncomplicated C not 6 ; run; ods select FitStatistics; proc logistic data = uti; freq count; class diagnosis treatment /param=ref; model respon

17、se = diagnosis|treatment; run; ods select FitStatistics GoodnessOfFit TypeIII OddsRatios; proc logistic data = uti; freq count; class diagnosis treatment; model response = diagnosis treatment / scale=none aggregate; run; *clodds：计算似然比的置信区间 clparm: 计算参数的置信区间 ; ods select ClparmPL CloddsPL; proc logis

18、tic data = uti; freq count; class diagnosis treatment; model response = diagnosis treatment / scale=none aggregate clodds=pl clparm=pl; run; *contrast:定制假设检验的方式，变量需要是矩阵形式; ods select ContrastTest ContrastEstimate; proc logistic data = uti; freq count; class diagnosis treatment /param=ref; model resp

19、onse = diagnosis treatment; contrast B versus A treatment -1 1 / estimate=exp; contrast A treatment 1 0; contrast joint test treatment 1 0, treatment 0 1; run;运行结果及说明：（部分）偏差和 Pearson 拟合优度统计量准则 值 自由度 值/自由度 Pr 卡方偏差 2.5147 2 1.2573 0.2844Pearson 2.7574 2 1.3787 0.25193 型效应分析效应 自由度 Wald 卡方 Pr 卡方diagnosi

20、s 1 10.2885 0.0013treatment 2 24.6219 卡方B versus A 1 8.6919 0.0032A 1 4.9020 0.0268joint test 2 24.6219 卡方B versus AEXP 1 2.6539 0.8786 0.05 1.3870 5.0778 8.6919 0.0032例 4 自变量连续有序变量（例如，年龄）的 Logistic 回归。代码：data coronary; input sex ecg age ca ; datalines; 0 0 28 0 1 0 42 1 0 1 46 0 1 1 45 0 0 0 34 0 1

21、 0 44 1 0 1 48 1 1 1 45 1 0 0 38 0 1 0 45 0 0 1 49 0 1 1 45 1 0 0 41 1 1 0 46 0 0 1 49 0 1 1 46 1 0 0 44 0 1 0 48 0 0 1 52 0 1 1 48 1 0 0 45 1 1 0 50 0 0 1 53 1 1 1 57 1 0 0 46 0 1 0 52 1 0 1 54 1 1 1 57 1 0 0 47 0 1 0 52 1 0 1 55 0 1 1 59 1 0 0 50 0 1 0 54 0 0 1 57 1 1 1 60 1 0 0 51 0 1 0 55 0 0 2

22、46 1 1 1 63 1 0 0 51 0 1 0 59 1 0 2 48 0 1 2 35 0 0 0 53 0 1 0 59 1 0 2 57 1 1 2 37 1 0 0 55 1 1 1 32 0 0 2 60 1 1 2 43 1 0 0 59 0 1 1 37 0 1 0 30 0 1 2 47 1 0 0 60 1 1 1 38 1 1 0 34 0 1 2 48 1 0 1 32 1 1 1 38 1 1 0 36 1 1 2 49 0 0 1 33 0 1 1 42 1 1 0 38 1 1 2 58 1 0 1 35 0 1 1 43 0 1 0 39 0 1 2 59

23、1 0 1 39 0 1 1 43 1 1 0 42 0 1 2 60 1 0 1 40 0 1 1 44 1 ; run; *selection用于选择逐步回归方法,包括forward,backward,stepwiseinclude:设定每个拟合模型中包含 model语句中列的因子的个数. units:可以设置自变量每次变化10 个单位，计算的调整的发生比率AOR; proc logistic data = coronary descending; model ca = sex ecg age ecg*ecg age*age sex*ecg sex*age ecg*age / select

24、ion=forward include=3 details lackfit; run; proc logistic descending; model ca=sex ecg age; units age=10; run;运行结果及说明（部分）： 最大似然估计值分析参数 自由度 估计值 标准误差Wald 卡方 Pr 卡方Intercept 1 -5.6418 1.8061 9.7572 0.0018sex 1 1.3564 0.5464 6.1616 0.0131ecg 1 0.8732 0.3843 5.1619 0.0231age 1 0.0929 0.0351 7.0003 0.0081经

25、过筛选之后的 3 个自变量及截距的估计，回归方程略。优比估计值效应 点估计值 95% Wald置信限sex 3.882 1.330 11.330ecg 2.395 1.127 5.086age 1.097 1.024 1.175残差卡方检验卡方 自由度 Pr 卡方2.3277 5 0.8022对条目合格的效应分析效应 自由度 评分卡方 Pr 卡方ecg*ecg 1 0.3766 0.5394age*age 1 0.7712 0.3798sex*ecg 1 0.0352 0.8513sex*age 1 0.0290 0.8647ecg*age 1 0.8825 0.3475上述 5 个自变量的卡

26、方检验 P 值=0.05 ，故不必加入回归方程。Hosmer 和 Lemeshow 检验的分区ca = 1 ca = 0 组 总计观测值 期望值 观测值 期望值1 8 2 1.02 6 6.982 8 1 1.80 7 6.203 8 3 2.59 5 5.414 8 3 3.42 5 4.585 8 4 4.07 4 3.936 9 6 5.38 3 3.627 9 4 5.97 5 3.038 8 7 5.99 1 2.019 8 7 6.98 1 1.0210 4 4 3.77 0 0.23Hosmer 和 Lemeshow 拟合优度检验卡方 自由度 Pr 卡方4.7766 8 0.78

27、12优比效应 单位 估计值优比效应 单位 估计值age 10.0000 2.531自变量年龄按 10 岁为单位变化，得到的调整优势比为：AOR=2.531.例 5 回归诊断。代码：data uti2; input diagnosis : $13. treatment $ response trials; datalines; complicated A 78 106 complicated B 101 112 complicated C 68 114 uncomplicated A 40 45 uncomplicated B 54 59 uncomplicated C 34 40 ; *INF

28、LUENCE诊断; proc logistic data=uti2; class diagnosis treatment / param=ref; model response/trials = diagnosis treatment / influence; run; proc logistic data=uti2; class diagnosis treatment / param=ref; model response/trials = diagnosis / scale=none aggregate=(treatment diagnosis) influence iplots; run

29、; 运行结果及说明（略）例 6 精确检验。代码：data liver; input time $ group $ status $ count ; datalines; early antidote severe 6 early antidote not 12 early control severe 6 early control not 2 delayed antidote severe 3 delayed antidote not 4 delayed control severe 3 delayed control not 0 late antidote severe 5 late an

30、tidote not 1 late control severe 6 late control not 0 ; *estimate=both,表示对第一个exact语句中指定的变量进行精确点估计. joint,表示对第二个exact 中time、group进行联合检验 ; proc logistic descending; freq count; class time(ref=early) group(ref=control) / param=ref; model status = time group / scale=none aggregate clparm=wald; exact Mod

31、el 1 intercept time group / estimate=both; exact Joint Test time group / joint; run;运行结果（部分）：最大似然估计值分析参数 自由度 估计值 标准误差Wald 卡方 Pr 卡方Intercept 1 1.4132 0.7970 3.1439 0.0762time delayed 1 0.7024 0.8344 0.7087 0.3999time late 1 2.5533 1.1667 4.7893 0.0286group antidote 1 -2.2170 0.8799 6.3480 0.0118优比估计值

32、效应 点估计值 95% Wald置信限time early-delayed 2.019 0.393 10.359time early-late 12.849 1.305 126.471group control-antidote 0.109 0.019 0.611参数估计和 Wald 置信区间参数 估计值 95% 置信限Intercept 1.4132 -0.1489 2.9754time delayed 0.7024 -0.9330 2.3378参数估计和 Wald 置信区间参数 估计值 95% 置信限time late 2.5533 0.2666 4.8400group antidote

33、-2.2170 -3.9417 -0.4924例 7 绘制 ROC 曲线。ROC 曲线是根据一系列不同的二分类方式（分界值或决定阈），以真阳性率（灵敏度）为纵坐标，假阳性率（1-特异度）为横坐标绘制的曲线。1. ROC 曲线能很容易地查出任意界限值时的对疾病的识别能力。2. 选择最佳的诊断界限值。ROC 曲线越靠近左上角,试验的准确性就越高。最靠近左上角的 ROC 曲线的点是错误最少的最好阈值，其假阳性和假阴性的总数最少。3. 两种或两种以上不同诊断试验对疾病识别能力的比较。在对同一种疾病的两种或两种以上诊断方法进行比较时，可将各试验的ROC 曲线绘制到同一坐标中，以直观地鉴别优劣，靠近左上角

34、的ROC 曲线所代表的受试者工作最准确。亦可通过分别计算各个试验的 ROC 曲线下的面积(AUC) 进行比较，哪一种试验的 AUC 最大，则哪一种试验的诊断价值最佳。代码：data Data1;input disease n age;datalines;0 14 250 20 350 19 457 18 556 12 6517 17 75;ods graphics on;proc logistic data=Data1 plots(only)=roc(id=obs);model disease/n=age / scale=none clparm=wald clodds=pl rsquare;u

35、nits age=10;effectplot;run;ods graphics off;运行结果（部分）：偏差和 Pearson 拟合优度统计量准则 值 自由度 值/自由度 Pr 卡方偏差 7.7756 4 1.9439 0.1002Pearson 6.6020 4 1.6505 0.1585最大似然估计值分析参数 自由度 估计值 标准误差Wald 卡方 Pr 卡方Intercept 1 -12.5016 2.5555 23.9317 卡方似然比 30.2480 10 0.0008评分 28.3738 10 0.0016Wald 25.6828 10 0.00423 型效应分析效应 自由度 W

36、ald 卡方 Pr 卡方School 4 14.5522 0.0057Program 2 10.4815 0.0053School*Program 4 1.7439 0.7827最大似然估计值分析参数 Style 自由度估计值 标准误差Wald 卡方Pr 卡方Intercept self 1 -0.8097 0.1488 29.5989 .0001Intercept team 1 -0.6585 0.1366 23.2449 .0001School 1 self 1 -0.8194 0.2281 12.9066 0.0003School 1 team 1 -0.2675 0.1881 2.02

37、33 0.1549School 2 self 1 0.2974 0.1919 2.4007 0.1213School 2 team 1 -0.1033 0.1898 0.2961 0.5863Program regular self 1 0.3985 0.1488 7.1684 0.0074Program regular team 1 0.3537 0.1366 6.7071 0.0096School*Program 1 regular self 1 0.2751 0.2281 1.4547 0.2278School*Program 1 regular team 1 0.1474 0.1881 0.6143 0.4332School*Program 2 regular self 1 -0.0998 0.1919 0.2702 0.6032School*Program 2 regular team 1 -0.0168 0.1898 0.0079 0.9293